FLEXPART.f90 12.2 KB
Newer Older
Matthias Langer's avatar
 
Matthias Langer committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
!                                                                     *
! This file is part of FLEXPART.                                      *
!                                                                     *
! FLEXPART is free software: you can redistribute it and/or modify    *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or   *
! (at your option) any later version.                                 *
!                                                                     *
! FLEXPART is distributed in the hope that it will be useful,         *
! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
! GNU General Public License for more details.                        *
!                                                                     *
! You should have received a copy of the GNU General Public License   *
! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
!**********************************************************************

program flexpart

  !*****************************************************************************
  !                                                                            *
  !     This is the Lagrangian Particle Dispersion Model FLEXPART.             *
  !     The main program manages the reading of model run specifications, etc. *
  !     All actual computing is done within subroutine timemanager.            *
  !                                                                            *
  !     Author: A. Stohl                                                       *
  !                                                                            *
  !     18 May 1996                                                            *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  !                                                                            *
  ! Constants:                                                                 *
  !                                                                            *
  !*****************************************************************************

  use point_mod
  use par_mod
  use com_mod
  use conv_mod

  implicit none

  integer :: i,j,ix,jy,inest
  integer :: idummy = -320
51
52
  character(len=256) :: inline_options  !pathfile, flexversion, arg2

Matthias Langer's avatar
 
Matthias Langer committed
53
54
55
56
57
58
59
60
61

  ! Generate a large number of random numbers
  !******************************************

  do i=1,maxrand-1,2
    call gasdev1(idummy,rannumb(i),rannumb(i+1))
  end do
  call gasdev1(idummy,rannumb(maxrand),rannumb(maxrand-1))

62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
  ! 
  flexversion='Version 9.1.8  (2013-12-08)'
  !verbosity=0
  ! Read the pathnames where input/output files are stored
  !*******************************************************

  inline_options='none'
  select case (iargc())
  case (2)
    call getarg(1,arg1)
    pathfile=arg1
    call getarg(2,arg2)
    inline_options=arg2
  case (1)
    call getarg(1,arg1)
    pathfile=arg1
    verbosity=0
    if (arg1(1:1).eq.'-') then
        write(pathfile,'(a11)') './pathnames'
        inline_options=arg1 
    endif
  case (0)
    write(pathfile,'(a11)') './pathnames'
    verbosity=0
  end select
  
    if (inline_options(1:1).eq.'-') then
      print*, 'inline options=', inline_options       
      if (trim(inline_options).eq.'-v'.or.trim(inline_options).eq.'-v1') then
         print*, 'verbose mode 1: additional information will be displayed'
         verbosity=1
      endif
      if (trim(inline_options).eq.'-v2') then
         print*, 'verbose mode 2: additional information will be displayed'
         verbosity=2
      endif
      if (trim(inline_options).eq.'-i') then
         print*, 'info mode: will provide run specific information and stop'
         verbosity=1
         info_flag=1
      endif
      if (trim(inline_options).eq.'-i2') then
         print*, 'info mode: will provide run specific information and stop'
         verbosity=2
         info_flag=1
      endif
    endif


Matthias Langer's avatar
 
Matthias Langer committed
111
112
  ! Print the GPL License statement
  !*******************************************************
113
  print*,'Welcome to FLEXPART', trim(flexversion)
Matthias Langer's avatar
 
Matthias Langer committed
114
115
  print*,'FLEXPART is free software released under the GNU Genera'// &
       'l Public License.'
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
            
  if (verbosity.gt.0) then
        WRITE(*,*) 'call readpaths'
  endif 
  call readpaths(pathfile)
   
 
  if (verbosity.gt.1) then !show clock info 
     !print*,'length(4)',length(4)
     !count=0,count_rate=1000
     CALL SYSTEM_CLOCK(count_clock0, count_rate, count_max)
     !WRITE(*,*) 'SYSTEM_CLOCK',count, count_rate, count_max
     !WRITE(*,*) 'SYSTEM_CLOCK, count_clock0', count_clock0
     !WRITE(*,*) 'SYSTEM_CLOCK, count_rate', count_rate
     !WRITE(*,*) 'SYSTEM_CLOCK, count_max', count_max
  endif
