#if ( HYBRID_COORD==1 ) # define mu(...) (c1h(k)*XXPCXX(__VA_ARGS__)) # define XXPCXX(...) mu(__VA_ARGS__) # define mub(...) (c1h(k)*XXPCBXX(__VA_ARGS__)+c2h(k)) # define XXPCBXX(...) mub(__VA_ARGS__) #endif module module_stoch !*********************************************************************** ! ! Purpose: Stochastic Perturbation Schemes ! Author : Judith Berner, NCAR (berner@ucar.edu) ! Date : Apr 2014 ! !*********************************************************************** ! ! The scheme introduces stochastic perturbations to the rotational wind ! components and to the potential temperature field. The stochastic ! perturbations are generated by independent autoregressive processes for ! each wavenumber and results in smooth spatially and temporally correlated patterns. ! Details of the scheme and its performance in a meso-scale WRF-ensemble ! system are available in: ! ! Berner, J., S.-Y. Ha, J. P. Hacker, A. Fournier and C. Snyder 2011: ! "Model uncertainty in a mesoscale ensemble prediction system: Stochastic ! versus multi-physics representations", 2011, Mon. Wea. Rev., 139, 1972-1995 ! http://journals.ametsoc.org/doi/abs/10.1175/2010MWR3595.1 ! ! Features: ! Dissipation: Dissipation rates are assumed constant in space and time ! Vertical structure: Supports two options for vertical structure: ! 0) constant ! 1) random phase ! ! Optional namelist parameters: ! stoch_force_opt = 0, 0, 0: No stochastic parameterization ! = 1, 1, 1: Use SKEB scheme ! skebs_vertstruc = 0, 0, 0: Constant vertical structure of random pattern generator ! = 1, 1, 1: Random phase vertical structure random pattern generator ! tot_backscat_psi : Total backscattered dissipation rate for streamfunction; Controls ! amplitude of rotational wind perturbations Default value is 1.0E-5 m2/s3. ! tot_backscat_t : Total backscattered dissipation rate for potential temperature; ! Controls amplitude of potential temperature perturbations. Default value is 1.0E-6 m2/s3. ! nens : Random seed for random number stream. This parameter needs to be different ! for each member in ensemble forecasts. Is a function of initial start time ! to ensure different random number streams for different forecasts. ! ztau_psi : Decorrelation time (in s) for streamfunction perturbations. ! Default is 10800s. Recommended value is 216000s. ! ztau_t : Decorrelation time (in s) for potential temperature perturbations. ! Default 10800s. Recommended value is 216000s. ! rexponent_psi : Spectral slope for streamfunction perturbations. Default is -1.83 ! for a kinetic-energy forcing spectrum with slope -5/3. ! rexponent_t : Spectral slope of potential temperature perturbations. Default is -1.83 ! for a potential energy forcing spectrum with slope -1.832. ! kminforc : Minimal forcing wavenumber in longitude for streamfunction perturbations. Default is 1. ! lminforc : Minimal forcing wavenumber in latitude for streamfunction perturbations. Default is 1. ! kminforc : Minimal forcing wavenumber in longitude for potential temperature perturbations. Default is 1. ! lminforct : Minimal forcing wavenumber in latitude for potential temperature perturbations. Default is 1. ! kmaxforc : Maximal forcing wavenumber in longitude for streamfunction perturbations. ! Default is maximal possible wavenumbers determined by number of gridpoints. ! lmaxforc : Maximal forcing wavenumber in latitude for streamfunction perturbations. ! Default is maximal possible wavenumbers determined by number of gridpoints. ! kmaxforct : Maximal forcing wavenumber in longitude for potential temperature perturbations. ! Default is maximal possible wavenumbers determined by number of gridpoints. ! lmaxforct : Maximal forcing wavenumber in latitude for potential temperature perturbations. ! Default is maximal possible wavenumbers determined by number of gridpoints. ! zsigma2_eps : Noise variance in autoregressive process defining streamfunction perturbations. ! zsigma2_eta : Noise variance in autoregressive process defining in potential temperature perturbations. !*********************************************************************** ! ------------------------------------------------------------------ !************** DECLARE FIELDS AND VARIABLES FOR STOCHASTIC BACKSCATTER ! ------------------------------------------------------------------ implicit none public :: SETUP_RAND_PERTURB, UPDATE_STOCH,& rand_pert_update INTEGER :: LMINFORC, LMAXFORC, KMINFORC, KMAXFORC, & & LMINFORCT, LMAXFORCT, KMINFORCT, KMAXFORCT REAL :: ALPH, ALPH_PSI, ALPH_T, TOT_BACKSCAT_PSI, TOT_BACKSCAT_T, REXPONENT_PSI,REXPONENT_T ! ----------Fields for spectral transform ----------- INTEGER :: LENSAV INTEGER,ALLOCATABLE:: wavenumber_k(:), wavenumber_l(:) REAL, ALLOCATABLE :: WSAVE1(:),WSAVE2(:), WSAVE1_ideal(:),WSAVE2_ideal(:) REAL, ALLOCATABLE :: rindarrayik(:),rindarrayil(:),rindarrayl(:),rindarrayk(:) ! --------- Others ------------------------------------------------- REAL, PARAMETER:: RPI= 3.141592653589793 !4.0*atan(1.0) REAL, PARAMETER:: CP= 1006.0 ! specific heat of dry air in J/(Kg*K)= m^2/(K* s^2) REAL, PARAMETER:: T0= 300.0 ! Reference temperature in K character(len=33) :: filenamesave save !======================================================================= contains !======================================================================= ! ------------------------------------------------------------------ !!************** INITIALIZE STOCHASTIC ROUTINES ***************************** ! ------------------------------------------------------------------ ! This subroutine drives the initialization of the stochastic schemes SUBROUTINE INITIALIZE_STOCH (grid, config_flags, & first_trip_for_this_domain, & ips, ipe, jps, jpe, kps, kpe, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte, & imsx, imex, jmsx, jmex, kmsx, kmex, & ipsx, ipex, jpsx, jpex, kpsx, kpex, & imsy, imey, jmsy, jmey, kmsy, kmey, & ipsy, ipey, jpsy, jpey, kpsy, kpey ) USE module_configure USE module_domain, ONLY : domain IMPLICIT NONE TYPE (grid_config_rec_type) :: config_flags TYPE ( domain ), INTENT(INOUT) :: grid INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe, & its, ite, jts, jte, kts, kte INTEGER , INTENT(IN) :: imsx,imex,jmsx,jmex,kmsx,kmex, & ipsx,ipex,jpsx,jpex,kpsx,kpex, & imsy,imey,jmsy,jmey,kmsy,kmey, & ipsy,ipey,jpsy,jpey,kpsy,kpey LOGICAL :: first_trip_for_this_domain INTEGER, PARAMETER :: maxspinup=100 INTEGER :: K IF ( first_trip_for_this_domain ) THEN grid%did_stoch = .FALSE. END IF IF ((( grid%id == 1) .AND. (.NOT. grid%did_stoch)) .AND. & (( grid%skebs_on== 1) .OR.( grid%sppt_on== 1) .OR. ( grid%rand_perturb_on== 1) .OR. & ( grid%spp_conv== 1) .OR. ( grid%spp_pbl== 1) .OR. ( grid%spp_mp .GE. 1) .OR. & ( grid%spp_lsm== 1)) ) THEN grid%did_stoch = .TRUE. IF (grid%skebs_on==1) then ! Initialize SKEBS ! Initialize streamfunction (1) if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then call rand_seed (config_flags, grid%ISEED_SKEBS, grid%iseedarr_skebs , 1, config_flags%seed_dim) else call read_write_stochrestart ('stoch_skebs_U',grid,config_flags,grid%lmax_ideal,grid%kmax_ideal, & grid%SPSTREAMFORCS,grid%SPSTREAMFORCC,2) call read_write_stochrestart ('stoch_skebs_T',grid,config_flags,grid%lmax_ideal,grid%kmax_ideal, & grid%SPTFORCS,grid%SPTFORCC,2) ! broadcast SPFORC and SPFORS to all processes #ifdef DM_PARALLEL CALL wrf_dm_bcast_real ( grid%SPSTREAMFORCS ,grid%kmax_ideal*grid%lmax_ideal ) CALL wrf_dm_bcast_real ( grid%SPSTREAMFORCC ,grid%kmax_ideal*grid%lmax_ideal ) CALL wrf_dm_bcast_real ( grid%SPTFORCS,grid%kmax_ideal*grid%lmax_ideal ) CALL wrf_dm_bcast_real ( grid%SPTFORCC,grid%kmax_ideal*grid%lmax_ideal ) #endif endif call SETUP_RAND_PERTURB('W', & grid%skebs_vertstruc,config_flags%restart, & grid%SPSTREAM_AMP, & grid%SPSTREAMFORCS,grid%SPSTREAMFORCC,grid%ALPH_PSI,& grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPUV, & grid%KMINFORCT,grid%KMAXFORCT, & grid%LMINFORCT,grid%LMAXFORCT, & grid%KMAXFORCTH,grid%LMAXFORCTH, & grid%time_step,grid%DX,grid%DY, & grid%stepstoch, & grid%gridpt_stddev_sppt, & grid%lengthscale_sppt, & grid%timescale_sppt, & grid%TOT_BACKSCAT_PSI,grid%ZTAU_PSI, & grid%REXPONENT_PSI, & grid%kmax_ideal, grid%lmax_ideal, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then do k = 1,maxspinup CALL RAND_PERT_UPDATE(grid,'U', & grid%SPSTREAMFORCS,grid%SPSTREAMFORCC, & grid%SPSTREAM_AMP,grid%ALPH_PSI, & ips, ipe, jps, jpe, kps, kpe, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & kts,kte, & imsx,imex,jmsx,jmex,kmsx,kmex, & ipsx,ipex,jpsx,jpex,kpsx,kpex, & imsy,imey,jmsy,jmey,kmsy,kmey, & ipsy,ipey,jpsy,jpey,kpsy,kpey, & grid%kmax_ideal, grid%lmax_ideal, & grid% num_stoch_levels,grid% num_stoch_levels, & grid% num_stoch_levels,grid% num_stoch_levels, & config_flags%restart, grid%iseedarr_skebs, & config_flags%seed_dim, & grid%DX,grid%DY,grid%skebs_vertstruc, & grid%ru_tendf_stoch, & grid%stddev_cutoff_sppt,grid%gridpt_stddev_sppt, & grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPUV ) CALL RAND_PERT_UPDATE(grid,'V', & grid%SPSTREAMFORCS,grid%SPSTREAMFORCC, & grid%SPSTREAM_AMP,grid%ALPH_PSI, & ips, ipe, jps, jpe, kps, kpe, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & kts,kte, & imsx,imex,jmsx,jmex,kmsx,kmex, & ipsx,ipex,jpsx,jpex,kpsx,kpex, & imsy,imey,jmsy,jmey,kmsy,kmey, & ipsy,ipey,jpsy,jpey,kpsy,kpey, & grid%kmax_ideal, grid%lmax_ideal, & grid% num_stoch_levels,grid% num_stoch_levels, & grid% num_stoch_levels,grid% num_stoch_levels, & config_flags%restart, grid%iseedarr_skebs, & config_flags%seed_dim, & grid%DX,grid%DY,grid%skebs_vertstruc, & grid%rv_tendf_stoch, & grid%stddev_cutoff_sppt,grid%gridpt_stddev_sppt, & grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT ) enddo endif ! Initialize potential temperature (2) call SETUP_RAND_PERTURB('T', & grid%skebs_vertstruc,config_flags%restart, & grid%SPT_AMP, & grid%SPTFORCS,grid%SPTFORCC,grid%ALPH_T, & grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT, & grid%KMINFORCT,grid%KMAXFORCT, & grid%LMINFORCT,grid%LMAXFORCT, & grid%KMAXFORCTH,grid%LMAXFORCTH, & grid%time_step,grid%DX,grid%DY, & grid%stepstoch, & grid%gridpt_stddev_sppt, & grid%lengthscale_sppt, & grid%timescale_sppt, & grid%TOT_BACKSCAT_T,grid%ZTAU_T, & grid%REXPONENT_T, & grid%kmax_ideal, grid%lmax_ideal, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then do k = 1,maxspinup CALL RAND_PERT_UPDATE(grid,'T', & grid%SPTFORCS,grid%SPTFORCC, & grid%SPT_AMP,grid%ALPH_T, & ips, ipe, jps, jpe, kps, kpe, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & kts,kte, & imsx,imex,jmsx,jmex,kmsx,kmex, & ipsx,ipex,jpsx,jpex,kpsx,kpex, & imsy,imey,jmsy,jmey,kmsy,kmey, & ipsy,ipey,jpsy,jpey,kpsy,kpey, & grid%kmax_ideal, grid%lmax_ideal, & grid%num_stoch_levels,grid%num_stoch_levels, & grid%num_stoch_levels,grid%num_stoch_levels, & config_flags%restart, grid%iseedarr_skebs, & config_flags%seed_dim, & grid%DX,grid%DY,grid%skebs_vertstruc, & grid%rt_tendf_stoch, & grid%stddev_cutoff_sppt,grid%gridpt_stddev_sppt, & grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPUV ) enddo endif ENDIF IF (grid%sppt_on==1) then ! Initialize SPPT (3) if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then call rand_seed (config_flags, grid%ISEED_SPPT, grid%iseedarr_sppt , 1, config_flags%seed_dim) else call read_write_stochrestart ('stoch_sppt_xx',grid,config_flags,grid%lmax_ideal,grid%kmax_ideal, & grid%SPPTFORCS,grid%SPPTFORCC,2) ! broadcast SPFORC and SPFORS to all processes #ifdef DM_PARALLEL CALL wrf_dm_bcast_real ( grid%SPPTFORCS,grid%kmax_ideal*grid%lmax_ideal ) CALL wrf_dm_bcast_real ( grid%SPPTFORCC,grid%kmax_ideal*grid%lmax_ideal ) #endif endif call SETUP_RAND_PERTURB('P', & grid%sppt_vertstruc,config_flags%restart, & grid%SPPT_AMP, & grid%SPPTFORCC,grid%SPPTFORCS,grid%ALPH_SPPT, & grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT, & grid%KMINFORCT,grid%KMAXFORCT, & grid%LMINFORCT,grid%LMAXFORCT, & grid%KMAXFORCTH,grid%LMAXFORCTH, & grid%time_step,grid%DX,grid%DY, & grid%stepstoch, & grid%gridpt_stddev_sppt, & grid%lengthscale_sppt, & grid%timescale_sppt, & grid%TOT_BACKSCAT_PSI,grid%ZTAU_PSI, & grid%REXPONENT_PSI, & grid%kmax_ideal, grid%lmax_ideal, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then do k = 1,maxspinup CALL RAND_PERT_UPDATE(grid,'T', & grid%SPPTFORCS,grid%SPPTFORCC, & grid%SPPT_AMP,grid%ALPH_SPPT, & ips, ipe, jps, jpe, kps, kpe, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & kts,kte, & imsx,imex,jmsx,jmex,kmsx,kmex, & ipsx,ipex,jpsx,jpex,kpsx,kpex, & imsy,imey,jmsy,jmey,kmsy,kmey, & ipsy,ipey,jpsy,jpey,kpsy,kpey, & grid%kmax_ideal, grid%lmax_ideal, & grid%num_stoch_levels,grid%num_stoch_levels, & grid%num_stoch_levels,grid%num_stoch_levels, & config_flags%restart, grid%iseedarr_sppt, & config_flags%seed_dim, & grid%DX,grid%DY,grid%sppt_vertstruc, & grid%rstoch, & grid%stddev_cutoff_sppt,grid%gridpt_stddev_sppt, & grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT ) enddo endif ENDIF ! Initialize RAND_PERTURB (4) IF (grid%rand_perturb_on==1) then if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then call rand_seed (config_flags, grid%ISEED_RAND_PERT, grid%iseedarr_rand_pert , 1, config_flags%seed_dim) else call read_write_stochrestart ('stoch_ranpert',grid,config_flags,grid%lmax_ideal,grid%kmax_ideal, & grid%SPFORCS,grid%SPFORCC,2) ! broadcast SPFORC and SPFORS to all processes #ifdef DM_PARALLEL CALL wrf_dm_bcast_real ( grid%SPFORCS , grid%kmax_ideal*grid%lmax_ideal ) CALL wrf_dm_bcast_real ( grid%SPFORCC , grid%kmax_ideal*grid%lmax_ideal ) #endif endif call SETUP_RAND_PERTURB('R', & grid%rand_pert_vertstruc,config_flags%restart, & grid%SP_AMP, & grid%SPFORCC,grid%SPFORCS,grid%ALPH_RAND, & grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT, & grid%KMINFORCT,grid%KMAXFORCT, & grid%LMINFORCT,grid%LMAXFORCT, & grid%KMAXFORCTH,grid%LMAXFORCTH, & grid%time_step,grid%DX,grid%DY, & grid%stepstoch, & grid%gridpt_stddev_rand_pert, & grid%lengthscale_rand_pert, & grid%timescale_rand_pert, & grid%TOT_BACKSCAT_PSI,grid%ZTAU_PSI, & grid%REXPONENT_PSI, & grid%kmax_ideal, grid%lmax_ideal, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then do k = 1,maxspinup CALL RAND_PERT_UPDATE(grid,'T', & grid%SPFORCS,grid%SPFORCC, & grid%SP_AMP,grid%ALPH_RAND, & ips, ipe, jps, jpe, kps, kpe, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & kts, kte, & imsx,imex,jmsx,jmex,kmsx,kmex, & ipsx,ipex,jpsx,jpex,kpsx,kpex, & imsy,imey,jmsy,jmey,kmsy,kmey, & ipsy,ipey,jpsy,jpey,kpsy,kpey, & grid%kmax_ideal, grid%lmax_ideal, & grid%num_stoch_levels,grid%num_stoch_levels, & grid%num_stoch_levels,grid%num_stoch_levels, & config_flags%restart, grid%iseedarr_rand_pert, & config_flags%seed_dim, & grid%DX,grid%DY,grid%rand_pert_vertstruc, & grid%RAND_PERT, & grid%gridpt_stddev_rand_pert, & grid%lengthscale_rand_pert, & grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT ) enddo ENDIF !rand_perturb_on ENDIF ! Initialize Stochastic Parameter Perturbations to convection scheme IF (grid%spp_conv==1) then if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then call rand_seed (config_flags, grid%iseed_spp_conv, grid%iseedarr_spp_conv , 1, config_flags%seed_dim) else call read_write_stochrestart ('stoch_spp_con',grid,config_flags,grid%lmax_ideal,grid%kmax_ideal, & grid%SPFORCS2,grid%SPFORCC2,2) ! broadcast SPFORC and SPFORS to all processes #ifdef DM_PARALLEL CALL wrf_dm_bcast_real ( grid%SPFORCS2 ,grid%kmax_ideal*grid%lmax_ideal ) CALL wrf_dm_bcast_real ( grid%SPFORCC2 ,grid%kmax_ideal*grid%lmax_ideal ) #endif endif call SETUP_RAND_PERTURB('R', & grid%vertstruc_spp_conv,config_flags%restart, & grid%SP_AMP2, & grid%SPFORCC2,grid%SPFORCS2,grid%ALPH_RAND2, & grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT, & grid%KMINFORCT,grid%KMAXFORCT, & grid%LMINFORCT,grid%LMAXFORCT, & grid%KMAXFORCTH,grid%LMAXFORCTH, & grid%time_step,grid%DX,grid%DY, & grid%stepstoch, & grid%gridpt_stddev_spp_conv, & grid%lengthscale_spp_conv, & grid%timescale_spp_conv, & grid%TOT_BACKSCAT_PSI,grid%ZTAU_PSI, & grid%REXPONENT_PSI, & grid%kmax_ideal, grid%lmax_ideal, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then do k = 1,maxspinup CALL RAND_PERT_UPDATE(grid,'T', & grid%SPFORCS2,grid%SPFORCC2, & grid%SP_AMP2,grid%ALPH_RAND2, & ips, ipe, jps, jpe, kps, kpe, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & kts, kte, & imsx,imex,jmsx,jmex,kmsx,kmex, & ipsx,ipex,jpsx,jpex,kpsx,kpex, & imsy,imey,jmsy,jmey,kmsy,kmey, & ipsy,ipey,jpsy,jpey,kpsy,kpey, & grid%kmax_ideal, grid%lmax_ideal, & grid%num_stoch_levels,grid%num_stoch_levels, & grid%num_stoch_levels,grid%num_stoch_levels, & config_flags%restart, grid%iseedarr_spp_conv, & config_flags%seed_dim, & grid%DX,grid%DY,grid%vertstruc_spp_conv, & grid%pattern_spp_conv, & grid%stddev_cutoff_spp_conv,grid%gridpt_stddev_spp_conv, & grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT ) enddo endif ENDIF ! Initialize Stochastic Parameter Peturbations (SPP) to PBL scheme IF (grid%spp_pbl==1) then if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then call rand_seed (config_flags, grid%iseed_spp_pbl, grid%iseedarr_spp_pbl , 1, config_flags%seed_dim) else call read_write_stochrestart ('stoch_spp_pbl',grid,config_flags,grid%lmax_ideal,grid%kmax_ideal, & grid%SPFORCS3,grid%SPFORCC3,2) ! broadcast SPFORC and SPFORS to all processes #ifdef DM_PARALLEL CALL wrf_dm_bcast_real ( grid%SPFORCS3,grid%kmax_ideal*grid%lmax_ideal ) CALL wrf_dm_bcast_real ( grid%SPFORCC3,grid%kmax_ideal*grid%lmax_ideal ) #endif endif call SETUP_RAND_PERTURB('R', & grid%vertstruc_spp_pbl,config_flags%restart, & grid%SP_AMP3, & grid%SPFORCC3,grid%SPFORCS3,grid%ALPH_RAND3, & grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT, & grid%KMINFORCT,grid%KMAXFORCT, & grid%LMINFORCT,grid%LMAXFORCT, & grid%KMAXFORCTH,grid%LMAXFORCTH, & grid%time_step,grid%DX,grid%DY, & grid%stepstoch, & grid%gridpt_stddev_spp_pbl, & grid%lengthscale_spp_pbl, & grid%timescale_spp_pbl, & grid%TOT_BACKSCAT_PSI,grid%ZTAU_PSI, & grid%REXPONENT_PSI, & grid%kmax_ideal, grid%lmax_ideal, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then do k = 1,maxspinup CALL RAND_PERT_UPDATE(grid,'T', & grid%SPFORCS3,grid%SPFORCC3, & grid%SP_AMP3,grid%ALPH_RAND3, & ips, ipe, jps, jpe, kps, kpe, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & kts, kte, & imsx,imex,jmsx,jmex,kmsx,kmex, & ipsx,ipex,jpsx,jpex,kpsx,kpex, & imsy,imey,jmsy,jmey,kmsy,kmey, & ipsy,ipey,jpsy,jpey,kpsy,kpey, & grid%kmax_ideal, grid%lmax_ideal, & grid%num_stoch_levels,grid%num_stoch_levels, & grid%num_stoch_levels,grid%num_stoch_levels, & config_flags%restart, grid%iseedarr_spp_pbl, & config_flags%seed_dim, & grid%DX,grid%DY,grid%vertstruc_spp_pbl, & grid%pattern_spp_pbl, & grid%stddev_cutoff_spp_pbl,grid%gridpt_stddev_spp_pbl, & grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT ) enddo endif ENDIF ! Initialize Stochastic Parameter Peturbations (SPP) to MP Thompson scheme IF (grid%spp_mp .GE. 1) then if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then call rand_seed (config_flags, grid%iseed_spp_mp, grid%iseedarr_spp_mp, 1, config_flags%seed_dim) else call read_write_stochrestart ('stoch_spp_mph',grid,config_flags,grid%lmax_ideal,grid%kmax_ideal, & grid%SPFORCS5,grid%SPFORCC5,2) ! broadcast SPFORC and SPFORS to all processes #ifdef DM_PARALLEL CALL wrf_dm_bcast_real ( grid%SPFORCS5,grid%kmax_ideal*grid%lmax_ideal ) CALL wrf_dm_bcast_real ( grid%SPFORCC5,grid%kmax_ideal*grid%lmax_ideal ) #endif endif call SETUP_RAND_PERTURB('R', & grid%vertstruc_spp_mp,config_flags%restart, & grid%SP_AMP5, & grid%SPFORCC5,grid%SPFORCS5,grid%ALPH_RAND5, & grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT, & grid%KMINFORCT,grid%KMAXFORCT, & grid%LMINFORCT,grid%LMAXFORCT, & grid%KMAXFORCTH,grid%LMAXFORCTH, & grid%time_step,grid%DX,grid%DY, & grid%stepstoch, & grid%gridpt_stddev_spp_mp, & grid%lengthscale_spp_mp, & grid%timescale_spp_mp, & grid%TOT_BACKSCAT_PSI,grid%ZTAU_PSI, & grid%REXPONENT_PSI, & grid%kmax_ideal, grid%lmax_ideal, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then do k = 1,maxspinup CALL RAND_PERT_UPDATE(grid,'T', & grid%SPFORCS5,grid%SPFORCC5, & grid%SP_AMP5,grid%ALPH_RAND5, & ips, ipe, jps, jpe, kps, kpe, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & kts, kte, & imsx,imex,jmsx,jmex,kmsx,kmex, & ipsx,ipex,jpsx,jpex,kpsx,kpex, & imsy,imey,jmsy,jmey,kmsy,kmey, & ipsy,ipey,jpsy,jpey,kpsy,kpey, & grid%kmax_ideal, grid%lmax_ideal, & grid%num_stoch_levels,grid%num_stoch_levels, & grid%num_stoch_levels,grid%num_stoch_levels, & config_flags%restart, grid%iseedarr_spp_mp, & config_flags%seed_dim, & grid%DX,grid%DY,grid%vertstruc_spp_mp, & grid%pattern_spp_mp, & grid%stddev_cutoff_spp_mp,grid%gridpt_stddev_spp_mp,& grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT ) enddo endif ENDIF ! Initialize Stochastic Parameter Peturbations (SPP) to LSM scheme IF (grid%spp_lsm==1) then if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then call rand_seed (config_flags, grid%iseed_spp_lsm, grid%iseedarr_spp_lsm , 1, config_flags%seed_dim) else call read_write_stochrestart ('stoch_spp_lsm',grid,config_flags,grid%lmax_ideal,grid%kmax_ideal, & grid%SPFORCS4,grid%SPFORCC4,2) ! broadcast SPFORC and SPFORS to all processes #ifdef DM_PARALLEL CALL wrf_dm_bcast_real ( grid%SPFORCS4,grid%kmax_ideal*grid%lmax_ideal ) CALL wrf_dm_bcast_real ( grid%SPFORCC4,grid%kmax_ideal*grid%lmax_ideal ) #endif endif call SETUP_RAND_PERTURB('R', & grid%vertstruc_spp_lsm,config_flags%restart, & grid%SP_AMP4, & grid%SPFORCC4,grid%SPFORCS4,grid%ALPH_RAND4, & grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT, & grid%KMINFORCT,grid%KMAXFORCT, & grid%LMINFORCT,grid%LMAXFORCT, & grid%KMAXFORCTH,grid%LMAXFORCTH, & grid%time_step,grid%DX,grid%DY, & grid%stepstoch, & grid%gridpt_stddev_spp_lsm, & grid%lengthscale_spp_lsm, & grid%timescale_spp_lsm, & grid%TOT_BACKSCAT_PSI,grid%ZTAU_PSI, & grid%REXPONENT_PSI, & grid%kmax_ideal, grid%lmax_ideal, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) if ((.not.config_flags%restart) .and. (.not.config_flags%hrrr_cycling)) then do k = 1,maxspinup CALL RAND_PERT_UPDATE(grid,'T', & grid%SPFORCS4,grid%SPFORCC4, & grid%SP_AMP4,grid%ALPH_RAND4, & ips, ipe, jps, jpe, kps, kpe, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & kts, kte, & imsx,imex,jmsx,jmex,kmsx,kmex, & ipsx,ipex,jpsx,jpex,kpsx,kpex, & imsy,imey,jmsy,jmey,kmsy,kmey, & ipsy,ipey,jpsy,jpey,kpsy,kpey, & grid%kmax_ideal, grid%lmax_ideal, & grid%num_stoch_levels,grid%num_stoch_levels, & grid%num_stoch_levels,grid%num_stoch_levels, & config_flags%restart, grid%iseedarr_spp_lsm, & config_flags%seed_dim, & grid%DX,grid%DY,grid%vertstruc_spp_lsm, & grid%pattern_spp_lsm, & grid%stddev_cutoff_spp_lsm,grid%gridpt_stddev_spp_lsm,& grid%VERTSTRUCC,grid%VERTSTRUCS,grid%VERTAMPT ) enddo endif ENDIF ENDIF ! skebs or sppt or rand_perturb or spp END SUBROUTINE INITIALIZE_STOCH ! --- END SETUP STOCHASTIC PERTURBATION SCHEMES ---------- subroutine SETUP_RAND_PERTURB( variable_in,& skebs_vertstruc,restart, & SP_AMP,SPFORCC,SPFORCS,ALPH, & VERTSTRUCC,VERTSTRUCS,VERTAMP, & KMINFORCT,KMAXFORCTH,LMINFORCT,LMAXFORCTH, & KMAXFORCT,LMAXFORCT, & itime_step,DX,DY, & stepstoch, & gridpt_stddev_rand_perturb, l_rand_perturb, & tau_rand_perturb, & TOT_BACKSCAT,ZTAU,REXPONENT, & kmax_ideal, lmax_ideal, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! General control LOGICAL :: restart REAL, PARAMETER :: RPI= 3.141592653589793 !4.0*atan(1.0) CHARACTER, INTENT(IN) :: variable_in ! W=SKEBS_PSI, T=SKEBS_T, P=SPPT, R=RAND_PERTURB CHARACTER :: variable ! Common to all schemes INTEGER , INTENT(IN) :: kmax_ideal, lmax_ideal INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte INTEGER :: IER,IK,IL,I,J,itime_step,stepstoch,skebs_vertstruc, & KMINFORCT,LMINFORCT,KMAXFORCT,LMAXFORCT,KMAXFORCTH,LMAXFORCTH, & KMAX,LMAX,LENSAV,ILEV REAL :: DX,DY,RY,RX,ALPH,RHOKLMAX,ZREF,RHOKL,EPS REAL, DIMENSION (lmax_ideal,kmax_ideal) :: SPFORCS,SPFORCC,SP_AMP REAL, DIMENSION (ims:ime,kms:kme,jms:jme) :: VERTSTRUCC,VERTSTRUCS REAL, DIMENSION (kms:kme) :: VERTAMP REAL, DIMENSION (ids:ide,jds:jde) :: ZCHI ! SPPT and perturb_rand specific REAL :: gridpt_stddev_rand_perturb,kappat,tau_rand_perturb,l_rand_perturb REAL, DIMENSION (ims:ime,jms:jme) :: var_sigma1 ! SKEBS specific REAL :: z,phi,ZGAMMAN,ZCONSTF0,TOT_BACKSCAT,ZTAU,REXPONENT,ZSIGMA2 LOGICAL :: is_print = .true. variable = variable_in ! --------- SETUP PARAMETERS --------------------------------------- KMAX=(jde-jds)+1 !NLAT LMAX=(ide-ids)+1 !NLON RY= KMAX*DY RX= LMAX*DX LENSAV= 4*(KMAX+LMAX)+INT(LOG(REAL(KMAX))) + INT(LOG(REAL(LMAX))) + 8 LENSAV=MAX(LENSAV,4*(KMAX_ideal+LMAX_ideal)+INT(LOG(REAL(KMAX_ideal))) + INT(LOG(REAL(LMAX_ideal))) + 8) ! --------- ALLOCATE FIELDS FOR FFTPACK---------------------------- IF ( ALLOCATED(WSAVE1) ) DEALLOCATE(WSAVE1) IF ( ALLOCATED(WSAVE2) ) DEALLOCATE(WSAVE2) IF ( ALLOCATED(WSAVE1_ideal) ) DEALLOCATE(WSAVE1_ideal) IF ( ALLOCATED(WSAVE2_ideal) ) DEALLOCATE(WSAVE2_ideal) IF ( ALLOCATED(rindarrayil) ) DEALLOCATE(rindarrayil) IF ( ALLOCATED(rindarrayik) ) DEALLOCATE(rindarrayik) IF ( ALLOCATED(rindarrayl) ) DEALLOCATE(rindarrayl) IF ( ALLOCATED(rindarrayk) ) DEALLOCATE(rindarrayk) IF ( ALLOCATED(WAVENUMBER_K)) DEALLOCATE(WAVENUMBER_K) IF ( ALLOCATED(WAVENUMBER_L)) DEALLOCATE(WAVENUMBER_L) ALLOCATE(WSAVE1(LENSAV),WSAVE2(LENSAV)) ALLOCATE(WSAVE1_ideal(LENSAV),WSAVE2_ideal(LENSAV)) ALLOCATE(rindarrayil(LMAX_ideal),rindarrayik(KMAX_ideal),rindarrayl(LMAX),rindarrayk(KMAX)) ALLOCATE (wavenumber_k(jds:jde),wavenumber_l(ids:ide)) ! -------- INITIALIZE arrays to zero ----------------------------- SP_AMP=0.0 ! -------- INITIALIZE FFTPACK ROUTINES ----------------------------- call CFFT1I (LMAX, WSAVE1, LENSAV, IER) if(ier.ne. 0) write(*,95) ier call CFFT1I (KMAX, WSAVE2, LENSAV, IER) if(ier.ne. 0) write(*,95) ier call CFFT1I (LMAX_ideal, WSAVE1_ideal, LENSAV, IER) if(ier.ne. 0) write(*,95) ier call CFFT1I (KMAX_ideal, WSAVE2_ideal, LENSAV, IER) if(ier.ne. 0) write(*,95) ier ! Index arrays for spline interpolations - should be moved to setup routine do i=1,LMAX_ideal rindarrayil(i)=i enddo do i=1,KMAX_ideal rindarrayik(i)=i enddo do i=1,LMAX rindarrayl(i) = 1. + real(LMAX_ideal-1)/(LMAX-1)*(i-1) enddo do i=1,KMAX rindarrayk(i) = 1. + real(KMAX_ideal-1)/(KMAX-1)*(i-1) enddo 95 format('error in cFFT2I= ',i5) call findindex( wavenumber_k, wavenumber_l, & KMAX_ideal, LMAX_ideal) ! set maximal perturbed wavenumber based on gridpoints in domain KMAXFORCT=min0(((ide-ids)+1)/2,((jde-jds)+1 )/2)-5 LMAXFORCT=KMAXFORCT if (KMAXFORCT > KMAXFORCTH) then KMAXFORCT=KMAXFORCTH endif if (LMAXFORCT > LMAXFORCTH) then LMAXFORCT=LMAXFORCTH endif ! -------------------------------------------------------------------------------------- ! ---------- INITIALIZE STOCHASTIC KINETIC-ENERGY BACKSCATTER SCHEME (SKEBS) ---------- ! -------------------------------------------------------------------------------------- ALPH = float(itime_step)/ZTAU/float(stepstoch) ! approximation of 1.-exp(-itime_step/ZTAU_PSI) ZSIGMA2=1./(12.0*ALPH) if (is_print) then IF (variable == 'W') then WRITE(*,'('' '')') WRITE(*,'('' =============================================='')') WRITE(*,'('' >> Initializing STREAMFUNCTION forcing pattern of << '')') WRITE(*,'('' >> stochastic kinetic-energy backscatter scheme << '')') WRITE(*,'('' Total backscattered energy, TOT_BACKSCAT_PSI '',E12.5)') TOT_BACKSCAT WRITE(*,'('' Exponent for energy spectra, REXPONENT_PSI ='',E12.5)') REXPONENT WRITE(*,'('' Minimal wavenumber of streamfunction forcing, LMINFORC ='',I10)') LMINFORCT WRITE(*,'('' Maximal wavenumber of streamfunction forcing, LMAXFORC ='',I10)') LMAXFORCT WRITE(*,'('' Minimal wavenumber of streamfunction forcing, KMINFORC ='',I10)') KMINFORCT WRITE(*,'('' Maximal wavenumber of streamfunction forcing, KMAXFORC ='',I10)') KMAXFORCT WRITE(*,'('' skebs_vertstruc '',I10)') skebs_vertstruc WRITE(*,'('' Time step: itime_step='',I10)') itime_step WRITE(*,'('' Decorrelation time of noise, ZTAU_PSI ='',E12.5)') ZTAU WRITE(*,'('' Variance of noise, ZSIGMA2_EPS ='',E12.5)') ZSIGMA2 WRITE(*,'('' Autoregressive parameter 1-ALPH_PSI ='',E12.5)') 1.-ALPH WRITE(*,'('' =============================================='')') ! Unit of SPSTREAM_AMP: sqrt(m^2/s^3 1/s m**2(p+1)) m**-2(p/2) = m^/s^2 * m**[(p+1)-p] = m^2/s^2 m ELSEIF (variable == 'T') then WRITE(*,'('' '')') WRITE(*,'('' =============================================='')') WRITE(*,'('' >> Initializing TEMPERATURE forcing pattern of << '')') WRITE(*,'('' >> stochastic kinetic-energy backscatter scheme << '')') WRITE(*,'('' Total backscattered energy, TOT_BACKSCAT_T '',E12.5)') TOT_BACKSCAT WRITE(*,'('' Exponent for energy spectra, REXPONENT_T ='',E12.5)') REXPONENT WRITE(*,'('' Minimal wavenumber of tempearature forcing, LMINFORC ='',I10)') LMINFORCT WRITE(*,'('' Maximal wavenumber of tempearature forcing, LMAXFORC ='',I10)') LMAXFORCT WRITE(*,'('' Minimal wavenumber of tempearature forcing, KMINFORC ='',I10)') KMINFORCT WRITE(*,'('' Maximal wavenumber of tempearature forcing, KMAXFORC ='',I10)') KMAXFORCT WRITE(*,'('' skebs_vertstruc '',I10)') skebs_vertstruc WRITE(*,'('' Decorrelation time of noise, ZTAU_T ='',E12.5)') ZTAU WRITE(*,'('' Variance of noise, ZSIGMA2_ETA ='',E12.5)') ZSIGMA2 WRITE(*,'('' Autoregressive parameter 1-ALPH_T ='',E12.5)') 1.-ALPH WRITE(*,'('' =============================================='')') endif IF ((variable == 'P') .or. (variable == 'R')) then kappat= L_rand_perturb**2 ! L^2= kappa*T, where L is a length scale in m; set to for L=100km phi = exp (-float(itime_step*stepstoch)/tau_rand_perturb) alph = 1.-phi endif ! -------------------------------------------------------------------------------------- ! ---------- INITIALIZE STOCHASTICALLY PERTURBED PHYSICAL TENDENCY SCHEME -------------- ! -------------------------------------------------------------------------------------- if (variable == 'P') then WRITE(*,'('' '')') WRITE(*,'('' =============================================='')') WRITE(*,'('' >> Initializing Stochastically Perturbed Physics Tendency scheme << '')') WRITE(*,'('' sppt_vertstruc '',I10)') skebs_vertstruc WRITE(*,'('' Decorrelation time of noise, Tau ='',E12.5)') tau_rand_perturb WRITE(*,'('' Autoregressive parameter Phi ='',E12.5)') phi WRITE(*,'('' Length Scale L'',E12.5)') l_rand_perturb WRITE(*,'('' Variance in gridpoint space'',E12.5)') gridpt_stddev_rand_perturb WRITE(*,'('' =============================================='')') endif ! variable ! -------------------------------------------------------------------------------------- ! -------------------- INITIALIZE RANDOM PERTUBATIONS ------------------------------- ! -------------------------------------------------------------------------------------- if (variable == 'R') then WRITE(*,'('' '')') WRITE(*,'('' =============================================='')') WRITE(*,'('' >> Initializing random perturbations << '')') WRITE(*,'('' rand_pert_vertstruc '',I10)') skebs_vertstruc WRITE(*,'('' Decorrelation time of noise, Tau ='',E12.5)') tau_rand_perturb WRITE(*,'('' Autoregressive parameter Phi ='',E12.5)') phi WRITE(*,'('' Length Scale L'',E12.5)') l_rand_perturb WRITE(*,'('' Variance in gridpoint space'',E12.5)') gridpt_stddev_rand_perturb WRITE(*,'('' =============================================='')') endif ! variable endif !is print ! -------------------------------------------------------------------------------------- ! Compute Normalization constants ! -------------------------------------------------------------------------------------- ZCHI = 0.0 ZGAMMAN = 0.0 ! Fill lower left quadrant of ZCHI. For this range the indeces IK,IL DO IK=jds-1,jde ! These are now wavenumbers DO IL=ids-1,ide if (((sqrt((IK/RY*IK/RY)+(IL/RX*IL/RX)).lt.((KMAXFORCT+0.5)/RX)).and.& (sqrt((IK/RY*IK/RY)+(IL/RX*IL/RX)).ge.((KMINFORCT-0.5)/RX))) .or. & ((sqrt((IK/RY*IK/RY)+(IL/RX*IL/RX)).lt.((LMAXFORCT+0.5)/RY)).and.& (sqrt((IK/RY*IK/RY)+(IL/RX*IL/RX)).ge.((LMINFORCT-0.5)/RY))))then if ((IK>0).or.(IL>0)) then if (variable == 'W') then ZCHI(IL+1,IK+1)=((IK/RY*IK/RY)+(IL/RX*IL/RX))**(REXPONENT/2.) ! SKEBS :U ZGAMMAN= ZGAMMAN + ((IK/RY*IK/RY)+(IL/RX*IL/RX))**(REXPONENT+1) else if (variable == 'T') then ZCHI(IL+1,IK+1)=((IK/RY*IK/RY)+(IL/RX*IL/RX))**(REXPONENT/2.) ! SKEBS :T ZGAMMAN= ZGAMMAN + ((IK/RY*IK/RY)+(IL/RX*IL/RX))**(REXPONENT) else if ((variable == 'P') .or. (variable == 'R')) then ZCHI(IL+1,IK+1)=exp( -2*RPI**2*kappat*((IK/RY*IK/RY)+(IL/RX*IL/RX)) ) !SPPT ZGAMMAN= ZGAMMAN + exp( -4*RPI**2*kappat*((IK/RY*IK/RY)+(IL/RX*IL/RX)) ) !SPPT endif endif endif enddo enddo ZGAMMAN=4.0*ZGAMMAN !account for all quadrants, although only one is Filled if (variable == 'W') then ZCONSTF0=SQRT(ALPH*TOT_BACKSCAT/(float(itime_step)*ZSIGMA2*ZGAMMAN))/(2*RPI) elseif (variable == 'T') then ! Mod by Nick to avoid precision problems !ZCONSTF0=SQRT(T0*ALPH*TOT_BACKSCAT/(float(itime_step)*cp*ZSIGMA2*ZGAMMAN)) ZCONSTF0=SQRT(T0*ALPH*TOT_BACKSCAT)/SQRT(float(itime_step)*cp*ZSIGMA2) ZCONSTF0=ZCONSTF0/SQRT(ZGAMMAN) elseif ((variable == 'P') .or. (variable == 'R')) then ZCONSTF0= gridpt_stddev_rand_perturb*sqrt((1.-phi**2)/(2.*ZGAMMAN)) endif ! -------------------------------------------------------------------------------------- ! Now the wavenumber-dependent amplitudes ! -------------------------------------------------------------------------------------- ! Note: There are symmetries and anti-symmetries to ensure real-valued back transforms ! Fill lower left quadrant of matrix of noise amplitudes for wavenumbers K=0,KMAX/2 ! SP_AMP=0.0 DO IL=1, (LMAX_ideal/2+1)! upper left DO IK = 1, (KMAX_ideal/2+1) SP_AMP(IL,IK)=ZCONSTF0*ZCHI(IL,IK) ENDDO ENDDO DO IL=1, (LMAX_ideal/2+1) DO IK = 1, (KMAX_ideal/2+1) SP_AMP(LMAX_ideal-IL+1,KMAX_ideal-IK+1)= SP_AMP(IL,IK) !lower right SP_AMP(IL,KMAX_ideal-IK+1)= SP_AMP(IL,IK) !lower right SP_AMP(LMAX_ideal-IL+1,IK)= SP_AMP(IL,IK) !lower right ENDDO ENDDO ! ----------------------------------------- ! ----------------------------------------- ! Array for vertical structure if desired VERTAMP=1.0 ! Define vertical amplitude here. IF (skebs_vertstruc==1) then VERTSTRUCC=0.0 VERTSTRUCS=0.0 RHOKLMAX= sqrt(KMAX**2/DY**2 + LMAX**2/DX**2) ZREF=32.0 DO ILEV=kts,kte DO IK=jts,jte DO IL=its,ite if (IL.le.(LMAX/2)) then RHOKL = sqrt((IK+1)**2/DY**2 + (IL+1)**2/DX**2) EPS = ((RHOKLMAX - RHOKL)/ RHOKLMAX) * (ILEV/ZREF) * RPI VERTSTRUCC(IL,ILEV,IK) = cos ( eps* (IL+1) ) VERTSTRUCS(IL,ILEV,IK) = sin ( eps* (IL+1) ) else RHOKL = sqrt((IK+1)**2/DY**2 + (LMAX-IL+2)**2/DX**2) EPS = ((RHOKLMAX - RHOKL)/ RHOKLMAX) * (ILEV/ZREF) * RPI VERTSTRUCC (IL,ILEV,IK) = cos ( eps* (LMAX-IL+2) ) VERTSTRUCS (IL,ILEV,IK) = - sin ( eps* (LMAX-IL+2) ) endif ENDDO ENDDO ENDDO ENDIF END subroutine SETUP_RAND_PERTURB ! ------------------------------------------------------------------ !************** UPDATE STOCHASTIC PATTERN IN WAVENUMBER SPACE********** ! ------------------------------------------------------------------ subroutine UPDATE_STOCH( & SPFORCS,SPFORCC,SP_AMP,ALPH, & restart,iseedarr,seed_dim, & kmax_ideal, lmax_ideal, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE REAL, DIMENSION(lmax_ideal,kmax_ideal) ,INTENT(INOUT) :: SPFORCS,SPFORCC,SP_AMP INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte INTEGER , INTENT(IN) :: seed_dim INTEGER, DIMENSION (seed_dim), INTENT(INOUT) :: iseedarr INTEGER , INTENT(IN) :: kmax_ideal, lmax_ideal REAL, DIMENSION(LMAX_ideal,KMAX_ideal) :: ZRANDNOSS,ZRANDNOSC INTEGER , ALLOCATABLE , DIMENSION(:) :: iseed REAL :: Z,ALPH REAL, PARAMETER :: thresh = 3.0 INTEGER ::IL, IK,LMAX,KMAX,J INTEGER :: how_many LOGICAL :: LGAUSS,RESTART KMAX=(jde-jds)+1 !NLAT LMAX=(ide-ids)+1 !NATX call random_seed(put=iseedarr) ! Pick the distribution of the noise LGAUSS=.true. IF (LGAUSS) then DO IK=1,KMAX_ideal DO IL=1,LMAX_ideal do call gauss_noise(z) if (abs(z)1)) then ! Upper half DO IL=1,LMAX_ideal SPFORCC(IL,IK) = (1.-ALPH)*SPFORCC(IL,IK) + SP_AMP(IL,IK) * ZRANDNOSC(IL,IK) SPFORCS(IL,IK) = (1.-ALPH)*SPFORCS(IL,IK) + SP_AMP(IL,IK) * ZRANDNOSS(IL,IK) ENDDO ELSEIF (IK==1) then DO IL=1,LMAX_ideal if ((IL.le.(LMAX_ideal/2+1))) then SPFORCC(IL,IK) = (1.-ALPH)*SPFORCC(IL,IK) + SP_AMP(IL,IK) * ZRANDNOSC(IL,IK) SPFORCS(IL,IK) = (1.-ALPH)*SPFORCS(IL,IK) + SP_AMP(IL,IK) * ZRANDNOSS(IL,IK) elseif ((IL.gt.(LMAX_ideal/2+1))) then SPFORCC(IL,IK) = (1.-ALPH)*SPFORCC(IL,IK) + SP_AMP(IL,IK) * ZRANDNOSC(LMAX_ideal-IL+2,IK) SPFORCS(IL,IK) = (1.-ALPH)*SPFORCS(IL,IK) - SP_AMP(IL,IK) * ZRANDNOSS(LMAX_ideal-IL+2,IK) endif ENDDO ENDIF ENDDO DO IK=1,KMAX_ideal if (IK.gt.(KMAX_ideal/2+1)) then ! Lower half DO IL=1,LMAX_ideal if (IL.le.(LMAX_ideal/2+1).and.(IL.gt.1)) then !lower left SPFORCC(IL,IK) = (1.-ALPH)* SPFORCC(IL,IK) + SP_AMP(IL,IK) * ZRANDNOSC(LMAX_ideal-IL+2,KMAX_ideal-IK+2) SPFORCS(IL,IK) = (1.-ALPH)* SPFORCS(IL,IK) - SP_AMP(IL,IK) * ZRANDNOSS(LMAX_ideal-IL+2,KMAX_ideal-IK+2) elseif (IL.eq.1) then !don't exceed index SPFORCC(IL,IK) = (1.-ALPH)* SPFORCC(IL,IK) + SP_AMP(IL,IK) * ZRANDNOSC( 1,KMAX_ideal-IK+2) SPFORCS(IL,IK) = (1.-ALPH)* SPFORCS(IL,IK) - SP_AMP(IL,IK) * ZRANDNOSS( 1,KMAX_ideal-IK+2) elseif (IL.gt.(LMAX_ideal/2+1)) then !lower right SPFORCC(IL,IK) = (1.-ALPH)* SPFORCC(IL,IK) + SP_AMP(IL,IK) * ZRANDNOSC(LMAX_ideal-IL+2,KMAX_ideal-IK+2) SPFORCS(IL,IK) = (1.-ALPH)* SPFORCS(IL,IK) - SP_AMP(IL,IK) * ZRANDNOSS(LMAX_ideal-IL+2,KMAX_ideal-IK+2) endif ENDDO endif ENDDO call random_seed(get=iseedarr) END subroutine UPDATE_STOCH ! ------------------------------------------------------------------ SUBROUTINE UPDATE_STOCH_TEN(ru_tendf,rv_tendf,t_tendf, & ru_tendf_stoch,rv_tendf_stoch,rt_tendf_stoch,& mu,mub,c1h,c2h, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte, & kte_stoch,kme_stoch ) IMPLICIT NONE INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte, & kte_stoch,kme_stoch REAL , DIMENSION(ims:ime , kms:kme, jms:jme),INTENT(INOUT) :: & ru_tendf, rv_tendf, t_tendf REAL , DIMENSION(ims:ime , kms:kme_stoch, jms:jme) :: & ru_tendf_stoch,rv_tendf_stoch,rt_tendf_stoch REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: mu,mub REAL , DIMENSION(kms:kme) , INTENT(IN) :: c1h,c2h INTEGER :: I,J,K,kh REAL :: dt,xm DO j = jts,MIN(jde-1,jte) DO k = kts,kte-1 kh=min(k,kte_stoch) DO i = its,ite ru_tendf(i,k,j) = ru_tendf(i,k,j) + ru_tendf_stoch(i,kh,j) * (mu(i,j)+mub(i,j)) ENDDO ENDDO ENDDO DO j = jts,jte DO k = kts,kte-1 kh=min(k,kte_stoch) DO i = its,MIN(ide-1,ite) rv_tendf(i,k,j) = rv_tendf(i,k,j) + rv_tendf_stoch(i,kh,j) * (mu(i,j)+mub(i,j)) ENDDO ENDDO ENDDO DO j = jts,MIN(jde-1,jte) DO k = kts,kte-1 kh=min(k,kte_stoch) DO i = its,MIN(ide-1,ite) t_tendf(i,k,j) = t_tendf(i,k,j) + rt_tendf_stoch(i,kh,j) * (mu(i,j)+mub(i,j)) ENDDO ENDDO ENDDO END SUBROUTINE UPDATE_STOCH_TEN ! ------------------------------------------------------------------ !!************** PERTURB PHYSICS TENDENCIES (except T) FOR SPPT ******************* ! ------------------------------------------------------------------ subroutine perturb_physics_tend(gridpt_stddev_sppt, & sppt_thresh_fact,rstoch, & ru_tendf,rv_tendf,t_tendf,moist_tend, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte, & kte_stoch,kme_stoch ) ! This subroutine add stochastic perturbations of the form ! ! rx_tendf(i,k,j) = rx_tendf(i,k,j)*(1.0 + rstoch(i,k,j)) ! ! to the tendencies of U, V, and Q. ! Since the temperature perturbations do not include the micro-physics ! tendencies at this point, the stochastic tendency perturbations to ! temperature are added in subroutine rk_addtend_dry of module module_em.F IMPLICIT NONE INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte, & kte_stoch,kme_stoch REAL , DIMENSION(ims:ime , kms:kme, jms:jme),INTENT(INOUT) :: & ru_tendf, rv_tendf, t_tendf,moist_tend REAL , DIMENSION(ims:ime,kms:kme_stoch, jms:jme),INTENT(INOUT) :: rstoch REAL :: gridpt_stddev_sppt ,thresh,sppt_thresh_fact INTEGER :: I,J,K,kh ! Here the random process at each gridpoint is capped if it exceeds a value thresh thresh=sppt_thresh_fact*gridpt_stddev_sppt DO j = jts,jte DO k = kts,min(kte-1,kte_stoch-1) DO i = its,ite ! rstoch(i,k,j)=MAX(MIN(rstoch(i,k,j),thresh),-1.*thresh)) if (rstoch(i,k,j).lt.-thresh) then rstoch(i,k,j)=-thresh endif if (rstoch(i,k,j).gt.thresh) then rstoch(i,k,j)=thresh endif ENDDO ENDDO ENDDO ! Perturb the tendencies of u,v,q,t. DO j = jts,MIN(jde-1,jte) DO k = kts,kte-1 kh = min( k, kte_stoch-1 ) DO i = its,ite ru_tendf(i,k,j) = ru_tendf(i,k,j)*(1.0 + rstoch(i,kh,j)) ENDDO ENDDO ENDDO DO j = jts,jte DO k = kts,kte-1 kh = min( k, kte_stoch-1 ) DO i = its,MIN(ide-1,ite) rv_tendf(i,k,j) = rv_tendf(i,k,j)*(1.0 + rstoch(i,kh,j)) ENDDO ENDDO ENDDO DO j = jts,MIN(jde-1,jte) DO k = kts,kte-1 kh = min( k, kte_stoch-1 ) DO i = its,MIN(ide-1,ite) moist_tend(i,k,j) = moist_tend(i,k,j)*(1.0 + rstoch(i,kh,j)) t_tendf (i,k,j) = t_tendf(i,k,j)*(1.0 + rstoch(i,kh,j)) ENDDO ENDDO ENDDO end subroutine perturb_physics_tend ! ------------------------------------------------------------------ !!************** UPDATE SPECTRAL PATTERN AND TRANSFORM TO GRIDPOINT SPACE*** ! ------------------------------------------------------------------ ! This subroutine evolves the spectral pattern and transforms it back to gridpoint space. SUBROUTINE RAND_PERT_UPDATE (grid, variable_in, & SPFORCS,SPFORCC,SP_AMP,ALPH_RAND, & ips, ipe, jps, jpe, kps, kpe, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & kts, kte, & imsx,imex,jmsx,jmex,kmsx,kmex, & ipsx,ipex,jpsx,jpex,kpsx,kpex, & imsy,imey,jmsy,jmey,kmsy,kmey, & ipsy,ipey,jpsy,jpey,kpsy,kpey, & kmax_ideal, lmax_ideal, & kpe_stoch,kde_stoch,kme_stoch,kte_stoch, & restart,iseedarr,seed_dim, & DX,DY,skebs_vertstruc, & RAND_PERT,thresh_fact,gridpt_stddev, & VERTSTRUCC,VERTSTRUCS,VERTAMP ) USE module_domain, ONLY : domain #ifdef DM_PARALLEL USE module_dm, ONLY : local_communicator, mytask, ntasks, ntasks_x, ntasks_y, local_communicator_periodic, & wrf_dm_maxval, wrf_err_message, local_communicator_x, local_communicator_y, data_order_xzy #endif IMPLICIT NONE TYPE ( domain ), INTENT(INOUT) :: grid INTEGER , INTENT(IN) :: kmax_ideal, lmax_ideal INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe, & kts, kte INTEGER , INTENT(IN) :: imsx,imex,jmsx,jmex,kmsx,kmex, & ipsx,ipex,jpsx,jpex,kpsx,kpex, & imsy,imey,jmsy,jmey,kmsy,kmey, & ipsy,ipey,jpsy,jpey,kpsy,kpey INTEGER , INTENT(IN) :: seed_dim INTEGER :: kpe_stoch,kde_stoch,kme_stoch,kte_stoch REAL , INTENT(IN) :: ALPH_RAND,dx,dy,thresh_fact,gridpt_stddev INTEGER , INTENT(IN) :: skebs_vertstruc CHARACTER, INTENT(IN) :: variable_in ! T, U, V ! T ! random field, T ! U ! first derivative of streamfunction with regard to y; for skebs: U ! V ! first derivative of streamfunction with regard to x; for skebs: V INTEGER, DIMENSION (seed_dim), INTENT(INOUT) :: iseedarr REAL, DIMENSION(ims:ime,kms:kme, jms:jme),INTENT(IN) :: VERTSTRUCC,VERTSTRUCS REAL, DIMENSION(lmax_ideal,kmax_ideal) ,INTENT(INOUT) :: SPFORCS,SPFORCC,SP_AMP REAL, DIMENSION(kms:kme ) ,INTENT(IN) :: VERTAMP REAL, DIMENSION(ims:ime,kms:kme_stoch, jms:jme) :: RAND_PERT REAL :: RY,RX ! Local Variabels INTEGER :: IK,IL,ILEV,NLON,NLAT,IJ,I,J,K,numbands COMPLEX, DIMENSION (LMAX_ideal) :: dummy_complex2 COMPLEX, DIMENSION (KMAX_ideal) :: dummy_complex3 COMPLEX, DIMENSION (ipsx:ipex) :: dummy_complex REAL :: thresh REAL, DIMENSION(LMAX_ideal,KMAX_ideal) :: ZFORCC,ZFORCS REAL, DIMENSION(ids:ide,jds:jde) :: ZFORCg INTEGER :: IER,LENWRK,KMAX,LMAX,LENSAV,JJ,II INTEGER :: its,ite,jts,jte,ind REAL, ALLOCATABLE :: WORK(:) CHARACTER (LEN=160) :: mess real, dimension (LMAX_ideal) :: b,c,d ! Assums LMAX_ideal=KMAX_ideal which is the case for homogenous fields LOGICAL :: RESTART CHARACTER :: variable variable = variable_in NLAT=(jde-jds)+1 !KMAX NLON=(ide-ids)+1 !LMAX KMAX = NLAT LMAX = NLON RY= NLAT*DY RX= NLON*DX LENWRK=2*KMAX*LMAX ALLOCATE(WORK(LENWRK)) LENSAV= 4*(KMAX+LMAX)+INT(LOG(REAL(KMAX))) + INT(LOG(REAL(LMAX))) + 8 ! Update the pattern generator by evolving each spectral coefficients as AR1 IF (variable .ne. 'V') THEN !T, random field, U, don't update for V CALL UPDATE_STOCH( & SPFORCS,SPFORCC,SP_AMP,ALPH_RAND, & restart,iseedarr,seed_dim, & kmax_ideal, lmax_ideal, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & kts,kte,kts,kte,kts, kte ) endif ! Take derivatives of streamfunction of for skebs ZFORCC=SPFORCC ZFORCS=SPFORCS IF (variable == 'U') THEN !U DO J = 1, KMAX_ideal DO I = 1, LMAX_ideal ZFORCC(I,J) = 2*RPI/RY* wavenumber_k(J) * ZFORCS(I,J) ZFORCS(I,J) = -2*RPI/RY* wavenumber_k(J) * ZFORCC(I,J) ENDDO ENDDO ELSEIF (variable == 'V') THEN !V DO J = 1, KMAX_ideal DO I = 1, LMAX_ideal ZFORCC(I,J) = -2*RPI/RX* wavenumber_l(I) * ZFORCS(I,J) ZFORCS(I,J) = 2*RPI/RX* wavenumber_l(I) * ZFORCC(I,J) ENDDO ENDDO endif DO J = 1, KMAX_ideal DO i = 1, LMAX_ideal dummy_complex2(i)=cmplx(ZFORCC(i,j),ZFORCS(i,j)) ENDDO CALL cFFT1B (LMAX_ideal, 1 ,dummy_complex2,LMAX_ideal, WSAVE1_ideal, LENSAV, WORK, LENWRK, IER) if (ier.ne.0) then WRITE(mess,FMT='(A)') 'error in cFFT1B in do_fftback_along_x, field U' CALL wrf_debug(0,mess) end if DO i = 1, LMAX_ideal ZFORCC(i,j)=real(dummy_complex2(i)) ZFORCS(i,j)=imag(dummy_complex2(i)) END DO END DO ! J DO I = 1, LMAX_ideal DO j = 1, KMAX_ideal dummy_complex3(j)=cmplx(ZFORCC(i,j),ZFORCS(i,j)) ENDDO CALL cFFT1B (KMAX_ideal, 1 ,dummy_complex3,KMAX_ideal, WSAVE2_ideal, LENSAV, WORK, LENWRK, IER) if (ier.ne.0) then WRITE(mess,FMT='(A)') 'error in cFFT1B in do_fftback_along_y, field U' CALL wrf_debug(0,mess) end if DO j = 1, KMAX_ideal ZFORCC(i,j)=real(dummy_complex3(j)) ZFORCS(i,j)=imag(dummy_complex3(j)) END DO END DO ! I ! Spline interpolation if ((LMAX_ideal.ne.ide).or.(KMAX_ideal.ne.jde)) then ! https://ww2.odu.edu/~agodunov/computing/programs/book2/Ch01/spline.f90 DO J = 1, KMAX_ideal DO I = 1, LMAX_ideal dummy_complex2(i)=cmplx(ZFORCC(i,j),0.0) ENDDO call spline (rindarrayil, real(dummy_complex2) , b, c, d, LMAX_ideal ) do i=1,LMAX ZFORCg(i,j)=ispline (rindarrayl(i),rindarrayil, real(dummy_complex2) , b, c, d, LMAX_ideal ) enddo ENDDO DO I = 1, LMAX DO J = 1, KMAX_ideal dummy_complex3(j)=cmplx(ZFORCg(i,j),0.0) ENDDO call spline (rindarrayik, real(dummy_complex3) , b, c, d, KMAX_ideal ) do J=1,KMAX ZFORCg(i,j)=ispline (rindarrayk(j),rindarrayik, real(dummy_complex3) , b, c, d, KMAX_ideal ) enddo ENDDO else DO J = 1, KMAX_ideal DO I = 1, LMAX_ideal ZFORCg(i,j)=ZFORCC(i,j) ENDDO ENDDO endif ! LMAX_ideal thresh=thresh_fact*gridpt_stddev RAND_PERT=0.0 !$OMP PARALLEL DO & !$OMP PRIVATE ( ij ) DO ij = 1 , grid%num_tiles DO k=kts,min(kte,grid%num_stoch_levels) DO I=grid%i_start(ij), grid%i_end(ij) DO j=grid%j_start(ij), grid%j_end(ij) RAND_PERT(I,K,J)=MAX(MIN(ZFORCg(I,J),thresh),-1.0*thresh) ENDDO ENDDO ENDDO ENDDO !OMP !$OMP END PARALLEL DO END SUBROUTINE RAND_PERT_UPDATE ! ------------------------------------------------------------------ !!************** UPDATE SPECTRAL PATTERN AND TRANSFORM TO GRIDPOINT SPACE*** ! ------------------------------------------------------------------ ! This subroutine evolves the spectral pattern and transforms it back to gridpoint space. subroutine read_write_stochrestart (filename1,grid,config_flags,lmax_ideal,kmax_ideal,SPFORS,SPFORC,nswitch) !lala USE module_configure , ONLY : grid_config_rec_type USE module_date_time USE module_timing USE module_utility USE module_domain #ifdef DM_PARALLEL USE module_dm #endif IMPLICIT NONE TYPE(domain) :: grid TYPE (grid_config_rec_type) , INTENT(IN) :: config_flags REAL, DIMENSION(lmax_ideal,kmax_ideal) ,INTENT(INOUT) :: SPFORC, SPFORS INTEGER , INTENT(IN) :: kmax_ideal,lmax_ideal,nswitch INTEGER :: ierr INTEGER , PARAMETER :: write_restart = 1 INTEGER , PARAMETER :: read_restart = 2 INTEGER , PARAMETER :: fid=19 TYPE(WRFU_Time) :: next_time, currentTime, startTime CHARACTER (LEN=256) :: current_date_char, date_string, next_timestr_char character(len=13) :: filename1 character(len=50) :: filename character(len=70) :: message CHARACTER (len=len_current_date) current_timestr, next_timestr logical itsopen LOGICAL, EXTERNAL :: wrf_dm_on_monitor if(wrf_dm_on_monitor()) then IF ( nswitch .eq. write_restart) then ! WRITE!!! CALL domain_clock_get( grid, current_timestr=current_date_char ) !write(filename,FMT='(A,A,A)'),trim(filename1),'_rst_d01_',TRIM(current_date_char) write(filename,FMT='(A,A,A)'),trim(filename1),'_rst_d01' inquire(unit=19, opened=itsopen) if ( itsopen ) then write(*,*) 'Its open already' else print*,'opening ', filename, ' for writing restart file for stochastic suite' OPEN(19, FILE=filename,FORM='UNFORMATTED',STATUS='REPLACE',IOSTAT=ierr) IF(ierr .NE. 0 ) THEN WRITE(message,FMT='(A)') & 'problem with writing restart files for stochastic physics suite' CALL wrf_error_fatal ( message ) END IF endif write(19),SPFORC,SPFORS close(19) else IF ( nswitch .eq. read_restart) then ! READ!!! CALL domain_clock_get( grid, current_timestr=current_date_char ) !write(filename,FMT='(A,A,A)'),trim(filename1),'_rst_d01_',TRIM(current_date_char) ! The clock has not yet been advanced, hence the time on the stoch restart file is ! istep too early. Want to add tstep to minute ! write(filename,FMT='(A,A)'),'stochrst_d01_','2000-01-24_12:12:00' write(filename,FMT='(A,A,A)'),trim(filename1),'_rst_d01' inquire(unit=19, opened=itsopen) if ( itsopen ) then write(*,*) 'Its open already' else print*,'RESTART run: opening ', filename, ' for reading' OPEN(19, FILE=filename,FORM='UNFORMATTED',STATUS='OLD',IOSTAT=ierr) IF(ierr .NE. 0 ) THEN WRITE(message,FMT='(A)') & 'problem with reading restart files for stochastic physics suite' CALL wrf_error_fatal ( message ) END IF end if read(19),SPFORC,SPFORS close(19) endif endif !wrf_dm_on_monitor END SUBROUTINE read_write_stochrestart ! ------------------------------------------------------------------ !!************** TRANSFORM FROM GRIDPOINT SPACE TO SPHERICAL HARMONICS ** ! ------------------------------------------------------------------ subroutine findindex( wavenumber_k, wavenumber_L, & KMAX_ideal, LMAX_ideal) IMPLICIT NONE INTEGER :: IK,IL,KMAX,LMAX,KMAX_ideal, LMAX_ideal INTEGER, DIMENSION (1:KMAX_ideal):: wavenumber_k INTEGER, DIMENSION (1:LMAX_ideal):: wavenumber_l KMAX=KMAX_ideal LMAX=LMAX_ideal !map wave numbers K,L to indeces IK, IL DO IK=1,KMAX/2+1 wavenumber_k(IK)=IK-1 ENDDO DO IK=KMAX,KMAX/2+2,-1 wavenumber_k(IK)=IK-KMAX-1 ENDDO DO IL=1,LMAX/2+1 wavenumber_l(IL)=IL-1 ENDDO DO IL=LMAX,LMAX/2+2,-1 wavenumber_l(IL)=IL-LMAX-1 ENDDO END subroutine findindex ! ------------------------------------------------------------------ subroutine gauss_noise(z) real :: z ! output real :: x,y,r, coeff ! INPUT ! [2.1] Get two uniform variate random numbers IL range 0 to 1: do call random_number( x ) call random_number( y ) ! [2.2] Transform to range -1 to 1 and calculate sum of squares: x = 2.0 * x - 1.0 y = 2.0 * y - 1.0 r = x * x + y * y if ( r > 0.0 .and. r < 1.0 ) exit end do ! ! [2.3] Use Box-Muller transformation to get normal deviates: coeff = sqrt( -2.0 * log(r) / r ) z = coeff * x end subroutine gauss_noise ! ------------------------------------------------------------------ SUBROUTINE rand_seed (config_flags, iseed1, iseedarr, seed_start, seed_dim ) USE module_configure IMPLICIT NONE ! ! Structure that contains run-time configuration (namelist) data for domain TYPE (grid_config_rec_type) :: config_flags ! ! Arguments INTEGER :: iseed1 , seed_start, seed_dim INTEGER, DIMENSION (seed_start:seed_dim), INTENT(OUT):: iseedarr ! Local integer*8 :: fctime, one_big integer :: i fctime = config_flags%start_year * ( config_flags%start_month*100+config_flags%start_day) + config_flags%start_hour one_big = 1 iseedarr=0 do i = seed_start,seed_dim iseedarr(i)=mod(fctime+iseed1*config_flags%nens*1000000,19211*one_big) enddo end SUBROUTINE rand_seed ! ------------------------------------------------------------------ subroutine spline (x, y, b, c, d, n) !====================================================================== ! Calculate the coefficients b(i), c(i), and d(i), i=1,2,...,n ! for cubic spline interpolation ! s(x) = y(i) + b(i)*(x-x(i)) + c(i)*(x-x(i))**2 + d(i)*(x-x(i))**3 ! for x(i) <= x <= x(i+1) ! Alex G: January 2010 !---------------------------------------------------------------------- ! input.. ! x = the arrays of data abscissas (in strictly increasing order) ! y = the arrays of data ordinates ! n = size of the arrays xi() and yi() (n>=2) ! output.. ! b, c, d = arrays of spline coefficients ! comments ... ! spline.f90 program is based on fortran version of program spline.f ! the accompanying function fspline can be used for interpolation !====================================================================== implicit none integer n real x(n), y(n), b(n), c(n), d(n) integer i, j, gap real h gap = n-1 ! check input if ( n < 2 ) return if ( n < 3 ) then b(1) = (y(2)-y(1))/(x(2)-x(1)) ! linear interpolation c(1) = 0. d(1) = 0. b(2) = b(1) c(2) = 0. d(2) = 0. return end if ! ! step 1: preparation ! d(1) = x(2) - x(1) c(2) = (y(2) - y(1))/d(1) do i = 2, gap d(i) = x(i+1) - x(i) b(i) = 2.0*(d(i-1) + d(i)) c(i+1) = (y(i+1) - y(i))/d(i) c(i) = c(i+1) - c(i) end do ! ! step 2: end conditions ! b(1) = -d(1) b(n) = -d(n-1) c(1) = 0.0 c(n) = 0.0 if(n /= 3) then c(1) = c(3)/(x(4)-x(2)) - c(2)/(x(3)-x(1)) c(n) = c(n-1)/(x(n)-x(n-2)) - c(n-2)/(x(n-1)-x(n-3)) c(1) = c(1)*d(1)**2/(x(4)-x(1)) c(n) = -c(n)*d(n-1)**2/(x(n)-x(n-3)) end if ! ! step 3: forward elimination ! do i = 2, n h = d(i-1)/b(i-1) b(i) = b(i) - h*d(i-1) c(i) = c(i) - h*c(i-1) end do ! ! step 4: back substitution ! c(n) = c(n)/b(n) do j = 1, gap i = n-j c(i) = (c(i) - d(i)*c(i+1))/b(i) end do ! ! step 5: compute spline coefficients ! b(n) = (y(n) - y(gap))/d(gap) + d(gap)*(c(gap) + 2.0*c(n)) do i = 1, gap b(i) = (y(i+1) - y(i))/d(i) - d(i)*(c(i+1) + 2.0*c(i)) d(i) = (c(i+1) - c(i))/d(i) c(i) = 3.*c(i) end do c(n) = 3.0*c(n) d(n) = d(n-1) end subroutine spline function ispline(u, x, y, b, c, d, n) !====================================================================== ! function ispline evaluates the cubic spline interpolation at point z ! ispline = y(i)+b(i)*(u-x(i))+c(i)*(u-x(i))**2+d(i)*(u-x(i))**3 ! where x(i) <= u <= x(i+1) !---------------------------------------------------------------------- ! input.. ! u = the abscissa at which the spline is to be evaluated ! x, y = the arrays of given data points ! b, c, d = arrays of spline coefficients computed by spline ! n = the number of data points ! output: ! ispline = interpolated value at point u !======================================================================= implicit none real ispline integer n real u, x(n), y(n), b(n), c(n), d(n) integer i, j, k real dx ! if u is ouside the x() interval take a boundary value (left or right) if(u <= x(1)) then ispline = y(1) return end if if(u >= x(n)) then ispline = y(n) return end if !* ! binary search for for i, such that x(i) <= u <= x(i+1) !* i = 1 j = n+1 do while (j > i+1) k = (i+j)/2 if(u < x(k)) then j=k else i=k end if end do !* ! evaluate spline interpolation !* dx = u - x(i) ispline = y(i) + dx*(b(i) + dx*(c(i) + dx*d(i))) end function ispline ! ------------------------------------------------------------------ end module module_stoch