FLEXPART.f90 13.5 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
35
36
37
38
39
40
41
42
43
44
45
!**********************************************************************
! 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                                                            *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  !                                                                            *
  ! Constants:                                                                 *
  !                                                                            *
  !*****************************************************************************

  use point_mod
  use par_mod
  use com_mod
  use conv_mod
46
47
  use netcdf_output_mod, only: writeheader_netcdf
  use random_mod, only: gasdev1
Matthias Langer's avatar
 
Matthias Langer committed
48
49
50
51
52

  implicit none

  integer :: i,j,ix,jy,inest
  integer :: idummy = -320
53
  character(len=256) :: inline_options  !pathfile, flexversion, arg2
54
55
56
57


  ! Initialize arrays in com_mod
  !*****************************
58
59
  call com_mod_allocate_part(maxpart)
  
60

Matthias Langer's avatar
 
Matthias Langer committed
61
62
63
64
65
66
67
68
  ! 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))

69

70
  ! FLEXPART version string
71
  flexversion_major = '10' ! Major version number, also used for species file names
72
  flexversion='Version '//trim(flexversion_major)//'.1beta (2016-11-02)'
73
74
  verbosity=0

75
76
77
  ! Read the pathnames where input/output files are stored
  !*******************************************************

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

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

140
  if (verbosity.gt.0) then
141
    write(*,*) 'call readcommand'
142
  endif
Matthias Langer's avatar
 
Matthias Langer committed
143
  call readcommand
144
  if (verbosity.gt.0) then
145
146
    write(*,*) '    ldirect=', ldirect
    write(*,*) '    ibdate,ibtime=',ibdate,ibtime
147
    write(*,*) '    iedate,ietime=', iedate,ietime
148
    if (verbosity.gt.1) then   
149
150
151
      CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
      write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
    endif     
152
  endif
Matthias Langer's avatar
 
Matthias Langer committed
153
154
155

  ! Read the age classes to be used
  !********************************
156
  if (verbosity.gt.0) then
157
    write(*,*) 'call readageclasses'
158
  endif
Matthias Langer's avatar
 
Matthias Langer committed
159
160
  call readageclasses

161
  if (verbosity.gt.1) then   
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
164
  endif     
Matthias Langer's avatar
 
Matthias Langer committed
165
166
167
168

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

169
  if (verbosity.gt.0) then
170
    write(*,*) 'call readavailable'
171
  endif  
Matthias Langer's avatar
 
Matthias Langer committed
172
173
  call readavailable

174
175
176
177
178
179
180
181
  ! 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
182
183
184
  ! Read the model grid specifications,
  ! both for the mother domain and eventual nests
  !**********************************************
185
186
 
  if (verbosity.gt.0) then
187
     write(*,*) 'call gridcheck'
188
  endif
189

Matthias Langer's avatar
 
Matthias Langer committed
190
  call gridcheck
191

192
  if (verbosity.gt.1) then   
193
194
    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
195
196
197
  endif      

  if (verbosity.gt.0) then
198
    write(*,*) 'call gridcheck_nests'
199
  endif  
Matthias Langer's avatar
 
Matthias Langer committed
200
201
202
203
204
  call gridcheck_nests

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

205
  if (verbosity.gt.0) then
206
    write(*,*) 'call readoutgrid'
207
208
  endif

Matthias Langer's avatar
 
Matthias Langer committed
209
210
  call readoutgrid

211
  if (nested_output.eq.1) then
212
    call readoutgrid_nest
213
    if (verbosity.gt.0) then
214
      write(*,*) '# readoutgrid_nest'
215
216
    endif
  endif
Matthias Langer's avatar
 
Matthias Langer committed
217
218
219
220

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

221
  if (verbosity.eq.1) then
222
     print*,'call readreceptors'
223
  endif
Matthias Langer's avatar
 
Matthias Langer committed
224
225
226
227
228
229
230
231
232
233
234
  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
  !***************************

