FLEXPART.f90 14.1 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
64
65
66
67
68
  ! flexversion='Version 9.2 beta (2014-07-01)'
  flexversion='Version 9.2.0.1 (2015-01-27)'
  ! default inlide options
  inline_options='none'
  !verbosity flags  defined in com_mod.f90
 
69
70
71
72
  ! Read the pathnames where input/output files are stored
  !*******************************************************

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

  ! inline options allow to fine tune the verbosity during run time
  ! e.g.: show compilation parameters or input variables, time execution     
96
  if (inline_options(1:1).eq.'-') then
97
98
99
100
101
102
103
104
105
106
107
108
109
   ! 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'
110
       verbosity=1
111
       if (inline_options(index_v+1:index_v+1).eq.'2') then
112
       verbosity=2
113
114
115
       endif
       print*, 'verbosity level=', verbosity !inline_options(index_v+1:index_v+1)
              
116
    endif
117
118
119
120
121
122
123
124
125
    !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   
    !if (trim(inline_options).eq.'-i') then
       index_v=index(inline_options,'i')
       print*, 'Info mode: provide run specific information and stop'
126
127
       verbosity=1
       info_flag=1
128
129
130
131
132
133
       !if (trim(inline_options).eq.'-i2') then
       if (inline_options(index_v+1:index_v+1).eq.'2') then
           print*, 'Including input files'
       !   verbosity=1
       info_flag=2
       endif
134
    endif
135
136
137
138
139
140
141
142
143
    !if (trim(inline_options).eq.'-i2') then
    !   print*, 'Info mode: provide more detailed run specific information and stop'
    !   verbosity=1
    !   info_flag=2
    !endif
    if (index(inline_options,'t').gt.0) then
       time_flag=1
       print*, 'timing execution: not implemented'
       !stop
144
    endif
145
146
147
148
149
150
    if (index(inline_options,'d').gt.0) then
       debug_flag=1
       print*, 'debug: not implemented'
       print*, 'debug_flag=', debug_flag
       !stop
    endif 
151
152
  endif
           
153
  if (verbosity.gt.0) then
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
    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
171
    write(*,*) 'call readpaths'
172
173
174
  endif 
  call readpaths(pathfile)
 
175
  !if (time_flag.gt.1) then !show clock info 
176
     !count=0,count_rate=1000
177
  CALL SYSTEM_CLOCK(count_clock0, count_rate, count_max)
178
179
180
181
     !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
182
  !endif
Matthias Langer's avatar
 
Matthias Langer committed
183
184
185
186

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

187
  if (verbosity.gt.0) then
188
    write(*,*) 'call readcommand'
189
  endif
Matthias Langer's avatar
 
Matthias Langer committed
190
  call readcommand
191
  if (verbosity.gt.0) then
192
193
    write(*,*) '    ldirect      =', ldirect
    write(*,*) '    ibdate,ibtime=', ibdate,ibtime
194
    write(*,*) '    iedate,ietime=', iedate,ietime
195
196
  endif
    if (time_flag.gt.0) then   
197
198
199
      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
200
201
202

  ! Read the age classes to be used
  !********************************
203
  if (verbosity.gt.0) then
204
    write(*,*) 'call readageclasses'
205
  endif
Matthias Langer's avatar
 
Matthias Langer committed
206
207
  call readageclasses

208
  if (time_flag.gt.1) then   
209
210
    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
211
  endif     
Matthias Langer's avatar
 
Matthias Langer committed
212
213
214
215

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

216
  if (verbosity.gt.0) then
217
    write(*,*) 'call readavailable'
218
  endif  
Matthias Langer's avatar
 
Matthias Langer committed
219
220
221
222
223
  call readavailable

  ! Read the model grid specifications,
  ! both for the mother domain and eventual nests
  !**********************************************
224
225
 
  if (verbosity.gt.0) then
226
     write(*,*) 'FLEXPART> call gridcheck'
