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

Wrapper PR for: Thompson inner loop, Thompson subcycling bugfix, fix snet bug for Noah LSM #702

Merged
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
2 changes: 1 addition & 1 deletion physics/GFS_radiation_surface.F90
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,7 @@ subroutine GFS_radiation_surface_run ( &
alnwf, facsf, facwf, &
semis_lnd, semis_ice, snoalb
character(len=3) , dimension(:), intent(in) :: lndp_var_list
real(kind=kind_phys), dimension(:), intent(in) :: albdvis_lnd, albdnir_lnd, &
real(kind=kind_phys), dimension(:), intent(inout) :: albdvis_lnd, albdnir_lnd, &
albivis_lnd, albinir_lnd
real(kind=kind_phys), dimension(:), intent(in) :: albdvis_ice, albdnir_ice, &
albivis_ice, albinir_ice
Expand Down
8 changes: 4 additions & 4 deletions physics/GFS_radiation_surface.meta
Original file line number Diff line number Diff line change
Expand Up @@ -411,7 +411,7 @@
dimensions = (horizontal_loop_extent)
type = real
kind = kind_phys
intent = in
intent = inout
optional = F
[albdnir_lnd]
standard_name = surface_albedo_direct_NIR_over_land
Expand All @@ -420,7 +420,7 @@
dimensions = (horizontal_loop_extent)
type = real
kind = kind_phys
intent = in
intent = inout
optional = F
[albivis_lnd]
standard_name = surface_albedo_diffuse_visible_over_land
Expand All @@ -429,7 +429,7 @@
dimensions = (horizontal_loop_extent)
type = real
kind = kind_phys
intent = in
intent = inout
optional = F
[albinir_lnd]
standard_name = surface_albedo_diffuse_NIR_over_land
Expand All @@ -438,7 +438,7 @@
dimensions = (horizontal_loop_extent)
type = real
kind = kind_phys
intent = in
intent = inout
optional = F
[albdvis_ice]
standard_name = surface_albedo_direct_visible_over_ice
Expand Down
58 changes: 50 additions & 8 deletions physics/module_mp_thompson.F90
Original file line number Diff line number Diff line change
Expand Up @@ -969,7 +969,7 @@ END SUBROUTINE thompson_init
SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
nwfa, nifa, nwfa2d, nifa2d, &
tt, th, pii, &
p, w, dz, dt_in, &
p, w, dz, dt_in, dt_inner, &
RAINNC, RAINNCV, &
SNOWNC, SNOWNCV, &
ICENC, ICENCV, &
Expand Down Expand Up @@ -1046,7 +1046,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: &
vt_dbz_wt
LOGICAL, INTENT(IN) :: first_time_step
REAL, INTENT(IN):: dt_in
REAL, INTENT(IN):: dt_in, dt_inner
! To support subcycling: current step and maximum number of steps
INTEGER, INTENT (IN) :: istep, nsteps
LOGICAL, INTENT (IN) :: reset_dBZ
Expand Down Expand Up @@ -1107,6 +1107,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
LOGICAL, OPTIONAL, INTENT(IN) :: diagflag
INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref
logical :: melti = .false.
INTEGER :: ndt, it

! CCPP error handling
character(len=*), optional, intent( out) :: errmsg
Expand Down Expand Up @@ -1236,7 +1237,30 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
! j_end = jte
! endif

dt = dt_in
! dt = dt_in
RAINNC(:,:) = 0.0
SNOWNC(:,:) = 0.0
ICENC(:,:) = 0.0
GRAUPELNC(:,:) = 0.0
pcp_ra(:,:) = 0.0
pcp_sn(:,:) = 0.0
pcp_gr(:,:) = 0.0
pcp_ic(:,:) = 0.0
ndt = max(nint(dt_in/dt_inner),1)
dt = dt_in/ndt
if(dt_in .le. dt_inner) dt= dt_in
if(nsteps>1 .and. ndt>1) then
if (present(errmsg) .and. present(errflg)) then
write(errmsg, '(a)') 'Logic error in mp_gt_driver: inner loop cannot be used with subcycling'
errflg = 1
return
else
write(*,'(a)') 'Warning: inner loop cannot be used with subcycling, resetting ndt=1'
ndt = 1
endif
endif

do it = 1, ndt

qc_max = 0.
qr_max = 0.
Expand Down Expand Up @@ -1412,10 +1436,10 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
tten1, qvten1, qrten1, qsten1, &
qgten1, qiten1, niten1, nrten1, ncten1, qcten1)