Matthias Langer's avatar
 
Matthias Langer committed
132
133
134
135
136


  ! Read the user specifications for the current model run
  !*******************************************************

137
138
139
  if (verbosity.gt.0) then
        WRITE(*,*) 'call readcommand'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
140
  call readcommand
141
142
143
144
145
146
147
148
149
  if (verbosity.gt.0) then
        WRITE(*,*) '    ldirect=', ldirect
        WRITE(*,*) '    ibdate,ibtime=',ibdate,ibtime
        WRITE(*,*) '    iedate,ietime=', iedate,ietime
        if (verbosity.gt.1) then   
                CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
                WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
        endif     
  endif
Matthias Langer's avatar
 
Matthias Langer committed
150
151
152
153


  ! Read the age classes to be used
  !********************************
154
155
156
  if (verbosity.gt.0) then
        WRITE(*,*) 'call readageclasses'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
157
158
  call readageclasses

159
160
161
162
163
164
  if (verbosity.gt.1) then   
                CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
                WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
  endif     
 

Matthias Langer's avatar
 
Matthias Langer committed
165
166
167
168

  ! Read, which wind fields are available within the modelling period
  !******************************************************************

169
170
171
  if (verbosity.gt.0) then
        WRITE(*,*) 'call readavailable'
  endif  
Matthias Langer's avatar
 
Matthias Langer committed
172
173
174
175
176
  call readavailable

  ! Read the model grid specifications,
  ! both for the mother domain and eventual nests
  !**********************************************
177
178
179
180
 
  if (verbosity.gt.0) then
     WRITE(*,*) 'call gridcheck'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
181
182

  call gridcheck
183
184
185
186
187
188
189
190
191
192

  if (verbosity.gt.1) then   
     CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
     WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
  endif      
  

  if (verbosity.gt.0) then
        WRITE(*,*) 'call gridcheck_nests'
  endif  
Matthias Langer's avatar
 
Matthias Langer committed
193
194
195
196
197
198
  call gridcheck_nests


  ! Read the output grid specifications
  !************************************

199
200
201
202
  if (verbosity.gt.0) then
        WRITE(*,*) 'call readoutgrid'
  endif

Matthias Langer's avatar
 
Matthias Langer committed
203
204
  call readoutgrid

205
206
207
208
209
210
  if (nested_output.eq.1) then
          call readoutgrid_nest
    if (verbosity.gt.0) then
        WRITE(*,*) '# readoutgrid_nest'
    endif
  endif
Matthias Langer's avatar
 
Matthias Langer committed
211
212
213
214

  ! Read the receptor points for which extra concentrations are to be calculated
  !*****************************************************************************

215
216
217
  if (verbosity.eq.1) then
     print*,'call readreceptors'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
218
219
220
221
222
223
224
225
226
227
228
  call readreceptors

  ! Read the physico-chemical species property table
  !*************************************************
  !SEC: now only needed SPECIES are read in readreleases.f
  !call readspecies


  ! Read the landuse inventory
  !***************************

229
230
231
  if (verbosity.gt.0) then
    print*,'call readlanduse'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
232
233
234
235
236
237
  call readlanduse


  ! Assign fractional cover of landuse classes to each ECMWF grid point
  !********************************************************************

238
239
240
  if (verbosity.gt.0) then
    print*,'call assignland'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
241
242
243
244
245
246
247
  call assignland



  ! Read the coordinates of the release locations
  !**********************************************

248
249
250
  if (verbosity.gt.0) then
    print*,'call readreleases'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
251
252
  call readreleases

253

Matthias Langer's avatar
 
Matthias Langer committed
254
255
256
  ! Read and compute surface resistances to dry deposition of gases
  !****************************************************************

257
258
259
  if (verbosity.gt.0) then
    print*,'call readdepo'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
260
261
262
263
264
  call readdepo

  ! Convert the release point coordinates from geografical to grid coordinates
  !***************************************************************************

265
266
267
268
  call coordtrafo  
  if (verbosity.gt.0) then
    print*,'call coordtrafo'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
269
270
271
272
273


  ! Initialize all particles to non-existent
  !*****************************************

