#if ( HYBRID_COORD==1 ) 
#  define mu(...) (c1h(k)*XXPCXX(__VA_ARGS__))
#  define XXPCXX(...) mu(__VA_ARGS__)

#  define mub(...) (c1h(k)*XXPCBXX(__VA_ARGS__)+c2h(k))
#  define XXPCBXX(...) mub(__VA_ARGS__)
#endif

module module_stoch
!***********************************************************************
!
! Purpose: Stochastic Perturbation Schemes
! Author : Judith Berner, NCAR  (berner@ucar.edu)
! Date   : Apr 2014
!
!***********************************************************************
!
! The scheme introduces stochastic perturbations to the rotational wind 
! components and to the potential temperature field. The stochastic 
! perturbations are generated by independent autoregressive processes for 
! each wavenumber and results in smooth spatially and temporally correlated patterns.
! Details of the scheme and its performance in a meso-scale WRF-ensemble 
! system are available in:
!
! Berner, J.,  S.-Y. Ha, J. P. Hacker, A. Fournier and C. Snyder 2011:
! "Model uncertainty in a mesoscale ensemble prediction system: Stochastic 
! versus multi-physics representations",  2011, Mon. Wea. Rev., 139, 1972-1995
! http://journals.ametsoc.org/doi/abs/10.1175/2010MWR3595.1
!
! Features: 
!     Dissipation: Dissipation rates are assumed constant in space and time
!     Vertical structure: Supports two options for vertical structure: 
!                         0) constant
!                         1) random phase
!
! Optional namelist parameters:
!   stoch_force_opt		=	0, 0, 0: No stochastic parameterization
!   			=	1, 1, 1: Use SKEB scheme
!   skebs_vertstruc	=	0, 0, 0: Constant vertical structure of random pattern generator
!   			= 	1, 1, 1: Random phase vertical structure random pattern generator
!   tot_backscat_psi	:	Total backscattered dissipation rate for streamfunction; Controls 
!                                   amplitude of rotational wind perturbations Default value is 1.0E-5 m2/s3.
!   tot_backscat_t		:	Total backscattered dissipation rate for potential temperature; 
!   				Controls amplitude of potential temperature perturbations. Default value is 1.0E-6 m2/s3.
!   nens			:	Random seed for random number stream. This parameter needs to be different 
!                                   for each member in ensemble forecasts. Is a function of initial start time 
!   				to ensure different random number streams for different forecasts.
!   ztau_psi		:	Decorrelation time (in s) for streamfunction perturbations.
!   				Default is 10800s.  Recommended value is 216000s.
!   ztau_t			:	Decorrelation time (in s) for potential temperature perturbations. 
!   				Default 10800s. Recommended value is 216000s.
!   rexponent_psi		:	Spectral slope for streamfunction perturbations. Default is -1.83 
!                                   for a kinetic-energy forcing spectrum with slope -5/3.
!   rexponent_t 		:	Spectral slope of potential temperature perturbations. Default is -1.83 
!   				for a potential energy forcing spectrum with slope -1.832.
!   kminforc		:	Minimal forcing wavenumber in longitude for streamfunction perturbations. Default is 1.
!   lminforc		:	Minimal forcing wavenumber in latitude for streamfunction perturbations. Default is 1.
!   kminforc		: 	Minimal forcing wavenumber in longitude for potential temperature perturbations. Default is 1.
!   lminforct		:	Minimal forcing wavenumber in latitude for potential temperature perturbations. Default is 1.
!   kmaxforc		:	Maximal forcing wavenumber in longitude for streamfunction perturbations. 
!   				Default is maximal possible wavenumbers determined by number of gridpoints.
!   lmaxforc		:	Maximal forcing wavenumber in latitude for streamfunction perturbations. 
!   				Default is maximal possible wavenumbers determined by number of gridpoints.
!   kmaxforct		:	Maximal forcing wavenumber in longitude for potential temperature perturbations. 
!   				Default is maximal possible wavenumbers determined by number of gridpoints.
!   lmaxforct		:	Maximal forcing wavenumber in latitude for potential temperature perturbations. 
!   				Default is maximal possible wavenumbers determined by number of gridpoints.
!   zsigma2_eps		:	Noise variance in autoregressive process defining streamfunction perturbations.
!   zsigma2_eta		:	Noise variance in autoregressive process defining in potential temperature perturbations.

!***********************************************************************
!     ------------------------------------------------------------------
!************** DECLARE FIELDS AND VARIABLES FOR STOCHASTIC BACKSCATTER
!     ------------------------------------------------------------------
      implicit none
      public ::  SETUP_RAND_PERTURB, UPDATE_STOCH,& 
                         rand_pert_update

      INTEGER :: LMINFORC, LMAXFORC, KMINFORC, KMAXFORC, & 
      &          LMINFORCT, LMAXFORCT, KMINFORCT, KMAXFORCT
      REAL    :: ALPH, ALPH_PSI, ALPH_T, TOT_BACKSCAT_PSI, TOT_BACKSCAT_T,  REXPONENT_PSI,REXPONENT_T

!     ----------Fields for spectral transform -----------

      INTEGER :: LENSAV
      INTEGER,ALLOCATABLE:: wavenumber_k(:), wavenumber_l(:)
      REAL, ALLOCATABLE :: WSAVE1(:),WSAVE2(:), WSAVE1_ideal(:),WSAVE2_ideal(:)
      REAL, ALLOCATABLE :: rindarrayik(:),rindarrayil(:),rindarrayl(:),rindarrayk(:)

!     --------- Others -------------------------------------------------
      REAL, PARAMETER:: RPI= 3.141592653589793 !4.0*atan(1.0) 
      REAL, PARAMETER:: CP= 1006.0 ! specific heat of dry air in J/(Kg*K)= m^2/(K* s^2)
      REAL, PARAMETER:: T0= 300.0 ! Reference temperature in K 
      character(len=33)                         :: filenamesave

      save


!=======================================================================
contains
!=======================================================================
!     ------------------------------------------------------------------
!!************** INITIALIZE STOCHASTIC ROUTINES *****************************
!     ------------------------------------------------------------------
! This subroutine drives the initialization of the stochastic schemes  

      SUBROUTINE INITIALIZE_STOCH  (grid, config_flags,                       &
                          first_trip_for_this_domain,                         & 
                          ips, ipe, jps, jpe, kps, kpe,                       &
                          ids, ide, jds, jde, kds, kde,                       &
                          ims, ime, jms, jme, kms, kme,                       &
                          its, ite, jts, jte, kts, kte,                       & 
                          imsx, imex, jmsx, jmex, kmsx, kmex,  &
                          ipsx, ipex, jpsx, jpex, kpsx, kpex,  &
                          imsy, imey, jmsy, jmey, kmsy, kmey,  &
                          ipsy, ipey, jpsy, jpey, kpsy, kpey )


    USE module_configure
    USE module_domain, ONLY : domain

      IMPLICIT NONE

      TYPE (grid_config_rec_type)            :: config_flags
      TYPE ( domain ), INTENT(INOUT)         :: grid

      INTEGER , INTENT(IN)     ::               ids, ide, jds, jde, kds, kde,   &
                                                ims, ime, jms, jme, kms, kme,   &
                                                ips, ipe, jps, jpe, kps, kpe,   & 
                                                its, ite, jts, jte, kts, kte
      INTEGER , INTENT(IN)     ::               imsx,imex,jmsx,jmex,kmsx,kmex, &
                                                ipsx,ipex,jpsx,jpex,kpsx,kpex, &
                                                imsy,imey,jmsy,jmey,kmsy,kmey, &
                                                ipsy,ipey,jpsy,jpey,kpsy,kpey

      LOGICAL                  ::               first_trip_for_this_domain
      INTEGER, PARAMETER       ::               maxspinup=100
      INTEGER                  ::               K 


   IF ( first_trip_for_this_domain ) THEN
     grid%did_stoch = .FALSE.
   END IF


   IF ((( grid%id == 1) .AND. (.NOT. grid%did_stoch)) .AND. &
       (( grid%skebs_on== 1) .OR.( grid%sppt_on== 1) .OR. ( grid%rand_perturb_on== 1) .OR. &
         ( grid%spp_conv== 1) .OR. ( grid%spp_pbl== 1) .OR. ( grid%spp_mp .GE. 1) .OR. &
         ( grid%spp_lsm== 1)) ) THEN

     grid%did_stoch = .TRUE.

     IF (grid%skebs_on==1) then

! Initialize SKEBS
!    Initialize streamfunction (1)
     if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then 
         call rand_seed (config_flags, grid%ISEED_SKEBS, grid%iseedarr_skebs , 1, config_flags%seed_dim)
     else
         call read_write_stochrestart ('stoch_skebs_U',grid,config_flags,grid%lmax_ideal,grid%kmax_ideal, & 
                                       grid%SPSTREAMFORCS,grid%SPSTREAMFORCC,2)
         call read_write_stochrestart ('stoch_skebs_T',grid,config_flags,grid%lmax_ideal,grid%kmax_ideal, & 
                                       grid%SPTFORCS,grid%SPTFORCC,2)
         ! broadcast SPFORC and SPFORS to all processes
         #ifdef DM_PARALLEL
           CALL wrf_dm_bcast_real ( grid%SPSTREAMFORCS ,grid%kmax_ideal*grid%lmax_ideal )
           CALL wrf_dm_bcast_real ( grid%SPSTREAMFORCC ,grid%kmax_ideal*grid%lmax_ideal )
           CALL wrf_dm_bcast_real ( grid%SPTFORCS,grid%kmax_ideal*grid%lmax_ideal )
           CALL wrf_dm_bcast_real ( grid%SPTFORCC,grid%kmax_ideal*grid%lmax_ideal )
         #endif
     endif
     call SETUP_RAND_PERTURB('W',                                         &
                       grid%skebs_vertstruc,config_flags%restart,         &
                       grid%SPSTREAM_AMP,                                 &
                       grid%SPSTREAMFORCS,grid%SPSTREAMFORCC,grid%ALPH_PSI,&
                       grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPUV,     &
                       grid%KMINFORCT,grid%KMAXFORCT,                     &
                       grid%LMINFORCT,grid%LMAXFORCT,                     &
                       grid%KMAXFORCTH,grid%LMAXFORCTH,                   &
                       grid%time_step,grid%DX,grid%DY,                    &
                       grid%stepstoch,                                     & 
                       grid%gridpt_stddev_sppt,                           &
                       grid%lengthscale_sppt,                             &
                       grid%timescale_sppt,                               &
                       grid%TOT_BACKSCAT_PSI,grid%ZTAU_PSI,               &
                       grid%REXPONENT_PSI,                                &
                       grid%kmax_ideal, grid%lmax_ideal,                  & 
                       ids, ide, jds, jde, kds, kde,                      &
                       ims, ime, jms, jme, kms, kme,                      &
                       its, ite, jts, jte, kts, kte                       )
     if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then 
        do k = 1,maxspinup
          CALL RAND_PERT_UPDATE(grid,'U',                                    &
                           grid%SPSTREAMFORCS,grid%SPSTREAMFORCC,             &
                           grid%SPSTREAM_AMP,grid%ALPH_PSI,                   &
                           ips, ipe, jps, jpe, kps, kpe,                      &
                           ids, ide, jds, jde, kds, kde,                      &
                           ims, ime, jms, jme, kms, kme,                      &
                           kts,kte,                                    &
                           imsx,imex,jmsx,jmex,kmsx,kmex,                     &
                           ipsx,ipex,jpsx,jpex,kpsx,kpex,                     &
                           imsy,imey,jmsy,jmey,kmsy,kmey,                     &
                           ipsy,ipey,jpsy,jpey,kpsy,kpey,                     &
                           grid%kmax_ideal, grid%lmax_ideal,                  &
                           grid% num_stoch_levels,grid% num_stoch_levels,     &
                           grid% num_stoch_levels,grid% num_stoch_levels,     &
                           config_flags%restart, grid%iseedarr_skebs,         &
                           config_flags%seed_dim,                              &
                           grid%DX,grid%DY,grid%skebs_vertstruc,              &
                           grid%ru_tendf_stoch,                               &
                           grid%stddev_cutoff_sppt,grid%gridpt_stddev_sppt, &
                           grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPUV     )
          CALL RAND_PERT_UPDATE(grid,'V',                                    &
                           grid%SPSTREAMFORCS,grid%SPSTREAMFORCC,             &
                           grid%SPSTREAM_AMP,grid%ALPH_PSI,                   &
                           ips, ipe, jps, jpe, kps, kpe,                      &
                           ids, ide, jds, jde, kds, kde,                      &
                           ims, ime, jms, jme, kms, kme,                      &
                           kts,kte,                                    &
                           imsx,imex,jmsx,jmex,kmsx,kmex,                     &
                           ipsx,ipex,jpsx,jpex,kpsx,kpex,                     &
                           imsy,imey,jmsy,jmey,kmsy,kmey,                     &
                           ipsy,ipey,jpsy,jpey,kpsy,kpey,                     &
                           grid%kmax_ideal, grid%lmax_ideal,                  &
                           grid% num_stoch_levels,grid% num_stoch_levels,     &
                           grid% num_stoch_levels,grid% num_stoch_levels,     &
                           config_flags%restart, grid%iseedarr_skebs,         &
                           config_flags%seed_dim,                              &
                           grid%DX,grid%DY,grid%skebs_vertstruc,              &
                           grid%rv_tendf_stoch,                               &
                           grid%stddev_cutoff_sppt,grid%gridpt_stddev_sppt, &
                           grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT      )

         enddo
       endif