227
  endif
Matthias Langer's avatar
 
Matthias Langer committed
228
  call gridcheck
229

230
  if (time_flag.gt.0) then   
231
232
    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
233
234
235
  endif      

  if (verbosity.gt.0) then
236
    write(*,*) 'call gridcheck_nests'
237
  endif  
Matthias Langer's avatar
 
Matthias Langer committed
238
239
240
241
242
  call gridcheck_nests

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

243
  if (verbosity.gt.0) then
244
    write(*,*) 'FLEXPART> call readoutgrid'
245
246
  endif

Matthias Langer's avatar
 
Matthias Langer committed
247
248
  call readoutgrid

249
  if (nested_output.eq.1) then
250
    call readoutgrid_nest
251
    if (verbosity.gt.0) then
252
      write(*,*) 'FLEXPART>  readoutgrid_nest'
253
254
    endif
  endif
Matthias Langer's avatar
 
Matthias Langer committed
255
256
257
258

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

259
  if (verbosity.eq.1) then
260
     print*,'FLEXPART> call readreceptors'
261
  endif
Matthias Langer's avatar
 
Matthias Langer committed
262
263
264
265
266
267
268
269
270
271
272
  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
  !***************************

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

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

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

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

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

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

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

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

305
306
307
308
  call coordtrafo  
  if (verbosity.gt.0) then
    print*,'call coordtrafo'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
309
310
311
312

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

313
314
315
  if (verbosity.gt.0) then
    print*,'Initialize all particles to non-existent'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
316
317
318
319
320
321
322
323
  do j=1,maxpart
    itra1(j)=-999999999
  end do

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

  if (ipin.eq.1) then
324
325
326
    if (verbosity.gt.0) then
      print*,'call readpartpositions'
    endif
Matthias Langer's avatar
 
Matthias Langer committed
327
328
    call readpartpositions
  else
329
    if (verbosity.gt.0) then
330
      print*,'set numpart=0, numparticlecount=0'
331
    endif    
Matthias Langer's avatar
 
Matthias Langer committed
332
333
334
335
336
337
338
339
    numpart=0
    numparticlecount=0
  endif

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

340
341
342
  if (verbosity.gt.0) then
    print*,'call outgrid_init'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
343
344
345
346
347
348
  call outgrid_init
  if (nested_output.eq.1) call outgrid_init_nest

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

349
350
  if (OHREA.eqv..TRUE.) then
    if (verbosity.gt.0) then
351
      print*,'call readOHfield'
352
    endif
353
    call readOHfield
354
  endif
Matthias Langer's avatar
 
Matthias Langer committed
355
356
357
358
359

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

360
361
362
363
  if (verbosity.gt.0) then
    print*,'call writeheader'
  endif

Matthias Langer's avatar
 
Matthias Langer committed
364
  call writeheader
365
  ! write header in ASCII format 
366
367
368
369
370
371
372
373
374
375
376
  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
377
378
379
380
381
382
  call openreceptors
  if ((iout.eq.4).or.(iout.eq.5)) call openouttraj

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

383
384
385
  if (verbosity.gt.0) then
    print*,'discretize release times'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
386
  do i=1,numpoint
387
388
    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
389
390
391
392
393
  end do

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

394
395
396
397
  if (verbosity.gt.0) then
    print*,'Initialize cloud-base mass fluxes for the convection scheme'
  endif

Matthias Langer's avatar
 
Matthias Langer committed
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
  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
  !********************************

414
415
416
417
418
419
420
421
  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
    print*, 'info only mode (stop before call timemanager)'
    stop
  endif
422
423
424
425
  if (verbosity.gt.0) then
     print*,'call timemanager'
  endif

Matthias Langer's avatar
 
Matthias Langer committed
426
427
  call timemanager

428
  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLEXPART MODEL RUN!'
Matthias Langer's avatar
 
Matthias Langer committed
429
430

end program flexpart