From a4058658d59533af8ed2065ab0c67d536379be65 Mon Sep 17 00:00:00 2001 From: Michael Barlage Date: Tue, 18 Apr 2023 13:46:21 -0400 Subject: [PATCH 1/5] add veg and bare qsfc to output --- physics/module_sf_noahmplsm.F90 | 9 +++++++-- physics/noahmpdrv.F90 | 16 +++++++++++++--- 2 files changed, 20 insertions(+), 5 deletions(-) diff --git a/physics/module_sf_noahmplsm.F90 b/physics/module_sf_noahmplsm.F90 index e28b5447e..db653c548 100644 --- a/physics/module_sf_noahmplsm.F90 +++ b/physics/module_sf_noahmplsm.F90 @@ -442,7 +442,8 @@ subroutine noahmp_sflx (parameters, & shg , shc , shb , evg , evb , ghv , & ! out : ghb , irg , irc , irb , tr , evc , & ! out : chleaf , chuc , chv2 , chb2 , fpice , pahv , & - pahg , pahb , pah , esnow , canhs , laisun , laisha , rb & + pahg , pahb , pah , esnow , canhs , laisun , & + laisha , rb , qsfcveg , qsfcbare & #ifdef CCPP ,errmsg, errflg) #else @@ -587,6 +588,8 @@ subroutine noahmp_sflx (parameters, & real (kind=kind_phys) , intent(out) :: rb !< leaf boundary layer resistance (s/m) real (kind=kind_phys) , intent(out) :: laisun !< sunlit leaf area index (m2/m2) real (kind=kind_phys) , intent(out) :: laisha !< shaded leaf area index (m2/m2) + real (kind=kind_phys) , intent(out) :: qsfcveg !< effective spec humid over vegetation + real (kind=kind_phys) , intent(out) :: qsfcbare !< effective spec humid over bare soil !jref:start; output real (kind=kind_phys) , intent(out) :: t2mv !< 2-m air temperature over vegetated part [k] @@ -847,7 +850,9 @@ subroutine noahmp_sflx (parameters, & emissi ,pah ,canhs, & shg,shc,shb,evg,evb,ghv,ghb,irg,irc,irb,tr,evc,chleaf,chuc,chv2,chb2 ) !out - qsfc = q1 ! + qsfcveg = eah*0.622/(sfcprs - 0.378*eah) + qsfcbare = qsfc + qsfc = q1 !jref:end #ifdef CCPP if (errflg /= 0) return diff --git a/physics/noahmpdrv.F90 b/physics/noahmpdrv.F90 index 85d3363b4..88e75637b 100644 --- a/physics/noahmpdrv.F90 +++ b/physics/noahmpdrv.F90 @@ -195,7 +195,10 @@ subroutine noahmpdrv_run & ch_bare_ground_2m_ccpp, & precip_adv_heat_veg_ccpp, & precip_adv_heat_grd_v_ccpp, & - precip_adv_heat_grd_b_ccpp ) + precip_adv_heat_grd_b_ccpp, & + spec_humid_sfc_veg_ccpp, & + spec_humid_sfc_bare_ccpp & + ) use machine , only : kind_phys use funcphys, only : fpvs @@ -436,6 +439,8 @@ subroutine noahmpdrv_run & real(kind=kind_phys), dimension(:) , intent(out), optional :: precip_adv_heat_veg_ccpp real(kind=kind_phys), dimension(:) , intent(out), optional :: precip_adv_heat_grd_v_ccpp real(kind=kind_phys), dimension(:) , intent(out), optional :: precip_adv_heat_grd_b_ccpp + real(kind=kind_phys), dimension(:) , intent(out), optional :: spec_humid_sfc_veg_ccpp + real(kind=kind_phys), dimension(:) , intent(out), optional :: spec_humid_sfc_bare_ccpp ! ! --- some new options, hard code for now @@ -613,6 +618,8 @@ subroutine noahmpdrv_run & real (kind=kind_phys) :: leaf_air_resistance ! out | leaf boundary layer resistance [s/m] 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 @@ -997,11 +1004,12 @@ subroutine noahmpdrv_run & ch_vegetated_2m ,ch_bare_ground_2m ,precip_frozen_frac , & precip_adv_heat_veg ,precip_adv_heat_grd_v ,precip_adv_heat_grd_b , & precip_adv_heat_total ,snow_sublimation ,canopy_heat_storage , & + lai_sunlit ,lai_shaded ,leaf_air_resistance , & #ifdef CCPP - lai_sunlit ,lai_shaded ,leaf_air_resistance , & + spec_humid_sfc_veg ,spec_humid_sfc_bare , & errmsg ,errflg ) #else - lai_sunlit ,lai_shaded ,leaf_air_resistance ) + spec_humid_sfc_veg ,spec_humid_sfc_bare ) #endif #ifdef CCPP @@ -1129,6 +1137,8 @@ subroutine noahmpdrv_run & if(present(precip_adv_heat_veg_ccpp )) precip_adv_heat_veg_ccpp (i) = precip_adv_heat_veg if(present(precip_adv_heat_grd_v_ccpp)) precip_adv_heat_grd_v_ccpp(i) = precip_adv_heat_grd_v if(present(precip_adv_heat_grd_b_ccpp)) precip_adv_heat_grd_b_ccpp(i) = precip_adv_heat_grd_b + if(present(spec_humid_sfc_veg_ccpp )) spec_humid_sfc_veg_ccpp (i) = spec_humid_sfc_veg + if(present(spec_humid_sfc_bare_ccpp )) spec_humid_sfc_bare_ccpp (i) = spec_humid_sfc_bare wslakexy (i) = lake_water ! not active fwetxy (i) = canopy_wet_fraction From 8f2078ccbd85a47d81d1157e03db6c320bfbf7ac Mon Sep 17 00:00:00 2001 From: Zhichang Guo Date: Wed, 19 Apr 2023 18:57:48 +0000 Subject: [PATCH 2/5] include wet leaf contribution factor --- physics/module_sf_noahmplsm.F90 | 30 ++++++++++++++++++++---------- physics/rte-rrtmgp | 2 +- 2 files changed, 21 insertions(+), 11 deletions(-) diff --git a/physics/module_sf_noahmplsm.F90 b/physics/module_sf_noahmplsm.F90 index db653c548..4fe616a13 100644 --- a/physics/module_sf_noahmplsm.F90 +++ b/physics/module_sf_noahmplsm.F90 @@ -3928,6 +3928,10 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys) :: t, tdc !kelvin to degree celsius with limit -50 to +50 +! local variables + real(kind=kind_phys) :: cvw ! 0. .and. fwet > 0.) then + if (tv > tfrz) then + wlcf = min(fwet,canliq*latheav/dt/evpot) + else + wlcf = min(fwet,canice*latheav/dt/evpot) + endif + else + wlcf= fwet + endif + cvw = wlcf*cvw ctw = (1.-fwet)*(laisune/(rb+rssun) + laishae/(rb+rssha)) cgw = 1./(rawg+rsurf) - cond = caw + cew + ctw + cgw + cond = caw + cvw + ctw + cgw aea = (eair*caw + estg*cgw) / cond - bea = (cew+ctw)/cond - cev = (1.-bea)*cew*rhoair*cpair/gammav ! barlage: change to vegetation v3.6 + bea = (cvw+ctw)/cond + cev = (1.-bea)*cvw*rhoair*cpair/gammav ! barlage: change to vegetation v3.6 ctr = (1.-bea)*ctw*rhoair*cpair/gammav ! evaluate surface fluxes with current temperature and solve for dts @@ -4195,13 +4210,8 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & irc = fveg*(air + cir*tv**4) shc = fveg*rhoair*cpair*cvh * ( tv-tah) - evc = fveg*rhoair*cpair*cew * (estv-eah) / gammav ! barlage: change to v in v3.6 + evc = fveg*rhoair*cpair*cvw * (estv-eah) / gammav ! barlage: change to v in v3.6 tr = fveg*rhoair*cpair*ctw * (estv-eah) / gammav - if (tv > tfrz) then - evc = min(canliq*latheav/dt,evc) ! barlage: add if block for canice in v3.6 - else - evc = min(canice*latheav/dt,evc) - end if ! canopy heat capacity hcv = 0.02*vaie*cwat + canliq*cwat/denh2o + canice*cice/denice !j/m2/k diff --git a/physics/rte-rrtmgp b/physics/rte-rrtmgp index 0dc54f5ec..7f01618c9 160000 --- a/physics/rte-rrtmgp +++ b/physics/rte-rrtmgp @@ -1 +1 @@ -Subproject commit 0dc54f5ecaeb1e1e342efd1e02d0bcd41737bde2 +Subproject commit 7f01618c92409658bddd3afa9acb004c608f6a0d From 5edc13128252db440dea9f31c797499bc302e85a Mon Sep 17 00:00:00 2001 From: Zhichang Guo Date: Wed, 19 Apr 2023 19:33:42 +0000 Subject: [PATCH 3/5] include wet leaf contribution factor --- physics/module_sf_noahmplsm.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/module_sf_noahmplsm.F90 b/physics/module_sf_noahmplsm.F90 index 4fe616a13..a460e8426 100644 --- a/physics/module_sf_noahmplsm.F90 +++ b/physics/module_sf_noahmplsm.F90 @@ -4194,7 +4194,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & else wlcf= fwet endif - cvw = wlcf*cvw + cvw = wlcf*cew ctw = (1.-fwet)*(laisune/(rb+rssun) + laishae/(rb+rssha)) cgw = 1./(rawg+rsurf) cond = caw + cvw + ctw + cgw From c1ce88b1ecef8879cf225be82f121b1c7d6275f1 Mon Sep 17 00:00:00 2001 From: Michael Barlage Date: Thu, 20 Apr 2023 11:40:44 -0400 Subject: [PATCH 4/5] remove cvw and wlcf --- physics/module_sf_noahmplsm.F90 | 21 ++++++++------------- physics/rte-rrtmgp | 2 +- 2 files changed, 9 insertions(+), 14 deletions(-) diff --git a/physics/module_sf_noahmplsm.F90 b/physics/module_sf_noahmplsm.F90 index a460e8426..0473b43f9 100644 --- a/physics/module_sf_noahmplsm.F90 +++ b/physics/module_sf_noahmplsm.F90 @@ -3928,10 +3928,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys) :: t, tdc !kelvin to degree celsius with limit -50 to +50 -! local variables - real(kind=kind_phys) :: cvw ! 0. .and. fwet > 0.) then if (tv > tfrz) then - wlcf = min(fwet,canliq*latheav/dt/evpot) + cew = min(fwet,canliq*latheav/dt/evpot) * vaie/rb else - wlcf = min(fwet,canice*latheav/dt/evpot) + cew = min(fwet,canice*latheav/dt/evpot) * vaie/rb endif else - wlcf= fwet + cew= fwet * vaie/rb endif - cvw = wlcf*cew ctw = (1.-fwet)*(laisune/(rb+rssun) + laishae/(rb+rssha)) cgw = 1./(rawg+rsurf) - cond = caw + cvw + ctw + cgw + cond = caw + cew + ctw + cgw aea = (eair*caw + estg*cgw) / cond - bea = (cvw+ctw)/cond - cev = (1.-bea)*cvw*rhoair*cpair/gammav ! barlage: change to vegetation v3.6 + bea = (cew+ctw)/cond + cev = (1.-bea)*cew*rhoair*cpair/gammav ! barlage: change to vegetation v3.6 ctr = (1.-bea)*ctw*rhoair*cpair/gammav ! evaluate surface fluxes with current temperature and solve for dts @@ -4210,7 +4205,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & irc = fveg*(air + cir*tv**4) shc = fveg*rhoair*cpair*cvh * ( tv-tah) - evc = fveg*rhoair*cpair*cvw * (estv-eah) / gammav ! barlage: change to v in v3.6 + evc = fveg*rhoair*cpair*cew * (estv-eah) / gammav ! barlage: change to v in v3.6 tr = fveg*rhoair*cpair*ctw * (estv-eah) / gammav ! canopy heat capacity diff --git a/physics/rte-rrtmgp b/physics/rte-rrtmgp index 7f01618c9..0dc54f5ec 160000 --- a/physics/rte-rrtmgp +++ b/physics/rte-rrtmgp @@ -1 +1 @@ -Subproject commit 7f01618c92409658bddd3afa9acb004c608f6a0d +Subproject commit 0dc54f5ecaeb1e1e342efd1e02d0bcd41737bde2 From 6b9c79a09583c9115cbfdeacd364fa6c6536ae3b Mon Sep 17 00:00:00 2001 From: Michael Barlage Date: Thu, 20 Apr 2023 11:48:52 -0400 Subject: [PATCH 5/5] add back second check --- physics/module_sf_noahmplsm.F90 | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/physics/module_sf_noahmplsm.F90 b/physics/module_sf_noahmplsm.F90 index 0473b43f9..ef6f99c44 100644 --- a/physics/module_sf_noahmplsm.F90 +++ b/physics/module_sf_noahmplsm.F90 @@ -4207,6 +4207,11 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & shc = fveg*rhoair*cpair*cvh * ( tv-tah) evc = fveg*rhoair*cpair*cew * (estv-eah) / gammav ! barlage: change to v in v3.6 tr = fveg*rhoair*cpair*ctw * (estv-eah) / gammav + if (tv > tfrz) then + evc = min(canliq*latheav/dt,evc) ! barlage: add if block for canice in v3.6 + else + evc = min(canice*latheav/dt,evc) + end if ! canopy heat capacity hcv = 0.02*vaie*cwat + canliq*cwat/denh2o + canice*cice/denice !j/m2/k