235
  if (verbosity.gt.0) then
236
    print*,'call readlanduse'
237
  endif
Matthias Langer's avatar
 
Matthias Langer committed
238
239
240
241
242
  call readlanduse

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

243
  if (verbosity.gt.0) then
244
    print*,'call assignland'
245
  endif
Matthias Langer's avatar
 
Matthias Langer committed
246
247
248
249
250
  call assignland

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

251
  if (verbosity.gt.0) then
252
    print*,'call readreleases'
253
  endif
Matthias Langer's avatar
 
Matthias Langer committed
254
255
256
257
258
  call readreleases

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

259
  if (verbosity.gt.0) then
260
    print*,'call readdepo'
261
  endif
Matthias Langer's avatar
 
Matthias Langer committed
262
263
264
265
266
  call readdepo

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

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

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

275
  if (verbosity.gt.0) then
276
    print*,'Initialize all particles to non-existent'
277
  endif
Matthias Langer's avatar
 
Matthias Langer committed
278
279
280
281
282
283
284
285
  do j=1,maxpart
    itra1(j)=-999999999
  end do

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

  if (ipin.eq.1) then
286
    if (verbosity.gt.0) then
287
      print*,'call readpartpositions'
288
    endif
Matthias Langer's avatar
 
Matthias Langer committed
289
290
    call readpartpositions
  else
291
    if (verbosity.gt.0) then
292
      print*,'numpart=0, numparticlecount=0'
293
    endif    
Matthias Langer's avatar
 
Matthias Langer committed
294
295
296
297
298
299
300
301
    numpart=0
    numparticlecount=0
  endif

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

302
  if (verbosity.gt.0) then
303
    print*,'call outgrid_init'
304
  endif
Matthias Langer's avatar
 
Matthias Langer committed
305
306
307
308
309
310
  call outgrid_init
  if (nested_output.eq.1) call outgrid_init_nest

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

311
312
  if (OHREA.eqv..TRUE.) then
    if (verbosity.gt.0) then
313
      print*,'call readOHfield'
314
    endif
315
    call readOHfield
316
  endif
Matthias Langer's avatar
 
Matthias Langer committed
317
318
319
320
321

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

322
323
324
325
326
327
328
329
330
331
332
333
334
335
  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

336
  if (verbosity.gt.0) then
337
    print*,'call writeheader'
338
339
  endif

Matthias Langer's avatar
 
Matthias Langer committed
340
  call writeheader
341
  ! FLEXPART 9.2 ticket ?? write header in ASCII format 
342
343
344
345
346
347
348
349
350
  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
351
    print*,'call openreceptors'
352
  endif
Matthias Langer's avatar
 
Matthias Langer committed
353
354
355
356
357
358
  call openreceptors
  if ((iout.eq.4).or.(iout.eq.5)) call openouttraj

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

359
  if (verbosity.gt.0) then
360
    print*,'discretize release times'
361
  endif
Matthias Langer's avatar
 
Matthias Langer committed
362
  do i=1,numpoint
363
364
    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
365
366
367
368
369
  end do

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

370
  if (verbosity.gt.0) then
371
    print*,'Initialize cloud-base mass fluxes for the convection scheme'
372
373
  endif

Matthias Langer's avatar
 
Matthias Langer committed
374
375
376
377
378
379
380
381
382
383
384
385
386
  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

387
388
389
390
391
392
393
394
395
396
397
  ! 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
398
399
400
  ! Calculate particle trajectories
  !********************************

401
  if (verbosity.gt.0) then
402
403
404
405
406
407
408
409
410
     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'
411
412
  endif

Matthias Langer's avatar
 
Matthias Langer committed
413
414
  call timemanager

415
! NIK 16.02.2005 
416
417
418
419
420
421
422
423
424
425
  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
  
426
427
  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
       &XPART MODEL RUN!'
Matthias Langer's avatar
 
Matthias Langer committed
428
429

end program flexpart