Skip to content

Commit

Permalink
Merge pull request #820 from JeffBeck-NOAA/feature/stoch_spp
Browse files Browse the repository at this point in the history 
Add support for Stochastically Perturbed Parameterizations (SPP) in FV3
  • Loading branch information
@grantfirl
grantfirl authored Feb 23, 2022
2 parents 3b2f061 + f77322b commit ff6395c
Show file tree
Hide file tree
Showing 18 changed files with 254 additions and 100 deletions.
23 changes: 22 additions & 1 deletion physics/GFS_rrtmg_pre.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, &
gasvmr_ccl4, gasvmr_cfc113, aerodp, clouds6, clouds7, clouds8, &
clouds9, cldsa, cldfra, cldfra2d, lwp_ex,iwp_ex, lwp_fc,iwp_fc, &
faersw1, faersw2, faersw3, faerlw1, faerlw2, faerlw3, alpha, &
errmsg, errflg)
spp_wts_rad, spp_rad, errmsg, errflg)

use machine, only: kind_phys

Expand Down Expand Up @@ -104,6 +104,9 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, &
uni_cld, effr_in, do_mynnedmf, &
lmfshal, lmfdeep2, pert_clds

integer, intent(in) :: spp_rad
real(kind_phys), intent(in) :: spp_wts_rad(:,:)

real(kind=kind_phys), intent(in) :: fhswr, fhlwr, solhr, sup, julian, sppt_amp
real(kind=kind_phys), intent(in) :: con_eps, epsm1, fvirt, rog, rocp, con_rd

Expand Down Expand Up @@ -1060,6 +1063,24 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, &
enddo
enddo

if ( spp_rad == 1 ) then
do k=1,lm
if (k < levs) then
do i=1,im
clouds3(i,k) = clouds3(i,k) - spp_wts_rad(i,k) * clouds3(i,k)
clouds5(i,k) = clouds5(i,k) - spp_wts_rad(i,k) * clouds5(i,k)
clouds9(i,k) = clouds9(i,k) - spp_wts_rad(i,k) * clouds9(i,k)
enddo
else
do i=1,im
clouds3(i,k) = clouds3(i,k) - spp_wts_rad(i,k) * clouds3(i,k)
clouds5(i,k) = clouds5(i,k) - spp_wts_rad(i,k) * clouds5(i,k)
clouds9(i,k) = clouds9(i,k) - spp_wts_rad(i,k) * clouds9(i,k)
enddo
endif
enddo
endif

! mg, sfc-perts
! --- scale random patterns for surface perturbations with
! perturbation size
Expand Down
17 changes: 16 additions & 1 deletion physics/GFS_rrtmg_pre.meta
View file
Original file line number Diff line number Diff line change
Expand Up @@ -532,7 +532,7 @@
[sfc_wts]
standard_name = surface_stochastic_weights_from_coupled_process
long_name = weights for stochastic surface physics perturbation
units = none
units = 1
dimensions = (horizontal_loop_extent,number_of_perturbed_land_surface_variables)
type = real
kind = kind_phys
Expand Down Expand Up @@ -1082,6 +1082,21 @@
type = real
kind = kind_phys
intent = out
[spp_wts_rad]
standard_name = spp_weights_for_radiation_scheme
long_name = spp weights for radiation scheme
units = 1
dimensions = (horizontal_loop_extent,vertical_layer_dimension)
type = real
kind = kind_phys
intent = in
[spp_rad]
standard_name = control_for_radiation_spp_perturbations
long_name = control for radiation spp perturbations
units = count
dimensions = ()
type = integer
intent = in
[errmsg]
standard_name = ccpp_error_message
long_name = error message for error handling in CCPP
Expand Down
2 changes: 1 addition & 1 deletion physics/GFS_surface_generic.meta
View file
Original file line number Diff line number Diff line change
Expand Up @@ -310,7 +310,7 @@
[sfc_wts]
standard_name = surface_stochastic_weights_from_coupled_process
long_name = weights for stochastic surface physics perturbation
units = none
units = 1
dimensions = (horizontal_loop_extent,number_of_perturbed_land_surface_variables)
type = real
kind = kind_phys
Expand Down
47 changes: 35 additions & 12 deletions physics/drag_suite.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -218,7 +218,8 @@ subroutine drag_suite_run( &
& do_gsl_drag_ls_bl, do_gsl_drag_ss, do_gsl_drag_tofd, &
& dtend, dtidx, index_of_process_orographic_gwd, &
& index_of_temperature, index_of_x_wind, &
& index_of_y_wind, ldiag3d, errmsg, errflg)
& index_of_y_wind, ldiag3d, &
& spp_wts_gwd, spp_gwd, errmsg, errflg)

