Skip to content

add semi-Lagrangian sedimentation of graupel into Thompson MP #770

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 5 commits into from
Nov 23, 2021
Merged

add semi-Lagrangian sedimentation of graupel into Thompson MP #770

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
044c30f
add semi-Lagrangian sedimentation of graupel into Thompson MP
RuiyuSun Nov 5, 2021
ef3bac3
to address reviews
RuiyuSun Nov 6, 2021
4db27ba
update the namelist variables to control semi-lagrangian sedimentation
RuiyuSun Nov 18, 2021
bd4710a
remove trailing whitespace
RuiyuSun Nov 19, 2021
35cf7eb
Merge remote-tracking branch 'ncar/main' into semi-lagrangian_sedi_gr…
RuiyuSun Nov 22, 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
153 changes: 92 additions & 61 deletions physics/module_mp_thompson.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -970,8 +970,7 @@ 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, dt_inner, &
sedi_semi, sedi_semi_update, &
sedi_semi_decfl, &
sedi_semi, decfl, &
RAINNC, RAINNCV, &
SNOWNC, SNOWNCV, &
ICENC, ICENCV, &
Expand Down Expand Up @@ -1049,7 +1048,8 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
vt_dbz_wt
LOGICAL, INTENT(IN) :: first_time_step
REAL, INTENT(IN):: dt_in, dt_inner
LOGICAL, INTENT(IN) :: sedi_semi, sedi_semi_update, sedi_semi_decfl
LOGICAL, INTENT(IN) :: sedi_semi
INTEGER, INTENT(IN) :: decfl
! To support subcycling: current step and maximum number of steps
INTEGER, INTENT (IN) :: istep, nsteps
LOGICAL, INTENT (IN) :: reset_dBZ
Expand Down Expand Up @@ -1425,7 +1425,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
#endif
rand1, rand2, rand3, &
kts, kte, dt, i, j, ext_diag, &
sedi_semi, sedi_semi_update, sedi_semi_decfl, &
sedi_semi, decfl, &
!vtsk1, txri1, txrc1, &
prw_vcdc1, prw_vcde1, &
tpri_inu1, tpri_ide1_d, tpri_ide1_s, tprs_ide1, &
Expand Down Expand Up @@ -1821,7 +1821,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
! Extended diagnostics, most arrays only
! allocated if ext_diag flag is .true.
ext_diag, &
sedi_semi, sedi_semi_update, sedi_semi_decfl, &
sedi_semi, decfl, &
!vtsk1, txri1, txrc1, &
prw_vcdc1, prw_vcde1, &
tpri_inu1, tpri_ide1_d, tpri_ide1_s, tprs_ide1, &
Expand Down Expand Up @@ -1851,7 +1851,8 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
REAL, INTENT(IN):: rand1, rand2, rand3
! Extended diagnostics, most arrays only allocated if ext_diag is true
LOGICAL, INTENT(IN) :: ext_diag
LOGICAL, INTENT(IN) :: sedi_semi, sedi_semi_update, sedi_semi_decfl
LOGICAL, INTENT(IN) :: sedi_semi
INTEGER, INTENT(IN) :: decfl
REAL, DIMENSION(:), INTENT(OUT):: &
!vtsk1, txri1, txrc1, &
prw_vcdc1, &
Expand Down Expand Up @@ -1907,13 +1908,12 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
prg_rcg, prg_ihm

DOUBLE PRECISION, PARAMETER:: zeroD0 = 0.0d0
REAL, PARAMETER :: decfl = 8.0
REAL :: dtcfl,rainsfc
REAL :: dtcfl,rainsfc,graulsfc
Comment thread
climbfuji marked this conversation as resolved.
INTEGER :: niter

