source: branches/flexpart91_hasod/src_parallel/initialize.f90 @ 8

Last change on this file since 8 was 8, checked in by hasod, 11 years ago

Added parallel version of Flexpart91

File size: 8.8 KB
Line 
1!**********************************************************************
2! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
3! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
4! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
5!                                                                     *
6! This file is part of FLEXPART.                                      *
7!                                                                     *
8! FLEXPART is free software: you can redistribute it and/or modify    *
9! it under the terms of the GNU General Public License as published by*
10! the Free Software Foundation, either version 3 of the License, or   *
11! (at your option) any later version.                                 *
12!                                                                     *
13! FLEXPART is distributed in the hope that it will be useful,         *
14! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
15! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
16! GNU General Public License for more details.                        *
17!                                                                     *
18! You should have received a copy of the GNU General Public License   *
19! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
20!**********************************************************************
21
22subroutine initialize(pj, itime,ldt,up,vp,wp, &
23       usigold,vsigold,wsigold,xt,yt,zt,icbt)
24  !                        i    i   o  o  o
25  !        o       o       o    i  i  i   o
26  !*****************************************************************************
27  !                                                                            *
28  !  Calculation of trajectories utilizing a zero-acceleration scheme. The time*
29  !  step is determined by the Courant-Friedrichs-Lewy (CFL) criterion. This   *
30  !  means that the time step must be so small that the displacement within    *
31  !  this time step is smaller than 1 grid distance. Additionally, a temporal  *
32  !  CFL criterion is introduced: the time step must be smaller than the time  *
33  !  interval of the wind fields used for interpolation.                       *
34  !  For random walk simulations, these are the only time step criteria.       *
35  !  For the other options, the time step is also limited by the Lagrangian    *
36  !  time scale.                                                               *
37  !                                                                            *
38  !     Author: A. Stohl                                                       *
39  !                                                                            *
40  !     16 December 1997                                                       *
41  !                                                                            *
42  !  Literature:                                                               *
43  !                                                                            *
44  !*****************************************************************************
45  !                                                                            *
46  ! Variables:                                                                 *
47  ! h [m]              Mixing height                                           *
48  ! lwindinterv [s]    time interval between two wind fields                   *
49  ! itime [s]          current temporal position                               *
50  ! ldt [s]            Suggested time step for next integration                *
51  ! ladvance [s]       Total integration time period                           *
52  ! rannumb(maxrand)   normally distributed random variables                   *
53  ! up,vp,wp           random velocities due to turbulence                     *
54  ! usig,vsig,wsig     uncertainties of wind velocities due to interpolation   *
55  ! usigold,vsigold,wsigold  like usig, etc., but for the last time step       *
56  ! xt,yt,zt           Next time step's spatial position of trajectory         *
57  !                                                                            *
58  !                                                                            *
59  ! Constants:                                                                 *
60  ! cfl                factor, by which the time step has to be smaller than   *
61  !                    the spatial CFL-criterion                               *
62  ! cflt               factor, by which the time step has to be smaller than   *
63  !                    the temporal CFL-criterion                              *
64  !                                                                            *
65  !*****************************************************************************
66
67  use par_mod
68  use com_mod
69  use interpol_mod
70  use hanna_mod
71  use random_mod, only : ran3_initialize
72  ! reads: ran3_initialize
73  ! writes:
74
75  implicit none
76
77  integer, intent(in)   :: pj
78  integer, intent(in)   :: itime
79  integer :: ldt,nrand
80  integer(kind=2) :: icbt
81  real :: zt,dz,dz1,dz2,up,vp,wp,usigold,vsigold,wsigold,ran3
82  real(kind=dp) :: xt,yt
83  save idummy
84
85  integer :: idummy = -7
86
87  icbt=1           ! initialize particle to "no reflection"
88
89  nrand=int(ran3_initialize(pj)*real(maxrand-1))+1
90
91
92  !******************************
93  ! 2. Interpolate necessary data
94  !******************************
95
96  ! Compute maximum mixing height around particle position
97  !*******************************************************
98
99  ix=int(xt)
100  jy=int(yt)
101  ixp=ix+1
102  jyp=jy+1
103
104  h=max(hmix(ix ,jy ,1,memind(1)), &
105       hmix(ixp,jy ,1,memind(1)), &
106       hmix(ix ,jyp,1,memind(1)), &
107       hmix(ixp,jyp,1,memind(1)), &
108       hmix(ix ,jy ,1,memind(2)), &
109       hmix(ixp,jy ,1,memind(2)), &
110       hmix(ix ,jyp,1,memind(2)), &
111       hmix(ixp,jyp,1,memind(2)))
112
113  zeta=zt/h
114
115
116  !*************************************************************
117  ! If particle is in the PBL, interpolate once and then make a
118  ! time loop until end of interval is reached
119  !*************************************************************
120
121  if (zeta.le.1.) then
122
123    call interpol_all(itime,real(xt),real(yt),zt)
124
125
126  ! Vertical interpolation of u,v,w,rho and drhodz
127  !***********************************************
128
129  ! Vertical distance to the level below and above current position
130  ! both in terms of (u,v) and (w) fields
131  !****************************************************************
132
133    dz1=zt-height(indz)
134    dz2=height(indzp)-zt
135    dz=1./(dz1+dz2)
136
137    u=(dz1*uprof(indzp)+dz2*uprof(indz))*dz
138    v=(dz1*vprof(indzp)+dz2*vprof(indz))*dz
139    w=(dz1*wprof(indzp)+dz2*wprof(indz))*dz
140
141
142  ! Compute the turbulent disturbances
143
144  ! Determine the sigmas and the timescales
145  !****************************************
146
147    if (turbswitch) then
148      call hanna(zt)
149    else
150      call hanna1(zt)
151    endif
152
153
154  ! Determine the new diffusivity velocities
155  !*****************************************
156
157    if (nrand+2.gt.maxrand) nrand=1
158    up=rannumb(nrand)*sigu
159    vp=rannumb(nrand+1)*sigv
160    wp=rannumb(nrand+2)
161    if (.not.turbswitch) wp=wp*sigw
162
163
164  ! Determine time step for next integration
165  !*****************************************
166
167    if (turbswitch) then
168      ldt=int(min(tlw,h/max(2.*abs(wp*sigw),1.e-5), &
169           0.5/abs(dsigwdz),600.)*ctl)
170    else
171      ldt=int(min(tlw,h/max(2.*abs(wp),1.e-5),600.)*ctl)
172    endif
173    ldt=max(ldt,mintime)
174
175
176    usig=(usigprof(indzp)+usigprof(indz))/2.
177    vsig=(vsigprof(indzp)+vsigprof(indz))/2.
178    wsig=(wsigprof(indzp)+wsigprof(indz))/2.
179
180  else
181
182
183
184  !**********************************************************
185  ! For all particles that are outside the PBL, make a single
186  ! time step. Only horizontal turbulent disturbances are
187  ! calculated. Vertical disturbances are reset.
188  !**********************************************************
189
190
191  ! Interpolate the wind
192  !*********************
193
194    call interpol_wind(itime,real(xt),real(yt),zt)
195
196
197  ! Compute everything for above the PBL
198
199  ! Assume constant turbulent perturbations
200  !****************************************
201
202    ldt=abs(lsynctime)
203
204    if (nrand+1.gt.maxrand) nrand=1
205    up=rannumb(nrand)*0.3
206    vp=rannumb(nrand+1)*0.3
207    nrand=nrand+2
208    wp=0.
209    sigw=0.
210
211  endif
212
213  !****************************************************************
214  ! Add mesoscale random disturbances
215  ! This is done only once for the whole lsynctime interval to save
216  ! computation time
217  !****************************************************************
218
219
220  ! It is assumed that the average interpolation error is 1/2 sigma
221  ! of the surrounding points, autocorrelation time constant is
222  ! 1/2 of time interval between wind fields
223  !****************************************************************
224
225  if (nrand+2.gt.maxrand) nrand=1
226  usigold=rannumb(nrand)*usig
227  vsigold=rannumb(nrand+1)*vsig
228  wsigold=rannumb(nrand+2)*wsig
229
230end subroutine initialize
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG