FLEXPART.f90 15.7 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
125
    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*, 'verbosity=2' 
        print*, 'Verbose mode 2: display more detailed information during run'
        verbosity=2
      endif
126
    endif
127
128
129
130
131
132
133

    
    !debug mode 
    iopt=index(inline_options,'d')
    if (iopt.gt.0) then
      debug_mode=.true.
      endif
134
    endif
135
136
137
138
139
140
141
142
143
144
145
146
147


    
!      stop
!    if (trim(inline_options).eq.'-v'.or.trim(inline_options).eq.'-v1') then
!       print*, 'Verbose mode 1: display detailed information during run'
!       verbosity=1
!    endif
!    if (trim(inline_options).eq.'-v2') then
!       print*, 'Verbose mode 2: display more detailed information during run'
!       verbosity=2
!    endif

148
149
    if (trim(inline_options).eq.'-i') then
       print*, 'Info mode: provide detailed run specific information and stop'
150
151
152
       verbosity=1
       info_flag=1
    endif
153
154
155
156
    if (trim(inline_options).eq.'-i2') then
       print*, 'Info mode: provide more detailed run specific information and stop'
       verbosity=2
       info_flag=1
157
158
159
    endif
  endif
           
160
  if (verbosity.gt.0) then
Ignacio Pisso's avatar
Ignacio Pisso committed
161
162
163
164
    print*, 'nxmax=',nxmax
    print*, 'nymax=',nymax
    print*, 'nzmax=',nzmax
    print*,'nxshift=',nxshift 
165
166
167
  endif
  
  if (verbosity.gt.0) then
168
    write(*,*) 'call readpaths'
169
  endif 
170
  call readpaths
171
 
172
173
  if (verbosity.gt.1) then !show clock info 
     !print*,'length(4)',length(4)
174
     !count=0,count_rate=1000
175
     CALL SYSTEM_CLOCK(count_clock0, count_rate, count_max)
176
177
178
179
     !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
180
  endif
Matthias Langer's avatar
 
Matthias Langer committed
181
182
183
184

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

185
  if (verbosity.gt.0) then
186
    write(*,*) 'call readcommand'
187
  endif
Matthias Langer's avatar
 
Matthias Langer committed
188
  call readcommand
189
  if (verbosity.gt.0) then
190
191
    write(*,*) '    ldirect=', ldirect
    write(*,*) '    ibdate,ibtime=',ibdate,ibtime
192
    write(*,*) '    iedate,ietime=', iedate,ietime
193
    if (verbosity.gt.1) then   
194
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
    endif     
197
  endif
198
199
200
201
202

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

Matthias Langer's avatar
 
Matthias Langer committed
203
204
205

  ! Read the age classes to be used
  !********************************
206
  if (verbosity.gt.0) then
207
    write(*,*) 'call readageclasses'
208
  endif
Matthias Langer's avatar
 
Matthias Langer committed
209
210
  call readageclasses

211
  if (verbosity.gt.1) then   
212
213
    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
214
  endif     
Matthias Langer's avatar
 
Matthias Langer committed
215
216
217
218

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

219
  if (verbosity.gt.0) then
220
    write(*,*) 'call readavailable'
221
  endif  
Matthias Langer's avatar
 
Matthias Langer committed
222
223
  call readavailable

224
225
  ! Detect metdata format
  !**********************
226
227
228
  if (verbosity.gt.0) then
    write(*,*) 'call detectformat'
  endif
229
230
231
232
233
234
235
236
237

  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'
238
    stop
239
240
241
242
  endif



243
244
245
246
247
248
249
250
  ! 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
251
252
253
  ! Read the model grid specifications,
  ! both for the mother domain and eventual nests
  !**********************************************
254
255
 
  if (verbosity.gt.0) then
256
     write(*,*) 'call gridcheck'
257
  endif
258

259
260
261
262
263
  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
    call gridcheck_ecmwf
  else
    call gridcheck_gfs
  end if
264

265
  if (verbosity.gt.1) then   
266
267
    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
268
269
270
  endif      

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

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

278
  if (verbosity.gt.0) then
279
    write(*,*) 'call readoutgrid'
280
281
  endif

Matthias Langer's avatar
 
Matthias Langer committed
282
283
  call readoutgrid

284
  if (nested_output.eq.1) then
285
    call readoutgrid_nest
