FLEXPART.f90 14.8 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
71
72
73


  ! Initialize arrays in com_mod
  !*****************************
74
75
  call com_mod_allocate_part(maxpart)
  
76

Matthias Langer's avatar
 
Matthias Langer committed
77
78
79
80
81
82
83
84
  ! 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))

85

86
  ! FLEXPART version string
87
  flexversion_major = '10' ! Major version number, also used for species file names
Espen Sollum's avatar
Espen Sollum committed
88
  flexversion='Version '//trim(flexversion_major)//'.2beta (2017-08-01)'
89
90
  verbosity=0

91
92
93
  ! Read the pathnames where input/output files are stored
  !*******************************************************

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

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

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

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

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

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

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

194
195
196
197
198
199
200
201
202
203
204
  ! 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'
205
    stop
206
207
208
209
  endif



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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

455
  call timemanager(metdata_format)
Matthias Langer's avatar
 
Matthias Langer committed
456

457
! NIK 16.02.2005 
458
459
460
461
462
463
464
465
466
467
  do i=1,nspec
    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(*,*) '**********************************************'
  end do
  
468
469
  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
       &XPART MODEL RUN!'
Matthias Langer's avatar
 
Matthias Langer committed
470
471

end program flexpart