REAL, DIMENSION(kts:kte):: temp, pres, qv
REAL, DIMENSION(kts:kte):: rc, ri, rr, rs, rg, ni, nr, nc, nwfa, nifa
REAL, DIMENSION(kts:kte):: rr_tmp,nr_tmp
REAL, DIMENSION(kts:kte):: rr_tmp, nr_tmp, rg_tmp
REAL, DIMENSION(kts:kte):: rho, rhof, rhof2
REAL, DIMENSION(kts:kte):: qvs, qvsi, delQvs
REAL, DIMENSION(kts:kte):: satw, sati, ssatw, ssati
Expand All @@ -1927,7 +1927,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &

REAL, DIMENSION(kts:kte):: sed_r, sed_s, sed_g, sed_i, sed_n,sed_c

REAL:: rgvm, delta_tp, orho, lfus2
REAL:: rgvm, delta_tp, orho, lfus2, orhodt
REAL, DIMENSION(5):: onstep
DOUBLE PRECISION:: N0_exp, N0_min, lam_exp, lamc, lamr, lamg
DOUBLE PRECISION:: lami, ilami, ilamc
Expand Down Expand Up @@ -3930,44 +3930,41 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
else !if(.not. sedi_semi)
niter = 1
dtcfl = dt
if(sedi_semi_decfl) then
niter = int(nstep/decfl) + 1
dtcfl = dt/niter
endif
niter = int(nstep/max(decfl,1)) + 1
dtcfl = dt/niter
do n = 1, niter
rr_tmp(:) = rr(:)
nr_tmp(:) = nr(:)
call nislfv_rain_ppm(kte,dzq,vtrk,rr,rainsfc,dtcfl,R1)
call nislfv_rain_ppm(kte,dzq,vtnrk,nr,vtr,dtcfl,R2)
call semi_lagrange_sedim(kte,dzq,vtrk,rr,rainsfc,dtcfl,R1)
call semi_lagrange_sedim(kte,dzq,vtnrk,nr,vtr,dtcfl,R2)
do k = kts, kte
qrten(k) = qrten(k) + (rr(k) - rr_tmp(k))/rho(k)/dt
nrten(k) = nrten(k) + (nr(k) - nr_tmp(k))/rho(k)/dt
orhodt = 1./(rho(k)*dt)
qrten(k) = qrten(k) + (rr(k) - rr_tmp(k)) * orhodt
nrten(k) = nrten(k) + (nr(k) - nr_tmp(k)) * orhodt
enddo
pptrain = pptrain + rainsfc

if(sedi_semi_update) then
do k = kte+1, kts, -1
vtrk(k) = 0.
vtnrk(k) = 0.
enddo
do k = kte, kts, -1
vtr = 0.
if (rr(k).gt. R1) then
lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
vtr = rhof(k)*av_r*crg(6)*org3 * lamr**cre(3) &
*((lamr+fv_r)**(-cre(6)))
vtrk(k) = vtr
do k = kte+1, kts, -1
vtrk(k) = 0.
vtnrk(k) = 0.
enddo
do k = kte, kts, -1
vtr = 0.
if (rr(k).gt. R1) then
lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
vtr = rhof(k)*av_r*crg(6)*org3 * lamr**cre(3) &
*((lamr+fv_r)**(-cre(6)))
vtrk(k) = vtr
! First below is technically correct:
! vtr = rhof(k)*av_r*crg(5)*org2 * lamr**cre(2) &
! *((lamr+fv_r)**(-cre(5)))
! Test: make number fall faster (but still slower than mass)
! Goal: less prominent size sorting
vtr = rhof(k)*av_r*crg(7)/crg(12) * lamr**cre(12) &
*((lamr+fv_r)**(-cre(7)))
vtnrk(k) = vtr
endif
enddo
endif ! if(sedi_semi_update)
vtr = rhof(k)*av_r*crg(7)/crg(12) * lamr**cre(12) &
*((lamr+fv_r)**(-cre(7)))
vtnrk(k) = vtr
endif
enddo
enddo
endif! if(.not. sedi_semi)
endif
Expand Down Expand Up @@ -4054,28 +4051,59 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &

