FLEXPART.f90 15 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
!**********************************************************************
! 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                                                            *
  !                                                                            *
  !*****************************************************************************
35
36
37
38
39
40
41
42
  ! 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
43
44
45
46
47
48
49
50
51
52
53
  !                                                                            *
  ! Variables:                                                                 *
  !                                                                            *
  ! Constants:                                                                 *
  !                                                                            *
  !*****************************************************************************

  use point_mod
  use par_mod
  use com_mod
  use conv_mod
54

55
  use random_mod, only: gasdev1
56
  use class_gribfile
Matthias Langer's avatar
 
Matthias Langer committed
57

58
59
60
61
#ifdef USE_NCF
  use netcdf_output_mod, only: writeheader_netcdf
#endif

Matthias Langer's avatar
 
Matthias Langer committed
62
63
64
65
  implicit none

  integer :: i,j,ix,jy,inest
  integer :: idummy = -320
66
  character(len=256) :: inline_options  !pathfile, flexversion, arg2
67
68
69
  integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN
  integer :: detectformat

70
  
Matthias Langer's avatar
 
Matthias Langer committed
71
72
73
74
75
76
77
78
  ! 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))

79

80
  ! FLEXPART version string
81
  flexversion_major = '10' ! Major version number, also used for species file names
Ignacio Pisso's avatar
Ignacio Pisso committed
82
  flexversion='Version '//trim(flexversion_major)//'.3 (2019-06-26)'
83
84
  verbosity=0

85
86
87
  ! Read the pathnames where input/output files are stored
  !*******************************************************

88
  inline_options='none'
89
  select case (iargc())
90
  case (2)
91
92
93
94
    call getarg(1,arg1)
    pathfile=arg1
    call getarg(2,arg2)
    inline_options=arg2
95
  case (1)
96
97
98
    call getarg(1,arg1)
    pathfile=arg1
    if (arg1(1:1).eq.'-') then
99
100
      write(pathfile,'(a11)') './pathnames'
      inline_options=arg1 
101
    endif
102
  case (0)
103
104
105
    write(pathfile,'(a11)') './pathnames'
  end select
  
Matthias Langer's avatar
 
Matthias Langer committed
106
107
  ! Print the GPL License statement
  !*******************************************************
108
  print*,'Welcome to FLEXPART ', trim(flexversion)
109
  print*,'FLEXPART is free software released under the GNU General Public License.'
110
 
111
112
113
  if (inline_options(1:1).eq.'-') then
    if (trim(inline_options).eq.'-v'.or.trim(inline_options).eq.'-v1') then
       print*, 'Verbose mode 1: display detailed information during run'
114
       verbosity=1
115
116
117
    endif
    if (trim(inline_options).eq.'-v2') then
       print*, 'Verbose mode 2: display more detailed information during run'
118
119
       verbosity=2
    endif
120
121
    if (trim(inline_options).eq.'-i') then
       print*, 'Info mode: provide detailed run specific information and stop'
122
123
124
       verbosity=1
       info_flag=1
    endif
125
126
127
128
    if (trim(inline_options).eq.'-i2') then
       print*, 'Info mode: provide more detailed run specific information and stop'
       verbosity=2
       info_flag=1
129
130
131
    endif
  endif
           
132
  if (verbosity.gt.0) then
Ignacio Pisso's avatar
Ignacio Pisso committed
133
134
135
136
    print*, 'nxmax=',nxmax
    print*, 'nymax=',nymax
    print*, 'nzmax=',nzmax
    print*,'nxshift=',nxshift 
137
    write(*,*) 'call readpaths'
138
  endif 
139
  call readpaths
140
 
141
142
  if (verbosity.gt.1) then !show clock info 
     !print*,'length(4)',length(4)
143
     !count=0,count_rate=1000
144
     CALL SYSTEM_CLOCK(count_clock0, count_rate, count_max)
145
146
147
148
     !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
149
  endif
Matthias Langer's avatar
 
Matthias Langer committed
150
151
152
153

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

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

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

Matthias Langer's avatar
 
Matthias Langer committed
172
173
174

  ! Read the age classes to be used
  !********************************
175
  if (verbosity.gt.0) then
176
    write(*,*) 'call readageclasses'
177
  endif
Matthias Langer's avatar
 
Matthias Langer committed
178
179
  call readageclasses

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

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

188
  if (verbosity.gt.0) then
189
    write(*,*) 'call readavailable'
190
  endif  
Matthias Langer's avatar
 
Matthias Langer committed
191
192
  call readavailable

193
194
195
196
197
198
199
200
201
202
203
  ! Detect metdata format
  !**********************

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



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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

361
#ifdef USE_NCF
362
363
364
365
366
367
368
369
370
371
372
373
374
  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
375
#endif
376

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

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

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

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

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

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

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

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

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

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

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

end program flexpart