FLEXPART.f90 15.3 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
  implicit none

64
  integer :: i,j,ix,jy,inest, iopt
Matthias Langer's avatar
 
Matthias Langer committed
65
  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
82
  flexversion='Version '//trim(flexversion_major)//'.4 (2019-07-23)'
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

  ! Ingest inline options
  !*******************************************************
114
  if (inline_options(1:1).eq.'-') then
115
116
117
118
119
120
121
122
123
124
    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
125
    endif
126
127
128
129
    !debug mode 
    iopt=index(inline_options,'d')
    if (iopt.gt.0) then
      debug_mode=.true.
130
    endif
131
132
    if (trim(inline_options).eq.'-i') then
       print*, 'Info mode: provide detailed run specific information and stop'
133
134
135
       verbosity=1
       info_flag=1
    endif
136
137
138
139
    if (trim(inline_options).eq.'-i2') then
       print*, 'Info mode: provide more detailed run specific information and stop'
       verbosity=2
       info_flag=1
140
141
142
    endif
  endif
           
143
  if (verbosity.gt.0) then
Ignacio Pisso's avatar
Ignacio Pisso committed
144
145
146
147
    print*, 'nxmax=',nxmax
    print*, 'nymax=',nymax
    print*, 'nzmax=',nzmax
    print*,'nxshift=',nxshift 
148
149
150
  endif
  
  if (verbosity.gt.0) then
151
    write(*,*) 'call readpaths'
152
  endif 
153
  call readpaths
154
 
155
156
  if (verbosity.gt.1) then !show clock info 
     !print*,'length(4)',length(4)
157
     !count=0,count_rate=1000
158
     CALL SYSTEM_CLOCK(count_clock0, count_rate, count_max)
159
160
161
162
     !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
163
  endif
Matthias Langer's avatar
 
Matthias Langer committed
164
165
166
167

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

168
  if (verbosity.gt.0) then
169
    write(*,*) 'call readcommand'
170
  endif
Matthias Langer's avatar
 
Matthias Langer committed
171
  call readcommand
172
  if (verbosity.gt.0) then
173
174
    write(*,*) '    ldirect=', ldirect
    write(*,*) '    ibdate,ibtime=',ibdate,ibtime
175
    write(*,*) '    iedate,ietime=', iedate,ietime
176
    if (verbosity.gt.1) then   
177
178
179
      CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
      write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
    endif     
180
  endif
181
182
183
184
185

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

Matthias Langer's avatar
 
Matthias Langer committed
186
187
188

  ! Read the age classes to be used
  !********************************
189
  if (verbosity.gt.0) then
190
    write(*,*) 'call readageclasses'
191
  endif
Matthias Langer's avatar
 
Matthias Langer committed
192
193
  call readageclasses

194
  if (verbosity.gt.1) then   
195
196
    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
197
  endif     
Matthias Langer's avatar
 
Matthias Langer committed
198
199
200
201

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

202
  if (verbosity.gt.0) then
203
    write(*,*) 'call readavailable'
204
  endif  
Matthias Langer's avatar
 
Matthias Langer committed
205
206
  call readavailable

207
208
  ! Detect metdata format
  !**********************
209
210
211
  if (verbosity.gt.0) then
    write(*,*) 'call detectformat'
  endif
212
213
214
215
216
217
218
219
220

  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'
221
    stop
222
223
224
225
  endif



226
227
228
229
230
231
232
233
  ! 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
234
235
236
  ! Read the model grid specifications,
  ! both for the mother domain and eventual nests
  !**********************************************
237
238
 
  if (verbosity.gt.0) then
239
     write(*,*) 'call gridcheck'
240
  endif
241

242
243
244
245
246
  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
    call gridcheck_ecmwf
  else
    call gridcheck_gfs
  end if
247

248
  if (verbosity.gt.1) then   
249
250
    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
251
252
253
  endif      

  if (verbosity.gt.0) then
254
    write(*,*) 'call gridcheck_nests'
255
  endif  
Matthias Langer's avatar
 
Matthias Langer committed
256
257
258
259
260
  call gridcheck_nests

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

261
  if (verbosity.gt.0) then
262
    write(*,*) 'call readoutgrid'
263
264
  endif

Matthias Langer's avatar
 
Matthias Langer committed
265
266
  call readoutgrid

267
  if (nested_output.eq.1) then
