dyn_em/module_stoch.F

   1
   2 module module_stoch
   3 !***********************************************************************
   4 !
   5 ! Purpose: Stochastic Perturbation Schemes
   6 ! Author : Judith Berner, NCAR  (berner@ucar.edu)
   7 ! Date   : Apr 2014
   8 !
   9 !***********************************************************************
  10 !
  11 ! The scheme introduces stochastic perturbations to the rotational wind
  12 ! components and to the potential temperature field. The stochastic
  13 ! perturbations are generated by independent autoregressive processes for
  14 ! each wavenumber and results in smooth spatially and temporally correlated patterns.
  15 ! Details of the scheme and its performance in a meso-scale WRF-ensemble
  16 ! system are available in:
  17 !
  18 ! Berner, J.,  S.-Y. Ha, J. P. Hacker, A. Fournier and C. Snyder 2011:
  19 ! "Model uncertainty in a mesoscale ensemble prediction system: Stochastic
  20 ! versus multi-physics representations",  2011, Mon. Wea. Rev., 139, 1972-1995
  21 ! http://journals.ametsoc.org/doi/abs/10.1175/2010MWR3595.1
  22 !
  23 ! Features:
  24 !     Dissipation: Dissipation rates are assumed constant in space and time
  25 !     Vertical structure: Supports two options for vertical structure:
  26 !                         0) constant
  27 !                         1) random phase
  28 !
  29 ! Optional namelist parameters:
  30 !   stoch_force_opt = 0, 0, 0: No stochastic parameterization
  31 !    = 1, 1, 1: Use SKEB scheme
  32 !   skebs_vertstruc = 0, 0, 0: Constant vertical structure of random pattern generator
  33 !    = 1, 1, 1: Random phase vertical structure random pattern generator
  34 !   tot_backscat_psi : Total backscattered dissipation rate for streamfunction; Controls
  35 !                                   amplitude of rotational wind perturbations Default value is 1.0E-5 m2/s3.
  36 !   tot_backscat_t : Total backscattered dissipation rate for potential temperature;
  37 !    Controls amplitude of potential temperature perturbations. Default value is 1.0E-6 m2/s3.
  38 !   nens : Random seed for random number stream. This parameter needs to be different
  39 !                                   for each member in ensemble forecasts. Is a function of initial start time
  40 !    to ensure different random number streams for different forecasts.
  41 !   ztau_psi : Decorrelation time (in s) for streamfunction perturbations.
  42 !    Default is 10800s.  Recommended value is 216000s.
  43 !   ztau_t : Decorrelation time (in s) for potential temperature perturbations.
  44 !    Default 10800s. Recommended value is 216000s.
  45 !   rexponent_psi : Spectral slope for streamfunction perturbations. Default is -1.83
  46 !                                   for a kinetic-energy forcing spectrum with slope -5/3.
  47 !   rexponent_t : Spectral slope of potential temperature perturbations. Default is -1.83
  48 !    for a potential energy forcing spectrum with slope -1.832.
  49 !   kminforc : Minimal forcing wavenumber in longitude for streamfunction perturbations. Default is 1.
  50 !   lminforc : Minimal forcing wavenumber in latitude for streamfunction perturbations. Default is 1.
  51 !   kminforc : Minimal forcing wavenumber in longitude for potential temperature perturbations. Default is 1.
  52 !   lminforct : Minimal forcing wavenumber in latitude for potential temperature perturbations. Default is 1.
  53 !   kmaxforc : Maximal forcing wavenumber in longitude for streamfunction perturbations.
  54 !    Default is maximal possible wavenumbers determined by number of gridpoints.
  55 !   lmaxforc : Maximal forcing wavenumber in latitude for streamfunction perturbations.
  56 !    Default is maximal possible wavenumbers determined by number of gridpoints.
  57 !   kmaxforct : Maximal forcing wavenumber in longitude for potential temperature perturbations.
  58 !    Default is maximal possible wavenumbers determined by number of gridpoints.
  59 !   lmaxforct : Maximal forcing wavenumber in latitude for potential temperature perturbations.
  60 !    Default is maximal possible wavenumbers determined by number of gridpoints.
  61 !   zsigma2_eps : Noise variance in autoregressive process defining streamfunction perturbations.
  62 !   zsigma2_eta : Noise variance in autoregressive process defining in potential temperature perturbations.
  63
  64 !***********************************************************************
  65 !     ------------------------------------------------------------------
  66 !************** DECLARE FIELDS AND VARIABLES FOR STOCHASTIC BACKSCATTER
  67 !     ------------------------------------------------------------------
  68       implicit none
  69       public ::  SETUP_RAND_PERTURB, UPDATE_STOCH,&
  70                          do_fftback_along_x,do_fftback_along_y,&
  71                          rand_pert_update
  72
  73       INTEGER :: LMINFORC, LMAXFORC, KMINFORC, KMAXFORC, &
  74       &          LMINFORCT, LMAXFORCT, KMINFORCT, KMAXFORCT
  75       REAL    :: ALPH, ALPH_PSI, ALPH_T, TOT_BACKSCAT_PSI, TOT_BACKSCAT_T,  REXPONENT_PSI,REXPONENT_T
  76
  77 !     ----------Fields for spectral transform -----------
  78
  79       INTEGER :: LENSAV
  80       INTEGER,ALLOCATABLE:: wavenumber_k(:), wavenumber_l(:)
  81       REAL, ALLOCATABLE :: WSAVE1(:),WSAVE2(:)
  82
  83 !     --------- Others -------------------------------------------------
  84       REAL, PARAMETER:: RPI= 3.141592653589793 !4.0*atan(1.0)
  85       REAL, PARAMETER:: CP= 1006.0 ! specific heat of dry air in J/(Kg*K)= m^2/(K* s^2)
  86       REAL, PARAMETER:: T0= 300.0 ! Reference temperature in K
  87
  88       save
  89
  90
  91 !=======================================================================
  92 contains
  93 !=======================================================================
  94 !     ------------------------------------------------------------------
  95 !!************** INITIALIZE STOCHASTIC ROUTINES *****************************
  96 !     ------------------------------------------------------------------
  97 ! This subroutine drives the initialization of the stochastic schemes
  98
  99       SUBROUTINE INITIALIZE_STOCH  (grid, config_flags,                       &
 100                           first_trip_for_this_domain,                         &
 101                           ips, ipe, jps, jpe, kps, kpe,                       &
 102                           ids, ide, jds, jde, kds, kde,                       &
 103                           ims, ime, jms, jme, kms, kme,                       &
 104                           its, ite, jts, jte, kts, kte,                       &
 105                           imsx, imex, jmsx, jmex, kmsx, kmex,  &
 106                           ipsx, ipex, jpsx, jpex, kpsx, kpex,  &
 107                           imsy, imey, jmsy, jmey, kmsy, kmey,  &
 108                           ipsy, ipey, jpsy, jpey, kpsy, kpey )
 109
 110
 111     USE module_configure
 112     USE module_domain, ONLY : domain
 113     USE module_wrf_error
 114 #ifdef DM_PARALLEL
 115     USE module_dm, ONLY : local_communicator, mytask, ntasks, ntasks_x, ntasks_y, local_communicator_periodic, &
 116                           wrf_dm_maxval, wrf_err_message, local_communicator_x, local_communicator_y, data_order_xzy
 117 #endif
 118
 119       IMPLICIT NONE
 120
 121       TYPE (grid_config_rec_type)            :: config_flags
 122       TYPE ( domain ), INTENT(INOUT)         :: grid
 123
 124       INTEGER , INTENT(IN)     ::               ids, ide, jds, jde, kds, kde,   &
 125                                                 ims, ime, jms, jme, kms, kme,   &
 126                                                 ips, ipe, jps, jpe, kps, kpe,   &
 127                                                 its, ite, jts, jte, kts, kte
 128       INTEGER , INTENT(IN)     ::               imsx,imex,jmsx,jmex,kmsx,kmex, &
 129                                                 ipsx,ipex,jpsx,jpex,kpsx,kpex, &
 130                                                 imsy,imey,jmsy,jmey,kmsy,kmey, &
 131                                                 ipsy,ipey,jpsy,jpey,kpsy,kpey
 132
 133       LOGICAL                  ::               first_trip_for_this_domain
 134       INTEGER                  ::               K, n, nfield, ierr, iseed, vertstruc, kvar
 135       real                     :: std, lengthscale, timescale, std_cutoff
 136       character (len = 10)    :: varname
 137       LOGICAL , EXTERNAL       :: wrf_dm_on_monitor
 138       INTEGER                  :: stochpert_unit
 139       INTEGER , PARAMETER      :: OPEN_OK = 0
 140       CHARACTER (len = 1024)   :: message
 141       logical, parameter :: DEBUG_MULTI_PERT = .true.
 142       REAL, DIMENSION(ims:ime,jms:jme) :: SPFORCS3d_loc, SPFORCC3d_loc, SP_AMP3d_loc
 143
 144
 145    IF ( first_trip_for_this_domain ) THEN
 146      grid%did_stoch = .FALSE.
 147    END IF
 148
 149    IF ((( grid%id == 1) .AND. (.NOT. grid%did_stoch)) .AND. &
 150        (( grid%skebs_on== 1) .OR.( grid%sppt_on== 1) .OR. ( grid%rand_perturb_on== 1) .OR. &
 151          ( grid%spp_conv== 1) .OR. ( grid%spp_pbl== 1) .OR. ( grid%spp_lsm== 1)) .or. &
 152         (config_flags%multi_perturb == 1)) THEN
 153
 154      grid%did_stoch = .TRUE.
 155
 156      IF (grid%skebs_on==1) then
 157
 158 ! Initialize SKEBS
 159 !    Initialize streamfunction (1)
 160      if ((.not.config_flags%restart) .or. (.not.config_flags%hrrr_cycling)) then
 161          call rand_seed (config_flags, grid%ISEED_SKEBS, grid%iseedarr_skebs , 1, config_flags%seed_dim)
 162      endif
 163      call SETUP_RAND_PERTURB('W',                                         &
 164                        grid%skebs_vertstruc,config_flags%restart,         &
 165                        grid%SPSTREAM_AMP,                                 &
 166                        grid%SPSTREAMFORCS,grid%SPSTREAMFORCC,grid%ALPH_PSI,&
 167                        grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPUV,     &
 168                        grid%KMINFORCT,grid%KMAXFORCT,                     &
 169                        grid%LMINFORCT,grid%LMAXFORCT,                     &
 170                        grid%KMAXFORCTH,grid%LMAXFORCTH,                   &
 171                        grid%time_step,grid%DX,grid%DY,                    &
 172                        grid%gridpt_stddev_sppt,                           &
 173                        grid%lengthscale_sppt,                             &
 174                        grid%timescale_sppt,                               &
 175                        grid%TOT_BACKSCAT_PSI,grid%ZTAU_PSI,               &
 176                        grid%REXPONENT_PSI,                                &
 177                        ids, ide, jds, jde, kds, kde,                      &
 178                        ims, ime, jms, jme, kms, kme,                      &
 179                        its, ite, jts, jte, kts, kte                       )
 180 !    Initialize potential temperature (2)
 181      call SETUP_RAND_PERTURB('T',                                         &
 182                        grid%skebs_vertstruc,config_flags%restart,     &
 183                        grid%SPT_AMP,                                      &
 184                        grid%SPTFORCS,grid%SPTFORCC,grid%ALPH_T,           &
 185                        grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT,     &
 186                        grid%KMINFORCT,grid%KMAXFORCT,                     &
 187                        grid%LMINFORCT,grid%LMAXFORCT,                     &
 188                        grid%KMAXFORCTH,grid%LMAXFORCTH,                   &
 189                        grid%time_step,grid%DX,grid%DY,                    &
 190                        grid%gridpt_stddev_sppt,                           &
 191                        grid%lengthscale_sppt,                             &
 192                        grid%timescale_sppt,                               &
 193                        grid%TOT_BACKSCAT_T,grid%ZTAU_T,                   &
 194                        grid%REXPONENT_T,                                  &
 195                        ids, ide, jds, jde, kds, kde,                      &
 196                        ims, ime, jms, jme, kms, kme,                      &
 197                        its, ite, jts, jte, kts, kte                       )
 198      ENDIF
 199
 200 IF (grid%sppt_on==1) then
 201 ! Initialize SPPT (3)
 202      if ((.not.config_flags%restart) .or. (.not.config_flags%hrrr_cycling)) then
 203          call rand_seed (config_flags, grid%ISEED_SPPT, grid%iseedarr_sppt  , 1, config_flags%seed_dim)
 204      endif
 205      call SETUP_RAND_PERTURB('P',                                         &
 206                        grid%sppt_vertstruc,config_flags%restart,          &
 207                        grid%SPPT_AMP,                                     &
 208                        grid%SPPTFORCC,grid%SPPTFORCS,grid%ALPH_SPPT,      &
 209                        grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT,     &
 210                        grid%KMINFORCT,grid%KMAXFORCT,                     &
 211                        grid%LMINFORCT,grid%LMAXFORCT,                     &
 212                        grid%KMAXFORCTH,grid%LMAXFORCTH,                   &
 213                        grid%time_step,grid%DX,grid%DY,                    &
 214                        grid%gridpt_stddev_sppt,                           &
 215                        grid%lengthscale_sppt,                             &
 216                        grid%timescale_sppt,                               &
 217                        grid%TOT_BACKSCAT_PSI,grid%ZTAU_PSI,               &
 218                        grid%REXPONENT_PSI,                                &
 219                        ids, ide, jds, jde, kds, kde,                      &
 220                        ims, ime, jms, jme, kms, kme,                      &
 221                        its, ite, jts, jte, kts, kte                       )
 222      ENDIF
 223
 224 ! Initialize RAND_PERTURB (4)
 225      IF (grid%rand_perturb_on==1) then
 226      if ((.not.config_flags%restart) .or. (.not.config_flags%hrrr_cycling)) then
 227          call rand_seed (config_flags, grid%ISEED_RAND_PERT, grid%iseedarr_rand_pert  , 1, config_flags%seed_dim)
 228      endif
 229      call SETUP_RAND_PERTURB('R',                                         &
 230                        grid%rand_pert_vertstruc,config_flags%restart,     &
 231                        grid%SP_AMP,                                       &
 232                        grid%SPFORCC,grid%SPFORCS,grid%ALPH_RAND,          &
 233                        grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT,     &
 234                        grid%KMINFORCT,grid%KMAXFORCT,                     &
 235                        grid%LMINFORCT,grid%LMAXFORCT,                     &
 236                        grid%KMAXFORCTH,grid%LMAXFORCTH,                   &
 237                        grid%time_step,grid%DX,grid%DY,                    &
 238                        grid%gridpt_stddev_rand_pert,                      &
 239                        grid%lengthscale_rand_pert,                        &
 240                        grid%timescale_rand_pert,                          &
 241                        grid%TOT_BACKSCAT_PSI,grid%ZTAU_PSI,               &
 242                        grid%REXPONENT_PSI,                                &
 243                        ids, ide, jds, jde, kds, kde,                      &
 244                        ims, ime, jms, jme, kms, kme,                      &
 245                        its, ite, jts, jte, kts, kte                       )
 246
 247      if (.not.config_flags%restart) then ! spin up
 248         do k = 1,10
 249            CALL RAND_PERT_UPDATE(grid,'R',                                     &
 250                            grid%SPFORCS,grid%SPFORCC,                          &
 251                            grid%SP_AMP,grid%ALPH_RAND,                         &
 252                            ips, ipe, jps, jpe, kps, kpe,                       &
 253                            ids, ide, jds, jde, kds, kde,                       &
 254                            ims, ime, jms, jme, kms, kme,                       &
 255                            kts, kte,                                           &
 256                            imsx,imex,jmsx,jmex,kmsx,kmex,                      &
 257                            ipsx,ipex,jpsx,jpex,kpsx,kpex,                      &
 258                            imsy,imey,jmsy,jmey,kmsy,kmey,                      &
 259                            ipsy,ipey,jpsy,jpey,kpsy,kpey,                      &
 260                            grid%num_stoch_levels,grid%num_stoch_levels,        &
 261                            grid%num_stoch_levels,grid%num_stoch_levels,        &
 262                            config_flags%restart, grid%iseedarr_rand_pert,      &
 263                            config_flags%seed_dim,                              &
 264                            grid%DX,grid%DY,grid%rand_pert_vertstruc,           &
 265                            grid%RAND_PERT,                                     &
 266                            grid%stddev_cutoff_rand_pert,                       &
 267                            grid%gridpt_stddev_rand_pert,                       &
 268                            grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT       )
 269          enddo
 270        ENDIF !rand_perturb_on
 271      ENDIF
 272
 273        ! Initialize MULTI_PERTURB
 274      If_multi_perturb: IF (config_flags%multi_perturb == 1) then
 275          ! Set frequency to update perturbations: spdt is in min and dt is in s
 276        grid%stepsp = nint (config_flags%spdt * 60.0 / grid%dt)
 277        grid%stepsp = max (grid%stepsp, 1)
 278
 279          ! Initialize perturbations
 280        stochpert_unit = 97
 281        IF ( wrf_dm_on_monitor() ) THEN
 282          OPEN(stochpert_unit, FILE='STOCHPERT.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
 283          IF ( ierr .NE. OPEN_OK ) THEN
 284            WRITE(message,FMT='(A)') &
 285            'module_stoch.F: INITIALIZE_STOCH: open failure for STOCHPERT.TBL'
 286            CALL wrf_error_fatal ( message )
 287          END IF
 288          REWIND(stochpert_unit)
 289        ENDIF
 290
 291        if (DEBUG_MULTI_PERT) print *, 'num_pert3d = ', num_pert3d
 292        if (DEBUG_MULTI_PERT) print *, 'num_pert_3d = ', config_flags%num_pert_3d
 293          ! Reads stoch table
 294          if ( wrf_dm_on_monitor() ) then
 295            READ (stochpert_unit,*)  ! Skip
 296            READ (stochpert_unit,*)  ! the
 297            READ (stochpert_unit,*)  ! header
 298          end if
 299        DO n = 2, config_flags%num_pert_3d
 300          IF ( wrf_dm_on_monitor() ) THEN
 301            read (stochpert_unit, *) nfield, varname, std, lengthscale, timescale, std_cutoff, &
 302                iseed, vertstruc
 303            select case (trim(varname))
 304              case ('ALBEDO')
 305                kvar = P_PALBEDO
 306
 307              case ('AOD')
 308                kvar = P_PAOD
 309
 310              case ('ANGSTROM')
 311                kvar = P_PANGSTROM
 312
 313              case ('ASSYMFAC')
 314                kvar = P_PASSYMFAC
 315
 316              case ('QVAPOR')
 317                kvar = P_PQVAPOR
 318
 319              case ('QCLOUD')
 320                kvar = P_PQCLOUD
 321
 322              case ('QICE')
 323                kvar = P_PQICE
 324
 325              case ('QSNOW')
 326                kvar = P_PQSNOW
 327
 328              case ('NI')
 329                kvar = P_PNI
 330
 331              case ('TH')
 332                kvar = P_PTH
 333
 334              case ('TKE')
 335                kvar = P_PTKE
 336
 337              case ('SMOIS')
 338                kvar = P_PSMOIS
 339
 340              case ('TSOIL')
 341                kvar = P_PTSOIL
 342
 343              case ('W')
 344                kvar = P_PW
 345
 346              case default
 347                WRITE(message,FMT='(A)') &
 348                    'module_stoch.F: Invalid entry in STOCHPERT.TBL'
 349                CALL wrf_error_fatal ( message )
 350
 351            end select
 352
 353            grid%gridpt_stddev_mult3d(kvar) = std
 354            grid%lengthscale_mult3d(kvar) = lengthscale
 355            grid%timescale_mult3d(kvar) = timescale
 356            grid%stddev_cutoff_mult3d(kvar) = std_cutoff
 357            grid%ISEED_MULT3D(kvar) = iseed
 358            grid%mult3d_vertstruc(kvar) = vertstruc
 359
 360            if (DEBUG_MULTI_PERT) print *, kvar, varname, grid%gridpt_stddev_mult3d(kvar), &
 361                              grid%lengthscale_mult3d(kvar), grid%timescale_mult3d(kvar), &
 362                              grid%stddev_cutoff_mult3d(kvar), grid%ISEED_MULT3D(kvar), &
 363                              grid%mult3d_vertstruc(kvar)
 364          ENDIF
 365            ! Note that num_pert_3d is the dimension of the 1d arrays from the table
 366            ! and num_pert3d is the number of allocated pert3d elements from the package (<= num_pert_3d)
 367          CALL wrf_dm_bcast_bytes (grid%gridpt_stddev_mult3d,   config_flags%num_pert_3d * RWORDSIZE )
 368          CALL wrf_dm_bcast_bytes (grid%lengthscale_mult3d,     config_flags%num_pert_3d * RWORDSIZE )
 369          CALL wrf_dm_bcast_bytes (grid%timescale_mult3d,       config_flags%num_pert_3d * RWORDSIZE )
 370          CALL wrf_dm_bcast_bytes (grid%stddev_cutoff_mult3d,   config_flags%num_pert_3d * RWORDSIZE )
 371          CALL wrf_dm_bcast_bytes (grid%ISEED_MULT3D,           config_flags%num_pert_3d * IWORDSIZE )
 372          CALL wrf_dm_bcast_bytes (grid%mult3d_vertstruc,       config_flags%num_pert_3d * IWORDSIZE )
 373        END DO
 374
 375        Loop_multi_perturb: DO n = 2, config_flags%num_pert_3d
 376          CALL contiguize_2d ( .TRUE. , grid%SPFORCS3d, SPFORCS3d_loc, n,       &
 377                               ips, ipe, jps, jpe, 1, config_flags%num_pert_3d, &
 378                               ids, ide, jds, jde, 1, config_flags%num_pert_3d, &
 379                               ims, ime, jms, jme, 1, config_flags%num_pert_3d  )
 380          CALL contiguize_2d ( .TRUE. , grid%SPFORCC3d, SPFORCC3d_loc, n,       &
 381                               ips, ipe, jps, jpe, 1, config_flags%num_pert_3d, &
 382                               ids, ide, jds, jde, 1, config_flags%num_pert_3d, &
 383                               ims, ime, jms, jme, 1, config_flags%num_pert_3d  )
 384          CALL contiguize_2d ( .TRUE. , grid%SP_AMP3d,  SP_AMP3d_loc,  n,       &
 385                               ips, ipe, jps, jpe, 1, config_flags%num_pert_3d, &
 386                               ids, ide, jds, jde, 1, config_flags%num_pert_3d, &
 387                               ims, ime, jms, jme, 1, config_flags%num_pert_3d  )
 388          if ((.not.config_flags%restart) .or. (.not.config_flags%hrrr_cycling)) then
 389            call rand_seed (config_flags, grid%ISEED_MULT3D(n), grid%iseedarr_mult3d(:,n)  , 1, config_flags%seed_dim)
 390          endif
 391          call SETUP_RAND_PERTURB('R',                                     &
 392                        grid%mult3d_vertstruc(n),config_flags%restart,     &
 393                        SP_AMP3d_loc,SPFORCC3d_loc,SPFORCS3d_loc,          &
 394                        grid%ALPH_RAND3d(n),                               &
 395                        grid%VERTSTRUCC3d(:,:,:,n),grid%VERTSTRUCS3d(:,:,:,n),grid%VERTAMPT3d(:,n),     &
 396                        grid%KMINFORCT,grid%KMAXFORCT,                     &
 397                        grid%LMINFORCT,grid%LMAXFORCT,                     &
 398                        grid%KMAXFORCTH,grid%LMAXFORCTH,                   &
 399                        grid%stepsp * grid%time_step,grid%DX,grid%DY,      &
 400                        grid%gridpt_stddev_mult3d(n),                      &
 401                        grid%lengthscale_mult3d(n),                        &
 402                        grid%timescale_mult3d(n),                          &
 403                        grid%TOT_BACKSCAT_PSI,grid%ZTAU_PSI,               &
 404                        grid%REXPONENT_PSI,                                &
 405                        ids, ide, jds, jde, kds, kde,                      &
 406                        ims, ime, jms, jme, kms, kme,                      &
 407                        its, ite, jts, jte, kts, kte                       )
 408
 409 !        if (.not.config_flags%restart) then ! spin up
 410 !          do k = 1,10
 411 !            CALL RAND_PERT_UPDATE(grid,'R',                                   &
 412 !                          SPFORCS3d_loc,SPFORCC3d_loc,SP_AMP3d_loc,           &
 413 !                          grid%ALPH_RAND3d(n),                                &
 414 !                          ips, ipe, jps, jpe, kps, kpe,                       &
 415 !                          ids, ide, jds, jde, kds, kde,                       &
 416 !                          ims, ime, jms, jme, kms, kme,                       &
 417 !                          kts, kte,                                           &
 418 !                          imsx,imex,jmsx,jmex,kmsx,kmex,                      &
 419 !                          ipsx,ipex,jpsx,jpex,kpsx,kpex,                      &
 420 !                          imsy,imey,jmsy,jmey,kmsy,kmey,                      &
 421 !                          ipsy,ipey,jpsy,jpey,kpsy,kpey,                      &
 422 !                          grid%num_stoch_levels,grid%num_stoch_levels,        &
 423 !                          grid%num_stoch_levels,grid%num_stoch_levels,        &
 424 !                          config_flags%restart, grid%iseedarr_mult3d(:,n),    &
 425 !                          config_flags%seed_dim,                              &
 426 !                          grid%DX,grid%DY,grid%mult3d_vertstruc(n),           &
 427 !                          grid%PERT3D(:,:,:,n),                               &
 428 !                          grid%stddev_cutoff_mult3d(n),                       &
 429 !                          grid%gridpt_stddev_mult3d(n),                       &
 430 !                          grid%VERTSTRUCC3d(:,:,:,n),grid%VERTSTRUCS3d(:,:,:,n),grid%VERTAMPT3d(:,n)       )
 431 !          enddo
 432 !        end if
 433          CALL contiguize_2d ( .FALSE., grid%SPFORCS3d, SPFORCS3d_loc, n,       &
 434                               ips, ipe, jps, jpe, 1, config_flags%num_pert_3d, &
 435                               ids, ide, jds, jde, 1, config_flags%num_pert_3d, &
 436                               ims, ime, jms, jme, 1, config_flags%num_pert_3d  )
 437          CALL contiguize_2d ( .FALSE., grid%SPFORCC3d, SPFORCC3d_loc, n,       &
 438                               ips, ipe, jps, jpe, 1, config_flags%num_pert_3d, &
 439                               ids, ide, jds, jde, 1, config_flags%num_pert_3d, &
 440                               ims, ime, jms, jme, 1, config_flags%num_pert_3d  )
 441          CALL contiguize_2d ( .FALSE., grid%SP_AMP3d,  SP_AMP3d_loc,  n,       &
 442                               ips, ipe, jps, jpe, 1, config_flags%num_pert_3d, &
 443                               ids, ide, jds, jde, 1, config_flags%num_pert_3d, &
 444                               ims, ime, jms, jme, 1, config_flags%num_pert_3d  )
 445        END DO Loop_multi_perturb
 446      ENDIF  If_multi_perturb
 447      IF ( wrf_dm_on_monitor() ) CLOSE (stochpert_unit)
 448
 449
 450 ! Initialize Stochastic Parameter Perturbations to convection scheme
 451      IF (grid%spp_conv==1) then
 452      if (.not.config_flags%restart) then ! set random number seed (else  iseedarray is read in from restart files)
 453          call rand_seed (config_flags, grid%iseed_spp_conv, grid%iseedarr_spp_conv , 1, config_flags%seed_dim)
 454      endif
 455      call SETUP_RAND_PERTURB('S',                                         &
 456                        grid%vertstruc_spp_conv,config_flags%restart,      &
 457                        grid%SP_AMP2,                                      &
 458                        grid%SPFORCC2,grid%SPFORCS2,grid%ALPH_RAND2,       &
 459                        grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT,     &
 460                        grid%KMINFORCT,grid%KMAXFORCT,                     &
 461                        grid%LMINFORCT,grid%LMAXFORCT,                     &
 462                        grid%KMAXFORCTH,grid%LMAXFORCTH,                   &
 463                        grid%time_step,grid%DX,grid%DY,                    &
 464                        grid%gridpt_stddev_spp_conv,                       &
 465                        grid%lengthscale_spp_conv,                         &
 466                        grid%timescale_spp_conv,                           &
 467                        grid%TOT_BACKSCAT_PSI,grid%ZTAU_PSI,               &
 468                        grid%REXPONENT_PSI,                                &
 469                        ids, ide, jds, jde, kds, kde,                      &
 470                        ims, ime, jms, jme, kms, kme,                      &
 471                        its, ite, jts, jte, kts, kte                       )
 472      ENDIF
 473 ! Initialize Stochastic Parameter Peturbations (SPP) to PBL scheme
 474      IF (grid%spp_pbl==1) then
 475      if (.not.config_flags%restart) then ! set random number seed (else  iseedarray is read in from restart files)
 476          call rand_seed (config_flags, grid%iseed_spp_pbl, grid%iseedarr_spp_pbl , 1, config_flags%seed_dim)
 477      endif
 478      call SETUP_RAND_PERTURB('Q',                                         &
 479                        grid%vertstruc_spp_pbl,config_flags%restart,       &
 480                        grid%SP_AMP3,                                      &
 481                        grid%SPFORCC3,grid%SPFORCS3,grid%ALPH_RAND3,       &
 482                        grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT,     &
 483                        grid%KMINFORCT,grid%KMAXFORCT,                     &
 484                        grid%LMINFORCT,grid%LMAXFORCT,                     &
 485                        grid%KMAXFORCTH,grid%LMAXFORCTH,                   &
 486                        grid%time_step,grid%DX,grid%DY,                    &
 487                        grid%gridpt_stddev_spp_pbl,                        &
 488                        grid%lengthscale_spp_pbl,                          &
 489                        grid%timescale_spp_pbl,                            &
 490                        grid%TOT_BACKSCAT_PSI,grid%ZTAU_PSI,               &
 491                        grid%REXPONENT_PSI,                                &
 492                        ids, ide, jds, jde, kds, kde,                      &
 493                        ims, ime, jms, jme, kms, kme,                      &
 494                        its, ite, jts, jte, kts, kte                       )
 495      ENDIF
 496 ! Initialize Stochastic Parameter Peturbations (SPP) to LSM scheme
 497      IF (grid%spp_lsm==1) then
 498      if (.not.config_flags%restart) then ! set random number seed (else  iseedarray is read in from restart files)
 499          call rand_seed (config_flags, grid%iseed_spp_lsm, grid%iseedarr_spp_lsm , 1, config_flags%seed_dim)
 500      endif
 501      call SETUP_RAND_PERTURB('O',                                         &
 502                        grid%vertstruc_spp_lsm,config_flags%restart,       &
 503                        grid%SP_AMP4,                                      &
 504                        grid%SPFORCC4,grid%SPFORCS4,grid%ALPH_RAND4,       &
 505                        grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT,     &
 506                        grid%KMINFORCT,grid%KMAXFORCT,                     &
 507                        grid%LMINFORCT,grid%LMAXFORCT,                     &
 508                        grid%KMAXFORCTH,grid%LMAXFORCTH,                   &
 509                        grid%time_step,grid%DX,grid%DY,                    &
 510                        grid%gridpt_stddev_spp_lsm,                        &
 511                        grid%lengthscale_spp_lsm,                          &
 512                        grid%timescale_spp_lsm,                            &
 513                        grid%TOT_BACKSCAT_PSI,grid%ZTAU_PSI,               &
 514                        grid%REXPONENT_PSI,                                &
 515                        ids, ide, jds, jde, kds, kde,                      &
 516                        ims, ime, jms, jme, kms, kme,                      &
 517                        its, ite, jts, jte, kts, kte                       )
 518      ENDIF
 519
 520      ENDIF ! skebs or sppt or rand_perturb or spp
 521
 522       END SUBROUTINE INITIALIZE_STOCH
 523
 524       !   --- END SETUP STOCHASTIC PERTURBATION SCHEMES ----------
 525
 526
 527       subroutine SETUP_RAND_PERTURB( variable_in,&
 528                        skebs_vertstruc,restart,                  &
 529                        SP_AMP,SPFORCC,SPFORCS,ALPH,                  &
 530                        VERTSTRUCC,VERTSTRUCS,VERTAMP,                &
 531                        KMINFORCT,KMAXFORCTH,LMINFORCT,LMAXFORCTH,    &
 532                        KMAXFORCT,LMAXFORCT,                          &
 533                        itime_step,DX,DY,                             &
 534                        gridpt_stddev_rand_perturb, l_rand_perturb,   &
 535                        tau_rand_perturb,                             &
 536                        TOT_BACKSCAT,ZTAU,REXPONENT,                  &
 537                        ids, ide, jds, jde, kds, kde,                 &
 538                        ims, ime, jms, jme, kms, kme,                 &
 539                        its, ite, jts, jte, kts, kte                  )
 540
 541
 542
 543       IMPLICIT NONE
 544
 545 ! General control
 546       LOGICAL                                   :: restart
 547       REAL, PARAMETER                           :: RPI= 3.141592653589793 !4.0*atan(1.0)
 548       CHARACTER, INTENT(IN)                     :: variable_in ! W=SKEBS_PSI, T=SKEBS_T, P=SPPT, R=RAND_PERTURB
 549       CHARACTER                                 :: variable
 550
 551 ! Common to all schemes
 552       INTEGER , INTENT(IN)                      :: ids, ide, jds, jde, kds, kde,   &
 553                                                    ims, ime, jms, jme, kms, kme,   &
 554                                                    its, ite, jts, jte, kts, kte
 555       INTEGER                                   :: IER,IK,IL,I,J,itime_step,skebs_vertstruc, &
 556                                                    KMINFORCT,LMINFORCT,KMAXFORCT,LMAXFORCT,KMAXFORCTH,LMAXFORCTH, &
 557                                                    KMAX,LMAX,LENSAV,ILEV
 558       REAL                                      :: DX,DY,RY,RX,ALPH,RHOKLMAX,ZREF,RHOKL,EPS
 559       REAL, DIMENSION (ims:ime,jms:jme)         :: SPFORCS,SPFORCC,SP_AMP
 560       REAL, DIMENSION (ims:ime,kms:kme,jms:jme) :: VERTSTRUCC,VERTSTRUCS
 561       REAL, DIMENSION (kms:kme)                 :: VERTAMP
 562       REAL, DIMENSION (ids:ide,jds:jde)         :: ZCHI
 563
 564 ! SPPT and perturb_rand specific
 565       REAL                                      :: gridpt_stddev_rand_perturb,kappat,tau_rand_perturb,l_rand_perturb
 566       REAL, DIMENSION (ims:ime,jms:jme)         :: var_sigma1
 567
 568
 569 ! SKEBS specific
 570       REAL                                      :: z,phi,ZGAMMAN,ZCONSTF0,TOT_BACKSCAT,ZTAU,REXPONENT,ZSIGMA2
 571       LOGICAL                                   :: is_print = .true.
 572
 573
 574       variable = variable_in
 575 !     --------- SETUP PARAMETERS ---------------------------------------
 576       KMAX=(jde-jds)+1 !NLAT
 577       LMAX=(ide-ids)+1 !NLON
 578       RY=  KMAX*DY
 579       RX=  LMAX*DX
 580       LENSAV= 4*(KMAX+LMAX)+INT(LOG(REAL(KMAX))) + INT(LOG(REAL(LMAX))) + 8
 581
 582 !     --------- ALLOCATE FIELDS FOR FFTPACK----------------------------
 583 !     --------- ALLOCATE FIELDS FOR FFTPACK----------------------------
 584       IF ( ALLOCATED(WSAVE1)      ) DEALLOCATE(WSAVE1)
 585       IF ( ALLOCATED(WSAVE2)      ) DEALLOCATE(WSAVE2)
 586       ALLOCATE(WSAVE1(LENSAV),WSAVE2(LENSAV))
 587
 588       IF ( ALLOCATED(WAVENUMBER_K)) DEALLOCATE(WAVENUMBER_K)
 589       IF ( ALLOCATED(WAVENUMBER_L)) DEALLOCATE(WAVENUMBER_L)
 590       ALLOCATE (wavenumber_k(jds:jde),wavenumber_l(ids:ide))
 591
 592 !     -------- INITIALIZE FFTPACK ROUTINES -----------------------------
 593       call CFFT1I (LMAX, WSAVE1, LENSAV, IER)
 594       if(ier.ne. 0) write(*,95) ier
 595
 596       call CFFT1I (KMAX, WSAVE2, LENSAV, IER)
 597       if(ier.ne. 0) write(*,95) ier
 598
 599       95 format('error in cFFT2I=  ',i5)
 600
 601       call findindex( wavenumber_k, wavenumber_l,             &
 602                       ids, ide, jds, jde, kds, kde,                &
 603                       ims, ime, jms, jme, kms, kme,                &
 604                       its, ite, jts, jte, kts, kte                 )
 605
 606 !    set maximal perturbed wavenumber based on gridpoints in domain
 607      KMAXFORCT=min0(((ide-ids)+1)/2,((jde-jds)+1 )/2)-5
 608      LMAXFORCT=KMAXFORCT
 609      if (KMAXFORCT > KMAXFORCTH) then
 610         KMAXFORCT=KMAXFORCTH
 611      endif
 612      if (LMAXFORCT > LMAXFORCTH) then
 613         LMAXFORCT=LMAXFORCTH
 614      endif
 615
 616
 617 !     --------------------------------------------------------------------------------------
 618 !     ---------- INITIALIZE STOCHASTIC KINETIC-ENERGY BACKSCATTER SCHEME  (SKEBS) ----------
 619 !     --------------------------------------------------------------------------------------
 620       ALPH   =  float(itime_step)/ZTAU ! approximation of 1.-exp(-itime_step/ZTAU_PSI)
 621       ZSIGMA2=1./(12.0*ALPH)
 622
 623       if (is_print) then
 624       IF (variable == 'W') then
 625       WRITE(*,'(''                                               '')')
 626       WRITE(*,'('' =============================================='')')
 627       WRITE(*,'('' >> Initializing STREAMFUNCTION forcing pattern of  << '')')
 628       WRITE(*,'('' >> stochastic kinetic-energy backscatter scheme << '')')
 629       WRITE(*,'('' Total backscattered energy, TOT_BACKSCAT_PSI '',E12.5)') TOT_BACKSCAT
 630       WRITE(*,'('' Exponent for energy spectra, REXPONENT_PSI ='',E12.5)') REXPONENT
 631       WRITE(*,'('' Minimal wavenumber of streamfunction forcing, LMINFORC ='',I10)') LMINFORCT
 632       WRITE(*,'('' Maximal wavenumber of streamfunction forcing, LMAXFORC ='',I10)') LMAXFORCT
 633       WRITE(*,'('' Minimal wavenumber of streamfunction forcing, KMINFORC ='',I10)') KMINFORCT
 634       WRITE(*,'('' Maximal wavenumber of streamfunction forcing, KMAXFORC ='',I10)') KMAXFORCT
 635       WRITE(*,'('' skebs_vertstruc                             '',I10)') skebs_vertstruc
 636       WRITE(*,'('' Time step: itime_step='',I10)') itime_step
 637       WRITE(*,'('' Decorrelation time of noise, ZTAU_PSI ='',E12.5)') ZTAU
 638       WRITE(*,'('' Variance of noise, ZSIGMA2_EPS  ='',E12.5)') ZSIGMA2
 639       WRITE(*,'('' Autoregressive parameter 1-ALPH_PSI ='',E12.5)') 1.-ALPH
 640       WRITE(*,'('' =============================================='')')
 641
 642 !     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
 643       ELSEIF (variable == 'T') then
 644       WRITE(*,'(''                                               '')')
 645       WRITE(*,'('' =============================================='')')
 646       WRITE(*,'('' >> Initializing TEMPERATURE forcing pattern of  << '')')
 647       WRITE(*,'('' >> stochastic kinetic-energy backscatter scheme << '')')
 648       WRITE(*,'('' Total backscattered energy, TOT_BACKSCAT_T   '',E12.5)') TOT_BACKSCAT
 649       WRITE(*,'('' Exponent for energy spectra, REXPONENT_T   ='',E12.5)') REXPONENT
 650       WRITE(*,'('' Minimal wavenumber of tempearature forcing, LMINFORC ='',I10)') LMINFORCT
 651       WRITE(*,'('' Maximal wavenumber of tempearature forcing, LMAXFORC ='',I10)') LMAXFORCT
 652       WRITE(*,'('' Minimal wavenumber of tempearature forcing, KMINFORC ='',I10)') KMINFORCT
 653       WRITE(*,'('' Maximal wavenumber of tempearature forcing, KMAXFORC ='',I10)') KMAXFORCT
 654       WRITE(*,'('' skebs_vertstruc                             '',I10)') skebs_vertstruc
 655       WRITE(*,'('' Decorrelation time of noise, ZTAU_T ='',E12.5)') ZTAU
 656       WRITE(*,'('' Variance of noise, ZSIGMA2_ETA  ='',E12.5)') ZSIGMA2
 657       WRITE(*,'('' Autoregressive parameter 1-ALPH_T ='',E12.5)') 1.-ALPH
 658       WRITE(*,'('' =============================================='')')
 659       endif
 660
 661      IF ((variable == 'P') .or. (variable == 'R') .or. (variable == 'S') .or. (variable == 'Q') .or. (variable == 'O')) then
 662       kappat= L_rand_perturb**2 ! L^2= kappa*T,  where L is a length scale in m; set to for L=100km
 663       phi = exp (-float(itime_step)/tau_rand_perturb)
 664       alph = 1.-phi
 665      endif
 666
 667 !     --------------------------------------------------------------------------------------
 668 !     ---------- INITIALIZE STOCHASTICALLY PERTURBED PHYSICAL TENDENCY SCHEME --------------
 669 !     --------------------------------------------------------------------------------------
 670      if (variable == 'P') then
 671       WRITE(*,'(''                                               '')')
 672       WRITE(*,'('' =============================================='')')
 673       WRITE(*,'('' >> Initializing Stochastically Perturbed Physics Tendency scheme << '')')
 674       WRITE(*,'('' sppt_vertstruc                             '',I10)') skebs_vertstruc
 675       WRITE(*,'('' Decorrelation time of noise, Tau ='',E12.5)') tau_rand_perturb
 676       WRITE(*,'('' Autoregressive parameter Phi ='',E12.5)') phi
 677       WRITE(*,'('' Length Scale L'',E12.5)') l_rand_perturb
 678       WRITE(*,'('' Variance in gridpoint space'',E12.5)') gridpt_stddev_rand_perturb
 679       WRITE(*,'('' =============================================='')')
 680       endif ! variable
 681
 682 !     --------------------------------------------------------------------------------------
 683 !     --------------------  INITIALIZE  RANDOM PERTUBATIONS  -------------------------------
 684 !     --------------------------------------------------------------------------------------
 685      if (variable == 'R') then
 686       WRITE(*,'(''                                               '')')
 687       WRITE(*,'('' =============================================='')')
 688       WRITE(*,'('' >> Initializing random perturbations << '')')
 689       WRITE(*,'('' rand_pert_vertstruc                             '',I10)') skebs_vertstruc
 690       WRITE(*,'('' Decorrelation time of noise, Tau ='',E12.5)') tau_rand_perturb
 691       WRITE(*,'('' Autoregressive parameter Phi ='',E12.5)') phi
 692       WRITE(*,'('' Length Scale L'',E12.5)') l_rand_perturb
 693       WRITE(*,'('' Variance in gridpoint space'',E12.5)') gridpt_stddev_rand_perturb
 694       WRITE(*,'('' =============================================='')')
 695      endif ! variable
 696
 697      if (variable == 'S') then
 698       WRITE(*,'(''                                               '')')
 699       WRITE(*,'('' =============================================='')')
 700       WRITE(*,'('' >> Initializing stochastic parameter perturbations for convection<< '')')
 701       WRITE(*,'('' rand_pert_vertstruc2                            '',I10)') skebs_vertstruc
 702       WRITE(*,'('' Decorrelation time of noise, Tau ='',E12.5)') tau_rand_perturb
 703       WRITE(*,'('' Autoregressive parameter Phi ='',E12.5)') phi
 704       WRITE(*,'('' Length Scale L'',E12.5)') l_rand_perturb
 705       WRITE(*,'('' Variance in gridpoint space'',E12.5)') gridpt_stddev_rand_perturb
 706       WRITE(*,'('' =============================================='')')
 707      endif ! variable
 708
 709      if (variable == 'Q') then
 710       WRITE(*,'(''                                               '')')
 711       WRITE(*,'('' =============================================='')')
 712       WRITE(*,'('' >> Initializing stochastic parameter perturbations for PBL<< '')')
 713       WRITE(*,'('' rand_pert_vertstruc3                            '',I10)') skebs_vertstruc
 714       WRITE(*,'('' Decorrelation time of noise, Tau ='',E12.5)') tau_rand_perturb
 715       WRITE(*,'('' Autoregressive parameter Phi ='',E12.5)') phi
 716       WRITE(*,'('' Length Scale L'',E12.5)') l_rand_perturb
 717       WRITE(*,'('' Variance in gridpoint space'',E12.5)') gridpt_stddev_rand_perturb
 718       WRITE(*,'('' =============================================='')')
 719      endif ! variable
 720
 721      if (variable == 'O') then
 722       WRITE(*,'(''                                               '')')
 723       WRITE(*,'('' =============================================='')')
 724       WRITE(*,'('' >> Initializing stochastic parameter perturbations for LSM<< '')')
 725       WRITE(*,'('' rand_pert_vertstruc4                            '',I10)') skebs_vertstruc
 726       WRITE(*,'('' Decorrelation time of noise, Tau ='',E12.5)') tau_rand_perturb
 727       WRITE(*,'('' Autoregressive parameter Phi ='',E12.5)') phi
 728       WRITE(*,'('' Length Scale L'',E12.5)') l_rand_perturb
 729       WRITE(*,'('' Variance in gridpoint space'',E12.5)') gridpt_stddev_rand_perturb
 730       WRITE(*,'('' =============================================='')')
 731      endif ! variable
 732
 733      endif !is print
 734 !     --------------------------------------------------------------------------------------
 735 !     Compute Normalization constants
 736 !     --------------------------------------------------------------------------------------
 737
 738      ZCHI    =  0.0
 739      ZGAMMAN  =  0.0
 740       ! Fill lower left quadrant of ZCHI. For this range the indeces IK,IL
 741       DO IK=jds-1,jde   ! These are now wavenumbers
 742       DO IL=ids-1,ide
 743       if (((sqrt((IK/RY*IK/RY)+(IL/RX*IL/RX)).lt.((KMAXFORCT+0.5)/RX)).and.&
 744            (sqrt((IK/RY*IK/RY)+(IL/RX*IL/RX)).ge.((KMINFORCT-0.5)/RX))) .or. &
 745           ((sqrt((IK/RY*IK/RY)+(IL/RX*IL/RX)).lt.((LMAXFORCT+0.5)/RY)).and.&
 746            (sqrt((IK/RY*IK/RY)+(IL/RX*IL/RX)).ge.((LMINFORCT-0.5)/RY))))then
 747         if ((IK>0).or.(IL>0)) then
 748           if (variable == 'W') then
 749             ZCHI(IL+1,IK+1)=((IK/RY*IK/RY)+(IL/RX*IL/RX))**(REXPONENT/2.)  ! SKEBS :U
 750             ZGAMMAN= ZGAMMAN + ((IK/RY*IK/RY)+(IL/RX*IL/RX))**(REXPONENT+1)
 751           else if (variable == 'T') then
 752             ZCHI(IL+1,IK+1)=((IK/RY*IK/RY)+(IL/RX*IL/RX))**(REXPONENT/2.)  ! SKEBS :T
 753             ZGAMMAN= ZGAMMAN + ((IK/RY*IK/RY)+(IL/RX*IL/RX))**(REXPONENT)
 754           else if ((variable == 'P') .or. (variable == 'R') .or.  (variable == 'S') .or. (variable == 'O') .or.  (variable == 'Q     ')) then
 755             ZCHI(IL+1,IK+1)=exp(  -2*RPI**2*kappat*((IK/RY*IK/RY)+(IL/RX*IL/RX)) ) !SPPT
 756             ZGAMMAN= ZGAMMAN + exp(  -4*RPI**2*kappat*((IK/RY*IK/RY)+(IL/RX*IL/RX)) ) !SPPT
 757           endif
 758         endif
 759       endif
 760       enddo
 761       enddo
 762       ZGAMMAN=4.0*ZGAMMAN !account for all quadrants, although only one is Filled
 763       if (variable == 'W') then
 764          ZCONSTF0=SQRT(ALPH*TOT_BACKSCAT/(float(itime_step)*ZSIGMA2*ZGAMMAN))/(2*RPI)
 765       elseif  (variable == 'T') then
 766          ZCONSTF0=SQRT(T0*ALPH*TOT_BACKSCAT/(float(itime_step)*cp*ZSIGMA2*ZGAMMAN))
 767       elseif ((variable == 'P') .or. (variable == 'R') .or.  (variable == 'S') .or. (variable == 'O') .or.  (variable == 'Q     ')) then
 768          ZCONSTF0= gridpt_stddev_rand_perturb*sqrt((1.-phi**2)/(2.*ZGAMMAN))
 769       endif
 770
 771
 772 !     --------------------------------------------------------------------------------------
 773 !     Now the wavenumber-dependent amplitudes
 774 !     --------------------------------------------------------------------------------------
 775 !     Note: There are symmetries and anti-symmetries to ensure real-valued back transforms
 776 !     Fill lower left quadrant of matrix of noise amplitudes for wavenumbers K=0,KMAX/2
 777       SP_AMP=0.0
 778       DO IK=jts,jte
 779       DO IL=its,ite
 780       if ((IL .le. (LMAX/2+1)) .and. (IK .le. (KMAX/2+1)) ) then
 781         SP_AMP(IL,IK)      = ZCONSTF0*ZCHI(IL,IK)
 782       endif
 783       ENDDO
 784       ENDDO
 785
 786       ! Fill other quadrants:
 787       ! Upper left quadrant
 788       DO IK=jts,jte
 789       DO IL=its,ite
 790       if ( (IL .gt. (LMAX/2+1)) .and. (IK .le. (KMAX/2+1)) ) then
 791         SP_AMP(IL,IK)      =  ZCONSTF0*ZCHI(LMAX-IL+2,IK)
 792       endif
 793       ENDDO
 794       ENDDO
 795
 796 !     Lower right quadrant
 797       DO IK=jts,jte
 798       DO IL=its,ite
 799       if ((IK .gt. (KMAX/2+1)) .and. (IL.le.LMAX/2) ) then
 800         SP_AMP(IL,IK)      = ZCONSTF0*ZCHI(IL,KMAX-IK+2)
 801       endif
 802       ENDDO
 803       ENDDO
 804
 805 !     Upper right quadrant
 806       DO IK=jts,jte
 807       DO IL=its,ite
 808       if ((IK .gt. (KMAX/2+1)) .and. (IL.gt.LMAX/2) ) then
 809         SP_AMP(IL,IK)      = ZCONSTF0*ZCHI(LMAX-IL+2,KMAX-IK+2)
 810       endif
 811       ENDDO
 812       ENDDO
 813
 814 !     -----------------------------------------
 815 !     Array for vertical structure if desired
 816       VERTAMP=1.0 ! Define vertical amplitude here.
 817
 818       IF (skebs_vertstruc==1) then
 819         VERTSTRUCC=0.0
 820         VERTSTRUCS=0.0
 821         RHOKLMAX= sqrt(KMAX**2/DY**2 + LMAX**2/DX**2)
 822         ZREF=32.0
 823         DO ILEV=kts,kte
 824           DO IK=jts,jte
 825             DO IL=its,ite
 826             if (IL.le.(LMAX/2)) then
 827               RHOKL   = sqrt((IK+1)**2/DY**2 + (IL+1)**2/DX**2)
 828               EPS     = ((RHOKLMAX - RHOKL)/ RHOKLMAX)  * (ILEV/ZREF) * RPI
 829               VERTSTRUCC(IL,ILEV,IK) =  cos ( eps* (IL+1) )
 830               VERTSTRUCS(IL,ILEV,IK) =  sin ( eps* (IL+1) )
 831              else
 832               RHOKL   = sqrt((IK+1)**2/DY**2 + (LMAX-IL+2)**2/DX**2)
 833               EPS     = ((RHOKLMAX - RHOKL)/ RHOKLMAX)  * (ILEV/ZREF) * RPI
 834               VERTSTRUCC (IL,ILEV,IK) =   cos ( eps* (LMAX-IL+2) )
 835               VERTSTRUCS (IL,ILEV,IK) = - sin ( eps* (LMAX-IL+2) )
 836             endif
 837             ENDDO
 838           ENDDO
 839         ENDDO
 840       ENDIF
 841
 842      END subroutine SETUP_RAND_PERTURB
 843
 844 !     ------------------------------------------------------------------
 845 !************** UPDATE STOCHASTIC PATTERN IN WAVENUMBER SPACE**********
 846 !     ------------------------------------------------------------------
 847
 848      subroutine UPDATE_STOCH(                                         &
 849                       SPFORCS,SPFORCC,SP_AMP,ALPH,                    &
 850                       restart,iseedarr,seed_dim,                      &
 851                       ids, ide, jds, jde, kds, kde,                   &
 852                       ims, ime, jms, jme, kms, kme,                   &
 853                       its, ite, jts, jte, kts, kte                    )
 854      IMPLICIT NONE
 855
 856      REAL, DIMENSION( ids:ide,jds:jde)      :: ZRANDNOSS,ZRANDNOSC
 857      REAL, DIMENSION (ims:ime,jms:jme)      :: SPFORCS,SPFORCC,SP_AMP
 858      INTEGER , INTENT(IN)     ::               ids, ide, jds, jde, kds, kde,   &
 859                                                ims, ime, jms, jme, kms, kme,   &
 860                                                its, ite, jts, jte, kts, kte
 861      INTEGER , INTENT(IN)     ::               seed_dim
 862
 863      INTEGER, DIMENSION (seed_dim) , INTENT(INOUT) :: iseedarr
 864      INTEGER , ALLOCATABLE , DIMENSION(:) :: iseed
 865      REAL :: Z,ALPH
 866      REAL, PARAMETER :: thresh = 3.0
 867      INTEGER ::IL, IK,LMAX,KMAX
 868      INTEGER :: how_many
 869      LOGICAL :: LGAUSS,RESTART
 870
 871      KMAX=(jde-jds)+1 !NLAT
 872      LMAX=(ide-ids)+1 !NATX
 873
 874            CALL random_seed(size=how_many)
 875            IF ( ALLOCATED(iseed)) DEALLOCATE(iseed)
 876            ALLOCATE(iseed(how_many))
 877            iseed=iseedarr(1:how_many)
 878            call random_seed(put=iseed(1:how_many))
 879
 880 !     Pick the distribution of the noise
 881 !     Random noise uses global indexes to ensure necessary symmetries and anti-symmetries
 882 !     of random forcing when run on multiple processors
 883      LGAUSS=.true.
 884      IF (LGAUSS) then
 885        DO IK=jds,jde
 886          DO IL=ids,ide
 887           do
 888            call gauss_noise(z)
 889            if (abs(z)<thresh) exit
 890           ENDDO
 891           ZRANDNOSS(IL,IK)=z
 892           do
 893            call gauss_noise(z)
 894            if (abs(z)<thresh) exit
 895           ENDDO
 896           ZRANDNOSC(IL,IK)=z
 897          ENDDO
 898        ENDDO
 899      ELSE
 900        DO IK=jds,jde
 901          DO IL=ids,ide
 902            CALL RANDOM_NUMBER(z)
 903            ZRANDNOSS(IL,IK)=z-0.5
 904            CALL RANDOM_NUMBER(z)
 905            ZRANDNOSC(IL,IK)=z-0.5
 906           ENDDO
 907         ENDDO
 908       ENDIF
 909
 910 !     Note: There are symmetries and anti-symmetries to ensure real-valued back transforms
 911 ! for symmetric part: left and right half axis symmetric
 912
 913       DO IK=jts,jte
 914       if ((IK.le.(KMAX/2+1)) .and. (IK>1)) then ! Upper half
 915         DO IL=its,ite
 916           SPFORCC(IL,IK)       = (1.-ALPH)*SPFORCC(IL,IK)      + SP_AMP(IL,IK)     * ZRANDNOSC(IL,IK)
 917           SPFORCS(IL,IK)       = (1.-ALPH)*SPFORCS(IL,IK)      + SP_AMP(IL,IK)     * ZRANDNOSS(IL,IK)
 918         ENDDO
 919       ELSEIF (IK==1) then
 920         DO IL=its,ite
 921         if ((IL.le.(LMAX/2+1))) then
 922           SPFORCC(IL,IK)       = (1.-ALPH)*SPFORCC(IL,IK)      + SP_AMP(IL,IK)     * ZRANDNOSC(IL,IK)
 923           SPFORCS(IL,IK)       = (1.-ALPH)*SPFORCS(IL,IK)      + SP_AMP(IL,IK)     * ZRANDNOSS(IL,IK)
 924         elseif ((IL.gt.(LMAX/2+1))) then
 925           SPFORCC(IL,IK)       = (1.-ALPH)*SPFORCC(IL,IK)      + SP_AMP(IL,IK)     * ZRANDNOSC(LMAX-IL+2,IK)
 926           SPFORCS(IL,IK)       = (1.-ALPH)*SPFORCS(IL,IK)      - SP_AMP(IL,IK)     * ZRANDNOSS(LMAX-IL+2,IK)
 927         endif
 928         ENDDO
 929       ENDIF
 930       ENDDO
 931
 932       DO IK=jts,jte
 933       if (IK.gt.(KMAX/2+1)) then ! Lower half
 934         DO IL=its,ite
 935           if (IL.le.(LMAX/2+1).and.(IL.gt.1)) then !lower left
 936            SPFORCC(IL,IK)      = (1.-ALPH)* SPFORCC(IL,IK)      + SP_AMP(IL,IK)      * ZRANDNOSC(LMAX-IL+2,KMAX-IK+2)
 937            SPFORCS(IL,IK)      = (1.-ALPH)* SPFORCS(IL,IK)      - SP_AMP(IL,IK)      * ZRANDNOSS(LMAX-IL+2,KMAX-IK+2)
 938           elseif (IL.eq.1) then !don't exceed index
 939            SPFORCC(IL,IK)      = (1.-ALPH)* SPFORCC(IL,IK)      + SP_AMP(IL,IK)      * ZRANDNOSC(        1,KMAX-IK+2)
 940            SPFORCS(IL,IK)      = (1.-ALPH)* SPFORCS(IL,IK)      - SP_AMP(IL,IK)      * ZRANDNOSS(        1,KMAX-IK+2)
 941           elseif (IL.gt.(LMAX/2+1)) then !lower right
 942            SPFORCC(IL,IK)      = (1.-ALPH)* SPFORCC(IL,IK)      + SP_AMP(IL,IK)      * ZRANDNOSC(LMAX-IL+2,KMAX-IK+2)
 943            SPFORCS(IL,IK)      = (1.-ALPH)* SPFORCS(IL,IK)      - SP_AMP(IL,IK)      * ZRANDNOSS(LMAX-IL+2,KMAX-IK+2)
 944           endif
 945         ENDDO
 946       endif
 947       ENDDO
 948
 949       call random_seed(get=iseed(1:how_many))
 950       iseedarr=0.0
 951       iseedarr(1:how_many)=iseed
 952
 953      END subroutine UPDATE_STOCH
 954 !     ------------------------------------------------------------------
 955       SUBROUTINE UPDATE_STOCH_TEN(ru_tendf,rv_tendf,t_tendf, &
 956                        ru_tendf_stoch,rv_tendf_stoch,rt_tendf_stoch,&
 957                        mu,mub,c1h,c2h,                              &
 958                        ids, ide, jds, jde, kds, kde,                &
 959                        ims, ime, jms, jme, kms, kme,                &
 960                        its, ite, jts, jte, kts, kte,                &
 961                        kte_stoch,kme_stoch                          )
 962
 963        IMPLICIT NONE
 964        INTEGER , INTENT(IN)        ::  ids, ide, jds, jde, kds, kde,   &
 965                                        ims, ime, jms, jme, kms, kme,   &
 966                                        its, ite, jts, jte, kts, kte,   &
 967                                        kte_stoch,kme_stoch
 968
 969        REAL , DIMENSION(ims:ime , kms:kme, jms:jme),INTENT(INOUT) :: &
 970                                        ru_tendf, rv_tendf, t_tendf
 971
 972        REAL , DIMENSION(ims:ime , kms:kme_stoch, jms:jme)           :: &
 973                       ru_tendf_stoch,rv_tendf_stoch,rt_tendf_stoch
 974
 975        REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: mu,mub
 976        REAL , DIMENSION(kms:kme) , INTENT(IN) :: c1h,c2h
 977
 978        INTEGER :: I,J,K,kh
 979        REAL  :: dt,xm
 980
 981
 982        DO j = jts,MIN(jde-1,jte)
 983          DO k = kts,kte-1
 984            kh=min(k,kte_stoch)
 985            DO i = its,ite
 986              ru_tendf(i,k,j) = ru_tendf(i,k,j) +  ru_tendf_stoch(i,kh,j)  * ((c1h(k)*mu(i,j))+(c1h(k)*mub(i,j)+c2h(k)))
 987            ENDDO
 988          ENDDO
 989        ENDDO
 990
 991        DO j = jts,jte
 992          DO k = kts,kte-1
 993            kh=min(k,kte_stoch)
 994            DO i = its,MIN(ide-1,ite)
 995              rv_tendf(i,k,j) = rv_tendf(i,k,j) +  rv_tendf_stoch(i,kh,j) *  ((c1h(k)*mu(i,j))+(c1h(k)*mub(i,j)+c2h(k)))
 996            ENDDO
 997          ENDDO
 998        ENDDO
 999
1000        DO j = jts,MIN(jde-1,jte)
1001          DO k = kts,kte-1
1002            kh=min(k,kte_stoch)
1003            DO i = its,MIN(ide-1,ite)
1004              t_tendf(i,k,j) = t_tendf(i,k,j) + rt_tendf_stoch(i,kh,j) * ((c1h(k)*mu(i,j))+(c1h(k)*mub(i,j)+c2h(k)))
1005            ENDDO
1006          ENDDO
1007        ENDDO
1008
1009        END SUBROUTINE UPDATE_STOCH_TEN
1010 !     ------------------------------------------------------------------
1011 !!************** PERTURB PHYSICS TENDENCIES (except T) FOR SPPT *******************
1012 !     ------------------------------------------------------------------
1013       subroutine perturb_physics_tend(gridpt_stddev_sppt,               &
1014                        sppt_thresh_fact,rstoch,                         &
1015                        ru_tendf,rv_tendf,t_tendf,moist_tend,            &
1016                        ids, ide, jds, jde, kds, kde,                    &
1017                        ims, ime, jms, jme, kms, kme,                    &
1018                        its, ite, jts, jte, kts, kte,                    &
1019                        kte_stoch,kme_stoch                               )
1020
1021 !      This subroutine add stochastic perturbations of the form
1022 !
1023 !                  rx_tendf(i,k,j)      = rx_tendf(i,k,j)*(1.0 + rstoch(i,k,j))
1024 !
1025 !       to the tendencies of  U, V, and Q.
1026 !       Since the temperature perturbations do not include the micro-physics
1027 !       tendencies at this point, the stochastic tendency perturbations to
1028 !       temperature are added in subroutine rk_addtend_dry of module module_em.F
1029
1030        IMPLICIT NONE
1031        INTEGER , INTENT(IN)        ::  ids, ide, jds, jde, kds, kde,   &
1032                                        ims, ime, jms, jme, kms, kme,   &
1033                                        its, ite, jts, jte, kts, kte,   &
1034                                        kte_stoch,kme_stoch
1035
1036        REAL , DIMENSION(ims:ime , kms:kme, jms:jme),INTENT(INOUT) ::   &
1037                                         ru_tendf, rv_tendf, t_tendf,moist_tend
1038        REAL , DIMENSION(ims:ime,kms:kme_stoch, jms:jme),INTENT(INOUT) :: rstoch
1039        REAL :: gridpt_stddev_sppt ,thresh,sppt_thresh_fact
1040
1041        INTEGER :: I,J,K,kh
1042
1043 ! Here the random process at each gridpoint is capped if it exceeds a value thresh
1044
1045        thresh=sppt_thresh_fact*gridpt_stddev_sppt
1046        DO j = jts,jte
1047          DO k = kts,min(kte-1,kte_stoch-1)
1048            DO i = its,ite
1049 !                rstoch(i,k,j)=MAX(MIN(rstoch(i,k,j),thresh),-1.*thresh))
1050              if (rstoch(i,k,j).lt.-thresh) then
1051                  rstoch(i,k,j)=-thresh
1052              endif
1053              if (rstoch(i,k,j).gt.thresh) then
1054                  rstoch(i,k,j)=thresh
1055              endif
1056            ENDDO
1057          ENDDO
1058        ENDDO
1059
1060 ! Perturb the tendencies of u,v,q,t.
1061        DO j = jts,MIN(jde-1,jte)
1062          DO k = kts,kte-1
1063          kh = min( k, kte_stoch-1 )
1064            DO i = its,ite
1065               ru_tendf(i,k,j)      = ru_tendf(i,k,j)*(1.0 + rstoch(i,kh,j))
1066            ENDDO
1067          ENDDO
1068        ENDDO
1069
1070        DO j = jts,jte
1071          DO k = kts,kte-1
1072          kh = min( k, kte_stoch-1 )
1073             DO i = its,MIN(ide-1,ite)
1074               rv_tendf(i,k,j)      = rv_tendf(i,k,j)*(1.0 + rstoch(i,kh,j))
1075            ENDDO
1076          ENDDO
1077        ENDDO
1078
1079        DO j = jts,MIN(jde-1,jte)
1080          DO k = kts,kte-1
1081          kh = min( k, kte_stoch-1 )
1082            DO i = its,MIN(ide-1,ite)
1083               moist_tend(i,k,j)    = moist_tend(i,k,j)*(1.0 + rstoch(i,kh,j))
1084               t_tendf   (i,k,j)    =    t_tendf(i,k,j)*(1.0 + rstoch(i,kh,j))
1085            ENDDO
1086          ENDDO
1087        ENDDO
1088
1089       end subroutine perturb_physics_tend
1090
1091 !     ------------------------------------------------------------------
1092 !!************** UPDATE SPECTRAL PATTERN AND TRANFORM GRIDPOINT SPACE***
1093 !     ------------------------------------------------------------------
1094 ! This subroutine evolves the spectral pattern and transforms it back to gridpoint space.
1095
1096
1097 !pajm
1098       SUBROUTINE RAND_PERT_UPDATE       (grid, variable_in,                   &
1099                           SPFORCS,SPFORCC,SP_AMP,ALPH_RAND,                   &
1100                           ips, ipe, jps, jpe, kps, kpe,                       &
1101                           ids, ide, jds, jde, kds, kde,                       &
1102                           ims, ime, jms, jme, kms, kme,                       &
1103                           kts, kte,                                           &
1104                           imsx,imex,jmsx,jmex,kmsx,kmex,                      &
1105                           ipsx,ipex,jpsx,jpex,kpsx,kpex,                      &
1106                           imsy,imey,jmsy,jmey,kmsy,kmey,                      &
1107                           ipsy,ipey,jpsy,jpey,kpsy,kpey,                      &
1108                           kpe_stoch,kde_stoch,kme_stoch,kte_stoch,            &
1109                           restart,iseedarr,seed_dim,                          &
1110                           DX,DY,skebs_vertstruc,                              &
1111                           RAND_PERT,thresh_fact,gridpt_stddev,                &
1112                           VERTSTRUCC,VERTSTRUCS,VERTAMP                       )
1113
1114
1115
1116     USE module_domain, ONLY : domain
1117 !pajm Do we need the following line USE module_state_description, ONLY : num_pert3d
1118 #ifdef DM_PARALLEL
1119     USE module_dm, ONLY : local_communicator, mytask, ntasks, ntasks_x, ntasks_y, local_communicator_periodic, &
1120                           wrf_dm_maxval, wrf_err_message, local_communicator_x, local_communicator_y, data_order_xzy
1121 #endif
1122
1123
1124       IMPLICIT NONE
1125
1126       TYPE ( domain ), INTENT(INOUT) :: grid
1127
1128
1129       INTEGER , INTENT(IN)     ::               ids, ide, jds, jde, kds, kde,   &
1130                                                 ims, ime, jms, jme, kms, kme,   &
1131                                                 ips, ipe, jps, jpe, kps, kpe,   &
1132                                                 kts, kte
1133       INTEGER , INTENT(IN)     ::               imsx,imex,jmsx,jmex,kmsx,kmex, &
1134                                                 ipsx,ipex,jpsx,jpex,kpsx,kpex, &
1135                                                 imsy,imey,jmsy,jmey,kmsy,kmey, &
1136                                                 ipsy,ipey,jpsy,jpey,kpsy,kpey
1137       INTEGER , INTENT(IN)     ::               seed_dim
1138       INTEGER                  ::               kpe_stoch,kde_stoch,kme_stoch,kte_stoch
1139
1140       REAL    , INTENT(IN)     ::               ALPH_RAND,dx,dy,thresh_fact,gridpt_stddev
1141       INTEGER , INTENT(IN)     ::               skebs_vertstruc
1142       CHARACTER, INTENT(IN)    ::               variable_in ! T, U, V
1143                                                ! T          ! random field, T
1144                                                ! U          ! first derivative of streamfunction with regard to y; for skebs: U
1145                                                ! V          ! first derivative of streamfunction with regard to x; for skebs: V
1146
1147       INTEGER, DIMENSION (seed_dim),            INTENT(INOUT) :: iseedarr
1148       REAL, DIMENSION(ims:ime,kms:kme, jms:jme),INTENT(IN)    :: VERTSTRUCC,VERTSTRUCS
1149       REAL, DIMENSION(ims:ime,jms:jme)         ,INTENT(INOUT) :: SPFORCS,SPFORCC,SP_AMP
1150       REAL, DIMENSION(kms:kme )                ,INTENT(IN)    :: VERTAMP
1151       REAL, DIMENSION(ims:ime,kms:kme_stoch, jms:jme)         :: RAND_PERT
1152       REAL                     ::               RY,RX
1153
1154
1155 ! Local Variabels
1156       INTEGER                  ::               IK,IL,ILEV,NLON,NLAT,IJ,I,J,K
1157       INTEGER                  ::               gridsp32y,gridsm32y,gridsp32x,gridsm32x,gridsp32 ,gridsm32
1158       INTEGER                  ::               gridep32y,gridem32y,gridep32x,gridem32x,gridep32 ,gridem32
1159       REAL                     ::               thresh
1160
1161       REAL, DIMENSION(ims:ime,kms:kme_stoch, jms:jme)         :: RAND_REAL, RAND_IMAG
1162       LOGICAL :: RESTART
1163       CHARACTER  :: variable
1164
1165       variable = variable_in
1166
1167       NLAT=(jde-jds)+1 !KMAX
1168       NLON=(ide-ids)+1 !LMAX
1169       RY=  NLAT*DY
1170       RX=  NLON*DX
1171
1172 ! Update the pattern generator by evolving each spectral coefficients as AR1
1173
1174       !$OMP PARALLEL DO   &
1175       !$OMP PRIVATE ( ij )
1176        DO ij = 1 , grid%num_tiles
1177              IF (variable .ne. 'V') THEN  !T, random field, U, don't update for V
1178               CALL UPDATE_STOCH( &
1179                          SPFORCS,SPFORCC,SP_AMP,ALPH_RAND,                   &
1180                          restart,iseedarr,seed_dim,                          &
1181                          ids, ide, jds, jde, kds, kde,                       &
1182                          ims, ime, jms, jme, kms, kme,                       &
1183                          grid%i_start(ij), grid%i_end(ij), grid%j_start(ij), grid%j_end(ij), kts, kte                        )
1184              endif
1185
1186 ! Put spectral coefficients in arrays RAND_REAL,RAND_IMAG
1187
1188               IF (variable == 'T')   THEN  ! T, rand
1189               DO IK=grid%j_start(ij), grid%j_end(ij)
1190                DO ILEV=kts,kte_stoch
1191                 DO IL=grid%i_start(ij),grid%i_end(ij)
1192                        grid%RAND_REAL(IL,ILEV,IK)  = SPFORCC(IL,IK)
1193                        grid%RAND_IMAG(IL,ILEV,IK)  = SPFORCS(IL,IK)
1194                   ENDDO
1195                 ENDDO
1196               ENDDO
1197
1198               ELSEIF (variable == 'U') THEN !U
1199               DO IK=grid%j_start(ij), grid%j_end(ij)
1200                DO ILEV=kts,kte_stoch
1201                 DO IL=grid%i_start(ij),grid%i_end(ij)
1202                        grid%RAND_REAL(IL,ILEV,IK)  =  2*RPI/RY*  wavenumber_k(IK) * SPFORCS(IL,IK)
1203                        grid%RAND_IMAG(IL,ILEV,IK)  = -2*RPI/RY*  wavenumber_k(IK) * SPFORCC(IL,IK)
1204                   ENDDO
1205                 ENDDO
1206               ENDDO
1207
1208               ELSEIF (variable == 'V') THEN !V
1209               DO IK=grid%j_start(ij), grid%j_end(ij)
1210                DO ILEV=kts,kte_stoch
1211                 DO IL=grid%i_start(ij),grid%i_end(ij)
1212                        grid%RAND_REAL(IL,ILEV,IK)  = -2*RPI/RX*  wavenumber_l(IL) * SPFORCS(IL,IK)
1213                        grid%RAND_IMAG(IL,ILEV,IK)  =  2*RPI/RX*  wavenumber_l(IL) * SPFORCC(IL,IK)
1214                   ENDDO
1215                 ENDDO
1216               ENDDO
1217               endif
1218
1219
1220 ! Apply vertical structure function
1221
1222            IF (skebs_vertstruc.ne.0) then
1223              DO ILEV=kts,kte_stoch
1224               DO IL=grid%i_start(ij),grid%i_end(ij)
1225                 DO IK=grid%j_start(ij), grid%j_end(ij)
1226                   grid%RAND_REAL(IL,ILEV,IK)  = VERTAMP(ILEV) * &
1227                        (grid%RAND_REAL(IL,ILEV,IK) * VERTSTRUCC(IL,ILEV,IK) - grid%RAND_IMAG(IL,ILEV,IK) * VERTSTRUCS(IL,ILEV,IK))
1228                   grid%RAND_IMAG(IL,ILEV,IK)  = VERTAMP(ILEV) * &
1229                        (grid%RAND_REAL(IL,ILEV,IK) * VERTSTRUCS(IL,ILEV,IK) + grid%RAND_IMAG(IL,ILEV,IK) * VERTSTRUCC(IL,ILEV,IK))
1230                 ENDDO
1231              ENDDO
1232            ENDDO
1233            ENDIF
1234         ENDDO
1235        !$OMP END PARALLEL DO
1236
1237 ! Transform spectral pattern to gridpoint space
1238
1239 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
1240
1241 ! Roll out into latitude bands and perform FFT along latitude bands
1242
1243 ! Save a copy of the indices as we might need them to change when
1244 ! doing the "thin" 3d arrays (where the "k" dimension is unity).
1245 ! These are the original Z-transposed and X-transposed k-dimensions.
1246
1247         gridsp32x=grid%sp32x
1248         gridsm32x=grid%sm32x
1249         gridep32x=grid%ep32x
1250         gridem32x=grid%em32x
1251         gridsp32 =grid%sp32
1252         gridsm32 =grid%sm32
1253         gridep32 =grid%ep32
1254         gridem32 =grid%em32
1255
1256 ! Set number of vertical levels to which ever is smaller: the full number
1257 ! of vertical levels, or the number of levels to be transformed into
1258 ! gridpoint space.
1259
1260         grid%sp32x=min(kpsx,grid%num_stoch_levels)
1261         grid%sm32x=min(kmsx,grid%num_stoch_levels)
1262         grid%ep32x=min(kpex,grid%num_stoch_levels)
1263         grid%em32x=min(kmex,grid%num_stoch_levels)
1264         grid%sp32 =min(kps ,grid%num_stoch_levels)
1265         grid%sm32 =min(kms ,grid%num_stoch_levels)
1266         grid%ep32 =min(kpe ,grid%num_stoch_levels)
1267         grid%em32 =min(kme ,grid%num_stoch_levels)
1268
1269 #include "XPOSE_RAND_REAL_z2x.inc"
1270 #include "XPOSE_RAND_IMAG_z2x.inc"
1271         call do_fftback_along_x(grid%RAND_REAL_xxx,grid%RAND_IMAG_xxx,                  &
1272                               ids,ide,jds,jde,                                          &
1273                               imsx,imex,jmsx,jmex,kmsx,min(kmex,grid%num_stoch_levels), &
1274                               ipsx,ipex,jpsx,jpex,kpsx,min(kpex,grid%num_stoch_levels))
1275 #include "XPOSE_RAND_REAL_x2z.inc"
1276 #include "XPOSE_RAND_IMAG_x2z.inc"
1277
1278 ! Roll out into longitude bands and perform FFT along longitude bands
1279
1280 ! Save a copy of the indices as we might need them to change when
1281 ! doing the "thin" 3d arrays (where the "k" dimension is unity).
1282 ! These are the original Y-transposed k-dimensions.
1283
1284         gridsp32y=grid%sp32y
1285         gridsm32y=grid%sm32y
1286         gridep32y=grid%ep32y
1287         gridem32y=grid%em32y
1288
1289 ! Again, set number of vertical levels to the min of the number of levels and the
1290 ! number of stochastic levels.
1291
1292         grid%sp32y=min(kpsy,grid%num_stoch_levels)
1293         grid%sm32y=min(kmsy,grid%num_stoch_levels)
1294         grid%ep32y=min(kpey,grid%num_stoch_levels)
1295         grid%em32y=min(kmey,grid%num_stoch_levels)
1296
1297 #include "XPOSE_RAND_REAL_z2y.inc"
1298 #include "XPOSE_RAND_IMAG_z2y.inc"
1299         call do_fftback_along_y(grid%RAND_REAL_yyy,grid%RAND_IMAG_yyy,                  &
1300                               ids,ide,jds,jde,                                          &
1301                               imsy,imey,jmsy,jmey,kmsy,min(kmey,grid%num_stoch_levels), &
1302                               ipsy,ipey,jpsy,jpey,kpsy,min(kpey,grid%num_stoch_levels))
1303 #include "XPOSE_RAND_REAL_y2z.inc"
1304 #include "XPOSE_RAND_IMAG_y2z.inc"
1305
1306 ! Put the original vertical "k" dimensions back.
1307
1308         grid%sp32x=gridsp32x
1309         grid%sm32x=gridsm32x
1310         grid%ep32x=gridep32x
1311         grid%em32x=gridem32x
1312         grid%sp32y=gridsp32y
1313         grid%sm32y=gridsm32y
1314         grid%ep32y=gridep32y
1315         grid%em32y=gridem32y
1316         grid%sp32 =gridsp32
1317         grid%sm32 =gridsm32
1318         grid%ep32 =gridep32
1319         grid%em32 =gridem32
1320
1321 #else
1322         call do_fftback_along_x(grid%RAND_REAL,grid%RAND_IMAG,                          &
1323                               ids,ide,jds,jde,                                          &
1324                               ims,ime,jms,jme,kms,min(kme,grid%num_stoch_levels),       &
1325                               ips,ipe,jps,jpe,kps,min(kpe,grid%num_stoch_levels))
1326         call do_fftback_along_y(grid%RAND_REAL,grid%RAND_IMAG,                          &
1327                               ids,ide,jds,jde,                                          &
1328                               ims,ime,jms,jme,kms,min(kme,grid%num_stoch_levels),       &
1329                               ips,ipe,jps,jpe,kps,min(kpe,grid%num_stoch_levels))
1330 #endif
1331
1332
1333       thresh=thresh_fact*gridpt_stddev
1334       !$OMP PARALLEL DO   &
1335       !$OMP PRIVATE ( ij )
1336         DO ij = 1 , grid%num_tiles
1337                DO k=kts,min(kte,grid%num_stoch_levels)
1338                  DO I=grid%i_start(ij), grid%i_end(ij)
1339                    DO j=grid%j_start(ij), grid%j_end(ij)
1340                      RAND_PERT(I,K,J)=grid%RAND_REAL(I,K,J)
1341                      RAND_PERT(I,K,J)=MAX(MIN(grid%RAND_REAL(I,K,J),thresh),-1.0*thresh)
1342                    ENDDO
1343                  ENDDO
1344                 ENDDO
1345         ENDDO
1346        !$OMP END PARALLEL DO
1347
1348
1349        END SUBROUTINE RAND_PERT_UPDATE
1350
1351 !     ------------------------------------------------------------------
1352 !!************** SUBROUTINE DO_FFTBACK_ALONG_X
1353 !     ------------------------------------------------------------------
1354        subroutine do_fftback_along_x(                            &
1355                                fieldc,fields,                    &
1356                                ids,ide,jds,jde,                  &
1357                                imsx,imex,jmsx,jmex,kmsx,kmex,    &
1358                                ipsx,ipex,jpsx,jpex,kpsx,kpex     )
1359        IMPLICIT NONE
1360
1361        INTEGER, INTENT(IN):: imsx,imex,jmsx,jmex,kmsx,kmex,    &
1362                              ipsx,ipex,jpsx,jpex,kpsx,kpex,    &
1363                              ids,ide,jds,jde
1364
1365
1366        REAL, DIMENSION    (imsx:imex, kmsx:kmex, jmsx:jmex) :: fieldc,fields
1367
1368        COMPLEX, DIMENSION (ipsx:ipex)            :: dummy_complex
1369        INTEGER                                   :: IER,LENWRK,KMAX,LMAX,I,J,K
1370        REAL, ALLOCATABLE                         :: WORK(:)
1371
1372        CHARACTER (LEN=160) :: mess
1373
1374
1375        KMAX=(jde-jds)+1
1376        LMAX=(ide-ids)+1
1377
1378        LENWRK=2*KMAX*LMAX
1379        ALLOCATE(WORK(LENWRK))
1380        LENSAV= 4*(KMAX+LMAX)+INT(LOG(REAL(KMAX))) + INT(LOG(REAL(LMAX))) + 8
1381
1382        DO k=kpsx,kpex
1383          DO j = jpsx, jpex
1384            DO i = ipsx, ipex
1385              dummy_complex(i)=cmplx(fieldc(i,k,j),fields(i,k,j))
1386            ENDDO
1387            CALL cFFT1B (LMAX, 1 ,dummy_complex,LMAX, WSAVE1, LENSAV, WORK, LENWRK, IER)
1388            if (ier.ne.0) then
1389               WRITE(mess,FMT='(A)') 'error in cFFT1B in do_fftback_along_x, field U'
1390               CALL wrf_debug(0,mess)
1391            end if
1392            DO i = ipsx, ipex
1393              fieldc(i,k,j)=real(dummy_complex(i))
1394              fields(i,k,j)=imag(dummy_complex(i))
1395            END DO
1396          END DO
1397        END DO
1398
1399        DEALLOCATE(WORK)
1400        end subroutine do_fftback_along_x
1401
1402 !!     ------------------------------------------------------------------
1403 !!!************** SUBROUTINE DO_FFTBACK_ALONG_Y
1404 !!     ------------------------------------------------------------------
1405        subroutine do_fftback_along_y(                            &
1406                                fieldc,fields,                    &
1407                                ids,ide,jds,jde,                  &
1408                                imsy,imey,jmsy,jmey,kmsy,kmey,    &
1409                                ipsy,ipey,jpsy,jpey,kpsy,kpey     )
1410        IMPLICIT NONE
1411
1412        INTEGER :: IER,LENWRK,KMAX,LMAX,I,J,K,skebs_vertstruc
1413
1414        INTEGER, INTENT(IN) :: imsy,imey,jmsy,jmey,kmsy,kmey,    &
1415                               ipsy,ipey,jpsy,jpey,kpsy,kpey,    &
1416                               ids,ide,jds,jde
1417
1418        REAL, DIMENSION    (imsy:imey, kmsy:kmey, jmsy:jmey) :: fieldc,fields
1419
1420        COMPLEX, DIMENSION (jpsy:jpey)            :: dummy_complex
1421        REAL, ALLOCATABLE :: WORK(:)
1422
1423        CHARACTER (LEN=160) :: mess
1424
1425        KMAX=(jde-jds)+1
1426        LMAX=(ide-ids)+1
1427        LENWRK=2*KMAX*LMAX
1428        ALLOCATE(WORK(LENWRK))
1429        LENSAV= 4*(KMAX+LMAX)+INT(LOG(REAL(KMAX))) + INT(LOG(REAL(LMAX))) + 8
1430
1431
1432
1433         DO k=kpsy,kpey
1434           DO i = ipsy, ipey
1435             DO j = jpsy,jpey
1436             dummy_complex(j)=cmplx(fieldc(i,k,j),fields(i,k,j))
1437             ENDDO
1438             CALL cFFT1B (KMAX, 1 ,dummy_complex,KMAX, WSAVE2, LENSAV, WORK, LENWRK, IER)
1439             if (ier.ne.0) then
1440                WRITE(mess,FMT='(A)') 'error in cFFT1B in do_fftback_along_y, field U'
1441                CALL wrf_debug(0,mess)
1442             end if
1443             DO j = jpsy, jpey
1444             fieldc(i,k,j)=real(dummy_complex(j))
1445             fields(i,k,j)=imag(dummy_complex(j))
1446             END DO
1447           END DO
1448         END DO ! k_start-k_end
1449
1450
1451        DEALLOCATE(WORK)
1452        end subroutine do_fftback_along_y
1453 !     ------------------------------------------------------------------
1454 !!************** TRANSFORM FROM GRIDPOILT SPACE TO SPHERICAL HARMONICS **
1455 !     ------------------------------------------------------------------
1456       subroutine findindex( wavenumber_k, wavenumber_L,             &
1457                        ids, ide, jds, jde, kds, kde,                &
1458                        ims, ime, jms, jme, kms, kme,                &
1459                        its, ite, jts, jte, kts, kte                 )
1460
1461       IMPLICIT NONE
1462       INTEGER :: IK,IL,KMAX,LMAX
1463       INTEGER, DIMENSION (jds:jde)::  wavenumber_k
1464       INTEGER, DIMENSION (ids:ide)::  wavenumber_l
1465       INTEGER , INTENT(IN)     ::  ids, ide, jds, jde, kds, kde,   &
1466                                    ims, ime, jms, jme, kms, kme,   &
1467                                    its, ite, jts, jte, kts, kte
1468       KMAX=(jde-jds)+1
1469       LMAX=(ide-ids)+1
1470
1471       !map wave numbers K,L to indeces IK, IL
1472       DO IK=1,KMAX/2+1
1473         wavenumber_k(IK)=IK-1
1474       ENDDO
1475       DO IK=KMAX,KMAX/2+2,-1
1476         wavenumber_k(IK)=IK-KMAX-1
1477       ENDDO
1478       DO IL=1,LMAX/2+1
1479         wavenumber_l(IL)=IL-1
1480       ENDDO
1481       DO IL=LMAX,LMAX/2+2,-1
1482         wavenumber_l(IL)=IL-LMAX-1
1483       ENDDO
1484
1485       END subroutine findindex
1486
1487 !     ------------------------------------------------------------------
1488      subroutine gauss_noise(z)
1489       real :: z                    ! output
1490       real :: x,y,r, coeff         ! INPUT
1491
1492 !  [2.1] Get two uniform variate random numbers IL range 0 to 1:
1493
1494       do
1495
1496       call random_number( x )
1497       call random_number( y )
1498
1499 !     [2.2] Transform to range -1 to 1 and calculate sum of squares:
1500
1501       x = 2.0 * x - 1.0
1502       y = 2.0 * y - 1.0
1503       r = x * x + y * y
1504
1505       if ( r > 0.0 .and. r < 1.0 ) exit
1506
1507       end do
1508 !
1509 !  [2.3] Use Box-Muller transformation to get normal deviates:
1510
1511       coeff = sqrt( -2.0 * log(r) / r )
1512       z = coeff * x
1513
1514      end subroutine gauss_noise
1515 !     ------------------------------------------------------------------
1516      SUBROUTINE rand_seed (config_flags, iseed1, iseedarr, seed_start, seed_dim )
1517      USE module_configure
1518      IMPLICIT NONE
1519 !
1520 !  Structure that contains run-time configuration (namelist) data for domain
1521       TYPE (grid_config_rec_type)                       :: config_flags
1522 !
1523 ! Arguments
1524      INTEGER             :: seed_start, seed_dim, iseed1
1525      INTEGER, DIMENSION (seed_start:seed_dim), INTENT(OUT)           :: iseedarr
1526
1527 ! Local
1528       integer*8          :: fctime, one_big
1529       integer            :: i
1530
1531       fctime = config_flags%start_year * ( config_flags%start_month*100+config_flags%start_day) + config_flags%start_hour
1532
1533       one_big = 1
1534       iseedarr=0.0
1535       if ( seed_dim-3 .lt. seed_start ) then
1536          do i = seed_start,seed_dim
1537             iseedarr(i  )= iseed1+config_flags%nens*1000000
1538          enddo
1539       else
1540          do i = seed_start,seed_dim-3,4
1541             iseedarr(i  )= iseed1+config_flags%nens*1000000
1542             iseedarr(i+1)= mod(fctime+iseed1*config_flags%nens*1000000,19211*one_big)
1543             iseedarr(i+2)= mod(fctime+iseed1*config_flags%nens*1000000,71209*one_big)
1544             iseedarr(i+3)= mod(fctime+iseed1*config_flags%nens*1000000,11279*one_big)
1545          enddo
1546       end if
1547
1548       end SUBROUTINE rand_seed
1549
1550 !     ------------------------------------------------------------------
1551
1552       SUBROUTINE contiguize_2d ( in_out, d3, d2, n ,               &
1553                                  ips, ipe, jps, jpe, kps, kpe,     &
1554                                  ids, ide, jds, jde, kds, kde,     &
1555                                  ims, ime, jms, jme, kms, kme      )
1556          IMPLICIT NONE
1557
1558          LOGICAL, INTENT(IN) :: in_out ! d3->d2 = T; d2->d3 = F
1559          INTEGER, INTENT(IN) :: n
1560          INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde,     &
1561                                 ims, ime, jms, jme, kms, kme,     &
1562                                 ips, ipe, jps, jpe, kps, kpe
1563
1564          REAL, INTENT(INOUT), DIMENSION(ims:ime,        jms:jme) :: d2
1565          REAL, INTENT(INOUT), DIMENSION(ims:ime,kms:kme,jms:jme) :: d3
1566
1567          INTEGER :: i, j
1568
1569          IF ( in_out ) THEN
1570             DO j = jms, jme
1571                DO i = ims, ime
1572                   d2(i,j) = d3(i,n,j)
1573                END DO
1574             END DO
1575          ELSE
1576             DO j = jms, jme
1577                DO i = ims, ime
1578                   d3(i,n,j) = d2(i,j)
1579                END DO
1580             END DO
1581          END IF
1582
1583       END SUBROUTINE contiguize_2d
1584
1585 !     ------------------------------------------------------------------
1586       end module module_stoch