274
275
276
  if (verbosity.gt.0) then
    print*,'Initialize all particles to non-existent'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
277
278
279
280
281
282
283
284
  do j=1,maxpart
    itra1(j)=-999999999
  end do

  ! For continuation of previous run, read in particle positions
  !*************************************************************

  if (ipin.eq.1) then
285
286
287
    if (verbosity.gt.0) then
      print*,'call readpartpositions'
    endif
Matthias Langer's avatar
 
Matthias Langer committed
288
289
    call readpartpositions
  else
290
291
292
    if (verbosity.gt.0) then
      print*,'numpart=0, numparticlecount=0'
    endif    
Matthias Langer's avatar
 
Matthias Langer committed
293
294
295
296
297
298
299
300
301
    numpart=0
    numparticlecount=0
  endif


  ! Calculate volume, surface area, etc., of all output grid cells
  ! Allocate fluxes and OHfield if necessary
  !***************************************************************

302
303
304
305

  if (verbosity.gt.0) then
    print*,'call outgrid_init'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
306
307
308
309
310
311
312
  call outgrid_init
  if (nested_output.eq.1) call outgrid_init_nest


  ! Read the OH field
  !******************

313
314
315
316
  if (OHREA.eqv..TRUE.) then
    if (verbosity.gt.0) then
       print*,'call readOHfield'
    endif
Matthias Langer's avatar
 
Matthias Langer committed
317
       call readOHfield
318
  endif
Matthias Langer's avatar
 
Matthias Langer committed
319
320
321
322
323

  ! Write basic information on the simulation to a file "header"
  ! and open files that are to be kept open throughout the simulation
  !******************************************************************

324
325
326
327
328

  if (verbosity.gt.0) then
    print*,'call writeheader'
  endif

Matthias Langer's avatar
 
Matthias Langer committed
329
  call writeheader
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
  ! FLEXPART 9.2 ticket ?? write header in ASCII format 
  call writeheader_txt
  !if (nested_output.eq.1) call writeheader_nest
  if (nested_output.eq.1.and.surf_only.ne.1) call writeheader_nest

  if (nested_output.eq.1.and.surf_only.eq.1) call writeheader_nest_surf
  if (nested_output.ne.1.and.surf_only.eq.1) call writeheader_surf



  !open(unitdates,file=path(2)(1:length(2))//'dates')

  if (verbosity.gt.0) then
    print*,'call openreceptors'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
345
346
347
348
349
350
351
  call openreceptors
  if ((iout.eq.4).or.(iout.eq.5)) call openouttraj


  ! Releases can only start and end at discrete times (multiples of lsynctime)
  !***************************************************************************

352
353
354
  if (verbosity.gt.0) then
    print*,'discretize release times'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
355
356
357
358
359
360
361
362
363
364
365
  do i=1,numpoint
    ireleasestart(i)=nint(real(ireleasestart(i))/ &
         real(lsynctime))*lsynctime
    ireleaseend(i)=nint(real(ireleaseend(i))/ &
         real(lsynctime))*lsynctime
  end do


  ! Initialize cloud-base mass fluxes for the convection scheme
  !************************************************************

366
367
368
369
  if (verbosity.gt.0) then
    print*,'Initialize cloud-base mass fluxes for the convection scheme'
  endif

Matthias Langer's avatar
 
Matthias Langer committed
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
  do jy=0,nymin1
    do ix=0,nxmin1
      cbaseflux(ix,jy)=0.
    end do
  end do
  do inest=1,numbnests
    do jy=0,nyn(inest)-1
      do ix=0,nxn(inest)-1
        cbasefluxn(ix,jy,inest)=0.
      end do
    end do
  end do


  ! Calculate particle trajectories
  !********************************

387
388
389
390
391
392
393
394
395
396
397
398
  if (verbosity.gt.0) then
     if (verbosity.gt.1) then   
       CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
       WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
     endif
     if (info_flag.eq.1) then
         print*, 'info only mode (stop)'    
         stop
     endif
     print*,'call timemanager'
  endif

Matthias Langer's avatar
 
Matthias Langer committed
399
400
401
402
403
404
405
  call timemanager


  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
       &XPART MODEL RUN!'

end program flexpart