Skip to content

Commit

Permalink
Merge pull request #702 from climbfuji/thompson_inner_loop_subcycling…
Browse files Browse the repository at this point in the history 
…_bugfix_noah_bugfix

Wrapper PR for: Thompson inner loop, Thompson subcycling bugfix, fix snet bug for Noah LSM
  • Loading branch information
@climbfuji
climbfuji authored Jul 23, 2021
2 parents 09faa73 + cdbe374 commit 610fbcb
Show file tree
Hide file tree
Showing 9 changed files with 178 additions and 35 deletions.
2 changes: 1 addition & 1 deletion physics/GFS_radiation_surface.F90
View file
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
View file
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
View file
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
View file
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
View file
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
View file
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
View file
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
View file
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

0 comments on commit 610fbcb

Please sign in to comment.