FLEXPART.f90 11.9 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
46
47
48
49
50
!**********************************************************************
! 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

  implicit none

  integer :: i,j,ix,jy,inest
  integer :: idummy = -320
51
52
  character(len=256) :: inline_options  !pathfile, flexversion, arg2

Matthias Langer's avatar
 
Matthias Langer committed
53
54
55
56
57
58
59
60
61

  ! 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))

62
63
64
65
  ! FLEXPART version string
  flexversion='Version 9.2 beta (2014-07-01)'
  verbosity=0

66
67
68
69
70
71
72
73
74
75
76
77
78
79
  ! Read the pathnames where input/output files are stored
  !*******************************************************

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

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

131
  if (verbosity.gt.0) then
132
    write(*,*) 'call readcommand'
133
  endif
Matthias Langer's avatar
 
Matthias Langer committed
134
  call readcommand
135
  if (verbosity.gt.0) then
136
137
138
139
140
141
142
    write(*,*) '    ldirect=', ldirect
    write(*,*) '    ibdate,ibtime=',ibdate,ibtime
    write(*,*) '    iedate,ietime=', iedate,ietime
    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     
143
  endif
Matthias Langer's avatar
 
Matthias Langer committed
144
145
146

  ! Read the age classes to be used
  !********************************
147
  if (verbosity.gt.0) then
148
    write(*,*) 'call readageclasses'
149
  endif
Matthias Langer's avatar
 
Matthias Langer committed
150
151
  call readageclasses

152
  if (verbosity.gt.1) then   
153
154
    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
155
  endif     
Matthias Langer's avatar
 
Matthias Langer committed
156
157
158
159

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

160
  if (verbosity.gt.0) then
161
    write(*,*) 'call readavailable'
162
  endif  
Matthias Langer's avatar
 
Matthias Langer committed
163
164
165
166
167
  call readavailable

  ! Read the model grid specifications,
  ! both for the mother domain and eventual nests
  !**********************************************
168
169
 
  if (verbosity.gt.0) then
170
     write(*,*) 'call gridcheck'
171
  endif
Matthias Langer's avatar
 
Matthias Langer committed
172
173

  call gridcheck
174
175

  if (verbosity.gt.1) then   
176
177
    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
178
179
180
  endif      

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

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

188
  if (verbosity.gt.0) then
189
    write(*,*) 'call readoutgrid'
190
191
  endif

Matthias Langer's avatar
 
Matthias Langer committed
192
193
  call readoutgrid

194
  if (nested_output.eq.1) then
195
    call readoutgrid_nest
196
    if (verbosity.gt.0) then
197
      write(*,*) '# readoutgrid_nest'
198
199
    endif
  endif
Matthias Langer's avatar
 
Matthias Langer committed
200
201
202
203

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

204
205
206
  if (verbosity.eq.1) then
     print*,'call readreceptors'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
207
208
209
210
211
212
213
214
215
216
217
  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
  !***************************

218
219
220
  if (verbosity.gt.0) then
    print*,'call readlanduse'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
221
222
223
224
225
  call readlanduse

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

226
227
228
  if (verbosity.gt.0) then
    print*,'call assignland'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
229
230
231
232
233
  call assignland

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

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

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

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

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

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

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

258
259
260
  if (verbosity.gt.0) then
    print*,'Initialize all particles to non-existent'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
261
262
263
264
265
266
267
268
  do j=1,maxpart
    itra1(j)=-999999999
  end do

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

  if (ipin.eq.1) then
269
270
271
    if (verbosity.gt.0) then
      print*,'call readpartpositions'
    endif
Matthias Langer's avatar
 
Matthias Langer committed
272
273
    call readpartpositions
  else
274
275
276
    if (verbosity.gt.0) then
      print*,'numpart=0, numparticlecount=0'
    endif    
Matthias Langer's avatar
 
Matthias Langer committed
277
278
279
280
281
282
283
284
    numpart=0
    numparticlecount=0
  endif

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

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

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

294
295
  if (OHREA.eqv..TRUE.) then
    if (verbosity.gt.0) then
296
      print*,'call readOHfield'
297
    endif
298
    call readOHfield
299
  endif
Matthias Langer's avatar
 
Matthias Langer committed
300
301
302
303
304

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

305
306
307
308
  if (verbosity.gt.0) then
    print*,'call writeheader'
  endif

Matthias Langer's avatar
 
Matthias Langer committed
309
  call writeheader
310
311
312
313
314
315
316
317
318
319
320
321
  ! FLEXPART 9.2 ticket ?? write header in ASCII format 
  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
    print*,'call openreceptors'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
322
323
324
325
326
327
  call openreceptors
  if ((iout.eq.4).or.(iout.eq.5)) call openouttraj

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

328
329
330
  if (verbosity.gt.0) then
    print*,'discretize release times'
  endif
Matthias Langer's avatar
 
Matthias Langer committed
331
  do i=1,numpoint
332
333
    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
334
335
336
337
338
  end do

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

339
340
341
342
  if (verbosity.gt.0) then
    print*,'Initialize cloud-base mass fluxes for the convection scheme'
  endif

Matthias Langer's avatar
 
Matthias Langer committed
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
  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

  ! Calculate particle trajectories
  !********************************

359
360
361
  if (verbosity.gt.0) then
     if (verbosity.gt.1) then   
       CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
362
       write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
363
364
     endif
     if (info_flag.eq.1) then
365
366
       print*, 'info only mode (stop)'    
       stop
367
368
369
370
     endif
     print*,'call timemanager'
  endif

Matthias Langer's avatar
 
Matthias Langer committed
371
372
  call timemanager

373
  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLEXPART MODEL RUN!'
Matthias Langer's avatar
 
Matthias Langer committed
374
375

end program flexpart