From 12aef874d4b2064f3489e1637e6b182b0d70d8e8 Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Thu, 9 May 2024 12:35:14 -0400 Subject: [PATCH 1/4] add canopy resistance to noahmp%model DDT and pass as argument to noahmpdrv_run --- drivers/nuopc/lnd_comp_driver.F90 | 2 +- drivers/nuopc/lnd_comp_types.F90 | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) 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 From 0cc58c037ced60b5eb8268030b383739e373aa09 Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Thu, 9 May 2024 19:17:53 +0000 Subject: [PATCH 2/4] update NOAHMP files from CCPP to coincide with ccpp-physics PR#205 --- drivers/ccpp/noahmpdrv.F90 | 22 ++++++++++++--- drivers/ccpp/noahmpdrv.meta | 8 ++++++ drivers/ccpp/sfc_diff.f | 54 ++++++++++++++++++++++++++----------- src/module_sf_noahmplsm.F90 | 10 ++++--- 4 files changed, 73 insertions(+), 21 deletions(-) diff --git a/drivers/ccpp/noahmpdrv.F90 b/drivers/ccpp/noahmpdrv.F90 index 6aff506..f2f6b69 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 - + + ! total stomatal/canopy resistance Based on Bonan et al. (2011) conductance (1/Rs) equation + 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) = 0.0 + else + rca(i) = ((1.0/(rs_sunlit+leaf_air_resistance)*lai_sunlit) + & + ((1.0/(rs_shaded+leaf_air_resistance))*lai_shaded)) + rca(i) = 1.0/rca(i) !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/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) From 1dbf99ab7910e8880e2d6f20f479a3858203d6e2 Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Wed, 15 May 2024 12:28:50 -0400 Subject: [PATCH 3/4] update noahmpdrv.F90 from ccpp/physics PR#205 --- drivers/ccpp/noahmpdrv.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/drivers/ccpp/noahmpdrv.F90 b/drivers/ccpp/noahmpdrv.F90 index f2f6b69..db3ae47 100644 --- a/drivers/ccpp/noahmpdrv.F90 +++ b/drivers/ccpp/noahmpdrv.F90 @@ -1063,14 +1063,14 @@ subroutine noahmpdrv_run & zorl (i) = z0_total * 100.0 ! convert to cm ztmax (i) = z0h_total - ! total stomatal/canopy resistance Based on Bonan et al. (2011) conductance (1/Rs) equation + !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) = 0.0 - else + 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) = 1.0/rca(i) !resistance + rca(i) = max((1.0/rca(i)),parameters%rsmin) !resistance end if smc (i,:) = soil_moisture_vol From d696ac55551e850e564640897de1d5159b3db5a3 Mon Sep 17 00:00:00 2001 From: Michael Barlage Date: Mon, 8 Jul 2024 18:22:02 +0000 Subject: [PATCH 4/4] update three relevant file for consistency with ccpp HR4 updates --- drivers/ccpp/noahmpdrv.F90 | 97 +++++++++++++++++++------------------ drivers/ccpp/noahmpdrv.meta | 43 ++++++++++++++++ src/module_sf_noahmplsm.F90 | 50 +++++++++---------- 3 files changed, 117 insertions(+), 73 deletions(-) diff --git a/drivers/ccpp/noahmpdrv.F90 b/drivers/ccpp/noahmpdrv.F90 index db3ae47..1313e9f 100644 --- a/drivers/ccpp/noahmpdrv.F90 +++ b/drivers/ccpp/noahmpdrv.F90 @@ -292,11 +292,11 @@ subroutine noahmpdrv_run & integer , intent(in) :: iyrlen ! year length [days] real(kind=kind_phys) , intent(in) :: julian ! julian day of year real(kind=kind_phys), dimension(:) , intent(in) :: garea ! area of the grid cell - real(kind=kind_phys), dimension(:) , intent(in) :: rainn_mp ! microphysics non-convective precipitation [mm] - real(kind=kind_phys), dimension(:) , intent(in) :: rainc_mp ! microphysics convective precipitation [mm] - real(kind=kind_phys), dimension(:) , intent(in) :: snow_mp ! microphysics snow [mm] - real(kind=kind_phys), dimension(:) , intent(in) :: graupel_mp ! microphysics graupel [mm] - real(kind=kind_phys), dimension(:) , intent(in) :: ice_mp ! microphysics ice/hail [mm] + real(kind=kind_phys), dimension(:) , intent(in), optional :: rainn_mp ! microphysics non-convective precipitation [mm] + real(kind=kind_phys), dimension(:) , intent(in), optional :: rainc_mp ! microphysics convective precipitation [mm] + real(kind=kind_phys), dimension(:) , intent(in), optional :: snow_mp ! microphysics snow [mm] + real(kind=kind_phys), dimension(:) , intent(in), optional :: graupel_mp ! microphysics graupel [mm] + real(kind=kind_phys), dimension(:) , intent(in), optional :: ice_mp ! microphysics ice/hail [mm] real(kind=kind_phys), dimension(:) , intent(in) :: rhonewsn1 ! precipitation ice density (kg/m^3) real(kind=kind_phys) , intent(in) :: con_hvap ! latent heat condensation [J/kg] real(kind=kind_phys) , intent(in) :: con_cp ! specific heat air [J/kg/K] @@ -334,40 +334,40 @@ subroutine noahmpdrv_run & real(kind=kind_phys), dimension(:) , intent(inout) :: fm101 ! MOS function for momentum evaulated @ 10 m real(kind=kind_phys), dimension(:) , intent(inout) :: fh21 ! MOS function for heat evaulated @ 2m - real(kind=kind_phys), dimension(:) , intent(inout) :: snowxy ! actual no. of snow layers - real(kind=kind_phys), dimension(:) , intent(inout) :: tvxy ! vegetation leaf temperature [K] - real(kind=kind_phys), dimension(:) , intent(inout) :: tgxy ! bulk ground surface temperature [K] - real(kind=kind_phys), dimension(:) , intent(inout) :: canicexy ! canopy-intercepted ice [mm] - real(kind=kind_phys), dimension(:) , intent(inout) :: canliqxy ! canopy-intercepted liquid water [mm] - real(kind=kind_phys), dimension(:) , intent(inout) :: eahxy ! canopy air vapor pressure [Pa] - real(kind=kind_phys), dimension(:) , intent(inout) :: tahxy ! canopy air temperature [K] - real(kind=kind_phys), dimension(:) , intent(inout) :: cmxy ! bulk momentum drag coefficient [m/s] - real(kind=kind_phys), dimension(:) , intent(inout) :: chxy ! bulk sensible heat exchange coefficient [m/s] - real(kind=kind_phys), dimension(:) , intent(inout) :: fwetxy ! wetted or snowed fraction of the canopy [-] - real(kind=kind_phys), dimension(:) , intent(inout) :: sneqvoxy ! snow mass at last time step[mm h2o] - real(kind=kind_phys), dimension(:) , intent(inout) :: alboldxy ! snow albedo at last time step [-] - real(kind=kind_phys), dimension(:) , intent(inout) :: qsnowxy ! snowfall on the ground [mm/s] - real(kind=kind_phys), dimension(:) , intent(inout) :: wslakexy ! lake water storage [mm] - real(kind=kind_phys), dimension(:) , intent(inout) :: zwtxy ! water table depth [m] - real(kind=kind_phys), dimension(:) , intent(inout) :: waxy ! water in the "aquifer" [mm] - real(kind=kind_phys), dimension(:) , intent(inout) :: wtxy ! groundwater storage [mm] - real(kind=kind_phys), dimension(:,lsnowl:), intent(inout) :: tsnoxy ! snow temperature [K] - real(kind=kind_phys), dimension(:,lsnowl:), intent(inout) :: zsnsoxy ! snow/soil layer depth [m] - real(kind=kind_phys), dimension(:,lsnowl:), intent(inout) :: snicexy ! snow layer ice [mm] - real(kind=kind_phys), dimension(:,lsnowl:), intent(inout) :: snliqxy ! snow layer liquid water [mm] - real(kind=kind_phys), dimension(:) , intent(inout) :: lfmassxy ! leaf mass [g/m2] - real(kind=kind_phys), dimension(:) , intent(inout) :: rtmassxy ! mass of fine roots [g/m2] - real(kind=kind_phys), dimension(:) , intent(inout) :: stmassxy ! stem mass [g/m2] - real(kind=kind_phys), dimension(:) , intent(inout) :: woodxy ! mass of wood (incl. woody roots) [g/m2] - real(kind=kind_phys), dimension(:) , intent(inout) :: stblcpxy ! stable carbon in deep soil [g/m2] - real(kind=kind_phys), dimension(:) , intent(inout) :: fastcpxy ! short-lived carbon, shallow soil [g/m2] - real(kind=kind_phys), dimension(:) , intent(inout) :: xlaixy ! leaf area index [m2/m2] - real(kind=kind_phys), dimension(:) , intent(inout) :: xsaixy ! stem area index [m2/m2] - real(kind=kind_phys), dimension(:) , intent(inout) :: taussxy ! snow age factor [-] - real(kind=kind_phys), dimension(:,:) , intent(inout) :: smoiseq ! eq volumetric soil moisture [m3/m3] - real(kind=kind_phys), dimension(:) , intent(inout) :: smcwtdxy ! soil moisture content in the layer to the water table when deep - real(kind=kind_phys), dimension(:) , intent(inout) :: deeprechxy ! recharge to the water table when deep - real(kind=kind_phys), dimension(:) , intent(inout) :: rechxy ! recharge to the water table + real(kind=kind_phys), dimension(:) , intent(inout), optional :: snowxy ! actual no. of snow layers + real(kind=kind_phys), dimension(:) , intent(inout), optional :: tvxy ! vegetation leaf temperature [K] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: tgxy ! bulk ground surface temperature [K] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: canicexy ! canopy-intercepted ice [mm] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: canliqxy ! canopy-intercepted liquid water [mm] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: eahxy ! canopy air vapor pressure [Pa] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: tahxy ! canopy air temperature [K] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: cmxy ! bulk momentum drag coefficient [m/s] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: chxy ! bulk sensible heat exchange coefficient [m/s] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: fwetxy ! wetted or snowed fraction of the canopy [-] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: sneqvoxy ! snow mass at last time step[mm h2o] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: alboldxy ! snow albedo at last time step [-] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: qsnowxy ! snowfall on the ground [mm/s] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: wslakexy ! lake water storage [mm] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: zwtxy ! water table depth [m] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: waxy ! water in the "aquifer" [mm] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: wtxy ! groundwater storage [mm] + real(kind=kind_phys), dimension(:,lsnowl:), intent(inout), optional :: tsnoxy ! snow temperature [K] + real(kind=kind_phys), dimension(:,lsnowl:), intent(inout), optional :: zsnsoxy ! snow/soil layer depth [m] + real(kind=kind_phys), dimension(:,lsnowl:), intent(inout), optional :: snicexy ! snow layer ice [mm] + real(kind=kind_phys), dimension(:,lsnowl:), intent(inout), optional :: snliqxy ! snow layer liquid water [mm] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: lfmassxy ! leaf mass [g/m2] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: rtmassxy ! mass of fine roots [g/m2] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: stmassxy ! stem mass [g/m2] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: woodxy ! mass of wood (incl. woody roots) [g/m2] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: stblcpxy ! stable carbon in deep soil [g/m2] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: fastcpxy ! short-lived carbon, shallow soil [g/m2] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: xlaixy ! leaf area index [m2/m2] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: xsaixy ! stem area index [m2/m2] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: taussxy ! snow age factor [-] + real(kind=kind_phys), dimension(:,:) , intent(inout), optional :: smoiseq ! eq volumetric soil moisture [m3/m3] + real(kind=kind_phys), dimension(:) , intent(inout), optional :: smcwtdxy ! soil moisture content in the layer to the water table when deep + real(kind=kind_phys), dimension(:) , intent(inout), optional :: deeprechxy ! recharge to the water table when deep + real(kind=kind_phys), dimension(:) , intent(inout), optional :: rechxy ! recharge to the water table real(kind=kind_phys), dimension(:) , intent(out) :: albdvis ! albedo - direct visible [fraction] real(kind=kind_phys), dimension(:) , intent(out) :: albdnir ! albedo - direct NIR [fraction] real(kind=kind_phys), dimension(:) , intent(out) :: albivis ! albedo - diffuse visible [fraction] @@ -395,12 +395,12 @@ subroutine noahmpdrv_run & real(kind=kind_phys), dimension(:) , intent(out) :: snohf ! snow/freezing-rain latent heat flux [W/m2] real(kind=kind_phys), dimension(:) , intent(out) :: smcwlt2 ! dry soil moisture threshold [m3/m3] real(kind=kind_phys), dimension(:) , intent(out) :: smcref2 ! soil moisture threshold [m3/m3] - real(kind=kind_phys), dimension(:) , intent(out) :: wet1 ! normalized surface soil saturated fraction - real(kind=kind_phys), dimension(:) , intent(out) :: t2mmp ! combined T2m from tiles - real(kind=kind_phys), dimension(:) , intent(out) :: q2mp ! combined q2m from tiles + real(kind=kind_phys), dimension(:) , intent(out), optional :: wet1 ! normalized surface soil saturated fraction + real(kind=kind_phys), dimension(:) , intent(out), optional :: t2mmp ! combined T2m from tiles + real(kind=kind_phys), dimension(:) , intent(out), optional :: 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) + real(kind=kind_phys), dimension(:) , intent(out), optional :: rca ! total canopy/stomatal resistance (s/m) character(len=*) , intent(out) :: errmsg integer , intent(out) :: errflg @@ -663,6 +663,7 @@ subroutine noahmpdrv_run & real (kind=kind_phys) :: precip_freeze_frac_in ! used for penman calculation real (kind=kind_phys) :: virtfac1 ! virtual factor + real (kind=kind_phys) :: tflux ! surface flux temp real (kind=kind_phys) :: tvs1 ! surface virtual temp real (kind=kind_phys) :: vptemp ! virtual potential temp @@ -944,7 +945,8 @@ subroutine noahmpdrv_run & t2mmp(i) = temperature_bare_2m q2mp(i) = spec_humidity_bare_2m - tskin(i) = temperature_ground + tskin(i) = temperature_radiative + tflux = temperature_ground surface_temperature = temperature_ground vegetation_fraction = vegetation_frac ch_vegetated = 0.0 @@ -1038,7 +1040,8 @@ subroutine noahmpdrv_run & q2mp(i) = spec_humidity_veg_2m * vegetation_fraction + & spec_humidity_bare_2m * (1-vegetation_fraction) - tskin(i) = surface_temperature + tskin(i) = temperature_radiative + tflux = surface_temperature endif ! glacial split ends @@ -1194,9 +1197,9 @@ subroutine noahmpdrv_run & endif if(thsfc_loc) then ! Use local potential temperature - tvs1 = tskin(i) * virtfac1 + tvs1 = tflux * virtfac1 else ! Use potential temperature referenced to 1000 hPa - tvs1 = tskin(i)/prsik1(i) * virtfac1 + tvs1 = tflux/prsik1(i) * virtfac1 endif z0_total = max(min(z0_total,forcing_height),1.0e-6) diff --git a/drivers/ccpp/noahmpdrv.meta b/drivers/ccpp/noahmpdrv.meta index 7fa13f3..7535500 100644 --- a/drivers/ccpp/noahmpdrv.meta +++ b/drivers/ccpp/noahmpdrv.meta @@ -516,6 +516,7 @@ type = real kind = kind_phys intent = in + optional = True [rainc_mp] standard_name = convective_precipitation_rate_on_previous_timestep long_name = convective precipitation rate from previous timestep @@ -524,6 +525,7 @@ type = real kind = kind_phys intent = in + optional = True [snow_mp] standard_name = snowfall_rate_on_previous_timestep long_name = snow precipitation rate from previous timestep @@ -532,6 +534,7 @@ type = real kind = kind_phys intent = in + optional = True [graupel_mp] standard_name = graupel_precipitation_rate_on_previous_timestep long_name = graupel precipitation rate from previous timestep @@ -540,6 +543,7 @@ type = real kind = kind_phys intent = in + optional = True [ice_mp] standard_name = ice_precipitation_rate_on_previous_timestep long_name = ice precipitation rate from previous timestep @@ -548,6 +552,7 @@ type = real kind = kind_phys intent = in + optional = True [rhonewsn1] standard_name = surface_frozen_precipitation_density long_name = density of precipitation ice @@ -840,6 +845,7 @@ type = real kind = kind_phys intent = inout + optional = True [tvxy] standard_name = canopy_temperature long_name = vegetation temperature @@ -848,6 +854,7 @@ type = real kind = kind_phys intent = inout + optional = True [tgxy] standard_name = ground_temperature long_name = ground temperature for noahmp @@ -856,6 +863,7 @@ type = real kind = kind_phys intent = inout + optional = True [canicexy] standard_name = canopy_intercepted_ice_mass long_name = canopy intercepted ice mass @@ -864,6 +872,7 @@ type = real kind = kind_phys intent = inout + optional = True [canliqxy] standard_name = canopy_intercepted_liquid_water long_name = canopy intercepted liquid water @@ -872,6 +881,7 @@ type = real kind = kind_phys intent = inout + optional = True [eahxy] standard_name = air_vapor_pressure_in_canopy long_name = canopy air vapor pressure @@ -880,6 +890,7 @@ type = real kind = kind_phys intent = inout + optional = True [tahxy] standard_name = air_temperature_in_canopy long_name = canopy air temperature @@ -888,6 +899,7 @@ type = real kind = kind_phys intent = inout + optional = True [cmxy] standard_name = surface_drag_coefficient_for_momentum_for_noahmp long_name = surface drag coefficient for momentum for noahmp @@ -896,6 +908,7 @@ type = real kind = kind_phys intent = inout + optional = True [chxy] standard_name = surface_drag_coefficient_for_heat_and_moisture_for_noahmp long_name = surface exchange coeff heat & moisture for noahmp @@ -904,6 +917,7 @@ type = real kind = kind_phys intent = inout + optional = True [fwetxy] standard_name = wet_canopy_area_fraction long_name = area fraction of canopy that is wetted/snowed @@ -912,6 +926,7 @@ type = real kind = kind_phys intent = inout + optional = True [sneqvoxy] standard_name = lwe_thickness_of_snowfall_amount_on_previous_timestep long_name = snow mass at previous time step @@ -920,6 +935,7 @@ type = real kind = kind_phys intent = inout + optional = True [alboldxy] standard_name = surface_albedo_assuming_deep_snow_on_previous_timestep long_name = snow albedo at previous time step @@ -928,6 +944,7 @@ type = real kind = kind_phys intent = inout + optional = True [qsnowxy] standard_name = lwe_snowfall_rate long_name = snow precipitation rate at surface @@ -936,6 +953,7 @@ type = real kind = kind_phys intent = inout + optional = True [wslakexy] standard_name = water_storage_in_lake long_name = lake water storage @@ -944,6 +962,7 @@ type = real kind = kind_phys intent = inout + optional = True [zwtxy] standard_name = water_table_depth long_name = water table depth @@ -952,6 +971,7 @@ type = real kind = kind_phys intent = inout + optional = True [waxy] standard_name = water_storage_in_aquifer long_name = water storage in aquifer @@ -960,6 +980,7 @@ type = real kind = kind_phys intent = inout + optional = True [wtxy] standard_name = water_storage_in_aquifer_and_saturated_soil long_name = water storage in aquifer and saturated soil @@ -968,6 +989,7 @@ type = real kind = kind_phys intent = inout + optional = True [tsnoxy] standard_name = temperature_in_surface_snow long_name = temperature_in_surface_snow @@ -976,6 +998,7 @@ type = real kind = kind_phys intent = inout + optional = True [zsnsoxy] standard_name = depth_from_snow_surface_at_bottom_interface long_name = depth from the top of the snow surface at the bottom of the layer @@ -984,6 +1007,7 @@ type = real kind = kind_phys intent = inout + optional = True [snicexy] standard_name = lwe_thickness_of_ice_in_surface_snow long_name = lwe_thickness_of_ice_in_surface_snow @@ -992,6 +1016,7 @@ type = real kind = kind_phys intent = inout + optional = True [snliqxy] standard_name = lwe_thickness_of_liquid_water_in_surface_snow long_name = snow layer liquid water @@ -1000,6 +1025,7 @@ type = real kind = kind_phys intent = inout + optional = True [lfmassxy] standard_name = leaf_mass_content long_name = leaf mass @@ -1008,6 +1034,7 @@ type = real kind = kind_phys intent = inout + optional = True [rtmassxy] standard_name = fine_root_mass_content long_name = fine root mass @@ -1016,6 +1043,7 @@ type = real kind = kind_phys intent = inout + optional = True [stmassxy] standard_name = stem_mass_content long_name = stem mass @@ -1024,6 +1052,7 @@ type = real kind = kind_phys intent = inout + optional = True [woodxy] standard_name = wood_mass_content long_name = wood mass including woody roots @@ -1032,6 +1061,7 @@ type = real kind = kind_phys intent = inout + optional = True [stblcpxy] standard_name = slow_soil_pool_mass_content_of_carbon long_name = stable carbon in deep soil @@ -1040,6 +1070,7 @@ type = real kind = kind_phys intent = inout + optional = True [fastcpxy] standard_name = fast_soil_pool_mass_content_of_carbon long_name = short-lived carbon in shallow soil @@ -1048,6 +1079,7 @@ type = real kind = kind_phys intent = inout + optional = True [xlaixy] standard_name = leaf_area_index long_name = leaf area index @@ -1056,6 +1088,7 @@ type = real kind = kind_phys intent = inout + optional = True [xsaixy] standard_name = stem_area_index long_name = stem area index @@ -1064,6 +1097,7 @@ type = real kind = kind_phys intent = inout + optional = True [taussxy] standard_name = dimensionless_age_of_surface_snow long_name = non-dimensional snow age @@ -1072,6 +1106,7 @@ type = real kind = kind_phys intent = inout + optional = True [smoiseq] standard_name = volumetric_equilibrium_soil_moisture long_name = equilibrium soil water content @@ -1080,6 +1115,7 @@ type = real kind = kind_phys intent = inout + optional = True [smcwtdxy] standard_name = volumetric_soil_moisture_between_soil_bottom_and_water_table long_name = soil water content between the bottom of the soil and the water table @@ -1088,6 +1124,7 @@ type = real kind = kind_phys intent = inout + optional = True [deeprechxy] standard_name = water_table_recharge_assuming_deep long_name = recharge to or from the water table when deep @@ -1096,6 +1133,7 @@ type = real kind = kind_phys intent = inout + optional = True [rechxy] standard_name = water_table_recharge_assuming_shallow long_name = recharge to or from the water table when shallow @@ -1104,6 +1142,7 @@ type = real kind = kind_phys intent = inout + optional = True [albdvis] standard_name = surface_albedo_direct_visible_over_land long_name = direct surface albedo visible band over land @@ -1328,6 +1367,7 @@ type = real kind = kind_phys intent = out + optional = True [t2mmp] standard_name = temperature_at_2m_from_noahmp long_name = 2 meter temperature from noahmp @@ -1336,6 +1376,7 @@ type = real kind = kind_phys intent = out + optional = True [q2mp] standard_name = specific_humidity_at_2m_from_noahmp long_name = 2 meter specific humidity from noahmp @@ -1344,6 +1385,7 @@ type = real kind = kind_phys intent = out + optional = True [zvfun] standard_name = function_of_surface_roughness_length_and_green_vegetation_fraction long_name = function of surface roughness length and green vegetation fraction @@ -1368,6 +1410,7 @@ type = real kind = kind_phys intent = out + optional = True [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/src/module_sf_noahmplsm.F90 b/src/module_sf_noahmplsm.F90 index 2472f11..a76a354 100644 --- a/src/module_sf_noahmplsm.F90 +++ b/src/module_sf_noahmplsm.F90 @@ -1989,7 +1989,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in real (kind=kind_phys), parameter :: mpe = 1.e-6 real (kind=kind_phys), parameter :: psiwlt = -150. !metric potential for wilting point (m) - real (kind=kind_phys), parameter :: z0 = 0.002 ! bare-soil roughness length (m) (i.e., under the canopy) + real (kind=kind_phys), parameter :: z0 = 0.015 ! bare-soil roughness length (m) (i.e., under the canopy) ! --------------------------------------------------------------------------------------------------- ! initialize fluxes from veg. fraction @@ -2629,10 +2629,10 @@ subroutine csnow (parameters,isnow ,nsnow ,nsoil ,snice ,snliq ,dzsnso ! thermal conductivity of snow do iz = isnow+1, 0 -! tksno(iz) = 3.2217e-6*bdsnoi(iz)**2. ! stieglitz(yen,1965) +! tksno(iz) = 3.2217e-6*bdsnoi(iz)**2. ! stieglitz(yen,1965) ! tksno(iz) = 2e-2+2.5e-6*bdsnoi(iz)*bdsnoi(iz) ! anderson, 1976 -! tksno(iz) = 0.35 ! constant - tksno(iz) = 2.576e-6*bdsnoi(iz)**2. + 0.074 ! verseghy (1991) + tksno(iz) = 0.35 ! constant +! tksno(iz) = 2.576e-6*bdsnoi(iz)**2. + 0.074 ! verseghy (1991) ! tksno(iz) = 2.22*(bdsnoi(iz)/1000.)**1.88 ! douvill(yen, 1981) enddo @@ -4056,11 +4056,6 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & end if -! prepare for longwave rad. - - air = -emv*(1.+(1.-emv)*(1.-emg))*lwdn - emv*emg*sb*tg**4 - cir = (2.-emv*(1.-emg))*emv*sb -! if(opt_sfc == 4) then gdx = sqrt(garea1) @@ -4207,6 +4202,11 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & end if end if +! prepare for longwave rad. + + air = -emv*(1.+(1.-emv)*(1.-emg))*lwdn - emv*emg*sb*tg**4 + cir = (2.-emv*(1.-emg))*emv*sb + ! prepare for sensible heat flux above veg. cah = 1./rahc @@ -4269,7 +4269,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & ! update vegetation surface temperature tv = tv + dtv -! tah = ata + bta*tv ! canopy air t; update here for consistency + tah = ata + bta*tv ! canopy air t; update here for consistency ! for computing m-o length in the next iteration h = rhoair*cpair*(tah - sfctmp) /rahc @@ -4282,15 +4282,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & qfx = (qsfc-qair)*rhoair*caw endif - - if (liter == 1) then - exit loop1 - endif - if (iter >= 5 .and. abs(dtv) <= 0.01 .and. liter == 0) then - liter = 1 - endif - - end do loop1 ! end stability iteration +! after canopy balance, do the under-canopy ground balance ! under-canopy fluxes and tg @@ -4300,8 +4292,6 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & cev = rhoair*cpair / (gammag*(rawg+rsurf)) ! barlage: change to ground v3.6 cgh = 2.*df(isnow+1)/dzsnso(isnow+1) - loop2: do iter = 1, niterg - t = tdc(tg) call esat(t, esatw, esati, dsatw, dsati) if (t .gt. 0.) then @@ -4327,7 +4317,14 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & gh = gh + cgh*dtg tg = tg + dtg - end do loop2 + if (liter == 1) then + exit loop1 + endif + if (iter >= 5 .and. abs(dtv) <= 0.01 .and. abs(dtg) <= 0.01 .and. liter == 0) then + liter = 1 ! if conditions are met, then do one final loop + endif + + end do loop1 ! tah = (cah*sfctmp + cvh*tv + cgh*tg)/(cah + cvh + cgh) @@ -5824,7 +5821,8 @@ subroutine thermalz0(parameters, fveg, z0m, z0mg, zlvl, if (opt_trs == z0heqz0m) then - z0m_out = exp(fveg * log(z0m) + (1.0 - fveg) * log(z0mg)) +! z0m_out = exp(fveg * log(z0m) + (1.0 - fveg) * log(z0mg)) + z0m_out = fveg * z0m + (1.0 - fveg) * z0mg z0h_out = z0m_out elseif (opt_trs == chen09) then @@ -5841,7 +5839,7 @@ subroutine thermalz0(parameters, fveg, z0m, z0mg, zlvl, endif z0h_out = exp( fveg * log(z0m * exp(-czil*0.4*258.2*sqrt(ustarx*z0m))) + & - (1.0 - fveg) * log(max(z0m/exp(kb_sigma_f0),1.0e-6)) ) + (1.0 - fveg) * log(max(z0mg/exp(kb_sigma_f0),1.0e-6)) ) elseif (opt_trs == tessel) then @@ -5880,7 +5878,7 @@ subroutine thermalz0(parameters, fveg, z0m, z0mg, zlvl, z0h_out = z0m_out - elseif (opt_trs == chen09 .or. opt_trs == tessel) then + elseif (opt_trs == tessel) then if (vegtyp <= 5) then z0h_out = z0m_out @@ -5888,7 +5886,7 @@ subroutine thermalz0(parameters, fveg, z0m, z0mg, zlvl, z0h_out = z0m_out * 0.01 endif - elseif (opt_trs == blumel99) then + elseif (opt_trs == chen09 .or. opt_trs == blumel99) then reyn = ustarx*z0m_out/viscosity ! Blumel99 eqn 36c if (reyn > 2.0) then