Updated land perturbation scheme

CMakeLists.txt

@@ -53,8 +53,10 @@ add_library(
     ./cellular_automata_global.F90
     ./cellular_automata_sgs.F90
     ./update_ca.F90
+    ./lndp_apply_perts.F90 
 )
 
 target_link_libraries(stochastic_physics sp::sp_d)
 target_link_libraries(stochastic_physics fms)
+target_link_libraries(stochastic_physics gfsphysics)
 

compns_stochy.F90

@@ -61,8 +61,8 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
       skeb,skeb_tau,skeb_vdof,skeb_lscale,iseed_skeb,skeb_vfilt,skeb_diss_smooth, &
       skeb_sigtop1,skeb_sigtop2,skebnorm,sppt_sigtop1,sppt_sigtop2,&
       shum_sigefold,spptint,shumint,skebint,skeb_npass,use_zmtnblck,new_lscale
-      namelist /nam_sfcperts/nsfcpert,pertz0,pertshc,pertzt,pertlai, & ! mg, sfcperts
-      pertvegf,pertalb,iseed_sfc,sfc_tau,sfc_lscale,sppt_land
+      namelist /nam_sfcperts/lndp_type,lndp_var_list, lndp_prt_list, iseed_lndp, & 
+      lndp_tau,lndp_lscale 
 
       rerth  =6.3712e+6      ! radius of earth (m)
       tol=0.01  ! tolerance for calculations
@@ -75,27 +75,30 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
       sppt             = -999.  ! stochastic physics tendency amplitude
       shum             = -999.  ! stochastic boundary layer spf hum amp
       skeb             = -999.  ! stochastic KE backscatter amplitude
-      ! mg, sfcperts
-      pertz0           = -999.  ! momentum roughness length amplitude
-      pertshc          = -999.  ! soil hydraulic conductivity amp
-      pertzt           = -999.  ! mom/heat roughness length amplitude
-      pertlai          = -999.  ! leaf area index amplitude
-      pertvegf         = -999.  ! vegetation fraction amplitude
-      pertalb          = -999.  ! albedo perturbations amplitude
+      lndp_var_list  = 'XXX'
+      lndp_prt_list  = -999.
 ! logicals
       do_sppt = .false.
       use_zmtnblck = .false.
       new_lscale = .false.
       do_shum = .false.
       do_skeb = .false.
-      ! mg, sfcperts
-      do_sfcperts = .false.
-      sppt_land = .false.
-      nsfcpert = 0
-! for sfcperts random patterns
-      sfc_lscale  = -999.       ! length scales
-      sfc_tau     = -999.       ! time scales
-      iseed_sfc   = 0           ! random seeds (if 0 use system clock)
+! C. Draper July 2020.
+! input land pert variables: 
+! LNDP_TYPE = 0
+! no explicit land perturbations
+! LNDP_Type =  1 
+! this is the initial land sfc pert scheme, introduced and tested for impact on GEFS forecasts.
+! see https://journals.ametsoc.org/doi/full/10.1175/MWR-D-18-0057.1
+! 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. 
+      lndp_type = 0 !
+      lndp_lscale  = -999.       ! length scales
+      lndp_tau     = -999.       ! time scales
+      iseed_lndp   = 0           ! random seeds (if 0 use system clock)
 ! for SKEB random patterns.
       skeb_vfilt       = 0
       skebint          = 0
@@ -209,11 +212,6 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
         iret=9
         return
       ENDIF
-! mg, sfcperts
-      IF (pertz0(1) > 0 .OR. pertshc(1) > 0 .OR. pertzt(1) > 0 .OR. &
-          pertlai(1) > 0 .OR. pertvegf(1) > 0 .OR. pertalb(1) > 0) THEN
-        do_sfcperts=.true.
-      ENDIF
 !calculate ntrunc if not supplied
      if (ntrunc .LT. 1) then  
         if (me==0) print*,'ntrunc not supplied, calculating'
