[production/RRFS.v1]Fixes for REFS restart reproducibility and crash under DEBUG

get_stochy_pattern.F90

@@ -566,7 +566,7 @@ subroutine write_pattern(rpattern,outlun,lev,np,varid1,varid2,slice_of_3d,iret)
    integer, intent(in) :: np,varid1,varid2
    logical, intent(in) :: slice_of_3d
    integer, intent(out) :: iret
-   real(kind_phys), allocatable  :: pattern2d(:)
+   real(kind_dbl_prec), allocatable  :: pattern2d(:)
    integer nm,nn,arrlen,isize,ierr
    integer,allocatable :: isave(:)
    include 'netcdf.inc'

lndp_apply_perts.F90

@@ -216,8 +216,13 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsm_noahmp, iopt_dveg
                 case('smc')
                     if (do_pert_state) then
                         p=5.
-                        min_bound = smcmin(stype(nb,i))
-                        max_bound = smcmax(stype(nb,i))
+                        if (stype(nb,i) .lt. 0.0001) then ! avoid 0 index
+                          min_bound = 0
+                          max_bound = 0
+                        else
+                          min_bound = smcmin(stype(nb,i))
+                          max_bound = smcmax(stype(nb,i))
+                        end if
 
                       ! with RUC LSM perturb smc only at time step = 2, as in HRRR
                         do k=1,lsoil
@@ -360,8 +365,16 @@ subroutine apply_pert(vname,pert,print_flag, state,ierr,p,vmin, vmax)
               print*, 'error, flat-top function requires min & max to be specified'
            endif
 
-           z = -1. + 2*(state  - vmin)/(vmax - vmin) ! flat-top function
-           state =  state  + pert*(1-abs(z**p))
+           if ((vmax-vmin) .lt. 1e-30) then ! avoid divided by zero
+             state = state  + pert
+           else
+             if (abs(state) .lt. 1e10) then ! exclude missing values
+               z = -1. + 2*(state  - vmin)/(vmax - vmin) ! flat-top function
+               state =  state  + pert*(1-abs(z**p))
+             else
+               state =  state
+             end if
+           end if
        else
            state =  state  + pert
        endif

stochy_data_mod.F90

@@ -418,13 +418,13 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret)
       if (is_rootpe()) then
          print *, 'Initialize random pattern for SPP-PERTS'
          if (stochini) then
-            ierr=NF90_INQ_VARID(stochlun,"spppert_seed", varid1)
+            ierr=NF90_INQ_VARID(stochlun,"spp_seed", varid1)
             if (ierr .NE. 0) then
                write(0,*) 'error inquring SPP-PERTS seed'
                iret = ierr
                return
             end if
-            ierr=NF90_INQ_VARID(stochlun,"ppcpert_spec", varid2)
+            ierr=NF90_INQ_VARID(stochlun,"spp_spec", varid2)
             if (ierr .NE. 0) then
                write(0,*) 'error inquring SPP-PERTS spec'
                iret = ierr
@@ -702,7 +702,7 @@ subroutine read_pattern(rpattern,jcapin,lunptn,k,np,varid1,varid2,slice_of_3d,ir
    type(random_pattern), intent(inout) :: rpattern
    integer, intent(in) :: lunptn,np,varid1,varid2,jcapin
    logical, intent(in) :: slice_of_3d
-   real(kind_phys),allocatable  :: pattern2d(:),pattern2din(:)
+   real(kind_dbl_prec),allocatable  :: pattern2d(:),pattern2din(:)
    real(kind_phys) :: stdevin,varin
    integer nm,nn,iret,ierr,isize,k,ndimspec2
    integer, allocatable :: isave(:)

stochy_patterngenerator.F90

@@ -220,7 +220,7 @@ end subroutine computevarspec
  subroutine computevarspec_r(rpattern,dataspec,var)
 !\callgraph
     ! compute globally integrated variance from spectral coefficients
-    real(kind_phys), intent(in) :: dataspec(2*ndimspec)
+    real(kind_dbl_prec), intent(in) :: dataspec(2*ndimspec)
     real(kind_phys), intent(out) ::  var
     type(random_pattern), intent(in) :: rpattern
     integer n
@@ -364,8 +364,8 @@ end subroutine setvarspect
 !! restarting from a higher-resolution pattern
  subroutine chgres_pattern(pattern2din,pattern2dout,ntruncin,ntruncout)
 !\callgraph
-   real(kind_phys), intent(in) :: pattern2din((ntruncin+1)*(ntruncin+2))
-   real(kind_phys), intent(out) :: pattern2dout((ntruncout+1)*(ntruncout+2))
+   real(kind_dbl_prec), intent(in) :: pattern2din((ntruncin+1)*(ntruncin+2))
+   real(kind_dbl_prec), intent(out) :: pattern2dout((ntruncout+1)*(ntruncout+2))
    integer, intent(in) :: ntruncin,ntruncout
    integer             :: m,n,nm,ndimsspecin,ndimsspecout
    integer,allocatable, dimension(:,:):: idxin