Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add 'rca' variable to noahmpdrv_run interface for CCPP/physics PR#205 #10

Merged
merged 3 commits into from
May 17, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 19 additions & 3 deletions drivers/ccpp/noahmpdrv.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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)
Expand Down
8 changes: 8 additions & 0 deletions drivers/ccpp/noahmpdrv.meta
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
54 changes: 39 additions & 15 deletions drivers/ccpp/sfc_diff.f
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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

Expand All @@ -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, &
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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),
Expand Down Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion drivers/nuopc/lnd_comp_driver.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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)

!----------------------
Expand Down
3 changes: 3 additions & 0 deletions drivers/nuopc/lnd_comp_types.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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
Expand Down
10 changes: 7 additions & 3 deletions src/module_sf_noahmplsm.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand Down
Loading