pcp_ra(i,j) = pptrain
pcp_sn(i,j) = pptsnow
pcp_gr(i,j) = pptgraul
pcp_ic(i,j) = pptice
pcp_ra(i,j) = pcp_ra(i,j) + pptrain
pcp_sn(i,j) = pcp_sn(i,j) + pptsnow
pcp_gr(i,j) = pcp_gr(i,j) + pptgraul
pcp_ic(i,j) = pcp_ic(i,j) + pptice
RAINNCV(i,j) = pptrain + pptsnow + pptgraul + pptice
RAINNC(i,j) = RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice
IF ( PRESENT(snowncv) .AND. PRESENT(snownc) ) THEN
Expand Down Expand Up @@ -1594,9 +1618,26 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
enddo
endif assign_extended_diagnostics

if (ndt>1 .and. it==ndt) then

SR(i,j) = (pcp_sn(i,j) + pcp_gr(i,j) + pcp_ic(i,j))/(RAINNC(i,j)+1.e-12)
RAINNCV(i,j) = RAINNC(i,j)
IF ( PRESENT (snowncv) ) THEN
SNOWNCV(i,j) = SNOWNC(i,j)
ENDIF
IF ( PRESENT (icencv) ) THEN
ICENCV(i,j) = ICENC(i,j)
ENDIF
IF ( PRESENT (graupelncv) ) THEN
GRAUPELNCV(i,j) = GRAUPELNC(i,j)
ENDIF
endif

! Diagnostic calculations only for last step
! if Thompson MP is called multiple times
last_step_only: IF (istep == nsteps) THEN
last_step_only: IF ((ndt>1 .and. it==ndt) .or. &
(nsteps>1 .and. istep==nsteps) .or. &
(nsteps==1 .and. ndt==1)) THEN

!> - Call calc_refl10cm()

Expand Down Expand Up @@ -1656,6 +1697,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
! 'ni: ', ni_max, '(', imax_ni, ',', jmax_ni, ',', kmax_ni, ')', &
! 'nr: ', nr_max, '(', imax_nr, ',', jmax_nr, ',', kmax_nr, ')'
! END DEBUG - GT
enddo ! end of nt loop

! These are always allocated
!deallocate (vtsk1)
Expand Down
6 changes: 5 additions & 1 deletion physics/module_sf_noahmplsm.f90
Original file line number Diff line number Diff line change
Expand Up @@ -1904,6 +1904,10 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in
! ground snow cover fraction [niu and yang, 2007, jgr]

