diff --git a/drivers/ccpp/noahmpdrv.F90 b/drivers/ccpp/noahmpdrv.F90 index 6aff506..db3ae47 100644 --- a/drivers/ccpp/noahmpdrv.F90 +++ b/drivers/ccpp/noahmpdrv.F90 @@ -157,7 +157,7 @@ subroutine noahmpdrv_run & sncovr1, qsurf, gflux, drain, evap, hflx, ep, runoff, & cmm, chh, evbs, evcw, sbsno, pah, ecan, etran, edir, snowc,& stm, snohf,smcwlt2, smcref2, wet1, t2mmp, q2mp,zvfun, & - ztmax, errmsg, errflg, & + ztmax, rca, errmsg, errflg, & canopy_heat_storage_ccpp, & rainfall_ccpp, & sw_absorbed_total_ccpp, & @@ -400,6 +400,8 @@ subroutine noahmpdrv_run & real(kind=kind_phys), dimension(:) , intent(out) :: q2mp ! combined q2m from tiles real(kind=kind_phys), dimension(:) , intent(out) :: zvfun ! real(kind=kind_phys), dimension(:) , intent(out) :: ztmax ! thermal roughness length + real(kind=kind_phys), dimension(:) , intent(out) :: rca ! total canopy/stomatal resistance (s/m) + character(len=*) , intent(out) :: errmsg integer , intent(out) :: errflg @@ -623,7 +625,7 @@ subroutine noahmpdrv_run & real (kind=kind_phys) :: canopy_heat_storage ! out | within-canopy heat [W/m2] real (kind=kind_phys) :: spec_humid_sfc_veg ! out | surface specific humidty over vegetation [kg/kg] real (kind=kind_phys) :: spec_humid_sfc_bare ! out | surface specific humidty over bare soil [kg/kg] - + real (kind=kind_phys) :: ustarx ! inout |surface friction velocity real (kind=kind_phys) :: prslkix ! in exner function real (kind=kind_phys) :: prsik1x ! in exner function @@ -948,6 +950,10 @@ subroutine noahmpdrv_run & ch_vegetated = 0.0 ch_bare_ground = ch_noahmp canopy_heat_storage = 0.0 + lai_sunlit = 0.0 + lai_shaded = 0.0 + rs_sunlit = 0.0 + rs_shaded = 0.0 else ! not glacier @@ -1056,7 +1062,17 @@ subroutine noahmpdrv_run & chxy (i) = ch_noahmp zorl (i) = z0_total * 100.0 ! convert to cm ztmax (i) = z0h_total - + + !LAI-scale canopy resistance based on weighted sunlit shaded fraction + if(rs_sunlit .le. 0.0 .or. rs_shaded .le. 0.0 .or. & + lai_sunlit .eq. 0.0 .or. lai_shaded .eq. 0.0) then + rca(i) = parameters%rsmax + else !calculate LAI-scale canopy conductance (1/Rs) + rca(i) = ((1.0/(rs_sunlit+leaf_air_resistance)*lai_sunlit) + & + ((1.0/(rs_shaded+leaf_air_resistance))*lai_shaded)) + rca(i) = max((1.0/rca(i)),parameters%rsmin) !resistance + end if + smc (i,:) = soil_moisture_vol slc (i,:) = soil_liquid_vol snowxy (i) = float(snow_levels) diff --git a/drivers/ccpp/noahmpdrv.meta b/drivers/ccpp/noahmpdrv.meta index 39eed14..7fa13f3 100644 --- a/drivers/ccpp/noahmpdrv.meta +++ b/drivers/ccpp/noahmpdrv.meta @@ -1360,6 +1360,14 @@ type = real kind = kind_phys intent = out +[rca] + standard_name = aerodynamic_resistance_in_canopy + long_name = canopy resistance + units = s m-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/drivers/ccpp/sfc_diff.f b/drivers/ccpp/sfc_diff.f index c5ed8bf..e5cb429 100644 --- a/drivers/ccpp/sfc_diff.f +++ b/drivers/ccpp/sfc_diff.f @@ -62,6 +62,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) & flag_iter,redrag, & !intent(in) & flag_lakefreeze, & !intent(in) & u10m,v10m,sfc_z0_type, & !hafs,z0 type !intent(in) + & u1,v1,usfco,vsfco,icplocn2atm, & & wet,dry,icy, & !intent(in) & thsfc_loc, & !intent(in) & tskin_wat, tskin_lnd, tskin_ice, & !intent(in) @@ -86,6 +87,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) integer, parameter :: kp = kind_phys integer, intent(in) :: im, ivegsrc integer, intent(in) :: sfc_z0_type ! option for calculating surface roughness length over ocean + integer, intent(in) :: icplocn2atm ! option for including ocean current in the computation of flux integer, dimension(:), intent(in) :: vegtype @@ -97,6 +99,8 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) logical, intent(in) :: thsfc_loc ! Flag for reference pressure in theta calculation real(kind=kind_phys), dimension(:), intent(in) :: u10m,v10m + real(kind=kind_phys), dimension(:), intent(in) :: u1,v1 + real(kind=kind_phys), dimension(:), intent(in) :: usfco,vsfco real(kind=kind_phys), intent(in) :: rvrdm1, eps, epsm1, grav real(kind=kind_phys), dimension(:), intent(in) :: & & ps,t1,q1,z1,garea,prsl1,prslki,prsik1,prslk1, & @@ -127,6 +131,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) ! locals ! integer i + real(kind=kind_phys) :: windrel ! real(kind=kind_phys) :: rat, tv1, thv1, restar, wind10m, & czilc, tem1, tem2, virtfac @@ -351,11 +356,29 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) & * virtfac endif - z0 = 0.01_kp * z0rl_wat(i) - z0max = max(zmin, min(z0,z1(i))) -! ustar_wat(i) = sqrt(grav * z0 / charnock) - wind10m = sqrt(u10m(i)*u10m(i)+v10m(i)*v10m(i)) + if (icplocn2atm == 0) then + wind10m=sqrt(u10m(i)*u10m(i)+v10m(i)*v10m(i)) + windrel=wind(i) + else if (icplocn2atm ==1) then + wind10m=sqrt((u10m(i)-usfco(i))**2+(v10m(i)-vsfco(i))**2) + windrel=sqrt((u1(i)-usfco(i))**2+(v1(i)-vsfco(i))**2) + endif + + if (sfc_z0_type == -1) then ! using wave model derived momentum roughness + tem1 = 0.11 * vis / ustar_wat(i) + z0 = tem1 + 0.01_kp * z0rl_wav(i) + if (redrag) then + z0max = max(min(z0, z0s_max),1.0e-7_kp) + else + z0max = max(min(z0,0.1_kp), 1.0e-7_kp) + endif + z0rl_wat(i) = 100.0_kp * z0max ! cm + else + z0 = 0.01_kp * z0rl_wat(i) + z0max = max(zmin, min(z0,z1(i))) + endif +! !** test xubin's new z0 ! ztmax = z0max @@ -385,7 +408,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) ! call stability ! --- inputs: - & (z1(i), zvfun(i), gdx, tv1, thv1, wind(i), + & (z1(i), zvfun(i), gdx, tv1, thv1, windrel, & z0max, ztmax_wat(i), tvs, grav, thsfc_loc, ! --- outputs: & rb_wat(i), fm_wat(i), fh_wat(i), fm10_wat(i), fh2_wat(i), @@ -425,17 +448,18 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) z0rl_wat(i) = 1.0e-4_kp endif - elseif (z0rl_wav(i) <= 1.0e-7_kp .or. & - & z0rl_wav(i) > 1.0_kp) then -! z0 = (charnock / grav) * ustar_wat(i) * ustar_wat(i) - tem1 = 0.11 * vis / ustar_wat(i) - z0 = tem1 + (charnock/grav)*ustar_wat(i)*ustar_wat(i) + elseif (z0rl_wav(i) <= 1.0e-7_kp .or. + & z0rl_wav(i) > 1.0_kp) then +! z0 = (charnock / grav) * ustar_wat(i) * ustar_wat(i) + tem1 = 0.11 * vis / ustar_wat(i) + z0 = tem1 + (charnock/grav)*ustar_wat(i)*ustar_wat(i) + + if (redrag) then + z0rl_wat(i) = 100.0_kp * max(min(z0, z0s_max),1.0e-7_kp) + else + z0rl_wat(i) = 100.0_kp * max(min(z0,0.1_kp), 1.0e-7_kp) + endif - if (redrag) then - z0rl_wat(i) = 100.0_kp * max(min(z0, z0s_max),1.0e-7_kp) - else - z0rl_wat(i) = 100.0_kp * max(min(z0,0.1_kp), 1.0e-7_kp) - endif endif endif ! end of if(open ocean) diff --git a/drivers/nuopc/lnd_comp_driver.F90 b/drivers/nuopc/lnd_comp_driver.F90 index eec9fe3..a788815 100644 --- a/drivers/nuopc/lnd_comp_driver.F90 +++ b/drivers/nuopc/lnd_comp_driver.F90 @@ -660,7 +660,7 @@ subroutine drv_run(gcomp, noahmp, rc) noahmp%model%snowc , noahmp%model%stm , & noahmp%model%snohf , noahmp%model%smcwlt2 , noahmp%model%smcref2 , & noahmp%model%wet1 , noahmp%model%t2mmp , noahmp%model%q2mp , & - noahmp%model%zvfun , noahmp%model%ztmax , & + noahmp%model%zvfun , noahmp%model%ztmax , noahmp%model%rca , & noahmp%static%errmsg , noahmp%static%errflg) !---------------------- diff --git a/drivers/nuopc/lnd_comp_types.F90 b/drivers/nuopc/lnd_comp_types.F90 index bec9fbd..bad5df7 100644 --- a/drivers/nuopc/lnd_comp_types.F90 +++ b/drivers/nuopc/lnd_comp_types.F90 @@ -202,6 +202,7 @@ module lnd_comp_types real(kind=r8), allocatable :: q2mp (:) ! combined q2m from tiles real(kind=r8), allocatable :: zvfun (:) ! some function of vegetation used for gfs stability real(kind=r8), allocatable :: ztmax (:) ! bounded surface roughness length for heat over land + real(kind=r8), allocatable :: rca (:) ! total canopy/stomatal resistance (s/m) real(kind=r8), allocatable :: rho (:) ! air density real(kind=r8), allocatable :: pores (:) ! max soil moisture for a given soil type for land surface model real(kind=r8), allocatable :: resid (:) ! min soil moisture for a given soil type for land surface model @@ -493,6 +494,7 @@ subroutine InitializeAllocate(this, begl, endl, km, lsnowl) allocate(this%model%q2mp (begl:endl)) allocate(this%model%zvfun (begl:endl)) allocate(this%model%ztmax (begl:endl)) + allocate(this%model%rca (begl:endl)) allocate(this%model%rho (begl:endl)) allocate(this%model%pores (30)) allocate(this%model%resid (30)) @@ -661,6 +663,7 @@ subroutine InitializeDefault(this) this%model%q2mp = 0.0_r8 this%model%zvfun = 0.0_r8 this%model%ztmax = 0.0_r8 + this%model%rca = 0.0_r8 this%model%rho = 0.0_r8 this%model%pores = 0.0_r8 this%model%resid = 0.0_r8 diff --git a/src/module_sf_noahmplsm.F90 b/src/module_sf_noahmplsm.F90 index 6abd59f..2472f11 100644 --- a/src/module_sf_noahmplsm.F90 +++ b/src/module_sf_noahmplsm.F90 @@ -2013,6 +2013,8 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in chuc = 0. chv2 = 0. rb = 0. + laisun = 0. + laisha = 0. cdmnv = 0.0 ezpdv = 0.0 @@ -2263,7 +2265,8 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in csigmaf1, & !out !jref:start qc ,qsfc ,psfc , & !in - q2v ,chv2, chleaf, chuc) !inout + q2v ,chv2 ,chleaf ,chuc , & + rb) !out ! new coupling code @@ -3712,7 +3715,8 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & t2mv ,psnsun ,psnsha ,canhs , & !out csigmaf1, & !out qc ,qsfc ,psfc , & !in - q2v ,cah2 ,chleaf ,chuc ) !inout + q2v ,cah2 ,chleaf ,chuc , & !inout + rb) !out ! -------------------------------------------------------------------------------------------------- ! use newton-raphson iteration to solve for vegetation (tv) and @@ -3836,6 +3840,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys), intent(out) :: chuc !< under canopy exchange coefficient real (kind=kind_phys), intent(out) :: canhs !< canopy heat storage change (w/m2) real (kind=kind_phys), intent(out) :: q2v !< + real (kind=kind_phys), intent(out) :: rb !< bulk leaf boundary layer resistance (s/m) real (kind=kind_phys) :: cah !< sensible heat conductance, canopy air to zlvl air (m/s) real (kind=kind_phys) :: u10v !< 10 m wind speed in eastward dir (m/s) real (kind=kind_phys) :: v10v !< 10 m wind speed in eastward dir (m/s) @@ -3852,7 +3857,6 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys) :: z0mo !roughness length for intermediate output only (m) real (kind=kind_phys) :: z0h !roughness length, sensible heat (m) real (kind=kind_phys) :: z0hg !roughness length, sensible heat (m) - real (kind=kind_phys) :: rb !bulk leaf boundary layer resistance (s/m) real (kind=kind_phys) :: ramc !aerodynamic resistance for momentum (s/m) real (kind=kind_phys) :: rahc !aerodynamic resistance for sensible heat (s/m) real (kind=kind_phys) :: rawc !aerodynamic resistance for water vapor (s/m)