if (ANY(L_qg .eqv. .true.)) then
nstep = NINT(1./onstep(4))
do n = 1, nstep
do k = kte, kts, -1
sed_g(k) = vtgk(k)*rg(k)
enddo
k = kte
odzq = 1./dzq(k)
orho = 1./rho(k)
qgten(k) = qgten(k) - sed_g(k)*odzq*onstep(4)*orho
rg(k) = MAX(R1, rg(k) - sed_g(k)*odzq*DT*onstep(4))
do k = ksed1(4), kts, -1
odzq = 1./dzq(k)
orho = 1./rho(k)
qgten(k) = qgten(k) + (sed_g(k+1)-sed_g(k)) &
*odzq*onstep(4)*orho
rg(k) = MAX(R1, rg(k) + (sed_g(k+1)-sed_g(k)) &
if(.not. sedi_semi) then
do n = 1, nstep
do k = kte, kts, -1
sed_g(k) = vtgk(k)*rg(k)
enddo
k = kte
odzq = 1./dzq(k)
orho = 1./rho(k)
qgten(k) = qgten(k) - sed_g(k)*odzq*onstep(4)*orho
rg(k) = MAX(R1, rg(k) - sed_g(k)*odzq*DT*onstep(4))
do k = ksed1(4), kts, -1
odzq = 1./dzq(k)
orho = 1./rho(k)
qgten(k) = qgten(k) + (sed_g(k+1)-sed_g(k)) &
*odzq*onstep(4)*orho
rg(k) = MAX(R1, rg(k) + (sed_g(k+1)-sed_g(k)) &
*odzq*DT*onstep(4))
enddo
enddo

if (rg(kts).gt.R1*10.) &
pptgraul = pptgraul + sed_g(kts)*DT*onstep(4)
enddo
endif
if (rg(kts).gt.R1*10.) &
pptgraul = pptgraul + sed_g(kts)*DT*onstep(4)
enddo
else ! if(.not. sedi_semi) then
niter = 1
dtcfl = dt
niter = int(nstep/max(decfl,1)) + 1
dtcfl = dt/niter

do n = 1, niter
rg_tmp(:) = rg(:)
call semi_lagrange_sedim(kte,dzq,vtgk,rg,graulsfc,dtcfl,R1)
do k = kts, kte
orhodt = 1./(rho(k)*dt)
qgten(k) = qgten(k) + (rg(k) - rg_tmp(k))*orhodt
enddo
pptgraul = pptgraul + graulsfc
do k = kte+1, kts, -1
vtgk(k) = 0.
enddo
do k = kte, kts, -1
vtg = 0.
if (rg(k).gt. R1) then
vtg = rhof(k)*av_g*cgg(6)*ogg3 * ilamg(k)**bv_g
if (temp(k).gt. T_0) then
vtgk(k) = MAX(vtg, vtrk(k))
else
vtgk(k) = vtg
endif
endif
enddo
enddo
endif ! if(.not. sedi_semi) then
endif

!+---+-----------------------------------------------------------------+
!> - Instantly melt any cloud ice into cloud water if above 0C and
Expand Down Expand Up @@ -6102,13 +6130,13 @@ subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, &
end subroutine calc_refl10cm
!
!-------------------------------------------------------------------
SUBROUTINE nislfv_rain_ppm(km,dzl,wwl,rql,precip,dt,R1)
SUBROUTINE semi_lagrange_sedim(km,dzl,wwl,rql,precip,dt,R1)
!-------------------------------------------------------------------
!
! for non-iteration semi-Lagrangain forward advection for cloud
! This routine is a semi-Lagrangain forward advection for hydrometeors
! with mass conservation and positive definite advection
! 2nd order interpolation with monotonic piecewise parabolic method
! this routine is under assumption of decfl < 1 for semi_Lagrangian
! 2nd order interpolation with monotonic piecewise parabolic method is used.
! This routine is under assumption of decfl < 1 for semi_Lagrangian
!
! dzl depth of model layer in meter
! wwl terminal velocity at model layer m/s
Expand All @@ -6118,6 +6146,9 @@ SUBROUTINE nislfv_rain_ppm(km,dzl,wwl,rql,precip,dt,R1)
!
! author: hann-ming henry juang <henry.juang@noaa.gov>
! implemented by song-you hong
! reference: Juang, H.-M., and S.-Y. Hong, 2010: Forward semi-Lagrangian advection
! with mass conservation and positive definiteness for falling
! hydrometeors. *Mon. Wea. Rev.*, *138*, 1778-1791
!
implicit none

Expand Down Expand Up @@ -6320,7 +6351,7 @@ SUBROUTINE nislfv_rain_ppm(km,dzl,wwl,rql,precip,dt,R1)
!
! ----------------------------------
!
END SUBROUTINE nislfv_rain_ppm
END SUBROUTINE semi_lagrange_sedim
!+---+-----------------------------------------------------------------+
!+---+-----------------------------------------------------------------+
!+---+-----------------------------------------------------------------+
Expand Down
12 changes: 4 additions & 8 deletions physics/mp_thompson.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -302,8 +302,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, &
is_aerosol_aware, nc, nwfa, nifa, &
nwfa2d, nifa2d, &
tgrs, prsl, phii, omega, &
sedi_semi, sedi_semi_update, &
sedi_semi_decfl, dtp, dt_inner, &
sedi_semi, decfl, dtp, dt_inner, &
first_time_step, istep, nsteps, &
prcp, rain, graupel, ice, snow, sr, &
refl_10cm, reset_dBZ, do_radar_ref, &
Expand Down Expand Up @@ -359,8 +358,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, &
real(kind_phys), intent(inout) :: refl_10cm(:,:)
logical, intent(in ) :: do_radar_ref
logical, intent(in) :: sedi_semi
logical, intent(in) :: sedi_semi_update
logical, intent(in) :: sedi_semi_decfl
integer, intent(in) :: decfl
! MPI and block information
integer, intent(in) :: blkno
integer, intent(in) :: mpicomm
Expand Down Expand Up @@ -618,8 +616,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, &
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=dtstep, dt_inner=dt_inner, &
sedi_semi=sedi_semi, sedi_semi_update=sedi_semi_update, &
sedi_semi_decfl=sedi_semi_decfl, &
sedi_semi=sedi_semi, decfl=decfl, &
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 @@ -658,8 +655,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, &
tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, dt_inner=dt_inner, &
sedi_semi=sedi_semi, sedi_semi_update=sedi_semi_update, &
sedi_semi_decfl=sedi_semi_decfl, &
sedi_semi=sedi_semi, decfl=decfl, &
rainnc=rain_mp, rainncv=delta_rain_mp, &
snownc=snow_mp, snowncv=delta_snow_mp, &
icenc=ice_mp, icencv=delta_ice_mp, &
Expand Down
17 changes: 5 additions & 12 deletions physics/mp_thompson.meta
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -474,19 +474,12 @@
dimensions = ()
type = logical
intent = in
[sedi_semi_update]
standard_name = flag_for_v_update_in_semi_Lagrangian_sedi
long_name = flag for v update in semi Lagrangian sedi of rain
units = flag
[decfl]
standard_name = deformed_CFL_factor
long_name = deformed CFL factor
units = count
dimensions = ()
type = logical
intent = in
[sedi_semi_decfl]
standard_name = flag_for_iteration_with_semi_Lagrangian_sedi
long_name = flag for interation with semi Lagrangian sedi of rain
units = flag
dimensions = ()
type = logical
type = integer
intent = in
[dtp]
standard_name = timestep_for_physics
Expand Down