fsno = 0.
if(snowh <= 1.e-6 .or. sneqv <= 1.e-3) then
snowh = 0.0
sneqv = 0.0
end if
if(snowh.gt.0.) then
bdsno = sneqv / snowh
fmelt = (bdsno/100.)**parameters%mfsno
Expand Down Expand Up @@ -7070,7 +7074,7 @@ subroutine snowh2o (parameters,nsnow ,nsoil ,dt ,qsnfro ,qsnsub , & !in
end if
end if

if(snowh <= 1.e-8 .or. sneqv <= 1.e-6) then
if(snowh <= 1.e-6 .or. sneqv <= 1.e-3) then
snowh = 0.0
sneqv = 0.0
end if
Expand Down
11 changes: 6 additions & 5 deletions physics/mp_thompson.F90
Original file line number Diff line number Diff line change
Expand Up @@ -334,7 +334,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, &
spechum, qc, qr, qi, qs, qg, ni, nr, &
is_aerosol_aware, nc, nwfa, nifa, &
nwfa2d, nifa2d, &
tgrs, prsl, phii, omega, &
tgrs, prsl, phii, omega, dt_inner, &
dtp, first_time_step, istep, nsteps, &
prcp, rain, graupel, ice, snow, sr, &
refl_10cm, reset_dBZ, do_radar_ref, &
Expand Down Expand Up @@ -389,6 +389,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, &
! Radar reflectivity
real(kind_phys), intent( out) :: refl_10cm(:,:)
logical, optional, intent(in ) :: do_radar_ref
real, intent(in ) :: dt_inner
! Cloud effective radii
real(kind_phys), optional, intent( out) :: re_cloud(:,:)
real(kind_phys), optional, intent( out) :: re_ice(:,:)
Expand Down Expand Up @@ -673,7 +674,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, &
if (do_effective_radii) then
call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, &
nc=nc, nwfa=nwfa, nifa=nifa, nwfa2d=nwfa2d, nifa2d=nifa2d, &
tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, &
tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, dt_inner=dt_inner, &
rainnc=rain_mp, rainncv=delta_rain_mp, &
snownc=snow_mp, snowncv=delta_snow_mp, &
icenc=ice_mp, icencv=delta_ice_mp, &
Expand Down Expand Up @@ -713,7 +714,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, &
else
call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, &
nc=nc, nwfa=nwfa, nifa=nifa, nwfa2d=nwfa2d, nifa2d=nifa2d, &
tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, &
tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, dt_inner=dt_inner, &
rainnc=rain_mp, rainncv=delta_rain_mp, &
snownc=snow_mp, snowncv=delta_snow_mp, &
icenc=ice_mp, icencv=delta_ice_mp, &
Expand Down Expand Up @@ -753,7 +754,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, &
else
if (do_effective_radii) then
call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, &
tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, &
tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, dt_inner=dt_inner, &
rainnc=rain_mp, rainncv=delta_rain_mp, &
snownc=snow_mp, snowncv=delta_snow_mp, &
icenc=ice_mp, icencv=delta_ice_mp, &
Expand Down Expand Up @@ -792,7 +793,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, &
qcten3=qcten3)
else
call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, &
tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, &
tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, dt_inner=dt_inner, &
rainnc=rain_mp, rainncv=delta_rain_mp, &
snownc=snow_mp, snowncv=delta_snow_mp, &
icenc=ice_mp, icencv=delta_ice_mp, &
Expand Down
9 changes: 9 additions & 0 deletions physics/mp_thompson.meta
Original file line number Diff line number Diff line change
Expand Up @@ -562,6 +562,15 @@
kind = kind_phys
intent = in
optional = F
[dt_inner]
standard_name = time_step_for_inner_loop
long_name = time step for inner loop
units = s
dimensions = ()
type = real
kind = kind_phys
intent = in
optional = F
[first_time_step]
standard_name = flag_for_first_time_step
long_name = flag for first time step for time integration loop (cold/warmstart)
Expand Down
8 changes: 7 additions & 1 deletion physics/radiation_surface.f
Original file line number Diff line number Diff line change
Expand Up @@ -411,13 +411,15 @@ subroutine setalb &
real (kind=kind_phys), dimension(:), intent(in) :: &
& slmsk, snowf, zorlf, coszf, tsknf, tairf, hprif, &
& alvsf, alnsf, alvwf, alnwf, facsf, facwf, fice, tisfc, &
& lsmalbdvis, lsmalbdnir, lsmalbivis, lsmalbinir, &
& icealbdvis, icealbdnir, icealbivis, icealbinir, &
& sncovr, sncovr_ice, snoalb, albPpert ! sfc-perts, mgehne
real (kind=kind_phys), intent(in) :: pertalb ! sfc-perts, mgehne
real (kind=kind_phys), intent(in) :: min_seaice
real (kind=kind_phys), dimension(:), intent(in) :: &
& fracl, fraco, fraci
real (kind=kind_phys), dimension(:),intent(inout) :: &
& lsmalbdvis, lsmalbdnir, lsmalbivis, lsmalbinir

logical, dimension(:), intent(in) :: &
& icy

