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.

Sign up for GitHub

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

Fixes and updates related to Noah-MP #52

Merged
merged 10 commits into from
Jul 1, 2024
Merged

Fixes and updates related to Noah-MP #52

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
5c02a7b
Fix issue with iems = 2
kaiyuan-cheng Jun 27, 2024
7748d51
Avoid using ice as hail for Noah-MP
kaiyuan-cheng Jun 27, 2024
4faccbd
Correct conversion of mixing ratio
kaiyuan-cheng Jun 24, 2024
6897c12
fix potential temperature
kaiyuan-cheng Jun 26, 2024
d7f792d
fix water table recharge
kaiyuan-cheng Jun 26, 2024
18cb8d4
Update snow initialization
kaiyuan-cheng Jun 25, 2024
26ed91a
fix wilting point
kaiyuan-cheng Jun 27, 2024
7e91421
clean up
kaiyuan-cheng Jun 27, 2024
2b5fb08
Clean up
kaiyuan-cheng Jun 25, 2024
caeff54
clean up
kaiyuan-cheng Jun 24, 2024
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
16 changes: 14 additions & 2 deletions FV3GFS/FV3GFS_io.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -877,6 +877,7 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, enforce_
real(kind=kind_phys) :: masslai, masssai,snd
real(kind=kind_phys) :: ddz,expon,aa,bb,smc,func,dfunc,dx
real(kind=kind_phys) :: bexp, smcmax, smcwlt,dwsat,dksat,psisat
real(kind=kind_phys) :: emg, emv

real(kind=kind_phys), dimension(-2:0) :: dzsno
real(kind=kind_phys), dimension(-2:4) :: dzsnso
Expand Down Expand Up @@ -1351,7 +1352,6 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, enforce_
Sfcprop(nb)%albdnir(ix) = 0.2
Sfcprop(nb)%albivis(ix) = 0.2
Sfcprop(nb)%albinir(ix) = 0.2
Sfcprop(nb)%emiss(ix) = 0.95

Sfcprop(nb)%waxy(ix) = 4900.0
Sfcprop(nb)%wtxy(ix) = Sfcprop(nb)%waxy(ix)
Expand Down Expand Up @@ -1394,6 +1394,10 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, enforce_

endif ! non urban ...

emv = 1. - exp(-(Sfcprop(nb)%xlaixy(ix)+Sfcprop(nb)%xsaixy(ix))/1.0) ! Ignore snow-buried sai and lai during the initialization
emg = 0.97*(1.-Sfcprop(nb)%sncovr(ix)) + 1.0*Sfcprop(nb)%sncovr(ix)
Sfcprop(nb)%emiss(ix) = Sfcprop(nb)%vfrac(ix) * ( emg*(1-emv) + emv + emv*(1-emv)*(1-emg) ) + (1-Sfcprop(nb)%vfrac(ix)) * emg

if ( vegtyp == isice_table ) then
do lsoil = 1,Model%lsoil
Sfcprop(nb)%stc(ix,lsoil) = min(Sfcprop(nb)%stc(ix,lsoil),min(Sfcprop(nb)%tg3(ix),263.15))
Expand All @@ -1405,7 +1409,15 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, enforce_
snd = Sfcprop(nb)%snowd(ix)/1000.0 ! go to m from snwdph

if (Sfcprop(nb)%weasd(ix) /= 0.0 .and. snd == 0.0 ) then
snd = Sfcprop(nb)%weasd(ix)/1000.0
snd = Sfcprop(nb)%weasd(ix)*0.005
Sfcprop(nb)%snowd(ix) = snd*1000.0
endif

! cap SNOW at 2000, maintain density
if (Sfcprop(nb)%weasd(ix) > 2000.0 ) then
snd = snd * 2000.0 / Sfcprop(nb)%weasd(ix)
Sfcprop(nb)%weasd(ix) = 2000.0
Sfcprop(nb)%snowd(ix) = snd*1000.0
endif

