/wrfv2_fire/dyn_em/module_stoch.F
FORTRAN Legacy | 933 lines | 645 code | 120 blank | 168 comment | 29 complexity | 0795edc22316699ab25b0b14d52c0d60 MD5 | raw file
Possible License(s): AGPL-1.0
- module module_stoch
- !***********************************************************************
- !
- ! Purpose: Stochastic kinetic-energy backscatter scheme (SKEB)
- ! Author : Judith Berner, NCAR (berner@ucar.edu)
- ! Date : Dec 2010
- ! Version: 1.0
- !
- !***********************************************************************
- !
- ! 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 2010:
- ! Model uncertainty in a mesoscale ensemble prediction system: Stochastic
- ! versus multi-physics representations, MWR, accepted
- ! (available through the AMS early online release)
- !
- ! Features:
- ! Version 1.0:
- ! 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_opt - 0) No stochastic parameterization
- ! - 1) Stochastic kinetic-energy backscatter scheme (SKEB)
- ! stoch_vertstruc_opt - 0) Constant vertical structure
- ! - 1) Random phase vertical structure
- ! tot_backscat_psi - Strength of streamfunction perturbations
- ! tot_backscat-t - Strength of potential temperature perturbations
- !
- !***********************************************************************
- ! ------------------------------------------------------------------
- !************** DECLARE FIELDS AND VARIABLES FOR STOCHASTIC BACKSCATTER
- ! ------------------------------------------------------------------
- implicit none
- public :: SETUP_STOCH, UPDATE_STOCH,do_fftback_along_x,do_fftback_along_y,&
- SP2GP_prep
- INTEGER :: LMINFORC, LMAXFORC, KMINFORC, KMAXFORC, &
- & LMINFORCT, LMAXFORCT, KMINFORCT, KMAXFORCT
- REAL :: ALPH, TOT_BACKSCAT_PSI, TOT_BACKSCAT_T, REXPONENT
- ! ----------Fields for spectral transform -----------
- INTEGER :: LENSAV
- INTEGER,ALLOCATABLE:: wavenumber_k(:), wavenumber_l(:)
- REAL, ALLOCATABLE :: WSAVE1(:),WSAVE2(:)
- ! --------- Others -------------------------------------------------
- REAL, PARAMETER:: RPI= 3.141592653589793 !4.0*atan(1.0)
- REAL, PARAMETER:: CP= 1006 ! specific heat of dry air in J/(Kg*K)= m^2/(K* s^2)
- save
- !=======================================================================
- contains
- !=======================================================================
- ! ------------------------------------------------------------------
- !!******** INITIALIZE STOCHASTIC KINETIC ENERGY BACKSCATTER (SKEB) *****
- ! ------------------------------------------------------------------
- subroutine SETUP_STOCH( &
- VERTSTRUCC,VERTSTRUCS, &
- SPT_AMP,SPSTREAM_AMP, &
- stoch_vertstruc_opt, &
- ISEED1,ISEED2,itime_step,DX,DY, &
- TOT_BACKSCAT_PSI,TOT_BACKSCAT_T, &
- ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte )
- IMPLICIT NONE
- INTEGER :: IER,IK,IL,iseed1,iseed2,I,J
- INTEGER :: itime_step,stoch_vertstruc_opt
- INTEGER :: KMAX,LMAX,LENSAV,ILEV
- INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte
- REAL :: DX,DY,RY,RX,RATIO_BACKSCAT,TOT_BACKSCAT_PSI,TOT_BACKSCAT_T
- REAL :: ZGAMMAN,ZTAU,ZCONSTF0,ZCONSTF0T,ZSIGMA2_EPS, RHOKLMAX,ZREF,RHOKL,EPS
- REAL, DIMENSION (ims:ime,kms:kme,jms:jme) :: VERTSTRUCC,VERTSTRUCS
- REAL, DIMENSION (ims:ime,jms:jme) :: SPSTREAM_AMP,SPT_AMP
- REAL, DIMENSION (ids:ide,jds:jde) :: ZCHI,ZCHIT
- LOGICAL :: is_print = .true.
- INTEGER , ALLOCATABLE , DIMENSION(:) :: iseed
- INTEGER :: how_many
- LOGICAL , EXTERNAL :: wrf_dm_on_monitor
-
- ! --------- SETUP PARAMETERS ---------------------------------------
- KMAX=(jde-jds)+1 !NLAT
- LMAX=(ide-ids)+1 !NLON
- RY= KMAX*DY
- RX= LMAX*DY
- LENSAV= 4*(KMAX+LMAX)+INT(LOG(REAL(KMAX))) + INT(LOG(REAL(LMAX))) + 8
- ! --------- ALLOCATE FIELDS FOR FFTPACK----------------------------
- ALLOCATE(WSAVE1(LENSAV),WSAVE2(LENSAV))
- ALLOCATE (wavenumber_k(jds:jde),wavenumber_l(ids:ide))
- ! -------- 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
- 95 format('error in cFFT2I= 'i5)
- call findindex( wavenumber_k, wavenumber_l, &
- ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte )
- ! ---------- INITIAIZE STOCHASTIC KINETIC ENERGY BACKSCATTER PARAMETERS-----------
- REXPONENT=-1.83 !produces 2(p+1) kinetic energy spectra % p=-11/6=1.83 => k=-5/3
- ! TOT_BACKSCAT_PSI = 2.0
- ! TOT_BACKSCAT_T = 4.8E-04 ! 2.E-06/240
- KMINFORC=0
- KMAXFORC=min0(40,KMAX/2)
- LMINFORC=KMINFORC
- LMAXFORC=KMAXFORC
- KMINFORCT=0
- KMAXFORCT=KMAXFORC
- LMINFORCT=KMINFORCT
- LMAXFORCT=KMAXFORCT
- ZTAU = 2.E04/12.
- ALPH = float(itime_step)/ZTAU ! approximation of 1.-exp(-itime_step/ZTAU)
- ZSIGMA2_EPS=1./(12.0*ALPH)
- ! Sanity checks for forcing wavenumber range
- IF (LMAXFORC>LMAX/2) then
- LMAXFORC=min0(40,LMAX/2)-1
- KMAXFORC=LMAXFORC
- ENDIF
- IF (LMAXFORCT>LMAX/2) then
- LMAXFORCT=min0(40,LMAX/2)-1
- KMAXFORCT=LMAXFORCT
- ENDIF
- IF ((LMINFORC>LMAXFORC).or.(KMINFORC>KMAXFORC)) then
- WRITE(*,'('' LMINFORC>LMAXFORC IN SETUP_STOCH.F90'')')
- STOP
- ENDIF
- IF ((KMAXFORC>KMAX/2).or.(LMAXFORC>LMAX/2)) then
- WRITE(*,'('' KMAXFORC>KMAX/2 IN SETUP_STOCH.F90'')')
- print*,KMAXFORC,KMAX/2
- STOP
- ENDIF
- IF ((KMINFORC.ne.LMINFORC).or.(KMAXFORC.ne.LMAXFORC)) then
- WRITE(*,'('' Forcing is non-homogenious in latitude and longitude'')')
- WRITE(*,'('' If this is what you want, comment *stop* IN SETUP_STOCH.F90'')')
- STOP
- ENDIF
- ! Output of stochastic settings
- if (is_print) then
- WRITE(*,'('' '')')
- WRITE(*,'('' =============================================='')')
- WRITE(*,'('' >> Initializing stochastic kinetic-energy backscatter scheme << '')')
- WRITE(*,'('' Total backscattered energy, TOT_BACKSCAT_PSI '',E12.5)') TOT_BACKSCAT_PSI
- WRITE(*,'('' Total backscattered temperature, TOT_BACKSCAT_T '',E12.5)') TOT_BACKSCAT_T
- WRITE(*,'('' Exponent for energy spectra, REXPONENT ='',E12.5)') REXPONENT
- WRITE(*,'('' Minimal wavenumber of streamfunction forcing, LMINFORC ='',I10)') LMINFORC
- WRITE(*,'('' Maximal wavenumber of streamfunction forcing, LMAXFORC ='',I10)') LMAXFORC
- WRITE(*,'('' Minimal wavenumber of streamfunction forcing, KMINFORC ='',I10)') KMINFORC
- WRITE(*,'('' Maximal wavenumber of streamfunction forcing, KMAXFORC ='',I10)') KMAXFORC
- WRITE(*,'('' Minimal wavenumber of temperature forcing, LMINFORCT ='',I10)') LMINFORCT
- WRITE(*,'('' Maximal wavenumber of temperature forcing, LMAXFORCT ='',I10)') LMAXFORCT
- WRITE(*,'('' stoch_vertstruc_opt '',I10)') stoch_vertstruc_opt
- WRITE(*,'('' Time step: itime_step='',I10)') itime_step
- WRITE(*,'('' Decorrelation time of noise, ZTAU ='',E12.5)') ZTAU
- WRITE(*,'('' Variance of noise, ZSIGMA2_EPS ='',E12.5)') ZSIGMA2_EPS
- WRITE(*,'('' Autoregressive parameter 1-ALPH ='',E12.5)') 1.-ALPH
- WRITE(*,'('' =============================================='')')
- endif
- ! ---------- INITIALIZE NOISE AMPLITUDES ----------------------------------
- ! Amplitudes for streamfunction and temperature perturbations: SPSTREAM_AMP , SPT_AMP
- ! 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
- ! First the constants:
- ZCHI = 0.0
- ! Fill lower left quadrant of ZCHI. For this range the indeces IK,IL
- do IK=KMINFORC,KMAXFORC
- do IL=LMINFORC,LMAXFORC
- if ((sqrt(float(IK*IK+IL*IL))).le.(KMAXFORC)) then
- if ((IK>0).or.(IL>0)) then
- ZCHI(IL+1,IK+1)=((IK/RY*IK/RY)+(IL/RX*IL/RX))**(REXPONENT/2.)
- endif
- endif
- enddo
- enddo
- ZGAMMAN = 0.0
- DO IK=KMINFORC,KMAXFORC
- DO IL=LMINFORC,LMAXFORC
- if (sqrt(float(IK*IK+IL*IL)).le.KMAXFORC) then
- if ((IK>0).or.(IL>0)) then
- ZGAMMAN= ZGAMMAN + ((IK/RY*IK/RY)+(IL/RX*IL/RX))**(REXPONENT+1)
- endif
- endif
- ENDDO
- ENDDO
- ZGAMMAN=4.0*ZGAMMAN !account for all quadrants, although only one is Filled
- ! A value TOT_BACKSCAT_PSI=xx m^2/S^3 means that in each gridbox and on average,
- ! a dissipation rate of D=TOT_BACKSCAT_PSI m^2/s^3 is backscattered onto the resolved streamfunction
- ! The resulting units for ZCONSTF0 are sqrt(m^2/(s^3*s)) = m^2/s^2, which is the unit of dpsi/dt
- ! Note, that the unit of IK has the unit m here.
- ZCONSTF0=SQRT(ALPH*TOT_BACKSCAT_PSI/(float(itime_step)*ZSIGMA2_EPS*ZGAMMAN))/(2*RPI)
- ZCHIT = 0.0
- ! Fill lower left quadrant of ZCHI. For this range the indeces IK,IL
- do IK=KMINFORCT,KMAXFORCT
- do IL=LMINFORCT,LMAXFORCT
- if ((sqrt(float(IK*IK+IL*IL))).le.(KMAXFORCT)) then
- if ((IK>0).or.(IL>0)) then
- ZCHIT(IL+1,IK+1)=((IK/RY*IK/RY)+(IL/RX*IL/RX))**(REXPONENT/2.)
- endif
- endif
- enddo
- enddo
- ZGAMMAN = 0.0
- DO IK=KMINFORCT,KMAXFORCT
- DO IL=LMINFORCT,LMAXFORCT
- if (sqrt(float(IK*IK+IL*IL)).le.KMAXFORC) then
- if ((IK>0).or.(IL>0)) then
- ZGAMMAN= ZGAMMAN + ((IK/RY*IK/RY)+(IL/RX*IL/RX))**(REXPONENT+1)
- endif
- endif
- ENDDO
- ENDDO
- ZGAMMAN=4.0*ZGAMMAN !account for all quadrants, although only one is Filled
- ! A value TOT_BACKSCAT_T= xx m^2/S^3 means that in each gridbox and on average,
- ! a dissipation rate of D=TOT_BACKSCAT_T m^2/s^3 is backscattered onto the resolved temperture pattern
- ! The resulting units for ZCONSTF0T are m^2/s^3* (K* s^2)/m^2 = K/s , which is the unit of dT/dt
- ZCONSTF0T=TOT_BACKSCAT_T /cp* SQRT(ALPH/(ZSIGMA2_EPS*ZGAMMAN))/(2*RPI)
- ! 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
- SPSTREAM_AMP=0.0
- SPT_AMP=0.0
- SPT_AMP=0.0
- DO IK=jts,jte
- DO IL=its,ite
- if ((IL .le. (LMAX/2+1)) .and. (IK .le. (KMAX/2+1)) ) then
- SPSTREAM_AMP(IL,IK) = ZCONSTF0 *ZCHI(IL,IK)
- SPT_AMP(IL,IK) = ZCONSTF0T*ZCHIT(IL,IK)
- endif
- ENDDO
- ENDDO
- ! Fill other quadrants:
- ! Upper left quadrant
- DO IK=jts,jte
- DO IL=its,ite
- if ( (IL .gt. (LMAX/2+1)) .and. (IK .le. (KMAX/2+1)) ) then
- SPSTREAM_AMP(IL,IK) = ZCONSTF0 *ZCHI(LMAX-IL+2,IK)
- SPT_AMP(IL,IK) = ZCONSTF0T*ZCHIT(LMAX-IL+2,IK)
- endif
- ENDDO
- ENDDO
- ! Lower right quadrant
- DO IK=jts,jte
- DO IL=its,ite
- if ((IK .gt. (KMAX/2+1)) .and. (IL.le.LMAX/2) ) then
- SPSTREAM_AMP(IL,IK) = ZCONSTF0 *ZCHI(IL,KMAX-IK+2)
- SPT_AMP(IL,IK) = ZCONSTF0T*ZCHIT(IL,KMAX-IK+2)
- endif
- ENDDO
- ENDDO
- ! Upper right quadrant
- DO IK=jts,jte
- DO IL=its,ite
- if ((IK .gt. (KMAX/2+1)) .and. (IL.gt.LMAX/2) ) then
- SPSTREAM_AMP(IL,IK) = ZCONSTF0 *ZCHI(LMAX-IL+2,KMAX-IK+2)
- SPT_AMP(IL,IK) = ZCONSTF0T*ZCHIT(LMAX-IL+2,KMAX-IK+2)
- endif
- ENDDO
- ENDDO
- ! Array for vertical structure if desired
- IF (stoch_vertstruc_opt>0) 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
- END IF
-
- ! Set seed for random number generator
- CALL random_seed(size=how_many)
- ALLOCATE(iseed(how_many))
- IF ( wrf_dm_on_monitor() ) THEN
- iseed=0
- iseed(1) = iseed1
- iseed(2) = iseed2
- call random_seed(put=iseed(1:how_many)) ! set random seed on monitor.
- call random_seed(get=iseed(1:how_many))
- END IF
- #ifdef DM_PARALLEL
- CALL wrf_dm_bcast_integer ( iseed , how_many )
- CALL random_seed(put=iseed(1:how_many)) ! set random seed on each proc
- #endif
- END subroutine SETUP_STOCH
- ! ------------------------------------------------------------------
- !!************** UPDATE STOCHASTIC PATTERN IN WAVENUMBER SPACE**********
- ! ------------------------------------------------------------------
- subroutine UPDATE_STOCH( &
- SPSTREAMFORCS,SPSTREAMFORCC,SPTFORCS,SPTFORCC, &
- SPT_AMP,SPSTREAM_AMP, &
- ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte )
- IMPLICIT NONE
- REAL, DIMENSION( ids:ide,jds:jde) :: ZRANDNOSS1,ZRANDNOSC1,ZRANDNOSS2,ZRANDNOSC2
- REAL, DIMENSION (ims:ime,jms:jme) :: SPSTREAMFORCS,SPSTREAMFORCC,SPTFORCS,SPTFORCc,SPSTREAM_AMP,SPT_AMP
- INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte
-
- REAL :: Z
- INTEGER ::IL, IK,LMAX,KMAX
- LOGICAL :: LGAUSS
- KMAX=(jde-jds)+1 !NLAT
- LMAX=(ide-ids)+1 !NATX
- ! Pick the distribution of the noise
- ! Random noise uses global indexes to ensure necessary symmetries and anti-symmetries
- ! of random forcing when run on multiple processors
- LGAUSS=.false.
- IF (LGAUSS) then
- DO IK=jds,jde
- DO IL=ids,ide
- call gauss_noise(z)
- ZRANDNOSS1(IL,IK)=z
- call gauss_noise(z)
- ZRANDNOSC1(IL,IK)=z
- call gauss_noise(z)
- ZRANDNOSS2(IL,IK)=z
- call gauss_noise(z)
- ZRANDNOSC2(IL,IK)=z
- ENDDO
- ENDDO
- ELSE
- DO IK=jds,jde
- DO IL=ids,ide
- CALL RANDOM_NUMBER(z)
- ZRANDNOSS1(IL,IK)=z-0.5
- CALL RANDOM_NUMBER(z)
- ZRANDNOSC1(IL,IK)=z-0.5
- CALL RANDOM_NUMBER(z)
- ZRANDNOSS2(IL,IK)=z-0.5
- CALL RANDOM_NUMBER(z)
- ZRANDNOSC2(IL,IK)=z-0.5
- ENDDO
- ENDDO
- ENDIF
- ! Note: There are symmetries and anti-symmetries to ensure real-valued back transforms
- DO IK=jts,jte
- if (IK.le.(KMAX/2)) then
- DO IL=its,ite
- SPSTREAMFORCC(IL,IK) = (1.-ALPH)*SPSTREAMFORCC(IL,IK) + SPSTREAM_AMP(IL,IK)*(ZRANDNOSC1(IL,IK))
- SPSTREAMFORCS(IL,IK) = (1.-ALPH)*SPSTREAMFORCS(IL,IK) + SPSTREAM_AMP(IL,IK)*(ZRANDNOSS1(IL,IK))
- SPTFORCC(IL,IK) = (1.-ALPH)*SPTFORCC(IL,IK) + SPT_AMP(IL,IK) *(ZRANDNOSC2(IL,IK))
- SPTFORCS(IL,IK) = (1.-ALPH)*SPTFORCS(IL,IK) + SPT_AMP(IL,IK) *(ZRANDNOSS2(IL,IK))
- ENDDO
- endif
- ENDDO
- DO IK=jts,jte
- if (IK.ge.(KMAX/2+1))then
- DO IL=its,ite
- if (IL>1) then
- SPSTREAMFORCC(IL,IK)= (1.-ALPH)* SPSTREAMFORCC(IL,IK) + &
- SPSTREAM_AMP(IL,IK) * ZRANDNOSC1(LMAX-IL+2,KMAX-IK+2)
- SPSTREAMFORCS(IL,IK)= -((1.-ALPH)*(-1.0*SPSTREAMFORCS(IL,IK))+ &
- SPSTREAM_AMP(IL,IK) * ZRANDNOSS1(LMAX-IL+2,KMAX-IK+2))
- SPTFORCC(IL,IK)= (1.-ALPH)* SPTFORCC(IL,IK) + &
- SPT_AMP(IL,IK) * ZRANDNOSC2(LMAX-IL+2,KMAX-IK+2)
- SPTFORCS(IL,IK)= -((1.-ALPH)*(-1.0*SPTFORCS(IL,IK))+ &
- SPT_AMP(IL,IK) * ZRANDNOSS2(LMAX-IL+2,KMAX-IK+2))
- else
- SPSTREAMFORCC(1,IK) = (1.-ALPH) * SPSTREAMFORCC(1,IK) + &
- SPSTREAM_AMP(1,IK) * ZRANDNOSC1(1,KMAX-IK+2)
- SPSTREAMFORCS(1,IK) = -((1.-ALPH)*(-1.0*SPSTREAMFORCS(1,IK))+ &
- SPSTREAM_AMP(1,IK) * ZRANDNOSS1(1,KMAX-IK+2))
- SPTFORCC(1,IK) = (1.-ALPH) * SPTFORCC(1,IK) + &
- SPT_AMP(1,IK) * ZRANDNOSC2(1,KMAX-IK+2)
- SPTFORCS(1,IK) = -((1.-ALPH)*(-1.0*SPTFORCS(1,IK))+ &
- SPT_AMP(1,IK) * ZRANDNOSS2(1,KMAX-IK+2))
- endif
- ENDDO
- endif
- ENDDO
-
- END subroutine UPDATE_STOCH
- ! ------------------------------------------------------------------
- !************** ADD STOCHASTIC TENDENCIES TO PHYSICAL TENDENCIES********
- ! ------------------------------------------------------------------
- SUBROUTINE CALCULATE_STOCH_TEN( &
- ru_tendf,rv_tendf,t_tendf, &
- GPUFORC,GPVFORC,GPTFORC, &
- ru_real,rv_real,rt_real, &
- mu,mub, &
- ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte, &
- dt)
- IMPLICIT NONE
- INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte
- REAL , DIMENSION(ims:ime , kms:kme, jms:jme),INTENT(INOUT) :: &
- ru_tendf, rv_tendf, t_tendf, &
- GPUFORC,GPVFORC,GPTFORC
- REAL , DIMENSION(ims:ime , kms:kme, jms:jme),INTENT(IN) :: &
- ru_real,rv_real,rt_real
- REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: mu,mub
- INTEGER :: I,J,K
- REAL :: dt
- DO j = jts,MIN(jde-1,jte)
- DO k = kts,kte-1
- DO i = its,ite
- GPUFORC(i,k,j)= ru_real(i,k,j)
- ENDDO
- ENDDO
- ENDDO
- DO j = jts,jte
- DO k = kts,kte-1
- DO i = its,MIN(ide-1,ite)
- GPVFORC(i,k,j)= rv_real(i,k,j)
- ENDDO
- ENDDO
- ENDDO
-
- DO j = jts,MIN(jde-1,jte)
- DO k = kts,kte-1
- DO i = its,MIN(ide-1,ite)
- GPTFORC(i,k,j)= rt_real(i,k,j)
- ENDDO
- ENDDO
- ENDDO
- END SUBROUTINE CALCULATE_STOCH_TEN
- ! ------------------------------------------------------------------
- SUBROUTINE UPDATE_STOCH_TEN(ru_tendf,rv_tendf,t_tendf, &
- GPUFORC,GPVFORC,GPTFORC, &
- mu,mub, &
- ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte, &
- dt )
- IMPLICIT NONE
- INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte
- REAL , DIMENSION(ims:ime , kms:kme, jms:jme),INTENT(INOUT) :: &
- ru_tendf, rv_tendf, t_tendf
- !REAL , DIMENSION(ims:ime , kms:kme, jms:jme),INTENT(IN) :: &
- REAL , DIMENSION(ims:ime , kms:kme, jms:jme) :: &
- GPUFORC,GPVFORC,GPTFORC
- REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: mu,mub
- INTEGER :: I,J,K
- REAL :: dt,xm
- DO j = jts,MIN(jde-1,jte)
- DO k = kts,kte-1
- DO i = its,ite
- ru_tendf(i,k,j) = ru_tendf(i,k,j) + GPUFORC(i,k,j) * (mu(i,j)+mub(i,j))
- ENDDO
- ENDDO
- ENDDO
- DO j = jts,jte
- DO i = its,MIN(ide-1,ite)
- DO k = kts,kte-1
- rv_tendf(i,k,j) = rv_tendf(i,k,j) + GPVFORC(i,k,j) * (mu(i,j)+mub(i,j))
- ENDDO
- ENDDO
- ENDDO
- DO j = jts,MIN(jde-1,jte)
- DO k = kts,kte-1
- DO i = its,MIN(ide-1,ite)
- t_tendf(i,k,j) = t_tendf(i,k,j) + GPTFORC(i,k,j) * (mu(i,j)+mub(i,j))
- ENDDO
- ENDDO
- ENDDO
- END SUBROUTINE UPDATE_STOCH_TEN
- ! ------------------------------------------------------------------
- !!************** TRANSFORM FROM SPHERICAL HARMONICS TO GRIDPOILT SPACE**
- ! ------------------------------------------------------------------
- subroutine SP2GP_prep( &
- SPSTREAMFORCS,SPSTREAMFORCC,SPTFORCS,SPTFORCC, &
- VERTSTRUCC,VERTSTRUCS, &
- RU_REAL,RV_REAL,RT_REAL, &
- RU_IMAG,RV_IMAG,RT_IMAG, &
- dx,dy,stoch_vertstruc_opt, &
- ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte )
- IMPLICIT NONE
- INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte
- REAL, DIMENSION (ims:ime , jms:jme) :: SPSTREAMFORCS,SPSTREAMFORCC,SPTFORCS,SPTFORCC
- REAL, DIMENSION (ims:ime , kms:kme, jms:jme) :: RU_REAL,RV_REAL,RT_REAL,RU_IMAG,RV_IMAG,RT_IMAG, &
- VERTSTRUCC,VERTSTRUCS
- INTEGER :: IK,IL,ILEV,NLAT,NLON,stoch_vertstruc_opt
- REAL :: dx,dy,RY,RX
- NLAT=(jde-jds)+1 !KMAX
- NLON=(ide-ids)+1 !LMAX
- RY= NLAT*DY
- RX= NLON*DX
- DO ILEV=kts,kte
- if (stoch_vertstruc_opt==0) then
- DO IL=its,ite
- DO IK=jts,jte
- rt_real(IL,ILEV,IK) = SPTFORCC(IL,IK)
- rt_imag(IL,ILEV,IK) = SPTFORCS(IL,IK)
- ru_real(IL,ILEV,IK) = 2*RPI/RY* wavenumber_k(IK) * SPSTREAMFORCS(IL,IK)
- ru_imag(IL,ILEV,IK) =-2*RPI/RY* wavenumber_k(IK) * SPSTREAMFORCC(IL,IK)
- rv_real(IL,ILEV,IK) =-2*RPI/RX* wavenumber_l(IL) * SPSTREAMFORCS(IL,IK)
- rv_imag(IL,ILEV,IK) = 2*RPI/RX* wavenumber_l(IL) * SPSTREAMFORCC(IL,IK)
- ENDDO
- ENDDO
- elseif (stoch_vertstruc_opt==1) then
-
- DO IL=its,ite
- DO IK=jts,jte
- rt_real(IL,ILEV,IK) = SPTFORCC(IL,IK)*VERTSTRUCC(IL,ILEV,IK) - SPTFORCS(IL,IK)*VERTSTRUCS(IL,ILEV,IK)
- rt_imag(IL,ILEV,IK) = SPTFORCC(IL,IK)*VERTSTRUCS(IL,ILEV,IK) + SPTFORCS(IL,IK)*VERTSTRUCC(IL,ILEV,IK)
- ru_real(IL,ILEV,IK) = 2*RPI/RY* wavenumber_k(IK) *&
- (+SPSTREAMFORCC(IL,IK)*VERTSTRUCS(IL,ILEV,IK) + SPSTREAMFORCS(IL,IK)*VERTSTRUCC(IL,ILEV,IK))
- ru_imag(IL,ILEV,IK) = 2*RPI/RY* wavenumber_k(IK) *&
- (-SPSTREAMFORCC(IL,IK)*VERTSTRUCC(IL,ILEV,IK) + SPSTREAMFORCS(IL,IK)*VERTSTRUCS(IL,ILEV,IK))
- rv_real(IL,ILEV,IK) = 2*RPI/RX* wavenumber_l(IL) *&
- (-SPSTREAMFORCC(IL,IK)*VERTSTRUCS(IL,ILEV,IK) - SPSTREAMFORCS(IL,IK)*VERTSTRUCC(IL,ILEV,IK))
- rv_imag(IL,ILEV,IK) = 2*RPI/RX* wavenumber_l(IL) *&
- (+SPSTREAMFORCC(IL,IK)*VERTSTRUCC(IL,ILEV,IK) - SPSTREAMFORCS(IL,IK)*VERTSTRUCS(IL,ILEV,IK))
-
- ENDDO
- ENDDO
- endif
- ENDDO !ILEV
-
- END subroutine SP2GP_prep
- ! ------------------------------------------------------------------
- !!************** SUBROUTINE DO_FFTBACK_ALONG_X
- ! ------------------------------------------------------------------
- subroutine do_fftback_along_x(fieldc_U_xxx,fields_U_xxx, &
- fieldc_V_xxx,fields_V_xxx, &
- fieldc_T_xxx,fields_T_xxx, &
- ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- ips, ipe, jps, jpe, kps, kpe, &
- imsx,imex,jmsx,jmex,kmsx,kmex, &
- ipsx,ipex,jpsx,jpex,kpsx,kpex, &
- imsy,imey,jmsy,jmey,kmsy,kmey, &
- ipsy,ipey,jpsy,jpey,kpsy,kpey, &
- k_start , k_end &
- )
- IMPLICIT NONE
-
- INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- ips, ipe, jps, jpe, kps, kpe, &
- imsx,imex,jmsx,jmex,kmsx,kmex, &
- ipsx,ipex,jpsx,jpex,kpsx,kpex, &
- imsy,imey,jmsy,jmey,kmsy,kmey, &
- ipsy,ipey,jpsy,jpey,kpsy,kpey, &
- k_start , k_end
-
- REAL, DIMENSION (imsx:imex, kmsx:kmex, jmsx:jmex) :: fieldc_U_xxx,fields_U_xxx, &
- fieldc_V_xxx,fields_V_xxx, &
- fieldc_T_xxx,fields_T_xxx
- COMPLEX, DIMENSION (ipsx:ipex) :: dummy_complex
- INTEGER :: IER,LENWRK,KMAX,LMAX,I,J,K
- REAL, ALLOCATABLE :: WORK(:)
- CHARACTER (LEN=160) :: mess
- KMAX=(jde-jds)+1
- LMAX=(ide-ids)+1
- LENWRK=2*KMAX*LMAX
- ALLOCATE(WORK(LENWRK))
- LENSAV= 4*(KMAX+LMAX)+INT(LOG(REAL(KMAX))) + INT(LOG(REAL(LMAX))) + 8
- DO k=kpsx,kpex
- DO j = jpsx, jpex
- DO i = ipsx, ipex
- dummy_complex(i)=cmplx(fieldc_U_xxx(i,k,j),fields_U_xxx(i,k,j))
- ENDDO
- CALL cFFT1B (LMAX, 1 ,dummy_complex,LMAX, WSAVE1, 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 = ipsx, ipex
- fieldc_U_xxx(i,k,j)=real(dummy_complex(i))
- fields_U_xxx(i,k,j)=imag(dummy_complex(i))
- END DO
- END DO
- END DO
- DO k=kpsx,kpex
- DO j = jpsx, jpex
- DO i = ipsx, ipex
- dummy_complex(i)=cmplx(fieldc_V_xxx(i,k,j),fields_V_xxx(i,k,j))
- ENDDO
- CALL cFFT1B (LMAX, 1 ,dummy_complex,LMAX, WSAVE1, LENSAV, WORK, LENWRK, IER)
- if (ier.ne.0) then
- WRITE(mess,FMT='(A)') 'error in cFFT1B in do_fftback_along_x, field V'
- CALL wrf_debug(0,mess)
- end if
- DO i = ipsx,ipex
- fieldc_V_xxx(i,k,j)=real(dummy_complex(i))
- fields_V_xxx(i,k,j)=imag(dummy_complex(i))
- END DO
- END DO
- END DO
- DO k=kpsx,kpex
- DO j = jpsx, jpex
- DO i = ipsx, ipex
- dummy_complex(i)=cmplx(fieldc_T_xxx(i,k,j),fields_T_xxx(i,k,j))
- ENDDO
- CALL cFFT1B (LMAX, 1 ,dummy_complex,LMAX, WSAVE1, LENSAV, WORK, LENWRK, IER)
- if (ier.ne.0) then
- WRITE(mess,FMT='(A)') 'error in cFFT1B in do_fftback_along_x, field T'
- CALL wrf_debug(0,mess)
- end if
- DO i = ipsx, ipex
- fieldc_T_xxx(i,k,j)=real(dummy_complex(i))
- fields_T_xxx(i,k,j)=imag(dummy_complex(i))
- END DO
- END DO
- END DO !
- DEALLOCATE(WORK)
- end subroutine do_fftback_along_x
- !! ------------------------------------------------------------------
- !!!************** SUBROUTINE DO_FFTBACK_ALONG_Y
- !! ------------------------------------------------------------------
- subroutine do_fftback_along_y(fieldc_U_yyy,fields_U_yyy, &
- fieldc_V_yyy,fields_V_yyy, &
- fieldc_T_yyy,fields_T_yyy, &
- ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- ips, ipe, jps, jpe, kps, kpe, &
- imsx,imex,jmsx,jmex,kmsx,kmex, &
- ipsx,ipex,jpsx,jpex,kpsx,kpex, &
- imsy,imey,jmsy,jmey,kmsy,kmey, &
- ipsy,ipey,jpsy,jpey,kpsy,kpey, &
- k_start , k_end &
- )
- IMPLICIT NONE
-
- INTEGER :: IER,LENWRK,KMAX,LMAX,I,J,K
-
- INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- ips, ipe, jps, jpe, kps, kpe, &
- imsx,imex,jmsx,jmex,kmsx,kmex, &
- ipsx,ipex,jpsx,jpex,kpsx,kpex, &
- imsy,imey,jmsy,jmey,kmsy,kmey, &
- ipsy,ipey,jpsy,jpey,kpsy,kpey, &
- k_start , k_end
-
- REAL, DIMENSION (imsy:imey, kmsy:kmey, jmsy:jmey) :: fieldc_U_yyy,fields_U_yyy, &
- fieldc_V_yyy,fields_V_yyy, &
- fieldc_T_yyy,fields_T_yyy
- COMPLEX, DIMENSION (jpsy:jpey) :: dummy_complex
- REAL, ALLOCATABLE :: WORK(:)
- CHARACTER (LEN=160) :: mess
- KMAX=(jde-jds)+1
- LMAX=(ide-ids)+1
- LENWRK=2*KMAX*LMAX
- ALLOCATE(WORK(LENWRK))
- LENSAV= 4*(KMAX+LMAX)+INT(LOG(REAL(KMAX))) + INT(LOG(REAL(LMAX))) + 8
- DO k=kpsy,kpey
- DO i = ipsy, ipey
- DO j = jpsy,jpey
- dummy_complex(j)=cmplx(fieldc_U_yyy(i,k,j),fields_U_yyy(i,k,j))
- ENDDO
- CALL cFFT1B (KMAX, 1 ,dummy_complex,KMAX, WSAVE2, 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 = jpsy, jpey
- fieldc_U_yyy(i,k,j)=real(dummy_complex(j))
- fields_U_yyy(i,k,j)=imag(dummy_complex(j))
- END DO
- END DO
- END DO ! k_start-k_end
- DO k=kpsy,kpey
- DO i = ipsy, ipey
- DO j = jpsy, jpey
- dummy_complex(j)=cmplx(fieldc_V_yyy(i,k,j),fields_V_yyy(i,k,j))
- ENDDO
- CALL cFFT1B (KMAX, 1 ,dummy_complex,KMAX, WSAVE2, LENSAV, WORK, LENWRK, IER)
- if (ier.ne.0) then
- WRITE(mess,FMT='(A)') 'error in cFFT1B in do_fftback_along_y, field V'
- CALL wrf_debug(0,mess)
- end if
- DO j = jpsy, jpey
- fieldc_V_yyy(i,k,j)=real(dummy_complex(j))
- fields_V_yyy(i,k,j)=imag(dummy_complex(j))
- END DO
- END DO
- END DO ! k_start-k_end
- DO k=kpsy,kpey
- DO i = ipsy, ipey
- DO j = jpsy,jpey
- dummy_complex(j)=cmplx(fieldc_T_yyy(i,k,j),fields_T_yyy(i,k,j))
- ENDDO
- CALL cFFT1B (KMAX, 1 ,dummy_complex,KMAX, WSAVE2, LENSAV, WORK, LENWRK, IER)
- if (ier.ne.0) then
- WRITE(mess,FMT='(A)') 'error in cFFT1B in do_fftback_along_y, field T'
- CALL wrf_debug(0,mess)
- end if
- DO j = jpsy,jpey
- fieldc_T_yyy(i,k,j)=real(dummy_complex(j))
- fields_T_yyy(i,k,j)=imag(dummy_complex(j))
- END DO
- END DO
- END DO ! k_start-k_end
- DEALLOCATE(WORK)
- end subroutine do_fftback_along_y
- ! ------------------------------------------------------------------
- !!************** TRANSFORM FROM GRIDPOILT SPACE TO SPHERICAL HARMONICS **
- ! ------------------------------------------------------------------
- subroutine findindex( wavenumber_k, wavenumber_L, &
- ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte )
- IMPLICIT NONE
- INTEGER :: IK,IL,KMAX,LMAX
- INTEGER, DIMENSION (jds:jde):: wavenumber_k
- INTEGER, DIMENSION (ids:ide):: wavenumber_l
- INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte
- KMAX=(jde-jds)+1
- LMAX=(ide-ids)+1
- !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, seed1, seed2,nens )
- USE module_configure
- IMPLICIT NONE
- !
- ! Structure that contains run-time configuration (namelist) data for domain
- TYPE (grid_config_rec_type) :: config_flags
- !
- ! Arguments
- INTEGER, INTENT(OUT) :: seed1, seed2
- INTEGER, INTENT(IN ) :: nens
- ! Local
- integer :: date_time(8)
- integer*8 :: yyyy,mmdd,newtime
- integer*8 :: ihr,isc,idiv
- integer*8 :: iphys
- character (len=10) :: real_clock(3), time
- !
- LOGICAL :: is_print = .false.
- !
- ! Read CPU time
- call date_and_time(real_clock(1), real_clock(2),&
- real_clock(3), date_time)
- read(real_clock(1),'(i4)') yyyy
- read(real_clock(1),'(4x,i4)') mmdd
- read(real_clock(2),'(i6,1x,i3)') ihr,isc ! ihr = hhmmss , isc - milliseconds of the minute
- newtime = yyyy*mmdd+ihr+isc
- if (is_print) then
- print*,'real_clock: ',real_clock
- print*,'real_clock(1): ',real_clock(1)
- print*,'real_clock(2): ',real_clock(2)
- print*,'real_clock(3): ',real_clock(3)
- print*,'date_time: ',date_time
- print*,'newtime: ',newtime
- !
- ! get a seed number (w/ physics options for each ensemble member)
- !
- print *,'physics_options:',&
- config_flags%sf_surface_physics,&
- config_flags%ra_lw_physics, &
- config_flags%ra_sw_physics, &
- config_flags%mp_physics, &
- config_flags%sf_sfclay_physics,&
- config_flags%bl_pbl_physics,&
- config_flags%cu_physics
- endif !isprint
- iphys = config_flags%sf_surface_physics * 1000000 + &
- config_flags%ra_lw_physics * 100000 + &
- config_flags%ra_sw_physics * 10000 + &
- config_flags%mp_physics * 1000 + &
- config_flags%sf_sfclay_physics * 100 + &
- config_flags%bl_pbl_physics * 10 + &
- config_flags%cu_physics
-
-
- idiv=2;
- seed1 = newtime+iphys+nens*1000000
- seed2 = mod(newtime+iphys+nens*1000000,idiv)
- if(is_print) print *,'Rand_seed (iphys/newtime/idiv):',iphys,newtime,idiv,nens
- end SUBROUTINE rand_seed
- end module module_stoch