!    Initialize potential temperature (2)
     call SETUP_RAND_PERTURB('T',                                         &
                       grid%skebs_vertstruc,config_flags%restart,     &
                       grid%SPT_AMP,                                      &
                       grid%SPTFORCS,grid%SPTFORCC,grid%ALPH_T,           &
                       grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT,     &
                       grid%KMINFORCT,grid%KMAXFORCT,                     &
                       grid%LMINFORCT,grid%LMAXFORCT,                     &
                       grid%KMAXFORCTH,grid%LMAXFORCTH,                   &
                       grid%time_step,grid%DX,grid%DY,                    &
                       grid%stepstoch,                                     & 
                       grid%gridpt_stddev_sppt,                           &
                       grid%lengthscale_sppt,                             &
                       grid%timescale_sppt,                               &
                       grid%TOT_BACKSCAT_T,grid%ZTAU_T,                   &
                       grid%REXPONENT_T,                                  &
                       grid%kmax_ideal, grid%lmax_ideal,                  & 
                       ids, ide, jds, jde, kds, kde,                      &
                       ims, ime, jms, jme, kms, kme,                      &
                       its, ite, jts, jte, kts, kte                       )
     if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then 
        do k = 1,maxspinup
           CALL RAND_PERT_UPDATE(grid,'T',                                     &
                          grid%SPTFORCS,grid%SPTFORCC,                        &
                          grid%SPT_AMP,grid%ALPH_T,                           &
                          ips, ipe, jps, jpe, kps, kpe,                       &
                          ids, ide, jds, jde, kds, kde,                       &
                          ims, ime, jms, jme, kms, kme,                       &
                          kts,kte,                                     &
                          imsx,imex,jmsx,jmex,kmsx,kmex,                      &
                          ipsx,ipex,jpsx,jpex,kpsx,kpex,                      &
                          imsy,imey,jmsy,jmey,kmsy,kmey,                      &
                          ipsy,ipey,jpsy,jpey,kpsy,kpey,                      &
                          grid%kmax_ideal, grid%lmax_ideal,                   &
                          grid%num_stoch_levels,grid%num_stoch_levels,        &
                          grid%num_stoch_levels,grid%num_stoch_levels,        &
                          config_flags%restart, grid%iseedarr_skebs,          &
                          config_flags%seed_dim,                              &
                          grid%DX,grid%DY,grid%skebs_vertstruc,               &
                          grid%rt_tendf_stoch,                                &
                          grid%stddev_cutoff_sppt,grid%gridpt_stddev_sppt, &
                          grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPUV      )

         enddo
       endif
     ENDIF

IF (grid%sppt_on==1) then
! Initialize SPPT (3)
     if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then 
         call rand_seed (config_flags, grid%ISEED_SPPT, grid%iseedarr_sppt  ,  1, config_flags%seed_dim)
     else
         call read_write_stochrestart ('stoch_sppt_xx',grid,config_flags,grid%lmax_ideal,grid%kmax_ideal, & 
                                       grid%SPPTFORCS,grid%SPPTFORCC,2)
         ! broadcast SPFORC and SPFORS to all processes
         #ifdef DM_PARALLEL
           CALL wrf_dm_bcast_real ( grid%SPPTFORCS,grid%kmax_ideal*grid%lmax_ideal )
           CALL wrf_dm_bcast_real ( grid%SPPTFORCC,grid%kmax_ideal*grid%lmax_ideal )
         #endif
     endif
     call SETUP_RAND_PERTURB('P',                                         &
                       grid%sppt_vertstruc,config_flags%restart,          &
                       grid%SPPT_AMP,                                     &
                       grid%SPPTFORCC,grid%SPPTFORCS,grid%ALPH_SPPT,      &
                       grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT,     &
                       grid%KMINFORCT,grid%KMAXFORCT,                     &
                       grid%LMINFORCT,grid%LMAXFORCT,                     &
                       grid%KMAXFORCTH,grid%LMAXFORCTH,                   &
                       grid%time_step,grid%DX,grid%DY,                    &
                       grid%stepstoch,                                     & 
                       grid%gridpt_stddev_sppt,                           &
                       grid%lengthscale_sppt,                             &
                       grid%timescale_sppt,                               &
                       grid%TOT_BACKSCAT_PSI,grid%ZTAU_PSI,               &
                       grid%REXPONENT_PSI,                                &
                       grid%kmax_ideal, grid%lmax_ideal,                  & 
                       ids, ide, jds, jde, kds, kde,                      &
                       ims, ime, jms, jme, kms, kme,                      &
                       its, ite, jts, jte, kts, kte                       )
     if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then 
        do k = 1,maxspinup
          CALL RAND_PERT_UPDATE(grid,'T',                                     &
                          grid%SPPTFORCS,grid%SPPTFORCC,                      &
                          grid%SPPT_AMP,grid%ALPH_SPPT,                       &
                          ips, ipe, jps, jpe, kps, kpe,                       &
                          ids, ide, jds, jde, kds, kde,                       &
                          ims, ime, jms, jme, kms, kme,                       &
                          kts,kte,                                     &
                          imsx,imex,jmsx,jmex,kmsx,kmex,                      &
                          ipsx,ipex,jpsx,jpex,kpsx,kpex,                      &
                          imsy,imey,jmsy,jmey,kmsy,kmey,                      &
                          ipsy,ipey,jpsy,jpey,kpsy,kpey,                      &
                          grid%kmax_ideal, grid%lmax_ideal,                   &
                          grid%num_stoch_levels,grid%num_stoch_levels,        &
                          grid%num_stoch_levels,grid%num_stoch_levels,        &
                          config_flags%restart, grid%iseedarr_sppt,           &
                          config_flags%seed_dim,                              &
                          grid%DX,grid%DY,grid%sppt_vertstruc,                &
                          grid%rstoch,                                        &
                          grid%stddev_cutoff_sppt,grid%gridpt_stddev_sppt, &
                          grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT       )
         enddo
       endif
     ENDIF

! Initialize RAND_PERTURB (4)
     IF (grid%rand_perturb_on==1) then
     if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then 
         call rand_seed (config_flags, grid%ISEED_RAND_PERT, grid%iseedarr_rand_pert  ,  1, config_flags%seed_dim)
     else

         call read_write_stochrestart ('stoch_ranpert',grid,config_flags,grid%lmax_ideal,grid%kmax_ideal, & 
                                       grid%SPFORCS,grid%SPFORCC,2)
         ! broadcast SPFORC and SPFORS to all processes
         #ifdef DM_PARALLEL
           CALL wrf_dm_bcast_real ( grid%SPFORCS , grid%kmax_ideal*grid%lmax_ideal )
           CALL wrf_dm_bcast_real ( grid%SPFORCC , grid%kmax_ideal*grid%lmax_ideal )
         #endif
     endif
     call SETUP_RAND_PERTURB('R',                                         &
                       grid%rand_pert_vertstruc,config_flags%restart,     &
                       grid%SP_AMP,                                       &
                       grid%SPFORCC,grid%SPFORCS,grid%ALPH_RAND,          &
                       grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT,     &
                       grid%KMINFORCT,grid%KMAXFORCT,                     &
                       grid%LMINFORCT,grid%LMAXFORCT,                     &
                       grid%KMAXFORCTH,grid%LMAXFORCTH,                   &
                       grid%time_step,grid%DX,grid%DY,                    &
                       grid%stepstoch,                                     & 
                       grid%gridpt_stddev_rand_pert,                      &
                       grid%lengthscale_rand_pert,                        &
                       grid%timescale_rand_pert,                          &
                       grid%TOT_BACKSCAT_PSI,grid%ZTAU_PSI,               &
                       grid%REXPONENT_PSI,                                &
                       grid%kmax_ideal, grid%lmax_ideal,                  & 
                       ids, ide, jds, jde, kds, kde,                      &
                       ims, ime, jms, jme, kms, kme,                      &
                       its, ite, jts, jte, kts, kte                       )

     if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then 
        do k = 1,maxspinup
           CALL RAND_PERT_UPDATE(grid,'T',                                     &
                           grid%SPFORCS,grid%SPFORCC,                          &
                           grid%SP_AMP,grid%ALPH_RAND,                         &
                           ips, ipe, jps, jpe, kps, kpe,                       &
                           ids, ide, jds, jde, kds, kde,                       &
                           ims, ime, jms, jme, kms, kme,                       &
                           kts, kte,                                           &
                           imsx,imex,jmsx,jmex,kmsx,kmex,                      &
                           ipsx,ipex,jpsx,jpex,kpsx,kpex,                      &
                           imsy,imey,jmsy,jmey,kmsy,kmey,                      &
                           ipsy,ipey,jpsy,jpey,kpsy,kpey,                      &
                           grid%kmax_ideal, grid%lmax_ideal,                  & 
                           grid%num_stoch_levels,grid%num_stoch_levels,        &
                           grid%num_stoch_levels,grid%num_stoch_levels,        &
                           config_flags%restart, grid%iseedarr_rand_pert,      &
                           config_flags%seed_dim,                              &
                           grid%DX,grid%DY,grid%rand_pert_vertstruc,           &
                           grid%RAND_PERT,                                     &
                           grid%gridpt_stddev_rand_pert,                      &
                           grid%lengthscale_rand_pert,                        &
                           grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT       )
         enddo
       ENDIF !rand_perturb_on
     ENDIF

! Initialize Stochastic Parameter Perturbations to convection scheme
     IF (grid%spp_conv==1) then
     if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then 
         call rand_seed (config_flags, grid%iseed_spp_conv, grid%iseedarr_spp_conv ,  1, config_flags%seed_dim)
     else
         call read_write_stochrestart ('stoch_spp_con',grid,config_flags,grid%lmax_ideal,grid%kmax_ideal, & 
                                       grid%SPFORCS2,grid%SPFORCC2,2)
         ! broadcast SPFORC and SPFORS to all processes
         #ifdef DM_PARALLEL
           CALL wrf_dm_bcast_real ( grid%SPFORCS2 ,grid%kmax_ideal*grid%lmax_ideal )
           CALL wrf_dm_bcast_real ( grid%SPFORCC2 ,grid%kmax_ideal*grid%lmax_ideal )
         #endif
     endif
     call SETUP_RAND_PERTURB('R',                                        &
                       grid%vertstruc_spp_conv,config_flags%restart,    &
                       grid%SP_AMP2,                                      &
                       grid%SPFORCC2,grid%SPFORCS2,grid%ALPH_RAND2,       &
                       grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT,     &
                       grid%KMINFORCT,grid%KMAXFORCT,                     &
                       grid%LMINFORCT,grid%LMAXFORCT,                     &
                       grid%KMAXFORCTH,grid%LMAXFORCTH,                   &
                       grid%time_step,grid%DX,grid%DY,                    &
                       grid%stepstoch,                                     & 
                       grid%gridpt_stddev_spp_conv,                     &
                       grid%lengthscale_spp_conv,                       &
                       grid%timescale_spp_conv,                         &
                       grid%TOT_BACKSCAT_PSI,grid%ZTAU_PSI,               &
                       grid%REXPONENT_PSI,                                &
                       grid%kmax_ideal, grid%lmax_ideal,                  & 
                       ids, ide, jds, jde, kds, kde,                      &
                       ims, ime, jms, jme, kms, kme,                      &
                       its, ite, jts, jte, kts, kte                       )
     if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then 
        do k = 1,maxspinup
           CALL RAND_PERT_UPDATE(grid,'T',                                     &
                           grid%SPFORCS2,grid%SPFORCC2,                        &
                           grid%SP_AMP2,grid%ALPH_RAND2,                       &
                           ips, ipe, jps, jpe, kps, kpe,                       &
                           ids, ide, jds, jde, kds, kde,                       &
                           ims, ime, jms, jme, kms, kme,                       &
                           kts, kte,                                           &
                           imsx,imex,jmsx,jmex,kmsx,kmex,                      &
                           ipsx,ipex,jpsx,jpex,kpsx,kpex,                      &
                           imsy,imey,jmsy,jmey,kmsy,kmey,                      &
                           ipsy,ipey,jpsy,jpey,kpsy,kpey,                      &
                           grid%kmax_ideal, grid%lmax_ideal,                   &
                           grid%num_stoch_levels,grid%num_stoch_levels,        &
                           grid%num_stoch_levels,grid%num_stoch_levels,        &
                           config_flags%restart, grid%iseedarr_spp_conv,       &
                           config_flags%seed_dim,                              &
                           grid%DX,grid%DY,grid%vertstruc_spp_conv,            &
                           grid%pattern_spp_conv,                              &
                           grid%stddev_cutoff_spp_conv,grid%gridpt_stddev_spp_conv, &
                           grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT       )
         enddo
       endif
     ENDIF

! Initialize Stochastic Parameter Peturbations (SPP) to PBL scheme
     IF (grid%spp_pbl==1) then
     if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then 
         call rand_seed (config_flags, grid%iseed_spp_pbl, grid%iseedarr_spp_pbl ,  1, config_flags%seed_dim)
     else
         call read_write_stochrestart ('stoch_spp_pbl',grid,config_flags,grid%lmax_ideal,grid%kmax_ideal, & 
                                       grid%SPFORCS3,grid%SPFORCC3,2)
         ! broadcast SPFORC and SPFORS to all processes
         #ifdef DM_PARALLEL
           CALL wrf_dm_bcast_real ( grid%SPFORCS3,grid%kmax_ideal*grid%lmax_ideal )
           CALL wrf_dm_bcast_real ( grid%SPFORCC3,grid%kmax_ideal*grid%lmax_ideal )
         #endif
     endif
     call SETUP_RAND_PERTURB('R',                                        &
                       grid%vertstruc_spp_pbl,config_flags%restart,    &
                       grid%SP_AMP3,                                      &
                       grid%SPFORCC3,grid%SPFORCS3,grid%ALPH_RAND3,       &
                       grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT,     &
                       grid%KMINFORCT,grid%KMAXFORCT,                     &
                       grid%LMINFORCT,grid%LMAXFORCT,                     &
                       grid%KMAXFORCTH,grid%LMAXFORCTH,                   &
                       grid%time_step,grid%DX,grid%DY,                    &
                       grid%stepstoch,                                     & 
                       grid%gridpt_stddev_spp_pbl,                     &
                       grid%lengthscale_spp_pbl,                       &
                       grid%timescale_spp_pbl,                         &
                       grid%TOT_BACKSCAT_PSI,grid%ZTAU_PSI,               &
                       grid%REXPONENT_PSI,                                &
                       grid%kmax_ideal, grid%lmax_ideal,                  & 
                       ids, ide, jds, jde, kds, kde,                      &
                       ims, ime, jms, jme, kms, kme,                      &
                       its, ite, jts, jte, kts, kte                       )
     if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then 
        do k = 1,maxspinup
           CALL RAND_PERT_UPDATE(grid,'T',                                     &
                           grid%SPFORCS3,grid%SPFORCC3,                        &
                           grid%SP_AMP3,grid%ALPH_RAND3,                       &
                           ips, ipe, jps, jpe, kps, kpe,                       &
                           ids, ide, jds, jde, kds, kde,                       &
                           ims, ime, jms, jme, kms, kme,                       &
                           kts, kte,                                           &
                           imsx,imex,jmsx,jmex,kmsx,kmex,                      &
                           ipsx,ipex,jpsx,jpex,kpsx,kpex,                      &
                           imsy,imey,jmsy,jmey,kmsy,kmey,                      &
                           ipsy,ipey,jpsy,jpey,kpsy,kpey,                      &
                           grid%kmax_ideal, grid%lmax_ideal,                   &
                           grid%num_stoch_levels,grid%num_stoch_levels,        &
                           grid%num_stoch_levels,grid%num_stoch_levels,        &
                           config_flags%restart, grid%iseedarr_spp_pbl,        &
                           config_flags%seed_dim,                              &
                           grid%DX,grid%DY,grid%vertstruc_spp_pbl,             &
                           grid%pattern_spp_pbl,                               &
                           grid%stddev_cutoff_spp_pbl,grid%gridpt_stddev_spp_pbl, &
                           grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT       )

         enddo
       endif
     ENDIF
! Initialize Stochastic Parameter Peturbations (SPP) to MP Thompson scheme

     IF (grid%spp_mp .GE. 1) then 
     if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then 
         call rand_seed (config_flags, grid%iseed_spp_mp, grid%iseedarr_spp_mp,  1, config_flags%seed_dim)
     else
         call read_write_stochrestart ('stoch_spp_mph',grid,config_flags,grid%lmax_ideal,grid%kmax_ideal, & 
                                       grid%SPFORCS5,grid%SPFORCC5,2)
         ! broadcast SPFORC and SPFORS to all processes
         #ifdef DM_PARALLEL
           CALL wrf_dm_bcast_real ( grid%SPFORCS5,grid%kmax_ideal*grid%lmax_ideal )
           CALL wrf_dm_bcast_real ( grid%SPFORCC5,grid%kmax_ideal*grid%lmax_ideal )
         #endif
     endif
     call SETUP_RAND_PERTURB('R',                                        &    
                       grid%vertstruc_spp_mp,config_flags%restart,    &    
                       grid%SP_AMP5,                                      &    
                       grid%SPFORCC5,grid%SPFORCS5,grid%ALPH_RAND5,       &    
                       grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT,     &    
                       grid%KMINFORCT,grid%KMAXFORCT,                     &    
                       grid%LMINFORCT,grid%LMAXFORCT,                     &    
                       grid%KMAXFORCTH,grid%LMAXFORCTH,                   &    
                       grid%time_step,grid%DX,grid%DY,                    &    
                       grid%stepstoch,                                     & 
                       grid%gridpt_stddev_spp_mp,                     &    
                       grid%lengthscale_spp_mp,                       &    
                       grid%timescale_spp_mp,                         &    
                       grid%TOT_BACKSCAT_PSI,grid%ZTAU_PSI,               &    
                       grid%REXPONENT_PSI,                                &    
                       grid%kmax_ideal, grid%lmax_ideal,                  & 
                       ids, ide, jds, jde, kds, kde,                      &    
                       ims, ime, jms, jme, kms, kme,                      &    
                       its, ite, jts, jte, kts, kte                       )    
     if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then 
        do k = 1,maxspinup
           CALL RAND_PERT_UPDATE(grid,'T',                                     &
                           grid%SPFORCS5,grid%SPFORCC5,                          &
                           grid%SP_AMP5,grid%ALPH_RAND5,                         &
                           ips, ipe, jps, jpe, kps, kpe,                       &
                           ids, ide, jds, jde, kds, kde,                       &
                           ims, ime, jms, jme, kms, kme,                       &
                           kts, kte,                                           &
                           imsx,imex,jmsx,jmex,kmsx,kmex,                      &
                           ipsx,ipex,jpsx,jpex,kpsx,kpex,                      &
                           imsy,imey,jmsy,jmey,kmsy,kmey,                      &
                           ipsy,ipey,jpsy,jpey,kpsy,kpey,                      &
                           grid%kmax_ideal, grid%lmax_ideal,                   &
                           grid%num_stoch_levels,grid%num_stoch_levels,        &
                           grid%num_stoch_levels,grid%num_stoch_levels,        &
                           config_flags%restart, grid%iseedarr_spp_mp,        &
                           config_flags%seed_dim,                              &
                           grid%DX,grid%DY,grid%vertstruc_spp_mp,             &
                           grid%pattern_spp_mp,                               &
                           grid%stddev_cutoff_spp_mp,grid%gridpt_stddev_spp_mp,&
                           grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT       )
         enddo
       endif 
     ENDIF

! Initialize Stochastic Parameter Peturbations (SPP) to LSM scheme
     IF (grid%spp_lsm==1) then
     if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then 
         call rand_seed (config_flags, grid%iseed_spp_lsm, grid%iseedarr_spp_lsm ,  1, config_flags%seed_dim)
     else
         call read_write_stochrestart ('stoch_spp_lsm',grid,config_flags,grid%lmax_ideal,grid%kmax_ideal, & 
                                       grid%SPFORCS4,grid%SPFORCC4,2)
         ! broadcast SPFORC and SPFORS to all processes
         #ifdef DM_PARALLEL
           CALL wrf_dm_bcast_real ( grid%SPFORCS4,grid%kmax_ideal*grid%lmax_ideal )
           CALL wrf_dm_bcast_real ( grid%SPFORCC4,grid%kmax_ideal*grid%lmax_ideal )
         #endif
     endif
     call SETUP_RAND_PERTURB('R',                                        &
                       grid%vertstruc_spp_lsm,config_flags%restart,    &
                       grid%SP_AMP4,                                      &
                       grid%SPFORCC4,grid%SPFORCS4,grid%ALPH_RAND4,       &
                       grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT,     &
                       grid%KMINFORCT,grid%KMAXFORCT,                     &
                       grid%LMINFORCT,grid%LMAXFORCT,                     &
                       grid%KMAXFORCTH,grid%LMAXFORCTH,                   &
                       grid%time_step,grid%DX,grid%DY,                    &
                       grid%stepstoch,                                     & 
                       grid%gridpt_stddev_spp_lsm,                     &
                       grid%lengthscale_spp_lsm,                       &
                       grid%timescale_spp_lsm,                         &
                       grid%TOT_BACKSCAT_PSI,grid%ZTAU_PSI,               &
                       grid%REXPONENT_PSI,                                &
                       grid%kmax_ideal, grid%lmax_ideal,                  & 
                       ids, ide, jds, jde, kds, kde,                      &
                       ims, ime, jms, jme, kms, kme,                      &
                       its, ite, jts, jte, kts, kte                       )
     if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then 
        do k = 1,maxspinup
           CALL RAND_PERT_UPDATE(grid,'T',                                     &
                           grid%SPFORCS4,grid%SPFORCC4,                          &
                           grid%SP_AMP4,grid%ALPH_RAND4,                         &
                           ips, ipe, jps, jpe, kps, kpe,                       &
                           ids, ide, jds, jde, kds, kde,                       &
                           ims, ime, jms, jme, kms, kme,                       &
                           kts, kte,                                           &
                           imsx,imex,jmsx,jmex,kmsx,kmex,                      &
                           ipsx,ipex,jpsx,jpex,kpsx,kpex,                      &
                           imsy,imey,jmsy,jmey,kmsy,kmey,                      &
                           ipsy,ipey,jpsy,jpey,kpsy,kpey,                      &
                           grid%kmax_ideal, grid%lmax_ideal,                   &
                           grid%num_stoch_levels,grid%num_stoch_levels,        &
                           grid%num_stoch_levels,grid%num_stoch_levels,        &
                           config_flags%restart, grid%iseedarr_spp_lsm,        &
                           config_flags%seed_dim,                              &
                           grid%DX,grid%DY,grid%vertstruc_spp_lsm,             &
                           grid%pattern_spp_lsm,                               &
                           grid%stddev_cutoff_spp_lsm,grid%gridpt_stddev_spp_lsm,&
                           grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT       )
         enddo
       endif 
     ENDIF

     ENDIF ! skebs or sppt or rand_perturb or spp
 
      END SUBROUTINE INITIALIZE_STOCH       

      !   --- END SETUP STOCHASTIC PERTURBATION SCHEMES ----------


      subroutine SETUP_RAND_PERTURB( variable_in,& 
                       skebs_vertstruc,restart,                  & 
                       SP_AMP,SPFORCC,SPFORCS,ALPH,                  & 
                       VERTSTRUCC,VERTSTRUCS,VERTAMP,                &
                       KMINFORCT,KMAXFORCTH,LMINFORCT,LMAXFORCTH,    & 
                       KMAXFORCT,LMAXFORCT,                          &
                       itime_step,DX,DY,                             &
                       stepstoch,                                    & 
                       gridpt_stddev_rand_perturb, l_rand_perturb,   & 
                       tau_rand_perturb,                             &
                       TOT_BACKSCAT,ZTAU,REXPONENT,                  & 
                       kmax_ideal, lmax_ideal,                       & 
                       ids, ide, jds, jde, kds, kde,                 &
                       ims, ime, jms, jme, kms, kme,                 &
                       its, ite, jts, jte, kts, kte                  )



      IMPLICIT NONE

! General control 
      LOGICAL                                   :: restart
      REAL, PARAMETER                           :: RPI= 3.141592653589793 !4.0*atan(1.0) 
      CHARACTER, INTENT(IN)                     :: variable_in ! W=SKEBS_PSI, T=SKEBS_T, P=SPPT, R=RAND_PERTURB
      CHARACTER                                 :: variable

! Common to all schemes
      INTEGER , INTENT(IN)                      :: kmax_ideal, lmax_ideal
      INTEGER , INTENT(IN)                      :: ids, ide, jds, jde, kds, kde,   &
                                                   ims, ime, jms, jme, kms, kme,   &
                                                   its, ite, jts, jte, kts, kte
      INTEGER                                   :: IER,IK,IL,I,J,itime_step,stepstoch,skebs_vertstruc, & 
                                                   KMINFORCT,LMINFORCT,KMAXFORCT,LMAXFORCT,KMAXFORCTH,LMAXFORCTH, & 
                                                   KMAX,LMAX,LENSAV,ILEV
      REAL                                      :: DX,DY,RY,RX,ALPH,RHOKLMAX,ZREF,RHOKL,EPS
      REAL, DIMENSION (lmax_ideal,kmax_ideal)    :: SPFORCS,SPFORCC,SP_AMP
      REAL, DIMENSION (ims:ime,kms:kme,jms:jme) :: VERTSTRUCC,VERTSTRUCS
      REAL, DIMENSION (kms:kme)                 :: VERTAMP
      REAL, DIMENSION (ids:ide,jds:jde)         :: ZCHI

! SPPT and perturb_rand specific 
      REAL                                      :: gridpt_stddev_rand_perturb,kappat,tau_rand_perturb,l_rand_perturb
      REAL, DIMENSION (ims:ime,jms:jme)         :: var_sigma1


! SKEBS specific
      REAL                                      :: z,phi,ZGAMMAN,ZCONSTF0,TOT_BACKSCAT,ZTAU,REXPONENT,ZSIGMA2
      LOGICAL                                   :: is_print = .true.


      variable = variable_in
!     --------- SETUP PARAMETERS ---------------------------------------
      KMAX=(jde-jds)+1 !NLAT
      LMAX=(ide-ids)+1 !NLON
      RY=  KMAX*DY
      RX=  LMAX*DX
      LENSAV= 4*(KMAX+LMAX)+INT(LOG(REAL(KMAX))) + INT(LOG(REAL(LMAX))) + 8
      LENSAV=MAX(LENSAV,4*(KMAX_ideal+LMAX_ideal)+INT(LOG(REAL(KMAX_ideal))) + INT(LOG(REAL(LMAX_ideal))) + 8)

!     --------- ALLOCATE FIELDS FOR FFTPACK----------------------------
      IF ( ALLOCATED(WSAVE1)      ) DEALLOCATE(WSAVE1)
      IF ( ALLOCATED(WSAVE2)      ) DEALLOCATE(WSAVE2)
      IF ( ALLOCATED(WSAVE1_ideal)      ) DEALLOCATE(WSAVE1_ideal)
      IF ( ALLOCATED(WSAVE2_ideal)      ) DEALLOCATE(WSAVE2_ideal)
      IF ( ALLOCATED(rindarrayil)      ) DEALLOCATE(rindarrayil)
      IF ( ALLOCATED(rindarrayik)      ) DEALLOCATE(rindarrayik)
      IF ( ALLOCATED(rindarrayl)      ) DEALLOCATE(rindarrayl)
      IF ( ALLOCATED(rindarrayk)      ) DEALLOCATE(rindarrayk)
      IF ( ALLOCATED(WAVENUMBER_K)) DEALLOCATE(WAVENUMBER_K)
      IF ( ALLOCATED(WAVENUMBER_L)) DEALLOCATE(WAVENUMBER_L)
      ALLOCATE(WSAVE1(LENSAV),WSAVE2(LENSAV))
      ALLOCATE(WSAVE1_ideal(LENSAV),WSAVE2_ideal(LENSAV))
      ALLOCATE(rindarrayil(LMAX_ideal),rindarrayik(KMAX_ideal),rindarrayl(LMAX),rindarrayk(KMAX))
      ALLOCATE (wavenumber_k(jds:jde),wavenumber_l(ids:ide))

!     -------- INITIALIZE arrays to zero -----------------------------
      SP_AMP=0.0

!     -------- INITIALIZE FFTPACK ROUTINES -----------------------------
      call CFFT1I (LMAX, WSAVE1, LENSAV, IER)
      if(ier.ne. 0) write(*,95) ier

      call CFFT1I (KMAX, WSAVE2, LENSAV, IER)
      if(ier.ne. 0) write(*,95) ier

      call CFFT1I (LMAX_ideal, WSAVE1_ideal, LENSAV, IER)
      if(ier.ne. 0) write(*,95) ier

      call CFFT1I (KMAX_ideal, WSAVE2_ideal, LENSAV, IER)
      if(ier.ne. 0) write(*,95) ier

! Index arrays for spline interpolations - should be moved to setup routine
      do i=1,LMAX_ideal
        rindarrayil(i)=i
      enddo
      do i=1,KMAX_ideal
        rindarrayik(i)=i
      enddo
      do i=1,LMAX
        rindarrayl(i) = 1. + real(LMAX_ideal-1)/(LMAX-1)*(i-1)
      enddo
      do i=1,KMAX
        rindarrayk(i) = 1. + real(KMAX_ideal-1)/(KMAX-1)*(i-1)
      enddo

      95 format('error in cFFT2I=  ',i5)

     call findindex( wavenumber_k, wavenumber_l,             &
                      KMAX_ideal, LMAX_ideal)

!    set maximal perturbed wavenumber based on gridpoints in domain
     KMAXFORCT=min0(((ide-ids)+1)/2,((jde-jds)+1 )/2)-5
     LMAXFORCT=KMAXFORCT
     if (KMAXFORCT > KMAXFORCTH) then
        KMAXFORCT=KMAXFORCTH
     endif
     if (LMAXFORCT > LMAXFORCTH) then
        LMAXFORCT=LMAXFORCTH
     endif


!     --------------------------------------------------------------------------------------
!     ---------- INITIALIZE STOCHASTIC KINETIC-ENERGY BACKSCATTER SCHEME  (SKEBS) ----------
!     --------------------------------------------------------------------------------------
      ALPH   =  float(itime_step)/ZTAU/float(stepstoch) ! approximation of 1.-exp(-itime_step/ZTAU_PSI)
      ZSIGMA2=1./(12.0*ALPH)

      if (is_print) then
      IF (variable == 'W') then
      WRITE(*,'(''                                               '')')
      WRITE(*,'('' =============================================='')')
      WRITE(*,'('' >> Initializing STREAMFUNCTION forcing pattern of  << '')')
      WRITE(*,'('' >> stochastic kinetic-energy backscatter scheme << '')')
      WRITE(*,'('' Total backscattered energy, TOT_BACKSCAT_PSI '',E12.5)') TOT_BACKSCAT
      WRITE(*,'('' Exponent for energy spectra, REXPONENT_PSI ='',E12.5)') REXPONENT
      WRITE(*,'('' Minimal wavenumber of streamfunction forcing, LMINFORC ='',I10)') LMINFORCT
      WRITE(*,'('' Maximal wavenumber of streamfunction forcing, LMAXFORC ='',I10)') LMAXFORCT
      WRITE(*,'('' Minimal wavenumber of streamfunction forcing, KMINFORC ='',I10)') KMINFORCT
      WRITE(*,'('' Maximal wavenumber of streamfunction forcing, KMAXFORC ='',I10)') KMAXFORCT
      WRITE(*,'('' skebs_vertstruc                             '',I10)') skebs_vertstruc
      WRITE(*,'('' Time step: itime_step='',I10)') itime_step
      WRITE(*,'('' Decorrelation time of noise, ZTAU_PSI ='',E12.5)') ZTAU
      WRITE(*,'('' Variance of noise, ZSIGMA2_EPS  ='',E12.5)') ZSIGMA2
      WRITE(*,'('' Autoregressive parameter 1-ALPH_PSI ='',E12.5)') 1.-ALPH
      WRITE(*,'('' =============================================='')')

