From 2d4d5ef59dc9f92bf59918b1e2cbcfc72474a93d Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Tue, 8 Mar 2022 16:30:40 +0000 Subject: [PATCH 1/5] Updated lndp scheme to work with Noah-MP. Cleaned up namelist options for lndp scheme. Updated max_n_var_lndp at Jeff Ator's request. --- docs/source/references.rst | 3 +- docs/source/users_guide.rst | 4 +- lndp_apply_perts.F90 | 207 +++++++++++++++++++++++++----------- main.doc | 2 +- stochy_namelist_def.F90 | 3 +- 5 files changed, 150 insertions(+), 69 deletions(-) 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..022caa3 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,12 +14,45 @@ 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, +! combined with setting tfactor_state below. +! +! If adding new parameters, need to check how/whether +! the parameters are updated. + + subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsm_noahmp, 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) @@ -28,8 +62,7 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & ! 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) :: 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,7 +70,7 @@ 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 + ! true = parameters have just been updated by global_cycle real(kind=kind_dbl_prec), intent(in) :: stype(:,:) real(kind=kind_dbl_prec), intent(in) :: smcmax(:) real(kind=kind_dbl_prec), intent(in) :: smcmin(:) @@ -63,9 +96,10 @@ 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 + 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, tmp_sic, pert + 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 +114,69 @@ 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' + ierr = 10 + return + 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(0) + ! no perts requested, do nothing + 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) @@ -141,59 +216,63 @@ 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. + soiltyp = int( stype(nb,i)+0.5 ) ! also need for maxsmc + min_bound = smcmin(soiltyp) + max_bound = smcmax(soiltyp) + + ! 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*tfactor_state + + 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 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 + pert = pert*tfactor_param !call apply_pert ('alvsf',pert,print_flag, alvsf(nb,i), ierr,p,min_bound, max_bound) 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) @@ -202,33 +281,33 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & !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 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 From 7f8f0b0a86fed435f6241dfd243706e83d88de9e Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Fri, 18 Mar 2022 21:06:15 +0000 Subject: [PATCH 2/5] added lndp_model_type option. --- compns_stochy.F90 | 83 ++++++++++++++++++++++++++------------------ lndp_apply_perts.F90 | 11 +++--- 2 files changed, 54 insertions(+), 40 deletions(-) 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/lndp_apply_perts.F90 b/lndp_apply_perts.F90 index 022caa3..770f3d1 100644 --- a/lndp_apply_perts.F90 +++ b/lndp_apply_perts.F90 @@ -167,8 +167,6 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsm_noahmp, lsoil,& do_pert_param=.true. tfactor_param=1. endif - case(0) - ! no perts requested, do nothing case default print*, & 'ERROR: unrecognised lndp_model_type option in lndp_apply_pert, exiting', trim(lndp_var_list(v)) @@ -197,11 +195,6 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsm_noahmp, 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 @@ -210,6 +203,10 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsm_noahmp, lsoil,& print_flag=.false. endif + if (print_flag) then + write(6,*) 'LNDPtmp', vfrac(nb,i) + endif + do v = 1,n_var_lndp select case (trim(lndp_var_list(v))) !================================================================= From b5d029750dd75ab134d386e601d7051296daf857 Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Fri, 18 Mar 2022 21:29:40 +0000 Subject: [PATCH 3/5] Added check for veg option for lndp scheme with Noah-Mp. --- lndp_apply_perts.F90 | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/lndp_apply_perts.F90 b/lndp_apply_perts.F90 index 770f3d1..01c7128 100644 --- a/lndp_apply_perts.F90 +++ b/lndp_apply_perts.F90 @@ -45,14 +45,13 @@ module lndp_apply_perts_mod ! 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, -! combined with setting tfactor_state below. +! The above cases are controlled by the lndp_model_type variable. ! -! If adding new parameters, need to check how/whether -! the parameters are updated. +! 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, lsoil,& - dtf, kdt, n_var_lndp, lndp_var_list, lndp_prt_list, & + 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) @@ -61,7 +60,7 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsm_noahmp, lsoil,& ! intent(in) integer, intent(in) :: blksz(:) - integer, intent(in) :: n_var_lndp, lsoil, kdt + 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(:) @@ -128,9 +127,17 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsm_noahmp, lsoil,& 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' + ' 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 ~=4 ) then + print*, & + 'ERROR: veg perturbations have not yet been coded for dveg options other than 4' + ierr = 10 + return end select enddo endif From 4412cd3e1c41e5c85ac6da2e27667ba54b6e6e2d Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Fri, 18 Mar 2022 22:12:14 +0000 Subject: [PATCH 4/5] minor bug fix. --- lndp_apply_perts.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/lndp_apply_perts.F90 b/lndp_apply_perts.F90 index 01c7128..2ce512e 100644 --- a/lndp_apply_perts.F90 +++ b/lndp_apply_perts.F90 @@ -133,11 +133,12 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsm_noahmp, iopt_dveg ierr = 10 return case ('veg') - if (iopt_dveg ~=4 ) then + 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 From c05fb060901d5a2fa65523d80e16c2881af2d9d2 Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Wed, 23 Mar 2022 16:16:39 +0000 Subject: [PATCH 5/5] removed unused albedo terms, updated soil type to be an integer, fixed precisiom issue in slc update. --- lndp_apply_perts.F90 | 48 +++++++++++++++++++------------------------- 1 file changed, 21 insertions(+), 27 deletions(-) diff --git a/lndp_apply_perts.F90 b/lndp_apply_perts.F90 index 2ce512e..91022f4 100644 --- a/lndp_apply_perts.F90 +++ b/lndp_apply_perts.F90 @@ -53,8 +53,7 @@ module lndp_apply_perts_mod 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 @@ -70,7 +69,7 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsm_noahmp, iopt_dveg real(kind=kind_dbl_prec), intent(in) :: xlat(:,:) logical, intent(in) :: param_update_flag ! true = parameters have just been updated by global_cycle - real(kind=kind_dbl_prec), intent(in) :: stype(:,:) + integer, intent(in) :: stype(:,:) real(kind=kind_dbl_prec), intent(in) :: smcmax(:) real(kind=kind_dbl_prec), intent(in) :: smcmin(:) @@ -80,12 +79,8 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsm_noahmp, iopt_dveg 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(:,:) @@ -94,10 +89,12 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsm_noahmp, iopt_dveg ! local integer :: nblks, print_i, print_nb, i, nb - integer :: this_im, v, soiltyp, k + !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 + 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 @@ -211,10 +208,6 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsm_noahmp, iopt_dveg print_flag=.false. endif - if (print_flag) then - write(6,*) 'LNDPtmp', vfrac(nb,i) - endif - do v = 1,n_var_lndp select case (trim(lndp_var_list(v))) !================================================================= @@ -223,25 +216,31 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsm_noahmp, iopt_dveg case('smc') if (do_pert_state) then p=5. - soiltyp = int( stype(nb,i)+0.5 ) ! also need for maxsmc - min_bound = smcmin(soiltyp) - max_bound = smcmax(soiltyp) + 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 - !store frozen soil moisture - tmp_sic= smc(nb,i,k) - slc(nb,i,k) + ! 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 - ! 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*tfactor_state - call apply_pert('smc',pert,print_flag, smc(nb,i,k),ierr,p,min_bound, max_bound) + 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) = smc(nb,i,k) - tmp_sic + slc(nb,i,k) = slc(nb,i,k) + tmp_smc - smc(nb,i,k) + smc(nb,i,k) = tmp_smc + enddo endif @@ -278,12 +277,8 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsm_noahmp, iopt_dveg pert = sfc_wts(nb,i,v)*lndp_prt_list(v) pert = pert*tfactor_param - !call apply_pert ('alvsf',pert,print_flag, alvsf(nb,i), ierr,p,min_bound, max_bound) 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 (do_pert_param) then @@ -366,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