if (vegtyp == 15) then ! land ice in MODIS/IGBP
Expand Down
3 changes: 2 additions & 1 deletion GFS_layer/GFS_physics_driver.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1468,6 +1468,7 @@ subroutine GFS_physics_driver &
if (dry(i)) then
Sfcprop%t2m(i) = t2mmp(i)
Sfcprop%q2m(i) = q2mp(i)
Sfcprop%emiss(i) = Radtend%semis(i)
endif
enddo
endif
Expand Down Expand Up @@ -3675,7 +3676,7 @@ subroutine GFS_physics_driver &
Sfcprop%drainncprv(:) = tem * (frain * rain1(:))
Sfcprop%dsnowprv(:) = tem * Diag%snow(:)
Sfcprop%dgraupelprv(:) = tem * Diag%graupel(:)
Sfcprop%diceprv(:) = tem * Diag%ice(:)
Sfcprop%diceprv(:) = 0.0
else
Sfcprop%draincprv(:) = 0.0
Sfcprop%drainncprv(:) = 0.0
Expand Down
2 changes: 1 addition & 1 deletion gsmphys/module_sf_noahmplsm.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -784,7 +784,7 @@ subroutine atm (parameters,sfcprs ,sfctmp ,q2 ,

!jref: seems like pair should be p1000mb??
pair = sfcprs ! atm bottom level pressure (pa)
thair = sfctmp * (sfcprs/pair)**(rair/cpair)
thair = sfctmp * (100000./pair)**(rair/cpair)

qair = q2 ! in wrf, driver converts to specific humidity

Expand Down
55 changes: 15 additions & 40 deletions gsmphys/sfc_noahmp_drv.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ subroutine noahmpdrv &
& iopt_inf,iopt_rad, iopt_alb, iopt_snf,iopt_tbot,iopt_stc, &
& iopt_gla, &
& xlatin,xcoszin, iyrlen, julian,imon, &
& rainn_mp,rainc_mp,snow_mp,graupel_mp,ice_mp, &
& rainn_mp,rainc_mp,snow_mp,graupel_mp,hail_mp, &

! --- in/outs:
& weasd, snwdph, tskin, tprcp, srflag, smc, stc, slc, &
Expand All @@ -34,15 +34,14 @@ subroutine noahmpdrv &
!
!
use machine , only : kind_phys
! use date_def, only : idate
use funcphys, only : fpvs
use physcons, only : con_g, con_hvap, con_cp, con_jcal, &
& con_eps, con_epsm1, con_fvirt, con_rd,con_hfus

use module_sf_noahmplsm
use module_sf_noahmp_glacier
use noahmp_tables, only : isice_table, co2_table, o2_table, &
& isurban_table,smcref_table,smcdry_table, &
& isurban_table,smcref_table,smcwlt_table, &
& smcmax_table,co2_table,o2_table, &
& saim_table,laim_table

Expand Down Expand Up @@ -77,7 +76,6 @@ subroutine noahmpdrv &
real(kind=kind_phys), save :: zsoil(4),sldpth(4)
data zsoil / -0.1, -0.4, -1.0, -2.0 /
data sldpth /0.1, 0.3, 0.6, 1.0 /
! data dzs /0.1, 0.3, 0.6, 1.0 /

!
! --- input:
Expand All @@ -92,7 +90,7 @@ subroutine noahmpdrv &
& t1, q1, sigmaf, dlwflx, dswsfc, snet, tg3, cm, &
& ch, prsl1, prslki, wind, shdmin, shdmax, &
& snoalb, zf, &
& rainn_mp,rainc_mp,snow_mp,graupel_mp,ice_mp
& rainn_mp,rainc_mp,snow_mp,graupel_mp,hail_mp

logical, dimension(im), intent(in) :: dry

Expand Down Expand Up @@ -133,7 +131,6 @@ subroutine noahmpdrv &

integer, dimension(im) :: jsnowxy
real (kind=kind_phys),dimension(im) :: snodep
real (kind=kind_phys),dimension(im,-2:4) :: tsnsoxy

! --- output:

Expand Down Expand Up @@ -174,7 +171,7 @@ subroutine noahmpdrv &
& flx1, flx2, flx3, ffrozp, lwdn, pc, prcp, ptu, q2, &
& q2sat, solnet, rc, rcs, rct, rcq, rcsoil, rsmin, &
& runoff1, runoff2, runoff3, sfcspd, sfcprs, sfctmp, &
& sfcems, sheat, shdfac, shdmin1d, shdmax1d, smcwlt, &
& sheat, shdfac, shdmin1d, shdmax1d, smcwlt, &
& smcdry, smcref, smcmax, sneqv, snoalb1d, snowh, &
& snomlt, sncovr, soilw, soilm, ssoil, tsea, th2, &
& xlai, zlvl, swdn, tem, psfc,fdown,t2v,tbot, qmelt
Expand Down Expand Up @@ -318,7 +315,8 @@ subroutine noahmpdrv &

do i = 1, im
if (flag_iter(i) .and. flag(i)) then
q0(i) = max(q1(i), 1.e-8) !* q1=specific humidity at level 1 (kg/kg)
q0(i) = q1(i) / (1. - q1(i)) ! convert to mixing ratio (q0)
q0(i) = max(q0(i), 1.e-8)
theta1(i) = t1(i) * prslki(i) !* adiabatic temp at level 1 (k)

tv1(i) = t1(i) * (1.0 + con_fvirt*q0(i))
Expand Down Expand Up @@ -387,7 +385,6 @@ subroutine noahmpdrv &
lwdn = dlwflx(i) !..downward lw flux at sfc in w/m2
swdn = dswsfc(i) !..downward sw flux at sfc in w/m2
solnet = snet(i) !..net sw rad flx (dn-up) at sfc in w/m2
sfcems = sfcemis(i)

sfctmp = t1(i)
sfcprs = prsl1(i)
Expand All @@ -397,7 +394,6 @@ subroutine noahmpdrv &
if (prcp > 0.0) then
if (ffrozp > 0.0) then ! rain/snow flag, one condition is enough?
snowng = .true.
qsnowxy(i) = ffrozp * prcp/10.0 !still use rho water?
else
if (sfctmp <= 275.15) frzgra = .true.
endif
Expand Down Expand Up @@ -450,7 +446,6 @@ subroutine noahmpdrv &
cmc = canopy(i)/1000. ! convert from mm to m
tsea = tsurf(i) ! clu_q2m_iter

snowh = snwdph(i) * 0.001 ! convert from mm to m
sneqv = weasd(i) * 0.001 ! convert from mm to m


Expand Down Expand Up @@ -493,11 +488,11 @@ subroutine noahmpdrv &
wtx = waxy(i)

do k = -2,0
tsnsoxy(i,k) = tsnoxy(i,k)
tsnsox(k) = tsnoxy(i,k)
enddo

do k = 1,4
tsnsoxy(i,k) = stc(i,k)
tsnsox(k) = stc(i,k)
enddo

do k = -2,0
Expand All @@ -515,7 +510,6 @@ subroutine noahmpdrv &

do k = -2, km
zsnsox(k) = zsnsoxy(i,k)
tsnsox(k) = tsnsoxy(i,k)
enddo

lfmassx = lfmassxy(i)
Expand All @@ -539,8 +533,8 @@ subroutine noahmpdrv &
enddo

smcwtdx = smcwtdxy(i)
rechx = rechxy(i)
deeprechx = deeprechxy(i)
rechx = 0.
deeprechx = 0.
!--
! the optional details for precip
!--
Expand All @@ -556,7 +550,7 @@ subroutine noahmpdrv &
pshcv = 0.
psnow = snow_mp(i)
pgrpl = graupel_mp(i)
phail = ice_mp(i)
phail = hail_mp(i)
!
!-- old
!
Expand All @@ -567,10 +561,6 @@ subroutine noahmpdrv &

snowh = snwdph(i) * 0.001 ! convert from mm to m

if (swe /= 0.0 .and. snowh == 0.0) then
snowh = 10.0 * swe /1000.0
endif

chx = chxy(i) ! maybe chxy
cmx = cmxy(i)

Expand Down Expand Up @@ -633,10 +623,6 @@ subroutine noahmpdrv &
xlaix = undefined
xsaix = undefined

smcwtdx = 0.0
rechx = 0.0
deeprechx = 0.0

do k = 1,4
smoiseqx(k) = smsoil(k)
enddo
Expand All @@ -655,7 +641,6 @@ subroutine noahmpdrv &
else
ice = 0

! write(*,*)'tsnsox(1)=',tsnsox,'tgx=',tgx
call noahmp_sflx (parameters ,&
& i , 1 , lat , iyrlen , julian , cosz ,& ! in : time/space-related
& delt , dx , dz8w , nsoil , zsoil , nsnow ,& ! in : model configuration
Expand Down Expand Up @@ -749,8 +734,8 @@ subroutine noahmpdrv &

taussxy (i) = taussx

rechxy (i) = rechx
deeprechxy(i) = deeprechx
rechxy (i) = rechxy(i) + rechx
deeprechxy(i) = deeprechxy(i) + deeprechx
smcwtdxy(i) = smcwtdx
smoiseq(i,1:4) = smoiseqx(1:4)

Expand All @@ -767,12 +752,9 @@ subroutine noahmpdrv &
weasd(i) = swe
snwdph(i) = snowh * 1000.0

! write(*,*) 'swe,snowh,can'
! write (*,*) swe,snowh*1000.0,canopy(i)
!
smcmax = smcmax_table(stype)
smcref = smcref_table(stype)
smcwlt = smcdry_table(stype)
smcwlt = smcwlt_table(stype)
!
! outs
!
Expand All @@ -796,12 +778,6 @@ subroutine noahmpdrv &
trans(i) = fctr
evap(i) = eta

! write(*,*) 'vtype, stype are',vtype,stype
! write(*,*) 'fsh,gflx,eta',fsh,ssoil,eta
! write(*,*) 'esnow,runsrf,runsub',esnow,runsrf,runsub
! write(*,*) 'evbs,evcw,trans',fgev,fcev,fctr
! write(*,*) 'snowc',fsno

tsurf(i) = trad
sfcemis(i) = emissi
if(albedo .gt. 0.0) then
Expand All @@ -816,11 +792,10 @@ subroutine noahmpdrv &
& 1.0*smsoil(4))*1000.0 ! unit conversion from m to kg m-2
!
snohf (i) = qsnbot * con_hfus ! only part of it but is diagnostic
! write(*,*) 'snohf',snohf(i)


fdown = fsa + lwdn
t2v = sfctmp * (1.0 + 0.61*q2)
! ssoil = -1.0 *ssoil

call penman (sfctmp,sfcprs,chx,t2v,th2,prcp,fdown,ssoil, &
& q2,q2sat,etp,snowng,frzgra,ffrozp,dqsdt2,emissi,fsno)
Expand Down