!     Unit of SPSTREAM_AMP: sqrt(m^2/s^3 1/s m**2(p+1)) m**-2(p/2) = m^/s^2 * m**[(p+1)-p] = m^2/s^2 m
      ELSEIF (variable == 'T') then
      WRITE(*,'(''                                               '')')
      WRITE(*,'('' =============================================='')')
      WRITE(*,'('' >> Initializing TEMPERATURE forcing pattern of  << '')')
      WRITE(*,'('' >> stochastic kinetic-energy backscatter scheme << '')')
      WRITE(*,'('' Total backscattered energy, TOT_BACKSCAT_T   '',E12.5)') TOT_BACKSCAT
      WRITE(*,'('' Exponent for energy spectra, REXPONENT_T   ='',E12.5)') REXPONENT
      WRITE(*,'('' Minimal wavenumber of tempearature forcing, LMINFORC ='',I10)') LMINFORCT
      WRITE(*,'('' Maximal wavenumber of tempearature forcing, LMAXFORC ='',I10)') LMAXFORCT
      WRITE(*,'('' Minimal wavenumber of tempearature forcing, KMINFORC ='',I10)') KMINFORCT
      WRITE(*,'('' Maximal wavenumber of tempearature forcing, KMAXFORC ='',I10)') KMAXFORCT
      WRITE(*,'('' skebs_vertstruc                             '',I10)') skebs_vertstruc
      WRITE(*,'('' Decorrelation time of noise, ZTAU_T ='',E12.5)') ZTAU
      WRITE(*,'('' Variance of noise, ZSIGMA2_ETA  ='',E12.5)') ZSIGMA2
      WRITE(*,'('' Autoregressive parameter 1-ALPH_T ='',E12.5)') 1.-ALPH
      WRITE(*,'('' =============================================='')')
      endif

     IF ((variable == 'P') .or. (variable == 'R')) then
      kappat= L_rand_perturb**2 ! L^2= kappa*T,  where L is a length scale in m; set to for L=100km 
      phi = exp (-float(itime_step*stepstoch)/tau_rand_perturb)
      alph = 1.-phi
     endif 

!     --------------------------------------------------------------------------------------
!     ---------- INITIALIZE STOCHASTICALLY PERTURBED PHYSICAL TENDENCY SCHEME --------------
!     --------------------------------------------------------------------------------------
     if (variable == 'P') then
      WRITE(*,'(''                                               '')')
      WRITE(*,'('' =============================================='')')
      WRITE(*,'('' >> Initializing Stochastically Perturbed Physics Tendency scheme << '')')
      WRITE(*,'('' sppt_vertstruc                             '',I10)') skebs_vertstruc
      WRITE(*,'('' Decorrelation time of noise, Tau ='',E12.5)') tau_rand_perturb
      WRITE(*,'('' Autoregressive parameter Phi ='',E12.5)') phi
      WRITE(*,'('' Length Scale L'',E12.5)') l_rand_perturb
      WRITE(*,'('' Variance in gridpoint space'',E12.5)') gridpt_stddev_rand_perturb
      WRITE(*,'('' =============================================='')')
      endif ! variable 

!     --------------------------------------------------------------------------------------
!     --------------------  INITIALIZE  RANDOM PERTUBATIONS  -------------------------------
!     --------------------------------------------------------------------------------------
     if (variable == 'R') then
      WRITE(*,'(''                                               '')')
      WRITE(*,'('' =============================================='')')
      WRITE(*,'('' >> Initializing random perturbations << '')')
      WRITE(*,'('' rand_pert_vertstruc                             '',I10)') skebs_vertstruc
      WRITE(*,'('' Decorrelation time of noise, Tau ='',E12.5)') tau_rand_perturb
      WRITE(*,'('' Autoregressive parameter Phi ='',E12.5)') phi
      WRITE(*,'('' Length Scale L'',E12.5)') l_rand_perturb
      WRITE(*,'('' Variance in gridpoint space'',E12.5)') gridpt_stddev_rand_perturb
      WRITE(*,'('' =============================================='')')
     endif ! variable 

     endif !is print
!     --------------------------------------------------------------------------------------
!     Compute Normalization constants       
!     --------------------------------------------------------------------------------------

     ZCHI    =  0.0
     ZGAMMAN  =  0.0
      ! Fill lower left quadrant of ZCHI. For this range the indeces IK,IL
      DO IK=jds-1,jde   ! These are now wavenumbers
      DO IL=ids-1,ide
      if (((sqrt((IK/RY*IK/RY)+(IL/RX*IL/RX)).lt.((KMAXFORCT+0.5)/RX)).and.&
           (sqrt((IK/RY*IK/RY)+(IL/RX*IL/RX)).ge.((KMINFORCT-0.5)/RX))) .or. &
          ((sqrt((IK/RY*IK/RY)+(IL/RX*IL/RX)).lt.((LMAXFORCT+0.5)/RY)).and.&
           (sqrt((IK/RY*IK/RY)+(IL/RX*IL/RX)).ge.((LMINFORCT-0.5)/RY))))then
        if ((IK>0).or.(IL>0)) then
          if (variable == 'W') then
            ZCHI(IL+1,IK+1)=((IK/RY*IK/RY)+(IL/RX*IL/RX))**(REXPONENT/2.)  ! SKEBS :U 
            ZGAMMAN= ZGAMMAN + ((IK/RY*IK/RY)+(IL/RX*IL/RX))**(REXPONENT+1)        
          else if (variable == 'T') then
            ZCHI(IL+1,IK+1)=((IK/RY*IK/RY)+(IL/RX*IL/RX))**(REXPONENT/2.)  ! SKEBS :T
            ZGAMMAN= ZGAMMAN + ((IK/RY*IK/RY)+(IL/RX*IL/RX))**(REXPONENT)          
          else if ((variable == 'P') .or. (variable == 'R')) then
            ZCHI(IL+1,IK+1)=exp(  -2*RPI**2*kappat*((IK/RY*IK/RY)+(IL/RX*IL/RX)) ) !SPPT
            ZGAMMAN= ZGAMMAN + exp(  -4*RPI**2*kappat*((IK/RY*IK/RY)+(IL/RX*IL/RX)) ) !SPPT
          endif
        endif
      endif
      enddo
      enddo
      ZGAMMAN=4.0*ZGAMMAN !account for all quadrants, although only one is Filled
      if (variable == 'W') then
         ZCONSTF0=SQRT(ALPH*TOT_BACKSCAT/(float(itime_step)*ZSIGMA2*ZGAMMAN))/(2*RPI)
      elseif  (variable == 'T') then     
         ! Mod by Nick to avoid precision problems
         !ZCONSTF0=SQRT(T0*ALPH*TOT_BACKSCAT/(float(itime_step)*cp*ZSIGMA2*ZGAMMAN))
         ZCONSTF0=SQRT(T0*ALPH*TOT_BACKSCAT)/SQRT(float(itime_step)*cp*ZSIGMA2)
         ZCONSTF0=ZCONSTF0/SQRT(ZGAMMAN)
      elseif ((variable == 'P') .or. (variable == 'R')) then
         ZCONSTF0= gridpt_stddev_rand_perturb*sqrt((1.-phi**2)/(2.*ZGAMMAN))
      endif

!     --------------------------------------------------------------------------------------
!     Now the wavenumber-dependent amplitudes
!     --------------------------------------------------------------------------------------
!     Note: There are symmetries and anti-symmetries to ensure real-valued back transforms
!     Fill lower left quadrant of matrix of noise amplitudes for wavenumbers K=0,KMAX/2
!
      SP_AMP=0.0
      DO IL=1, (LMAX_ideal/2+1)! upper left
        DO IK = 1, (KMAX_ideal/2+1)
          SP_AMP(IL,IK)=ZCONSTF0*ZCHI(IL,IK)
        ENDDO
      ENDDO
      DO IL=1, (LMAX_ideal/2+1) 
        DO IK = 1, (KMAX_ideal/2+1)
          SP_AMP(LMAX_ideal-IL+1,KMAX_ideal-IK+1)= SP_AMP(IL,IK) !lower right
          SP_AMP(IL,KMAX_ideal-IK+1)= SP_AMP(IL,IK) !lower right
          SP_AMP(LMAX_ideal-IL+1,IK)= SP_AMP(IL,IK) !lower right
        ENDDO
      ENDDO

!     -----------------------------------------


!     -----------------------------------------
!     Array for vertical structure if desired
      VERTAMP=1.0 ! Define vertical amplitude here.

      IF (skebs_vertstruc==1) then
        VERTSTRUCC=0.0
        VERTSTRUCS=0.0
        RHOKLMAX= sqrt(KMAX**2/DY**2 + LMAX**2/DX**2)
        ZREF=32.0
        DO ILEV=kts,kte
          DO IK=jts,jte
            DO IL=its,ite
            if (IL.le.(LMAX/2)) then
              RHOKL   = sqrt((IK+1)**2/DY**2 + (IL+1)**2/DX**2)
              EPS     = ((RHOKLMAX - RHOKL)/ RHOKLMAX)  * (ILEV/ZREF) * RPI
              VERTSTRUCC(IL,ILEV,IK) =  cos ( eps* (IL+1) )
              VERTSTRUCS(IL,ILEV,IK) =  sin ( eps* (IL+1) )
             else
              RHOKL   = sqrt((IK+1)**2/DY**2 + (LMAX-IL+2)**2/DX**2)
              EPS     = ((RHOKLMAX - RHOKL)/ RHOKLMAX)  * (ILEV/ZREF) * RPI
              VERTSTRUCC (IL,ILEV,IK) =   cos ( eps* (LMAX-IL+2) )
              VERTSTRUCS (IL,ILEV,IK) = - sin ( eps* (LMAX-IL+2) )
            endif
            ENDDO
          ENDDO
        ENDDO
      ENDIF

     END subroutine SETUP_RAND_PERTURB

