diff --git a/compns_stochy.F90 b/compns_stochy.F90 index 15b68bd..8d464c1 100644 --- a/compns_stochy.F90 +++ b/compns_stochy.F90 @@ -63,8 +63,8 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) shum_sigefold,spptint,shumint,skebint,skeb_npass,use_zmtnblck,new_lscale, & epbl,epbl_lscale,epbl_tau,iseed_epbl, & ocnsppt,ocnsppt_lscale,ocnsppt_tau,iseed_ocnsppt - namelist /nam_sfcperts/lndp_type,lndp_var_list, lndp_prt_list, iseed_lndp, & - lndp_tau,lndp_lscale + namelist /nam_sfcperts/lndp_type,lndp_model_type, lndp_var_list, lndp_prt_list, & + iseed_lndp, lndp_tau,lndp_lscale ! For SPP physics parameterization perterbations namelist /nam_sppperts/spp_var_list, spp_prt_list, iseed_spp, & spp_tau,spp_lscale,spp_sigtop1, spp_sigtop2,spp_stddev_cutoff @@ -101,9 +101,18 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) ! perturbations are assigned once at the start of the forecast ! LNDP_TYPE = 2 ! this is the newer land pert scheme, introduced and tested for impact on UFS/GDAS cycling stsyem -! perturbations are assigned at each time step (for state variables), or each time parameters are updated -! and the perturbations evolve over time. +! see https://journals.ametsoc.org/view/journals/hydr/22/8/JHM-D-21-0016.1.xml lndp_type = 0 ! +! LNDP_MODEL_TYPE +! integer indicating the model type for applying perturbations for lndp_type=2 scheme. +! 1 - global model +! (cycling of prognostic variables between DA cycles, +! parameters may be periodically updated during forecast) +! 2 - regional model +! (short foreast, prognotic variables re-initialized each forecast, +! parameters not updated during forecast) +! 3 - special case to apply perturbations only at start of forecast. + lndp_model_type = 0 ! lndp_lscale = -999. ! length scales lndp_tau = -999. ! time scales iseed_lndp = 0 ! random seeds (if 0 use system clock) @@ -203,26 +212,32 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) endif ENDIF ! compute frequencty to estimate dissipation timescale - IF (skebint == 0.) skebint=deltim - nsskeb=nint(skebint/deltim) ! skebint in seconds - IF(nsskeb<=0 .or. abs(nsskeb-skebint/deltim)>tol) THEN - WRITE(0,*) "SKEB interval is invalid",skebint - iret=9 - return - ENDIF - IF (spptint == 0.) spptint=deltim - nssppt=nint(spptint/deltim) ! spptint in seconds - IF(nssppt<=0 .or. abs(nssppt-spptint/deltim)>tol) THEN - WRITE(0,*) "SPPT interval is invalid",spptint - iret=9 - return - ENDIF - IF (shumint == 0.) shumint=deltim - nsshum=nint(shumint/deltim) ! shumint in seconds - IF(nsshum<=0 .or. abs(nsshum-shumint/deltim)>tol) THEN - WRITE(0,*) "SHUM interval is invalid",shumint - iret=9 - return + IF (do_skeb) THEN + IF (skebint == 0.) skebint=deltim + nsskeb=nint(skebint/deltim) ! skebint in seconds + IF(nsskeb<=0 .or. abs(nsskeb-skebint/deltim)>tol) THEN + WRITE(0,*) "SKEB interval is invalid",skebint + iret=9 + return + ENDIF + ENDIF + IF (do_sppt) THEN + IF (spptint == 0.) spptint=deltim + nssppt=nint(spptint/deltim) ! spptint in seconds + IF(nssppt<=0 .or. abs(nssppt-spptint/deltim)>tol) THEN + WRITE(0,*) "SPPT interval is invalid",spptint + iret=9 + return + ENDIF + ENDIF + IF (do_shum) THEN + IF (shumint == 0.) shumint=deltim + nsshum=nint(shumint/deltim) ! shumint in seconds + IF(nsshum<=0 .or. abs(nsshum-shumint/deltim)>tol) THEN + WRITE(0,*) "SHUM interval is invalid",shumint + iret=9 + return + ENDIF ENDIF !calculate ntrunc if not supplied if (ntrunc .LT. 1) then @@ -270,13 +285,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) cycle else n_var_lndp=n_var_lndp+1 - lndp_var_list( n_var_lndp) = lndp_var_list(k) ! for lndp_type==2: - ! for state variables, unit is pert per hour - ! for parmaters, no time dimension in unit - ! since perturbations do not accumulate - ! (i.e., global_cycle overwrites the paramaters - ! each time it's called, so any previous perturbations - ! are lost). + lndp_var_list( n_var_lndp) = lndp_var_list(k) lndp_prt_list( n_var_lndp) = lndp_prt_list(k) endif enddo @@ -290,7 +299,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) if (lndp_type==1) then if (me==0) print*, & - 'lndp_type=1, land perturbations will be applied to selected paramaters, using older scheme designed for S2S fcst spread' + 'lndp_type=1, land perturbations will be applied to selected paramaters, using older scheme designed for S2S fcst spread with the noah LSM' ! sanity-check requested input do k =1,n_var_lndp select case (lndp_var_list(k)) @@ -305,6 +314,8 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) elseif(lndp_type==2) then if (me==0) print*, & 'land perturbations will be applied to selected paramaters, using newer scheme designed for DA ens spread' + ! check requested parameters have been coded. + ! note, Noah-MP specific checks will be done later (since need to know lsm type) do k =1,n_var_lndp select case (lndp_var_list(k)) case('vgf','smc','stc','alb', 'sal','emi','zol') @@ -315,6 +326,11 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) return end select enddo + if ( (lndp_model_type < 1) .or. (lndp_model_type > 3) ) then + print*, 'ERROR: for lndp_type=2, must have lndp_model_type = 1,2,3' + iret = 10 + return + endif endif case default @@ -367,6 +383,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) print *, ' do_shum : ', do_shum print *, ' do_skeb : ', do_skeb print *, ' lndp_type : ', lndp_type + print *, ' lndp_model_type : ', lndp_model_type if (lndp_type .NE. 0) print *, ' n_var_lndp : ', n_var_lndp print *, ' do_spp : ', do_spp print *, ' n_var_spp : ', n_var_spp @@ -429,7 +446,7 @@ subroutine compns_stochy_ocn (deltim,iret) epbl,epbl_lscale,epbl_tau,iseed_epbl, & ocnsppt,ocnsppt_lscale,ocnsppt_tau,iseed_ocnsppt - namelist /nam_sfcperts/lndp_type,lndp_var_list, lndp_prt_list, iseed_lndp, & + namelist /nam_sfcperts/lndp_type,lndp_model_type,lndp_var_list, lndp_prt_list, iseed_lndp, & lndp_tau,lndp_lscale diff --git a/docs/source/references.rst b/docs/source/references.rst index ac2a261..f2d9f43 100644 --- a/docs/source/references.rst +++ b/docs/source/references.rst @@ -3,7 +3,8 @@ References Berner, J., G. Shutts, M. Leutbecher, and T. Palmer, 2009: A spectral stochastic kinetic energy backscatter scheme and its impact on flow- dependent predictability in the ECMWF ensemble prediction system. J. Atmos. Sci., 66, 603–626, `doi:10.1175/2008JAS2677.1 `_ -Gehne, M., T. Hamill, G. Bates, P. Pegion, W. Kolczynski 2019: Land-surface parameter and state perturbations in the Global Ensemble Forecast System. Mon. Wea. Rev. 147, 1319–1340 `doi:10.1175/MWR-D-18-0057.1 `_ +Draper, C, 2021: Accounting for Land Model Uncertainty in Numerical Weather Prediction Ensemble Systems: Toward Ensemble-Based Coupled Land–Atmosphere Data Assimilation, J Hydromet, 22(8), 2089-2104, ``_ + Palmer, T. N., R. Buizza, F. Doblas-Reyes, T. Jung, M. Leutbecher, G. J. Shutts,M. Steinheimer, and A.Weisheimer, 2009: Stochastic parametrization and model uncertainty. ECMWF Tech. Memo. 598, 42 pp `doi:10.21957/ps8gbwbdv `_ diff --git a/docs/source/users_guide.rst b/docs/source/users_guide.rst index 0d25818..6e03919 100644 --- a/docs/source/users_guide.rst +++ b/docs/source/users_guide.rst @@ -2,7 +2,7 @@ Users Guide ================================================== The stochastic physics currently only works with the UFS-atmosphere model -Currently, 3 stochastic schemes are used operationally at NCEP/EMC: Stochastic Kinetic Energy Backscatter (SKEB; Berner et al., 2009), Stochastically Perturbed Physics Tendencies (SPPT; Palmer et al., 2009), and Specific Humidity perturbations (SHUM), which is inspired by Tompkins and Berner, 2008. In addition there is the ability to perturb certain land model/surface parameters (Gehne et al, 2019), and a cellular automata scheme (Bengtsson et al. 2019) which interacts directly with the convective parameterization. +Currently, 3 stochastic schemes are used operationally at NCEP/EMC: Stochastic Kinetic Energy Backscatter (SKEB; Berner et al., 2009), Stochastically Perturbed Physics Tendencies (SPPT; Palmer et al., 2009), and Specific Humidity perturbations (SHUM), which is inspired by Tompkins and Berner, 2008. In addition there is the ability to perturb certain land model/surface parameters (Draper, 2021), and a cellular automata scheme (Bengtsson et al. 2019) which interacts directly with the convective parameterization. SKEB adds wind perturbations to model state. Perturbations are random in space/time, but amplitude is determined by a smoothed dissipation estimate provided by the dynamical core. Addresses errors in the dynamics - more active in the mid-latitudes @@ -11,7 +11,7 @@ SPPT multiplies the physics tendencies by a random number O [0,2] before updatin SHUM multiply the low-level specific humidity by a small random number each time-step. It attempts to address missing physics (cold pools, gust fronts), most active in convective regions -Land surface perturbations allow for land surface parameters such as Albedo, Soil Hydraulic Conductivity, LAI, and roughness lengths to vary in space. Addresses error in the land model and land-atmosphere interactions. +Land surface perturbations can be applied to the land model parameters and land model prognostic variables. The scheme is intended to address errors in the land model and land-atmosphere interactions. Due to the model’s numerics, any stochastic perturbation needs to be correlated in space and time in order to have the desired effect of upscale growth of the perturbations. This is achieved by creating a random pattern that has a specified decorrelation length-scale and is a first order auto-regressive process AR(1) in time with a specified decorrelation time-scale. (The CA random pattern generator also satisfies this condition) diff --git a/lndp_apply_perts.F90 b/lndp_apply_perts.F90 index bd3ef9f..91022f4 100644 --- a/lndp_apply_perts.F90 +++ b/lndp_apply_perts.F90 @@ -1,6 +1,7 @@ module lndp_apply_perts_mod use kinddef, only : kind_dbl_prec + use stochy_namelist_def implicit none @@ -13,23 +14,53 @@ module lndp_apply_perts_mod !==================================================================== ! lndp_apply_perts !==================================================================== -! Driver for applying perturbations to sprecified land states or parameters +! Driver for applying perturbations to specified land states or parameters, +! following the LNDP_TYPE=2 scheme. ! Draper, July 2020. - subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & - dtf, kdt, lndp_each_step, & - n_var_lndp, lndp_var_list, lndp_prt_list, & +! Can select perturbations by specifying lndp_var_list and lndp_pert_list + +! Notes on how the perturbations are added. +! 1. Model prognostic variables. +! If running a long forecast or cycling DA system (as in global UFS), +! perturbing the prognostic variables only at the initial conditions will +! have very limited impact, and they should instead be perturbed at every time step. +! In this case, the pertrubations should be specified as a rate (pert/hr) +! to avoid the ensemble spread being dependent on the model time step. +! +! For a short forecast (~days, as in regional HRRR), can see impact from +! perturbing only the initial conditions. In this case, the perturbation +! is specified as an absolute value (not a rate). +! +! 2. Model parameters: +! The timing of how to perturb the parameters depends on how / whether +! the parameters are updated over time. For the UFS global system, global_cycle +! is periodically called to update the parameters (controlled by FHCYC). +! Each time it's called global_cycle overwrites most of the +! prior parameters (overwriting any perturbations applied to those +! parameters). Hence, the perturbations are applied only immediately after global_cycle +! has been called, and the parameters are not applied as a rate (since they +! don't accumulate). +! For the regional models, FHCYC is 0, and the global_cycle is not called, so +! can perturb parameters every time step. Hence, need to specify the perturbations +! as a rate. +! +! The above cases are controlled by the lndp_model_type variable. +! +! If adding perturbations to new parameters, need to check how/whether +! the parameters are updated by the model. + + subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsm_noahmp, iopt_dveg, & + lsoil, dtf, kdt, n_var_lndp, lndp_var_list, lndp_prt_list, & sfc_wts, xlon, xlat, stype, smcmax, smcmin, param_update_flag, & - smc, slc, stc, vfrac, alvsf, alnsf, alvwf, alnwf, facsf, facwf, & - snoalb, semis, zorll, ierr) + smc, slc, stc, vfrac, alnsf, alnwf, snoalb, semis, zorll, ierr) implicit none ! intent(in) integer, intent(in) :: blksz(:) - integer, intent(in) :: n_var_lndp, lsoil, kdt - logical, intent(in) :: lndp_each_step - integer, intent(in) :: lsm, lsm_noah, lsm_ruc + integer, intent(in) :: n_var_lndp, lsoil, kdt, iopt_dveg + integer, intent(in) :: lsm, lsm_noah, lsm_ruc, lsm_noahmp character(len=3), intent(in) :: lndp_var_list(:) real(kind=kind_dbl_prec), intent(in) :: lndp_prt_list(:) real(kind=kind_dbl_prec), intent(in) :: dtf @@ -37,8 +68,8 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & real(kind=kind_dbl_prec), intent(in) :: xlon(:,:) real(kind=kind_dbl_prec), intent(in) :: xlat(:,:) logical, intent(in) :: param_update_flag - ! true = parameters have been updated, apply perts - real(kind=kind_dbl_prec), intent(in) :: stype(:,:) + ! true = parameters have just been updated by global_cycle + integer, intent(in) :: stype(:,:) real(kind=kind_dbl_prec), intent(in) :: smcmax(:) real(kind=kind_dbl_prec), intent(in) :: smcmin(:) @@ -48,12 +79,8 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & real(kind=kind_dbl_prec), intent(inout) :: stc(:,:,:) real(kind=kind_dbl_prec), intent(inout) :: vfrac(:,:) real(kind=kind_dbl_prec), intent(inout) :: snoalb(:,:) - real(kind=kind_dbl_prec), intent(inout) :: alvsf(:,:) real(kind=kind_dbl_prec), intent(inout) :: alnsf(:,:) - real(kind=kind_dbl_prec), intent(inout) :: alvwf(:,:) real(kind=kind_dbl_prec), intent(inout) :: alnwf(:,:) - real(kind=kind_dbl_prec), intent(inout) :: facsf(:,:) - real(kind=kind_dbl_prec), intent(inout) :: facwf(:,:) real(kind=kind_dbl_prec), intent(inout) :: semis(:,:) real(kind=kind_dbl_prec), intent(inout) :: zorll(:,:) @@ -62,10 +89,13 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & ! local integer :: nblks, print_i, print_nb, i, nb - integer :: this_im, v, soiltyp, k - logical :: print_flag + !integer :: this_im, v, soiltyp, k + integer :: this_im, v, k + logical :: print_flag, do_pert_state, do_pert_param - real(kind=kind_dbl_prec) :: p, min_bound, max_bound, tmp_sic, pert, factor + real(kind=kind_dbl_prec) :: p, min_bound, max_bound, pert + real(kind=kind_dbl_prec) :: tmp_smc + real(kind=kind_dbl_prec) :: conv_hr2tstep, tfactor_state, tfactor_param real(kind=kind_dbl_prec), dimension(lsoil) :: zslayer, smc_vertscale, stc_vertscale ! decrease in applied pert with depth @@ -80,28 +110,76 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & ierr = 0 - if (lsm/=lsm_noah .and. lsm/=lsm_ruc) then - write(6,*) 'ERROR: lndp_apply_pert assumes LSM is Noah or RUC,', & + if (lsm/=lsm_noah .and. lsm/=lsm_ruc .and. lsm/=lsm_noahmp) then + write(6,*) 'ERROR: lndp_apply_pert assumes LSM is Noah, Noah-MP, or RUC,', & ' may need to adapt variable names for a different LSM' ierr=10 return endif - !write (0,*) 'Input to lndp_apply_pert' - !write (0,*) 'lsm, lsoil, lsm_ruc, lsoil_lsm =', lsm, lsoil, lsm_ruc, lsoil_lsm - !write (0,*) 'zs_lsm =', zs_lsm - !write (0,*) 'n_var_lndp, lndp_var_list =', n_var_lndp, lndp_var_list - !write (0,*) 'smcmin =',smcmin - - ! lndp_prt_list input is per hour, factor converts to per timestep - ! Do conversion only when variables are perturbed at every time step - if(lndp_each_step) then - factor = dtf/3600. - else - factor = 1. + if (lsm==lsm_noahmp) then + do v = 1,n_var_lndp + select case (trim(lndp_var_list(v))) + case ('alb','sal','emi','zol') + print*, & + 'ERROR: lndp_prt_list option in lndp_apply_pert', trim(lndp_var_list(v)) , & + ' has not been checked for Noah-MP. Please check how the parameter is set/updated ', & + ' before applying. Note: in Noah-MP many variables that have traditionally been', & + ' externally specified parameters are now prognostic. Also, many parameters are', & + ' set at each Noah-MP model call from data tables' + ierr = 10 + return + case ('veg') + if (iopt_dveg .NE. 4 ) then + print*, & + 'ERROR: veg perturbations have not yet been coded for dveg options other than 4' + ierr = 10 + return + endif + end select + enddo endif - if (lsm == lsm_noah) then + ! for perturbations applied as a rate, lndp_prt_list input is per hour. Converts to per timestep + conv_hr2tstep = dtf/3600. ! conversion factor from per hour to per tstep. + + ! determine whether updating state variables and/or parameters + + do_pert_state=.false. + do_pert_param=.false. + + select case (lndp_model_type) + case(1) ! global, perturb states every time step (pert applied as a rate) + ! perturb parameters only when they've been update (pert is not a rate) + do_pert_state=.true. + tfactor_state=conv_hr2tstep + if (param_update_flag) then + do_pert_param=.true. + tfactor_param=1. + endif + case(2) ! regional, perturb states only at first time step (pert is not a rate) + ! perurb parameters at every time step (pert is a rate) + if ( kdt == 2 ) then + do_pert_state=.true. + tfactor_state=1. + endif + do_pert_param = .true. + tfactor_param = conv_hr2tstep + case(3) ! special case to apply perturbations at initial time step only (pert is not a rate) + if ( kdt == 2 ) then + do_pert_state=.true. + tfactor_state=1. + do_pert_param=.true. + tfactor_param=1. + endif + case default + print*, & + 'ERROR: unrecognised lndp_model_type option in lndp_apply_pert, exiting', trim(lndp_var_list(v)) + ierr = 10 + return + end select + + if (lsm == lsm_noah .or. lsm == lsm_noahmp) then do k = 1, lsoil zslayer(k) = zs_noah(k) smc_vertscale(k) = smc_vertscale_noah(k) @@ -122,11 +200,6 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & do nb =1,nblks do i = 1, blksz(nb) - !if ( nint(Sfcprop(nb)%slmsk(i)) .NE. 1) cycle ! skip if not land - - !if ( ((isot == 1) .and. (soiltyp == 16)) & - ! .or.( (isot == 0) .and. (soiltyp == 9 )) ) cycle ! skip if land-ice - if ( smc(nb,i,1) .EQ. 1.) cycle ! skip non-soil (land-ice and non-land) ! set printing if ( (i==print_i) .and. (nb==print_nb) ) then @@ -141,94 +214,100 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & ! State updates - performed every cycle !================================================================= case('smc') - p=5. - soiltyp = int( stype(nb,i)+0.5 ) ! also need for maxsmc - min_bound = smcmin(soiltyp) - max_bound = smcmax(soiltyp) - - if ((lsm /= lsm_ruc) .or. (lsm == lsm_ruc .and. kdt == 2)) then - ! with RUC LSM perturb smc only at time step = 2, as in HRRR - do k=1,lsoil - !store frozen soil moisture - tmp_sic= smc(nb,i,k) - slc(nb,i,k) - - ! perturb total soil moisture - ! factor of sldepth*1000 converts from mm to m3/m3 - ! lndp_prt_list(v) = 0.3 in input.nml - pert = sfc_wts(nb,i,v)*smc_vertscale(k)*lndp_prt_list(v)/(zslayer(k)*1000.) - pert = pert*dtf/3600. ! lndp_prt_list input is per hour, convert to per timestep - ! (necessary for state vars only) - call apply_pert('smc',pert,print_flag, smc(nb,i,k),ierr,p,min_bound, max_bound) - - ! assign all of applied pert to the liquid soil moisture - slc(nb,i,k) = smc(nb,i,k) - tmp_sic - enddo - endif + if (do_pert_state) then + p=5. + min_bound = smcmin(stype(nb,i)) + max_bound = smcmax(stype(nb,i)) + + ! with RUC LSM perturb smc only at time step = 2, as in HRRR + do k=1,lsoil + ! apply perts to a copy of smc, retain original smc + ! for later update to liquid soil moisture. + ! note: previously we were saving the ice water content + ! (smc-slc) and subtracting this from the perturbed smc to + ! get the perturbed slc. This was introducing small errors in the slc + ! when passing back to the calling program, I think due to precision issues, + ! as the ice content is typically zero. Clara Draper, March, 2022. + tmp_smc = smc(nb,i,k) + + ! perturb total soil moisture + ! factor of sldepth*1000 converts from mm to m3/m3 + pert = sfc_wts(nb,i,v)*smc_vertscale(k)*lndp_prt_list(v)/(zslayer(k)*1000.) + pert = pert*tfactor_state + + call apply_pert('smc',pert,print_flag, tmp_smc,ierr,p,min_bound, max_bound) + + ! assign all of applied pert to the liquid soil moisture + slc(nb,i,k) = slc(nb,i,k) + tmp_smc - smc(nb,i,k) + smc(nb,i,k) = tmp_smc + + enddo + endif case('stc') + if (do_pert_state) then + do k=1,lsoil + pert = sfc_wts(nb,i,v)*stc_vertscale(k)*lndp_prt_list(v) + pert = tfactor_state + call apply_pert('stc',pert,print_flag, stc(nb,i,k),ierr) + enddo + endif - do k=1,lsoil - pert = sfc_wts(nb,i,v)*stc_vertscale(k)*lndp_prt_list(v) - pert = pert*dtf/3600. ! lndp_prt_list input is per hour, convert to per timestep - ! (necessary for state vars only) - call apply_pert('stc',pert,print_flag, stc(nb,i,k),ierr) - enddo !================================================================= - ! Parameter updates - only if param_update_flag = TRUE + ! Parameter updates !================================================================= + + ! are all of the params below included in noah? + case('vgf') ! vegetation fraction - if (param_update_flag .or. lndp_each_step) then + if (do_pert_param) then p =5. min_bound=0. max_bound=1. pert = sfc_wts(nb,i,v)*lndp_prt_list(v) - pert = pert*factor + pert = pert*tfactor_param call apply_pert ('vfrac',pert,print_flag, vfrac(nb,i), ierr,p,min_bound, max_bound) endif case('alb') ! albedo - if (param_update_flag .or. lndp_each_step) then + if (do_pert_param) then p =5. min_bound=0.0 max_bound=0.4 pert = sfc_wts(nb,i,v)*lndp_prt_list(v) - pert = pert*factor - !call apply_pert ('alvsf',pert,print_flag, alvsf(nb,i), ierr,p,min_bound, max_bound) + pert = pert*tfactor_param call apply_pert ('alnsf',pert,print_flag, alnsf(nb,i), ierr,p,min_bound, max_bound) - !call apply_pert ('alvwf',pert,print_flag, alvwf(nb,i), ierr,p,min_bound, max_bound) call apply_pert ('alnwf',pert,print_flag, alnwf(nb,i), ierr,p,min_bound, max_bound) - !call apply_pert ('facsf',pert,print_flag, facsf(nb,i), ierr,p,min_bound, max_bound) - !call apply_pert ('facwf',pert,print_flag, facwf(nb,i), ierr,p,min_bound, max_bound) endif case('sal') ! snow albedo - if (param_update_flag .or. lndp_each_step) then + if (do_pert_param) then p =5. min_bound=0.3 max_bound=0.85 pert = sfc_wts(nb,i,v)*lndp_prt_list(v) - pert = pert*factor + pert = pert*tfactor_param call apply_pert ('snoalb',pert,print_flag, snoalb(nb,i), ierr,p,min_bound, max_bound) endif case('emi') ! emissivity - if (param_update_flag .or. lndp_each_step) then + if (do_pert_param) then p =5. min_bound=0.8 max_bound=1. pert = sfc_wts(nb,i,v)*lndp_prt_list(v) - pert = pert*factor + pert = pert*tfactor_param call apply_pert ('semis',pert,print_flag, semis(nb,i), ierr,p,min_bound, max_bound) endif case('zol') ! land roughness length - if (param_update_flag .or. lndp_each_step) then + if (do_pert_param) then p =5. min_bound=0. max_bound=300. pert = sfc_wts(nb,i,v)*lndp_prt_list(v) - pert = pert*factor + pert = pert*tfactor_param call apply_pert ('zol',pert,print_flag, zorll(nb,i), ierr,p,min_bound, max_bound) endif case default @@ -282,12 +361,11 @@ subroutine apply_pert(vname,pert,print_flag, state,ierr,p,vmin, vmax) z = -1. + 2*(state - vmin)/(vmax - vmin) ! flat-top function state = state + pert*(1-abs(z**p)) else - state = state + pert + state = state + pert endif if (present(vmax)) state = min( state , vmax ) if (present(vmin)) state = max( state , vmin ) - !state = max( min( state , vmax ), vmin ) if ( print_flag ) then write(*,*) 'LNDP - applying lndp to ',vname, ', final value', state diff --git a/main.doc b/main.doc index eeedbef..dfcdb1e 100644 --- a/main.doc +++ b/main.doc @@ -6,7 +6,7 @@ * *The stochastic physics currently only works with the UFS-atmosphere model * - *Currently, 3 stochastic schemes are used operationally at NCEP/EMC: Stochastic Kinetic Energy Backscatter (SKEB; Berner et al., 2009), Stochastically Perturbed Physics Tendencies (SPPT; Palmer et al., 2009), and Specific Humidity perturbations (SHUM), which is inspired by Tompkins and Berner, 2008. In addition there is the ability to perturb certain land model/surface parameters (Gehne et al, 2019), and a cellular automata scheme (Bengtsson et al. 20XX) which interacts directly with the convective parameterization. + *Currently, 3 stochastic schemes are used operationally at NCEP/EMC: Stochastic Kinetic Energy Backscatter (SKEB; Berner et al., 2009), Stochastically Perturbed Physics Tendencies (SPPT; Palmer et al., 2009), and Specific Humidity perturbations (SHUM), which is inspired by Tompkins and Berner, 2008. In addition there is the ability to perturb certain land model/surface parameters (Draper, 2021), and a cellular automata scheme (Bengtsson et al. 20XX) which interacts directly with the convective parameterization. * *SKEB adds wind perturbations to model state. Perturbations are random in space/time, but amplitude is determined by a smoothed dissipation estimate provided by the dynamical core. *Addresses errors in the dynamics - more active in the mid-latitudes diff --git a/stochy_namelist_def.F90 b/stochy_namelist_def.F90 index f5a9373..1afc5d0 100644 --- a/stochy_namelist_def.F90 +++ b/stochy_namelist_def.F90 @@ -10,7 +10,7 @@ module stochy_namelist_def implicit none public - integer, parameter :: max_n_var_lndp = 6 ! must match value used in GFS_typedefs + integer, parameter :: max_n_var_lndp = 20 ! must match value used in GFS_typedefs integer, parameter :: max_n_var_spp = 6 ! must match value used in GFS_typedefs integer nssppt,nsshum,nsepbl,nsocnsppt,nsskeb,lon_s,lat_s,ntrunc @@ -37,6 +37,7 @@ module stochy_namelist_def integer n_var_lndp integer(kind=kind_dbl_prec),dimension(5) ::iseed_lndp integer lndp_type + integer lndp_model_type character(len=3), dimension(max_n_var_lndp) :: lndp_var_list real(kind=kind_dbl_prec), dimension(max_n_var_lndp) :: lndp_prt_list