u_profile_val_func Subroutine

    subroutine u_profile_val_func(z,L,u_val,z_val_in,z_pbl,z0,u,u_star0,u_pbl)

        implicit none

        real, intent(in) :: z,L,u_val,z_val_in,z_pbl,z0
        real, intent(out) :: u,u_star0,u_pbl
        real a,b,p,kappa
        parameter (a=16.,b=5.,p=-0.25,kappa=0.4)
        real z_l,z_val
        real phim,phim_i
        real phim_val,phim_i_val
        !real phim_pbl,phih_pbl,phim_i_pbl,phih_i_pbl

        !If the input height is above the boundary layer then set the height to pbl height and calculate
        z_val = z_val_in
        if (z_val.ge.z_pbl) z_val = z_pbl

        z_l=0.4*z_pbl

        call phim_func(z,L,phim,phim_i)
        call phim_func(z_val,L,phim_val,phim_i_val)

        !call phi_func(z_pbl,L,phim_pbl,phih_pbl,phim_i_pbl,phih_i_pbl)

        if (L.ge.0.) then
            u_star0=u_val*kappa/(log(z_val/z0)-phim_i_val+kappa*z_val/z_l*(1-z_val/2./z_pbl)-z_val/z_pbl*(1+b*z_val/2./L))
            !u_pbl=u_star0/kappa*(log(z_pbl/z0)-phim_i_pbl+kappa*z_pbl/z_l*(1-z_pbl/2./z_pbl)-z_pbl/z_pbl*(1+b*z_pbl/2./L))
            u=u_star0/kappa*(log(z/z0)-phim_i+kappa*z/z_l*(1-z/2./z_pbl)-z/z_pbl*(1+b*z/2./L))
        else
            u_star0=u_val*kappa/(log(z_val/z0)-phim_i_val+kappa*z_val/z_l*(1-z_val/2./z_pbl)-1./z_pbl*((a*z_val-L)*phim_val+L)/a/(p+1))
            !u_pbl=u_star0/kappa*(log(z_pbl/z0)-phim_i_pbl+kappa*z_pbl/z_l*(1-z_pbl/2./z_pbl)-1./z_pbl*((a*z_pbl-L)*phim_pbl+L)/a/(p+1))
            u=u_star0/kappa*(log(z/z0)-phim_i+kappa*z/z_l*(1-z/2./z_pbl)-1./z_pbl*((a*z-L)*phim+L)/a/(p+1))
        endif

        !u_pbl not used and so not calculated here
        u_pbl=u
        if (z.ge.z_pbl) then
            u=u_pbl
        endif

    end subroutine u_profile_val_func