FLEXPART.f90 14 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
  character(len=256) :: inline_options  !pathfile, flexversion, arg2
52
  integer :: index_v
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
  ! FLEXPART version string
63
  ! flexversion='Version 9.2 beta (2014-07-01)'
64
65
  !flexversion='Version 9.2.0.1 (2015-01-27)'
  flexversion='Version 9.2.0.2 (2015-03-01)'
66
67
68
69
  ! default inlide options
  inline_options='none'
  !verbosity flags  defined in com_mod.f90
 
70
71
72
73
  ! Read the pathnames where input/output files are stored
  !*******************************************************

  select case (iargc())
74
  case (2) !2 parameters: pathfile and inline options
75
76
77
78
    call getarg(1,arg1)
    pathfile=arg1
    call getarg(2,arg2)
    inline_options=arg2
79
  case (1) !1 parameter pathfiel or inline options
80
81
82
    call getarg(1,arg1)
    pathfile=arg1
    if (arg1(1:1).eq.'-') then
83
84
      write(pathfile,'(a11)') './pathnames'
      inline_options=arg1 
85
    endif
86
  case (0) !default behavior
87
88
89
    write(pathfile,'(a11)') './pathnames'
  end select
  
Matthias Langer's avatar
 
Matthias Langer committed
90
91
  ! Print the GPL License statement
  !*******************************************************
92
  print*,'Welcome to FLEXPART ', trim(flexversion)
93
  print*,'FLEXPART is free software released under the GNU General Public License.'
94
95
96

  ! inline options allow to fine tune the verbosity during run time
  ! e.g.: show compilation parameters or input variables, time execution     
97
  if (inline_options(1:1).eq.'-') then
98
99
100
101
102
103
104
105
106
107
108
109
110
   ! if (index(inline_options,'v').gt.0) then
   !    print*, 'verbose mode'
   !    verbosity=1
   !    index_v=index(inline_options,'v')
   !    if (inline_options(index_v+1:index_v+1).eq.'2') then
   !    verbosity=2
   !    endif         
   ! endif   
 
    !if (trim(inline_options).eq.'-v'.or.trim(inline_options).eq.'-v1') then
    if (index(inline_options,'v').gt.0) then
       index_v=index(inline_options,'v')
       print*, 'Verbose mode: display  additional information during run'
111
       verbosity=1
112
       if (inline_options(index_v+1:index_v+1).eq.'2') then
113
       verbosity=2
114
115
       endif
       print*, 'verbosity level=', verbosity !inline_options(index_v+1:index_v+1)
116
    endif
117
118
119
120
121
122
123
    !iif (trim(inline_options).eq.'-v2') then
    !   print*, 'Verbose mode 2: display more detailed information during run'
    !   verbosity=2
    !endif

    if (index(inline_options,'i').gt.0) then   
       index_v=index(inline_options,'i')
124
       print*, 'Info mode: provide compile and run specific information, then stop'
125
126
       verbosity=1
       info_flag=1
127
128
129
       if (inline_options(index_v+1:index_v+1).eq.'2') then
       info_flag=2
       endif
130
    endif
131
132
    if (index(inline_options,'t').gt.0) then
       time_flag=1
133
       print*, 'timing execution activated'
134
       !stop
135
    endif
136
137
    if (index(inline_options,'d').gt.0) then
       debug_flag=1
138
       print*, 'debug messages activated'
139
       print*, 'debug_flag=', debug_flag
140
       !these messages give additional info on top on verbose mode
141
    endif 
142
143
  endif
           
144
  if (verbosity.gt.0) then
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
    print*, 'FLEXPART>******************************'
    print*, 'FLEXPART>* verbosity level:', verbosity
    print*, 'FLEXPART>* info only:      ', info_flag
    print*, 'FLEXPART>* time execution: ', time_flag
    print*, 'FLEXPART>******************************'
    
    print*, 'FLEXPART> parameters from par_mod'    
    print*, 'FLEXPART> nxmax=  ', nxmax
    print*, 'FLEXPART> nymax=  ', nymax
    print*, 'FLEXPART> nuvzmax=', nuvzmax
    print*, 'FLEXPART> nwzmax= ', nwzmax
    print*, 'FLEXPART> nzmax=  ', nzmax
    print*, 'FLEXPART> nxshift=', nxshift
    print*, 'FLEXPART> maxpart=', maxpart
    print*, 'FLEXPART> maxspec=', maxspec 

    if (info_flag.eq.1) stop
