Skip to content

adding inner loop to Thompson microphysics #685

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

Merged
climbfuji merged 6 commits into from
Jul 23, 2021
Merged

adding inner loop to Thompson microphysics #685

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
dd04dc3
adding inner loop to Thompson microphysics
RuiyuSun Jun 27, 2021
7178cb2
change the logic
RuiyuSun Jun 29, 2021
f7450c2
remove extra i
RuiyuSun Jul 1, 2021
581f5a6
Bugfix in Thompson MP, pass correct timestep to core routine
climbfuji Jun 26, 2021
b498507
Remove test code from physics/mp_thompson.F90
DomHeinzeller Jun 28, 2021
9e1a47a
Merge pull request #1 from climbfuji/innerloop_ruiyu_dom_mods
RuiyuSun Jul 1, 2021
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
97 changes: 88 additions & 9 deletions physics/module_mp_thompson.F90
View file
Open in desktop
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 @@ -1029,7 +1029,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 All @@ -1056,6 +1056,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 @@ -1141,7 +1142,25 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
! j_end = jte
! endif

dt = dt_in
!->rsun minor time interval
! 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) then
write(*,*) 'WARNING: innerloop cant not be used with sybcycle'
ndt = 1
endif
do it = 1, ndt
!<-rsun

qc_max = 0.
qr_max = 0.
Expand Down Expand Up @@ -1260,10 +1279,10 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
rand1, rand2, rand3, &
kts, kte, dt, i, j)

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 All @@ -1287,8 +1306,6 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
ENDIF
SR(i,j) = (pptsnow + pptgraul + pptice)/(RAINNCV(i,j)+1.e-12)



!..Reset lowest model level to initial state aerosols (fake sfc source).
!.. Changed 13 May 2013 to fake emissions in which nwfa2d is aerosol
!.. number tendency (number per kg per second).
Expand Down Expand Up @@ -1396,6 +1413,66 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
endif
enddo

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

!> - Call calc_refl10cm()

diagflag_presenti: IF ( PRESENT (diagflag) ) THEN

@climbfuji climbfuji Jun 28, 2021

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The code duplication in lines 1432-1472 can and should be avoided if I see this correctly by using the appropriate combination of conditions.

if (diagflag .and. do_radar_ref == 1) then
!
! Only set melti to true at the output times
if (reset_dBZ) then
melti=.true.
else
melti=.false.
endif
!
if (present(vt_dbz_wt)) then
call calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, &
t1d, p1d, dBZ, rand1, kts, kte, i, j, &
melti, vt_dbz_wt(i,:,j), &
first_time_step)
else
call calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, &
t1d, p1d, dBZ, rand1, kts, kte, i, j, &
melti)
end if
do k = kts, kte
refl_10cm(i,k,j) = MAX(-35., dBZ(k))
enddo
endif
ENDIF diagflag_presenti

IF (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) THEN
do k = kts, kte
re_qc1d(k) = re_qc_min
re_qi1d(k) = re_qi_min
re_qs1d(k) = re_qs_min
enddo
!> - Call calc_effectrad()
call calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, &
re_qc1d, re_qi1d, re_qs1d, kts, kte)
do k = kts, kte
re_cloud(i,k,j) = MAX(re_qc_min, MIN(re_qc1d(k), re_qc_max))
re_ice(i,k,j) = MAX(re_qi_min, MIN(re_qi1d(k), re_qi_max))
re_snow(i,k,j) = MAX(re_qs_min, MIN(re_qs1d(k), re_qs_max))
enddo
ENDIF

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

! Diagnostic calculations only for last step
! if Thompson MP is called multiple times
last_step_only: IF (istep == nsteps) THEN
Expand Down Expand Up @@ -1444,6 +1521,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
enddo
ENDIF
ENDIF last_step_only
endif !if (ndt > 1 .and. it == ndt) then

enddo i_loop
enddo j_loop
Expand All @@ -1458,6 +1536,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

END SUBROUTINE mp_gt_driver
!> @}
Expand Down
11 changes: 6 additions & 5 deletions physics/mp_thompson.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -320,7 +320,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 @@ -374,6 +374,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 @@ -565,7 +566,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=dtp, 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 All @@ -586,7 +587,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=dtp, 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 All @@ -607,7 +608,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=dtp, 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 All @@ -627,7 +628,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, &
first_time_step=first_time_step, errmsg=errmsg, errflg=errflg)
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=dtp, 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
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -545,6 +545,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