!     ------------------------------------------------------------------
!************** UPDATE STOCHASTIC PATTERN IN WAVENUMBER SPACE**********
!     ------------------------------------------------------------------

     subroutine UPDATE_STOCH(                                         & 
                      SPFORCS,SPFORCC,SP_AMP,ALPH,                    &
                      restart,iseedarr,seed_dim,                      &
                      kmax_ideal, lmax_ideal,                         &
                      ids, ide, jds, jde, kds, kde,                   &
                      ims, ime, jms, jme, kms, kme,                   &
                      its, ite, jts, jte, kts, kte                    )
     IMPLICIT NONE

     REAL, DIMENSION(lmax_ideal,kmax_ideal)   ,INTENT(INOUT) :: SPFORCS,SPFORCC,SP_AMP
     INTEGER , INTENT(IN)     ::               ids, ide, jds, jde, kds, kde,   &
                                               ims, ime, jms, jme, kms, kme,   &
                                               its, ite, jts, jte, kts, kte
     INTEGER , INTENT(IN)     ::               seed_dim
   
     INTEGER, DIMENSION (seed_dim),             INTENT(INOUT) :: iseedarr
     INTEGER , INTENT(IN)                      :: kmax_ideal, lmax_ideal
     REAL, DIMENSION(LMAX_ideal,KMAX_ideal) :: ZRANDNOSS,ZRANDNOSC 
     INTEGER , ALLOCATABLE , DIMENSION(:) :: iseed
     REAL :: Z,ALPH
     REAL, PARAMETER :: thresh = 3.0
     INTEGER ::IL, IK,LMAX,KMAX,J
     INTEGER :: how_many
     LOGICAL :: LGAUSS,RESTART

     KMAX=(jde-jds)+1 !NLAT
     LMAX=(ide-ids)+1 !NATX

     call random_seed(put=iseedarr)

!    Pick the distribution of the noise

     LGAUSS=.true.
     IF (LGAUSS) then
       DO IK=1,KMAX_ideal
         DO IL=1,LMAX_ideal 
          do
           call gauss_noise(z)
           if (abs(z)<thresh) exit
          ENDDO
          ZRANDNOSS(IL,IK)=z
          do
           call gauss_noise(z)
           if (abs(z)<thresh) exit
          ENDDO
          ZRANDNOSC(IL,IK)=z
         ENDDO
       ENDDO
     ELSE
       DO IK=1,KMAX_ideal
         DO IL=1,LMAX_ideal 
           CALL RANDOM_NUMBER(z)
           ZRANDNOSS(IL,IK)=z-0.5
           CALL RANDOM_NUMBER(z)
           ZRANDNOSC(IL,IK)=z-0.5
          ENDDO
        ENDDO
      ENDIF

!     Note: There are symmetries and anti-symmetries to ensure real-valued back transforms
! for symmetric part: left and right half axis symmetric

        DO IK=1,KMAX_ideal
        if ((IK.le.(KMAX_ideal/2+1)) .and. (IK>1)) then ! Upper half
          DO IL=1,LMAX_ideal
            SPFORCC(IL,IK)       = (1.-ALPH)*SPFORCC(IL,IK)      + SP_AMP(IL,IK)     * ZRANDNOSC(IL,IK)  
            SPFORCS(IL,IK)       = (1.-ALPH)*SPFORCS(IL,IK)      + SP_AMP(IL,IK)     * ZRANDNOSS(IL,IK)  
          ENDDO
        ELSEIF (IK==1) then
          DO IL=1,LMAX_ideal
          if ((IL.le.(LMAX_ideal/2+1))) then
            SPFORCC(IL,IK)       = (1.-ALPH)*SPFORCC(IL,IK)      + SP_AMP(IL,IK)     * ZRANDNOSC(IL,IK)  
            SPFORCS(IL,IK)       = (1.-ALPH)*SPFORCS(IL,IK)      + SP_AMP(IL,IK)     * ZRANDNOSS(IL,IK)  
          elseif ((IL.gt.(LMAX_ideal/2+1))) then
            SPFORCC(IL,IK)       = (1.-ALPH)*SPFORCC(IL,IK)      + SP_AMP(IL,IK)     * ZRANDNOSC(LMAX_ideal-IL+2,IK)  
            SPFORCS(IL,IK)       = (1.-ALPH)*SPFORCS(IL,IK)      - SP_AMP(IL,IK)     * ZRANDNOSS(LMAX_ideal-IL+2,IK)  
          endif
          ENDDO
        ENDIF
        ENDDO
        DO IK=1,KMAX_ideal
        if (IK.gt.(KMAX_ideal/2+1)) then ! Lower half
          DO IL=1,LMAX_ideal
           if (IL.le.(LMAX_ideal/2+1).and.(IL.gt.1)) then !lower left 
             SPFORCC(IL,IK)      = (1.-ALPH)* SPFORCC(IL,IK)      + SP_AMP(IL,IK)      * ZRANDNOSC(LMAX_ideal-IL+2,KMAX_ideal-IK+2)
             SPFORCS(IL,IK)      = (1.-ALPH)* SPFORCS(IL,IK)      - SP_AMP(IL,IK)      * ZRANDNOSS(LMAX_ideal-IL+2,KMAX_ideal-IK+2)
            elseif (IL.eq.1) then !don't exceed index
             SPFORCC(IL,IK)      = (1.-ALPH)* SPFORCC(IL,IK)      + SP_AMP(IL,IK)      * ZRANDNOSC(        1,KMAX_ideal-IK+2)
             SPFORCS(IL,IK)      = (1.-ALPH)* SPFORCS(IL,IK)      - SP_AMP(IL,IK)      * ZRANDNOSS(        1,KMAX_ideal-IK+2)
            elseif (IL.gt.(LMAX_ideal/2+1)) then !lower right
             SPFORCC(IL,IK)      = (1.-ALPH)* SPFORCC(IL,IK)      + SP_AMP(IL,IK)      * ZRANDNOSC(LMAX_ideal-IL+2,KMAX_ideal-IK+2)
             SPFORCS(IL,IK)      = (1.-ALPH)* SPFORCS(IL,IK)      - SP_AMP(IL,IK)      * ZRANDNOSS(LMAX_ideal-IL+2,KMAX_ideal-IK+2)
            endif
          ENDDO
        endif
        ENDDO
     
      call random_seed(get=iseedarr)

     END subroutine UPDATE_STOCH
!     ------------------------------------------------------------------
      SUBROUTINE UPDATE_STOCH_TEN(ru_tendf,rv_tendf,t_tendf, &
                       ru_tendf_stoch,rv_tendf_stoch,rt_tendf_stoch,& 
                       mu,mub,c1h,c2h,                              & 
                       ids, ide, jds, jde, kds, kde,                &
                       ims, ime, jms, jme, kms, kme,                &
                       its, ite, jts, jte, kts, kte,                &
                       kte_stoch,kme_stoch                          )

       IMPLICIT NONE
       INTEGER , INTENT(IN)        ::  ids, ide, jds, jde, kds, kde,   &
                                       ims, ime, jms, jme, kms, kme,   &
                                       its, ite, jts, jte, kts, kte,   & 
                                       kte_stoch,kme_stoch

       REAL , DIMENSION(ims:ime , kms:kme, jms:jme),INTENT(INOUT) :: &
                                       ru_tendf, rv_tendf, t_tendf

       REAL , DIMENSION(ims:ime , kms:kme_stoch, jms:jme)           :: &
                      ru_tendf_stoch,rv_tendf_stoch,rt_tendf_stoch 

       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: mu,mub
       REAL , DIMENSION(kms:kme) , INTENT(IN) :: c1h,c2h

       INTEGER :: I,J,K,kh
       REAL  :: dt,xm


   
       DO j = jts,MIN(jde-1,jte)
         DO k = kts,kte-1
           kh=min(k,kte_stoch)
           DO i = its,ite 
             ru_tendf(i,k,j) = ru_tendf(i,k,j) +  ru_tendf_stoch(i,kh,j)  * (mu(i,j)+mub(i,j)) 
           ENDDO
         ENDDO
       ENDDO

       DO j = jts,jte
         DO k = kts,kte-1
           kh=min(k,kte_stoch)
           DO i = its,MIN(ide-1,ite)
             rv_tendf(i,k,j) = rv_tendf(i,k,j) +  rv_tendf_stoch(i,kh,j) *  (mu(i,j)+mub(i,j)) 
           ENDDO
         ENDDO
       ENDDO

       DO j = jts,MIN(jde-1,jte)
         DO k = kts,kte-1
           kh=min(k,kte_stoch)
           DO i = its,MIN(ide-1,ite)
             t_tendf(i,k,j) = t_tendf(i,k,j) + rt_tendf_stoch(i,kh,j) * (mu(i,j)+mub(i,j))
           ENDDO
         ENDDO
       ENDDO

       END SUBROUTINE UPDATE_STOCH_TEN
!     ------------------------------------------------------------------
!!************** PERTURB PHYSICS TENDENCIES (except T) FOR SPPT *******************
!     ------------------------------------------------------------------
      subroutine perturb_physics_tend(gridpt_stddev_sppt,               & 
                       sppt_thresh_fact,rstoch,                         & 
                       ru_tendf,rv_tendf,t_tendf,moist_tend,            &
                       ids, ide, jds, jde, kds, kde,                    &
                       ims, ime, jms, jme, kms, kme,                    &
                       its, ite, jts, jte, kts, kte,                    & 
                       kte_stoch,kme_stoch                               )

!      This subroutine add stochastic perturbations of the form 
!
!                  rx_tendf(i,k,j)      = rx_tendf(i,k,j)*(1.0 + rstoch(i,k,j))
!
!       to the tendencies of  U, V, and Q. 
!       Since the temperature perturbations do not include the micro-physics
!       tendencies at this point, the stochastic tendency perturbations to 
!       temperature are added in subroutine rk_addtend_dry of module module_em.F

       IMPLICIT NONE
       INTEGER , INTENT(IN)        ::  ids, ide, jds, jde, kds, kde,   &
                                       ims, ime, jms, jme, kms, kme,   &
                                       its, ite, jts, jte, kts, kte,   & 
                                       kte_stoch,kme_stoch

       REAL , DIMENSION(ims:ime , kms:kme, jms:jme),INTENT(INOUT) ::   &
                                        ru_tendf, rv_tendf, t_tendf,moist_tend
       REAL , DIMENSION(ims:ime,kms:kme_stoch, jms:jme),INTENT(INOUT) :: rstoch
       REAL :: gridpt_stddev_sppt ,thresh,sppt_thresh_fact

       INTEGER :: I,J,K,kh

! Here the random process at each gridpoint is capped if it exceeds a value thresh 

       thresh=sppt_thresh_fact*gridpt_stddev_sppt
       DO j = jts,jte
         DO k = kts,min(kte-1,kte_stoch-1)
           DO i = its,ite 
!                rstoch(i,k,j)=MAX(MIN(rstoch(i,k,j),thresh),-1.*thresh))
             if (rstoch(i,k,j).lt.-thresh) then
                 rstoch(i,k,j)=-thresh
             endif
             if (rstoch(i,k,j).gt.thresh) then
                 rstoch(i,k,j)=thresh
             endif
           ENDDO
         ENDDO
       ENDDO