@@ -224,6 +222,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
            if (shum(k).GT.0) l_min=min(shum_lscale(k),l_min)
            if (skeb(k).GT.0) l_min=min(skeb_lscale(k),l_min)
        enddo
+       if (lndp_type.GT.0) l_min=min(lndp_lscale(1),l_min)
        !ntrunc=1.5*circ/l_min
        ntrunc=circ/l_min
        if (me==0) print*,'ntrunc calculated from l_min',l_min,ntrunc
@@ -241,7 +240,76 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
         lon_s=lon_s*2
         if (me==0) print*,'gaussian grid not set, defining here',lon_s,lat_s
      endif
+
+! 
+! land perts  - parse nml input
 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+
+     select case (lndp_type)
+     case (0) 
+        if (me==0) print*, & 
+           'no land perturbations selected'
+     case (1,2) 
+        ! count requested pert variables
+        n_var_lndp= 0
+        do k =1,size(lndp_var_list)
+            if  ( (lndp_var_list(k) .EQ. 'XXX') .or. (lndp_prt_list(k) .LE. 0.) ) then
+               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_prt_list( n_var_lndp) = lndp_prt_list(k) 
+            endif
+        enddo 
+        if (n_var_lndp > max_n_var_lndp) then 
+              print*, 'ERROR: land perturbation requested for too many parameters', & 
+                       'increase max_n_var_lndp'
+              iret = 10 
+              return
+        endif
+
+
+        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' 
+           !  sanity-check requested input
+           do k =1,n_var_lndp
+               select case (lndp_var_list(k))
+               case('rz0','rzt','shc','lai','vgf','alb') 
+                   if (me==0) print*, 'land perturbation will be applied to ', lndp_var_list(k)
+               case default
+                  print*, 'ERROR: land perturbation requested for unknown parameter', lndp_var_list(k)
+                  iret = 10 
+                  return
+               end select 
+           enddo
+        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'
+               do k =1,n_var_lndp
+                   select case (lndp_var_list(k))
+                   case('vgf','smc','stc') 
+                       if (me==0) print*, 'land perturbation will be applied to ', lndp_var_list(k)
+                   case default
+                      print*, 'ERROR: land perturbation requested for new parameter - will need to be coded in lndp_apply_pert', lndp_var_list(k)
+                      iret = 10 
+                      return
+                   end select 
+               enddo
+        endif
+
+     case default 
+        if (me==0) print*, & 
+         'lndp_type out of range, set to 0 (none), 1 (for fcst spread), 2 (for cycling DA spread)'
+         iret = 10 
+         return 
+     end select 
 !
 !  All checks are successful.
 !
@@ -250,7 +318,8 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
          print *, ' do_sppt : ', do_sppt
          print *, ' do_shum : ', do_shum
          print *, ' do_skeb : ', do_skeb
-         print *, ' do_sfcperts : ', do_sfcperts
+         print *, ' lndp_type : ', lndp_type
+         if (lndp_type .NE. 0) print *, ' n_var_lndp : ', n_var_lndp
       endif
       iret = 0
 !

get_stochy_pattern.F90

@@ -6,11 +6,11 @@ module get_stochy_pattern_mod
  use spectral_layout_mod, only : ipt_lats_node_a, lat1s_a, lats_dim_a,      &
                                  lats_node_a, lon_dim_a, len_trie_ls,       &
                                  len_trio_ls, ls_dim, nodes, stochy_la2ga
- use stochy_namelist_def, only : nsfcpert, ntrunc, stochini
+ use stochy_namelist_def, only : n_var_lndp, ntrunc, stochini
  use stochy_data_mod, only : gg_lats, gg_lons, inttyp, nskeb, nshum, nsppt, &
                              rad2deg, rnlat, rpattern_sfc, rpattern_skeb,   &
                              rpattern_shum, rpattern_sppt, skebu_save,      &
