FLEXPART.f90 14 KB
Newer Older
1
2
! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
! SPDX-License-Identifier: GPL-3.0-or-later
3

Matthias Langer's avatar
 
Matthias Langer committed
4
5
6
7
8
9
10
11
12
13
14
15
16
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                                                            *
  !                                                                            *
  !*****************************************************************************
17
18
19
20
21
22
23
24
  ! 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
25
26
27
28
29
30
31
32
33
34
35
  !                                                                            *
  ! Variables:                                                                 *
  !                                                                            *
  ! Constants:                                                                 *
  !                                                                            *
  !*****************************************************************************

  use point_mod
  use par_mod
  use com_mod
  use conv_mod
36

37
  use random_mod, only: gasdev1
38
  use class_gribfile
Matthias Langer's avatar
 
Matthias Langer committed
39

40
41
42
43
#ifdef USE_NCF
  use netcdf_output_mod, only: writeheader_netcdf
#endif

Matthias Langer's avatar
 
Matthias Langer committed
44
45
  implicit none

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

52
  
Matthias Langer's avatar
 
Matthias Langer committed
53
54
55
56
57
58
59
60
  ! 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))

61

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

67
68
69
  ! Read the pathnames where input/output files are stored
  !*******************************************************

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

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

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

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

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

Matthias Langer's avatar
 
Matthias Langer committed
168
169
170

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

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

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

184
  if (verbosity.gt.0) then
185
    write(*,*) 'call readavailable'
186
  endif  
Matthias Langer's avatar
 
Matthias Langer committed
187
188
  call readavailable

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

  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'
203
    stop
204
205
206
207
  endif



208
209
210
211
212
213
214
215
  ! 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
216
217
218
  ! Read the model grid specifications,
  ! both for the mother domain and eventual nests
  !**********************************************
219
220
 
  if (verbosity.gt.0) then
221
     write(*,*) 'call gridcheck'
222
  endif
223

224
225
226
227
228
  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
    call gridcheck_ecmwf
  else
    call gridcheck_gfs
  end if
229

230
  if (verbosity.gt.1) 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(*,*) '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(*,*) '# 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*,'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*,'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*,'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*,'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
  if (verbosity.gt.0) then
298
    print*,'call readdepo'
299
  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
  call coordtrafo  
  if (verbosity.gt.0) then
307
    print*,'call coordtrafo'
308
  endif
Matthias Langer's avatar
 
Matthias Langer committed
309
310
311
312

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

313
  if (verbosity.gt.0) then
314
    print*,'Initialize all particles to non-existent'
315
  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
    if (verbosity.gt.0) then
325
      print*,'call readpartpositions'
326
    endif
Matthias Langer's avatar
 
Matthias Langer committed
327
328
    call readpartpositions
  else
329
    if (verbosity.gt.0) then
330
      print*,'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
  if (verbosity.gt.0) then
341
    print*,'call outgrid_init'
342
  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
#ifdef USE_NCF
361
362
363
364
365
366
367
368
369
370
371
372
373
  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
374
#endif
375

376
  if (verbosity.gt.0) then
377
    print*,'call writeheader'
378
379
  endif

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

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

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

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

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

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

427
428
429
  ! Inform whether output kernel is used or not
  !*********************************************
  if (lroot) then
430
    if (.not.lusekerneloutput) then
431
432
433
434
435
436
437
      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
438
439
440
  ! Calculate particle trajectories
  !********************************

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

453
  if (verbosity.gt.0) write (*,*) 'timemanager> call wetdepo'
454
  call timemanager(metdata_format)
455
 
Matthias Langer's avatar
 
Matthias Langer committed
456

457
  if (verbosity.gt.0) then
458
! NIK 16.02.2005 
459
460
461
462
463
464
465
466
467
468
469
470
    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
471
  
472
473
  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
       &XPART MODEL RUN!'
Matthias Langer's avatar
 
Matthias Langer committed
474
475

end program flexpart