! Perturb the tendencies of u,v,q,t.
       DO j = jts,MIN(jde-1,jte)
         DO k = kts,kte-1
         kh = min( k, kte_stoch-1 ) 
           DO i = its,ite 
              ru_tendf(i,k,j)      = ru_tendf(i,k,j)*(1.0 + rstoch(i,kh,j))
           ENDDO
         ENDDO
       ENDDO

       DO j = jts,jte
         DO k = kts,kte-1
         kh = min( k, kte_stoch-1 ) 
            DO i = its,MIN(ide-1,ite)
              rv_tendf(i,k,j)      = rv_tendf(i,k,j)*(1.0 + rstoch(i,kh,j))
           ENDDO
         ENDDO
       ENDDO

       DO j = jts,MIN(jde-1,jte)
         DO k = kts,kte-1
         kh = min( k, kte_stoch-1 ) 
           DO i = its,MIN(ide-1,ite)
              moist_tend(i,k,j)    = moist_tend(i,k,j)*(1.0 + rstoch(i,kh,j))
              t_tendf   (i,k,j)    =    t_tendf(i,k,j)*(1.0 + rstoch(i,kh,j))
           ENDDO
         ENDDO
       ENDDO

      end subroutine perturb_physics_tend

!     ------------------------------------------------------------------
!!************** UPDATE SPECTRAL PATTERN AND TRANSFORM TO GRIDPOINT SPACE***
!     ------------------------------------------------------------------
! This subroutine evolves the spectral pattern and transforms it back to gridpoint space.


      SUBROUTINE RAND_PERT_UPDATE       (grid, variable_in,                   &
                          SPFORCS,SPFORCC,SP_AMP,ALPH_RAND,                   &
                          ips, ipe, jps, jpe, kps, kpe,                       &
                          ids, ide, jds, jde, kds, kde,                       &
                          ims, ime, jms, jme, kms, kme,                       &
                          kts, kte,                                           &
                          imsx,imex,jmsx,jmex,kmsx,kmex,                      &
                          ipsx,ipex,jpsx,jpex,kpsx,kpex,                      &
                          imsy,imey,jmsy,jmey,kmsy,kmey,                      &
                          ipsy,ipey,jpsy,jpey,kpsy,kpey,                      &
                          kmax_ideal, lmax_ideal,                             &
                          kpe_stoch,kde_stoch,kme_stoch,kte_stoch,            &
                          restart,iseedarr,seed_dim,                          &
                          DX,DY,skebs_vertstruc,                          &
                          RAND_PERT,thresh_fact,gridpt_stddev,                &
                          VERTSTRUCC,VERTSTRUCS,VERTAMP                       )



    USE module_domain, ONLY : domain
#ifdef DM_PARALLEL
    USE module_dm, ONLY : local_communicator, mytask, ntasks, ntasks_x, ntasks_y, local_communicator_periodic, &
                          wrf_dm_maxval, wrf_err_message, local_communicator_x, local_communicator_y, data_order_xzy
#endif


      IMPLICIT NONE

      TYPE ( domain ), INTENT(INOUT) :: grid

 
      INTEGER , INTENT(IN)     ::               kmax_ideal, lmax_ideal
      INTEGER , INTENT(IN)     ::               ids, ide, jds, jde, kds, kde,   &
                                                ims, ime, jms, jme, kms, kme,   &
                                                ips, ipe, jps, jpe, kps, kpe,   & 
                                                kts, kte                       
      INTEGER , INTENT(IN)     ::               imsx,imex,jmsx,jmex,kmsx,kmex, &
                                                ipsx,ipex,jpsx,jpex,kpsx,kpex, &
                                                imsy,imey,jmsy,jmey,kmsy,kmey, &
                                                ipsy,ipey,jpsy,jpey,kpsy,kpey
      INTEGER , INTENT(IN)     ::               seed_dim
      INTEGER                  ::               kpe_stoch,kde_stoch,kme_stoch,kte_stoch 

      REAL    , INTENT(IN)     ::               ALPH_RAND,dx,dy,thresh_fact,gridpt_stddev
      INTEGER , INTENT(IN)     ::               skebs_vertstruc
      CHARACTER, INTENT(IN)    ::               variable_in ! T, U, V                
                                               ! T          ! random field, T
                                               ! U          ! first derivative of streamfunction with regard to y; for skebs: U
                                               ! V          ! first derivative of streamfunction with regard to x; for skebs: V

      INTEGER, DIMENSION (seed_dim),             INTENT(INOUT) :: iseedarr
      REAL, DIMENSION(ims:ime,kms:kme, jms:jme),INTENT(IN)    :: VERTSTRUCC,VERTSTRUCS
      REAL, DIMENSION(lmax_ideal,kmax_ideal)   ,INTENT(INOUT) :: SPFORCS,SPFORCC,SP_AMP
      REAL, DIMENSION(kms:kme )                ,INTENT(IN)    :: VERTAMP
      REAL, DIMENSION(ims:ime,kms:kme_stoch, jms:jme)         :: RAND_PERT  
      REAL                     ::               RY,RX


! Local Variabels 

      INTEGER                  ::               IK,IL,ILEV,NLON,NLAT,IJ,I,J,K,numbands
      COMPLEX, DIMENSION (LMAX_ideal)            :: dummy_complex2   
      COMPLEX, DIMENSION (KMAX_ideal)            :: dummy_complex3   
      COMPLEX, DIMENSION (ipsx:ipex)            :: dummy_complex
      REAL                     ::               thresh
      REAL, DIMENSION(LMAX_ideal,KMAX_ideal)    :: ZFORCC,ZFORCS
      REAL, DIMENSION(ids:ide,jds:jde)          :: ZFORCg
      INTEGER                  ::               IER,LENWRK,KMAX,LMAX,LENSAV,JJ,II
      INTEGER                  ::                its,ite,jts,jte,ind
      REAL, ALLOCATABLE        :: WORK(:)
      CHARACTER (LEN=160)      :: mess
      real, dimension (LMAX_ideal) ::  b,c,d
                                              ! Assums LMAX_ideal=KMAX_ideal which is the case for homogenous fields


      LOGICAL :: RESTART
      CHARACTER  :: variable

      variable = variable_in

      NLAT=(jde-jds)+1 !KMAX
      NLON=(ide-ids)+1 !LMAX
      KMAX = NLAT
      LMAX = NLON
      RY=  NLAT*DY
      RX=  NLON*DX

      LENWRK=2*KMAX*LMAX
      ALLOCATE(WORK(LENWRK))
      LENSAV= 4*(KMAX+LMAX)+INT(LOG(REAL(KMAX))) + INT(LOG(REAL(LMAX))) + 8


! Update the pattern generator by evolving each spectral coefficients as AR1

              IF (variable .ne. 'V') THEN  !T, random field, U, don't update for V 
               CALL UPDATE_STOCH( &
                          SPFORCS,SPFORCC,SP_AMP,ALPH_RAND,                   &
                          restart,iseedarr,seed_dim,                          &
                          kmax_ideal, lmax_ideal,                             &
                          ids, ide, jds, jde, kds, kde,                       &
                          ims, ime, jms, jme, kms, kme,                       &
			  kts,kte,kts,kte,kts, kte                        )      
              endif
   
         
! Take derivatives of streamfunction of for skebs
             ZFORCC=SPFORCC
             ZFORCS=SPFORCS

              IF (variable == 'U') THEN !U
               DO J = 1, KMAX_ideal
                 DO I = 1,  LMAX_ideal
                       ZFORCC(I,J)  =  2*RPI/RY*  wavenumber_k(J) * ZFORCS(I,J)
                       ZFORCS(I,J)  = -2*RPI/RY*  wavenumber_k(J) * ZFORCC(I,J)
                 ENDDO
               ENDDO
              ELSEIF (variable == 'V') THEN !V
               DO J = 1, KMAX_ideal
                 DO I = 1,  LMAX_ideal
                       ZFORCC(I,J)  = -2*RPI/RX*  wavenumber_l(I) * ZFORCS(I,J)
                       ZFORCS(I,J)  =  2*RPI/RX*  wavenumber_l(I) * ZFORCC(I,J)
                 ENDDO
               ENDDO
              endif

         DO J = 1, KMAX_ideal
           DO i = 1,  LMAX_ideal
             dummy_complex2(i)=cmplx(ZFORCC(i,j),ZFORCS(i,j))
           ENDDO
           CALL cFFT1B (LMAX_ideal, 1 ,dummy_complex2,LMAX_ideal, WSAVE1_ideal, LENSAV, WORK, LENWRK, IER) 
           if (ier.ne.0) then 
              WRITE(mess,FMT='(A)') 'error in cFFT1B in do_fftback_along_x, field U'
              CALL wrf_debug(0,mess)
           end if
           DO i = 1,  LMAX_ideal 
             ZFORCC(i,j)=real(dummy_complex2(i))
             ZFORCS(i,j)=imag(dummy_complex2(i))
           END DO
         END DO ! J

          DO I = 1, LMAX_ideal 
           DO j = 1, KMAX_ideal 
            dummy_complex3(j)=cmplx(ZFORCC(i,j),ZFORCS(i,j))
           ENDDO
            CALL cFFT1B (KMAX_ideal, 1 ,dummy_complex3,KMAX_ideal, WSAVE2_ideal, LENSAV, WORK, LENWRK, IER) 
            if (ier.ne.0) then 
               WRITE(mess,FMT='(A)') 'error in cFFT1B in do_fftback_along_y, field U'
               CALL wrf_debug(0,mess)
            end if
            DO j =  1, KMAX_ideal 
                 ZFORCC(i,j)=real(dummy_complex3(j))
                 ZFORCS(i,j)=imag(dummy_complex3(j))
            END DO
          END DO ! I

       ! Spline interpolation
         if ((LMAX_ideal.ne.ide).or.(KMAX_ideal.ne.jde))  then
            ! https://ww2.odu.edu/~agodunov/computing/programs/book2/Ch01/spline.f90
             DO J = 1, KMAX_ideal 
              DO I = 1,  LMAX_ideal
               dummy_complex2(i)=cmplx(ZFORCC(i,j),0.0)
              ENDDO
              call spline (rindarrayil, real(dummy_complex2) , b, c, d,  LMAX_ideal )
              do i=1,LMAX
                ZFORCg(i,j)=ispline (rindarrayl(i),rindarrayil, real(dummy_complex2) , b, c, d,  LMAX_ideal )
              enddo
             ENDDO
             DO I = 1, LMAX
               DO J = 1, KMAX_ideal
                 dummy_complex3(j)=cmplx(ZFORCg(i,j),0.0)
               ENDDO
               call spline (rindarrayik, real(dummy_complex3) , b, c, d,  KMAX_ideal )
               do J=1,KMAX
                 ZFORCg(i,j)=ispline (rindarrayk(j),rindarrayik, real(dummy_complex3) , b, c, d,  KMAX_ideal )
               enddo
             ENDDO
         else
          DO J = 1, KMAX_ideal 
            DO I = 1,  LMAX_ideal
              ZFORCg(i,j)=ZFORCC(i,j)
            ENDDO
           ENDDO
         endif ! LMAX_ideal

      thresh=thresh_fact*gridpt_stddev
      RAND_PERT=0.0
      !$OMP PARALLEL DO   &
        !$OMP PRIVATE ( ij )
        DO ij = 1 , grid%num_tiles
               DO k=kts,min(kte,grid%num_stoch_levels)
                 DO I=grid%i_start(ij), grid%i_end(ij)
                   DO j=grid%j_start(ij), grid%j_end(ij)
                    RAND_PERT(I,K,J)=MAX(MIN(ZFORCg(I,J),thresh),-1.0*thresh)
                   ENDDO
                 ENDDO
               ENDDO
         ENDDO !OMP
        !$OMP END PARALLEL DO

       END SUBROUTINE RAND_PERT_UPDATE
