FLEXPART.f90 13.9 KB
Newer Older
Matthias Langer's avatar
 
Matthias Langer committed
1
2
3
4
5
6
7
8
9
10
11
12
13
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                                                            *
  !                                                                            *
  !*****************************************************************************
14
15
16
17
18
19
20
21
  ! Changes:                                                                   *
  !   Unified ECMWF and GFS builds                                             *
  !   Marian Harustak, 12.5.2017                                               *
  !     - Added detection of metdata format using gributils routines           *
  !     - Distinguished calls to ecmwf/gfs gridcheck versions based on         *
  !       detected metdata format                                              *
  !     - Passed metdata format down to timemanager                            *
  !*****************************************************************************
Matthias Langer's avatar
 
Matthias Langer committed
22
23
24
25
26
27
28
29
30
31
32
  !                                                                            *
  ! Variables:                                                                 *
  !                                                                            *
  ! Constants:                                                                 *
  !                                                                            *
  !*****************************************************************************

  use point_mod
  use par_mod
  use com_mod
  use conv_mod
33

34
  use random_mod, only: gasdev1
35
  use class_gribfile
Matthias Langer's avatar
 
Matthias Langer committed
36

37
38
39
40
#ifdef USE_NCF
  use netcdf_output_mod, only: writeheader_netcdf
#endif

Matthias Langer's avatar
 
Matthias Langer committed
41
42
  implicit none

43
  integer :: i,j,ix,jy,inest, iopt
Matthias Langer's avatar
 
Matthias Langer committed
44
  integer :: idummy = -320
45
  character(len=256) :: inline_options  !pathfile, flexversion, arg2
46
47
48
  integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN
  integer :: detectformat

49
  
Matthias Langer's avatar
 
Matthias Langer committed
50
51
52
53
54
55
56
57
  ! 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))

58

59
  ! FLEXPART version string
60
  flexversion_major = '10' ! Major version number, also used for species file names
61
  flexversion='Version '//trim(flexversion_major)//'.4 (2019-11-12)'
62
63
  verbosity=0

64
65
66
  ! Read the pathnames where input/output files are stored
  !*******************************************************

67
  inline_options='none'
68
  select case (iargc())
69
  case (2)
70
71
72
73
    call getarg(1,arg1)
    pathfile=arg1
    call getarg(2,arg2)
    inline_options=arg2
74
  case (1)
75
76
77
    call getarg(1,arg1)
    pathfile=arg1
    if (arg1(1:1).eq.'-') then
78
79
      write(pathfile,'(a11)') './pathnames'
      inline_options=arg1 
80
    endif
81
  case (0)
82
83
84
    write(pathfile,'(a11)') './pathnames'
  end select
  
Matthias Langer's avatar
 
Matthias Langer committed
85
86
  ! Print the GPL License statement
  !*******************************************************
87
  print*,'Welcome to FLEXPART ', trim(flexversion)
88
  print*,'FLEXPART is free software released under the GNU General Public License.'
89
 
90
91
92

  ! Ingest inline options
  !*******************************************************
93
  if (inline_options(1:1).eq.'-') then
94
95
96
97
98
99
100
101
102
103
    print*,'inline_options:',inline_options
    !verbose mode
    iopt=index(inline_options,'v') 
    if (iopt.gt.0) then
      verbosity=1
      !print*, iopt, inline_options(iopt+1:iopt+1)
      if  (trim(inline_options(iopt+1:iopt+1)).eq.'2') then
        print*, 'Verbose mode 2: display more detailed information during run'
        verbosity=2
      endif
104
    endif
105
106
107
108
    !debug mode 
    iopt=index(inline_options,'d')
    if (iopt.gt.0) then
      debug_mode=.true.
109
    endif
110
111
    if (trim(inline_options).eq.'-i') then
       print*, 'Info mode: provide detailed run specific information and stop'
112
113
114
       verbosity=1
       info_flag=1
    endif
115
116
117
118
    if (trim(inline_options).eq.'-i2') then
       print*, 'Info mode: provide more detailed run specific information and stop'
       verbosity=2
       info_flag=1
119
120
121
    endif
  endif
           
122
  if (verbosity.gt.0) then
Ignacio Pisso's avatar
Ignacio Pisso committed
123
124
125
126
    print*, 'nxmax=',nxmax
    print*, 'nymax=',nymax
    print*, 'nzmax=',nzmax
    print*,'nxshift=',nxshift 
127
128
129
  endif
  
  if (verbosity.gt.0) then
130
    write(*,*) 'call readpaths'
131
  endif 
132
  call readpaths
133
 
134
135
  if (verbosity.gt.1) then !show clock info 
     !print*,'length(4)',length(4)
136
     !count=0,count_rate=1000
137
     CALL SYSTEM_CLOCK(count_clock0, count_rate, count_max)
138
139
140
141
     !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
142
  endif
Matthias Langer's avatar
 
Matthias Langer committed
143
144
145
146

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

147
  if (verbosity.gt.0) then
148
    write(*,*) 'call readcommand'
149
  endif
Matthias Langer's avatar
 
Matthias Langer committed
150
  call readcommand
151
  if (verbosity.gt.0) then
152
153
    write(*,*) '    ldirect=', ldirect
    write(*,*) '    ibdate,ibtime=',ibdate,ibtime
154
    write(*,*) '    iedate,ietime=', iedate,ietime
155
    if (verbosity.gt.1) then   
156
157
158
      CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
      write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
    endif     
159
  endif
160
161
162
163
164

  ! Initialize arrays in com_mod
  !*****************************
  call com_mod_allocate_part(maxpart)

Matthias Langer's avatar
 
Matthias Langer committed
165
166
167

  ! Read the age classes to be used
  !********************************