-                             skebv_save, skeb_vwts, skeb_vpts, wlon
+                             skebv_save, skeb_vwts, skeb_vpts, wlon, nlndp
  use stochy_gg_def, only : coslat_a
  use stochy_patterngenerator_mod, only: random_pattern, ndimspec,           &
                                         patterngenerator_advance
@@ -25,7 +25,7 @@ module get_stochy_pattern_mod
  private
 
  public  get_random_pattern_fv3,get_random_pattern_fv3_vect
- public  get_random_pattern_sfc_fv3
+ public  get_random_pattern_fv3_sfc
  public  dump_patterns
  logical :: first_call=.true.
  contains
@@ -102,16 +102,17 @@ end subroutine get_random_pattern_fv3
 
 !>@brief The subroutine 'get_random_pattern_fv3_sfc' converts spherical harmonics to the gaussian grid then interpolates to the cubed-sphere grid once
 !>@details This subroutine is for a 2-D (lat-lon) scalar field
-subroutine get_random_pattern_sfc_fv3(rpattern,npatterns,&
+subroutine get_random_pattern_fv3_sfc(rpattern,npatterns,&
            gis_stochy,xlat,xlon,blksz,nblks,maxlen,pattern_3d)
 !\callgraph
 
 ! generate a random pattern for stochastic physics
  implicit none
- type(random_pattern), intent(inout)  :: rpattern(npatterns)
- type(stochy_internal_state), target  :: gis_stochy
+ type(random_pattern), intent(inout) :: rpattern(npatterns)
+ type(stochy_internal_state), target :: gis_stochy
  real(kind=kind_dbl_prec), intent(in) :: xlat(:,:),xlon(:,:)
  integer,intent(in)                   :: npatterns,blksz(:),nblks,maxlen
+ real(kind=kind_dbl_prec), intent(out) :: pattern_3d(nblks,maxlen,n_var_lndp)
 
  integer i,j,l,lat,ierr,n,nn,k,nt
  real(kind=kind_dbl_prec), dimension(lonf,gis_stochy%lats_node_a,1):: wrk2d
@@ -123,16 +124,16 @@ subroutine get_random_pattern_sfc_fv3(rpattern,npatterns,&
  real (kind=kind_dbl_prec)   glolal(lonf,gis_stochy%lats_node_a)
  integer kmsk0(lonf,gis_stochy%lats_node_a),len
  real(kind=kind_dbl_prec) :: globalvar,globalvar0
- real(kind=kind_dbl_prec) :: pattern_3d(nblks,maxlen,nsfcpert)
  real(kind=kind_dbl_prec) :: pattern_1d(maxlen)
  real(kind=kind_dbl_prec), allocatable, dimension(:,:) :: rslmsk
  integer :: blk
 
- do k=1,nsfcpert
+ do k=1,n_var_lndp
    kmsk0 = 0
    glolal = 0.
    do n=1,npatterns
-     if (is_master()) print *, 'Random pattern for SFC-PERTS in get_random_pattern_sfc_fv3: k, min, max ',k,minval(rpattern_sfc(n)%spec_o(:,:,k)), maxval(rpattern_sfc(n)%spec_o(:,:,k))
+     call patterngenerator_advance(rpattern(n),k,.false.)
+     if (is_master()) print *, 'Random pattern for LNDP PERTS in get_random_pattern_fv3_sfc: k, min, max ',k,minval(rpattern_sfc(n)%spec_o(:,:,k)), maxval(rpattern_sfc(n)%spec_o(:,:,k))
      call scalarspect_to_gaugrid(                       &
          rpattern(n)%spec_e(:,:,k),rpattern(n)%spec_o(:,:,k),wrk2d,&
          gis_stochy%ls_node,gis_stochy%ls_nodes,gis_stochy%max_ls_nodes,&
@@ -151,7 +152,7 @@ subroutine get_random_pattern_sfc_fv3(rpattern,npatterns,&
    enddo
 
    call mp_reduce_sum(workg,lonf,latg)
-   if (is_master()) print *, 'workg after mp_reduce_sum for SFC-PERTS in get_random_pattern_sfc_fv3: k, min, max ',k,minval(workg), maxval(workg)
+   if (is_master()) print *, 'workg after mp_reduce_sum for LNDP PERTS in get_random_pattern_fv3_sfc: k, min, max ',k,minval(workg), maxval(workg)
 
 ! interpolate to cube grid
 
@@ -166,13 +167,13 @@ subroutine get_random_pattern_sfc_fv3(rpattern,npatterns,&
       pattern_3d(blk,:,k)=pattern_1d(:)
       end associate
    enddo
-   if (is_master()) print *, '3D pattern for SFC-PERTS in get_random_pattern_sfc_fv3: k, min, max ',k,minval(pattern_3d(:,:,k)), maxval(pattern_3d(:,:,k))
+   if (is_master()) print *, '3D pattern for LNDP PERTS in get_random_pattern_fv3_sfc: k, min, max ',k,minval(pattern_3d(:,:,k)), maxval(pattern_3d(:,:,k))
    deallocate(rslmsk)
    deallocate(workg)
 
- enddo  ! loop over k, nsfcpert
+ enddo  ! loop over k, n_var_lndp
 
-end subroutine get_random_pattern_sfc_fv3
+end subroutine get_random_pattern_fv3_sfc
 
 
 !>@brief The subroutine 'get_random_pattern_fv3_vect' converts spherical harmonics to a vector on gaussian grid then interpolates to the cubed-sphere grid 
@@ -186,15 +187,15 @@ subroutine get_random_pattern_fv3_vect(rpattern,npatterns,&
  type(stochy_internal_state), target :: gis_stochy
  integer,              intent(in)    :: levs
  type(random_pattern), intent(inout) :: rpattern(npatterns)
+ real(kind=kind_dbl_prec), intent(out)  :: upattern_3d(nblks,maxlen,levs)
+ real(kind=kind_dbl_prec), intent(out) :: vpattern_3d(nblks,maxlen,levs)
 
  real(kind=kind_evod), dimension(len_trie_ls,2,1) ::  vrtspec_e,divspec_e
  real(kind=kind_evod), dimension(len_trio_ls,2,1) ::  vrtspec_o,divspec_o
  real(kind=kind_dbl_prec), intent(in) :: xlat(:,:),xlon(:,:)
  integer,intent(in)                   :: npatterns,blksz(:),nblks,maxlen
 
  integer :: blk,len
- real(kind=kind_dbl_prec) :: upattern_3d(nblks,maxlen,levs)
- real(kind=kind_dbl_prec) :: vpattern_3d(nblks,maxlen,levs)
  real(kind=kind_dbl_prec) :: pattern_1d(maxlen)
  real(kind=kind_dbl_prec), allocatable, dimension(:,:) :: rslmsk
  integer i,j,l,lat,ierr,n,nn,k,nt
@@ -393,7 +394,7 @@ subroutine dump_patterns(sfile)
     integer :: stochlun,k,n
     stochlun=99
     if (is_master()) then
-       if (nsppt > 0 .OR. nshum > 0 .OR. nskeb > 0) then
+       if (nsppt > 0 .OR. nshum > 0 .OR. nskeb > 0 .OR. nlndp > 0 ) then
           OPEN(stochlun,file=sfile,form='unformatted')
           print*,'open ',sfile,' for output'
        endif
@@ -415,6 +416,13 @@ subroutine dump_patterns(sfile)
        enddo
        enddo
     endif
+    if (nlndp > 0) then
+       do n=1,nlndp
+       do k=1,n_var_lndp
+          call write_pattern(rpattern_sfc(n),k,stochlun)
+       enddo
+       enddo
+    endif
     close(stochlun)
  end subroutine dump_patterns
 !>@brief The subroutine 'write_patterns' writes out a single pattern and the seed associated with the random number sequence