PageRenderTime 50ms CodeModel.GetById 13ms RepoModel.GetById 1ms app.codeStats 0ms

/wrfv2_fire/dyn_em/module_stoch.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 933 lines | 645 code | 120 blank | 168 comment | 29 complexity | 0795edc22316699ab25b0b14d52c0d60 MD5 | raw file
Possible License(s): AGPL-1.0
  1. module module_stoch
  2. !***********************************************************************
  3. !
  4. ! Purpose: Stochastic kinetic-energy backscatter scheme (SKEB)
  5. ! Author : Judith Berner, NCAR (berner@ucar.edu)
  6. ! Date : Dec 2010
  7. ! Version: 1.0
  8. !
  9. !***********************************************************************
  10. !
  11. ! The scheme introduces stochastic perturbations to the rotational wind
  12. ! components and to the potential temperature field. The stochastic
  13. ! perturbations are generated by independent autoregressive processes for
  14. ! each wavenumber and results in smooth spatially and temporally correlated patterns.
  15. ! Details of the scheme and its performance in a meso-scale WRF-ensemble
  16. ! system are available in:
  17. !
  18. ! Berner, J., S.-Y. Ha, J. P. Hacker, A. Fournier and C. Snyder 2010:
  19. ! Model uncertainty in a mesoscale ensemble prediction system: Stochastic
  20. ! versus multi-physics representations, MWR, accepted
  21. ! (available through the AMS early online release)
  22. !
  23. ! Features:
  24. ! Version 1.0:
  25. ! Dissipation: Dissipation rates are assumed constant in space and time
  26. ! Vertical structure: Supports two options for vertical structure:
  27. ! 0) constant
  28. ! 1) random phase
  29. !
  30. ! Optional namelist parameters:
  31. ! stoch_opt - 0) No stochastic parameterization
  32. ! - 1) Stochastic kinetic-energy backscatter scheme (SKEB)
  33. ! stoch_vertstruc_opt - 0) Constant vertical structure
  34. ! - 1) Random phase vertical structure
  35. ! tot_backscat_psi - Strength of streamfunction perturbations
  36. ! tot_backscat-t - Strength of potential temperature perturbations
  37. !
  38. !***********************************************************************
  39. ! ------------------------------------------------------------------
  40. !************** DECLARE FIELDS AND VARIABLES FOR STOCHASTIC BACKSCATTER
  41. ! ------------------------------------------------------------------
  42. implicit none
  43. public :: SETUP_STOCH, UPDATE_STOCH,do_fftback_along_x,do_fftback_along_y,&
  44. SP2GP_prep
  45. INTEGER :: LMINFORC, LMAXFORC, KMINFORC, KMAXFORC, &
  46. & LMINFORCT, LMAXFORCT, KMINFORCT, KMAXFORCT
  47. REAL :: ALPH, TOT_BACKSCAT_PSI, TOT_BACKSCAT_T, REXPONENT
  48. ! ----------Fields for spectral transform -----------
  49. INTEGER :: LENSAV
  50. INTEGER,ALLOCATABLE:: wavenumber_k(:), wavenumber_l(:)
  51. REAL, ALLOCATABLE :: WSAVE1(:),WSAVE2(:)
  52. ! --------- Others -------------------------------------------------
  53. REAL, PARAMETER:: RPI= 3.141592653589793 !4.0*atan(1.0)
  54. REAL, PARAMETER:: CP= 1006 ! specific heat of dry air in J/(Kg*K)= m^2/(K* s^2)
  55. save
  56. !=======================================================================
  57. contains
  58. !=======================================================================
  59. ! ------------------------------------------------------------------
  60. !!******** INITIALIZE STOCHASTIC KINETIC ENERGY BACKSCATTER (SKEB) *****
  61. ! ------------------------------------------------------------------
  62. subroutine SETUP_STOCH( &
  63. VERTSTRUCC,VERTSTRUCS, &
  64. SPT_AMP,SPSTREAM_AMP, &
  65. stoch_vertstruc_opt, &
  66. ISEED1,ISEED2,itime_step,DX,DY, &
  67. TOT_BACKSCAT_PSI,TOT_BACKSCAT_T, &
  68. ids, ide, jds, jde, kds, kde, &
  69. ims, ime, jms, jme, kms, kme, &
  70. its, ite, jts, jte, kts, kte )
  71. IMPLICIT NONE
  72. INTEGER :: IER,IK,IL,iseed1,iseed2,I,J
  73. INTEGER :: itime_step,stoch_vertstruc_opt
  74. INTEGER :: KMAX,LMAX,LENSAV,ILEV
  75. INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
  76. ims, ime, jms, jme, kms, kme, &
  77. its, ite, jts, jte, kts, kte
  78. REAL :: DX,DY,RY,RX,RATIO_BACKSCAT,TOT_BACKSCAT_PSI,TOT_BACKSCAT_T
  79. REAL :: ZGAMMAN,ZTAU,ZCONSTF0,ZCONSTF0T,ZSIGMA2_EPS, RHOKLMAX,ZREF,RHOKL,EPS
  80. REAL, DIMENSION (ims:ime,kms:kme,jms:jme) :: VERTSTRUCC,VERTSTRUCS
  81. REAL, DIMENSION (ims:ime,jms:jme) :: SPSTREAM_AMP,SPT_AMP
  82. REAL, DIMENSION (ids:ide,jds:jde) :: ZCHI,ZCHIT
  83. LOGICAL :: is_print = .true.
  84. INTEGER , ALLOCATABLE , DIMENSION(:) :: iseed
  85. INTEGER :: how_many
  86. LOGICAL , EXTERNAL :: wrf_dm_on_monitor
  87. ! --------- SETUP PARAMETERS ---------------------------------------
  88. KMAX=(jde-jds)+1 !NLAT
  89. LMAX=(ide-ids)+1 !NLON
  90. RY= KMAX*DY
  91. RX= LMAX*DY
  92. LENSAV= 4*(KMAX+LMAX)+INT(LOG(REAL(KMAX))) + INT(LOG(REAL(LMAX))) + 8
  93. ! --------- ALLOCATE FIELDS FOR FFTPACK----------------------------
  94. ALLOCATE(WSAVE1(LENSAV),WSAVE2(LENSAV))
  95. ALLOCATE (wavenumber_k(jds:jde),wavenumber_l(ids:ide))
  96. ! -------- INITIALIZE FFTPACK ROUTINES -----------------------------
  97. call CFFT1I (LMAX, WSAVE1, LENSAV, IER)
  98. if(ier.ne. 0) write(*,95) ier
  99. call CFFT1I (KMAX, WSAVE2, LENSAV, IER)
  100. if(ier.ne. 0) write(*,95) ier
  101. 95 format('error in cFFT2I= 'i5)
  102. call findindex( wavenumber_k, wavenumber_l, &
  103. ids, ide, jds, jde, kds, kde, &
  104. ims, ime, jms, jme, kms, kme, &
  105. its, ite, jts, jte, kts, kte )
  106. ! ---------- INITIAIZE STOCHASTIC KINETIC ENERGY BACKSCATTER PARAMETERS-----------
  107. REXPONENT=-1.83 !produces 2(p+1) kinetic energy spectra % p=-11/6=1.83 => k=-5/3
  108. ! TOT_BACKSCAT_PSI = 2.0
  109. ! TOT_BACKSCAT_T = 4.8E-04 ! 2.E-06/240
  110. KMINFORC=0
  111. KMAXFORC=min0(40,KMAX/2)
  112. LMINFORC=KMINFORC
  113. LMAXFORC=KMAXFORC
  114. KMINFORCT=0
  115. KMAXFORCT=KMAXFORC
  116. LMINFORCT=KMINFORCT
  117. LMAXFORCT=KMAXFORCT
  118. ZTAU = 2.E04/12.
  119. ALPH = float(itime_step)/ZTAU ! approximation of 1.-exp(-itime_step/ZTAU)
  120. ZSIGMA2_EPS=1./(12.0*ALPH)
  121. ! Sanity checks for forcing wavenumber range
  122. IF (LMAXFORC>LMAX/2) then
  123. LMAXFORC=min0(40,LMAX/2)-1
  124. KMAXFORC=LMAXFORC
  125. ENDIF
  126. IF (LMAXFORCT>LMAX/2) then
  127. LMAXFORCT=min0(40,LMAX/2)-1
  128. KMAXFORCT=LMAXFORCT
  129. ENDIF
  130. IF ((LMINFORC>LMAXFORC).or.(KMINFORC>KMAXFORC)) then
  131. WRITE(*,'('' LMINFORC>LMAXFORC IN SETUP_STOCH.F90'')')
  132. STOP
  133. ENDIF
  134. IF ((KMAXFORC>KMAX/2).or.(LMAXFORC>LMAX/2)) then
  135. WRITE(*,'('' KMAXFORC>KMAX/2 IN SETUP_STOCH.F90'')')
  136. print*,KMAXFORC,KMAX/2
  137. STOP
  138. ENDIF
  139. IF ((KMINFORC.ne.LMINFORC).or.(KMAXFORC.ne.LMAXFORC)) then
  140. WRITE(*,'('' Forcing is non-homogenious in latitude and longitude'')')
  141. WRITE(*,'('' If this is what you want, comment *stop* IN SETUP_STOCH.F90'')')
  142. STOP
  143. ENDIF
  144. ! Output of stochastic settings
  145. if (is_print) then
  146. WRITE(*,'('' '')')
  147. WRITE(*,'('' =============================================='')')
  148. WRITE(*,'('' >> Initializing stochastic kinetic-energy backscatter scheme << '')')
  149. WRITE(*,'('' Total backscattered energy, TOT_BACKSCAT_PSI '',E12.5)') TOT_BACKSCAT_PSI
  150. WRITE(*,'('' Total backscattered temperature, TOT_BACKSCAT_T '',E12.5)') TOT_BACKSCAT_T
  151. WRITE(*,'('' Exponent for energy spectra, REXPONENT ='',E12.5)') REXPONENT
  152. WRITE(*,'('' Minimal wavenumber of streamfunction forcing, LMINFORC ='',I10)') LMINFORC
  153. WRITE(*,'('' Maximal wavenumber of streamfunction forcing, LMAXFORC ='',I10)') LMAXFORC
  154. WRITE(*,'('' Minimal wavenumber of streamfunction forcing, KMINFORC ='',I10)') KMINFORC
  155. WRITE(*,'('' Maximal wavenumber of streamfunction forcing, KMAXFORC ='',I10)') KMAXFORC
  156. WRITE(*,'('' Minimal wavenumber of temperature forcing, LMINFORCT ='',I10)') LMINFORCT
  157. WRITE(*,'('' Maximal wavenumber of temperature forcing, LMAXFORCT ='',I10)') LMAXFORCT
  158. WRITE(*,'('' stoch_vertstruc_opt '',I10)') stoch_vertstruc_opt
  159. WRITE(*,'('' Time step: itime_step='',I10)') itime_step
  160. WRITE(*,'('' Decorrelation time of noise, ZTAU ='',E12.5)') ZTAU
  161. WRITE(*,'('' Variance of noise, ZSIGMA2_EPS ='',E12.5)') ZSIGMA2_EPS
  162. WRITE(*,'('' Autoregressive parameter 1-ALPH ='',E12.5)') 1.-ALPH
  163. WRITE(*,'('' =============================================='')')
  164. endif
  165. ! ---------- INITIALIZE NOISE AMPLITUDES ----------------------------------
  166. ! Amplitudes for streamfunction and temperature perturbations: SPSTREAM_AMP , SPT_AMP
  167. ! 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
  168. ! First the constants:
  169. ZCHI = 0.0
  170. ! Fill lower left quadrant of ZCHI. For this range the indeces IK,IL
  171. do IK=KMINFORC,KMAXFORC
  172. do IL=LMINFORC,LMAXFORC
  173. if ((sqrt(float(IK*IK+IL*IL))).le.(KMAXFORC)) then
  174. if ((IK>0).or.(IL>0)) then
  175. ZCHI(IL+1,IK+1)=((IK/RY*IK/RY)+(IL/RX*IL/RX))**(REXPONENT/2.)
  176. endif
  177. endif
  178. enddo
  179. enddo
  180. ZGAMMAN = 0.0
  181. DO IK=KMINFORC,KMAXFORC
  182. DO IL=LMINFORC,LMAXFORC
  183. if (sqrt(float(IK*IK+IL*IL)).le.KMAXFORC) then
  184. if ((IK>0).or.(IL>0)) then
  185. ZGAMMAN= ZGAMMAN + ((IK/RY*IK/RY)+(IL/RX*IL/RX))**(REXPONENT+1)
  186. endif
  187. endif
  188. ENDDO
  189. ENDDO
  190. ZGAMMAN=4.0*ZGAMMAN !account for all quadrants, although only one is Filled
  191. ! A value TOT_BACKSCAT_PSI=xx m^2/S^3 means that in each gridbox and on average,
  192. ! a dissipation rate of D=TOT_BACKSCAT_PSI m^2/s^3 is backscattered onto the resolved streamfunction
  193. ! The resulting units for ZCONSTF0 are sqrt(m^2/(s^3*s)) = m^2/s^2, which is the unit of dpsi/dt
  194. ! Note, that the unit of IK has the unit m here.
  195. ZCONSTF0=SQRT(ALPH*TOT_BACKSCAT_PSI/(float(itime_step)*ZSIGMA2_EPS*ZGAMMAN))/(2*RPI)
  196. ZCHIT = 0.0
  197. ! Fill lower left quadrant of ZCHI. For this range the indeces IK,IL
  198. do IK=KMINFORCT,KMAXFORCT
  199. do IL=LMINFORCT,LMAXFORCT
  200. if ((sqrt(float(IK*IK+IL*IL))).le.(KMAXFORCT)) then
  201. if ((IK>0).or.(IL>0)) then
  202. ZCHIT(IL+1,IK+1)=((IK/RY*IK/RY)+(IL/RX*IL/RX))**(REXPONENT/2.)
  203. endif
  204. endif
  205. enddo
  206. enddo
  207. ZGAMMAN = 0.0
  208. DO IK=KMINFORCT,KMAXFORCT
  209. DO IL=LMINFORCT,LMAXFORCT
  210. if (sqrt(float(IK*IK+IL*IL)).le.KMAXFORC) then
  211. if ((IK>0).or.(IL>0)) then
  212. ZGAMMAN= ZGAMMAN + ((IK/RY*IK/RY)+(IL/RX*IL/RX))**(REXPONENT+1)
  213. endif
  214. endif
  215. ENDDO
  216. ENDDO
  217. ZGAMMAN=4.0*ZGAMMAN !account for all quadrants, although only one is Filled
  218. ! A value TOT_BACKSCAT_T= xx m^2/S^3 means that in each gridbox and on average,
  219. ! a dissipation rate of D=TOT_BACKSCAT_T m^2/s^3 is backscattered onto the resolved temperture pattern
  220. ! The resulting units for ZCONSTF0T are m^2/s^3* (K* s^2)/m^2 = K/s , which is the unit of dT/dt
  221. ZCONSTF0T=TOT_BACKSCAT_T /cp* SQRT(ALPH/(ZSIGMA2_EPS*ZGAMMAN))/(2*RPI)
  222. ! Now the wavenumber-dependent amplitudes
  223. ! Note: There are symmetries and anti-symmetries to ensure real-valued back transforms
  224. ! Fill lower left quadrant of matrix of noise amplitudes for wavenumbers K=0,KMAX/2
  225. SPSTREAM_AMP=0.0
  226. SPT_AMP=0.0
  227. SPT_AMP=0.0
  228. DO IK=jts,jte
  229. DO IL=its,ite
  230. if ((IL .le. (LMAX/2+1)) .and. (IK .le. (KMAX/2+1)) ) then
  231. SPSTREAM_AMP(IL,IK) = ZCONSTF0 *ZCHI(IL,IK)
  232. SPT_AMP(IL,IK) = ZCONSTF0T*ZCHIT(IL,IK)
  233. endif
  234. ENDDO
  235. ENDDO
  236. ! Fill other quadrants:
  237. ! Upper left quadrant
  238. DO IK=jts,jte
  239. DO IL=its,ite
  240. if ( (IL .gt. (LMAX/2+1)) .and. (IK .le. (KMAX/2+1)) ) then
  241. SPSTREAM_AMP(IL,IK) = ZCONSTF0 *ZCHI(LMAX-IL+2,IK)
  242. SPT_AMP(IL,IK) = ZCONSTF0T*ZCHIT(LMAX-IL+2,IK)
  243. endif
  244. ENDDO
  245. ENDDO
  246. ! Lower right quadrant
  247. DO IK=jts,jte
  248. DO IL=its,ite
  249. if ((IK .gt. (KMAX/2+1)) .and. (IL.le.LMAX/2) ) then
  250. SPSTREAM_AMP(IL,IK) = ZCONSTF0 *ZCHI(IL,KMAX-IK+2)
  251. SPT_AMP(IL,IK) = ZCONSTF0T*ZCHIT(IL,KMAX-IK+2)
  252. endif
  253. ENDDO
  254. ENDDO
  255. ! Upper right quadrant
  256. DO IK=jts,jte
  257. DO IL=its,ite
  258. if ((IK .gt. (KMAX/2+1)) .and. (IL.gt.LMAX/2) ) then
  259. SPSTREAM_AMP(IL,IK) = ZCONSTF0 *ZCHI(LMAX-IL+2,KMAX-IK+2)
  260. SPT_AMP(IL,IK) = ZCONSTF0T*ZCHIT(LMAX-IL+2,KMAX-IK+2)
  261. endif
  262. ENDDO
  263. ENDDO
  264. ! Array for vertical structure if desired
  265. IF (stoch_vertstruc_opt>0) then
  266. VERTSTRUCC=0.0
  267. VERTSTRUCS=0.0
  268. RHOKLMAX= sqrt(KMAX**2/DY**2 + LMAX**2/DX**2)
  269. ZREF=32.0
  270. DO ILEV=kts,kte
  271. DO IK=jts,jte
  272. DO IL=its,ite
  273. if (IL.le.(LMAX/2)) then
  274. RHOKL = sqrt((IK+1)**2/DY**2 + (IL+1)**2/DX**2)
  275. EPS = ((RHOKLMAX - RHOKL)/ RHOKLMAX) * (ILEV/ZREF) * RPI
  276. VERTSTRUCC(IL,ILEV,IK) = cos ( eps* (IL+1) )
  277. VERTSTRUCS(IL,ILEV,IK) = sin ( eps* (IL+1) )
  278. else
  279. RHOKL = sqrt((IK+1)**2/DY**2 + (LMAX-IL+2)**2/DX**2)
  280. EPS = ((RHOKLMAX - RHOKL)/ RHOKLMAX) * (ILEV/ZREF) * RPI
  281. VERTSTRUCC (IL,ILEV,IK) = cos ( eps* (LMAX-IL+2) )
  282. VERTSTRUCS (IL,ILEV,IK) = - sin ( eps* (LMAX-IL+2) )
  283. endif
  284. ENDDO
  285. ENDDO
  286. ENDDO
  287. END IF
  288. ! Set seed for random number generator
  289. CALL random_seed(size=how_many)
  290. ALLOCATE(iseed(how_many))
  291. IF ( wrf_dm_on_monitor() ) THEN
  292. iseed=0
  293. iseed(1) = iseed1
  294. iseed(2) = iseed2
  295. call random_seed(put=iseed(1:how_many)) ! set random seed on monitor.
  296. call random_seed(get=iseed(1:how_many))
  297. END IF
  298. #ifdef DM_PARALLEL
  299. CALL wrf_dm_bcast_integer ( iseed , how_many )
  300. CALL random_seed(put=iseed(1:how_many)) ! set random seed on each proc
  301. #endif
  302. END subroutine SETUP_STOCH
  303. ! ------------------------------------------------------------------
  304. !!************** UPDATE STOCHASTIC PATTERN IN WAVENUMBER SPACE**********
  305. ! ------------------------------------------------------------------
  306. subroutine UPDATE_STOCH( &
  307. SPSTREAMFORCS,SPSTREAMFORCC,SPTFORCS,SPTFORCC, &
  308. SPT_AMP,SPSTREAM_AMP, &
  309. ids, ide, jds, jde, kds, kde, &
  310. ims, ime, jms, jme, kms, kme, &
  311. its, ite, jts, jte, kts, kte )
  312. IMPLICIT NONE
  313. REAL, DIMENSION( ids:ide,jds:jde) :: ZRANDNOSS1,ZRANDNOSC1,ZRANDNOSS2,ZRANDNOSC2
  314. REAL, DIMENSION (ims:ime,jms:jme) :: SPSTREAMFORCS,SPSTREAMFORCC,SPTFORCS,SPTFORCc,SPSTREAM_AMP,SPT_AMP
  315. INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
  316. ims, ime, jms, jme, kms, kme, &
  317. its, ite, jts, jte, kts, kte
  318. REAL :: Z
  319. INTEGER ::IL, IK,LMAX,KMAX
  320. LOGICAL :: LGAUSS
  321. KMAX=(jde-jds)+1 !NLAT
  322. LMAX=(ide-ids)+1 !NATX
  323. ! Pick the distribution of the noise
  324. ! Random noise uses global indexes to ensure necessary symmetries and anti-symmetries
  325. ! of random forcing when run on multiple processors
  326. LGAUSS=.false.
  327. IF (LGAUSS) then
  328. DO IK=jds,jde
  329. DO IL=ids,ide
  330. call gauss_noise(z)
  331. ZRANDNOSS1(IL,IK)=z
  332. call gauss_noise(z)
  333. ZRANDNOSC1(IL,IK)=z
  334. call gauss_noise(z)
  335. ZRANDNOSS2(IL,IK)=z
  336. call gauss_noise(z)
  337. ZRANDNOSC2(IL,IK)=z
  338. ENDDO
  339. ENDDO
  340. ELSE
  341. DO IK=jds,jde
  342. DO IL=ids,ide
  343. CALL RANDOM_NUMBER(z)
  344. ZRANDNOSS1(IL,IK)=z-0.5
  345. CALL RANDOM_NUMBER(z)
  346. ZRANDNOSC1(IL,IK)=z-0.5
  347. CALL RANDOM_NUMBER(z)
  348. ZRANDNOSS2(IL,IK)=z-0.5
  349. CALL RANDOM_NUMBER(z)
  350. ZRANDNOSC2(IL,IK)=z-0.5
  351. ENDDO
  352. ENDDO
  353. ENDIF
  354. ! Note: There are symmetries and anti-symmetries to ensure real-valued back transforms
  355. DO IK=jts,jte
  356. if (IK.le.(KMAX/2)) then
  357. DO IL=its,ite
  358. SPSTREAMFORCC(IL,IK) = (1.-ALPH)*SPSTREAMFORCC(IL,IK) + SPSTREAM_AMP(IL,IK)*(ZRANDNOSC1(IL,IK))
  359. SPSTREAMFORCS(IL,IK) = (1.-ALPH)*SPSTREAMFORCS(IL,IK) + SPSTREAM_AMP(IL,IK)*(ZRANDNOSS1(IL,IK))
  360. SPTFORCC(IL,IK) = (1.-ALPH)*SPTFORCC(IL,IK) + SPT_AMP(IL,IK) *(ZRANDNOSC2(IL,IK))
  361. SPTFORCS(IL,IK) = (1.-ALPH)*SPTFORCS(IL,IK) + SPT_AMP(IL,IK) *(ZRANDNOSS2(IL,IK))
  362. ENDDO
  363. endif
  364. ENDDO
  365. DO IK=jts,jte
  366. if (IK.ge.(KMAX/2+1))then
  367. DO IL=its,ite
  368. if (IL>1) then
  369. SPSTREAMFORCC(IL,IK)= (1.-ALPH)* SPSTREAMFORCC(IL,IK) + &
  370. SPSTREAM_AMP(IL,IK) * ZRANDNOSC1(LMAX-IL+2,KMAX-IK+2)
  371. SPSTREAMFORCS(IL,IK)= -((1.-ALPH)*(-1.0*SPSTREAMFORCS(IL,IK))+ &
  372. SPSTREAM_AMP(IL,IK) * ZRANDNOSS1(LMAX-IL+2,KMAX-IK+2))
  373. SPTFORCC(IL,IK)= (1.-ALPH)* SPTFORCC(IL,IK) + &
  374. SPT_AMP(IL,IK) * ZRANDNOSC2(LMAX-IL+2,KMAX-IK+2)
  375. SPTFORCS(IL,IK)= -((1.-ALPH)*(-1.0*SPTFORCS(IL,IK))+ &
  376. SPT_AMP(IL,IK) * ZRANDNOSS2(LMAX-IL+2,KMAX-IK+2))
  377. else
  378. SPSTREAMFORCC(1,IK) = (1.-ALPH) * SPSTREAMFORCC(1,IK) + &
  379. SPSTREAM_AMP(1,IK) * ZRANDNOSC1(1,KMAX-IK+2)
  380. SPSTREAMFORCS(1,IK) = -((1.-ALPH)*(-1.0*SPSTREAMFORCS(1,IK))+ &
  381. SPSTREAM_AMP(1,IK) * ZRANDNOSS1(1,KMAX-IK+2))
  382. SPTFORCC(1,IK) = (1.-ALPH) * SPTFORCC(1,IK) + &
  383. SPT_AMP(1,IK) * ZRANDNOSC2(1,KMAX-IK+2)
  384. SPTFORCS(1,IK) = -((1.-ALPH)*(-1.0*SPTFORCS(1,IK))+ &
  385. SPT_AMP(1,IK) * ZRANDNOSS2(1,KMAX-IK+2))
  386. endif
  387. ENDDO
  388. endif
  389. ENDDO
  390. END subroutine UPDATE_STOCH
  391. ! ------------------------------------------------------------------
  392. !************** ADD STOCHASTIC TENDENCIES TO PHYSICAL TENDENCIES********
  393. ! ------------------------------------------------------------------
  394. SUBROUTINE CALCULATE_STOCH_TEN( &
  395. ru_tendf,rv_tendf,t_tendf, &
  396. GPUFORC,GPVFORC,GPTFORC, &
  397. ru_real,rv_real,rt_real, &
  398. mu,mub, &
  399. ids, ide, jds, jde, kds, kde, &
  400. ims, ime, jms, jme, kms, kme, &
  401. its, ite, jts, jte, kts, kte, &
  402. dt)
  403. IMPLICIT NONE
  404. INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
  405. ims, ime, jms, jme, kms, kme, &
  406. its, ite, jts, jte, kts, kte
  407. REAL , DIMENSION(ims:ime , kms:kme, jms:jme),INTENT(INOUT) :: &
  408. ru_tendf, rv_tendf, t_tendf, &
  409. GPUFORC,GPVFORC,GPTFORC
  410. REAL , DIMENSION(ims:ime , kms:kme, jms:jme),INTENT(IN) :: &
  411. ru_real,rv_real,rt_real
  412. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: mu,mub
  413. INTEGER :: I,J,K
  414. REAL :: dt
  415. DO j = jts,MIN(jde-1,jte)
  416. DO k = kts,kte-1
  417. DO i = its,ite
  418. GPUFORC(i,k,j)= ru_real(i,k,j)
  419. ENDDO
  420. ENDDO
  421. ENDDO
  422. DO j = jts,jte
  423. DO k = kts,kte-1
  424. DO i = its,MIN(ide-1,ite)
  425. GPVFORC(i,k,j)= rv_real(i,k,j)
  426. ENDDO
  427. ENDDO
  428. ENDDO
  429. DO j = jts,MIN(jde-1,jte)
  430. DO k = kts,kte-1
  431. DO i = its,MIN(ide-1,ite)
  432. GPTFORC(i,k,j)= rt_real(i,k,j)
  433. ENDDO
  434. ENDDO
  435. ENDDO
  436. END SUBROUTINE CALCULATE_STOCH_TEN
  437. ! ------------------------------------------------------------------
  438. SUBROUTINE UPDATE_STOCH_TEN(ru_tendf,rv_tendf,t_tendf, &
  439. GPUFORC,GPVFORC,GPTFORC, &
  440. mu,mub, &
  441. ids, ide, jds, jde, kds, kde, &
  442. ims, ime, jms, jme, kms, kme, &
  443. its, ite, jts, jte, kts, kte, &
  444. dt )
  445. IMPLICIT NONE
  446. INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
  447. ims, ime, jms, jme, kms, kme, &
  448. its, ite, jts, jte, kts, kte
  449. REAL , DIMENSION(ims:ime , kms:kme, jms:jme),INTENT(INOUT) :: &
  450. ru_tendf, rv_tendf, t_tendf
  451. !REAL , DIMENSION(ims:ime , kms:kme, jms:jme),INTENT(IN) :: &
  452. REAL , DIMENSION(ims:ime , kms:kme, jms:jme) :: &
  453. GPUFORC,GPVFORC,GPTFORC
  454. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: mu,mub
  455. INTEGER :: I,J,K
  456. REAL :: dt,xm
  457. DO j = jts,MIN(jde-1,jte)
  458. DO k = kts,kte-1
  459. DO i = its,ite
  460. ru_tendf(i,k,j) = ru_tendf(i,k,j) + GPUFORC(i,k,j) * (mu(i,j)+mub(i,j))
  461. ENDDO
  462. ENDDO
  463. ENDDO
  464. DO j = jts,jte
  465. DO i = its,MIN(ide-1,ite)
  466. DO k = kts,kte-1
  467. rv_tendf(i,k,j) = rv_tendf(i,k,j) + GPVFORC(i,k,j) * (mu(i,j)+mub(i,j))
  468. ENDDO
  469. ENDDO
  470. ENDDO
  471. DO j = jts,MIN(jde-1,jte)
  472. DO k = kts,kte-1
  473. DO i = its,MIN(ide-1,ite)
  474. t_tendf(i,k,j) = t_tendf(i,k,j) + GPTFORC(i,k,j) * (mu(i,j)+mub(i,j))
  475. ENDDO
  476. ENDDO
  477. ENDDO
  478. END SUBROUTINE UPDATE_STOCH_TEN
  479. ! ------------------------------------------------------------------
  480. !!************** TRANSFORM FROM SPHERICAL HARMONICS TO GRIDPOILT SPACE**
  481. ! ------------------------------------------------------------------
  482. subroutine SP2GP_prep( &
  483. SPSTREAMFORCS,SPSTREAMFORCC,SPTFORCS,SPTFORCC, &
  484. VERTSTRUCC,VERTSTRUCS, &
  485. RU_REAL,RV_REAL,RT_REAL, &
  486. RU_IMAG,RV_IMAG,RT_IMAG, &
  487. dx,dy,stoch_vertstruc_opt, &
  488. ids, ide, jds, jde, kds, kde, &
  489. ims, ime, jms, jme, kms, kme, &
  490. its, ite, jts, jte, kts, kte )
  491. IMPLICIT NONE
  492. INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
  493. ims, ime, jms, jme, kms, kme, &
  494. its, ite, jts, jte, kts, kte
  495. REAL, DIMENSION (ims:ime , jms:jme) :: SPSTREAMFORCS,SPSTREAMFORCC,SPTFORCS,SPTFORCC
  496. REAL, DIMENSION (ims:ime , kms:kme, jms:jme) :: RU_REAL,RV_REAL,RT_REAL,RU_IMAG,RV_IMAG,RT_IMAG, &
  497. VERTSTRUCC,VERTSTRUCS
  498. INTEGER :: IK,IL,ILEV,NLAT,NLON,stoch_vertstruc_opt
  499. REAL :: dx,dy,RY,RX
  500. NLAT=(jde-jds)+1 !KMAX
  501. NLON=(ide-ids)+1 !LMAX
  502. RY= NLAT*DY
  503. RX= NLON*DX
  504. DO ILEV=kts,kte
  505. if (stoch_vertstruc_opt==0) then
  506. DO IL=its,ite
  507. DO IK=jts,jte
  508. rt_real(IL,ILEV,IK) = SPTFORCC(IL,IK)
  509. rt_imag(IL,ILEV,IK) = SPTFORCS(IL,IK)
  510. ru_real(IL,ILEV,IK) = 2*RPI/RY* wavenumber_k(IK) * SPSTREAMFORCS(IL,IK)
  511. ru_imag(IL,ILEV,IK) =-2*RPI/RY* wavenumber_k(IK) * SPSTREAMFORCC(IL,IK)
  512. rv_real(IL,ILEV,IK) =-2*RPI/RX* wavenumber_l(IL) * SPSTREAMFORCS(IL,IK)
  513. rv_imag(IL,ILEV,IK) = 2*RPI/RX* wavenumber_l(IL) * SPSTREAMFORCC(IL,IK)
  514. ENDDO
  515. ENDDO
  516. elseif (stoch_vertstruc_opt==1) then
  517. DO IL=its,ite
  518. DO IK=jts,jte
  519. rt_real(IL,ILEV,IK) = SPTFORCC(IL,IK)*VERTSTRUCC(IL,ILEV,IK) - SPTFORCS(IL,IK)*VERTSTRUCS(IL,ILEV,IK)
  520. rt_imag(IL,ILEV,IK) = SPTFORCC(IL,IK)*VERTSTRUCS(IL,ILEV,IK) + SPTFORCS(IL,IK)*VERTSTRUCC(IL,ILEV,IK)
  521. ru_real(IL,ILEV,IK) = 2*RPI/RY* wavenumber_k(IK) *&
  522. (+SPSTREAMFORCC(IL,IK)*VERTSTRUCS(IL,ILEV,IK) + SPSTREAMFORCS(IL,IK)*VERTSTRUCC(IL,ILEV,IK))
  523. ru_imag(IL,ILEV,IK) = 2*RPI/RY* wavenumber_k(IK) *&
  524. (-SPSTREAMFORCC(IL,IK)*VERTSTRUCC(IL,ILEV,IK) + SPSTREAMFORCS(IL,IK)*VERTSTRUCS(IL,ILEV,IK))
  525. rv_real(IL,ILEV,IK) = 2*RPI/RX* wavenumber_l(IL) *&
  526. (-SPSTREAMFORCC(IL,IK)*VERTSTRUCS(IL,ILEV,IK) - SPSTREAMFORCS(IL,IK)*VERTSTRUCC(IL,ILEV,IK))
  527. rv_imag(IL,ILEV,IK) = 2*RPI/RX* wavenumber_l(IL) *&
  528. (+SPSTREAMFORCC(IL,IK)*VERTSTRUCC(IL,ILEV,IK) - SPSTREAMFORCS(IL,IK)*VERTSTRUCS(IL,ILEV,IK))
  529. ENDDO
  530. ENDDO
  531. endif
  532. ENDDO !ILEV
  533. END subroutine SP2GP_prep
  534. ! ------------------------------------------------------------------
  535. !!************** SUBROUTINE DO_FFTBACK_ALONG_X
  536. ! ------------------------------------------------------------------
  537. subroutine do_fftback_along_x(fieldc_U_xxx,fields_U_xxx, &
  538. fieldc_V_xxx,fields_V_xxx, &
  539. fieldc_T_xxx,fields_T_xxx, &
  540. ids, ide, jds, jde, kds, kde, &
  541. ims, ime, jms, jme, kms, kme, &
  542. ips, ipe, jps, jpe, kps, kpe, &
  543. imsx,imex,jmsx,jmex,kmsx,kmex, &
  544. ipsx,ipex,jpsx,jpex,kpsx,kpex, &
  545. imsy,imey,jmsy,jmey,kmsy,kmey, &
  546. ipsy,ipey,jpsy,jpey,kpsy,kpey, &
  547. k_start , k_end &
  548. )
  549. IMPLICIT NONE
  550. INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
  551. ims, ime, jms, jme, kms, kme, &
  552. ips, ipe, jps, jpe, kps, kpe, &
  553. imsx,imex,jmsx,jmex,kmsx,kmex, &
  554. ipsx,ipex,jpsx,jpex,kpsx,kpex, &
  555. imsy,imey,jmsy,jmey,kmsy,kmey, &
  556. ipsy,ipey,jpsy,jpey,kpsy,kpey, &
  557. k_start , k_end
  558. REAL, DIMENSION (imsx:imex, kmsx:kmex, jmsx:jmex) :: fieldc_U_xxx,fields_U_xxx, &
  559. fieldc_V_xxx,fields_V_xxx, &
  560. fieldc_T_xxx,fields_T_xxx
  561. COMPLEX, DIMENSION (ipsx:ipex) :: dummy_complex
  562. INTEGER :: IER,LENWRK,KMAX,LMAX,I,J,K
  563. REAL, ALLOCATABLE :: WORK(:)
  564. CHARACTER (LEN=160) :: mess
  565. KMAX=(jde-jds)+1
  566. LMAX=(ide-ids)+1
  567. LENWRK=2*KMAX*LMAX
  568. ALLOCATE(WORK(LENWRK))
  569. LENSAV= 4*(KMAX+LMAX)+INT(LOG(REAL(KMAX))) + INT(LOG(REAL(LMAX))) + 8
  570. DO k=kpsx,kpex
  571. DO j = jpsx, jpex
  572. DO i = ipsx, ipex
  573. dummy_complex(i)=cmplx(fieldc_U_xxx(i,k,j),fields_U_xxx(i,k,j))
  574. ENDDO
  575. CALL cFFT1B (LMAX, 1 ,dummy_complex,LMAX, WSAVE1, LENSAV, WORK, LENWRK, IER)
  576. if (ier.ne.0) then
  577. WRITE(mess,FMT='(A)') 'error in cFFT1B in do_fftback_along_x, field U'
  578. CALL wrf_debug(0,mess)
  579. end if
  580. DO i = ipsx, ipex
  581. fieldc_U_xxx(i,k,j)=real(dummy_complex(i))
  582. fields_U_xxx(i,k,j)=imag(dummy_complex(i))
  583. END DO
  584. END DO
  585. END DO
  586. DO k=kpsx,kpex
  587. DO j = jpsx, jpex
  588. DO i = ipsx, ipex
  589. dummy_complex(i)=cmplx(fieldc_V_xxx(i,k,j),fields_V_xxx(i,k,j))
  590. ENDDO
  591. CALL cFFT1B (LMAX, 1 ,dummy_complex,LMAX, WSAVE1, LENSAV, WORK, LENWRK, IER)
  592. if (ier.ne.0) then
  593. WRITE(mess,FMT='(A)') 'error in cFFT1B in do_fftback_along_x, field V'
  594. CALL wrf_debug(0,mess)
  595. end if
  596. DO i = ipsx,ipex
  597. fieldc_V_xxx(i,k,j)=real(dummy_complex(i))
  598. fields_V_xxx(i,k,j)=imag(dummy_complex(i))
  599. END DO
  600. END DO
  601. END DO
  602. DO k=kpsx,kpex
  603. DO j = jpsx, jpex
  604. DO i = ipsx, ipex
  605. dummy_complex(i)=cmplx(fieldc_T_xxx(i,k,j),fields_T_xxx(i,k,j))
  606. ENDDO
  607. CALL cFFT1B (LMAX, 1 ,dummy_complex,LMAX, WSAVE1, LENSAV, WORK, LENWRK, IER)
  608. if (ier.ne.0) then
  609. WRITE(mess,FMT='(A)') 'error in cFFT1B in do_fftback_along_x, field T'
  610. CALL wrf_debug(0,mess)
  611. end if
  612. DO i = ipsx, ipex
  613. fieldc_T_xxx(i,k,j)=real(dummy_complex(i))
  614. fields_T_xxx(i,k,j)=imag(dummy_complex(i))
  615. END DO
  616. END DO
  617. END DO !
  618. DEALLOCATE(WORK)
  619. end subroutine do_fftback_along_x
  620. !! ------------------------------------------------------------------
  621. !!!************** SUBROUTINE DO_FFTBACK_ALONG_Y
  622. !! ------------------------------------------------------------------
  623. subroutine do_fftback_along_y(fieldc_U_yyy,fields_U_yyy, &
  624. fieldc_V_yyy,fields_V_yyy, &
  625. fieldc_T_yyy,fields_T_yyy, &
  626. ids, ide, jds, jde, kds, kde, &
  627. ims, ime, jms, jme, kms, kme, &
  628. ips, ipe, jps, jpe, kps, kpe, &
  629. imsx,imex,jmsx,jmex,kmsx,kmex, &
  630. ipsx,ipex,jpsx,jpex,kpsx,kpex, &
  631. imsy,imey,jmsy,jmey,kmsy,kmey, &
  632. ipsy,ipey,jpsy,jpey,kpsy,kpey, &
  633. k_start , k_end &
  634. )
  635. IMPLICIT NONE
  636. INTEGER :: IER,LENWRK,KMAX,LMAX,I,J,K
  637. INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
  638. ims, ime, jms, jme, kms, kme, &
  639. ips, ipe, jps, jpe, kps, kpe, &
  640. imsx,imex,jmsx,jmex,kmsx,kmex, &
  641. ipsx,ipex,jpsx,jpex,kpsx,kpex, &
  642. imsy,imey,jmsy,jmey,kmsy,kmey, &
  643. ipsy,ipey,jpsy,jpey,kpsy,kpey, &
  644. k_start , k_end
  645. REAL, DIMENSION (imsy:imey, kmsy:kmey, jmsy:jmey) :: fieldc_U_yyy,fields_U_yyy, &
  646. fieldc_V_yyy,fields_V_yyy, &
  647. fieldc_T_yyy,fields_T_yyy
  648. COMPLEX, DIMENSION (jpsy:jpey) :: dummy_complex
  649. REAL, ALLOCATABLE :: WORK(:)
  650. CHARACTER (LEN=160) :: mess
  651. KMAX=(jde-jds)+1
  652. LMAX=(ide-ids)+1
  653. LENWRK=2*KMAX*LMAX
  654. ALLOCATE(WORK(LENWRK))
  655. LENSAV= 4*(KMAX+LMAX)+INT(LOG(REAL(KMAX))) + INT(LOG(REAL(LMAX))) + 8
  656. DO k=kpsy,kpey
  657. DO i = ipsy, ipey
  658. DO j = jpsy,jpey
  659. dummy_complex(j)=cmplx(fieldc_U_yyy(i,k,j),fields_U_yyy(i,k,j))
  660. ENDDO
  661. CALL cFFT1B (KMAX, 1 ,dummy_complex,KMAX, WSAVE2, LENSAV, WORK, LENWRK, IER)
  662. if (ier.ne.0) then
  663. WRITE(mess,FMT='(A)') 'error in cFFT1B in do_fftback_along_y, field U'
  664. CALL wrf_debug(0,mess)
  665. end if
  666. DO j = jpsy, jpey
  667. fieldc_U_yyy(i,k,j)=real(dummy_complex(j))
  668. fields_U_yyy(i,k,j)=imag(dummy_complex(j))
  669. END DO
  670. END DO
  671. END DO ! k_start-k_end
  672. DO k=kpsy,kpey
  673. DO i = ipsy, ipey
  674. DO j = jpsy, jpey
  675. dummy_complex(j)=cmplx(fieldc_V_yyy(i,k,j),fields_V_yyy(i,k,j))
  676. ENDDO
  677. CALL cFFT1B (KMAX, 1 ,dummy_complex,KMAX, WSAVE2, LENSAV, WORK, LENWRK, IER)
  678. if (ier.ne.0) then
  679. WRITE(mess,FMT='(A)') 'error in cFFT1B in do_fftback_along_y, field V'
  680. CALL wrf_debug(0,mess)
  681. end if
  682. DO j = jpsy, jpey
  683. fieldc_V_yyy(i,k,j)=real(dummy_complex(j))
  684. fields_V_yyy(i,k,j)=imag(dummy_complex(j))
  685. END DO
  686. END DO
  687. END DO ! k_start-k_end
  688. DO k=kpsy,kpey
  689. DO i = ipsy, ipey
  690. DO j = jpsy,jpey
  691. dummy_complex(j)=cmplx(fieldc_T_yyy(i,k,j),fields_T_yyy(i,k,j))
  692. ENDDO
  693. CALL cFFT1B (KMAX, 1 ,dummy_complex,KMAX, WSAVE2, LENSAV, WORK, LENWRK, IER)
  694. if (ier.ne.0) then
  695. WRITE(mess,FMT='(A)') 'error in cFFT1B in do_fftback_along_y, field T'
  696. CALL wrf_debug(0,mess)
  697. end if
  698. DO j = jpsy,jpey
  699. fieldc_T_yyy(i,k,j)=real(dummy_complex(j))
  700. fields_T_yyy(i,k,j)=imag(dummy_complex(j))
  701. END DO
  702. END DO
  703. END DO ! k_start-k_end
  704. DEALLOCATE(WORK)
  705. end subroutine do_fftback_along_y
  706. ! ------------------------------------------------------------------
  707. !!************** TRANSFORM FROM GRIDPOILT SPACE TO SPHERICAL HARMONICS **
  708. ! ------------------------------------------------------------------
  709. subroutine findindex( wavenumber_k, wavenumber_L, &
  710. ids, ide, jds, jde, kds, kde, &
  711. ims, ime, jms, jme, kms, kme, &
  712. its, ite, jts, jte, kts, kte )
  713. IMPLICIT NONE
  714. INTEGER :: IK,IL,KMAX,LMAX
  715. INTEGER, DIMENSION (jds:jde):: wavenumber_k
  716. INTEGER, DIMENSION (ids:ide):: wavenumber_l
  717. INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
  718. ims, ime, jms, jme, kms, kme, &
  719. its, ite, jts, jte, kts, kte
  720. KMAX=(jde-jds)+1
  721. LMAX=(ide-ids)+1
  722. !map wave numbers K,L to indeces IK, IL
  723. DO IK=1,KMAX/2+1
  724. wavenumber_k(IK)=IK-1
  725. ENDDO
  726. DO IK=KMAX,KMAX/2+2,-1
  727. wavenumber_k(IK)=IK-KMAX-1
  728. ENDDO
  729. DO IL=1,LMAX/2+1
  730. wavenumber_l(IL)=IL-1
  731. ENDDO
  732. DO IL=LMAX,LMAX/2+2,-1
  733. wavenumber_l(IL)=IL-LMAX-1
  734. ENDDO
  735. END subroutine findindex
  736. ! ------------------------------------------------------------------
  737. subroutine gauss_noise(z)
  738. real :: z ! output
  739. real :: x,y,r, coeff ! INPUT
  740. ! [2.1] Get two uniform variate random numbers IL range 0 to 1:
  741. do
  742. call random_number( x )
  743. call random_number( y )
  744. ! [2.2] Transform to range -1 to 1 and calculate sum of squares:
  745. x = 2.0 * x - 1.0
  746. y = 2.0 * y - 1.0
  747. r = x * x + y * y
  748. if ( r > 0.0 .and. r < 1.0 ) exit
  749. end do
  750. !
  751. ! [2.3] Use Box-Muller transformation to get normal deviates:
  752. coeff = sqrt( -2.0 * log(r) / r )
  753. z = coeff * x
  754. end subroutine gauss_noise
  755. ! ------------------------------------------------------------------
  756. SUBROUTINE rand_seed (config_flags, seed1, seed2,nens )
  757. USE module_configure
  758. IMPLICIT NONE
  759. !
  760. ! Structure that contains run-time configuration (namelist) data for domain
  761. TYPE (grid_config_rec_type) :: config_flags
  762. !
  763. ! Arguments
  764. INTEGER, INTENT(OUT) :: seed1, seed2
  765. INTEGER, INTENT(IN ) :: nens
  766. ! Local
  767. integer :: date_time(8)
  768. integer*8 :: yyyy,mmdd,newtime
  769. integer*8 :: ihr,isc,idiv
  770. integer*8 :: iphys
  771. character (len=10) :: real_clock(3), time
  772. !
  773. LOGICAL :: is_print = .false.
  774. !
  775. ! Read CPU time
  776. call date_and_time(real_clock(1), real_clock(2),&
  777. real_clock(3), date_time)
  778. read(real_clock(1),'(i4)') yyyy
  779. read(real_clock(1),'(4x,i4)') mmdd
  780. read(real_clock(2),'(i6,1x,i3)') ihr,isc ! ihr = hhmmss , isc - milliseconds of the minute
  781. newtime = yyyy*mmdd+ihr+isc
  782. if (is_print) then
  783. print*,'real_clock: ',real_clock
  784. print*,'real_clock(1): ',real_clock(1)
  785. print*,'real_clock(2): ',real_clock(2)
  786. print*,'real_clock(3): ',real_clock(3)
  787. print*,'date_time: ',date_time
  788. print*,'newtime: ',newtime
  789. !
  790. ! get a seed number (w/ physics options for each ensemble member)
  791. !
  792. print *,'physics_options:',&
  793. config_flags%sf_surface_physics,&
  794. config_flags%ra_lw_physics, &
  795. config_flags%ra_sw_physics, &
  796. config_flags%mp_physics, &
  797. config_flags%sf_sfclay_physics,&
  798. config_flags%bl_pbl_physics,&
  799. config_flags%cu_physics
  800. endif !isprint
  801. iphys = config_flags%sf_surface_physics * 1000000 + &
  802. config_flags%ra_lw_physics * 100000 + &
  803. config_flags%ra_sw_physics * 10000 + &
  804. config_flags%mp_physics * 1000 + &
  805. config_flags%sf_sfclay_physics * 100 + &
  806. config_flags%bl_pbl_physics * 10 + &
  807. config_flags%cu_physics
  808. idiv=2;
  809. seed1 = newtime+iphys+nens*1000000
  810. seed2 = mod(newtime+iphys+nens*1000000,idiv)
  811. if(is_print) print *,'Rand_seed (iphys/newtime/idiv):',iphys,newtime,idiv,nens
  812. end SUBROUTINE rand_seed
  813. end module module_stoch