! ********************************************************************
! -----> I M P L E M E N T A T I O N V E R S I O N <----------
Expand Down Expand Up @@ -365,6 +366,11 @@ subroutine drag_suite_run( &
real(kind=kind_phys), dimension(im,km) :: zl ! = PHIL/g

!SPP
real(kind=kind_phys), dimension(im) :: var_stoch, varss_stoch, &
varmax_ss_stoch, varmax_fd_stoch
real(kind=kind_phys), intent(in) :: spp_wts_gwd(:,:)
integer, intent(in) :: spp_gwd

real(kind=kind_phys), dimension(im) :: rstoch

!Output:
Expand Down Expand Up @@ -595,6 +601,23 @@ subroutine drag_suite_run( &
endif
enddo

! SPP, if spp_gwd is 0, no perturbations are applied.
if ( spp_gwd==1 ) then
do i = its,im
var_stoch(i) = var(i) + var(i)*0.75*spp_wts_gwd(i,1)
varss_stoch(i) = varss(i) + varss(i)*0.75*spp_wts_gwd(i,1)
varmax_ss_stoch(i) = varmax_ss + varmax_ss*0.75*spp_wts_gwd(i,1)
varmax_fd_stoch(i) = varmax_fd + varmax_fd*0.75*spp_wts_gwd(i,1)
enddo
else
do i = its,im
var_stoch(i) = var(i)
varss_stoch(i) = varss(i)
varmax_ss_stoch(i) = varmax_ss
varmax_fd_stoch(i) = varmax_fd
enddo
endif

!--- calculate length of grid for flow-blocking drag
!
do i=1,im
Expand Down Expand Up @@ -711,7 +734,7 @@ subroutine drag_suite_run( &
! determine reference level: maximum of 2*var and pbl heights
!
do i = its,im
zlowtop(i) = 2. * var(i)
zlowtop(i) = 2. * var_stoch(i)
enddo
!
do i = its,im
Expand Down Expand Up @@ -867,7 +890,7 @@ subroutine drag_suite_run( &
!
ldrag(i) = ldrag(i) .or. bnv2(i,1).le.0.0
ldrag(i) = ldrag(i) .or. ulow(i).eq.1.0
ldrag(i) = ldrag(i) .or. var(i) .le. 0.0
ldrag(i) = ldrag(i) .or. var_stoch(i) .le. 0.0
!
! set all ri low level values to the low level value
!
Expand All @@ -877,7 +900,7 @@ subroutine drag_suite_run( &
!
if (.not.ldrag(i)) then
bnv(i) = sqrt( bnv2(i,1) )
fr(i) = bnv(i) * rulow(i) * 2. * var(i) * od(i)
fr(i) = bnv(i) * rulow(i) * 2. * var_stoch(i) * od(i)
fr(i) = min(fr(i),frmax)
xn(i) = ubar(i) * rulow(i)
yn(i) = vbar(i) * rulow(i)
Expand Down Expand Up @@ -961,7 +984,7 @@ subroutine drag_suite_run( &
exit
ENDIF
enddo
if((xland(i)-1.5).le.0. .and. 2.*varss(i).le.hpbl(i))then
if((xland(i)-1.5).le.0. .and. 2.*varss_stoch(i).le.hpbl(i))then
if(br1(i).gt.0. .and. thvx(i,kpbl2)-thvx(i,kts) > 0.)then
cleff_ss = sqrt(dxy(i)**2 + dxyp(i)**2) ! WRF
! cleff_ss = 3. * max(dx(i),cleff_ss)
Expand All @@ -980,8 +1003,8 @@ subroutine drag_suite_run( &
!tauwavex0=0.5*XNBV*xlinv(i)*(2*MIN(varss(i),75.))**2*ro(i,kts)*u1(i,kpbl(i))
!tauwavex0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*u1(i,kpbl2)
!tauwavex0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*u1(i,3)
var_temp = MIN(varss(i),varmax_ss) + &
MAX(0.,beta_ss*(varss(i)-varmax_ss))
var_temp = MIN(varss_stoch(i),varmax_ss_stoch(i)) + &
MAX(0.,beta_ss*(varss_stoch(i)-varmax_ss_stoch(i)))
! Note: This is a semi-implicit treatment of the time differencing
var_temp2 = 0.5*XNBV*xlinv(i)*(2.*var_temp)**2*ro(i,kvar) ! this is greater than zero
tauwavex0=var_temp2*u1(i,kvar)/(1.+var_temp2*deltim)
Expand All @@ -995,8 +1018,8 @@ subroutine drag_suite_run( &
!tauwavey0=0.5*XNBV*xlinv(i)*(2*MIN(varss(i),75.))**2*ro(i,kts)*v1(i,kpbl(i))
!tauwavey0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*v1(i,kpbl2)
!tauwavey0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*v1(i,3)
var_temp = MIN(varss(i),varmax_ss) + &
MAX(0.,beta_ss*(varss(i)-varmax_ss))
var_temp = MIN(varss_stoch(i),varmax_ss_stoch(i)) + &
MAX(0.,beta_ss*(varss_stoch(i)-varmax_ss_stoch(i)))
! Note: This is a semi-implicit treatment of the time differencing
var_temp2 = 0.5*XNBV*xlinv(i)*(2.*var_temp)**2*ro(i,kvar) ! this is greater than zero
tauwavey0=var_temp2*v1(i,kvar)/(1.+var_temp2*deltim)
Expand Down Expand Up @@ -1060,16 +1083,16 @@ subroutine drag_suite_run( &

IF ((xland(i)-1.5) .le. 0.) then
!(IH*kflt**n1)**-1 = (0.00102*0.00035**-1.9)**-1 = 0.00026615161
var_temp = MIN(varss(i),varmax_fd) + &
MAX(0.,beta_fd*(varss(i)-varmax_fd))
var_temp = MIN(varss_stoch(i),varmax_fd_stoch(i)) + &
MAX(0.,beta_fd*(varss_stoch(i)-varmax_fd_stoch(i)))
var_temp = MIN(var_temp, 250.)
a1=0.00026615161*var_temp**2
! a1=0.00026615161*MIN(varss(i),varmax)**2
! a1=0.00026615161*(0.5*varss(i))**2
! k1**(n1-n2) = 0.003**(-1.9 - -2.8) = 0.003**0.9 = 0.005363
a2=a1*0.005363
! Revise e-folding height based on PBL height and topographic std. dev. -- M. Toy 3/12/2018
H_efold = max(2*varss(i),hpbl(i))
H_efold = max(2*varss_stoch(i),hpbl(i))
H_efold = min(H_efold,1500.)
DO k=kts,km
wsp=SQRT(u1(i,k)**2 + v1(i,k)**2)
Expand Down
15 changes: 15 additions & 0 deletions physics/drag_suite.meta
View file
Original file line number Diff line number Diff line change
Expand Up @@ -624,6 +624,21 @@
dimensions = ()
type = logical
intent = in
[spp_wts_gwd]
standard_name = spp_weights_for_gravity_wave_drag_scheme
long_name = spp weights for gravity wave drag scheme
units = 1
dimensions = (horizontal_loop_extent,vertical_layer_dimension)
type = real
kind = kind_phys
intent = in
[spp_gwd]
standard_name = control_for_gravity_wave_drag_spp_perturbations
long_name = control for gravity wave drag spp perturbations
units = count
dimensions = ()
type = integer
intent = in
[errmsg]
standard_name = ccpp_error_message
long_name = error message for error handling in CCPP
Expand Down
21 changes: 12 additions & 9 deletions physics/module_MYNNPBL_wrapper.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -108,7 +108,7 @@ SUBROUTINE mynnedmf_wrapper_run( &
& icloud_bl, do_mynnsfclay, &
& imp_physics, imp_physics_gfdl, &
& imp_physics_thompson, imp_physics_wsm6, &
& ltaerosol, lprnt, huge, errmsg, errflg )
& ltaerosol, spp_wts_pbl, spp_pbl, lprnt, huge, errmsg, errflg )

! should be moved to inside the mynn:
use machine , only : kind_phys
Expand Down Expand Up @@ -209,7 +209,8 @@ SUBROUTINE mynnedmf_wrapper_run( &
& bl_mynn_output, &
& grav_settling, &
& imp_physics, imp_physics_wsm6, &
& imp_physics_thompson, imp_physics_gfdl
& imp_physics_thompson, imp_physics_gfdl, &
& spp_pbl

!TENDENCY DIAGNOSTICS
real(kind=kind_phys), intent(inout), optional :: dtend(:,:,:)
Expand All @@ -220,8 +221,8 @@ SUBROUTINE mynnedmf_wrapper_run( &

!MISC CONFIGURATION OPTIONS
INTEGER, PARAMETER :: &
& spp_pbl=0, &
& bl_mynn_mixscalars=1
& bl_mynn_mixscalars=1, &
& levflag=2
REAL, PARAMETER :: &
& closure=2.6 !2.5, 2.6 or 3.0
LOGICAL :: &
Expand All @@ -231,7 +232,7 @@ SUBROUTINE mynnedmf_wrapper_run( &
LOGICAL, PARAMETER :: cycling = .false.
INTEGER, PARAMETER :: param_first_scalar = 1
INTEGER :: &
& p_qc, p_qr, p_qi, p_qs, p_qg, p_qnc, p_qni
& p_qc, p_qr, p_qi, p_qs, p_qg, p_qnc, p_qni

!MYNN-1D
REAL(kind=kind_phys), intent(in) :: delt, dtf
Expand Down Expand Up @@ -275,15 +276,17 @@ SUBROUTINE mynnedmf_wrapper_run( &
& Tsq, Qsq, Cov, exch_h, exch_m
real(kind=kind_phys), dimension(:), intent(in) :: xmu
real(kind=kind_phys), dimension(:,:), intent(in) :: htrsw, htrlw
! spp_wts_pbl only allocated if spp_pbl == 1
real(kind_phys), dimension(:,:), intent(in) :: spp_wts_pbl

!LOCAL
real(kind=kind_phys), dimension(im,levs) :: &
& sqv,sqc,sqi,qnc,qni,ozone,qnwfa,qnifa, &
& dz, w, p, rho, th, qv, &
& RUBLTEN, RVBLTEN, RTHBLTEN, RQVBLTEN, &
& RQCBLTEN, RQNCBLTEN, RQIBLTEN, RQNIBLTEN, &
& RQNWFABLTEN, RQNIFABLTEN, &
& dqke,qWT,qSHEAR,qBUOY,qDISS, &
& pattern_spp_pbl
& dqke,qWT,qSHEAR,qBUOY,qDISS
real(kind=kind_phys), allocatable :: old_ozone(:,:)

!MYNN-CHEM arrays
Expand Down Expand Up @@ -525,9 +528,9 @@ SUBROUTINE mynnedmf_wrapper_run( &
! qi(i,k)=qi(i,k)/(1.0 - qvsh(i,k))
rho(i,k)=prsl(i,k)/(r_d*t3d(i,k))
w(i,k) = -omega(i,k)/(rho(i,k)*g)
pattern_spp_pbl(i,k)=0.0
enddo
enddo

do i=1,im
if (slmsk(i)==1. .or. slmsk(i)==2.) then !sea/land/ice mask (=0/1/2) in FV3
xland(i)=1.0 !but land/water = (1/2) in SFCLAY_mynn
Expand Down Expand Up @@ -703,7 +706,7 @@ SUBROUTINE mynnedmf_wrapper_run( &
& ,det_thl3D=det_thl,det_sqv3D=det_sqv &
& ,nupdraft=nupdraft,maxMF=maxMF & !output
& ,ktop_plume=ktop_plume & !output
& ,spp_pbl=spp_pbl,pattern_spp_pbl=pattern_spp_pbl & !input
& ,spp_pbl=spp_pbl,pattern_spp_pbl=spp_wts_pbl & !input
& ,RTHRATEN=htrlw & !input
& ,FLAG_QI=flag_qi,FLAG_QNI=flag_qni & !input
& ,FLAG_QC=flag_qc,FLAG_QNC=flag_qnc & !input
Expand Down
15 changes: 15 additions & 0 deletions physics/module_MYNNPBL_wrapper.meta
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1257,6 +1257,21 @@
dimensions = ()
type = logical
intent = in
[spp_wts_pbl]
standard_name = spp_weights_for_pbl_scheme
long_name = spp weights for pbl scheme
units = 1
dimensions = (horizontal_loop_extent,vertical_layer_dimension)
type = real
kind = kind_phys
intent = in
[spp_pbl]
standard_name = control_for_pbl_spp_perturbations
long_name = control for pbl spp perturbations
units = count
dimensions = ()
type = integer
intent = in
[lprnt]
standard_name = flag_print
long_name = control flag for diagnostic print out
Expand Down
11 changes: 7 additions & 4 deletions physics/module_MYNNSFC_wrapper.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,6 +88,7 @@ SUBROUTINE mynnsfc_wrapper_run( &
& FLHC, FLQC, &
& U10, V10, TH2, T2, Q2, &
& wstar, CHS2, CQS2, &
& spp_wts_sfc, spp_sfc, &
! & CP, G, ROVCP, R, XLV, &
! & SVP1, SVP2, SVP3, SVPT0, &
! & EP1,EP2,KARMAN, &
Expand Down Expand Up @@ -143,7 +144,6 @@ SUBROUTINE mynnsfc_wrapper_run( &

!MISC CONFIGURATION OPTIONS
INTEGER, PARAMETER :: &
& spp_pbl = 0, &
& isftcflx = 0, & !control: 0
& iz0tlnd = 0, & !control: 0
& isfflx = 1
Expand All @@ -155,12 +155,15 @@ SUBROUTINE mynnsfc_wrapper_run( &
integer, intent(in) :: ivegsrc
integer, intent(in) :: sfc_z0_type ! option for calculating surface roughness length over ocean
logical, intent(in) :: redrag ! reduced drag coeff. flag for high wind over sea (j.han)
integer, intent(in) :: spp_sfc ! flag for using SPP perturbations

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

!Input data
integer, dimension(:), intent(in) :: vegtype
real(kind=kind_phys), dimension(:), intent(in) :: &
& sigmaf,shdmax,z0pert,ztpert
real(kind_phys), dimension(:,:), intent(in) :: spp_wts_sfc

real(kind=kind_phys), dimension(:,:), &
& intent(in) :: phii
Expand Down Expand Up @@ -207,7 +210,7 @@ SUBROUTINE mynnsfc_wrapper_run( &
& cpm, qgh, qfx, snowh_wat

real(kind=kind_phys), dimension(im,levs) :: &
& pattern_spp_pbl, dz, th, qv
& dz, th, qv

!MYNN-1D
INTEGER :: k, i
Expand All @@ -228,7 +231,6 @@ SUBROUTINE mynnsfc_wrapper_run( &
! endif

! prep MYNN-only variables
pattern_spp_pbl(:,:) = 0
dz(:,:) = 0
th(:,:) = 0
qv(:,:) = 0
Expand All @@ -243,6 +245,7 @@ SUBROUTINE mynnsfc_wrapper_run( &
qv(i,k)=qvsh(i,k)/(1.0 - qvsh(i,k))
enddo
enddo

do i=1,im
if (slmsk(i)==1. .or. slmsk(i)==2.)then !sea/land/ice mask (=0/1/2) in FV3
xland(i)=1.0 !but land/water = (1/2) in SFCLAY_mynn
Expand Down Expand Up @@ -336,7 +339,7 @@ SUBROUTINE mynnsfc_wrapper_run( &
QGH=qgh,QSFC=qsfc, &
U10=u10,V10=v10,TH2=th2,T2=t2,Q2=q2, &
GZ1OZ0=GZ1OZ0,WSPD=wspd,wstar=wstar, &
spp_pbl=spp_pbl,pattern_spp_pbl=pattern_spp_pbl, &
spp_sfc=spp_sfc,pattern_spp_sfc=spp_wts_sfc, &
ids=1,ide=im, jds=1,jde=1, kds=1,kde=levs, &
ims=1,ime=im, jms=1,jme=1, kms=1,kme=levs, &
its=1,ite=im, jts=1,jte=1, kts=1,kte=levs, &
Expand Down
Loading

0 comments on commit ff6395c

Please sign in to comment.