From b99666221163823d433e9b05677e0c77cb5f9140 Mon Sep 17 00:00:00 2001 From: Laurie Carson Date: Thu, 11 Nov 2021 09:33:48 -0700 Subject: [PATCH 01/28] Add SPP option to several physics parameterizations --- physics/GFS_rrtmg_pre.F90 | 26 ++++++++++++++++- physics/GFS_rrtmg_pre.meta | 17 +++++++++++ physics/drag_suite.F90 | 23 ++++++++++++++- physics/drag_suite.meta | 17 +++++++++++ physics/module_MYNNPBL_wrapper.F90 | 20 ++++++++----- physics/module_MYNNPBL_wrapper.meta | 17 +++++++++++ physics/module_MYNNSFC_wrapper.F90 | 16 ++++++---- physics/module_MYNNSFC_wrapper.meta | 17 +++++++++++ physics/mp_thompson.F90 | 45 +++++++++++++++++++---------- physics/mp_thompson.meta | 17 +++++++++++ physics/unified_ugwp.F90 | 8 +++-- physics/unified_ugwp.meta | 17 +++++++++++ 12 files changed, 209 insertions(+), 31 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index dbea66985..8bb4f2073 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -35,7 +35,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & gasvmr_o2, gasvmr_co, gasvmr_cfc11, gasvmr_cfc12, gasvmr_cfc22, & gasvmr_ccl4, gasvmr_cfc113, aerodp, clouds6, clouds7, clouds8, & clouds9, cldsa, cldfra, faersw1, faersw2, faersw3, faerlw1, faerlw2, & - faerlw3, alpha, errmsg, errflg) + faerlw3, alpha, spp_wts_rad, do_spp, errmsg, errflg) use machine, only: kind_phys @@ -102,6 +102,9 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & uni_cld, effr_in, do_mynnedmf, & lmfshal, lmfdeep2, pert_clds + logical, optional, intent(in) :: do_spp + 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 @@ -1086,6 +1089,27 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & enddo enddo +! --- add spp + if ( do_spp .ne. 0 ) then + + do k=1,lm + if (k < levs) then + do i=1,im + effrl_inout(i,k) = effrl_inout(i,k) - spp_wts_rad(i,k) * effrl_inout(i,k) + effri_inout(i,k) = effri_inout(i,k) - spp_wts_rad(i,k) * effri_inout(i,k) + effrs_inout(i,k) = effrs_inout(i,k) - spp_wts_rad(i,k) * effrs_inout(i,k) + enddo + else + do i=1,im + effrl_inout(i,k) = effrl_inout(i,k) - spp_wts_rad(i,levs) * effrl_inout(i,k) + effri_inout(i,k) = effri_inout(i,k) - spp_wts_rad(i,levs) * effri_inout(i,k) + effrs_inout(i,k) = effrs_inout(i,k) - spp_wts_rad(i,levs) * effrs_inout(i,k) + enddo + endif + enddo + + endif + ! mg, sfc-perts ! --- scale random patterns for surface perturbations with ! perturbation size diff --git a/physics/GFS_rrtmg_pre.meta b/physics/GFS_rrtmg_pre.meta index 48ddc586d..247cad3b5 100644 --- a/physics/GFS_rrtmg_pre.meta +++ b/physics/GFS_rrtmg_pre.meta @@ -1169,6 +1169,23 @@ kind = kind_phys intent = out optional = F +[spp_wts_rad] + standard_name = weights_for_stochastic_spp_rad_perturbation + long_name = weights for stochastic spp rad perturbation + units = none + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[do_spp] + standard_name = flag_for_stochastic_spp_option + long_name = flag for stochastic spp option + units = flag + dimensions = () + type = logical + intent = in + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/drag_suite.F90 b/physics/drag_suite.F90 index 9b110d689..cf5fb6616 100644 --- a/physics/drag_suite.F90 +++ b/physics/drag_suite.F90 @@ -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, do_spp, errmsg, errflg) ! ******************************************************************** ! -----> I M P L E M E N T A T I O N V E R S I O N <---------- @@ -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(:,:) + logical, intent(in) :: do_spp + real(kind=kind_phys), dimension(im) :: rstoch !Output: @@ -602,6 +608,21 @@ subroutine drag_suite_run( & end if ! if (me==master) print *,"in Drag Suite, ss_taper:",ss_taper +! SPP, if do_spp is false, no perturbations are applied. +if ( do_spp ) 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 + var_stoch = var + varss_stoch = varss + varmax_ss_stoch = varmax_ss + varmax_fd_stoch = varmax_fd +endif + !--- calculate length of grid for flow-blocking drag ! delx = dx(1) diff --git a/physics/drag_suite.meta b/physics/drag_suite.meta index a36fc1e82..acb9fe9a2 100644 --- a/physics/drag_suite.meta +++ b/physics/drag_suite.meta @@ -703,6 +703,23 @@ type = logical intent = in optional = F +[spp_wts_gwd] + standard_name = weights_for_stochastic_spp_gwd_perturbation + long_name = weights for stochastic spp gwd perturbation + units = none + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[do_spp] + standard_name = flag_for_stochastic_spp_option + long_name = flag for stochastic spp option + units = flag + dimensions = () + type = logical + intent = in + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/module_MYNNPBL_wrapper.F90 b/physics/module_MYNNPBL_wrapper.F90 index 4b034f588..32cf6044d 100644 --- a/physics/module_MYNNPBL_wrapper.F90 +++ b/physics/module_MYNNPBL_wrapper.F90 @@ -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, errmsg, errflg ) + & ltaerosol, spp_wts_pbl, do_spp, lprnt, errmsg, errflg ) ! should be moved to inside the mynn: use machine , only : kind_phys @@ -195,6 +195,7 @@ SUBROUTINE mynnedmf_wrapper_run( & ! NAMELIST OPTIONS (INPUT): LOGICAL, INTENT(IN) :: bl_mynn_tkeadvect, ltaerosol, & lprnt, do_mynnsfclay, & + do_spp, & flag_for_pbl_generic_tend INTEGER, INTENT(IN) :: & & bl_mynn_cloudpdf, & @@ -221,7 +222,6 @@ SUBROUTINE mynnedmf_wrapper_run( & !MISC CONFIGURATION OPTIONS INTEGER, PARAMETER :: & - & spp_pbl=0, & & bl_mynn_mixscalars=1, & & levflag=2 LOGICAL :: & @@ -231,7 +231,8 @@ 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 + & spp_pbl, & + & p_qc, p_qr, p_qi, p_qs, p_qg, p_qnc, p_qni !MYNN-1D REAL(kind=kind_phys), intent(in) :: delt, dtf @@ -275,6 +276,9 @@ 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 do_spp == .true. + 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, & @@ -282,8 +286,7 @@ SUBROUTINE mynnedmf_wrapper_run( & & 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 @@ -525,9 +528,12 @@ 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 + if ( do_spp ) then + spp_pbl=1 + endif + 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 @@ -708,7 +714,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 diff --git a/physics/module_MYNNPBL_wrapper.meta b/physics/module_MYNNPBL_wrapper.meta index e88975aff..04d3fdec6 100644 --- a/physics/module_MYNNPBL_wrapper.meta +++ b/physics/module_MYNNPBL_wrapper.meta @@ -1427,6 +1427,23 @@ type = logical intent = in optional = F +[spp_wts_pbl] + standard_name = weights_for_stochastic_spp_pbl_perturbation + long_name = weights for stochastic spp pbl perturbation + units = none + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[do_spp] + standard_name = flag_for_stochastic_spp_option + long_name = flag for stochastic spp option + units = flag + dimensions = () + type = logical + intent = in + optional = F [lprnt] standard_name = flag_print long_name = control flag for diagnostic print out diff --git a/physics/module_MYNNSFC_wrapper.F90 b/physics/module_MYNNSFC_wrapper.F90 index 271ca5a24..ac060a783 100644 --- a/physics/module_MYNNSFC_wrapper.F90 +++ b/physics/module_MYNNSFC_wrapper.F90 @@ -88,6 +88,7 @@ SUBROUTINE mynnsfc_wrapper_run( & & FLHC, FLQC, & & U10, V10, TH2, T2, Q2, & & wstar, CHS2, CQS2, & + & spp_wts_sfc, do_spp, & ! & CP, G, ROVCP, R, XLV, & ! & SVP1, SVP2, SVP3, SVPT0, & ! & EP1,EP2,KARMAN, & @@ -134,7 +135,6 @@ SUBROUTINE mynnsfc_wrapper_run( & !MISC CONFIGURATION OPTIONS INTEGER, PARAMETER :: & - & spp_pbl = 0, & & isftcflx = 0, & !control: 0 & iz0tlnd = 0, & !control: 0 & isfflx = 1 @@ -145,12 +145,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) + logical, intent(in) :: do_spp ! 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 @@ -197,10 +200,10 @@ SUBROUTINE mynnsfc_wrapper_run( & & cpm, qgh, qfx, qsfc_ruc, snowh_wat real(kind=kind_phys), dimension(im,levs) :: & - & pattern_spp_pbl, dz, th, qv + & dz, th, qv !MYNN-1D - INTEGER :: k, i + INTEGER :: k, i, spp_sfc INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE, & & IMS,IME,JMS,JME,KMS,KME, & & ITS,ITE,JTS,JTE,KTS,KTE @@ -224,9 +227,12 @@ SUBROUTINE mynnsfc_wrapper_run( & th(i,k)=t3d(i,k)/exner(i,k) !qc(i,k)=MAX(qgrs(i,k,ntcw),0.0) qv(i,k)=qvsh(i,k)/(1.0 - qvsh(i,k)) - pattern_spp_pbl(i,k)=0.0 enddo enddo + if ( do_spp ) then + spp_sfc=1 + endif + 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 @@ -321,7 +327,7 @@ SUBROUTINE mynnsfc_wrapper_run( & QGH=qgh,QSFC=qsfc,QSFC_RUC=qsfc_ruc, & 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_pbl=spp_pbl,pattern_spp_pbl=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 ) diff --git a/physics/module_MYNNSFC_wrapper.meta b/physics/module_MYNNSFC_wrapper.meta index b91a026e3..d3471be4f 100644 --- a/physics/module_MYNNSFC_wrapper.meta +++ b/physics/module_MYNNSFC_wrapper.meta @@ -931,6 +931,23 @@ kind = kind_phys intent = inout optional = F +[spp_wts_sfc] + standard_name = weights_for_stochastic_spp_sfc_perturbation + long_name = weights for stochastic spp sfc perturbation + units = none + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[do_spp] + standard_name = flag_for_stochastic_spp_option + long_name = flag for stochastic spp option + units = flag + dimensions = () + type = logical + intent = in + optional = F [lprnt] standard_name = flag_print long_name = control flag for diagnostic print out diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index c31d90b09..f43cfd6ae 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -341,6 +341,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & re_cloud, re_ice, re_snow, & mpicomm, mpirank, mpiroot, blkno, & ext_diag, diag3d, reset_diag3d, & + spp_wts_mp, do_spp, & errmsg, errflg) implicit none @@ -436,12 +437,23 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & integer :: has_reqc integer :: has_reqi integer :: has_reqs - ! DH* 2020-06-05 hardcode these values for not using random perturbations, - ! hasn't been tested yet with this version of module_mp_thompson.F90 - integer, parameter :: rand_perturb_on = 0 integer, parameter :: kme_stoch = 1 - !real(kind_phys) :: rand_pert(1:ncol,1:kme_stoch) - ! *DH 2020-06-05 + logical, intent(in ) :: do_spp + real(kind_phys), intent(in) :: spp_wts_mp(:,:) + !+---+-----------------------------------------------------------------+ + !gthompsn 21Mar2018 + ! Setting spp_mp to 1 gives graupel Y-intercept pertubations (2^0) + ! 2 gives cloud water distribution gamma shape parameter + ! perturbations (2^1) + ! 4 gives CCN & IN activation perturbations (2^2) + ! 3 gives both 1+2 + ! 5 gives both 1+4 + ! 6 gives both 2+4 + ! 7 gives all 1+2+4 + ! For now (22Mar2018), standard deviation should be only 0.25 and cut-off at 1.5 + ! in order to constrain the various perturbations from being too extreme. + !+---+-----------------------------------------------------------------+ + integer :: spp_mp ! Dimensions used in mp_gt_driver integer :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & @@ -499,6 +511,13 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & return end if + ! Set stochastic physics selection to apply all perturbations + if ( do_spp ) then + spp_mp=7 + else + spp_mp=0 + endif + ! Set reduced time step if subcycling is used if (nsteps>1) then dtstep = dtp/real(nsteps, kind=kind_phys) @@ -683,10 +702,8 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & diagflag=diagflag, do_radar_ref=do_radar_ref_mp, & re_cloud=re_cloud, re_ice=re_ice, re_snow=re_snow, & has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & - rand_perturb_on=rand_perturb_on, kme_stoch=kme_stoch, & - ! DH* 2020-06-05 not passing this optional argument, see - ! comment in module_mp_thompson.F90 / mp_gt_driver - !rand_pert=rand_pert, & + rand_perturb_on=spp_mp, kme_stoch=kme_stoch, & + rand_pert=spp_wts_mp, & ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, & ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, & its=its, ite=ite, jts=jts, jte=jte, kts=kts, kte=kte, & @@ -722,10 +739,8 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & refl_10cm=refl_10cm, & diagflag=diagflag, do_radar_ref=do_radar_ref_mp, & has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & - rand_perturb_on=rand_perturb_on, kme_stoch=kme_stoch, & - ! DH* 2020-06-05 not passing this optional argument, see - ! comment in module_mp_thompson.F90 / mp_gt_driver - !rand_pert=rand_pert, & + rand_perturb_on=spp_mp, kme_stoch=kme_stoch, & + rand_pert=spp_wts_mp, & ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, & ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, & its=its, ite=ite, jts=jts, jte=jte, kts=kts, kte=kte, & @@ -763,7 +778,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & diagflag=diagflag, do_radar_ref=do_radar_ref_mp, & re_cloud=re_cloud, re_ice=re_ice, re_snow=re_snow, & has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & - rand_perturb_on=rand_perturb_on, kme_stoch=kme_stoch, & + rand_perturb_on=spp_mp, kme_stoch=kme_stoch, & ! DH* 2020-06-05 not passing this optional argument, see ! comment in module_mp_thompson.F90 / mp_gt_driver !rand_pert=rand_pert, & @@ -801,7 +816,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & refl_10cm=refl_10cm, & diagflag=diagflag, do_radar_ref=do_radar_ref_mp, & has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & - rand_perturb_on=rand_perturb_on, kme_stoch=kme_stoch, & + rand_perturb_on=spp_mp, kme_stoch=kme_stoch, & ! DH* 2020-06-05 not passing this optional argument, see ! comment in module_mp_thompson.F90 / mp_gt_driver !rand_pert=rand_pert, & diff --git a/physics/mp_thompson.meta b/physics/mp_thompson.meta index ab00e6524..ec7d30af5 100644 --- a/physics/mp_thompson.meta +++ b/physics/mp_thompson.meta @@ -758,6 +758,23 @@ type = logical intent = in optional = F +[spp_wts_mp] + standard_name = weights_for_stochastic_spp_mp_perturbation + long_name = weights for stochastic spp mp perturbation + units = none + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[do_spp] + standard_name = flag_for_stochastic_spp_option + long_name = flag for stochastic spp option + units = flag + dimensions = () + type = logical + intent = in + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/unified_ugwp.F90 b/physics/unified_ugwp.F90 index da79ecde8..491a4168c 100644 --- a/physics/unified_ugwp.F90 +++ b/physics/unified_ugwp.F90 @@ -218,7 +218,7 @@ subroutine unified_ugwp_run(me, master, im, levs, ntrac, dtp, fhzero, kdt, index_of_process_orographic_gwd, index_of_process_nonorographic_gwd, & ldiag3d, lssav, flag_for_gwd_generic_tend, do_ugwp_v0, do_ugwp_v0_orog_only, & do_ugwp_v0_nst_only, do_gsl_drag_ls_bl, do_gsl_drag_ss, do_gsl_drag_tofd, & - gwd_opt, errmsg, errflg) + gwd_opt, spp_wts_gwd, do_spp, errmsg, errflg) implicit none @@ -296,6 +296,9 @@ subroutine unified_ugwp_run(me, master, im, levs, ntrac, dtp, fhzero, kdt, do_gsl_drag_ls_bl, do_gsl_drag_ss, & do_gsl_drag_tofd + real(kind=kind_phys), intent(in) :: spp_wts_gwd(:,:) + logical, intent(in) :: do_spp + character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -343,7 +346,8 @@ subroutine unified_ugwp_run(me, master, im, levs, ntrac, dtp, fhzero, kdt, 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,do_spp, & + errmsg, errflg) ! ! put zeros due to xy GSL-drag style: dtaux2d_bl,dtauy2d_bl,dtaux2d_ss.......dusfc_ls,dvsfc_ls ! diff --git a/physics/unified_ugwp.meta b/physics/unified_ugwp.meta index b2f35e45f..2857d7012 100644 --- a/physics/unified_ugwp.meta +++ b/physics/unified_ugwp.meta @@ -1330,6 +1330,23 @@ type = integer intent = in optional = F +[spp_wts_gwd] + standard_name = weights_for_stochastic_spp_gwd_perturbation + long_name = weights for stochastic spp gwd perturbation + units = none + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[do_spp] + standard_name = flag_for_stochastic_spp_option + long_name = flag for stochastic spp option + units = flag + dimensions = () + type = logical + intent = in + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP From 094737206c34ed5462c155e482eed13a46ca1eac Mon Sep 17 00:00:00 2001 From: Laurie Carson Date: Mon, 22 Nov 2021 09:34:50 -0700 Subject: [PATCH 02/28] bug fixes --- physics/module_MYNNSFC_wrapper.F90 | 2 +- physics/module_mp_thompson.F90 | 18 ++++-------------- physics/mp_thompson.F90 | 8 ++------ 3 files changed, 7 insertions(+), 21 deletions(-) diff --git a/physics/module_MYNNSFC_wrapper.F90 b/physics/module_MYNNSFC_wrapper.F90 index ac060a783..d8c6bd38b 100644 --- a/physics/module_MYNNSFC_wrapper.F90 +++ b/physics/module_MYNNSFC_wrapper.F90 @@ -327,7 +327,7 @@ SUBROUTINE mynnsfc_wrapper_run( & QGH=qgh,QSFC=qsfc,QSFC_RUC=qsfc_ruc, & U10=u10,V10=v10,TH2=th2,T2=t2,Q2=q2, & GZ1OZ0=GZ1OZ0,WSPD=wspd,wstar=wstar, & - spp_pbl=spp_pbl,pattern_spp_pbl=spp_wts_sfc, & + spp_pbl=spp_sfc,pattern_spp_pbl=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 ) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index f05aa8ba2..90b70a66f 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -1119,23 +1119,13 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & ! No need to test for every subcycling step test_only_once: if (first_time_step .and. istep==1) then - ! DH* 2020-06-05: The stochastic perturbations code was retrofitted - ! from a newer version of the Thompson MP scheme, but it has not been - ! tested yet. - if (rand_perturb_on .ne. 0) then - errmsg = 'Logic error in mp_gt_driver: the stochastic perturbations code ' // & - 'has not been tested yet with this version of the Thompson scheme' + ! Activate this code when removing the guard above + if (rand_perturb_on .ne. 0 .and. .not. present(rand_pert)) then + errmsg = 'Logic error in mp_gt_driver: random perturbations are on, ' // & + 'but optional argument rand_pert is not present' errflg = 1 return end if - ! Activate this code when removing the guard above - !if (rand_perturb_on .ne. 0 .and. .not. present(rand_pert)) then - ! errmsg = 'Logic error in mp_gt_driver: random perturbations are on, ' // & - ! 'but optional argument rand_pert is not present' - ! errflg = 1 - ! return - !end if - ! *DH 2020-06-05 if ( (present(tt) .and. (present(th) .or. present(pii))) .or. & (.not.present(tt) .and. .not.(present(th) .and. present(pii))) ) then diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index f43cfd6ae..5ae812d32 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -779,9 +779,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & re_cloud=re_cloud, re_ice=re_ice, re_snow=re_snow, & has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & rand_perturb_on=spp_mp, kme_stoch=kme_stoch, & - ! DH* 2020-06-05 not passing this optional argument, see - ! comment in module_mp_thompson.F90 / mp_gt_driver - !rand_pert=rand_pert, & + rand_pert=spp_wts_mp, & ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, & ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, & its=its, ite=ite, jts=jts, jte=jte, kts=kts, kte=kte, & @@ -817,9 +815,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & diagflag=diagflag, do_radar_ref=do_radar_ref_mp, & has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & rand_perturb_on=spp_mp, kme_stoch=kme_stoch, & - ! DH* 2020-06-05 not passing this optional argument, see - ! comment in module_mp_thompson.F90 / mp_gt_driver - !rand_pert=rand_pert, & + rand_pert=spp_wts_mp, & ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, & ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, & its=its, ite=ite, jts=jts, jte=jte, kts=kts, kte=kte, & From 9fcd6da96fe74d6fd24b638d1743d527528781ef Mon Sep 17 00:00:00 2001 From: Laurie Carson Date: Tue, 23 Nov 2021 12:28:46 -0700 Subject: [PATCH 03/28] checkout dtc/ccpp branch (not main), for rte-rrtmgp --- physics/rte-rrtmgp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/rte-rrtmgp b/physics/rte-rrtmgp index 9c01f0e22..56c549450 160000 --- a/physics/rte-rrtmgp +++ b/physics/rte-rrtmgp @@ -1 +1 @@ -Subproject commit 9c01f0e225f493b5b993aee709a59bec1c445e86 +Subproject commit 56c549450787cc7f592a66b4005a299244056568 From 12f92fb9bb3b260e4395a79fa46fc791a64d95e1 Mon Sep 17 00:00:00 2001 From: Laurie Carson Date: Tue, 23 Nov 2021 12:33:35 -0700 Subject: [PATCH 04/28] remove optional keyword from metadata files --- physics/GFS_rrtmg_pre.meta | 3 --- physics/drag_suite.meta | 3 --- physics/module_MYNNPBL_wrapper.meta | 3 --- physics/module_MYNNSFC_wrapper.meta | 3 --- physics/mp_thompson.meta | 3 --- physics/unified_ugwp.meta | 3 --- 6 files changed, 18 deletions(-) diff --git a/physics/GFS_rrtmg_pre.meta b/physics/GFS_rrtmg_pre.meta index a4e5a0d62..7049ffed7 100644 --- a/physics/GFS_rrtmg_pre.meta +++ b/physics/GFS_rrtmg_pre.meta @@ -1034,7 +1034,6 @@ type = real kind = kind_phys intent = out - optional = F [spp_wts_rad] standard_name = weights_for_stochastic_spp_rad_perturbation long_name = weights for stochastic spp rad perturbation @@ -1043,7 +1042,6 @@ type = real kind = kind_phys intent = in - optional = F [do_spp] standard_name = flag_for_stochastic_spp_option long_name = flag for stochastic spp option @@ -1051,7 +1049,6 @@ dimensions = () type = logical intent = in - optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/drag_suite.meta b/physics/drag_suite.meta index df27b51af..4ef257f95 100644 --- a/physics/drag_suite.meta +++ b/physics/drag_suite.meta @@ -624,7 +624,6 @@ dimensions = () type = logical intent = in - optional = F [spp_wts_gwd] standard_name = weights_for_stochastic_spp_gwd_perturbation long_name = weights for stochastic spp gwd perturbation @@ -633,7 +632,6 @@ type = real kind = kind_phys intent = in - optional = F [do_spp] standard_name = flag_for_stochastic_spp_option long_name = flag for stochastic spp option @@ -641,7 +639,6 @@ dimensions = () type = logical intent = in - optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/module_MYNNPBL_wrapper.meta b/physics/module_MYNNPBL_wrapper.meta index 8efeea11b..db05bf4d8 100644 --- a/physics/module_MYNNPBL_wrapper.meta +++ b/physics/module_MYNNPBL_wrapper.meta @@ -1264,7 +1264,6 @@ dimensions = () type = logical intent = in - optional = F [spp_wts_pbl] standard_name = weights_for_stochastic_spp_pbl_perturbation long_name = weights for stochastic spp pbl perturbation @@ -1273,7 +1272,6 @@ type = real kind = kind_phys intent = in - optional = F [do_spp] standard_name = flag_for_stochastic_spp_option long_name = flag for stochastic spp option @@ -1281,7 +1279,6 @@ dimensions = () type = logical intent = in - optional = F [lprnt] standard_name = flag_print long_name = control flag for diagnostic print out diff --git a/physics/module_MYNNSFC_wrapper.meta b/physics/module_MYNNSFC_wrapper.meta index 901256efe..91f55ac74 100644 --- a/physics/module_MYNNSFC_wrapper.meta +++ b/physics/module_MYNNSFC_wrapper.meta @@ -827,7 +827,6 @@ type = real kind = kind_phys intent = inout - optional = F [spp_wts_sfc] standard_name = weights_for_stochastic_spp_sfc_perturbation long_name = weights for stochastic spp sfc perturbation @@ -836,7 +835,6 @@ type = real kind = kind_phys intent = in - optional = F [do_spp] standard_name = flag_for_stochastic_spp_option long_name = flag for stochastic spp option @@ -844,7 +842,6 @@ dimensions = () type = logical intent = in - optional = F [lprnt] standard_name = flag_print long_name = control flag for diagnostic print out diff --git a/physics/mp_thompson.meta b/physics/mp_thompson.meta index 245c68da5..52daa1c42 100644 --- a/physics/mp_thompson.meta +++ b/physics/mp_thompson.meta @@ -638,7 +638,6 @@ dimensions = () type = logical intent = in - optional = F [spp_wts_mp] standard_name = weights_for_stochastic_spp_mp_perturbation long_name = weights for stochastic spp mp perturbation @@ -647,7 +646,6 @@ type = real kind = kind_phys intent = in - optional = F [do_spp] standard_name = flag_for_stochastic_spp_option long_name = flag for stochastic spp option @@ -655,7 +653,6 @@ dimensions = () type = logical intent = in - optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/unified_ugwp.meta b/physics/unified_ugwp.meta index eff265ba6..5ed4e9370 100644 --- a/physics/unified_ugwp.meta +++ b/physics/unified_ugwp.meta @@ -1178,7 +1178,6 @@ dimensions = () type = integer intent = in - optional = F [spp_wts_gwd] standard_name = weights_for_stochastic_spp_gwd_perturbation long_name = weights for stochastic spp gwd perturbation @@ -1187,7 +1186,6 @@ type = real kind = kind_phys intent = in - optional = F [do_spp] standard_name = flag_for_stochastic_spp_option long_name = flag for stochastic spp option @@ -1195,7 +1193,6 @@ dimensions = () type = logical intent = in - optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP From a851742e606313644316e1d57fa23b6cfad4d71b Mon Sep 17 00:00:00 2001 From: Laurie Carson Date: Tue, 23 Nov 2021 12:37:31 -0700 Subject: [PATCH 05/28] specific commit hash for rte-rrtmgp --- physics/rte-rrtmgp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/rte-rrtmgp b/physics/rte-rrtmgp index 56c549450..9c51cb7c3 160000 --- a/physics/rte-rrtmgp +++ b/physics/rte-rrtmgp @@ -1 +1 @@ -Subproject commit 56c549450787cc7f592a66b4005a299244056568 +Subproject commit 9c51cb7c3e227c9e84c2bff29ce4f438c7a54ae6 From d5414d2a3ac26ec49c639c29c18ab9d8db12dc46 Mon Sep 17 00:00:00 2001 From: Laurie Carson Date: Tue, 23 Nov 2021 14:02:05 -0700 Subject: [PATCH 06/28] Add spp args to gwd call, when called directly (vs unified gwd) --- physics/ugwpv1_gsldrag.F90 | 9 +++++++-- physics/ugwpv1_gsldrag.meta | 15 +++++++++++++++ 2 files changed, 22 insertions(+), 2 deletions(-) diff --git a/physics/ugwpv1_gsldrag.F90 b/physics/ugwpv1_gsldrag.F90 index 71193ed88..72cc5e6f5 100644 --- a/physics/ugwpv1_gsldrag.F90 +++ b/physics/ugwpv1_gsldrag.F90 @@ -321,7 +321,8 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ntrac, lonr, dtp, fhzero,kd zogw, zlwb, zobl, zngw, dusfcg, dvsfcg, dudt, dvdt, dtdt, rdxzb, & dtend, dtidx, index_of_x_wind, index_of_y_wind, index_of_temperature, & index_of_process_orographic_gwd, index_of_process_nonorographic_gwd, & - lprnt, ipr, errmsg, errflg) + lprnt, ipr, spp_wts_gwd, do_spp, errmsg, errflg) + ! !######################################################################## ! Attention New Arrays and Names must be ADDED inside @@ -436,6 +437,9 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ntrac, lonr, dtp, fhzero,kd real(kind=kind_phys), intent(out), dimension(:) :: rdxzb ! for stoch phys. mtb-level + real(kind=kind_phys), intent(in) :: spp_wts_gwd(:,:) + logical, intent(in) :: do_spp + character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -558,7 +562,8 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ntrac, lonr, dtp, fhzero,kd 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,do_spp, & + errmsg, errflg) ! ! dusfcg = du_ogwcol + du_oblcol + du_osscol + du_ofdcol ! diff --git a/physics/ugwpv1_gsldrag.meta b/physics/ugwpv1_gsldrag.meta index 64d6b0d64..f4cdf685b 100644 --- a/physics/ugwpv1_gsldrag.meta +++ b/physics/ugwpv1_gsldrag.meta @@ -1082,6 +1082,21 @@ dimensions = () type = integer intent = in +[spp_wts_gwd] + standard_name = weights_for_stochastic_spp_gwd_perturbation + long_name = weights for stochastic spp gwd perturbation + units = none + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[do_spp] + standard_name = flag_for_stochastic_spp_option + long_name = flag for stochastic spp option + units = flag + dimensions = () + type = logical + intent = in [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP From 5523ff7d2870392b0afb7605975652b149273a5d Mon Sep 17 00:00:00 2001 From: jeff beck Date: Tue, 28 Dec 2021 04:14:14 +0000 Subject: [PATCH 07/28] Add missing drag_suite.F90 SPP code --- physics/drag_suite.F90 | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/physics/drag_suite.F90 b/physics/drag_suite.F90 index 7c20a9cf0..cd45577fd 100644 --- a/physics/drag_suite.F90 +++ b/physics/drag_suite.F90 @@ -732,7 +732,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 @@ -888,7 +888,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 ! @@ -898,7 +898,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) @@ -982,7 +982,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) @@ -1001,8 +1001,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) + & + MAX(0.,beta_ss*(varss_stoch(i)-varmax_ss_stoch)) ! 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) @@ -1016,8 +1016,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) + & + MAX(0.,beta_ss*(varss_stoch(i)-varmax_ss_stoch)) ! 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) @@ -1081,8 +1081,8 @@ 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) + & + MAX(0.,beta_fd*(varss_stoch(i)-varmax_fd_stoch)) var_temp = MIN(var_temp, 250.) a1=0.00026615161*var_temp**2 ! a1=0.00026615161*MIN(varss(i),varmax)**2 @@ -1090,7 +1090,7 @@ subroutine drag_suite_run( & ! 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) From c2e8bbf25b4b0dc1c0cb4208a1c35c8d6d8ae3e4 Mon Sep 17 00:00:00 2001 From: jeff beck Date: Tue, 28 Dec 2021 04:30:31 +0000 Subject: [PATCH 08/28] Requested changes from code review. --- physics/GFS_rrtmg_pre.F90 | 2 +- physics/mp_thompson.F90 | 6 ++++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index 893a0a002..5254aa843 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -104,7 +104,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & uni_cld, effr_in, do_mynnedmf, & lmfshal, lmfdeep2, pert_clds - logical, optional, intent(in) :: do_spp + logical, intent(in) :: do_spp real(kind_phys), intent(in) :: spp_wts_rad(:,:) real(kind=kind_phys), intent(in) :: fhswr, fhlwr, solhr, sup, julian, sppt_amp diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index 9e263ea19..6843a6cd0 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -373,6 +373,10 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & ! CCPP error handling character(len=*), intent( out) :: errmsg integer, intent( out) :: errflg + + ! SPP + logical, intent(in) :: do_spp + real(kind_phys), intent(in) :: spp_wts_mp(:,:) ! Local variables @@ -403,8 +407,6 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & integer, parameter :: has_reqi = 0 integer, parameter :: has_reqs = 0 integer, parameter :: kme_stoch = 1 - logical, intent(in ) :: do_spp - real(kind_phys), intent(in) :: spp_wts_mp(:,:) integer :: spp_mp ! Dimensions used in mp_gt_driver integer :: ids,ide, jds,jde, kds,kde, & From 307a507f24829210b88e7abf7e6872cdacc50462 Mon Sep 17 00:00:00 2001 From: jeff beck Date: Tue, 28 Dec 2021 18:40:43 +0000 Subject: [PATCH 09/28] Fix dimension-related bug. --- physics/drag_suite.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/physics/drag_suite.F90 b/physics/drag_suite.F90 index cd45577fd..df3ecb254 100644 --- a/physics/drag_suite.F90 +++ b/physics/drag_suite.F90 @@ -1001,8 +1001,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_stoch(i),varmax_ss_stoch) + & - MAX(0.,beta_ss*(varss_stoch(i)-varmax_ss_stoch)) + 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) @@ -1016,8 +1016,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_stoch(i),varmax_ss_stoch) + & - MAX(0.,beta_ss*(varss_stoch(i)-varmax_ss_stoch)) + 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) @@ -1081,8 +1081,8 @@ 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_stoch(i),varmax_fd_stoch) + & - MAX(0.,beta_fd*(varss_stoch(i)-varmax_fd_stoch)) + 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 From e44b9aa4bc7985fe1d9d9e7602a266e5544d3ce1 Mon Sep 17 00:00:00 2001 From: jeff beck Date: Wed, 12 Jan 2022 20:23:29 +0000 Subject: [PATCH 10/28] Changes for PR modification requests --- physics/GFS_rrtmg_pre.F90 | 6 +-- physics/GFS_rrtmg_pre.meta | 14 +++---- physics/drag_suite.F90 | 8 ++-- physics/drag_suite.meta | 14 +++---- physics/module_MYNNPBL_wrapper.F90 | 12 ++---- physics/module_MYNNPBL_wrapper.meta | 14 +++---- physics/module_MYNNSFC_wrapper.F90 | 11 ++---- physics/module_MYNNSFC_wrapper.meta | 14 +++---- physics/module_mp_thompson.F90 | 2 +- physics/module_sf_mynn.F90 | 60 ++++++++++++++--------------- physics/mp_thompson.F90 | 16 ++++---- physics/mp_thompson.meta | 14 +++---- physics/ugwpv1_gsldrag.F90 | 6 +-- physics/ugwpv1_gsldrag.meta | 14 +++---- physics/unified_ugwp.F90 | 6 +-- physics/unified_ugwp.meta | 14 +++---- 16 files changed, 109 insertions(+), 116 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index b1f56c21f..3ad790614 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -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, & - spp_wts_rad, do_spp, errmsg, errflg) + spp_wts_rad, spp_rad, errmsg, errflg) use machine, only: kind_phys @@ -104,7 +104,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & uni_cld, effr_in, do_mynnedmf, & lmfshal, lmfdeep2, pert_clds - logical, intent(in) :: do_spp + 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 @@ -1080,7 +1080,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & enddo ! --- add spp - if ( do_spp .ne. 0 ) then + if ( spp_rad==1 ) then do k=1,lm if (k < levs) then diff --git a/physics/GFS_rrtmg_pre.meta b/physics/GFS_rrtmg_pre.meta index f0568458f..c9aa88aae 100644 --- a/physics/GFS_rrtmg_pre.meta +++ b/physics/GFS_rrtmg_pre.meta @@ -1083,19 +1083,19 @@ kind = kind_phys intent = out [spp_wts_rad] - standard_name = weights_for_stochastic_spp_rad_perturbation - long_name = weights for stochastic spp rad perturbation + standard_name = weights_for_stochastic_spp_rad_perturbations + long_name = weights for stochastic spp rad perturbations units = none dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = in -[do_spp] - standard_name = flag_for_stochastic_spp_option - long_name = flag for stochastic spp option - units = flag +[spp_rad] + standard_name = control_for_radiation_spp_perturbations + long_name = control for radiation spp perturbations + units = count dimensions = () - type = logical + type = integer intent = in [errmsg] standard_name = ccpp_error_message diff --git a/physics/drag_suite.F90 b/physics/drag_suite.F90 index df3ecb254..919484e87 100644 --- a/physics/drag_suite.F90 +++ b/physics/drag_suite.F90 @@ -219,7 +219,7 @@ subroutine drag_suite_run( & & dtend, dtidx, index_of_process_orographic_gwd, & & index_of_temperature, index_of_x_wind, & & index_of_y_wind, ldiag3d, & - & spp_wts_gwd, do_spp, errmsg, errflg) + & 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 <---------- @@ -369,7 +369,7 @@ subroutine drag_suite_run( & 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(:,:) - logical, intent(in) :: do_spp + integer, intent(in) :: spp_gwd real(kind=kind_phys), dimension(im) :: rstoch @@ -601,8 +601,8 @@ subroutine drag_suite_run( & endif enddo -! SPP, if do_spp is false, no perturbations are applied. -if ( do_spp ) then +! 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) diff --git a/physics/drag_suite.meta b/physics/drag_suite.meta index 7f2c50237..e5ed1e9f3 100644 --- a/physics/drag_suite.meta +++ b/physics/drag_suite.meta @@ -625,19 +625,19 @@ type = logical intent = in [spp_wts_gwd] - standard_name = weights_for_stochastic_spp_gwd_perturbation - long_name = weights for stochastic spp gwd perturbation + standard_name = weights_for_stochastic_spp_gwd_perturbations + long_name = weights for stochastic spp gwd perturbations units = none dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = in -[do_spp] - standard_name = flag_for_stochastic_spp_option - long_name = flag for stochastic spp option - units = flag +[spp_gwd] + standard_name = control_for_gravity_wave_drag_spp_perturbations + long_name = control for gravity wave drag spp perturbations + units = count dimensions = () - type = logical + type = integer intent = in [errmsg] standard_name = ccpp_error_message diff --git a/physics/module_MYNNPBL_wrapper.F90 b/physics/module_MYNNPBL_wrapper.F90 index 0aa2648f3..7c0ba1ba4 100644 --- a/physics/module_MYNNPBL_wrapper.F90 +++ b/physics/module_MYNNPBL_wrapper.F90 @@ -108,7 +108,7 @@ SUBROUTINE mynnedmf_wrapper_run( & & icloud_bl, do_mynnsfclay, & & imp_physics, imp_physics_gfdl, & & imp_physics_thompson, imp_physics_wsm6, & - & ltaerosol, spp_wts_pbl, do_spp, 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 @@ -195,7 +195,6 @@ SUBROUTINE mynnedmf_wrapper_run( & ! NAMELIST OPTIONS (INPUT): LOGICAL, INTENT(IN) :: bl_mynn_tkeadvect, ltaerosol, & lprnt, do_mynnsfclay, & - do_spp, & flag_for_pbl_generic_tend INTEGER, INTENT(IN) :: & & bl_mynn_cloudpdf, & @@ -210,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(:,:,:) @@ -232,7 +232,6 @@ SUBROUTINE mynnedmf_wrapper_run( & LOGICAL, PARAMETER :: cycling = .false. INTEGER, PARAMETER :: param_first_scalar = 1 INTEGER :: & - & spp_pbl, & & p_qc, p_qr, p_qi, p_qs, p_qg, p_qnc, p_qni !MYNN-1D @@ -277,7 +276,7 @@ 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 do_spp == .true. + ! spp_wts_pbl only allocated if spp_pbl == 1 real(kind_phys), dimension(:,:), intent(in) :: spp_wts_pbl !LOCAL @@ -531,9 +530,6 @@ SUBROUTINE mynnedmf_wrapper_run( & w(i,k) = -omega(i,k)/(rho(i,k)*g) enddo enddo - if ( do_spp ) then - spp_pbl=1 - endif do i=1,im if (slmsk(i)==1. .or. slmsk(i)==2.) then !sea/land/ice mask (=0/1/2) in FV3 diff --git a/physics/module_MYNNPBL_wrapper.meta b/physics/module_MYNNPBL_wrapper.meta index 0d024498b..8f313bcbc 100644 --- a/physics/module_MYNNPBL_wrapper.meta +++ b/physics/module_MYNNPBL_wrapper.meta @@ -1258,19 +1258,19 @@ type = logical intent = in [spp_wts_pbl] - standard_name = weights_for_stochastic_spp_pbl_perturbation - long_name = weights for stochastic spp pbl perturbation + standard_name = weights_for_stochastic_spp_pbl_perturbations + long_name = weights for stochastic spp pbl perturbations units = none dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = in -[do_spp] - standard_name = flag_for_stochastic_spp_option - long_name = flag for stochastic spp option - units = flag +[spp_pbl] + standard_name = control_for_planetary_boundary_layer_spp_perturbations + long_name = control for planetary boundary layer spp perturbations + units = count dimensions = () - type = logical + type = integer intent = in [lprnt] standard_name = flag_print diff --git a/physics/module_MYNNSFC_wrapper.F90 b/physics/module_MYNNSFC_wrapper.F90 index 4032f8b30..e49891a18 100644 --- a/physics/module_MYNNSFC_wrapper.F90 +++ b/physics/module_MYNNSFC_wrapper.F90 @@ -88,7 +88,7 @@ SUBROUTINE mynnsfc_wrapper_run( & & FLHC, FLQC, & & U10, V10, TH2, T2, Q2, & & wstar, CHS2, CQS2, & - & spp_wts_sfc, do_spp, & + & spp_wts_sfc, spp_sfc, & ! & CP, G, ROVCP, R, XLV, & ! & SVP1, SVP2, SVP3, SVPT0, & ! & EP1,EP2,KARMAN, & @@ -155,7 +155,7 @@ 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) - logical, intent(in) :: do_spp ! flag for using SPP perturbations + integer, intent(in) :: spp_sfc ! flag for using SPP perturbations real(kind=kind_phys), intent(in) :: delt @@ -213,7 +213,7 @@ SUBROUTINE mynnsfc_wrapper_run( & & dz, th, qv !MYNN-1D - INTEGER :: k, i, spp_sfc + INTEGER :: k, i INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE, & & IMS,IME,JMS,JME,KMS,KME, & & ITS,ITE,JTS,JTE,KTS,KTE @@ -239,9 +239,6 @@ SUBROUTINE mynnsfc_wrapper_run( & qv(i,k)=qvsh(i,k)/(1.0 - qvsh(i,k)) enddo enddo - if ( do_spp ) then - spp_sfc=1 - endif do i=1,im if (slmsk(i)==1. .or. slmsk(i)==2.)then !sea/land/ice mask (=0/1/2) in FV3 @@ -336,7 +333,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_sfc,pattern_spp_pbl=spp_wts_sfc, & + 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 ) diff --git a/physics/module_MYNNSFC_wrapper.meta b/physics/module_MYNNSFC_wrapper.meta index 99468fd7c..c1ef9fe69 100644 --- a/physics/module_MYNNSFC_wrapper.meta +++ b/physics/module_MYNNSFC_wrapper.meta @@ -835,19 +835,19 @@ kind = kind_phys intent = inout [spp_wts_sfc] - standard_name = weights_for_stochastic_spp_sfc_perturbation - long_name = weights for stochastic spp sfc perturbation + standard_name = weights_for_stochastic_spp_sfc_perturbations + long_name = weights for stochastic spp sfc perturbations units = none dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = in -[do_spp] - standard_name = flag_for_stochastic_spp_option - long_name = flag for stochastic spp option - units = flag +[spp_sfc] + standard_name = control_for_surface_layer_spp_perturbations + long_name = control for surface layer spp perturbations + units = count dimensions = () - type = logical + type = integer intent = in [lprnt] standard_name = flag_print diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 5b975aaad..788a3010b 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -1290,7 +1290,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & !+---+-----------------------------------------------------------------+ !..Introduce stochastic parameter perturbations by creating as many scalar rand1, rand2, ... !.. variables as needed to perturb different pieces of microphysics. gthompsn 21Mar2018 -! Setting spp_mp to 1 gives graupel Y-intercept pertubations (2^0) +! Setting spp_mp_opt to 1 gives graupel Y-intercept pertubations (2^0) ! 2 gives cloud water distribution gamma shape parameter perturbations (2^1) ! 4 gives CCN & IN activation perturbations (2^2) ! 3 gives both 1+2 diff --git a/physics/module_sf_mynn.F90 b/physics/module_sf_mynn.F90 index c98dc2169..f2bc66e52 100644 --- a/physics/module_sf_mynn.F90 +++ b/physics/module_sf_mynn.F90 @@ -161,7 +161,7 @@ SUBROUTINE SFCLAY_mynn( & QGH,QSFC, & U10,V10,TH2,T2,Q2, & GZ1OZ0,WSPD,WSTAR, & - spp_pbl,pattern_spp_pbl, & + spp_sfc,pattern_spp_sfc, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) @@ -271,7 +271,7 @@ SUBROUTINE SFCLAY_mynn( & !NAMELIST/CONFIGURATION OPTIONS: INTEGER, INTENT(IN) :: ISFFLX, LSM, LSM_RUC INTEGER, OPTIONAL, INTENT(IN) :: ISFTCFLX, IZ0TLND - INTEGER, OPTIONAL, INTENT(IN) :: spp_pbl, psi_opt + INTEGER, OPTIONAL, INTENT(IN) :: spp_sfc, psi_opt 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) @@ -293,7 +293,7 @@ SUBROUTINE SFCLAY_mynn( & th3d,pi3d REAL, DIMENSION( ims:ime, kms:kme), OPTIONAL, & - INTENT(IN) :: pattern_spp_pbl + INTENT(IN) :: pattern_spp_sfc !=================================== ! 2D VARIABLES !=================================== @@ -396,8 +396,8 @@ SUBROUTINE SFCLAY_mynn( & QC1D(i)=QC3D(i,kts) P1D(i) =P3D(i,kts) T1D(i) =T3D(i,kts) - if (spp_pbl==1) then - rstoch1D(i)=pattern_spp_pbl(i,kts) + if (spp_sfc==1) then + rstoch1D(i)=pattern_spp_sfc(i,kts) else rstoch1D(i)=0.0 endif @@ -462,7 +462,7 @@ SUBROUTINE SFCLAY_mynn( & HFLX,HFX,QFLX,QFX,LH,FLHC,FLQC, & QGH,QSFC,U10,V10,TH2,T2,Q2, & GZ1OZ0,WSPD,wstar, & - spp_pbl,rstoch1D, & + spp_sfc,rstoch1D, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte & @@ -510,7 +510,7 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, & QGH,QSFC, & U10,V10,TH2,T2,Q2, & GZ1OZ0,WSPD,wstar, & - spp_pbl,rstoch1D, & + spp_sfc,rstoch1D, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte & @@ -537,7 +537,7 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, & !----------------------------- INTEGER, INTENT(IN) :: ISFFLX INTEGER, OPTIONAL, INTENT(IN ) :: ISFTCFLX, IZ0TLND - INTEGER, INTENT(IN) :: spp_pbl, psi_opt + INTEGER, INTENT(IN) :: spp_sfc, psi_opt 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) @@ -1111,7 +1111,7 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, & endif !-end wave model check ! add stochastic perturbation of ZNT - if (spp_pbl==1) then + if (spp_sfc==1) then ZNTstoch_wat(I) = MAX(ZNT_wat(I) + ZNT_wat(I)*1.0*rstoch1D(i), 1e-6) else ZNTstoch_wat(I) = ZNT_wat(I) @@ -1140,29 +1140,29 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, & IF ( ISFTCFLX .EQ. 0 ) THEN IF (COARE_OPT .EQ. 3.0) THEN CALL fairall_etal_2003(ZT_wat(i),ZQ_wat(i),restar,UST_wat(i),visc,& - rstoch1D(i),spp_pbl) + rstoch1D(i),spp_sfc) ELSE !presumably, this will be published soon, but hasn't yet CALL fairall_etal_2014(ZT_wat(i),ZQ_wat(i),restar,UST_wat(i),visc,& - rstoch1D(i),spp_pbl) + rstoch1D(i),spp_sfc) ENDIF ELSEIF ( ISFTCFLX .EQ. 1 ) THEN IF (COARE_OPT .EQ. 3.0) THEN CALL fairall_etal_2003(ZT_wat(i),ZQ_wat(i),restar,UST_wat(i),visc,& - rstoch1D(i),spp_pbl) + rstoch1D(i),spp_sfc) ELSE CALL fairall_etal_2014(ZT_wat(i),ZQ_wat(i),restar,UST_wat(i),visc,& - rstoch1D(i),spp_pbl) + rstoch1D(i),spp_sfc) ENDIF ELSEIF ( ISFTCFLX .EQ. 2 ) THEN CALL garratt_1992(ZT_wat(i),ZQ_wat(i),ZNTstoch_wat(i),restar,2.0) ELSEIF ( ISFTCFLX .EQ. 3 ) THEN IF (COARE_OPT .EQ. 3.0) THEN CALL fairall_etal_2003(ZT_wat(i),ZQ_wat(i),restar,UST_wat(i),visc,& - rstoch1D(i),spp_pbl) + rstoch1D(i),spp_sfc) ELSE CALL fairall_etal_2014(ZT_wat(i),ZQ_wat(i),restar,UST_wat(i),visc,& - rstoch1D(i),spp_pbl) + rstoch1D(i),spp_sfc) ENDIF ELSEIF ( ISFTCFLX .EQ. 4 ) THEN !GFS zt formulation @@ -1173,10 +1173,10 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, & !DEFAULT TO COARE 3.0/3.5 IF (COARE_OPT .EQ. 3.0) THEN CALL fairall_etal_2003(ZT_wat(i),ZQ_wat(i),restar,UST_wat(i),visc,& - rstoch1D(i),spp_pbl) + rstoch1D(i),spp_sfc) ELSE CALL fairall_etal_2014(ZT_wat(i),ZQ_wat(i),restar,UST_wat(i),visc,& - rstoch1D(i),spp_pbl) + rstoch1D(i),spp_sfc) ENDIF ENDIF IF (debug_code > 1) THEN @@ -1201,7 +1201,7 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, & endif ! add stochastic perturbaction of ZNT - if (spp_pbl==1) then + if (spp_sfc==1) then ZNTstoch_lnd(I) = MAX(ZNT_lnd(I) + ZNT_lnd(I)*1.0*rstoch1D(i), 1e-6) else ZNTstoch_lnd(I) = ZNT_lnd(I) @@ -1222,7 +1222,7 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, & IF ( PRESENT(IZ0TLND) ) THEN IF ( IZ0TLND .LE. 1 ) THEN CALL zilitinkevich_1995(ZNTstoch_lnd(i),ZT_lnd(i),ZQ_lnd(i),restar,& - UST_lnd(I),KARMAN,1.0,IZ0TLND,spp_pbl,rstoch1D(i)) + UST_lnd(I),KARMAN,1.0,IZ0TLND,spp_sfc,rstoch1D(i)) ELSEIF ( IZ0TLND .EQ. 2 ) THEN CALL Yang_2008(ZNTSTOCH_lnd(i),ZT_lnd(i),ZQ_lnd(i),UST_lnd(i),MOL(I),& qstar(I),restar,visc) @@ -1237,7 +1237,7 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, & ELSE !DEFAULT TO ZILITINKEVICH CALL zilitinkevich_1995(ZNTSTOCH_lnd(i),ZT_lnd(i),ZQ_lnd(i),restar,& - UST_lnd(I),KARMAN,1.0,0,spp_pbl,rstoch1D(i)) + UST_lnd(I),KARMAN,1.0,0,spp_sfc,rstoch1D(i)) ENDIF ENDIF IF (ZNTstoch_lnd(i) < 1E-8 .OR. Zt_lnd(i) < 1E-10) THEN @@ -1263,7 +1263,7 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, & IF (icy(I)) THEN ! add stochastic perturbaction of ZNT - if (spp_pbl==1) then + if (spp_sfc==1) then ZNTstoch_ice(I) = MAX(ZNT_ice(I) + ZNT_ice(I)*1.0*rstoch1D(i), 1e-6) else ZNTstoch_ice(I) = ZNT_ice(I) @@ -2246,7 +2246,7 @@ END SUBROUTINE SFCLAY1D_mynn !! to work with the Noah LSM and may be specific for that !! LSM only. Tests with RUC LSM showed no improvements. SUBROUTINE zilitinkevich_1995(Z_0,Zt,Zq,restar,ustar,KARMAN,& - & landsea,IZ0TLND2,spp_pbl,rstoch) + & landsea,IZ0TLND2,spp_sfc,rstoch) IMPLICIT NONE REAL, INTENT(IN) :: Z_0,restar,ustar,KARMAN,landsea @@ -2255,7 +2255,7 @@ SUBROUTINE zilitinkevich_1995(Z_0,Zt,Zq,restar,ustar,KARMAN,& REAL :: CZIL !=0.100 in Chen et al. (1997) !=0.075 in Zilitinkevich (1995) !=0.500 in Lemone et al. (2008) - INTEGER, INTENT(IN) :: spp_pbl + INTEGER, INTENT(IN) :: spp_sfc REAL, INTENT(IN) :: rstoch @@ -2296,7 +2296,7 @@ SUBROUTINE zilitinkevich_1995(Z_0,Zt,Zq,restar,ustar,KARMAN,& ! stochastically perturb thermal and moisture roughness length. ! currently set to half the amplitude: - if (spp_pbl==1) then + if (spp_sfc==1) then Zt = Zt + Zt * 0.5 * rstoch Zt = MAX(Zt, 0.0001) Zq = Zt @@ -2461,11 +2461,11 @@ END SUBROUTINE garratt_1992 !!(1992, p. 102), is available for flows with Ren < 2. !! !!This is for use over water only. - SUBROUTINE fairall_etal_2003(Zt,Zq,Ren,ustar,visc,rstoch,spp_pbl) + SUBROUTINE fairall_etal_2003(Zt,Zq,Ren,ustar,visc,rstoch,spp_sfc) IMPLICIT NONE REAL, INTENT(IN) :: Ren,ustar,visc,rstoch - INTEGER, INTENT(IN):: spp_pbl + INTEGER, INTENT(IN):: spp_sfc REAL, INTENT(OUT) :: Zt,Zq IF (Ren .le. 2.) then @@ -2484,7 +2484,7 @@ SUBROUTINE fairall_etal_2003(Zt,Zq,Ren,ustar,visc,rstoch,spp_pbl) ENDIF - if (spp_pbl==1) then + if (spp_sfc==1) then Zt = Zt + Zt * 0.5 * rstoch Zq = Zt endif @@ -2505,18 +2505,18 @@ END SUBROUTINE fairall_etal_2003 !! COARE 3.5/4.0 formulation, empirically derived from COARE and HEXMAX data !! [Fairall et al. (2014? coming soon, not yet published as of July 2014)]. !! This is for use over water only. - SUBROUTINE fairall_etal_2014(Zt,Zq,Ren,ustar,visc,rstoch,spp_pbl) + SUBROUTINE fairall_etal_2014(Zt,Zq,Ren,ustar,visc,rstoch,spp_sfc) IMPLICIT NONE REAL, INTENT(IN) :: Ren,ustar,visc,rstoch - INTEGER, INTENT(IN):: spp_pbl + INTEGER, INTENT(IN):: spp_sfc REAL, INTENT(OUT) :: Zt,Zq !Zt = (5.5e-5)*(Ren**(-0.60)) Zt = MIN(1.6E-4, 5.8E-5/(Ren**0.72)) Zq = Zt - IF (spp_pbl ==1) THEN + IF (spp_sfc ==1) THEN Zt = MAX(Zt + Zt*0.5*rstoch,2.0e-9) Zq = MAX(Zt + Zt*0.5*rstoch,2.0e-9) ELSE diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index 6843a6cd0..5316eb181 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -308,7 +308,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & refl_10cm, reset_dBZ, do_radar_ref, & mpicomm, mpirank, mpiroot, blkno, & ext_diag, diag3d, reset_diag3d, & - spp_wts_mp, do_spp, & + spp_wts_mp, spp_mp, & errmsg, errflg) implicit none @@ -375,7 +375,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & integer, intent( out) :: errflg ! SPP - logical, intent(in) :: do_spp + integer, intent(in) :: spp_mp real(kind_phys), intent(in) :: spp_wts_mp(:,:) ! Local variables @@ -407,7 +407,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & integer, parameter :: has_reqi = 0 integer, parameter :: has_reqs = 0 integer, parameter :: kme_stoch = 1 - integer :: spp_mp + integer :: spp_mp_opt ! Dimensions used in mp_gt_driver integer :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & @@ -466,10 +466,10 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & end if ! Set stochastic physics selection to apply all perturbations - if ( do_spp ) then - spp_mp=7 + if ( spp_mp ) then + spp_mp_opt=7 else - spp_mp=0 + spp_mp_opt=0 endif ! Set reduced time step if subcycling is used @@ -633,7 +633,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & refl_10cm=refl_10cm, & diagflag=diagflag, do_radar_ref=do_radar_ref_mp, & has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & - rand_perturb_on=spp_mp, kme_stoch=kme_stoch, & + rand_perturb_on=spp_mp_opt, kme_stoch=kme_stoch, & rand_pert=spp_wts_mp, & ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, & ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, & @@ -670,7 +670,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & refl_10cm=refl_10cm, & diagflag=diagflag, do_radar_ref=do_radar_ref_mp, & has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & - rand_perturb_on=spp_mp, kme_stoch=kme_stoch, & + rand_perturb_on=spp_mp_opt, kme_stoch=kme_stoch, & rand_pert=spp_wts_mp, & ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, & ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, & diff --git a/physics/mp_thompson.meta b/physics/mp_thompson.meta index 52daa1c42..896c01873 100644 --- a/physics/mp_thompson.meta +++ b/physics/mp_thompson.meta @@ -639,19 +639,19 @@ type = logical intent = in [spp_wts_mp] - standard_name = weights_for_stochastic_spp_mp_perturbation - long_name = weights for stochastic spp mp perturbation + standard_name = weights_for_stochastic_spp_mp_perturbations + long_name = weights for stochastic spp mp perturbations units = none dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = in -[do_spp] - standard_name = flag_for_stochastic_spp_option - long_name = flag for stochastic spp option - units = flag +[spp_mp] + standard_name = control_for_microphysics_spp_perturbations + long_name = control for microphysics spp perturbations + units = count dimensions = () - type = logical + type = integer intent = in [errmsg] standard_name = ccpp_error_message diff --git a/physics/ugwpv1_gsldrag.F90 b/physics/ugwpv1_gsldrag.F90 index 72cc5e6f5..3c4f3400f 100644 --- a/physics/ugwpv1_gsldrag.F90 +++ b/physics/ugwpv1_gsldrag.F90 @@ -321,7 +321,7 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ntrac, lonr, dtp, fhzero,kd zogw, zlwb, zobl, zngw, dusfcg, dvsfcg, dudt, dvdt, dtdt, rdxzb, & dtend, dtidx, index_of_x_wind, index_of_y_wind, index_of_temperature, & index_of_process_orographic_gwd, index_of_process_nonorographic_gwd, & - lprnt, ipr, spp_wts_gwd, do_spp, errmsg, errflg) + lprnt, ipr, spp_wts_gwd, spp_gwd, errmsg, errflg) ! !######################################################################## @@ -438,7 +438,7 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ntrac, lonr, dtp, fhzero,kd real(kind=kind_phys), intent(out), dimension(:) :: rdxzb ! for stoch phys. mtb-level real(kind=kind_phys), intent(in) :: spp_wts_gwd(:,:) - logical, intent(in) :: do_spp + logical, intent(in) :: spp_gwd character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -562,7 +562,7 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ntrac, lonr, dtp, fhzero,kd 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, spp_wts_gwd,do_spp, & + index_of_y_wind, ldiag3d, spp_wts_gwd, spp_gwd, & errmsg, errflg) ! ! dusfcg = du_ogwcol + du_oblcol + du_osscol + du_ofdcol diff --git a/physics/ugwpv1_gsldrag.meta b/physics/ugwpv1_gsldrag.meta index edce8c201..97113d20c 100644 --- a/physics/ugwpv1_gsldrag.meta +++ b/physics/ugwpv1_gsldrag.meta @@ -1083,19 +1083,19 @@ type = integer intent = in [spp_wts_gwd] - standard_name = weights_for_stochastic_spp_gwd_perturbation - long_name = weights for stochastic spp gwd perturbation + standard_name = weights_for_stochastic_spp_gwd_perturbations + long_name = weights for stochastic spp gwd perturbations units = none dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = in -[do_spp] - standard_name = flag_for_stochastic_spp_option - long_name = flag for stochastic spp option - units = flag +[spp_gwd] + standard_name = control_for_gravity_wave_drag_spp_perturbations + long_name = control for gravity wave drag spp perturbations + units = count dimensions = () - type = logical + type = integer intent = in [errmsg] standard_name = ccpp_error_message diff --git a/physics/unified_ugwp.F90 b/physics/unified_ugwp.F90 index 491a4168c..9e93bd5fc 100644 --- a/physics/unified_ugwp.F90 +++ b/physics/unified_ugwp.F90 @@ -218,7 +218,7 @@ subroutine unified_ugwp_run(me, master, im, levs, ntrac, dtp, fhzero, kdt, index_of_process_orographic_gwd, index_of_process_nonorographic_gwd, & ldiag3d, lssav, flag_for_gwd_generic_tend, do_ugwp_v0, do_ugwp_v0_orog_only, & do_ugwp_v0_nst_only, do_gsl_drag_ls_bl, do_gsl_drag_ss, do_gsl_drag_tofd, & - gwd_opt, spp_wts_gwd, do_spp, errmsg, errflg) + gwd_opt, spp_wts_gwd, spp_gwd, errmsg, errflg) implicit none @@ -297,7 +297,7 @@ subroutine unified_ugwp_run(me, master, im, levs, ntrac, dtp, fhzero, kdt, do_gsl_drag_tofd real(kind=kind_phys), intent(in) :: spp_wts_gwd(:,:) - logical, intent(in) :: do_spp + integer, intent(in) :: spp_gwd character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -346,7 +346,7 @@ subroutine unified_ugwp_run(me, master, im, levs, ntrac, dtp, fhzero, kdt, 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, spp_wts_gwd,do_spp, & + index_of_y_wind, ldiag3d, spp_wts_gwd, spp_gwd, & errmsg, errflg) ! ! put zeros due to xy GSL-drag style: dtaux2d_bl,dtauy2d_bl,dtaux2d_ss.......dusfc_ls,dvsfc_ls diff --git a/physics/unified_ugwp.meta b/physics/unified_ugwp.meta index d919e0b9e..46e79ea33 100644 --- a/physics/unified_ugwp.meta +++ b/physics/unified_ugwp.meta @@ -1179,19 +1179,19 @@ type = integer intent = in [spp_wts_gwd] - standard_name = weights_for_stochastic_spp_gwd_perturbation - long_name = weights for stochastic spp gwd perturbation + standard_name = weights_for_stochastic_spp_gwd_perturbations + long_name = weights for stochastic spp gwd perturbations units = none dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = in -[do_spp] - standard_name = flag_for_stochastic_spp_option - long_name = flag for stochastic spp option - units = flag +[spp_gwd] + standard_name = control_for_gravity_wave_drag_spp_perturbations + long_name = control for gravity wave drag spp perturbations + units = count dimensions = () - type = logical + type = integer intent = in [errmsg] standard_name = ccpp_error_message From 8d599bc0608c30675568cc81424d16905646944c Mon Sep 17 00:00:00 2001 From: jeff beck Date: Thu, 13 Jan 2022 02:45:27 +0000 Subject: [PATCH 11/28] Fix bug in SPP implementation --- physics/GFS_rrtmg_pre.F90 | 34 ++++++++++++++++------------------ 1 file changed, 16 insertions(+), 18 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index 3ad790614..c63221075 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -1080,24 +1080,22 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & enddo ! --- add spp - if ( spp_rad==1 ) then - - do k=1,lm - if (k < levs) then - do i=1,im - effrl_inout(i,k) = effrl_inout(i,k) - spp_wts_rad(i,k) * effrl_inout(i,k) - effri_inout(i,k) = effri_inout(i,k) - spp_wts_rad(i,k) * effri_inout(i,k) - effrs_inout(i,k) = effrs_inout(i,k) - spp_wts_rad(i,k) * effrs_inout(i,k) - enddo - else - do i=1,im - effrl_inout(i,k) = effrl_inout(i,k) - spp_wts_rad(i,levs) * effrl_inout(i,k) - effri_inout(i,k) = effri_inout(i,k) - spp_wts_rad(i,levs) * effri_inout(i,k) - effrs_inout(i,k) = effrs_inout(i,k) - spp_wts_rad(i,levs) * effrs_inout(i,k) - enddo - endif - enddo - + if ( spp_rad == 1 ) then + do k=1,lm + if (k < levs) then + do i=1,im + effrl_inout(i,k) = effrl_inout(i,k) - spp_wts_rad(i,k) * effrl_inout(i,k) + effri_inout(i,k) = effri_inout(i,k) - spp_wts_rad(i,k) * effri_inout(i,k) + effrs_inout(i,k) = effrs_inout(i,k) - spp_wts_rad(i,k) * effrs_inout(i,k) + enddo + else + do i=1,im + effrl_inout(i,k) = effrl_inout(i,k) - spp_wts_rad(i,levs) * effrl_inout(i,k) + effri_inout(i,k) = effri_inout(i,k) - spp_wts_rad(i,levs) * effri_inout(i,k) + effrs_inout(i,k) = effrs_inout(i,k) - spp_wts_rad(i,levs) * effrs_inout(i,k) + enddo + endif + enddo endif ! mg, sfc-perts From 6f898244276bf1f3e52f6fc084d4eb7b1963c1c9 Mon Sep 17 00:00:00 2001 From: jeff beck Date: Fri, 14 Jan 2022 01:51:44 +0000 Subject: [PATCH 12/28] Perturb cloud* instead of effr*_inout --- physics/GFS_rrtmg_pre.F90 | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index c63221075..f4f413d0d 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -1079,20 +1079,19 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & enddo enddo -! --- add spp if ( spp_rad == 1 ) then do k=1,lm if (k < levs) then do i=1,im - effrl_inout(i,k) = effrl_inout(i,k) - spp_wts_rad(i,k) * effrl_inout(i,k) - effri_inout(i,k) = effri_inout(i,k) - spp_wts_rad(i,k) * effri_inout(i,k) - effrs_inout(i,k) = effrs_inout(i,k) - spp_wts_rad(i,k) * effrs_inout(i,k) + 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 - effrl_inout(i,k) = effrl_inout(i,k) - spp_wts_rad(i,levs) * effrl_inout(i,k) - effri_inout(i,k) = effri_inout(i,k) - spp_wts_rad(i,levs) * effri_inout(i,k) - effrs_inout(i,k) = effrs_inout(i,k) - spp_wts_rad(i,levs) * effrs_inout(i,k) + 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 From 9c90b4702eff5bc1c0ff924620e3b8a0b06448dd Mon Sep 17 00:00:00 2001 From: jeff beck Date: Fri, 14 Jan 2022 16:25:08 +0000 Subject: [PATCH 13/28] Change spp_gwd from logical to integer --- physics/ugwpv1_gsldrag.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/ugwpv1_gsldrag.F90 b/physics/ugwpv1_gsldrag.F90 index 3c4f3400f..41290ed68 100644 --- a/physics/ugwpv1_gsldrag.F90 +++ b/physics/ugwpv1_gsldrag.F90 @@ -438,7 +438,7 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ntrac, lonr, dtp, fhzero,kd real(kind=kind_phys), intent(out), dimension(:) :: rdxzb ! for stoch phys. mtb-level real(kind=kind_phys), intent(in) :: spp_wts_gwd(:,:) - logical, intent(in) :: spp_gwd + integer, intent(in) :: spp_gwd character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg From a7b987594372c4446103c3afa0310b5b98c45723 Mon Sep 17 00:00:00 2001 From: jeff beck Date: Fri, 14 Jan 2022 19:41:26 +0000 Subject: [PATCH 14/28] MYNN SFC perturbation pattern name fix --- physics/module_MYNNSFC_wrapper.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/module_MYNNSFC_wrapper.F90 b/physics/module_MYNNSFC_wrapper.F90 index 852694329..c767d863a 100644 --- a/physics/module_MYNNSFC_wrapper.F90 +++ b/physics/module_MYNNSFC_wrapper.F90 @@ -231,7 +231,7 @@ SUBROUTINE mynnsfc_wrapper_run( & ! endif ! prep MYNN-only variables - pattern_spp_pbl(:,:) = 0 + pattern_spp_sfc(:,:) = 0 dz(:,:) = 0 th(:,:) = 0 qv(:,:) = 0 From dd419a3db2a9a1b825287c39da42610b1546dcb9 Mon Sep 17 00:00:00 2001 From: jeff beck Date: Fri, 14 Jan 2022 19:54:52 +0000 Subject: [PATCH 15/28] Remove initialization of pattern_spp_sfc=0 --- physics/module_MYNNSFC_wrapper.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/physics/module_MYNNSFC_wrapper.F90 b/physics/module_MYNNSFC_wrapper.F90 index c767d863a..150a66472 100644 --- a/physics/module_MYNNSFC_wrapper.F90 +++ b/physics/module_MYNNSFC_wrapper.F90 @@ -231,7 +231,6 @@ SUBROUTINE mynnsfc_wrapper_run( & ! endif ! prep MYNN-only variables - pattern_spp_sfc(:,:) = 0 dz(:,:) = 0 th(:,:) = 0 qv(:,:) = 0 From d466eb1fee4ffef9821247973d7aed16d5c8801a Mon Sep 17 00:00:00 2001 From: jeff beck Date: Sun, 16 Jan 2022 03:53:54 +0000 Subject: [PATCH 16/28] Updates to standard names and units --- physics/GFS_rrtmg_pre.meta | 6 +++--- physics/drag_suite.meta | 6 +++--- physics/module_MYNNPBL_wrapper.meta | 10 +++++----- physics/module_MYNNSFC_wrapper.meta | 6 +++--- physics/mp_thompson.meta | 6 +++--- physics/ugwpv1_gsldrag.meta | 6 +++--- physics/unified_ugwp.meta | 6 +++--- 7 files changed, 23 insertions(+), 23 deletions(-) diff --git a/physics/GFS_rrtmg_pre.meta b/physics/GFS_rrtmg_pre.meta index c9aa88aae..675fc095a 100644 --- a/physics/GFS_rrtmg_pre.meta +++ b/physics/GFS_rrtmg_pre.meta @@ -1083,9 +1083,9 @@ kind = kind_phys intent = out [spp_wts_rad] - standard_name = weights_for_stochastic_spp_rad_perturbations - long_name = weights for stochastic spp rad perturbations - units = none + 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 diff --git a/physics/drag_suite.meta b/physics/drag_suite.meta index e5ed1e9f3..ae474d1ee 100644 --- a/physics/drag_suite.meta +++ b/physics/drag_suite.meta @@ -625,9 +625,9 @@ type = logical intent = in [spp_wts_gwd] - standard_name = weights_for_stochastic_spp_gwd_perturbations - long_name = weights for stochastic spp gwd perturbations - units = none + 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 diff --git a/physics/module_MYNNPBL_wrapper.meta b/physics/module_MYNNPBL_wrapper.meta index 8f313bcbc..c8c1ca106 100644 --- a/physics/module_MYNNPBL_wrapper.meta +++ b/physics/module_MYNNPBL_wrapper.meta @@ -1258,16 +1258,16 @@ type = logical intent = in [spp_wts_pbl] - standard_name = weights_for_stochastic_spp_pbl_perturbations - long_name = weights for stochastic spp pbl perturbations - units = none + 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_planetary_boundary_layer_spp_perturbations - long_name = control for planetary boundary layer spp perturbations + standard_name = control_for_pbl_spp_perturbations + long_name = control for pbl spp perturbations units = count dimensions = () type = integer diff --git a/physics/module_MYNNSFC_wrapper.meta b/physics/module_MYNNSFC_wrapper.meta index c1ef9fe69..ec0c183ed 100644 --- a/physics/module_MYNNSFC_wrapper.meta +++ b/physics/module_MYNNSFC_wrapper.meta @@ -835,9 +835,9 @@ kind = kind_phys intent = inout [spp_wts_sfc] - standard_name = weights_for_stochastic_spp_sfc_perturbations - long_name = weights for stochastic spp sfc perturbations - units = none + standard_name = spp_weights_for_surface_layer_scheme + long_name = spp weights for surface layer scheme + units = 1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys diff --git a/physics/mp_thompson.meta b/physics/mp_thompson.meta index 896c01873..ddddffb7f 100644 --- a/physics/mp_thompson.meta +++ b/physics/mp_thompson.meta @@ -639,9 +639,9 @@ type = logical intent = in [spp_wts_mp] - standard_name = weights_for_stochastic_spp_mp_perturbations - long_name = weights for stochastic spp mp perturbations - units = none + standard_name = spp_weights_for_microphysics_scheme + long_name = spp weights for microphysics scheme + units = 1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys diff --git a/physics/ugwpv1_gsldrag.meta b/physics/ugwpv1_gsldrag.meta index 34487e8a3..b0ca7c7fb 100644 --- a/physics/ugwpv1_gsldrag.meta +++ b/physics/ugwpv1_gsldrag.meta @@ -1083,9 +1083,9 @@ type = integer intent = in [spp_wts_gwd] - standard_name = weights_for_stochastic_spp_gwd_perturbations - long_name = weights for stochastic spp gwd perturbations - units = none + 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 diff --git a/physics/unified_ugwp.meta b/physics/unified_ugwp.meta index 181f0bd55..26fc6d6f4 100644 --- a/physics/unified_ugwp.meta +++ b/physics/unified_ugwp.meta @@ -1179,9 +1179,9 @@ type = integer intent = in [spp_wts_gwd] - standard_name = weights_for_stochastic_spp_gwd_perturbations - long_name = weights for stochastic spp gwd perturbations - units = none + 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 From f685e2011fa14017571bc51540789e0c7932d21c Mon Sep 17 00:00:00 2001 From: jeff beck Date: Tue, 18 Jan 2022 20:47:55 +0000 Subject: [PATCH 17/28] Fix varmax field dimensions --- physics/drag_suite.F90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/physics/drag_suite.F90 b/physics/drag_suite.F90 index 919484e87..31fb4fd50 100644 --- a/physics/drag_suite.F90 +++ b/physics/drag_suite.F90 @@ -606,8 +606,8 @@ subroutine drag_suite_run( & 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) + varmax_ss_stoch = varmax_ss + varmax_ss*0.75*spp_wts_gwd(i,1) + varmax_fd_stoch = varmax_fd + varmax_fd*0.75*spp_wts_gwd(i,1) enddo else var_stoch = var @@ -1001,8 +1001,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_stoch(i),varmax_ss_stoch(i)) + & - MAX(0.,beta_ss*(varss_stoch(i)-varmax_ss_stoch(i))) + var_temp = MIN(varss_stoch(i),varmax_ss_stoch) + & + MAX(0.,beta_ss*(varss_stoch(i)-varmax_ss_stoch)) ! 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) @@ -1016,8 +1016,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_stoch(i),varmax_ss_stoch(i)) + & - MAX(0.,beta_ss*(varss_stoch(i)-varmax_ss_stoch(i))) + var_temp = MIN(varss_stoch(i),varmax_ss_stoch) + & + MAX(0.,beta_ss*(varss_stoch(i)-varmax_ss_stoch)) ! 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) @@ -1081,8 +1081,8 @@ 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_stoch(i),varmax_fd_stoch(i)) + & - MAX(0.,beta_fd*(varss_stoch(i)-varmax_fd_stoch(i))) + var_temp = MIN(varss_stoch(i),varmax_fd_stoch) + & + MAX(0.,beta_fd*(varss_stoch(i)-varmax_fd_stoch)) var_temp = MIN(var_temp, 250.) a1=0.00026615161*var_temp**2 ! a1=0.00026615161*MIN(varss(i),varmax)**2 From 73eb0f35b3ee2a47f821fe808c50c052207ca9fc Mon Sep 17 00:00:00 2001 From: jeff beck Date: Thu, 20 Jan 2022 04:25:34 +0000 Subject: [PATCH 18/28] Add dimensions to SPP variables --- drag_suite.F90 | 1381 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1381 insertions(+) create mode 100755 drag_suite.F90 diff --git a/drag_suite.F90 b/drag_suite.F90 new file mode 100755 index 000000000..7fea98b13 --- /dev/null +++ b/drag_suite.F90 @@ -0,0 +1,1381 @@ +!> \File drag_suite.F90 +!! This file is the parameterization of orographic gravity wave +!! drag, mountain blocking, and form drag. + +!> This module contains the CCPP-compliant orographic gravity wave dray scheme. + module drag_suite + + contains + + subroutine drag_suite_init(gwd_opt, errmsg, errflg) + + integer, intent(in) :: gwd_opt + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 + + ! Consistency checks + if (gwd_opt/=3 .and. gwd_opt/=33) then + write(errmsg,'(*(a))') "Logic error: namelist choice of gravity wave & + & drag is different from drag_suite scheme" + errflg = 1 + return + end if + end subroutine drag_suite_init + +! \defgroup GFS_ogwd GFS Orographic Gravity Wave Drag +!> \defgroup gfs_drag_suite GFS drag_suite Main +!! \brief This subroutine includes orographic gravity wave drag, mountain +!! blocking, and form drag. +!! +!> The time tendencies of zonal and meridional wind are altered to +!! include the effect of mountain induced gravity wave drag from +!! subgrid scale orography including convective breaking, shear +!! breaking and the presence of critical levels. +!! +!> \section arg_table_drag_suite_run Argument Table +!! \htmlinclude drag_suite_run.html +!! +!> \section gen_drag_suite GFS Orographic GWD Scheme General Algorithm +!! -# Calculate subgrid mountain blocking +!! -# Calculate orographic wave drag +!! +!! The NWP model gravity wave drag (GWD) scheme in the GFS has two +!! main components: how the surface stress is computed, and then how +!! that stress is distributed over a vertical column where it may +!! interact with the models momentum. Each of these depends on the +!! large scale environmental atmospheric state and assumptions about +!! the sub-grid scale processes. In Alpert GWD (1987) based on linear, +!! two-dimensional non-rotating, stably stratified flow over a mountain ridge, +!! sub-grid scale gravity wave motions are assumed which propagate away +!! from the mountain. Described in Alpert (1987), the flux measured over +!! a "low level" vertically averaged layer, in the atmosphere defines a base +!! level flux. "Low level" was taken to be the first 1/3 of the troposphere +!! in the 1987 implementation. This choice was meant to encompass a thick +!! low layer for vertical averages of the environmental (large scale) flow +!! quantities. The vertical momentum flux or gravity wave stress in a +!! grid box due to a single mountain is given as in Pierrehumbert, (1987) (PH): +!! +!! \f$ \tau = \frac {\rho \: U^{3}\: G(F_{r})} {\Delta X \; N } \f$ +!! +!! emetic \f$ \Delta X \f$ is a grid increment, N is the Brunt Viasala frequency +!! +!! +!! \f$ N(\sigma) = \frac{-g \: \sigma \: +!! \frac{\partial\Theta}{\partial\sigma}}{\Theta \:R \:T} \f$ +!! +!! The environmental variables are calculated from a mass weighted vertical +!! average over a base layer. G(Fr) is a monotonically increasing +!! function of Froude number, +!! +!! \f$ F_{r} = \frac{N h^{'}}{U} \f$ +!! +!! where U is the wind speed calculated as a mass weighted vertical average in +!! the base layer, and h', is the vertical displacement caused by the orography +!! variance. An effective mountain length for the gravity wave processes, +!! +!! \f$ l^{*} = \frac{\Delta X}{m} \f$ +!! +!! where m is the number of mountains in a grid box, can then +!! be defined to obtain the form of the base level stress +!! +!! +!! \f$ \tau = \frac {\rho \: U^{3} \: G(F_{r})} {N \;l^{*}} \f$ +!! +!! giving the stress induced from the surface in a model grid box. +!! PH gives the form for the function G(Fr) as +!! +!! +!! \f$ G(F_{r}) = \bar{G}\frac{F^{2}_{r}}{F^{2}_{r}\: + \:a^{2}} \f$ +!! +!! Where \f$ \bar{G} \f$ is an order unity non-dimensional saturation +!! flux set to 1 and 'a' is a function of the mountain aspect ratio also +!!set to 1 in the 1987 implementation of the GFS GWD. Typical values of +!! U=10m/s, N=0.01 1/s, l*=100km, and a=1, gives a flux of 1 Pascal and +!! if this flux is made to go to zero linearly with height then the +!! decelerations would be about 10/m/s/day which is consistent with +!! observations in PH. +!! +!! +!! In Kim, Moorthi, Alpert's (1998, 2001) GWD currently in GFS operations, +!! the GWD scheme has the same physical basis as in Alpert (1987) with the addition +!! of enhancement factors for the amplitude, G, and mountain shape details +!! in G(Fr) to account for effects from the mountain blocking. A factor, +!! E m', is an enhancement factor on the stress in the Alpert '87 scheme. +!! The E ranges from no enhancement to an upper limit of 3, E=E(OA)[1-3], +!! and is a function of OA, the Orographic Asymmetry defined in KA (1995) as +!! +!! Orographic Asymmetry (OA) = \f$ \frac{ \bar{x} \; - \; +!! \sum\limits_{j=1}^{N_{b}} x_{j} \; n_{j} }{\sigma_{x}} \f$ +!! +!! where Nb is the total number of bottom blocks in the mountain barrier, +!! \f$ \sigma_{x} \f$ is the standard deviation of the horizontal distance defined by +!! +!! \f$ \sigma_{x} = \sqrt{ \frac{\sum\limits_{j=1}^{N_{b}} +!! \; (x_{j} \; - \; \bar{x} )^2}{N_{x}} } \f$ +!! +!! +!! where Nx is the number of grid intervals for the large scale domain being +!! considered. So the term, E(OA)m'/ \f$ \Delta X \f$ in Kim's scheme represents +!! a multiplier on G shown in Alpert's eq (1), where m' is the number of mountains +!! in a sub-grid scale box. Kim increased the complexity of m' making it a +!! function of the fractional area of the sub-grid mountain and the asymmetry +!! and convexity statistics which are found from running a gravity wave +!! model for a large number of cases: +!! +!! \f$ m^{'} = C_{m} \Delta X \left[ \frac{1 \; + \; +!! \sum\limits_{x} L_{h} }{\Delta X} \right]^{OA+1} \f$ +!! +!! Where, according to Kim, \f$ \sum \frac{L_{h}}{\Delta X} \f$ is +!! the fractional area covered by the subgrid-scale orography higher than +!! a critical height \f$ h_{c} = Fr_{c} U_{0}/N_{0} \f$ , over the +!! "low level" vertically averaged layer, for a grid box with the interval +!! \f$ \Delta X \f$. Each \f$ L_{n}\f$ is the width of a segment of +!! orography intersection at the critical height: +!! +!! \f$ Fr_{0} = \frac{N_{0} \; h^{'}}{U_{0}} \f$ +!! +!! \f$ G^{'}(OC,Fr_{0}) = \frac{Fr_{0}^{2}}{Fr_{0}^{2} \; + \; a^{2}} \f$ +!! +!! \f$ a^{2} = \frac{C_{G}}{OC} \f$ +!! +!! \f$ E(OA, Fr_{0}) = (OA \; + \; 2)^{\delta} \f$ and \f$ \delta +!! \; = \; \frac{C_{E} \; Fr_{0}}{Fr_{c}} \f$ where \f$ Fr_{c} \f$ +!! is as in Alpert. +!! +!! +!! This represents a closed scheme, somewhat empirical adjustments +!! to the original scheme to calculate the surface stress. +!! +!! Momentum is deposited by the sub-grid scale gravity waves break due +!! to the presence of convective mixing assumed to occur when the +!! minimum Richardson number: +!! +!! Orographic Convexity (OC) = \f$ \frac{ \sum\limits_{j=1}^{N_{x}} +!! \; (h_{j} \; - \; \bar{h})^4 }{N_{x} \;\sigma_{h}^4} \f$ , +!! and where \f$ \sigma_{h} = \sqrt{ \frac{\sum\limits_{j=1}^{N_{x}} +!! \; (h_{j} \; - \; \bar{h} )^2}{N_{x}} } \f$ +!! +!! This represents a closed scheme, somewhat empirical adjustments +!! to the original scheme to calculate the surface stress. +!! +!! Momentum is deposited by the sub-grid scale gravity waves break due +!! to the presence of convective mixing assumed to occur when +!! the minimum Richardson number: +!! +!! \f$ Ri_{m} = \frac{Ri(1 \; - \; Fr)}{(1 \; + \; \sqrt{Ri}Fr)^2} \f$ +!! +!! Is less than 1/4 Or if critical layers are encountered in a layer +!! the the momentum flux will vanish. The critical layer is defined +!! when the base layer wind becomes perpendicular to the environmental +!! wind. Otherwise, wave breaking occurs at a level where the amplification +!! of the wave causes the local Froude number or similarly a truncated +!! (first term of the) Scorer parameter, to be reduced below a critical +!! value by the saturation hypothesis (Lindzen,). This is done through +!! eq 1 which can be written as +!! +!! \f$ \tau = \rho U N k h^{'2} \f$ +!! +!! For small Froude number this is discretized in the vertical so at each +!! level the stress is reduced by ratio of the Froude or truncated Scorer +!! parameter, \f$ \frac{U^{2}}{N^{2}} = \frac{N \tau_{l-1}}{\rho U^{3} k} \f$ , +!! where the stress is from the layer below beginning with that found near +!! the surface. The respective change in momentum is applied in +!! that layer building up from below. +!! +!! An amplitude factor is part of the calibration of this scheme which is +!! a function of the model resolution and the vertical diffusion. This +!! is because the vertical diffusion and the GWD account encompass +!! similar physical processes. Thus, one needs to run the model over +!! and over for various amplitude factors for GWD and vertical diffusion. +!! +!! In addition, there is also mountain blocking from lift and frictional +!! forces. Improved integration between how the GWD is calculated and +!! the mountain blocking of wind flow around sub-grid scale orography +!! is underway at NCEP. The GFS already has convectively forced GWD +!! an independent process. The next step is to test +!! +!> \section det_drag_suite GFS Orographic GWD Scheme Detailed Algorithm +!> @{ + subroutine drag_suite_run( & + & IM,KM,dvdt,dudt,dtdt,U1,V1,T1,Q1,KPBL, & + & PRSI,DEL,PRSL,PRSLK,PHII,PHIL,DELTIM,KDT, & + & var,oc1,oa4,ol4, & + & varss,oc1ss,oa4ss,ol4ss, & + & THETA,SIGMA,GAMMA,ELVMAX, & + & dtaux2d_ls,dtauy2d_ls,dtaux2d_bl,dtauy2d_bl, & + & dtaux2d_ss,dtauy2d_ss,dtaux2d_fd,dtauy2d_fd, & + & dusfc,dvsfc, & + & dusfc_ls,dvsfc_ls,dusfc_bl,dvsfc_bl, & + & dusfc_ss,dvsfc_ss,dusfc_fd,dvsfc_fd, & + & slmsk,br1,hpbl, & + & g, cp, rd, rv, fv, pi, imx, cdmbgwd, me, master, & + & lprnt, ipr, rdxzb, dx, gwd_opt, & + & 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, & + & 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 <---------- +! +! ----- This code ----- +!begin WRF code + +! this code handles the time tendencies of u v due to the effect of mountain +! induced gravity wave drag from sub-grid scale orography. this routine +! not only treats the traditional upper-level wave breaking due to mountain +! variance (alpert 1988), but also the enhanced lower-tropospheric wave +! breaking due to mountain convexity and asymmetry (kim and arakawa 1995). +! thus, in addition to the terrain height data in a model grid box, +! additional 10-2d topographic statistics files are needed, including +! orographic standard deviation (var), convexity (oc1), asymmetry (oa4) +! and ol (ol4). these data sets are prepared based on the 30 sec usgs orography +! hong (1999). the current scheme was implmented as in hong et al.(2008) +! +! Originally coded by song-you hong and young-joon kim and implemented by song-you hong +! +! program history log: +! 2014-10-01 Hyun-Joo Choi (from KIAPS) flow-blocking drag of kim and doyle +! with blocked height by dividing streamline theory +! 2017-04-06 Joseph Olson (from Gert-Jan Steeneveld) added small-scale +! orographic grabity wave drag: +! 2017-09-15 Joseph Olson, with some bug fixes from Michael Toy: added the +! topographic form drag of Beljaars et al. (2004, QJRMS) +! Activation of each component is done by specifying the integer-parameters +! (defined below) to 0: inactive or 1: active +! gwd_opt_ls = 0 or 1: large-scale +! gwd_opt_bl = 0 or 1: blocking drag +! gwd_opt_ss = 0 or 1: small-scale gravity wave drag +! gwd_opt_fd = 0 or 1: topographic form drag +! 2017-09-25 Michael Toy (from NCEP GFS model) added dissipation heating +! gsd_diss_ht_opt = 0: dissipation heating off +! gsd_diss_ht_opt = 1: dissipation heating on +! 2020-08-25 Michael Toy changed logic control for drag component selection +! for CCPP. +! Namelist options: +! do_gsl_drag_ls_bl - logical flag for large-scale GWD + blocking +! do_gsl_drag_ss - logical flag for small-scale GWD +! do_gsl_drag_tofd - logical flag for turbulent form drag +! Compile-time options (same as before): +! gwd_opt_ls = 0 or 1: large-scale GWD +! gwd_opt_bl = 0 or 1: blocking drag +! +! References: +! Hong et al. (2008), wea. and forecasting +! Kim and Doyle (2005), Q. J. R. Meteor. Soc. +! Kim and Arakawa (1995), j. atmos. sci. +! Alpert et al. (1988), NWP conference. +! Hong (1999), NCEP office note 424. +! Steeneveld et al (2008), JAMC +! Tsiringakis et al. (2017), Q. J. R. Meteor. Soc. +! Beljaars et al. (2004), Q. J. R. Meteor. Soc. +! +! notice : comparible or lower resolution orography files than model resolution +! are desirable in preprocess (wps) to prevent weakening of the drag +!------------------------------------------------------------------------------- +! +! input +! dudt (im,km) non-lin tendency for u wind component +! dvdt (im,km) non-lin tendency for v wind component +! u1(im,km) zonal wind / sqrt(rcl) m/sec at t0-dt +! v1(im,km) meridional wind / sqrt(rcl) m/sec at t0-dt +! t1(im,km) temperature deg k at t0-dt +! q1(im,km) specific humidity at t0-dt +! deltim time step secs +! del(km) positive increment of pressure across layer (pa) +! KPBL(IM) is the index of the top layer of the PBL +! ipr & lprnt for diagnostics +! +! output +! dudt, dvdt wind tendency due to gwdo +! dTdt +! +!------------------------------------------------------------------------------- + +!end wrf code +!----------------------------------------------------------------------C +! USE +! ROUTINE IS CALLED FROM CCPP (AFTER CALLING PBL SCHEMES) +! +! PURPOSE +! USING THE GWD PARAMETERIZATIONS OF PS-GLAS AND PH- +! GFDL TECHNIQUE. THE TIME TENDENCIES OF U V +! ARE ALTERED TO INCLUDE THE EFFECT OF MOUNTAIN INDUCED +! GRAVITY WAVE DRAG FROM SUB-GRID SCALE OROGRAPHY INCLUDING +! CONVECTIVE BREAKING, SHEAR BREAKING AND THE PRESENCE OF +! CRITICAL LEVELS +! +! +! ******************************************************************** + USE MACHINE , ONLY : kind_phys + implicit none + + ! Interface variables + integer, intent(in) :: im, km, imx, kdt, ipr, me, master + integer, intent(in) :: gwd_opt + logical, intent(in) :: lprnt + integer, intent(in) :: KPBL(:) + real(kind=kind_phys), intent(in) :: deltim, G, CP, RD, RV, cdmbgwd(:) + real(kind=kind_phys), intent(inout) :: dtend(:,:,:) + logical, intent(in) :: ldiag3d + integer, intent(in) :: dtidx(:,:), index_of_temperature, & + & index_of_process_orographic_gwd, index_of_x_wind, index_of_y_wind + + integer :: kpblmax + integer, parameter :: ims=1, kms=1, its=1, kts=1 + real(kind=kind_phys), intent(in) :: fv, pi + real(kind=kind_phys) :: rcl, cdmb + real(kind=kind_phys) :: g_inv + + real(kind=kind_phys), intent(inout) :: & + & dudt(:,:),dvdt(:,:), & + & dtdt(:,:) + real(kind=kind_phys), intent(out) :: rdxzb(:) + real(kind=kind_phys), intent(in) :: & + & u1(:,:),v1(:,:), & + & t1(:,:),q1(:,:), & + & PHII(:,:),prsl(:,:), & + & prslk(:,:),PHIL(:,:) + real(kind=kind_phys), intent(in) :: prsi(:,:), & + & del(:,:) + real(kind=kind_phys), intent(in) :: var(:),oc1(:), & + & oa4(:,:),ol4(:,:), & + & dx(:) + real(kind=kind_phys), intent(in) :: varss(:),oc1ss(:), & + & oa4ss(:,:),ol4ss(:,:) + real(kind=kind_phys), intent(in) :: THETA(:),SIGMA(:), & + & GAMMA(:),ELVMAX(:) + +! added for small-scale orographic wave drag + real(kind=kind_phys), dimension(im,km) :: utendwave,vtendwave,thx,thvx + real(kind=kind_phys), intent(in) :: br1(:), & + & hpbl(:), & + & slmsk(:) + real(kind=kind_phys), dimension(im) :: govrth,xland + !real(kind=kind_phys), dimension(im,km) :: dz2 + real(kind=kind_phys) :: tauwavex0,tauwavey0, & + & XNBV,density,tvcon,hpbl2 + integer :: kpbl2,kvar + !real(kind=kind_phys), dimension(im,km+1) :: zq ! = PHII/g + 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: + real(kind=kind_phys), intent(out) :: & + & dusfc(:), dvsfc(:) +!Output (optional): + real(kind=kind_phys), intent(out) :: & + & dusfc_ls(:),dvsfc_ls(:), & + & dusfc_bl(:),dvsfc_bl(:), & + & dusfc_ss(:),dvsfc_ss(:), & + & dusfc_fd(:),dvsfc_fd(:) + real(kind=kind_phys), intent(out) :: & + & dtaux2d_ls(:,:),dtauy2d_ls(:,:), & + & dtaux2d_bl(:,:),dtauy2d_bl(:,:), & + & dtaux2d_ss(:,:),dtauy2d_ss(:,:), & + & dtaux2d_fd(:,:),dtauy2d_fd(:,:) + +!Misc arrays + real(kind=kind_phys), dimension(im,km) :: dtaux2d, dtauy2d + +!------------------------------------------------------------------------- +! Flags to regulate the activation of specific components of drag suite: +! Each component is tapered off automatically as a function of dx, so best to +! keep them activated (.true.). + logical, intent(in) :: & + do_gsl_drag_ls_bl, & ! large-scale gravity wave drag and blocking + do_gsl_drag_ss, & ! small-scale gravity wave drag (Steeneveld et al. 2008) + do_gsl_drag_tofd ! form drag (Beljaars et al. 2004, QJRMS) + +! Additional flags + integer, parameter :: & + gwd_opt_ls = 1, & ! large-scale gravity wave drag + gwd_opt_bl = 1, & ! blocking drag + gsd_diss_ht_opt = 0 + +! Parameters for bounding the scale-adaptive variability: +! Small-scale GWD + turbulent form drag + real(kind=kind_phys), parameter :: dxmin_ss = 1000., & + & dxmax_ss = 12000. ! min,max range of tapering (m) +! Large-scale GWD + blocking + real(kind=kind_phys), parameter :: dxmin_ls = 3000., & + & dxmax_ls = 13000. ! min,max range of tapering (m) + real(kind=kind_phys), dimension(im) :: ss_taper, ls_taper ! small- and large-scale tapering factors (-) +! +! Variables for limiting topographic standard deviation (var) + real(kind=kind_phys), parameter :: varmax_ss = 50., & + varmax_fd = 150., & + beta_ss = 0.1, & + beta_fd = 0.2 + real(kind=kind_phys) :: var_temp, var_temp2 + +! added Beljaars orographic form drag + real(kind=kind_phys), dimension(im,km) :: utendform,vtendform + real(kind=kind_phys) :: a1,a2,wsp + real(kind=kind_phys) :: H_efold + +! critical richardson number for wave breaking : ! larger drag with larger value + real(kind=kind_phys), parameter :: ric = 0.25 + real(kind=kind_phys), parameter :: dw2min = 1. + real(kind=kind_phys), parameter :: rimin = -100. + real(kind=kind_phys), parameter :: bnv2min = 1.0e-5 + real(kind=kind_phys), parameter :: efmin = 0.0 + real(kind=kind_phys), parameter :: efmax = 10.0 + real(kind=kind_phys), parameter :: xl = 4.0e4 + real(kind=kind_phys), parameter :: critac = 1.0e-5 + real(kind=kind_phys), parameter :: gmax = 1. + real(kind=kind_phys), parameter :: veleps = 1.0 + real(kind=kind_phys), parameter :: factop = 0.5 + real(kind=kind_phys), parameter :: frc = 1.0 + real(kind=kind_phys), parameter :: ce = 0.8 + real(kind=kind_phys), parameter :: cg = 0.5 + integer,parameter :: kpblmin = 2 + +! +! local variables +! + integer :: i,j,k,lcap,lcapp1,nwd,idir, & + klcap,kp1 +! + real(kind=kind_phys) :: rcs,csg,fdir,cleff,cleff_ss,cs, & + rcsks,wdir,ti,rdz,tem2,dw2,shr2, & + bvf2,rdelks,wtkbj,tem,gfobnv,hd,fro, & + rim,temc,tem1,efact,temv,dtaux,dtauy, & + dtauxb,dtauyb,eng0,eng1 +! + logical :: ldrag(im),icrilv(im), & + flag(im),kloop1(im) +! + real(kind=kind_phys) :: taub(im),taup(im,km+1), & + xn(im),yn(im), & + ubar(im),vbar(im), & + fr(im),ulow(im), & + rulow(im),bnv(im), & + oa(im),ol(im), & + oass(im),olss(im), & + roll(im),dtfac(im), & + brvf(im),xlinv(im), & + delks(im),delks1(im), & + bnv2(im,km),usqj(im,km), & + taud_ls(im,km),taud_bl(im,km), & + ro(im,km), & + vtk(im,km),vtj(im,km), & + zlowtop(im),velco(im,km-1), & + coefm(im),coefm_ss(im) +! + integer :: kbl(im),klowtop(im) + integer,parameter :: mdir=8 + !integer :: nwdir(mdir) + !data nwdir/6,7,5,8,2,3,1,4/ + integer, parameter :: nwdir(8) = (/6,7,5,8,2,3,1,4/) +! +! variables for flow-blocking drag +! + real(kind=kind_phys),parameter :: frmax = 10. + real(kind=kind_phys),parameter :: olmin = 1.0e-5 + real(kind=kind_phys),parameter :: odmin = 0.1 + real(kind=kind_phys),parameter :: odmax = 10. + real(kind=kind_phys),parameter :: erad = 6371.315e+3 + integer :: komax(im) + integer :: kblk + real(kind=kind_phys) :: cd + real(kind=kind_phys) :: zblk,tautem + real(kind=kind_phys) :: pe,ke + real(kind=kind_phys) :: delx,dely + real(kind=kind_phys) :: dxy4(im,4),dxy4p(im,4) + real(kind=kind_phys) :: dxy(im),dxyp(im) + real(kind=kind_phys) :: ol4p(4),olp(im),od(im) + real(kind=kind_phys) :: taufb(im,km+1) + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + integer :: udtend, vdtend, Tdtend + + ! Calculate inverse of gravitational acceleration + g_inv = 1./G + + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 + + ! Initialize local variables + var_temp2 = 0. + udtend = -1 + vdtend = -1 + Tdtend = -1 + + if(ldiag3d) then + udtend = dtidx(index_of_x_wind,index_of_process_orographic_gwd) + vdtend = dtidx(index_of_y_wind,index_of_process_orographic_gwd) + Tdtend = dtidx(index_of_temperature,index_of_process_orographic_gwd) + endif + +!-------------------------------------------------------------------- +! SCALE-ADPTIVE PARAMETER FROM GFS GWD SCHEME +!-------------------------------------------------------------------- +! parameter (cdmb = 1.0) ! non-dim sub grid mtn drag Amp (*j*) +! non-dim sub grid mtn drag Amp (*j*) +! cdmb = 1.0/float(IMX/192) +! cdmb = 192.0/float(IMX) + cdmb = 4.0 * 192.0/float(IMX) + if (cdmbgwd(1) >= 0.0) cdmb = cdmb * cdmbgwd(1) + +!>-# Orographic Gravity Wave Drag Section + kpblmax = km / 2 ! maximum pbl height : # of vertical levels / 2 +! +! Scale cleff between IM=384*2 and 192*2 for T126/T170 and T62 +! + if (imx > 0) then +! cleff = 1.0E-5 * SQRT(FLOAT(IMX)/384.0) +! cleff = 1.0E-5 * SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF! +! cleff = 0.5E-5 * SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF! +! cleff = 1.0E-5 * SQRT(FLOAT(IMX)/192)/float(IMX/192) +! cleff = 1.0E-5 / SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF! + cleff = 0.5E-5 / SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF! +! hmhj for ndsl +! jw cleff = 0.1E-5 / SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF! +! cleff = 2.0E-5 * SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF! +! cleff = 2.5E-5 * SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF! + endif + if (cdmbgwd(2) >= 0.0) cleff = cleff * cdmbgwd(2) +!-------------------------------------------------------------------- +! END SCALE-ADPTIVE PARAMETER SECTION +!-------------------------------------------------------------------- +! +!---- constants +! + rcl = 1. + rcs = sqrt(rcl) + cs = 1. / sqrt(rcl) + csg = cs * g + lcap = km + lcapp1 = lcap + 1 + fdir = mdir / (2.0*pi) + + 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 this module + else + xland(i)=2.0 + endif + RDXZB(i) = 0.0 + enddo + +!--- calculate scale-aware tapering factors +do i=1,im + if ( dx(i) .ge. dxmax_ls ) then + ls_taper(i) = 1. + else + if ( dx(i) .le. dxmin_ls) then + ls_taper(i) = 0. + else + ls_taper(i) = 0.5 * ( SIN(pi*(dx(i)-0.5*(dxmax_ls+dxmin_ls))/ & + (dxmax_ls-dxmin_ls)) + 1. ) + endif + endif +enddo + +do i=1,im + if ( dx(i) .ge. dxmax_ss ) then + ss_taper(i) = 1. + else + if ( dx(i) .le. dxmin_ss) then + ss_taper(i) = 0. + else + ss_taper(i) = dxmax_ss * (1. - dxmin_ss/dx(i))/(dxmax_ss-dxmin_ss) + endif + 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 + delx = dx(i) + dely = dx(i) + dxy4(i,1) = delx + dxy4(i,2) = dely + dxy4(i,3) = sqrt(delx*delx + dely*dely) + dxy4(i,4) = dxy4(i,3) + dxy4p(i,1) = dxy4(i,2) + dxy4p(i,2) = dxy4(i,1) + dxy4p(i,3) = dxy4(i,4) + dxy4p(i,4) = dxy4(i,3) +enddo +! +!-----initialize arrays +! + dtaux = 0.0 + dtauy = 0.0 + do i = its,im + klowtop(i) = 0 + kbl(i) = 0 + enddo +! + do i = its,im + xn(i) = 0.0 + yn(i) = 0.0 + ubar (i) = 0.0 + vbar (i) = 0.0 + roll (i) = 0.0 + taub (i) = 0.0 + oa(i) = 0.0 + ol(i) = 0.0 + oass(i) = 0.0 + olss(i) = 0.0 + ulow (i) = 0.0 + dtfac(i) = 1.0 + rstoch(i) = 0.0 + ldrag(i) = .false. + icrilv(i) = .false. + flag(i) = .true. + enddo + + do k = kts,km + do i = its,im + usqj(i,k) = 0.0 + bnv2(i,k) = 0.0 + vtj(i,k) = 0.0 + vtk(i,k) = 0.0 + taup(i,k) = 0.0 + taud_ls(i,k) = 0.0 + taud_bl(i,k) = 0.0 + dtaux2d(i,k) = 0.0 + dtauy2d(i,k) = 0.0 + enddo + enddo +! + if ( (gwd_opt == 33).or.(gwd_opt == 22) ) then + do i = its,im + dusfc_ls(i) = 0.0 + dvsfc_ls(i) = 0.0 + dusfc_bl(i) = 0.0 + dvsfc_bl(i) = 0.0 + dusfc_ss(i) = 0.0 + dvsfc_ss(i) = 0.0 + dusfc_fd(i) = 0.0 + dvsfc_fd(i) = 0.0 + enddo + do k = kts,km + do i = its,im + dtaux2d_ls(i,k)= 0.0 + dtauy2d_ls(i,k)= 0.0 + dtaux2d_bl(i,k)= 0.0 + dtauy2d_bl(i,k)= 0.0 + dtaux2d_ss(i,k)= 0.0 + dtauy2d_ss(i,k)= 0.0 + dtaux2d_fd(i,k)= 0.0 + dtauy2d_fd(i,k)= 0.0 + enddo + enddo + endif + + do i = its,im + taup(i,km+1) = 0.0 + xlinv(i) = 1.0/xl + dusfc(i) = 0.0 + dvsfc(i) = 0.0 + enddo +! +! initialize array for flow-blocking drag +! + taufb(1:im,1:km+1) = 0.0 + komax(1:im) = 0 +! + do k = kts,km + do i = its,im + vtj(i,k) = t1(i,k) * (1.+fv*q1(i,k)) + vtk(i,k) = vtj(i,k) / prslk(i,k) + ro(i,k) = 1./rd * prsl(i,k) / vtj(i,k) ! density kg/m**3 + enddo + enddo +! +! calculate mid-layer height (zl), interface height (zq), and layer depth (dz2). +! + !zq=0. + do k = kts,km + do i = its,im + !zq(i,k+1) = PHII(i,k+1)*g_inv + !dz2(i,k) = (PHII(i,k+1)-PHII(i,k))*g_inv + zl(i,k) = PHIL(i,k)*g_inv + enddo + enddo +! +! determine reference level: maximum of 2*var and pbl heights +! + do i = its,im + zlowtop(i) = 2. * var_stoch(i) + enddo +! + do i = its,im + kloop1(i) = .true. + enddo +! + do k = kts+1,km + do i = its,im + if(kloop1(i).and.zl(i,k)-zl(i,1).ge.zlowtop(i)) then + klowtop(i) = k+1 + kloop1(i) = .false. + endif + enddo + enddo +! + do i = its,im + kbl(i) = max(kpbl(i), klowtop(i)) + kbl(i) = max(min(kbl(i),kpblmax),kpblmin) + enddo +! +! determine the level of maximum orographic height +! + ! komax(:) = kbl(:) + komax(:) = klowtop(:) - 1 ! modification by NOAA/GSD March 2018 +! + do i = its,im + delks(i) = 1.0 / (prsi(i,1) - prsi(i,kbl(i))) + delks1(i) = 1.0 / (prsl(i,1) - prsl(i,kbl(i))) + enddo +! +! compute low level averages within pbl +! + do k = kts,kpblmax + do i = its,im + if (k.lt.kbl(i)) then + rcsks = rcs * del(i,k) * delks(i) + rdelks = del(i,k) * delks(i) + ubar(i) = ubar(i) + rcsks * u1(i,k) ! pbl u mean + vbar(i) = vbar(i) + rcsks * v1(i,k) ! pbl v mean + roll(i) = roll(i) + rdelks * ro(i,k) ! ro mean + endif + enddo + enddo +! +! figure out low-level horizontal wind direction +! +! nwd 1 2 3 4 5 6 7 8 +! wd w s sw nw e n ne se +! + do i = its,im + wdir = atan2(ubar(i),vbar(i)) + pi + idir = mod(nint(fdir*wdir),mdir) + 1 + nwd = nwdir(idir) + oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(i,mod(nwd-1,4)+1) + ol(i) = ol4(i,mod(nwd-1,4)+1) + ! Repeat for small-scale gwd + oass(i) = (1-2*int( (nwd-1)/4 )) * oa4ss(i,mod(nwd-1,4)+1) + olss(i) = ol4ss(i,mod(nwd-1,4)+1) + +! +!----- compute orographic width along (ol) and perpendicular (olp) +!----- the direction of wind +! + ol4p(1) = ol4(i,2) + ol4p(2) = ol4(i,1) + ol4p(3) = ol4(i,4) + ol4p(4) = ol4(i,3) + olp(i) = ol4p(mod(nwd-1,4)+1) +! +!----- compute orographic direction (horizontal orographic aspect ratio) +! + od(i) = olp(i)/max(ol(i),olmin) + od(i) = min(od(i),odmax) + od(i) = max(od(i),odmin) +! +!----- compute length of grid in the along(dxy) and cross(dxyp) wind directions +! + dxy(i) = dxy4(i,MOD(nwd-1,4)+1) + dxyp(i) = dxy4p(i,MOD(nwd-1,4)+1) + enddo +! +! END INITIALIZATION; BEGIN GWD CALCULATIONS: +! +IF ( (do_gsl_drag_ls_bl).and. & + ((gwd_opt_ls .EQ. 1).or.(gwd_opt_bl .EQ. 1)) ) then + + do i=its,im + + if ( ls_taper(i).GT.1.E-02 ) then + +! +!--- saving richardson number in usqj for migwdi +! + do k = kts,km-1 + ti = 2.0 / (t1(i,k)+t1(i,k+1)) + rdz = 1./(zl(i,k+1) - zl(i,k)) + tem1 = u1(i,k) - u1(i,k+1) + tem2 = v1(i,k) - v1(i,k+1) + dw2 = rcl*(tem1*tem1 + tem2*tem2) + shr2 = max(dw2,dw2min) * rdz * rdz + bvf2 = g*(g/cp+rdz*(vtj(i,k+1)-vtj(i,k))) * ti + usqj(i,k) = max(bvf2/shr2,rimin) + bnv2(i,k) = 2.0*g*rdz*(vtk(i,k+1)-vtk(i,k))/(vtk(i,k+1)+vtk(i,k)) + bnv2(i,k) = max( bnv2(i,k), bnv2min ) + enddo +! +!----compute the "low level" or 1/3 wind magnitude (m/s) +! + ulow(i) = max(sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i)), 1.0) + rulow(i) = 1./ulow(i) +! + do k = kts,km-1 + velco(i,k) = (0.5*rcs) * ((u1(i,k)+u1(i,k+1)) * ubar(i) & + + (v1(i,k)+v1(i,k+1)) * vbar(i)) + velco(i,k) = velco(i,k) * rulow(i) + if ((velco(i,k).lt.veleps) .and. (velco(i,k).gt.0.)) then + velco(i,k) = veleps + endif + enddo +! +! no drag when critical level in the base layer +! + ldrag(i) = velco(i,1).le.0. +! +! no drag when velco.lt.0 +! + do k = kpblmin,kpblmax + if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. velco(i,k).le.0. + enddo +! +! no drag when bnv2.lt.0 +! + do k = kts,kpblmax + if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. bnv2(i,k).lt.0. + enddo +! +!-----the low level weighted average ri is stored in usqj(1,1; im) +!-----the low level weighted average n**2 is stored in bnv2(1,1; im) +!---- this is called bnvl2 in phys_gwd_alpert_sub not bnv2 +!---- rdelks (del(k)/delks) vert ave factor so we can * instead of / +! + wtkbj = (prsl(i,1)-prsl(i,2)) * delks1(i) + bnv2(i,1) = wtkbj * bnv2(i,1) + usqj(i,1) = wtkbj * usqj(i,1) +! + do k = kpblmin,kpblmax + if (k .lt. kbl(i)) then + rdelks = (prsl(i,k)-prsl(i,k+1)) * delks1(i) + bnv2(i,1) = bnv2(i,1) + bnv2(i,k) * rdelks + usqj(i,1) = usqj(i,1) + usqj(i,k) * rdelks + endif + enddo +! + 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_stoch(i) .le. 0.0 +! +! set all ri low level values to the low level value +! + do k = kpblmin,kpblmax + if (k .lt. kbl(i)) usqj(i,k) = usqj(i,1) + enddo +! + if (.not.ldrag(i)) then + bnv(i) = sqrt( bnv2(i,1) ) + 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) + endif +! +! compute the base level stress and store it in taub +! calculate enhancement factor, number of mountains & aspect +! ratio const. use simplified relationship between standard +! deviation & critical hgt + + if (.not. ldrag(i)) then + efact = (oa(i) + 2.) ** (ce*fr(i)/frc) + efact = min( max(efact,efmin), efmax ) +!!!!!!! cleff (effective grid length) is highly tunable parameter +!!!!!!! the bigger (smaller) value produce weaker (stronger) wave drag +!WRF cleff = sqrt(dxy(i)**2. + dxyp(i)**2.) +!WRF cleff = 3. * max(dx(i),cleff) + coefm(i) = (1. + ol(i)) ** (oa(i)+1.) +!WRF xlinv(i) = coefm(i) / cleff + xlinv(i) = coefm(i) * cleff + tem = fr(i) * fr(i) * oc1(i) + gfobnv = gmax * tem / ((tem + cg)*bnv(i)) + if ( gwd_opt_ls .NE. 0 ) then + taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i) & + * ulow(i) * gfobnv * efact + else ! We've gotten what we need for the blocking scheme + taub(i) = 0.0 + end if + else + taub(i) = 0.0 + xn(i) = 0.0 + yn(i) = 0.0 + endif + + endif ! (ls_taper(i).GT.1.E-02) + + enddo ! do i=its,im + +ENDIF ! (do_gsl_drag_ls_bl).and.((gwd_opt_ls .EQ. 1).or.(gwd_opt_bl .EQ. 1)) + +!========================================================= +! add small-scale wavedrag for stable boundary layer +!========================================================= + XNBV=0. + tauwavex0=0. + tauwavey0=0. + density=1.2 + utendwave=0. + vtendwave=0. +! +IF ( do_gsl_drag_ss ) THEN + + do i=its,im + + if ( ss_taper(i).GT.1.E-02 ) then + ! + ! calculating potential temperature + ! + do k = kts,km + thx(i,k) = t1(i,k)/prslk(i,k) + enddo + ! + do k = kts,km + tvcon = (1.+fv*q1(i,k)) + thvx(i,k) = thx(i,k)*tvcon + enddo + + hpbl2 = hpbl(i)+10. + kpbl2 = kpbl(i) + !kvar = MIN(kpbl, k-level of var) + kvar = 1 + do k=kts+1,MAX(kpbl(i),kts+1) +! IF (zl(i,k)>2.*var(i) .or. zl(i,k)>2*varmax) then + IF (zl(i,k)>300.) then + kpbl2 = k + IF (k == kpbl(i)) then + hpbl2 = hpbl(i)+10. + ELSE + hpbl2 = zl(i,k)+10. + ENDIF + exit + ENDIF + enddo + 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) +! cleff_ss = 10. * max(dxmax_ss,cleff_ss) + cleff_ss = 0.1 * max(dxmax_ss,cleff_ss) ! WRF +! cleff_ss = 0.1 * 12000. + coefm_ss(i) = (1. + olss(i)) ** (oass(i)+1.) + xlinv(i) = coefm_ss(i) / cleff_ss + !govrth(i)=g/(0.5*(thvx(i,kpbl(i))+thvx(i,kts))) + govrth(i)=g/(0.5*(thvx(i,kpbl2)+thvx(i,kts))) + !XNBV=sqrt(govrth(i)*(thvx(i,kpbl(i))-thvx(i,kts))/hpbl(i)) + XNBV=sqrt(govrth(i)*(thvx(i,kpbl2)-thvx(i,kts))/hpbl2) +! + !if(abs(XNBV/u1(i,kpbl(i))).gt.xlinv(i))then + if(abs(XNBV/u1(i,kpbl2)).gt.xlinv(i))then + !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_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) + tauwavex0=tauwavex0*ss_taper(i) + else + tauwavex0=0. + endif +! + !if(abs(XNBV/v1(i,kpbl(i))).gt.xlinv(i))then + if(abs(XNBV/v1(i,kpbl2)).gt.xlinv(i))then + !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_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) + tauwavey0=tauwavey0*ss_taper(i) + else + tauwavey0=0. + endif + + do k=kts,kpbl(i) !MIN(kpbl2+1,km-1) +!original + !utendwave(i,k)=-1.*tauwavex0*2.*max((1.-zl(i,k)/hpbl(i)),0.)/hpbl(i) + !vtendwave(i,k)=-1.*tauwavey0*2.*max((1.-zl(i,k)/hpbl(i)),0.)/hpbl(i) +!new + utendwave(i,k)=-1.*tauwavex0*2.*max((1.-zl(i,k)/hpbl2),0.)/hpbl2 + vtendwave(i,k)=-1.*tauwavey0*2.*max((1.-zl(i,k)/hpbl2),0.)/hpbl2 +!mod-to be used in HRRRv3/RAPv4 + !utendwave(i,k)=-1.*tauwavex0 * max((1.-zl(i,k)/hpbl2),0.)**2 + !vtendwave(i,k)=-1.*tauwavey0 * max((1.-zl(i,k)/hpbl2),0.)**2 + enddo + endif + endif + + do k = kts,km + dudt(i,k) = dudt(i,k) + utendwave(i,k) + dvdt(i,k) = dvdt(i,k) + vtendwave(i,k) + dusfc(i) = dusfc(i) + utendwave(i,k) * del(i,k) + dvsfc(i) = dvsfc(i) + vtendwave(i,k) * del(i,k) + enddo + if(udtend>0) then + dtend(i,kts:km,udtend) = dtend(i,kts:km,udtend) + utendwave(i,kts:km)*deltim + endif + if(vdtend>0) then + dtend(i,kts:km,vdtend) = dtend(i,kts:km,vdtend) + vtendwave(i,kts:km)*deltim + endif + if ( (gwd_opt == 33).or.(gwd_opt == 22) ) then + do k = kts,km + dusfc_ss(i) = dusfc_ss(i) + utendwave(i,k) * del(i,k) + dvsfc_ss(i) = dvsfc_ss(i) + vtendwave(i,k) * del(i,k) + dtaux2d_ss(i,k) = utendwave(i,k) + dtauy2d_ss(i,k) = vtendwave(i,k) + enddo + endif + + endif ! if (ss_taper(i).GT.1.E-02) + + enddo ! i=its,im + +ENDIF ! if (do_gsl_drag_ss) + +!================================================================ +! Topographic Form Drag from Beljaars et al. (2004, QJRMS, equ. 16): +!================================================================ +IF ( do_gsl_drag_tofd ) THEN + + do i=its,im + + if ( ss_taper(i).GT.1.E-02 ) then + + utendform=0. + vtendform=0. + + 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_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_stoch(i),hpbl(i)) + H_efold = min(H_efold,1500.) + DO k=kts,km + wsp=SQRT(u1(i,k)**2 + v1(i,k)**2) + ! alpha*beta*Cmd*Ccorr*2.109 = 12.*1.*0.005*0.6*2.109 = 0.0759 + var_temp = 0.0759*EXP(-(zl(i,k)/H_efold)**1.5)*a2* & + zl(i,k)**(-1.2)*ss_taper(i) ! this is greater than zero + ! Note: This is a semi-implicit treatment of the time differencing + ! per Beljaars et al. (2004, QJRMS) + utendform(i,k) = - var_temp*wsp*u1(i,k)/(1. + var_temp*deltim*wsp) + vtendform(i,k) = - var_temp*wsp*v1(i,k)/(1. + var_temp*deltim*wsp) + !IF(zl(i,k) > 4000.) exit + ENDDO + ENDIF + + do k = kts,km + dudt(i,k) = dudt(i,k) + utendform(i,k) + dvdt(i,k) = dvdt(i,k) + vtendform(i,k) + dusfc(i) = dusfc(i) + utendform(i,k) * del(i,k) + dvsfc(i) = dvsfc(i) + vtendform(i,k) * del(i,k) + enddo + if(udtend>0) then + dtend(i,kts:km,udtend) = dtend(i,kts:km,udtend) + utendform(i,kts:km)*deltim + endif + if(vdtend>0) then + dtend(i,kts:km,vdtend) = dtend(i,kts:km,vdtend) + vtendform(i,kts:km)*deltim + endif + if ( (gwd_opt == 33).or.(gwd_opt == 22) ) then + do k = kts,km + dtaux2d_fd(i,k) = utendform(i,k) + dtauy2d_fd(i,k) = vtendform(i,k) + dusfc_fd(i) = dusfc_fd(i) + utendform(i,k) * del(i,k) + dvsfc_fd(i) = dvsfc_fd(i) + vtendform(i,k) * del(i,k) + enddo + endif + + endif ! if (ss_taper(i).GT.1.E-02) + + enddo ! i=its,im + +ENDIF ! if (do_gsl_drag_tofd) +!======================================================= +! More for the large-scale gwd component +IF ( (do_gsl_drag_ls_bl).and.(gwd_opt_ls .EQ. 1) ) THEN + + do i=its,im + + if ( ls_taper(i).GT.1.E-02 ) then + +! +! now compute vertical structure of the stress. + do k = kts,kpblmax + if (k .le. kbl(i)) taup(i,k) = taub(i) + enddo +! + do k = kpblmin, km-1 ! vertical level k loop! + kp1 = k + 1 +! +! unstablelayer if ri < ric +! unstable layer if upper air vel comp along surf vel <=0 (crit lay) +! at (u-c)=0. crit layer exists and bit vector should be set (.le.) +! + if (k .ge. kbl(i)) then + icrilv(i) = icrilv(i) .or. ( usqj(i,k) .lt. ric) & + .or. (velco(i,k) .le. 0.0) + brvf(i) = max(bnv2(i,k),bnv2min) ! brunt-vaisala frequency squared + brvf(i) = sqrt(brvf(i)) ! brunt-vaisala frequency + endif +! + if (k .ge. kbl(i) .and. (.not. ldrag(i))) then + if (.not.icrilv(i) .and. taup(i,k) .gt. 0.0 ) then + temv = 1.0 / velco(i,k) + tem1 = coefm(i)/dxy(i)*(ro(i,kp1)+ro(i,k))*brvf(i)* & + velco(i,k)*0.5 + hd = sqrt(taup(i,k) / tem1) + fro = brvf(i) * hd * temv +! +! rim is the minimum-richardson number by shutts (1985) + tem2 = sqrt(usqj(i,k)) + tem = 1. + tem2 * fro + rim = usqj(i,k) * (1.-fro) / (tem * tem) +! +! check stability to employ the 'saturation hypothesis' +! of lindzen (1981) except at tropospheric downstream regions +! + if (rim .le. ric) then ! saturation hypothesis! + if ((oa(i) .le. 0.).or.(kp1 .ge. kpblmin )) then + temc = 2.0 + 1.0 / tem2 + hd = velco(i,k) * (2.*sqrt(temc)-temc) / brvf(i) + taup(i,kp1) = tem1 * hd * hd + endif + else ! no wavebreaking! + taup(i,kp1) = taup(i,k) + endif + endif + endif + enddo +! + if(lcap.lt.km) then + do klcap = lcapp1,km + taup(i,klcap) = prsi(i,klcap) / prsi(i,lcap) * taup(i,lcap) + enddo + endif + + endif ! if ( ls_taper(i).GT.1.E-02 ) + + enddo ! do i=its,im + +ENDIF ! (do_gsl_drag_ls_bl).and.(gwd_opt_ls .EQ. 1) +!=============================================================== +!COMPUTE BLOCKING COMPONENT +!=============================================================== +IF ( (do_gsl_drag_ls_bl) .and. (gwd_opt_bl .EQ. 1) ) THEN + + do i=its,im + + if ( ls_taper(i).GT.1.E-02 ) then + + if (.not.ldrag(i)) then +! +!------- determine the height of flow-blocking layer +! + kblk = 0 + pe = 0.0 + do k = km, kpblmin, -1 + if(kblk.eq.0 .and. k.le.komax(i)) then + pe = pe + bnv2(i,k)*(zl(i,komax(i))-zl(i,k))* & + del(i,k)/g/ro(i,k) + ke = 0.5*((rcs*u1(i,k))**2.+(rcs*v1(i,k))**2.) +! +!---------- apply flow-blocking drag when pe >= ke +! + if(pe.ge.ke) then + kblk = k + kblk = min(kblk,kbl(i)) + zblk = zl(i,kblk)-zl(i,kts) + RDXZB(i) = real(k,kind=kind_phys) + endif + endif + enddo + if(kblk.ne.0) then +! +!--------- compute flow-blocking stress +! + cd = max(2.0-1.0/od(i),0.0) + taufb(i,kts) = 0.5 * roll(i) * coefm(i) / & + max(dxmax_ls,dxy(i))**2 * cd * dxyp(i) * & + olp(i) * zblk * ulow(i)**2 + tautem = taufb(i,kts)/float(kblk-kts) + do k = kts+1, kblk + taufb(i,k) = taufb(i,k-1) - tautem + enddo +! +!----------sum orographic GW stress and flow-blocking stress +! + ! taup(i,:) = taup(i,:) + taufb(i,:) ! Keep taup and taufb separate for now + endif + + endif ! if (.not.ldrag(i)) + + endif ! if ( ls_taper(i).GT.1.E-02 ) + + enddo ! do i=its,im + +ENDIF ! IF ( (do_gsl_drag_ls_bl) .and. (gwd_opt_bl .EQ. 1) ) +!=========================================================== +IF ( (do_gsl_drag_ls_bl) .and. & + (gwd_opt_ls .EQ. 1 .OR. gwd_opt_bl .EQ. 1) ) THEN + + do i=its,im + + if ( ls_taper(i) .GT. 1.E-02 ) then + +! +! calculate - (g)*d(tau)/d(pressure) and deceleration terms dtaux, dtauy +! + do k = kts,km + taud_ls(i,k) = 1. * (taup(i,k+1) - taup(i,k)) * csg / del(i,k) + taud_bl(i,k) = 1. * (taufb(i,k+1) - taufb(i,k)) * csg / del(i,k) + enddo +! +! limit de-acceleration (momentum deposition ) at top to 1/2 value +! the idea is some stuff must go out the 'top' + do klcap = lcap,km + taud_ls(i,klcap) = taud_ls(i,klcap) * factop + taud_bl(i,klcap) = taud_bl(i,klcap) * factop + enddo +! +! if the gravity wave drag would force a critical line +! in the lower ksmm1 layers during the next deltim timestep, +! then only apply drag until that critical line is reached. +! + do k = kts,kpblmax-1 + if (k .le. kbl(i)) then + if ((taud_ls(i,k)+taud_bl(i,k)).ne.0.) & + dtfac(i) = min(dtfac(i),abs(velco(i,k) & + /(deltim*rcs*(taud_ls(i,k)+taud_bl(i,k))))) + endif + enddo +! + do k = kts,km + taud_ls(i,k) = taud_ls(i,k)*dtfac(i)* ls_taper(i) *(1.-rstoch(i)) + taud_bl(i,k) = taud_bl(i,k)*dtfac(i)* ls_taper(i) *(1.-rstoch(i)) + + dtaux = taud_ls(i,k) * xn(i) + dtauy = taud_ls(i,k) * yn(i) + dtauxb = taud_bl(i,k) * xn(i) + dtauyb = taud_bl(i,k) * yn(i) + + !add blocking and large-scale contributions to tendencies + dudt(i,k) = dtaux + dtauxb + dudt(i,k) + dvdt(i,k) = dtauy + dtauyb + dvdt(i,k) + + if ( gsd_diss_ht_opt .EQ. 1 ) then + ! Calculate dissipation heating + ! Initial kinetic energy (at t0-dt) + eng0 = 0.5*( (rcs*u1(i,k))**2. + (rcs*v1(i,k))**2. ) + ! Kinetic energy after wave-breaking/flow-blocking + eng1 = 0.5*( (rcs*(u1(i,k)+(dtaux+dtauxb)*deltim))**2 + & + (rcs*(v1(i,k)+(dtauy+dtauyb)*deltim))**2 ) + ! Modify theta tendency + dtdt(i,k) = dtdt(i,k) + max((eng0-eng1),0.0)/cp/deltim + if ( Tdtend>0 ) then + dtend(i,k,Tdtend) = dtend(i,k,Tdtend) + max((eng0-eng1),0.0)/cp + endif + endif + + dusfc(i) = dusfc(i) + taud_ls(i,k)*xn(i)*del(i,k) + & + taud_bl(i,k)*xn(i)*del(i,k) + dvsfc(i) = dvsfc(i) + taud_ls(i,k)*yn(i)*del(i,k) + & + taud_bl(i,k)*yn(i)*del(i,k) + if(udtend>0) then + dtend(i,k,udtend) = dtend(i,k,udtend) + (taud_ls(i,k) * & + xn(i) + taud_bl(i,k) * xn(i)) * deltim + endif + if(vdtend>0) then + dtend(i,k,vdtend) = dtend(i,k,vdtend) + (taud_ls(i,k) * & + yn(i) + taud_bl(i,k) * yn(i)) * deltim + endif + + enddo + + ! Finalize dusfc and dvsfc diagnostics + dusfc(i) = (-1./g*rcs) * dusfc(i) + dvsfc(i) = (-1./g*rcs) * dvsfc(i) + + if ( (gwd_opt == 33).or.(gwd_opt == 22) ) then + do k = kts,km + dtaux2d_ls(i,k) = taud_ls(i,k) * xn(i) + dtauy2d_ls(i,k) = taud_ls(i,k) * yn(i) + dtaux2d_bl(i,k) = taud_bl(i,k) * xn(i) + dtauy2d_bl(i,k) = taud_bl(i,k) * yn(i) + dusfc_ls(i) = dusfc_ls(i) + dtaux2d_ls(i,k) * del(i,k) + dvsfc_ls(i) = dvsfc_ls(i) + dtauy2d_ls(i,k) * del(i,k) + dusfc_bl(i) = dusfc_bl(i) + dtaux2d_bl(i,k) * del(i,k) + dvsfc_bl(i) = dvsfc_bl(i) + dtauy2d_bl(i,k) * del(i,k) + enddo + endif + + endif ! if ( ls_taper(i) .GT. 1.E-02 ) + + enddo ! do i=its,im + +ENDIF ! (do_gsl_drag_ls_bl).and.(gwd_opt_ls.EQ.1 .OR. gwd_opt_bl.EQ.1) + +if ( (gwd_opt == 33).or.(gwd_opt == 22) ) then + ! Finalize dusfc and dvsfc diagnostics + do i = its,im + dusfc_ls(i) = (-1./g*rcs) * dusfc_ls(i) + dvsfc_ls(i) = (-1./g*rcs) * dvsfc_ls(i) + dusfc_bl(i) = (-1./g*rcs) * dusfc_bl(i) + dvsfc_bl(i) = (-1./g*rcs) * dvsfc_bl(i) + dusfc_ss(i) = (-1./g*rcs) * dusfc_ss(i) + dvsfc_ss(i) = (-1./g*rcs) * dvsfc_ss(i) + dusfc_fd(i) = (-1./g*rcs) * dusfc_fd(i) + dvsfc_fd(i) = (-1./g*rcs) * dvsfc_fd(i) + enddo +endif +! + return + end subroutine drag_suite_run +!------------------------------------------------------------------- +! + subroutine drag_suite_finalize() + end subroutine drag_suite_finalize + + end module drag_suite From 96502941f925197d6278caa3419dcdfe0bd96553 Mon Sep 17 00:00:00 2001 From: jeff beck Date: Thu, 20 Jan 2022 04:34:35 +0000 Subject: [PATCH 19/28] Update dimensions of SPP fields --- physics/drag_suite.F90 | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/physics/drag_suite.F90 b/physics/drag_suite.F90 index 31fb4fd50..0f843387f 100644 --- a/physics/drag_suite.F90 +++ b/physics/drag_suite.F90 @@ -606,14 +606,16 @@ subroutine drag_suite_run( & 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 = varmax_ss + varmax_ss*0.75*spp_wts_gwd(i,1) - varmax_fd_stoch = varmax_fd + varmax_fd*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 - var_stoch = var - varss_stoch = varss - varmax_ss_stoch = varmax_ss - varmax_fd_stoch = varmax_fd + 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 @@ -1001,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_stoch(i),varmax_ss_stoch) + & - MAX(0.,beta_ss*(varss_stoch(i)-varmax_ss_stoch)) + 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) @@ -1016,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_stoch(i),varmax_ss_stoch) + & - MAX(0.,beta_ss*(varss_stoch(i)-varmax_ss_stoch)) + 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) From 63c5f96260aca15fc3757f358437a317b48243e0 Mon Sep 17 00:00:00 2001 From: jeff beck Date: Thu, 20 Jan 2022 04:35:55 +0000 Subject: [PATCH 20/28] Update dimensions for two SPP fields --- physics/drag_suite.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/drag_suite.F90 b/physics/drag_suite.F90 index 0f843387f..7fea98b13 100644 --- a/physics/drag_suite.F90 +++ b/physics/drag_suite.F90 @@ -610,7 +610,7 @@ subroutine drag_suite_run( & varmax_fd_stoch(i) = varmax_fd + varmax_fd*0.75*spp_wts_gwd(i,1) enddo else - do i=its,im + do i = its,im var_stoch(i) = var(i) varss_stoch(i) = varss(i) varmax_ss_stoch(i) = varmax_ss @@ -1083,8 +1083,8 @@ 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_stoch(i),varmax_fd_stoch) + & - MAX(0.,beta_fd*(varss_stoch(i)-varmax_fd_stoch)) + 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 From 070a9cdea3980de5d936786d41d260efcc3644a5 Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Fri, 11 Feb 2022 14:24:25 -0700 Subject: [PATCH 21/28] change units of surface_stochastic_weights_from_coupled_process from none to 1 to match change in fv3atm --- physics/GFS_rrtmg_pre.meta | 2 +- physics/GFS_surface_generic.meta | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/physics/GFS_rrtmg_pre.meta b/physics/GFS_rrtmg_pre.meta index 37cd8d17a..1eac8a571 100644 --- a/physics/GFS_rrtmg_pre.meta +++ b/physics/GFS_rrtmg_pre.meta @@ -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 diff --git a/physics/GFS_surface_generic.meta b/physics/GFS_surface_generic.meta index 6ad2953a6..28c88c5ea 100644 --- a/physics/GFS_surface_generic.meta +++ b/physics/GFS_surface_generic.meta @@ -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 From 308b0b38347d5460351f91b8648646e15a90f220 Mon Sep 17 00:00:00 2001 From: jeff beck Date: Wed, 16 Feb 2022 16:45:21 +0000 Subject: [PATCH 22/28] Remove misplaced file. --- drag_suite.F90 | 1381 ------------------------------------------------ 1 file changed, 1381 deletions(-) delete mode 100755 drag_suite.F90 diff --git a/drag_suite.F90 b/drag_suite.F90 deleted file mode 100755 index 7fea98b13..000000000 --- a/drag_suite.F90 +++ /dev/null @@ -1,1381 +0,0 @@ -!> \File drag_suite.F90 -!! This file is the parameterization of orographic gravity wave -!! drag, mountain blocking, and form drag. - -!> This module contains the CCPP-compliant orographic gravity wave dray scheme. - module drag_suite - - contains - - subroutine drag_suite_init(gwd_opt, errmsg, errflg) - - integer, intent(in) :: gwd_opt - character(len=*), intent(out) :: errmsg - integer, intent(out) :: errflg - - - ! Initialize CCPP error handling variables - errmsg = '' - errflg = 0 - - ! Consistency checks - if (gwd_opt/=3 .and. gwd_opt/=33) then - write(errmsg,'(*(a))') "Logic error: namelist choice of gravity wave & - & drag is different from drag_suite scheme" - errflg = 1 - return - end if - end subroutine drag_suite_init - -! \defgroup GFS_ogwd GFS Orographic Gravity Wave Drag -!> \defgroup gfs_drag_suite GFS drag_suite Main -!! \brief This subroutine includes orographic gravity wave drag, mountain -!! blocking, and form drag. -!! -!> The time tendencies of zonal and meridional wind are altered to -!! include the effect of mountain induced gravity wave drag from -!! subgrid scale orography including convective breaking, shear -!! breaking and the presence of critical levels. -!! -!> \section arg_table_drag_suite_run Argument Table -!! \htmlinclude drag_suite_run.html -!! -!> \section gen_drag_suite GFS Orographic GWD Scheme General Algorithm -!! -# Calculate subgrid mountain blocking -!! -# Calculate orographic wave drag -!! -!! The NWP model gravity wave drag (GWD) scheme in the GFS has two -!! main components: how the surface stress is computed, and then how -!! that stress is distributed over a vertical column where it may -!! interact with the models momentum. Each of these depends on the -!! large scale environmental atmospheric state and assumptions about -!! the sub-grid scale processes. In Alpert GWD (1987) based on linear, -!! two-dimensional non-rotating, stably stratified flow over a mountain ridge, -!! sub-grid scale gravity wave motions are assumed which propagate away -!! from the mountain. Described in Alpert (1987), the flux measured over -!! a "low level" vertically averaged layer, in the atmosphere defines a base -!! level flux. "Low level" was taken to be the first 1/3 of the troposphere -!! in the 1987 implementation. This choice was meant to encompass a thick -!! low layer for vertical averages of the environmental (large scale) flow -!! quantities. The vertical momentum flux or gravity wave stress in a -!! grid box due to a single mountain is given as in Pierrehumbert, (1987) (PH): -!! -!! \f$ \tau = \frac {\rho \: U^{3}\: G(F_{r})} {\Delta X \; N } \f$ -!! -!! emetic \f$ \Delta X \f$ is a grid increment, N is the Brunt Viasala frequency -!! -!! -!! \f$ N(\sigma) = \frac{-g \: \sigma \: -!! \frac{\partial\Theta}{\partial\sigma}}{\Theta \:R \:T} \f$ -!! -!! The environmental variables are calculated from a mass weighted vertical -!! average over a base layer. G(Fr) is a monotonically increasing -!! function of Froude number, -!! -!! \f$ F_{r} = \frac{N h^{'}}{U} \f$ -!! -!! where U is the wind speed calculated as a mass weighted vertical average in -!! the base layer, and h', is the vertical displacement caused by the orography -!! variance. An effective mountain length for the gravity wave processes, -!! -!! \f$ l^{*} = \frac{\Delta X}{m} \f$ -!! -!! where m is the number of mountains in a grid box, can then -!! be defined to obtain the form of the base level stress -!! -!! -!! \f$ \tau = \frac {\rho \: U^{3} \: G(F_{r})} {N \;l^{*}} \f$ -!! -!! giving the stress induced from the surface in a model grid box. -!! PH gives the form for the function G(Fr) as -!! -!! -!! \f$ G(F_{r}) = \bar{G}\frac{F^{2}_{r}}{F^{2}_{r}\: + \:a^{2}} \f$ -!! -!! Where \f$ \bar{G} \f$ is an order unity non-dimensional saturation -!! flux set to 1 and 'a' is a function of the mountain aspect ratio also -!!set to 1 in the 1987 implementation of the GFS GWD. Typical values of -!! U=10m/s, N=0.01 1/s, l*=100km, and a=1, gives a flux of 1 Pascal and -!! if this flux is made to go to zero linearly with height then the -!! decelerations would be about 10/m/s/day which is consistent with -!! observations in PH. -!! -!! -!! In Kim, Moorthi, Alpert's (1998, 2001) GWD currently in GFS operations, -!! the GWD scheme has the same physical basis as in Alpert (1987) with the addition -!! of enhancement factors for the amplitude, G, and mountain shape details -!! in G(Fr) to account for effects from the mountain blocking. A factor, -!! E m', is an enhancement factor on the stress in the Alpert '87 scheme. -!! The E ranges from no enhancement to an upper limit of 3, E=E(OA)[1-3], -!! and is a function of OA, the Orographic Asymmetry defined in KA (1995) as -!! -!! Orographic Asymmetry (OA) = \f$ \frac{ \bar{x} \; - \; -!! \sum\limits_{j=1}^{N_{b}} x_{j} \; n_{j} }{\sigma_{x}} \f$ -!! -!! where Nb is the total number of bottom blocks in the mountain barrier, -!! \f$ \sigma_{x} \f$ is the standard deviation of the horizontal distance defined by -!! -!! \f$ \sigma_{x} = \sqrt{ \frac{\sum\limits_{j=1}^{N_{b}} -!! \; (x_{j} \; - \; \bar{x} )^2}{N_{x}} } \f$ -!! -!! -!! where Nx is the number of grid intervals for the large scale domain being -!! considered. So the term, E(OA)m'/ \f$ \Delta X \f$ in Kim's scheme represents -!! a multiplier on G shown in Alpert's eq (1), where m' is the number of mountains -!! in a sub-grid scale box. Kim increased the complexity of m' making it a -!! function of the fractional area of the sub-grid mountain and the asymmetry -!! and convexity statistics which are found from running a gravity wave -!! model for a large number of cases: -!! -!! \f$ m^{'} = C_{m} \Delta X \left[ \frac{1 \; + \; -!! \sum\limits_{x} L_{h} }{\Delta X} \right]^{OA+1} \f$ -!! -!! Where, according to Kim, \f$ \sum \frac{L_{h}}{\Delta X} \f$ is -!! the fractional area covered by the subgrid-scale orography higher than -!! a critical height \f$ h_{c} = Fr_{c} U_{0}/N_{0} \f$ , over the -!! "low level" vertically averaged layer, for a grid box with the interval -!! \f$ \Delta X \f$. Each \f$ L_{n}\f$ is the width of a segment of -!! orography intersection at the critical height: -!! -!! \f$ Fr_{0} = \frac{N_{0} \; h^{'}}{U_{0}} \f$ -!! -!! \f$ G^{'}(OC,Fr_{0}) = \frac{Fr_{0}^{2}}{Fr_{0}^{2} \; + \; a^{2}} \f$ -!! -!! \f$ a^{2} = \frac{C_{G}}{OC} \f$ -!! -!! \f$ E(OA, Fr_{0}) = (OA \; + \; 2)^{\delta} \f$ and \f$ \delta -!! \; = \; \frac{C_{E} \; Fr_{0}}{Fr_{c}} \f$ where \f$ Fr_{c} \f$ -!! is as in Alpert. -!! -!! -!! This represents a closed scheme, somewhat empirical adjustments -!! to the original scheme to calculate the surface stress. -!! -!! Momentum is deposited by the sub-grid scale gravity waves break due -!! to the presence of convective mixing assumed to occur when the -!! minimum Richardson number: -!! -!! Orographic Convexity (OC) = \f$ \frac{ \sum\limits_{j=1}^{N_{x}} -!! \; (h_{j} \; - \; \bar{h})^4 }{N_{x} \;\sigma_{h}^4} \f$ , -!! and where \f$ \sigma_{h} = \sqrt{ \frac{\sum\limits_{j=1}^{N_{x}} -!! \; (h_{j} \; - \; \bar{h} )^2}{N_{x}} } \f$ -!! -!! This represents a closed scheme, somewhat empirical adjustments -!! to the original scheme to calculate the surface stress. -!! -!! Momentum is deposited by the sub-grid scale gravity waves break due -!! to the presence of convective mixing assumed to occur when -!! the minimum Richardson number: -!! -!! \f$ Ri_{m} = \frac{Ri(1 \; - \; Fr)}{(1 \; + \; \sqrt{Ri}Fr)^2} \f$ -!! -!! Is less than 1/4 Or if critical layers are encountered in a layer -!! the the momentum flux will vanish. The critical layer is defined -!! when the base layer wind becomes perpendicular to the environmental -!! wind. Otherwise, wave breaking occurs at a level where the amplification -!! of the wave causes the local Froude number or similarly a truncated -!! (first term of the) Scorer parameter, to be reduced below a critical -!! value by the saturation hypothesis (Lindzen,). This is done through -!! eq 1 which can be written as -!! -!! \f$ \tau = \rho U N k h^{'2} \f$ -!! -!! For small Froude number this is discretized in the vertical so at each -!! level the stress is reduced by ratio of the Froude or truncated Scorer -!! parameter, \f$ \frac{U^{2}}{N^{2}} = \frac{N \tau_{l-1}}{\rho U^{3} k} \f$ , -!! where the stress is from the layer below beginning with that found near -!! the surface. The respective change in momentum is applied in -!! that layer building up from below. -!! -!! An amplitude factor is part of the calibration of this scheme which is -!! a function of the model resolution and the vertical diffusion. This -!! is because the vertical diffusion and the GWD account encompass -!! similar physical processes. Thus, one needs to run the model over -!! and over for various amplitude factors for GWD and vertical diffusion. -!! -!! In addition, there is also mountain blocking from lift and frictional -!! forces. Improved integration between how the GWD is calculated and -!! the mountain blocking of wind flow around sub-grid scale orography -!! is underway at NCEP. The GFS already has convectively forced GWD -!! an independent process. The next step is to test -!! -!> \section det_drag_suite GFS Orographic GWD Scheme Detailed Algorithm -!> @{ - subroutine drag_suite_run( & - & IM,KM,dvdt,dudt,dtdt,U1,V1,T1,Q1,KPBL, & - & PRSI,DEL,PRSL,PRSLK,PHII,PHIL,DELTIM,KDT, & - & var,oc1,oa4,ol4, & - & varss,oc1ss,oa4ss,ol4ss, & - & THETA,SIGMA,GAMMA,ELVMAX, & - & dtaux2d_ls,dtauy2d_ls,dtaux2d_bl,dtauy2d_bl, & - & dtaux2d_ss,dtauy2d_ss,dtaux2d_fd,dtauy2d_fd, & - & dusfc,dvsfc, & - & dusfc_ls,dvsfc_ls,dusfc_bl,dvsfc_bl, & - & dusfc_ss,dvsfc_ss,dusfc_fd,dvsfc_fd, & - & slmsk,br1,hpbl, & - & g, cp, rd, rv, fv, pi, imx, cdmbgwd, me, master, & - & lprnt, ipr, rdxzb, dx, gwd_opt, & - & 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, & - & 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 <---------- -! -! ----- This code ----- -!begin WRF code - -! this code handles the time tendencies of u v due to the effect of mountain -! induced gravity wave drag from sub-grid scale orography. this routine -! not only treats the traditional upper-level wave breaking due to mountain -! variance (alpert 1988), but also the enhanced lower-tropospheric wave -! breaking due to mountain convexity and asymmetry (kim and arakawa 1995). -! thus, in addition to the terrain height data in a model grid box, -! additional 10-2d topographic statistics files are needed, including -! orographic standard deviation (var), convexity (oc1), asymmetry (oa4) -! and ol (ol4). these data sets are prepared based on the 30 sec usgs orography -! hong (1999). the current scheme was implmented as in hong et al.(2008) -! -! Originally coded by song-you hong and young-joon kim and implemented by song-you hong -! -! program history log: -! 2014-10-01 Hyun-Joo Choi (from KIAPS) flow-blocking drag of kim and doyle -! with blocked height by dividing streamline theory -! 2017-04-06 Joseph Olson (from Gert-Jan Steeneveld) added small-scale -! orographic grabity wave drag: -! 2017-09-15 Joseph Olson, with some bug fixes from Michael Toy: added the -! topographic form drag of Beljaars et al. (2004, QJRMS) -! Activation of each component is done by specifying the integer-parameters -! (defined below) to 0: inactive or 1: active -! gwd_opt_ls = 0 or 1: large-scale -! gwd_opt_bl = 0 or 1: blocking drag -! gwd_opt_ss = 0 or 1: small-scale gravity wave drag -! gwd_opt_fd = 0 or 1: topographic form drag -! 2017-09-25 Michael Toy (from NCEP GFS model) added dissipation heating -! gsd_diss_ht_opt = 0: dissipation heating off -! gsd_diss_ht_opt = 1: dissipation heating on -! 2020-08-25 Michael Toy changed logic control for drag component selection -! for CCPP. -! Namelist options: -! do_gsl_drag_ls_bl - logical flag for large-scale GWD + blocking -! do_gsl_drag_ss - logical flag for small-scale GWD -! do_gsl_drag_tofd - logical flag for turbulent form drag -! Compile-time options (same as before): -! gwd_opt_ls = 0 or 1: large-scale GWD -! gwd_opt_bl = 0 or 1: blocking drag -! -! References: -! Hong et al. (2008), wea. and forecasting -! Kim and Doyle (2005), Q. J. R. Meteor. Soc. -! Kim and Arakawa (1995), j. atmos. sci. -! Alpert et al. (1988), NWP conference. -! Hong (1999), NCEP office note 424. -! Steeneveld et al (2008), JAMC -! Tsiringakis et al. (2017), Q. J. R. Meteor. Soc. -! Beljaars et al. (2004), Q. J. R. Meteor. Soc. -! -! notice : comparible or lower resolution orography files than model resolution -! are desirable in preprocess (wps) to prevent weakening of the drag -!------------------------------------------------------------------------------- -! -! input -! dudt (im,km) non-lin tendency for u wind component -! dvdt (im,km) non-lin tendency for v wind component -! u1(im,km) zonal wind / sqrt(rcl) m/sec at t0-dt -! v1(im,km) meridional wind / sqrt(rcl) m/sec at t0-dt -! t1(im,km) temperature deg k at t0-dt -! q1(im,km) specific humidity at t0-dt -! deltim time step secs -! del(km) positive increment of pressure across layer (pa) -! KPBL(IM) is the index of the top layer of the PBL -! ipr & lprnt for diagnostics -! -! output -! dudt, dvdt wind tendency due to gwdo -! dTdt -! -!------------------------------------------------------------------------------- - -!end wrf code -!----------------------------------------------------------------------C -! USE -! ROUTINE IS CALLED FROM CCPP (AFTER CALLING PBL SCHEMES) -! -! PURPOSE -! USING THE GWD PARAMETERIZATIONS OF PS-GLAS AND PH- -! GFDL TECHNIQUE. THE TIME TENDENCIES OF U V -! ARE ALTERED TO INCLUDE THE EFFECT OF MOUNTAIN INDUCED -! GRAVITY WAVE DRAG FROM SUB-GRID SCALE OROGRAPHY INCLUDING -! CONVECTIVE BREAKING, SHEAR BREAKING AND THE PRESENCE OF -! CRITICAL LEVELS -! -! -! ******************************************************************** - USE MACHINE , ONLY : kind_phys - implicit none - - ! Interface variables - integer, intent(in) :: im, km, imx, kdt, ipr, me, master - integer, intent(in) :: gwd_opt - logical, intent(in) :: lprnt - integer, intent(in) :: KPBL(:) - real(kind=kind_phys), intent(in) :: deltim, G, CP, RD, RV, cdmbgwd(:) - real(kind=kind_phys), intent(inout) :: dtend(:,:,:) - logical, intent(in) :: ldiag3d - integer, intent(in) :: dtidx(:,:), index_of_temperature, & - & index_of_process_orographic_gwd, index_of_x_wind, index_of_y_wind - - integer :: kpblmax - integer, parameter :: ims=1, kms=1, its=1, kts=1 - real(kind=kind_phys), intent(in) :: fv, pi - real(kind=kind_phys) :: rcl, cdmb - real(kind=kind_phys) :: g_inv - - real(kind=kind_phys), intent(inout) :: & - & dudt(:,:),dvdt(:,:), & - & dtdt(:,:) - real(kind=kind_phys), intent(out) :: rdxzb(:) - real(kind=kind_phys), intent(in) :: & - & u1(:,:),v1(:,:), & - & t1(:,:),q1(:,:), & - & PHII(:,:),prsl(:,:), & - & prslk(:,:),PHIL(:,:) - real(kind=kind_phys), intent(in) :: prsi(:,:), & - & del(:,:) - real(kind=kind_phys), intent(in) :: var(:),oc1(:), & - & oa4(:,:),ol4(:,:), & - & dx(:) - real(kind=kind_phys), intent(in) :: varss(:),oc1ss(:), & - & oa4ss(:,:),ol4ss(:,:) - real(kind=kind_phys), intent(in) :: THETA(:),SIGMA(:), & - & GAMMA(:),ELVMAX(:) - -! added for small-scale orographic wave drag - real(kind=kind_phys), dimension(im,km) :: utendwave,vtendwave,thx,thvx - real(kind=kind_phys), intent(in) :: br1(:), & - & hpbl(:), & - & slmsk(:) - real(kind=kind_phys), dimension(im) :: govrth,xland - !real(kind=kind_phys), dimension(im,km) :: dz2 - real(kind=kind_phys) :: tauwavex0,tauwavey0, & - & XNBV,density,tvcon,hpbl2 - integer :: kpbl2,kvar - !real(kind=kind_phys), dimension(im,km+1) :: zq ! = PHII/g - 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: - real(kind=kind_phys), intent(out) :: & - & dusfc(:), dvsfc(:) -!Output (optional): - real(kind=kind_phys), intent(out) :: & - & dusfc_ls(:),dvsfc_ls(:), & - & dusfc_bl(:),dvsfc_bl(:), & - & dusfc_ss(:),dvsfc_ss(:), & - & dusfc_fd(:),dvsfc_fd(:) - real(kind=kind_phys), intent(out) :: & - & dtaux2d_ls(:,:),dtauy2d_ls(:,:), & - & dtaux2d_bl(:,:),dtauy2d_bl(:,:), & - & dtaux2d_ss(:,:),dtauy2d_ss(:,:), & - & dtaux2d_fd(:,:),dtauy2d_fd(:,:) - -!Misc arrays - real(kind=kind_phys), dimension(im,km) :: dtaux2d, dtauy2d - -!------------------------------------------------------------------------- -! Flags to regulate the activation of specific components of drag suite: -! Each component is tapered off automatically as a function of dx, so best to -! keep them activated (.true.). - logical, intent(in) :: & - do_gsl_drag_ls_bl, & ! large-scale gravity wave drag and blocking - do_gsl_drag_ss, & ! small-scale gravity wave drag (Steeneveld et al. 2008) - do_gsl_drag_tofd ! form drag (Beljaars et al. 2004, QJRMS) - -! Additional flags - integer, parameter :: & - gwd_opt_ls = 1, & ! large-scale gravity wave drag - gwd_opt_bl = 1, & ! blocking drag - gsd_diss_ht_opt = 0 - -! Parameters for bounding the scale-adaptive variability: -! Small-scale GWD + turbulent form drag - real(kind=kind_phys), parameter :: dxmin_ss = 1000., & - & dxmax_ss = 12000. ! min,max range of tapering (m) -! Large-scale GWD + blocking - real(kind=kind_phys), parameter :: dxmin_ls = 3000., & - & dxmax_ls = 13000. ! min,max range of tapering (m) - real(kind=kind_phys), dimension(im) :: ss_taper, ls_taper ! small- and large-scale tapering factors (-) -! -! Variables for limiting topographic standard deviation (var) - real(kind=kind_phys), parameter :: varmax_ss = 50., & - varmax_fd = 150., & - beta_ss = 0.1, & - beta_fd = 0.2 - real(kind=kind_phys) :: var_temp, var_temp2 - -! added Beljaars orographic form drag - real(kind=kind_phys), dimension(im,km) :: utendform,vtendform - real(kind=kind_phys) :: a1,a2,wsp - real(kind=kind_phys) :: H_efold - -! critical richardson number for wave breaking : ! larger drag with larger value - real(kind=kind_phys), parameter :: ric = 0.25 - real(kind=kind_phys), parameter :: dw2min = 1. - real(kind=kind_phys), parameter :: rimin = -100. - real(kind=kind_phys), parameter :: bnv2min = 1.0e-5 - real(kind=kind_phys), parameter :: efmin = 0.0 - real(kind=kind_phys), parameter :: efmax = 10.0 - real(kind=kind_phys), parameter :: xl = 4.0e4 - real(kind=kind_phys), parameter :: critac = 1.0e-5 - real(kind=kind_phys), parameter :: gmax = 1. - real(kind=kind_phys), parameter :: veleps = 1.0 - real(kind=kind_phys), parameter :: factop = 0.5 - real(kind=kind_phys), parameter :: frc = 1.0 - real(kind=kind_phys), parameter :: ce = 0.8 - real(kind=kind_phys), parameter :: cg = 0.5 - integer,parameter :: kpblmin = 2 - -! -! local variables -! - integer :: i,j,k,lcap,lcapp1,nwd,idir, & - klcap,kp1 -! - real(kind=kind_phys) :: rcs,csg,fdir,cleff,cleff_ss,cs, & - rcsks,wdir,ti,rdz,tem2,dw2,shr2, & - bvf2,rdelks,wtkbj,tem,gfobnv,hd,fro, & - rim,temc,tem1,efact,temv,dtaux,dtauy, & - dtauxb,dtauyb,eng0,eng1 -! - logical :: ldrag(im),icrilv(im), & - flag(im),kloop1(im) -! - real(kind=kind_phys) :: taub(im),taup(im,km+1), & - xn(im),yn(im), & - ubar(im),vbar(im), & - fr(im),ulow(im), & - rulow(im),bnv(im), & - oa(im),ol(im), & - oass(im),olss(im), & - roll(im),dtfac(im), & - brvf(im),xlinv(im), & - delks(im),delks1(im), & - bnv2(im,km),usqj(im,km), & - taud_ls(im,km),taud_bl(im,km), & - ro(im,km), & - vtk(im,km),vtj(im,km), & - zlowtop(im),velco(im,km-1), & - coefm(im),coefm_ss(im) -! - integer :: kbl(im),klowtop(im) - integer,parameter :: mdir=8 - !integer :: nwdir(mdir) - !data nwdir/6,7,5,8,2,3,1,4/ - integer, parameter :: nwdir(8) = (/6,7,5,8,2,3,1,4/) -! -! variables for flow-blocking drag -! - real(kind=kind_phys),parameter :: frmax = 10. - real(kind=kind_phys),parameter :: olmin = 1.0e-5 - real(kind=kind_phys),parameter :: odmin = 0.1 - real(kind=kind_phys),parameter :: odmax = 10. - real(kind=kind_phys),parameter :: erad = 6371.315e+3 - integer :: komax(im) - integer :: kblk - real(kind=kind_phys) :: cd - real(kind=kind_phys) :: zblk,tautem - real(kind=kind_phys) :: pe,ke - real(kind=kind_phys) :: delx,dely - real(kind=kind_phys) :: dxy4(im,4),dxy4p(im,4) - real(kind=kind_phys) :: dxy(im),dxyp(im) - real(kind=kind_phys) :: ol4p(4),olp(im),od(im) - real(kind=kind_phys) :: taufb(im,km+1) - - character(len=*), intent(out) :: errmsg - integer, intent(out) :: errflg - - integer :: udtend, vdtend, Tdtend - - ! Calculate inverse of gravitational acceleration - g_inv = 1./G - - ! Initialize CCPP error handling variables - errmsg = '' - errflg = 0 - - ! Initialize local variables - var_temp2 = 0. - udtend = -1 - vdtend = -1 - Tdtend = -1 - - if(ldiag3d) then - udtend = dtidx(index_of_x_wind,index_of_process_orographic_gwd) - vdtend = dtidx(index_of_y_wind,index_of_process_orographic_gwd) - Tdtend = dtidx(index_of_temperature,index_of_process_orographic_gwd) - endif - -!-------------------------------------------------------------------- -! SCALE-ADPTIVE PARAMETER FROM GFS GWD SCHEME -!-------------------------------------------------------------------- -! parameter (cdmb = 1.0) ! non-dim sub grid mtn drag Amp (*j*) -! non-dim sub grid mtn drag Amp (*j*) -! cdmb = 1.0/float(IMX/192) -! cdmb = 192.0/float(IMX) - cdmb = 4.0 * 192.0/float(IMX) - if (cdmbgwd(1) >= 0.0) cdmb = cdmb * cdmbgwd(1) - -!>-# Orographic Gravity Wave Drag Section - kpblmax = km / 2 ! maximum pbl height : # of vertical levels / 2 -! -! Scale cleff between IM=384*2 and 192*2 for T126/T170 and T62 -! - if (imx > 0) then -! cleff = 1.0E-5 * SQRT(FLOAT(IMX)/384.0) -! cleff = 1.0E-5 * SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF! -! cleff = 0.5E-5 * SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF! -! cleff = 1.0E-5 * SQRT(FLOAT(IMX)/192)/float(IMX/192) -! cleff = 1.0E-5 / SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF! - cleff = 0.5E-5 / SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF! -! hmhj for ndsl -! jw cleff = 0.1E-5 / SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF! -! cleff = 2.0E-5 * SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF! -! cleff = 2.5E-5 * SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF! - endif - if (cdmbgwd(2) >= 0.0) cleff = cleff * cdmbgwd(2) -!-------------------------------------------------------------------- -! END SCALE-ADPTIVE PARAMETER SECTION -!-------------------------------------------------------------------- -! -!---- constants -! - rcl = 1. - rcs = sqrt(rcl) - cs = 1. / sqrt(rcl) - csg = cs * g - lcap = km - lcapp1 = lcap + 1 - fdir = mdir / (2.0*pi) - - 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 this module - else - xland(i)=2.0 - endif - RDXZB(i) = 0.0 - enddo - -!--- calculate scale-aware tapering factors -do i=1,im - if ( dx(i) .ge. dxmax_ls ) then - ls_taper(i) = 1. - else - if ( dx(i) .le. dxmin_ls) then - ls_taper(i) = 0. - else - ls_taper(i) = 0.5 * ( SIN(pi*(dx(i)-0.5*(dxmax_ls+dxmin_ls))/ & - (dxmax_ls-dxmin_ls)) + 1. ) - endif - endif -enddo - -do i=1,im - if ( dx(i) .ge. dxmax_ss ) then - ss_taper(i) = 1. - else - if ( dx(i) .le. dxmin_ss) then - ss_taper(i) = 0. - else - ss_taper(i) = dxmax_ss * (1. - dxmin_ss/dx(i))/(dxmax_ss-dxmin_ss) - endif - 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 - delx = dx(i) - dely = dx(i) - dxy4(i,1) = delx - dxy4(i,2) = dely - dxy4(i,3) = sqrt(delx*delx + dely*dely) - dxy4(i,4) = dxy4(i,3) - dxy4p(i,1) = dxy4(i,2) - dxy4p(i,2) = dxy4(i,1) - dxy4p(i,3) = dxy4(i,4) - dxy4p(i,4) = dxy4(i,3) -enddo -! -!-----initialize arrays -! - dtaux = 0.0 - dtauy = 0.0 - do i = its,im - klowtop(i) = 0 - kbl(i) = 0 - enddo -! - do i = its,im - xn(i) = 0.0 - yn(i) = 0.0 - ubar (i) = 0.0 - vbar (i) = 0.0 - roll (i) = 0.0 - taub (i) = 0.0 - oa(i) = 0.0 - ol(i) = 0.0 - oass(i) = 0.0 - olss(i) = 0.0 - ulow (i) = 0.0 - dtfac(i) = 1.0 - rstoch(i) = 0.0 - ldrag(i) = .false. - icrilv(i) = .false. - flag(i) = .true. - enddo - - do k = kts,km - do i = its,im - usqj(i,k) = 0.0 - bnv2(i,k) = 0.0 - vtj(i,k) = 0.0 - vtk(i,k) = 0.0 - taup(i,k) = 0.0 - taud_ls(i,k) = 0.0 - taud_bl(i,k) = 0.0 - dtaux2d(i,k) = 0.0 - dtauy2d(i,k) = 0.0 - enddo - enddo -! - if ( (gwd_opt == 33).or.(gwd_opt == 22) ) then - do i = its,im - dusfc_ls(i) = 0.0 - dvsfc_ls(i) = 0.0 - dusfc_bl(i) = 0.0 - dvsfc_bl(i) = 0.0 - dusfc_ss(i) = 0.0 - dvsfc_ss(i) = 0.0 - dusfc_fd(i) = 0.0 - dvsfc_fd(i) = 0.0 - enddo - do k = kts,km - do i = its,im - dtaux2d_ls(i,k)= 0.0 - dtauy2d_ls(i,k)= 0.0 - dtaux2d_bl(i,k)= 0.0 - dtauy2d_bl(i,k)= 0.0 - dtaux2d_ss(i,k)= 0.0 - dtauy2d_ss(i,k)= 0.0 - dtaux2d_fd(i,k)= 0.0 - dtauy2d_fd(i,k)= 0.0 - enddo - enddo - endif - - do i = its,im - taup(i,km+1) = 0.0 - xlinv(i) = 1.0/xl - dusfc(i) = 0.0 - dvsfc(i) = 0.0 - enddo -! -! initialize array for flow-blocking drag -! - taufb(1:im,1:km+1) = 0.0 - komax(1:im) = 0 -! - do k = kts,km - do i = its,im - vtj(i,k) = t1(i,k) * (1.+fv*q1(i,k)) - vtk(i,k) = vtj(i,k) / prslk(i,k) - ro(i,k) = 1./rd * prsl(i,k) / vtj(i,k) ! density kg/m**3 - enddo - enddo -! -! calculate mid-layer height (zl), interface height (zq), and layer depth (dz2). -! - !zq=0. - do k = kts,km - do i = its,im - !zq(i,k+1) = PHII(i,k+1)*g_inv - !dz2(i,k) = (PHII(i,k+1)-PHII(i,k))*g_inv - zl(i,k) = PHIL(i,k)*g_inv - enddo - enddo -! -! determine reference level: maximum of 2*var and pbl heights -! - do i = its,im - zlowtop(i) = 2. * var_stoch(i) - enddo -! - do i = its,im - kloop1(i) = .true. - enddo -! - do k = kts+1,km - do i = its,im - if(kloop1(i).and.zl(i,k)-zl(i,1).ge.zlowtop(i)) then - klowtop(i) = k+1 - kloop1(i) = .false. - endif - enddo - enddo -! - do i = its,im - kbl(i) = max(kpbl(i), klowtop(i)) - kbl(i) = max(min(kbl(i),kpblmax),kpblmin) - enddo -! -! determine the level of maximum orographic height -! - ! komax(:) = kbl(:) - komax(:) = klowtop(:) - 1 ! modification by NOAA/GSD March 2018 -! - do i = its,im - delks(i) = 1.0 / (prsi(i,1) - prsi(i,kbl(i))) - delks1(i) = 1.0 / (prsl(i,1) - prsl(i,kbl(i))) - enddo -! -! compute low level averages within pbl -! - do k = kts,kpblmax - do i = its,im - if (k.lt.kbl(i)) then - rcsks = rcs * del(i,k) * delks(i) - rdelks = del(i,k) * delks(i) - ubar(i) = ubar(i) + rcsks * u1(i,k) ! pbl u mean - vbar(i) = vbar(i) + rcsks * v1(i,k) ! pbl v mean - roll(i) = roll(i) + rdelks * ro(i,k) ! ro mean - endif - enddo - enddo -! -! figure out low-level horizontal wind direction -! -! nwd 1 2 3 4 5 6 7 8 -! wd w s sw nw e n ne se -! - do i = its,im - wdir = atan2(ubar(i),vbar(i)) + pi - idir = mod(nint(fdir*wdir),mdir) + 1 - nwd = nwdir(idir) - oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(i,mod(nwd-1,4)+1) - ol(i) = ol4(i,mod(nwd-1,4)+1) - ! Repeat for small-scale gwd - oass(i) = (1-2*int( (nwd-1)/4 )) * oa4ss(i,mod(nwd-1,4)+1) - olss(i) = ol4ss(i,mod(nwd-1,4)+1) - -! -!----- compute orographic width along (ol) and perpendicular (olp) -!----- the direction of wind -! - ol4p(1) = ol4(i,2) - ol4p(2) = ol4(i,1) - ol4p(3) = ol4(i,4) - ol4p(4) = ol4(i,3) - olp(i) = ol4p(mod(nwd-1,4)+1) -! -!----- compute orographic direction (horizontal orographic aspect ratio) -! - od(i) = olp(i)/max(ol(i),olmin) - od(i) = min(od(i),odmax) - od(i) = max(od(i),odmin) -! -!----- compute length of grid in the along(dxy) and cross(dxyp) wind directions -! - dxy(i) = dxy4(i,MOD(nwd-1,4)+1) - dxyp(i) = dxy4p(i,MOD(nwd-1,4)+1) - enddo -! -! END INITIALIZATION; BEGIN GWD CALCULATIONS: -! -IF ( (do_gsl_drag_ls_bl).and. & - ((gwd_opt_ls .EQ. 1).or.(gwd_opt_bl .EQ. 1)) ) then - - do i=its,im - - if ( ls_taper(i).GT.1.E-02 ) then - -! -!--- saving richardson number in usqj for migwdi -! - do k = kts,km-1 - ti = 2.0 / (t1(i,k)+t1(i,k+1)) - rdz = 1./(zl(i,k+1) - zl(i,k)) - tem1 = u1(i,k) - u1(i,k+1) - tem2 = v1(i,k) - v1(i,k+1) - dw2 = rcl*(tem1*tem1 + tem2*tem2) - shr2 = max(dw2,dw2min) * rdz * rdz - bvf2 = g*(g/cp+rdz*(vtj(i,k+1)-vtj(i,k))) * ti - usqj(i,k) = max(bvf2/shr2,rimin) - bnv2(i,k) = 2.0*g*rdz*(vtk(i,k+1)-vtk(i,k))/(vtk(i,k+1)+vtk(i,k)) - bnv2(i,k) = max( bnv2(i,k), bnv2min ) - enddo -! -!----compute the "low level" or 1/3 wind magnitude (m/s) -! - ulow(i) = max(sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i)), 1.0) - rulow(i) = 1./ulow(i) -! - do k = kts,km-1 - velco(i,k) = (0.5*rcs) * ((u1(i,k)+u1(i,k+1)) * ubar(i) & - + (v1(i,k)+v1(i,k+1)) * vbar(i)) - velco(i,k) = velco(i,k) * rulow(i) - if ((velco(i,k).lt.veleps) .and. (velco(i,k).gt.0.)) then - velco(i,k) = veleps - endif - enddo -! -! no drag when critical level in the base layer -! - ldrag(i) = velco(i,1).le.0. -! -! no drag when velco.lt.0 -! - do k = kpblmin,kpblmax - if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. velco(i,k).le.0. - enddo -! -! no drag when bnv2.lt.0 -! - do k = kts,kpblmax - if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. bnv2(i,k).lt.0. - enddo -! -!-----the low level weighted average ri is stored in usqj(1,1; im) -!-----the low level weighted average n**2 is stored in bnv2(1,1; im) -!---- this is called bnvl2 in phys_gwd_alpert_sub not bnv2 -!---- rdelks (del(k)/delks) vert ave factor so we can * instead of / -! - wtkbj = (prsl(i,1)-prsl(i,2)) * delks1(i) - bnv2(i,1) = wtkbj * bnv2(i,1) - usqj(i,1) = wtkbj * usqj(i,1) -! - do k = kpblmin,kpblmax - if (k .lt. kbl(i)) then - rdelks = (prsl(i,k)-prsl(i,k+1)) * delks1(i) - bnv2(i,1) = bnv2(i,1) + bnv2(i,k) * rdelks - usqj(i,1) = usqj(i,1) + usqj(i,k) * rdelks - endif - enddo -! - 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_stoch(i) .le. 0.0 -! -! set all ri low level values to the low level value -! - do k = kpblmin,kpblmax - if (k .lt. kbl(i)) usqj(i,k) = usqj(i,1) - enddo -! - if (.not.ldrag(i)) then - bnv(i) = sqrt( bnv2(i,1) ) - 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) - endif -! -! compute the base level stress and store it in taub -! calculate enhancement factor, number of mountains & aspect -! ratio const. use simplified relationship between standard -! deviation & critical hgt - - if (.not. ldrag(i)) then - efact = (oa(i) + 2.) ** (ce*fr(i)/frc) - efact = min( max(efact,efmin), efmax ) -!!!!!!! cleff (effective grid length) is highly tunable parameter -!!!!!!! the bigger (smaller) value produce weaker (stronger) wave drag -!WRF cleff = sqrt(dxy(i)**2. + dxyp(i)**2.) -!WRF cleff = 3. * max(dx(i),cleff) - coefm(i) = (1. + ol(i)) ** (oa(i)+1.) -!WRF xlinv(i) = coefm(i) / cleff - xlinv(i) = coefm(i) * cleff - tem = fr(i) * fr(i) * oc1(i) - gfobnv = gmax * tem / ((tem + cg)*bnv(i)) - if ( gwd_opt_ls .NE. 0 ) then - taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i) & - * ulow(i) * gfobnv * efact - else ! We've gotten what we need for the blocking scheme - taub(i) = 0.0 - end if - else - taub(i) = 0.0 - xn(i) = 0.0 - yn(i) = 0.0 - endif - - endif ! (ls_taper(i).GT.1.E-02) - - enddo ! do i=its,im - -ENDIF ! (do_gsl_drag_ls_bl).and.((gwd_opt_ls .EQ. 1).or.(gwd_opt_bl .EQ. 1)) - -!========================================================= -! add small-scale wavedrag for stable boundary layer -!========================================================= - XNBV=0. - tauwavex0=0. - tauwavey0=0. - density=1.2 - utendwave=0. - vtendwave=0. -! -IF ( do_gsl_drag_ss ) THEN - - do i=its,im - - if ( ss_taper(i).GT.1.E-02 ) then - ! - ! calculating potential temperature - ! - do k = kts,km - thx(i,k) = t1(i,k)/prslk(i,k) - enddo - ! - do k = kts,km - tvcon = (1.+fv*q1(i,k)) - thvx(i,k) = thx(i,k)*tvcon - enddo - - hpbl2 = hpbl(i)+10. - kpbl2 = kpbl(i) - !kvar = MIN(kpbl, k-level of var) - kvar = 1 - do k=kts+1,MAX(kpbl(i),kts+1) -! IF (zl(i,k)>2.*var(i) .or. zl(i,k)>2*varmax) then - IF (zl(i,k)>300.) then - kpbl2 = k - IF (k == kpbl(i)) then - hpbl2 = hpbl(i)+10. - ELSE - hpbl2 = zl(i,k)+10. - ENDIF - exit - ENDIF - enddo - 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) -! cleff_ss = 10. * max(dxmax_ss,cleff_ss) - cleff_ss = 0.1 * max(dxmax_ss,cleff_ss) ! WRF -! cleff_ss = 0.1 * 12000. - coefm_ss(i) = (1. + olss(i)) ** (oass(i)+1.) - xlinv(i) = coefm_ss(i) / cleff_ss - !govrth(i)=g/(0.5*(thvx(i,kpbl(i))+thvx(i,kts))) - govrth(i)=g/(0.5*(thvx(i,kpbl2)+thvx(i,kts))) - !XNBV=sqrt(govrth(i)*(thvx(i,kpbl(i))-thvx(i,kts))/hpbl(i)) - XNBV=sqrt(govrth(i)*(thvx(i,kpbl2)-thvx(i,kts))/hpbl2) -! - !if(abs(XNBV/u1(i,kpbl(i))).gt.xlinv(i))then - if(abs(XNBV/u1(i,kpbl2)).gt.xlinv(i))then - !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_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) - tauwavex0=tauwavex0*ss_taper(i) - else - tauwavex0=0. - endif -! - !if(abs(XNBV/v1(i,kpbl(i))).gt.xlinv(i))then - if(abs(XNBV/v1(i,kpbl2)).gt.xlinv(i))then - !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_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) - tauwavey0=tauwavey0*ss_taper(i) - else - tauwavey0=0. - endif - - do k=kts,kpbl(i) !MIN(kpbl2+1,km-1) -!original - !utendwave(i,k)=-1.*tauwavex0*2.*max((1.-zl(i,k)/hpbl(i)),0.)/hpbl(i) - !vtendwave(i,k)=-1.*tauwavey0*2.*max((1.-zl(i,k)/hpbl(i)),0.)/hpbl(i) -!new - utendwave(i,k)=-1.*tauwavex0*2.*max((1.-zl(i,k)/hpbl2),0.)/hpbl2 - vtendwave(i,k)=-1.*tauwavey0*2.*max((1.-zl(i,k)/hpbl2),0.)/hpbl2 -!mod-to be used in HRRRv3/RAPv4 - !utendwave(i,k)=-1.*tauwavex0 * max((1.-zl(i,k)/hpbl2),0.)**2 - !vtendwave(i,k)=-1.*tauwavey0 * max((1.-zl(i,k)/hpbl2),0.)**2 - enddo - endif - endif - - do k = kts,km - dudt(i,k) = dudt(i,k) + utendwave(i,k) - dvdt(i,k) = dvdt(i,k) + vtendwave(i,k) - dusfc(i) = dusfc(i) + utendwave(i,k) * del(i,k) - dvsfc(i) = dvsfc(i) + vtendwave(i,k) * del(i,k) - enddo - if(udtend>0) then - dtend(i,kts:km,udtend) = dtend(i,kts:km,udtend) + utendwave(i,kts:km)*deltim - endif - if(vdtend>0) then - dtend(i,kts:km,vdtend) = dtend(i,kts:km,vdtend) + vtendwave(i,kts:km)*deltim - endif - if ( (gwd_opt == 33).or.(gwd_opt == 22) ) then - do k = kts,km - dusfc_ss(i) = dusfc_ss(i) + utendwave(i,k) * del(i,k) - dvsfc_ss(i) = dvsfc_ss(i) + vtendwave(i,k) * del(i,k) - dtaux2d_ss(i,k) = utendwave(i,k) - dtauy2d_ss(i,k) = vtendwave(i,k) - enddo - endif - - endif ! if (ss_taper(i).GT.1.E-02) - - enddo ! i=its,im - -ENDIF ! if (do_gsl_drag_ss) - -!================================================================ -! Topographic Form Drag from Beljaars et al. (2004, QJRMS, equ. 16): -!================================================================ -IF ( do_gsl_drag_tofd ) THEN - - do i=its,im - - if ( ss_taper(i).GT.1.E-02 ) then - - utendform=0. - vtendform=0. - - 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_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_stoch(i),hpbl(i)) - H_efold = min(H_efold,1500.) - DO k=kts,km - wsp=SQRT(u1(i,k)**2 + v1(i,k)**2) - ! alpha*beta*Cmd*Ccorr*2.109 = 12.*1.*0.005*0.6*2.109 = 0.0759 - var_temp = 0.0759*EXP(-(zl(i,k)/H_efold)**1.5)*a2* & - zl(i,k)**(-1.2)*ss_taper(i) ! this is greater than zero - ! Note: This is a semi-implicit treatment of the time differencing - ! per Beljaars et al. (2004, QJRMS) - utendform(i,k) = - var_temp*wsp*u1(i,k)/(1. + var_temp*deltim*wsp) - vtendform(i,k) = - var_temp*wsp*v1(i,k)/(1. + var_temp*deltim*wsp) - !IF(zl(i,k) > 4000.) exit - ENDDO - ENDIF - - do k = kts,km - dudt(i,k) = dudt(i,k) + utendform(i,k) - dvdt(i,k) = dvdt(i,k) + vtendform(i,k) - dusfc(i) = dusfc(i) + utendform(i,k) * del(i,k) - dvsfc(i) = dvsfc(i) + vtendform(i,k) * del(i,k) - enddo - if(udtend>0) then - dtend(i,kts:km,udtend) = dtend(i,kts:km,udtend) + utendform(i,kts:km)*deltim - endif - if(vdtend>0) then - dtend(i,kts:km,vdtend) = dtend(i,kts:km,vdtend) + vtendform(i,kts:km)*deltim - endif - if ( (gwd_opt == 33).or.(gwd_opt == 22) ) then - do k = kts,km - dtaux2d_fd(i,k) = utendform(i,k) - dtauy2d_fd(i,k) = vtendform(i,k) - dusfc_fd(i) = dusfc_fd(i) + utendform(i,k) * del(i,k) - dvsfc_fd(i) = dvsfc_fd(i) + vtendform(i,k) * del(i,k) - enddo - endif - - endif ! if (ss_taper(i).GT.1.E-02) - - enddo ! i=its,im - -ENDIF ! if (do_gsl_drag_tofd) -!======================================================= -! More for the large-scale gwd component -IF ( (do_gsl_drag_ls_bl).and.(gwd_opt_ls .EQ. 1) ) THEN - - do i=its,im - - if ( ls_taper(i).GT.1.E-02 ) then - -! -! now compute vertical structure of the stress. - do k = kts,kpblmax - if (k .le. kbl(i)) taup(i,k) = taub(i) - enddo -! - do k = kpblmin, km-1 ! vertical level k loop! - kp1 = k + 1 -! -! unstablelayer if ri < ric -! unstable layer if upper air vel comp along surf vel <=0 (crit lay) -! at (u-c)=0. crit layer exists and bit vector should be set (.le.) -! - if (k .ge. kbl(i)) then - icrilv(i) = icrilv(i) .or. ( usqj(i,k) .lt. ric) & - .or. (velco(i,k) .le. 0.0) - brvf(i) = max(bnv2(i,k),bnv2min) ! brunt-vaisala frequency squared - brvf(i) = sqrt(brvf(i)) ! brunt-vaisala frequency - endif -! - if (k .ge. kbl(i) .and. (.not. ldrag(i))) then - if (.not.icrilv(i) .and. taup(i,k) .gt. 0.0 ) then - temv = 1.0 / velco(i,k) - tem1 = coefm(i)/dxy(i)*(ro(i,kp1)+ro(i,k))*brvf(i)* & - velco(i,k)*0.5 - hd = sqrt(taup(i,k) / tem1) - fro = brvf(i) * hd * temv -! -! rim is the minimum-richardson number by shutts (1985) - tem2 = sqrt(usqj(i,k)) - tem = 1. + tem2 * fro - rim = usqj(i,k) * (1.-fro) / (tem * tem) -! -! check stability to employ the 'saturation hypothesis' -! of lindzen (1981) except at tropospheric downstream regions -! - if (rim .le. ric) then ! saturation hypothesis! - if ((oa(i) .le. 0.).or.(kp1 .ge. kpblmin )) then - temc = 2.0 + 1.0 / tem2 - hd = velco(i,k) * (2.*sqrt(temc)-temc) / brvf(i) - taup(i,kp1) = tem1 * hd * hd - endif - else ! no wavebreaking! - taup(i,kp1) = taup(i,k) - endif - endif - endif - enddo -! - if(lcap.lt.km) then - do klcap = lcapp1,km - taup(i,klcap) = prsi(i,klcap) / prsi(i,lcap) * taup(i,lcap) - enddo - endif - - endif ! if ( ls_taper(i).GT.1.E-02 ) - - enddo ! do i=its,im - -ENDIF ! (do_gsl_drag_ls_bl).and.(gwd_opt_ls .EQ. 1) -!=============================================================== -!COMPUTE BLOCKING COMPONENT -!=============================================================== -IF ( (do_gsl_drag_ls_bl) .and. (gwd_opt_bl .EQ. 1) ) THEN - - do i=its,im - - if ( ls_taper(i).GT.1.E-02 ) then - - if (.not.ldrag(i)) then -! -!------- determine the height of flow-blocking layer -! - kblk = 0 - pe = 0.0 - do k = km, kpblmin, -1 - if(kblk.eq.0 .and. k.le.komax(i)) then - pe = pe + bnv2(i,k)*(zl(i,komax(i))-zl(i,k))* & - del(i,k)/g/ro(i,k) - ke = 0.5*((rcs*u1(i,k))**2.+(rcs*v1(i,k))**2.) -! -!---------- apply flow-blocking drag when pe >= ke -! - if(pe.ge.ke) then - kblk = k - kblk = min(kblk,kbl(i)) - zblk = zl(i,kblk)-zl(i,kts) - RDXZB(i) = real(k,kind=kind_phys) - endif - endif - enddo - if(kblk.ne.0) then -! -!--------- compute flow-blocking stress -! - cd = max(2.0-1.0/od(i),0.0) - taufb(i,kts) = 0.5 * roll(i) * coefm(i) / & - max(dxmax_ls,dxy(i))**2 * cd * dxyp(i) * & - olp(i) * zblk * ulow(i)**2 - tautem = taufb(i,kts)/float(kblk-kts) - do k = kts+1, kblk - taufb(i,k) = taufb(i,k-1) - tautem - enddo -! -!----------sum orographic GW stress and flow-blocking stress -! - ! taup(i,:) = taup(i,:) + taufb(i,:) ! Keep taup and taufb separate for now - endif - - endif ! if (.not.ldrag(i)) - - endif ! if ( ls_taper(i).GT.1.E-02 ) - - enddo ! do i=its,im - -ENDIF ! IF ( (do_gsl_drag_ls_bl) .and. (gwd_opt_bl .EQ. 1) ) -!=========================================================== -IF ( (do_gsl_drag_ls_bl) .and. & - (gwd_opt_ls .EQ. 1 .OR. gwd_opt_bl .EQ. 1) ) THEN - - do i=its,im - - if ( ls_taper(i) .GT. 1.E-02 ) then - -! -! calculate - (g)*d(tau)/d(pressure) and deceleration terms dtaux, dtauy -! - do k = kts,km - taud_ls(i,k) = 1. * (taup(i,k+1) - taup(i,k)) * csg / del(i,k) - taud_bl(i,k) = 1. * (taufb(i,k+1) - taufb(i,k)) * csg / del(i,k) - enddo -! -! limit de-acceleration (momentum deposition ) at top to 1/2 value -! the idea is some stuff must go out the 'top' - do klcap = lcap,km - taud_ls(i,klcap) = taud_ls(i,klcap) * factop - taud_bl(i,klcap) = taud_bl(i,klcap) * factop - enddo -! -! if the gravity wave drag would force a critical line -! in the lower ksmm1 layers during the next deltim timestep, -! then only apply drag until that critical line is reached. -! - do k = kts,kpblmax-1 - if (k .le. kbl(i)) then - if ((taud_ls(i,k)+taud_bl(i,k)).ne.0.) & - dtfac(i) = min(dtfac(i),abs(velco(i,k) & - /(deltim*rcs*(taud_ls(i,k)+taud_bl(i,k))))) - endif - enddo -! - do k = kts,km - taud_ls(i,k) = taud_ls(i,k)*dtfac(i)* ls_taper(i) *(1.-rstoch(i)) - taud_bl(i,k) = taud_bl(i,k)*dtfac(i)* ls_taper(i) *(1.-rstoch(i)) - - dtaux = taud_ls(i,k) * xn(i) - dtauy = taud_ls(i,k) * yn(i) - dtauxb = taud_bl(i,k) * xn(i) - dtauyb = taud_bl(i,k) * yn(i) - - !add blocking and large-scale contributions to tendencies - dudt(i,k) = dtaux + dtauxb + dudt(i,k) - dvdt(i,k) = dtauy + dtauyb + dvdt(i,k) - - if ( gsd_diss_ht_opt .EQ. 1 ) then - ! Calculate dissipation heating - ! Initial kinetic energy (at t0-dt) - eng0 = 0.5*( (rcs*u1(i,k))**2. + (rcs*v1(i,k))**2. ) - ! Kinetic energy after wave-breaking/flow-blocking - eng1 = 0.5*( (rcs*(u1(i,k)+(dtaux+dtauxb)*deltim))**2 + & - (rcs*(v1(i,k)+(dtauy+dtauyb)*deltim))**2 ) - ! Modify theta tendency - dtdt(i,k) = dtdt(i,k) + max((eng0-eng1),0.0)/cp/deltim - if ( Tdtend>0 ) then - dtend(i,k,Tdtend) = dtend(i,k,Tdtend) + max((eng0-eng1),0.0)/cp - endif - endif - - dusfc(i) = dusfc(i) + taud_ls(i,k)*xn(i)*del(i,k) + & - taud_bl(i,k)*xn(i)*del(i,k) - dvsfc(i) = dvsfc(i) + taud_ls(i,k)*yn(i)*del(i,k) + & - taud_bl(i,k)*yn(i)*del(i,k) - if(udtend>0) then - dtend(i,k,udtend) = dtend(i,k,udtend) + (taud_ls(i,k) * & - xn(i) + taud_bl(i,k) * xn(i)) * deltim - endif - if(vdtend>0) then - dtend(i,k,vdtend) = dtend(i,k,vdtend) + (taud_ls(i,k) * & - yn(i) + taud_bl(i,k) * yn(i)) * deltim - endif - - enddo - - ! Finalize dusfc and dvsfc diagnostics - dusfc(i) = (-1./g*rcs) * dusfc(i) - dvsfc(i) = (-1./g*rcs) * dvsfc(i) - - if ( (gwd_opt == 33).or.(gwd_opt == 22) ) then - do k = kts,km - dtaux2d_ls(i,k) = taud_ls(i,k) * xn(i) - dtauy2d_ls(i,k) = taud_ls(i,k) * yn(i) - dtaux2d_bl(i,k) = taud_bl(i,k) * xn(i) - dtauy2d_bl(i,k) = taud_bl(i,k) * yn(i) - dusfc_ls(i) = dusfc_ls(i) + dtaux2d_ls(i,k) * del(i,k) - dvsfc_ls(i) = dvsfc_ls(i) + dtauy2d_ls(i,k) * del(i,k) - dusfc_bl(i) = dusfc_bl(i) + dtaux2d_bl(i,k) * del(i,k) - dvsfc_bl(i) = dvsfc_bl(i) + dtauy2d_bl(i,k) * del(i,k) - enddo - endif - - endif ! if ( ls_taper(i) .GT. 1.E-02 ) - - enddo ! do i=its,im - -ENDIF ! (do_gsl_drag_ls_bl).and.(gwd_opt_ls.EQ.1 .OR. gwd_opt_bl.EQ.1) - -if ( (gwd_opt == 33).or.(gwd_opt == 22) ) then - ! Finalize dusfc and dvsfc diagnostics - do i = its,im - dusfc_ls(i) = (-1./g*rcs) * dusfc_ls(i) - dvsfc_ls(i) = (-1./g*rcs) * dvsfc_ls(i) - dusfc_bl(i) = (-1./g*rcs) * dusfc_bl(i) - dvsfc_bl(i) = (-1./g*rcs) * dvsfc_bl(i) - dusfc_ss(i) = (-1./g*rcs) * dusfc_ss(i) - dvsfc_ss(i) = (-1./g*rcs) * dvsfc_ss(i) - dusfc_fd(i) = (-1./g*rcs) * dusfc_fd(i) - dvsfc_fd(i) = (-1./g*rcs) * dvsfc_fd(i) - enddo -endif -! - return - end subroutine drag_suite_run -!------------------------------------------------------------------- -! - subroutine drag_suite_finalize() - end subroutine drag_suite_finalize - - end module drag_suite From def0e78c4c6e94ea32e9f526e1408ae31e2467db Mon Sep 17 00:00:00 2001 From: jeff beck Date: Thu, 17 Feb 2022 20:10:50 +0000 Subject: [PATCH 23/28] Update spp_mp to equal 7 in if statement. --- physics/mp_thompson.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index 01874db40..7c76ea933 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -495,7 +495,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & end if ! Set stochastic physics selection to apply all perturbations - if ( spp_mp ) then + if ( spp_mp==7 ) then spp_mp_opt=7 else spp_mp_opt=0 From 372febe9ab1c26a0116fba27b2d85f70ec165ab5 Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Mon, 21 Feb 2022 11:39:15 -0700 Subject: [PATCH 24/28] edit module_mp_thompson.F90 to remove optional keyword for rand_pert (AKA spp_wts_mp); change to assumed-shape and reduce rank to match rank of input (i,k) instead of (i,j,k) --- physics/module_mp_thompson.F90 | 14 ++++---------- 1 file changed, 4 insertions(+), 10 deletions(-) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 9a25f0315..c23b6d1d8 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -1026,7 +1026,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):: & re_cloud, re_ice, re_snow INTEGER, INTENT(IN) :: rand_perturb_on, kme_stoch - REAL, DIMENSION(ims:ime,kms:kme_stoch,jms:jme), INTENT(IN), OPTIONAL:: & + REAL, DIMENSION(:,:), INTENT(IN) :: & rand_pert INTEGER, INTENT(IN):: has_reqc, has_reqi, has_reqs @@ -1123,12 +1123,6 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & ! No need to test for every subcycling step test_only_once: if (first_time_step .and. istep==1) then ! Activate this code when removing the guard above - if (rand_perturb_on .ne. 0 .and. .not. present(rand_pert)) then - errmsg = 'Logic error in mp_gt_driver: random perturbations are on, ' // & - 'but optional argument rand_pert is not present' - errflg = 1 - return - end if if ( (present(tt) .and. (present(th) .or. present(pii))) .or. & (.not.present(tt) .and. .not.(present(th) .and. present(pii))) ) then @@ -1294,11 +1288,11 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & rand2 = 0.0 rand3 = 0.0 if (rand_perturb_on .ne. 0) then - if (MOD(rand_perturb_on,2) .ne. 0) rand1 = rand_pert(i,1,j) + if (MOD(rand_perturb_on,2) .ne. 0) rand1 = rand_pert(i,1) m = RSHIFT(ABS(rand_perturb_on),1) - if (MOD(m,2) .ne. 0) rand2 = rand_pert(i,1,j)*2. + if (MOD(m,2) .ne. 0) rand2 = rand_pert(i,1)*2. m = RSHIFT(ABS(rand_perturb_on),2) - if (MOD(m,2) .ne. 0) rand3 = 0.25*(rand_pert(i,1,j)+ABS(min_rand)) + if (MOD(m,2) .ne. 0) rand3 = 0.25*(rand_pert(i,1)+ABS(min_rand)) m = RSHIFT(ABS(rand_perturb_on),3) endif !+---+-----------------------------------------------------------------+ From eef1e235a14cf41ddd0a34c388b1eb9da795d605 Mon Sep 17 00:00:00 2001 From: jeff beck Date: Mon, 21 Feb 2022 19:10:23 +0000 Subject: [PATCH 25/28] Loop over i,j instead of i,1 for rand_pert field. --- physics/module_mp_thompson.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index c23b6d1d8..2478f828c 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -1288,11 +1288,11 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & rand2 = 0.0 rand3 = 0.0 if (rand_perturb_on .ne. 0) then - if (MOD(rand_perturb_on,2) .ne. 0) rand1 = rand_pert(i,1) + if (MOD(rand_perturb_on,2) .ne. 0) rand1 = rand_pert(i,j) m = RSHIFT(ABS(rand_perturb_on),1) - if (MOD(m,2) .ne. 0) rand2 = rand_pert(i,1)*2. + if (MOD(m,2) .ne. 0) rand2 = rand_pert(i,j)*2. m = RSHIFT(ABS(rand_perturb_on),2) - if (MOD(m,2) .ne. 0) rand3 = 0.25*(rand_pert(i,1)+ABS(min_rand)) + if (MOD(m,2) .ne. 0) rand3 = 0.25*(rand_pert(i,j)+ABS(min_rand)) m = RSHIFT(ABS(rand_perturb_on),3) endif !+---+-----------------------------------------------------------------+ From ac7cde77cc1a983b500570864e4b12db8364e325 Mon Sep 17 00:00:00 2001 From: jeff beck Date: Mon, 21 Feb 2022 19:23:11 +0000 Subject: [PATCH 26/28] Revert dimensions changes to rand_pert. --- physics/module_mp_thompson.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 2478f828c..c23b6d1d8 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -1288,11 +1288,11 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & rand2 = 0.0 rand3 = 0.0 if (rand_perturb_on .ne. 0) then - if (MOD(rand_perturb_on,2) .ne. 0) rand1 = rand_pert(i,j) + if (MOD(rand_perturb_on,2) .ne. 0) rand1 = rand_pert(i,1) m = RSHIFT(ABS(rand_perturb_on),1) - if (MOD(m,2) .ne. 0) rand2 = rand_pert(i,j)*2. + if (MOD(m,2) .ne. 0) rand2 = rand_pert(i,1)*2. m = RSHIFT(ABS(rand_perturb_on),2) - if (MOD(m,2) .ne. 0) rand3 = 0.25*(rand_pert(i,j)+ABS(min_rand)) + if (MOD(m,2) .ne. 0) rand3 = 0.25*(rand_pert(i,1)+ABS(min_rand)) m = RSHIFT(ABS(rand_perturb_on),3) endif !+---+-----------------------------------------------------------------+ From 14ca01effa097f0932f6132b66defa5b63eab5f6 Mon Sep 17 00:00:00 2001 From: jeff beck Date: Mon, 21 Feb 2022 19:23:11 +0000 Subject: [PATCH 27/28] Revert dimension changes to rand_pert. --- physics/module_mp_thompson.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 2478f828c..c23b6d1d8 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -1288,11 +1288,11 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & rand2 = 0.0 rand3 = 0.0 if (rand_perturb_on .ne. 0) then - if (MOD(rand_perturb_on,2) .ne. 0) rand1 = rand_pert(i,j) + if (MOD(rand_perturb_on,2) .ne. 0) rand1 = rand_pert(i,1) m = RSHIFT(ABS(rand_perturb_on),1) - if (MOD(m,2) .ne. 0) rand2 = rand_pert(i,j)*2. + if (MOD(m,2) .ne. 0) rand2 = rand_pert(i,1)*2. m = RSHIFT(ABS(rand_perturb_on),2) - if (MOD(m,2) .ne. 0) rand3 = 0.25*(rand_pert(i,j)+ABS(min_rand)) + if (MOD(m,2) .ne. 0) rand3 = 0.25*(rand_pert(i,1)+ABS(min_rand)) m = RSHIFT(ABS(rand_perturb_on),3) endif !+---+-----------------------------------------------------------------+ From 88d0dd3939d3d1f65adb0ad4e9b0b51803ae1271 Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Mon, 21 Feb 2022 15:33:34 -0700 Subject: [PATCH 28/28] change optional and explicitly-shaped SPP arrays to non-optional and assumed-shape in module_bl_mynn and module_sf_mynn --- physics/module_bl_mynn.F90 | 7 ++++--- physics/module_sf_mynn.F90 | 4 ++-- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/physics/module_bl_mynn.F90 b/physics/module_bl_mynn.F90 index a492e50e0..3b0150e9e 100644 --- a/physics/module_bl_mynn.F90 +++ b/physics/module_bl_mynn.F90 @@ -4669,9 +4669,10 @@ SUBROUTINE mynn_bl_driver( & LOGICAL :: INITIALIZE_QKE ! Stochastic fields - INTEGER, INTENT(IN) ::spp_pbl - REAL, DIMENSION( ims:ime, kms:kme), INTENT(IN),OPTIONAL ::pattern_spp_pbl - REAL, DIMENSION(KTS:KTE) :: rstoch_col + INTEGER, INTENT(IN) :: spp_pbl + !GJF: this array must be assumed-shape since it's conditionally-allocated + REAL, DIMENSION(:,:), INTENT(IN) :: pattern_spp_pbl + REAL, DIMENSION(KTS:KTE) :: rstoch_col ! Substepping TKE INTEGER :: nsub diff --git a/physics/module_sf_mynn.F90 b/physics/module_sf_mynn.F90 index d1b3ce340..5f227750a 100644 --- a/physics/module_sf_mynn.F90 +++ b/physics/module_sf_mynn.F90 @@ -295,7 +295,8 @@ SUBROUTINE SFCLAY_mynn( & U3D,V3D, & th3d,pi3d - REAL, DIMENSION( ims:ime, kms:kme), OPTIONAL, & + !GJF: This array must be assumed-shape since it is conditionally-allocated + REAL, DIMENSION( :,: ), & INTENT(IN) :: pattern_spp_sfc !=================================== ! 2D VARIABLES @@ -3765,4 +3766,3 @@ REAL function psih_unstable(zolf,psi_opt) !======================================================================== END MODULE module_sf_mynn -