162
    write(*,*) 'call readpaths'
163
164
165
  endif 
  call readpaths(pathfile)
 
166
  !if (time_flag.gt.1) then !show clock info 
167
     !count=0,count_rate=1000
168
  CALL SYSTEM_CLOCK(count_clock0, count_rate, count_max)
169
170
171
172
     !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
173
  !endif
Matthias Langer's avatar
 
Matthias Langer committed
174
175
176
177

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

178
  if (verbosity.gt.0) then
179
    write(*,*) 'FLEXPART> call readcommand'
180
  endif
Matthias Langer's avatar
 
Matthias Langer committed
181
  call readcommand
182
  if (verbosity.gt.0) then
183
184
    write(*,*) '    ldirect      =', ldirect
    write(*,*) '    ibdate,ibtime=', ibdate,ibtime
185
    write(*,*) '    iedate,ietime=', iedate,ietime
186
187
  endif
    if (time_flag.gt.0) then   
188
189
190
      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
191
192
193

  ! Read the age classes to be used
  !********************************
194
  if (verbosity.gt.0) then
195
    write(*,*) 'FLEXPART> call readageclasses'
196
  endif
Matthias Langer's avatar
 
Matthias Langer committed
197
198
  call readageclasses

199
  if (time_flag.gt.1) then   
200
201
    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
202
  endif     
Matthias Langer's avatar
 
Matthias Langer committed
203
204
205
206

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

207
  if (verbosity.gt.0) then
208
    write(*,*) 'FLEXPART> call readavailable'
209
  endif  
Matthias Langer's avatar
 
Matthias Langer committed
210
211
212
213
214
  call readavailable

  ! Read the model grid specifications,
  ! both for the mother domain and eventual nests
  !**********************************************
215
216
 
  if (verbosity.gt.0) then
217
     write(*,*) 'FLEXPART> call gridcheck'
218
  endif
Matthias Langer's avatar
 
Matthias Langer committed
219
  call gridcheck
220

221
  if (time_flag.gt.0) then   
222
223
    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
224
225
226
  endif      

  if (verbosity.gt.0) then
227
    write(*,*) 'FLEXPART> call gridcheck_nests'
228
  endif  
Matthias Langer's avatar
 
Matthias Langer committed
229
230
231
232
233
  call gridcheck_nests

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

234
  if (verbosity.gt.0) then
235
    write(*,*) 'FLEXPART> call readoutgrid'
236
237
  endif

Matthias Langer's avatar
 
Matthias Langer committed
238
239
  call readoutgrid

240
  if (nested_output.eq.1) then
241
    call readoutgrid_nest
242
    if (verbosity.gt.0) then
243
      write(*,*) 'FLEXPART>  readoutgrid_nest'
244
245
    endif
  endif
Matthias Langer's avatar
 
Matthias Langer committed
246
247
248
249

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

250
  if (verbosity.eq.1) then
251
     print*,'FLEXPART> call readreceptors'
252
  endif
Matthias Langer's avatar
 
Matthias Langer committed
253
254
255
256
257
258
259
260
261
262
263
  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
  !***************************

264
  if (verbosity.gt.0) then
265
    print*,'FLEXPART> call readlanduse'
266
  endif
Matthias Langer's avatar
 
Matthias Langer committed
267
268
269
270
271
  call readlanduse

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

272
  if (verbosity.gt.0) then
273
    print*,'FLEXPART> call assignland'
274
  endif
Matthias Langer's avatar
 
Matthias Langer committed
275
276
277
278
279
  call assignland

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

280
  if (verbosity.gt.0) then