168
  if (verbosity.gt.0) then
169
    write(*,*) 'call readageclasses'
170
  endif
Matthias Langer's avatar
 
Matthias Langer committed
171
172
  call readageclasses

173
  if (verbosity.gt.1) then   
174
175
    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
176
  endif     
Matthias Langer's avatar
 
Matthias Langer committed
177
178
179
180

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

181
  if (verbosity.gt.0) then
182
    write(*,*) 'call readavailable'
183
  endif  
Matthias Langer's avatar
 
Matthias Langer committed
184
185
  call readavailable

186
187
  ! Detect metdata format
  !**********************
188
189
190
  if (verbosity.gt.0) then
    write(*,*) 'call detectformat'
  endif
191
192
193
194
195
196
197
198
199

  metdata_format = detectformat()

  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
    print *,'ECMWF metdata detected'
  elseif (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
    print *,'NCEP metdata detected'
  else
    print *,'Unknown metdata format'
200
    stop
201
202
203
204
  endif



205
206
207
208
209
210
211
212
  ! If nested wind fields are used, allocate arrays
  !************************************************

  if (verbosity.gt.0) then
    write(*,*) 'call com_mod_allocate_nests'
  endif
  call com_mod_allocate_nests

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

221
222
223
224
225
  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
    call gridcheck_ecmwf
  else
    call gridcheck_gfs
  end if
226

227
  if (verbosity.gt.1) then   
228
229
    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
230
231
232
  endif      

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

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

240
  if (verbosity.gt.0) then
241
    write(*,*) 'call readoutgrid'
242
243
  endif

Matthias Langer's avatar
 
Matthias Langer committed
244
245
  call readoutgrid

246
  if (nested_output.eq.1) then
247
    call readoutgrid_nest
248
    if (verbosity.gt.0) then
249
      write(*,*) '# readoutgrid_nest'
250
251
    endif
  endif
Matthias Langer's avatar
 
Matthias Langer committed
252
253
254
255

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

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

270
  if (verbosity.gt.0) then
271
    print*,'call readlanduse'
272
  endif
Matthias Langer's avatar
 
Matthias Langer committed
273
274
275
276
277
  call readlanduse

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

278
  if (verbosity.gt.0) then
279
    print*,'call assignland'
280
  endif
Matthias Langer's avatar
 
Matthias Langer committed
281
282
283
284
285
  call assignland

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

286
  if (verbosity.gt.0) then
287
    print*,'call readreleases'
288
  endif
Matthias Langer's avatar
 
Matthias Langer committed
289
290
291
292
293
  call readreleases

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

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

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

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

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

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

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

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

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

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

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

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

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

357
#ifdef USE_NCF
358
359
360
361
362
363
364
365
366
367
368
369
370
  if (lnetcdfout.eq.1) then 
    call writeheader_netcdf(lnest=.false.)
  else 
    call writeheader
  end if

  if (nested_output.eq.1) then
    if (lnetcdfout.eq.1) then
      call writeheader_netcdf(lnest=.true.)
    else
      call writeheader_nest
    endif
  endif
371
#endif
372

373
  if (verbosity.gt.0) then
374
    print*,'call writeheader'
375
376
  endif

Matthias Langer's avatar
 
Matthias Langer committed
377
  call writeheader
378
  ! FLEXPART 9.2 ticket ?? write header in ASCII format 
379
380
381
382
383
384
385
386
387
  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
388
    print*,'call openreceptors'
389
  endif
Matthias Langer's avatar
 
Matthias Langer committed
390
391
392
393
394
395
  call openreceptors
  if ((iout.eq.4).or.(iout.eq.5)) call openouttraj

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

396
  if (verbosity.gt.0) then
397
    print*,'discretize release times'
398
  endif
Matthias Langer's avatar
 
Matthias Langer committed
399
  do i=1,numpoint
400
401
    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
402
403
404
405
406
  end do

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

407
  if (verbosity.gt.0) then
408
    print*,'Initialize cloud-base mass fluxes for the convection scheme'
409
410
  endif

Matthias Langer's avatar
 
Matthias Langer committed
411
412
413
414
415
416
417
418
419
420
421
422
423
  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

424
425
426
  ! Inform whether output kernel is used or not
  !*********************************************
  if (lroot) then
427
    if (.not.lusekerneloutput) then
428
429
430
431
432
433
434
      write(*,*) "Concentrations are calculated without using kernel"
    else
      write(*,*) "Concentrations are calculated using kernel"
    end if
  end if


Matthias Langer's avatar
 
Matthias Langer committed
435
436
437
  ! Calculate particle trajectories
  !********************************

438
  if (verbosity.gt.0) then
439
440
441
442
443
444
445
446
447
     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'
448
449
  endif

450
  if (verbosity.gt.0) write (*,*) 'timemanager> call wetdepo'
451
  call timemanager(metdata_format)
452
 
Matthias Langer's avatar
 
Matthias Langer committed
453

454
  if (verbosity.gt.0) then
455
! NIK 16.02.2005 
456
457
458
459
460
461
462
463
464
465
466
467
    do i=1,nspec
      if (tot_inc_count(i).gt.0) then
         write(*,*) '**********************************************'
         write(*,*) 'Scavenging statistics for species ', species(i), ':'
         write(*,*) 'Total number of occurences of below-cloud scavenging', &
           & tot_blc_count(i)
         write(*,*) 'Total number of occurences of in-cloud    scavenging', &
           & tot_inc_count(i)
         write(*,*) '**********************************************'
      endif
    end do
  endif
468
  
469
470
  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
       &XPART MODEL RUN!'
Matthias Langer's avatar
 
Matthias Langer committed
471
472

end program flexpart