source: trunk/src/initialize.f90 @ 28

Last change on this file since 28 was 4, checked in by mlanger, 11 years ago
File size: 8.7 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(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
72  implicit none
73
74  integer :: itime
75  integer :: ldt,nrand
76  integer(kind=2) :: icbt
77  real :: zt,dz,dz1,dz2,up,vp,wp,usigold,vsigold,wsigold,ran3
78  real(kind=dp) :: xt,yt
79  save idummy
80
81  integer :: idummy = -7
82
83  icbt=1           ! initialize particle to "no reflection"
84
85  nrand=int(ran3(idummy)*real(maxrand-1))+1
86
87
88  !******************************
89  ! 2. Interpolate necessary data
90  !******************************
91
92  ! Compute maximum mixing height around particle position
93  !*******************************************************
94
95  ix=int(xt)
96  jy=int(yt)
97  ixp=ix+1
98  jyp=jy+1
99
100  h=max(hmix(ix ,jy ,1,memind(1)), &
101       hmix(ixp,jy ,1,memind(1)), &
102       hmix(ix ,jyp,1,memind(1)), &
103       hmix(ixp,jyp,1,memind(1)), &
104       hmix(ix ,jy ,1,memind(2)), &
105       hmix(ixp,jy ,1,memind(2)), &
106       hmix(ix ,jyp,1,memind(2)), &
107       hmix(ixp,jyp,1,memind(2)))
108
109  zeta=zt/h
110
111
112  !*************************************************************
113  ! If particle is in the PBL, interpolate once and then make a
114  ! time loop until end of interval is reached
115  !*************************************************************
116
117  if (zeta.le.1.) then
118
119    call interpol_all(itime,real(xt),real(yt),zt)
120
121
122  ! Vertical interpolation of u,v,w,rho and drhodz
123  !***********************************************
124
125  ! Vertical distance to the level below and above current position
126  ! both in terms of (u,v) and (w) fields
127  !****************************************************************
128
129    dz1=zt-height(indz)
130    dz2=height(indzp)-zt
131    dz=1./(dz1+dz2)
132
133    u=(dz1*uprof(indzp)+dz2*uprof(indz))*dz
134    v=(dz1*vprof(indzp)+dz2*vprof(indz))*dz
135    w=(dz1*wprof(indzp)+dz2*wprof(indz))*dz
136
137
138  ! Compute the turbulent disturbances
139
140  ! Determine the sigmas and the timescales
141  !****************************************
142
143    if (turbswitch) then
144      call hanna(zt)
145    else
146      call hanna1(zt)
147    endif
148
149
150  ! Determine the new diffusivity velocities
151  !*****************************************
152
153    if (nrand+2.gt.maxrand) nrand=1
154    up=rannumb(nrand)*sigu
155    vp=rannumb(nrand+1)*sigv
156    wp=rannumb(nrand+2)
157    if (.not.turbswitch) wp=wp*sigw
158
159
160  ! Determine time step for next integration
161  !*****************************************
162
163    if (turbswitch) then
164      ldt=int(min(tlw,h/max(2.*abs(wp*sigw),1.e-5), &
165           0.5/abs(dsigwdz),600.)*ctl)
166    else
167      ldt=int(min(tlw,h/max(2.*abs(wp),1.e-5),600.)*ctl)
168    endif
169    ldt=max(ldt,mintime)
170
171
172    usig=(usigprof(indzp)+usigprof(indz))/2.
173    vsig=(vsigprof(indzp)+vsigprof(indz))/2.
174    wsig=(wsigprof(indzp)+wsigprof(indz))/2.
175
176  else
177
178
179
180  !**********************************************************
181  ! For all particles that are outside the PBL, make a single
182  ! time step. Only horizontal turbulent disturbances are
183  ! calculated. Vertical disturbances are reset.
184  !**********************************************************
185
186
187  ! Interpolate the wind
188  !*********************
189
190    call interpol_wind(itime,real(xt),real(yt),zt)
191
192
193  ! Compute everything for above the PBL
194
195  ! Assume constant turbulent perturbations
196  !****************************************
197
198    ldt=abs(lsynctime)
199
200    if (nrand+1.gt.maxrand) nrand=1
201    up=rannumb(nrand)*0.3
202    vp=rannumb(nrand+1)*0.3
203    nrand=nrand+2
204    wp=0.
205    sigw=0.
206
207  endif
208
209  !****************************************************************
210  ! Add mesoscale random disturbances
211  ! This is done only once for the whole lsynctime interval to save
212  ! computation time
213  !****************************************************************
214
215
216  ! It is assumed that the average interpolation error is 1/2 sigma
217  ! of the surrounding points, autocorrelation time constant is
218  ! 1/2 of time interval between wind fields
219  !****************************************************************
220
221  if (nrand+2.gt.maxrand) nrand=1
222  usigold=rannumb(nrand)*usig
223  vsigold=rannumb(nrand+1)*vsig
224  wsigold=rannumb(nrand+2)*wsig
225
226end subroutine initialize
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG