Merge pull request #56 from ClaraDraper-NOAA/feature/lndp_noahmp2

Feature/lndp noahmp2 : add land perturbation scheme for Noah-MP fractional veg
NOAA-PSL · Apr 27, 2022 · 5232846 · 5232846
2 parents 2b1240f + 0e409df
commit 5232846
Show file tree

Hide file tree

Showing 6 changed files with 219 additions and 122 deletions.
diff --git 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
@@ -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 <https://journals.ametsoc.org/doi/full/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 <https://journals.ametsoc.org/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, `<https://journals.ametsoc.org/view/journals/hydr/22/8/JHM-D-21-0016.1.xml>`_
+
 
 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 <https://www.ecmwf.int/node/11577>`_
 

diff --git 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)