Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_stoch.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 36.2 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format