281
    print*,'FLEXPART> call readreleases'
282
  endif
Matthias Langer's avatar
 
Matthias Langer committed
283
284
285
286
287
  call readreleases

  ! Read and compute surface resistances to dry deposition of gases
  !****************************************************************

288
  if (verbosity.gt.0) then
289
    print*,'FLEXPART> call readdepo'
290
  endif
Matthias Langer's avatar
 
Matthias Langer committed
291
292
293
294
295
  call readdepo

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

296
297
  call coordtrafo  
  if (verbosity.gt.0) then
298
    print*,'FLEXPART> call coordtrafo'
299
  endif
Matthias Langer's avatar
 
Matthias Langer committed
300
301
302
303

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

304
  if (verbosity.gt.0) then
305
    print*,'FLEXPART> Initialize all particles to non-existent'
306
  endif
Matthias Langer's avatar
 
Matthias Langer committed
307
308
309
310
311
312
313
314
  do j=1,maxpart
    itra1(j)=-999999999
  end do

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

  if (ipin.eq.1) then
315
    if (verbosity.gt.0) then
316
      print*,'FLEXPART> call readpartpositions'
317
    endif
Matthias Langer's avatar
 
Matthias Langer committed
318
319
    call readpartpositions
  else
320
    if (verbosity.gt.0) then
321
      print*,'FLEXPART> set numpart=0, numparticlecount=0'
322
    endif    
Matthias Langer's avatar
 
Matthias Langer committed
323
324
325
326
327
328
329
330
    numpart=0
    numparticlecount=0
  endif

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

331
  if (verbosity.gt.0) then
332
    print*,'FLEXPART> call outgrid_init'
333
  endif
Matthias Langer's avatar
 
Matthias Langer committed
334
335
336
337
338
339
  call outgrid_init
  if (nested_output.eq.1) call outgrid_init_nest

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

340
341
  if (OHREA.eqv..TRUE.) then
    if (verbosity.gt.0) then
342
      print*,'FLEXPART> call readOHfield'
343
    endif
344
    call readOHfield
345
  endif
Matthias Langer's avatar
 
Matthias Langer committed
346
347
348
349
350

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

351
  if (verbosity.gt.0) then
352
    print*,'FLEXPART> call variuos writeheader routines'
353
354
  endif

Matthias Langer's avatar
 
Matthias Langer committed
355
  call writeheader
356
  ! write header in ASCII format 
357
358
359
360
361
362
363
364
365
  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
366
    print*,'FLEXPART> call openreceptors'
367
  endif
Matthias Langer's avatar
 
Matthias Langer committed
368
369
370
371
372
373
  call openreceptors
  if ((iout.eq.4).or.(iout.eq.5)) call openouttraj

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

374
  if (verbosity.gt.0) then
375
    print*,'FLEXPART> discretize release times'
376
  endif
Matthias Langer's avatar
 
Matthias Langer committed
377
  do i=1,numpoint
378
379
    ireleasestart(i)=nint(real(ireleasestart(i))/real(lsynctime))*lsynctime
    ireleaseend(i)=nint(real(ireleaseend(i))/real(lsynctime))*lsynctime
Matthias Langer's avatar
 
Matthias Langer committed
380
381
382
383
384
  end do

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

385
  if (verbosity.gt.0) then
386
    print*,'FLEXPART> Initialize cloud-base mass fluxes for the convection scheme'
387
388
  endif

Matthias Langer's avatar
 
Matthias Langer committed
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
  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
  !********************************

405
406
407
408
409
  if (time_flag.gt.0) 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.2) then
410
    print*, 'FLEXPART> info only mode (stop before call timemanager)'
411
412
    stop
  endif
413
  if (verbosity.gt.0) then
414
     print*,'FLEXPART> call timemanager'
415
416
  endif

Matthias Langer's avatar
 
Matthias Langer committed
417
418
  call timemanager

419
  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLEXPART MODEL RUN!'
Matthias Langer's avatar
 
Matthias Langer committed
420
421

end program flexpart