268
    call readoutgrid_nest
269
    if (verbosity.gt.0) then
270
      write(*,*) '# readoutgrid_nest'
271
272
    endif
  endif
Matthias Langer's avatar
 
Matthias Langer committed
273
274
275
276

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

277
  if (verbosity.eq.1) then
278
     print*,'call readreceptors'
279
  endif
Matthias Langer's avatar
 
Matthias Langer committed
280
281
282
283
284
285
286
287
288
289
290
  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
  !***************************

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

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

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

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

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

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

315
  if (verbosity.gt.0) then
316
    print*,'call readdepo'
317
  endif
Matthias Langer's avatar
 
Matthias Langer committed
318
319
320
321
322
  call readdepo

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

323
324
  call coordtrafo  
  if (verbosity.gt.0) then
325
    print*,'call coordtrafo'
326
  endif
Matthias Langer's avatar
 
Matthias Langer committed
327
328
329
330

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

331
  if (verbosity.gt.0) then
332
    print*,'Initialize all particles to non-existent'
333
  endif
Matthias Langer's avatar
 
Matthias Langer committed
334
335
336
337
338
339
340
341
  do j=1,maxpart
    itra1(j)=-999999999
  end do

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

  if (ipin.eq.1) then
342
    if (verbosity.gt.0) then
343
      print*,'call readpartpositions'
344
    endif
Matthias Langer's avatar
 
Matthias Langer committed
345
346
    call readpartpositions
  else
347
    if (verbosity.gt.0) then
348
      print*,'numpart=0, numparticlecount=0'
349
    endif    
Matthias Langer's avatar
 
Matthias Langer committed
350
351
352
353
354
355
356
357
    numpart=0
    numparticlecount=0
  endif

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

358
  if (verbosity.gt.0) then
359
    print*,'call outgrid_init'
360
  endif
Matthias Langer's avatar
 
Matthias Langer committed
361
362
363
364
365
366
  call outgrid_init
  if (nested_output.eq.1) call outgrid_init_nest

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

367
368
  if (OHREA.eqv..TRUE.) then
    if (verbosity.gt.0) then
369
      print*,'call readOHfield'
370
    endif
371
    call readOHfield
372
  endif
Matthias Langer's avatar
 
Matthias Langer committed
373
374
375
376
377

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

378
#ifdef USE_NCF
379
380
381
382
383
384
385
386
387
388
389
390
391
  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
392
#endif
393

394
  if (verbosity.gt.0) then
395
    print*,'call writeheader'
396
397
  endif

Matthias Langer's avatar
 
Matthias Langer committed
398
  call writeheader
399
  ! FLEXPART 9.2 ticket ?? write header in ASCII format 
400
401
402
403
404
405
406
407
408
  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
409
    print*,'call openreceptors'
410
  endif
Matthias Langer's avatar
 
Matthias Langer committed
411
412
413
414
415
416
  call openreceptors
  if ((iout.eq.4).or.(iout.eq.5)) call openouttraj

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

417
  if (verbosity.gt.0) then
418
    print*,'discretize release times'
419
  endif
Matthias Langer's avatar
 
Matthias Langer committed
420
  do i=1,numpoint
421
422
    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
423
424
425
426
427
  end do

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

428
  if (verbosity.gt.0) then
429
    print*,'Initialize cloud-base mass fluxes for the convection scheme'
430
431
  endif

Matthias Langer's avatar
 
Matthias Langer committed
432
433
434
435
436
437
438
439
440
441
442
443
444
  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

445
446
447
  ! Inform whether output kernel is used or not
  !*********************************************
  if (lroot) then
448
    if (.not.lusekerneloutput) then
449
450
451
452
453
454
455
      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
456
457
458
  ! Calculate particle trajectories
  !********************************

459
  if (verbosity.gt.0) then
460
461
462
463
464
465
466
467
468
     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'
469
470
  endif

471
  if (verbosity.gt.0) write (*,*) 'timemanager> call wetdepo'
472
  call timemanager(metdata_format)
473
 
Matthias Langer's avatar
 
Matthias Langer committed
474

475
  if (verbosity.gt.0) then
476
! NIK 16.02.2005 
477
478
479
480
481
482
483
484
485
486
487
488
    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
489
  
490
491
  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
       &XPART MODEL RUN!'
Matthias Langer's avatar
 
Matthias Langer committed
492
493

end program flexpart