!     ------------------------------------------------------------------
!!************** UPDATE SPECTRAL PATTERN AND TRANSFORM TO GRIDPOINT SPACE***
!     ------------------------------------------------------------------
! This subroutine evolves the spectral pattern and transforms it back to gridpoint space.
        subroutine read_write_stochrestart (filename1,grid,config_flags,lmax_ideal,kmax_ideal,SPFORS,SPFORC,nswitch)  
!lala
      USE module_configure , ONLY : grid_config_rec_type
      USE module_date_time
      USE module_timing
      USE module_utility
      USE module_domain
#ifdef DM_PARALLEL
      USE module_dm
#endif

      IMPLICIT NONE
      TYPE(domain)                               :: grid 
      TYPE (grid_config_rec_type) , INTENT(IN)   :: config_flags

      REAL, DIMENSION(lmax_ideal,kmax_ideal)   ,INTENT(INOUT) :: SPFORC, SPFORS
      INTEGER , INTENT(IN)                      :: kmax_ideal,lmax_ideal,nswitch
      INTEGER                                   :: ierr 

      INTEGER , PARAMETER                       :: write_restart = 1
      INTEGER , PARAMETER                       :: read_restart = 2
      INTEGER , PARAMETER                       :: fid=19
      TYPE(WRFU_Time)                           :: next_time, currentTime, startTime
      CHARACTER (LEN=256)                       :: current_date_char, date_string, next_timestr_char
      character(len=13)                         :: filename1
      character(len=50)                         :: filename
      character(len=70)                         :: message 
      CHARACTER (len=len_current_date) current_timestr,  next_timestr
      logical itsopen 


      LOGICAL, EXTERNAL      :: wrf_dm_on_monitor

     if(wrf_dm_on_monitor()) then
      IF ( nswitch .eq. write_restart) then ! WRITE!!!

         CALL domain_clock_get( grid, current_timestr=current_date_char )
         !write(filename,FMT='(A,A,A)'),trim(filename1),'_rst_d01_',TRIM(current_date_char) 
         write(filename,FMT='(A,A,A)'),trim(filename1),'_rst_d01'
         inquire(unit=19, opened=itsopen) 
         if ( itsopen ) then
            write(*,*) 'Its open already'
         else
          print*,'opening ', filename, ' for writing restart file for stochastic suite'
          OPEN(19, FILE=filename,FORM='UNFORMATTED',STATUS='REPLACE',IOSTAT=ierr)
          IF(ierr .NE. 0 ) THEN
            WRITE(message,FMT='(A)') &
            'problem with writing restart files for stochastic physics suite'
             CALL wrf_error_fatal ( message )
          END IF
         endif
         write(19),SPFORC,SPFORS
         close(19)

      else IF ( nswitch .eq. read_restart) then ! READ!!!
        CALL domain_clock_get( grid, current_timestr=current_date_char )
         !write(filename,FMT='(A,A,A)'),trim(filename1),'_rst_d01_',TRIM(current_date_char) 
         ! The clock has not yet been advanced, hence the time on the stoch restart file is
         ! istep too early. Want to add tstep to minute
         ! write(filename,FMT='(A,A)'),'stochrst_d01_','2000-01-24_12:12:00'
         write(filename,FMT='(A,A,A)'),trim(filename1),'_rst_d01'

         inquire(unit=19, opened=itsopen) 
         if ( itsopen ) then
           write(*,*) 'Its open already'
         else
           print*,'RESTART run: opening ', filename, ' for reading'
           OPEN(19, FILE=filename,FORM='UNFORMATTED',STATUS='OLD',IOSTAT=ierr)
         IF(ierr .NE. 0 ) THEN
           WRITE(message,FMT='(A)') &
           'problem with reading restart files for stochastic physics suite'
            CALL wrf_error_fatal ( message )
         END IF
         end if
         read(19),SPFORC,SPFORS
         close(19)
      endif
     endif  !wrf_dm_on_monitor

      END SUBROUTINE read_write_stochrestart
!     ------------------------------------------------------------------
!!************** TRANSFORM FROM GRIDPOINT SPACE TO SPHERICAL HARMONICS **
!     ------------------------------------------------------------------
            subroutine findindex( wavenumber_k, wavenumber_L,             &    
                       KMAX_ideal, LMAX_ideal) 

      IMPLICIT NONE 
      INTEGER :: IK,IL,KMAX,LMAX,KMAX_ideal, LMAX_ideal
      INTEGER, DIMENSION (1:KMAX_ideal)::  wavenumber_k
      INTEGER, DIMENSION (1:LMAX_ideal)::  wavenumber_l

      KMAX=KMAX_ideal  
      LMAX=LMAX_ideal 

      !map wave numbers K,L to indeces IK, IL
      DO IK=1,KMAX/2+1 
        wavenumber_k(IK)=IK-1 
      ENDDO
      DO IK=KMAX,KMAX/2+2,-1
        wavenumber_k(IK)=IK-KMAX-1 
      ENDDO
      DO IL=1,LMAX/2+1 
        wavenumber_l(IL)=IL-1
      ENDDO
      DO IL=LMAX,LMAX/2+2,-1
        wavenumber_l(IL)=IL-LMAX-1 
      ENDDO

      END subroutine findindex 
!     ------------------------------------------------------------------
     subroutine gauss_noise(z)
      real :: z                    ! output
      real :: x,y,r, coeff         ! INPUT

!  [2.1] Get two uniform variate random numbers IL range 0 to 1:

      do

      call random_number( x )
      call random_number( y )

!     [2.2] Transform to range -1 to 1 and calculate sum of squares:

      x = 2.0 * x - 1.0
      y = 2.0 * y - 1.0
      r = x * x + y * y

      if ( r > 0.0 .and. r < 1.0 ) exit

      end do
!
!  [2.3] Use Box-Muller transformation to get normal deviates:

      coeff = sqrt( -2.0 * log(r) / r )
      z = coeff * x

     end subroutine gauss_noise
!     ------------------------------------------------------------------
     SUBROUTINE rand_seed (config_flags, iseed1, iseedarr, seed_start, seed_dim )
     USE module_configure
     IMPLICIT NONE
!
!  Structure that contains run-time configuration (namelist) data for domain
      TYPE (grid_config_rec_type)                       :: config_flags
!
! Arguments
     INTEGER                                              :: iseed1 , seed_start, seed_dim
     INTEGER, DIMENSION (seed_start:seed_dim), INTENT(OUT):: iseedarr

! Local
      integer*8          :: fctime, one_big
      integer            :: i 

      fctime = config_flags%start_year * ( config_flags%start_month*100+config_flags%start_day) + config_flags%start_hour

      one_big = 1
      iseedarr=0
      do i = seed_start,seed_dim
         iseedarr(i)=mod(fctime+iseed1*config_flags%nens*1000000,19211*one_big)
      enddo

      end SUBROUTINE rand_seed
! ------------------------------------------------------------------
   subroutine spline (x, y, b, c, d, n)
!======================================================================
!  Calculate the coefficients b(i), c(i), and d(i), i=1,2,...,n
!  for cubic spline interpolation
!  s(x) = y(i) + b(i)*(x-x(i)) + c(i)*(x-x(i))**2 + d(i)*(x-x(i))**3
!  for  x(i) <= x <= x(i+1)
!  Alex G: January 2010
!----------------------------------------------------------------------
!  input..
!  x = the arrays of data abscissas (in strictly increasing order)
!  y = the arrays of data ordinates
!  n = size of the arrays xi() and yi() (n>=2)
!  output..
!  b, c, d  = arrays of spline coefficients
!  comments ...
!  spline.f90 program is based on fortran version of program spline.f
!  the accompanying function fspline can be used for interpolation
!======================================================================
implicit none
integer n
real x(n), y(n), b(n), c(n), d(n)
integer i, j, gap
real h

gap = n-1
! check input
if ( n < 2 ) return
if ( n < 3 ) then
  b(1) = (y(2)-y(1))/(x(2)-x(1))   ! linear interpolation
  c(1) = 0.
  d(1) = 0.
  b(2) = b(1)
  c(2) = 0.
  d(2) = 0.
  return
end if
!
! step 1: preparation
!
d(1) = x(2) - x(1)
c(2) = (y(2) - y(1))/d(1)
do i = 2, gap
  d(i) = x(i+1) - x(i)
  b(i) = 2.0*(d(i-1) + d(i))
  c(i+1) = (y(i+1) - y(i))/d(i)
  c(i) = c(i+1) - c(i)
end do
!
! step 2: end conditions 
!
b(1) = -d(1)
b(n) = -d(n-1)
c(1) = 0.0
c(n) = 0.0
if(n /= 3) then
  c(1) = c(3)/(x(4)-x(2)) - c(2)/(x(3)-x(1))
  c(n) = c(n-1)/(x(n)-x(n-2)) - c(n-2)/(x(n-1)-x(n-3))
  c(1) = c(1)*d(1)**2/(x(4)-x(1))
  c(n) = -c(n)*d(n-1)**2/(x(n)-x(n-3))
end if
!
! step 3: forward elimination 
!
do i = 2, n
  h = d(i-1)/b(i-1)
  b(i) = b(i) - h*d(i-1)
  c(i) = c(i) - h*c(i-1)
end do
!
! step 4: back substitution
!
c(n) = c(n)/b(n)
do j = 1, gap
  i = n-j
  c(i) = (c(i) - d(i)*c(i+1))/b(i)
end do
!
! step 5: compute spline coefficients
!
b(n) = (y(n) - y(gap))/d(gap) + d(gap)*(c(gap) + 2.0*c(n))
do i = 1, gap
  b(i) = (y(i+1) - y(i))/d(i) - d(i)*(c(i+1) + 2.0*c(i))
  d(i) = (c(i+1) - c(i))/d(i)
  c(i) = 3.*c(i)
end do
c(n) = 3.0*c(n)
d(n) = d(n-1)
end subroutine spline

  function ispline(u, x, y, b, c, d, n)
!======================================================================
! function ispline evaluates the cubic spline interpolation at point z
! ispline = y(i)+b(i)*(u-x(i))+c(i)*(u-x(i))**2+d(i)*(u-x(i))**3
! where  x(i) <= u <= x(i+1)
!----------------------------------------------------------------------
! input..
! u       = the abscissa at which the spline is to be evaluated
! x, y    = the arrays of given data points
! b, c, d = arrays of spline coefficients computed by spline
! n       = the number of data points
! output:
! ispline = interpolated value at point u
!=======================================================================
implicit none
real ispline
integer n
real  u, x(n), y(n), b(n), c(n), d(n)
integer i, j, k
real dx

! if u is ouside the x() interval take a boundary value (left or right)
if(u <= x(1)) then
  ispline = y(1)
  return
end if
if(u >= x(n)) then
  ispline = y(n)
  return
end if

!*
!  binary search for for i, such that x(i) <= u <= x(i+1)
!*
i = 1
j = n+1
do while (j > i+1)
  k = (i+j)/2
  if(u < x(k)) then
    j=k
    else
    i=k
   end if
end do
!*
!  evaluate spline interpolation
!*
dx = u - x(i)
ispline = y(i) + dx*(b(i) + dx*(c(i) + dx*d(i)))
end function ispline
!     ------------------------------------------------------------------
      end module module_stoch