Expand Down Expand Up @@ -537,6 +539,10 @@ subroutine setalb &
alndnd = alnwf(i)*flnd + snoalb(i) * fsno
alndvb = ab2bm *flnd + snoalb(i) * fsno
alndvd = alvwf(i)*flnd + snoalb(i) * fsno
lsmalbdnir(i) = min(0.99,max(0.01,alndnb))
lsmalbinir(i) = min(0.99,max(0.01,alndnd))
lsmalbdvis(i) = min(0.99,max(0.01,alndvb))
lsmalbivis(i) = min(0.99,max(0.01,alndvd))
else
!-- fill in values for land albedo
alndnb = 0.
Expand Down
30 changes: 24 additions & 6 deletions physics/sfc_drv.f
Original file line number Diff line number Diff line change
Expand Up @@ -94,10 +94,12 @@ end subroutine lsm_noah_finalize
! call sfc_drv !
! --- inputs: !
! ( im, km, ps, t1, q1, soiltyp, vegtype, sigmaf, !
! sfcemis, dlwflx, dswsfc, snet, delt, tg3, cm, ch, !
! sfcemis, dlwflx, dswsfc, delt, tg3, cm, ch, !
! prsl1, prslki, zf, land, wind, slopetyp, !
! shdmin, shdmax, snoalb, sfalb, flag_iter, flag_guess, !
! lheatstrg, isot, ivegsrc, !
! albdvis_lnd, albdnir_lnd, albivis_lnd, albinir_lnd, !
! adjvisbmd, adjnirbmd, adjvisdfd, adjnirdfd, !
! --- in/outs: !
! weasd, snwdph, tskin, tprcp, srflag, smc, stc, slc, !
! canopy, trans, tsurf, zorl, !
Expand Down Expand Up @@ -133,7 +135,6 @@ end subroutine lsm_noah_finalize
! sfcemis - real, sfc lw emissivity ( fraction ) im !
! dlwflx - real, total sky sfc downward lw flux ( w/m**2 ) im !
! dswflx - real, total sky sfc downward sw flux ( w/m**2 ) im !
! snet - real, total sky sfc netsw flx into ground(w/m**2) im !
! delt - real, time interval (second) 1 !
! tg3 - real, deep soil temperature (k) im !
! cm - real, surface exchange coeff for momentum (m/s) im !
Expand All @@ -148,6 +149,14 @@ end subroutine lsm_noah_finalize
! shdmax - real, max fractnl cover of green veg (not used) im !
! snoalb - real, upper bound on max albedo over deep snow im !
! sfalb - real, mean sfc diffused sw albedo (fractional) im !
! albdvis_lnd - real, sfc vis direct sw albedo over land im !
! albdnir_lnd - real, sfc nir direct sw albedo over land im !
! albivis_lnd - real, sfc vis diffused sw albedo over land im !
! albinir_lnd - real, sfc nir diffused sw albedo over land im !
! adjvisbmd - real, t adj sfc uv+vis-beam sw dnward flux im !
! adjnirbmd - real, t adj sfc nir-beam sw dnward flux im !
! adjvisdfd - real, t adj sfc uv+vis-diff sw dnward flux im !
! adjnirdfd - real, t adj sfc nir-diff sw dnward flux im !
! flag_iter- logical, im !
! flag_guess-logical, im !
! lheatstrg- logical, flag for canopy heat storage 1 !
Expand Down Expand Up @@ -203,11 +212,13 @@ end subroutine lsm_noah_finalize
subroutine lsm_noah_run &
& ( im, km, grav, cp, hvap, rd, eps, epsm1, rvrdm1, ps, & ! --- inputs:
& t1, q1, soiltyp, vegtype, sigmaf, &
& sfcemis, dlwflx, dswsfc, snet, delt, tg3, cm, ch, &
& sfcemis, dlwflx, dswsfc, delt, tg3, cm, ch, &
& prsl1, prslki, zf, land, wind, slopetyp, &
& shdmin, shdmax, snoalb, sfalb, flag_iter, flag_guess, &
& lheatstrg, isot, ivegsrc, &
& bexppert, xlaipert, vegfpert,pertvegf, & ! sfc perts, mgehne
& albdvis_lnd, albdnir_lnd, albivis_lnd, albinir_lnd, &
& adjvisbmd, adjnirbmd, adjvisdfd, adjnirdfd, &
! --- in/outs:
& weasd, snwdph, tskin, tprcp, srflag, smc, stc, slc, &
& canopy, trans, tsurf, zorl, &
Expand Down Expand Up @@ -247,10 +258,12 @@ subroutine lsm_noah_run &
integer, dimension(:), intent(in) :: soiltyp, vegtype, slopetyp

real (kind=kind_phys), dimension(:), intent(in) :: ps, &
& t1, q1, sigmaf, sfcemis, dlwflx, dswsfc, snet, tg3, cm, &
& t1, q1, sigmaf, sfcemis, dlwflx, dswsfc, tg3, cm, &
& ch, prsl1, prslki, wind, shdmin, shdmax, &
& snoalb, sfalb, zf, &
& bexppert, xlaipert, vegfpert
& bexppert, xlaipert, vegfpert, &
& albdvis_lnd, albdnir_lnd, albivis_lnd, albinir_lnd, &
& adjvisbmd, adjnirbmd, adjvisdfd, adjnirdfd

real (kind=kind_phys), intent(in) :: delt

Expand Down Expand Up @@ -400,7 +413,12 @@ subroutine lsm_noah_run &

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
!net sw rad flx (dn-up) at sfc in w/m2
solnet = adjvisbmd(i)*(1-albdvis_lnd(i)) &
& +adjnirbmd(i)*(1-albdnir_lnd(i)) &
& +adjvisdfd(i)*(1-albivis_lnd(i)) &
& +adjnirdfd(i)*(1-albinir_lnd(i))

sfcems = sfcemis(i)

sfcprs = prsl1(i)
Expand Down
Loading