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