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