286
    if (verbosity.gt.0) then
287
      write(*,*) '# readoutgrid_nest'
288
289
    endif
  endif
Matthias Langer's avatar
 
Matthias Langer committed
290
291
292
293

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

294
  if (verbosity.eq.1) then
295
     print*,'call readreceptors'
296
  endif
Matthias Langer's avatar
 
Matthias Langer committed
297
298
299
300
301
302
303
304
305
306
307
  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
  !***************************

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

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

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

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

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

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

332
  if (verbosity.gt.0) then
333
    print*,'call readdepo'
334
  endif
Matthias Langer's avatar
 
Matthias Langer committed
335
336
337
338
339
  call readdepo

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

340
341
  call coordtrafo  
  if (verbosity.gt.0) then
342
    print*,'call coordtrafo'
343
  endif
Matthias Langer's avatar
 
Matthias Langer committed
344
345
346
347

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

348
  if (verbosity.gt.0) then
349
    print*,'Initialize all particles to non-existent'
350
  endif
Matthias Langer's avatar
 
Matthias Langer committed
351
352
353
354
355
356
357
358
  do j=1,maxpart
    itra1(j)=-999999999
  end do

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

  if (ipin.eq.1) then
359
    if (verbosity.gt.0) then
360
      print*,'call readpartpositions'
361
    endif
Matthias Langer's avatar
 
Matthias Langer committed
362
363
    call readpartpositions
  else
364
    if (verbosity.gt.0) then
365
      print*,'numpart=0, numparticlecount=0'
366
    endif    
Matthias Langer's avatar
 
Matthias Langer committed
367
368
369
370
371
372
373
374
    numpart=0
    numparticlecount=0
  endif

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

375
  if (verbosity.gt.0) then
376
    print*,'call outgrid_init'
377
  endif
Matthias Langer's avatar
 
Matthias Langer committed
378
379
380
381
382
383
  call outgrid_init
  if (nested_output.eq.1) call outgrid_init_nest

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

384
385
  if (OHREA.eqv..TRUE.) then
    if (verbosity.gt.0) then
386
      print*,'call readOHfield'
387
    endif
388
    call readOHfield
389
  endif
Matthias Langer's avatar
 
Matthias Langer committed
390
391
392
393
394

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

395
#ifdef USE_NCF
396
397
398
399
400
401
402
403
404
405
406
407
408
  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
409
#endif
410

411
  if (verbosity.gt.0) then
412
    print*,'call writeheader'
413
414
  endif

Matthias Langer's avatar
 
Matthias Langer committed
415
  call writeheader
416
  ! FLEXPART 9.2 ticket ?? write header in ASCII format 
417
418
419
420
421
422
423
424
425
  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
426
    print*,'call openreceptors'
427
  endif
Matthias Langer's avatar
 
Matthias Langer committed
428
429
430
431
432
433
  call openreceptors
  if ((iout.eq.4).or.(iout.eq.5)) call openouttraj

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

434
  if (verbosity.gt.0) then
435
    print*,'discretize release times'
436
  endif
Matthias Langer's avatar
 
Matthias Langer committed
437
  do i=1,numpoint
438
439
    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
440
441
442
443
444
  end do

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

445
  if (verbosity.gt.0) then
446
    print*,'Initialize cloud-base mass fluxes for the convection scheme'
447
448
  endif

Matthias Langer's avatar
 
Matthias Langer committed
449
450
451
452
453
454
455
456
457
458
459
460
461
  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

462
463
464
  ! Inform whether output kernel is used or not
  !*********************************************
  if (lroot) then
465
    if (.not.lusekerneloutput) then
466
467
468
469
470
471
472
      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
473
474
475
  ! Calculate particle trajectories
  !********************************

476
  if (verbosity.gt.0) then
477
478
479
480
481
482
483
484
485
     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'
486
487
  endif

488
  if (verbosity.gt.0) write (*,*) 'timemanager> call wetdepo'
489
  call timemanager(metdata_format)
490
 
Matthias Langer's avatar
 
Matthias Langer committed
491

492
  if (verbosity.gt.0) then
493
! NIK 16.02.2005 
494
495
496
497
498
499
500
501
502
503
504
505
    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
506
  
507
508
  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
       &XPART MODEL RUN!'
Matthias Langer's avatar
 
Matthias Langer committed
509
510

end program flexpart