FLEXPART.f90 14.6 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 netcdf_output_mod, only: writeheader_netcdf
  use random_mod, only: gasdev1
56
  use class_gribfile
Matthias Langer's avatar
 
Matthias Langer committed
57
58
59
60
61

  implicit none

  integer :: i,j,ix,jy,inest
  integer :: idummy = -320
62
  character(len=256) :: inline_options  !pathfile, flexversion, arg2
63
64
65
  integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN
  integer :: detectformat

66
67
68
69


  ! Initialize arrays in com_mod
  !*****************************
70
71
  call com_mod_allocate_part(maxpart)
  
72

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

81

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

87
88
89
  ! Read the pathnames where input/output files are stored
  !*******************************************************

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

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

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

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

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

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

181
  if (verbosity.gt.0) then
182
    write(*,*) 'call readavailable'
183
  endif  
Matthias Langer's avatar
 
Matthias Langer committed
184
185
  call readavailable

186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
  ! 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'
    return
  endif



202
203
204
205
206
207
208
209
  ! 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
210
211
212
  ! Read the model grid specifications,
  ! both for the mother domain and eventual nests
  !**********************************************
213
214
 
  if (verbosity.gt.0) then
215
     write(*,*) 'call gridcheck'
216
  endif
217

218
219
220
221
222
  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
    call gridcheck_ecmwf
  else
    call gridcheck_gfs
  end if
223

224
  if (verbosity.gt.1) then   
225
226
    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
227
228
229
  endif      

  if (verbosity.gt.0) then
230
    write(*,*) 'call gridcheck_nests'
231
  endif  
Matthias Langer's avatar
 
Matthias Langer committed
232
233
234
235
236
  call gridcheck_nests

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

237
  if (verbosity.gt.0) then
238
    write(*,*) 'call readoutgrid'
239
240
  endif

Matthias Langer's avatar
 
Matthias Langer committed
241
242
  call readoutgrid

243
  if (nested_output.eq.1) then
244
    call readoutgrid_nest
245
    if (verbosity.gt.0) then
246
      write(*,*) '# readoutgrid_nest'
247
248
    endif
  endif
Matthias Langer's avatar
 
Matthias Langer committed
249
250
251
252

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

253
  if (verbosity.eq.1) then
254
     print*,'call readreceptors'
255
  endif
Matthias Langer's avatar
 
Matthias Langer committed
256
257
258
259
260
261
262
263
264
265
266
  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
  !***************************

267
  if (verbosity.gt.0) then
268
    print*,'call readlanduse'
269
  endif
Matthias Langer's avatar
 
Matthias Langer committed
270
271
272
273
274
  call readlanduse

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

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

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

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

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

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

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

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

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

307
  if (verbosity.gt.0) then
308
    print*,'Initialize all particles to non-existent'
309
  endif
Matthias Langer's avatar
 
Matthias Langer committed
310
311
312
313
314
315
316
317
  do j=1,maxpart
    itra1(j)=-999999999
  end do

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

  if (ipin.eq.1) then
318
    if (verbosity.gt.0) then
319
      print*,'call readpartpositions'
320
    endif
Matthias Langer's avatar
 
Matthias Langer committed
321
322
    call readpartpositions
  else
323
    if (verbosity.gt.0) then
324
      print*,'numpart=0, numparticlecount=0'
325
    endif    
Matthias Langer's avatar
 
Matthias Langer committed
326
327
328
329
330
331
332
333
    numpart=0
    numparticlecount=0
  endif

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

334
  if (verbosity.gt.0) then
335
    print*,'call outgrid_init'
336
  endif
Matthias Langer's avatar
 
Matthias Langer committed
337
338
339
340
341
342
  call outgrid_init
  if (nested_output.eq.1) call outgrid_init_nest

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

343
344
  if (OHREA.eqv..TRUE.) then
    if (verbosity.gt.0) then
345
      print*,'call readOHfield'
346
    endif
347
    call readOHfield
348
  endif
Matthias Langer's avatar
 
Matthias Langer committed
349
350
351
352
353

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

354
355
356
357
358
359
360
361
362
363
364
365
366
367
  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

368
  if (verbosity.gt.0) then
369
    print*,'call writeheader'
370
371
  endif

Matthias Langer's avatar
 
Matthias Langer committed
372
  call writeheader
373
  ! FLEXPART 9.2 ticket ?? write header in ASCII format 
374
375
376
377
378
379
380
381
382
  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
383
    print*,'call openreceptors'
384
  endif
Matthias Langer's avatar
 
Matthias Langer committed
385
386
387
388
389
390
  call openreceptors
  if ((iout.eq.4).or.(iout.eq.5)) call openouttraj

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

391
  if (verbosity.gt.0) then
392
    print*,'discretize release times'
393
  endif
Matthias Langer's avatar
 
Matthias Langer committed
394
  do i=1,numpoint
395
396
    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
397
398
399
400
401
  end do

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

402
  if (verbosity.gt.0) then
403
    print*,'Initialize cloud-base mass fluxes for the convection scheme'
404
405
  endif

Matthias Langer's avatar
 
Matthias Langer committed
406
407
408
409
410
411
412
413
414
415
416
417
418
  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

419
420
421
422
423
424
425
426
427
428
429
  ! Inform whether output kernel is used or not
  !*********************************************
  if (lroot) then
    if (lnokernel) then
      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
430
431
432
  ! Calculate particle trajectories
  !********************************

433
  if (verbosity.gt.0) then
434
435
436
437
438
439
440
441
442
     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'
443
444
  endif

445
  call timemanager(metdata_format)
Matthias Langer's avatar
 
Matthias Langer committed
446

447
! NIK 16.02.2005 
448
449
450
451
452
453
454
455
456
457
  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
  
458
459
  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
       &XPART MODEL RUN!'
Matthias Langer's avatar
 
Matthias Langer committed
460
461

end program flexpart