Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

initialize.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 11 years ago
sources for flexwrf v3.1
File size: 13.1 KB

Line
1	!***********************************************************************
2	!* Copyright 2012,2013 *
3	!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine, *
4	!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
5	!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso, *
6	!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann, *
7	!* *
8	!* This file is part of FLEXPART WRF *
9	!* *
10	!* FLEXPART is free software: you can redistribute it and/or modify *
11	!* it under the terms of the GNU General Public License as published by*
12	!* the Free Software Foundation, either version 3 of the License, or *
13	!* (at your option) any later version. *
14	!* *
15	!* FLEXPART is distributed in the hope that it will be useful, *
16	!* but WITHOUT ANY WARRANTY; without even the implied warranty of *
17	!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
18	!* GNU General Public License for more details. *
19	!* *
20	!* You should have received a copy of the GNU General Public License *
21	!* along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
22	!***********************************************************************
23	subroutine initialize(itime,ldt,up,vp,wp, &
24	usigold,vsigold,wsigold,xt,yt,zt,icbt, &
25	ngrid,depoindicator,indzindicator,cpt2,ompid,myid,n_threads,mts )
26	! uprof,vprof,wprof, usigprof,vsigprof,wsigprof, &
27	! rhoprof,rhogradprof, tkeprof,pttprof, &
28	! u,v,w,usig,vsig,wsig,pvi, &
29	!! p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2, &
30	! ix,jy,ixp,jyp,ngrid,indz,indzp, depoindicator, &
31	! indzindicator, &
32	! ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
33	! sigw,dsigwdz,dsigw2dz,cpt,ompid)
34
35	! i i o o o
36	! o o o i i i o
37	!*******************************************************************************
38	! *
39	! Calculation of trajectories utilizing a zero-acceleration scheme. *
40	! The time step is determined by the Courant-Friedrichs-Lewy (CFL) criterion. *
41	! This means that the time step must be so small that the displacement within *
42	! this time step is smaller than 1 grid distance. Additionally, a temporal *
43	! CFL criterion is introduced: the time step must be smaller than the time *
44	! interval of the wind fields used for interpolation. *
45	! For random walk simulations, these are the only time step criteria. *
46	! For the other options, the time step is also limited by the Lagrangian time *
47	! scale. *
48	! *
49	! Author: A. Stohl *
50	! *
51	! 16 December 1997 *
52	! *
53	! Literature: *
54	! *
55	!*******************************************************************************
56	! *
57	! Variables: *
58	! h [m] Mixing height *
59	! lwindinterv [s] time interval between two wind fields *
60	! itime [s] current temporal position *
61	! ldt [s] Suggested time step for next integration *
62	! ladvance [s] Total integration time period *
63	! rannumb(maxrand) normally distributed random variables *
64	! up,vp,wp random velocities due to turbulence *
65	! usig,vsig,wsig uncertainties of wind velocities due to interpolation *
66	! usigold,vsigold,wsigold like usig, etc., but for the last time step *
67	! xt,yt,zt Next time step's spatial position of trajectory *
68	! *
69	! *
70	! Constants: *
71	! cfl factor, by which the time step has to be smaller than the *
72	! spatial CFL-criterion *
73	! cflt factor, by which the time step has to be smaller than the *
74	! temporal CFL-criterion *
75	! *
76	!*******************************************************************************
77	! 12 JUNE 2007 W. Wang
78	! use WRF TKE option to compute turbulence
79	! Mar 2012: J. Brioude modification to handle openmp. *
80	! Jan 2013 M. Cassiani modification to use CBL scheme
81	!*******************************************************************************
82	use par_mod
83	use com_mod
84	use mt_stream
85
86	! use interpol_mod
87	! use hanna_mod
88	! use ran_mod
89	implicit none
90
91	integer :: itime
92	integer :: ldt,nrand,ompid
93	!OMP_GET_THREAD_NUM
94	integer(kind=2) :: icbt
95	real :: zt,dz,dz1,dz2,up,vp,wp,usigold,vsigold,wsigold,ran3
96	real(kind=dp) :: xt,yt
97	! save idummy
98
99	integer :: idummy = -7
100
101	real :: uprof(nzmax),vprof(nzmax),wprof(nzmax)
102	real :: usigprof(nzmax),vsigprof(nzmax),wsigprof(nzmax)
103	real :: rhoprof(nzmax),rhogradprof(nzmax)
104	real :: tkeprof(nzmax),pttprof(nzmax)
105	real :: u,v,w,usig,vsig,wsig,pvi,mu,mv
106
107	real :: p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2
108	integer :: ix,jy,ixp,jyp,ngrid,indz,indzp,cpt2,maxrand2
109	logical :: depoindicator(maxspec)
110	logical :: indzindicator(nzmax)
111
112	real :: ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw
113	real :: sigw,dsigwdz,dsigw2dz
114
115	real :: dcas,dcas1,dcas2 !modified by by mc, random number needed in initialize_cbl_vel
116	integer :: myid,n_threads !added by mc for parallel random number generation
117	integer(4) :: rannum
118	real(4) :: real_rannum
119	type (mt_state) :: mts (0: MAX_STREAM)
120
121
122	idummy=7
123	icbt=1 ! initialize particle to "no reflection"
124
125	if (newrandomgen.eq.0) then
126	cpt2=cpt2+1
127	! cpt=cpt+1000367
128	cpt2=mod(cpt2,maxrand)+1;
129
130	! nrand=int(ran3(idummy,inext,inextp,ma,iff)*real(maxrand-1))+1
131	! nrand=int(ran3(idummy)*real(maxrand-1))+1
132	nrand=cpt2+ompid*maxrand
133	maxrand2=maxrandomp
134	else
135	!mc
136	! rannum=genrand_int32(mts(ompid+1+(myid*n_threads))) !integer random
137	! number at 32 bit resolution
138	rannum=genrand_int32(mts(ompid+1)) !integer random number at 32 bit resolution
139	real_rannum = sngl(0.5_DP + 0.2328306e-9_DP * rannum) !conversion to single precision 32bit real between 0-1
140	nrand=int(real_rannum*real(maxrand-1))+1
141	maxrand2=maxrand
142	endif
143
144	!******************************
145	! 2. Interpolate necessary data
146	!******************************
147
148	! Compute maximum mixing height around particle position
149	!*******************************************************
150
151	ix=int(xt)
152	jy=int(yt)
153	ixp=ix+1
154	jyp=jy+1
155
156	h=max(hmix(ix ,jy ,1,memind(1)), &
157	hmix(ixp,jy ,1,memind(1)), &
158	hmix(ix ,jyp,1,memind(1)), &
159	hmix(ixp,jyp,1,memind(1)), &
160	hmix(ix ,jy ,1,memind(2)), &
161	hmix(ixp,jy ,1,memind(2)), &
162	hmix(ix ,jyp,1,memind(2)), &
163	hmix(ixp,jyp,1,memind(2)))
164
165	! JB
166	zeta=zt/h
167
168	!*************************************************************
169	! If particle is in the PBL, interpolate once and then make a
170	! time loop until end of interval is reached
171	!*************************************************************
172
173	if (zeta.le.1.) then
174
175	call interpol_all(itime,real(xt),real(yt),zt, &
176	uprof,vprof,wprof, usigprof,vsigprof,wsigprof, &
177	rhoprof,rhogradprof, tkeprof,pttprof, &
178	u,v,w,usig,vsig,wsig,pvi, &
179	p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2, &
180	ix,jy,ixp,jyp,ngrid,indz,indzp, depoindicator, &
181	indzindicator, &
182	ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
183	sigw,dsigwdz,dsigw2dz,mu,mv)
184
185
186	! Vertical interpolation of u,v,w,rho and drhodz
187	!***********************************************
188
189	! Vertical distance to the level below and above current position
190	! both in terms of (u,v) and (w) fields
191	!****************************************************************
192
193	dz1=zt-height(indz)
194	dz2=height(indzp)-zt
195	dz=1./(dz1+dz2)
196
197	u=(dz1uprof(indzp)+dz2uprof(indz))*dz
198	v=(dz1vprof(indzp)+dz2vprof(indz))*dz
199	w=(dz1wprof(indzp)+dz2wprof(indz))*dz
200
201
202	! Compute the turbulent disturbances
203
204	! Determine the sigmas and the timescales
205	!****************************************
206	! FLEXPART WRF
207	! write(,)'initial.f','turb_option=',turb_option
208	! write(,)'turb_option_mytke=',turb_option_mytke
209	if (turb_option .eq. turb_option_mytke) then
210	! write(,)'initial.f'
211	call tke_partition_my(zt, &
212	ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
213	sigw,dsigwdz,dsigw2dz,uprof,vprof,tkeprof,pttprof,indz,indzp)
214	elseif (turb_option .eq. turb_option_tke) then
215	call tke_partition_hanna(zt, &
216	ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
217	sigw,dsigwdz,dsigw2dz,uprof,vprof,tkeprof,pttprof,indz,indzp)
218	else
219
220	if (turbswitch) then
221	call hanna(zt, &
222	ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
223	sigw,dsigwdz,dsigw2dz)
224
225	else
226	call hanna1(zt, &
227	ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
228	sigw,dsigwdz,dsigw2dz)
229	endif
230	endif
231
232	! Determine the new diffusivity velocities
233	!*****************************************
234
235	if (nrand+2.gt.maxrand2) nrand=1
236	up=rannumb(nrand)*sigu
237	vp=rannumb(nrand+1)*sigv
238	wp=rannumb(nrand+2)
239	nrand=nrand+2 ! added by mc: it was missing previously a bug I think here and in original flexpart
240
241	if (.not.turbswitch) then ! modified by mc
242	wp=wp*sigw
243	else if (cblflag.eq.1) then ! modified by mc
244	if (-h/ol.gt.5) then !unstable conditions from -h/ol >5
245	!if (ol.lt.0.) then
246	!if (ol.gt.0.) then !by mc : gt.0 is only for test the correct is lt.0^M
247	dcas=uniform_rannumb(nrand) !uniform^M
248	dcas1=rannumb(nrand) !gaussian^M
249	nrand=nrand+3
250	call initialize_cbl_vel(idummy,zt,ust,wst,h,sigw,wp,dcas,dcas1,ol)
251	else
252	wp=wp*sigw
253	end if
254	end if
255	! Determine time step for next integration
256	!*****************************************
257
258	if (turbswitch) then
259	ldt=int(min(tlw,h/max(2.abs(wpsigw),1.e-5), &
260	0.5/abs(dsigwdz),600.)*ctl)
261	else
262	ldt=int(min(tlw,h/max(2.abs(wp),1.e-5),600.)ctl)
263	endif
264	ldt=max(ldt,mintime)
265
266
267	usig=(usigprof(indzp)+usigprof(indz))/2.
268	vsig=(vsigprof(indzp)+vsigprof(indz))/2.
269	wsig=(wsigprof(indzp)+wsigprof(indz))/2.
270
271	else
272
273
274
275	!**********************************************************
276	! For all particles that are outside the PBL, make a single
277	! time step. Only horizontal turbulent disturbances are
278	! calculated. Vertical disturbances are reset.
279	!**********************************************************
280
281
282	! Interpolate the wind
283	!*********************
284
285	call interpol_wind(itime,real(xt),real(yt),zt, &
286	uprof,vprof,wprof, usigprof,vsigprof,wsigprof, &
287	rhoprof,rhogradprof, tkeprof,pttprof, &
288	u,v,w,usig,vsig,wsig,pvi, &
289	p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2, &
290	ix,jy,ixp,jyp,ngrid,indz,indzp, depoindicator, &
291	indzindicator,mu,mv)
292
293
294	! Compute everything for above the PBL
295
296	! Assume constant turbulent perturbations
297	!****************************************
298
299	ldt=abs(lsynctime)
300
301	if (nrand+1.gt.maxrand2) nrand=1
302	up=rannumb(nrand)*0.3
303	vp=rannumb(nrand+1)*0.3
304	nrand=nrand+2
305	wp=0.
306	sigw=0.
307
308	endif
309
310	!****************************************************************
311	! Add mesoscale random disturbances
312	! This is done only once for the whole lsynctime interval to save
313	! computation time
314	!****************************************************************
315
316
317	! It is assumed that the average interpolation error is 1/2 sigma
318	! of the surrounding points, autocorrelation time constant is
319	! 1/2 of time interval between wind fields
320	!****************************************************************
321
322	if (nrand+2.gt.maxrand2) nrand=1
323	usigold=rannumb(nrand)*usig
324	vsigold=rannumb(nrand+1)*vsig
325	wsigold=rannumb(nrand+2)*wsig
326
327	end subroutine initialize
328

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

Download in other formats:

Original Format