FLEXPART.f90 13.2 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
  ! FLEXPART version string
70
71
  flexversion_major = '10' ! Major version number, also used for species file names
  flexversion='Version '//trim(flexversion_major)//'.0beta (2015-05-01)'
72
73
  verbosity=0

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

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

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

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

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

160
  if (verbosity.gt.1) then   
161
162
    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
163
  endif     
Matthias Langer's avatar
 
Matthias Langer committed
164
165
166
167

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

317
318
319
320
321
  !! testing !!
  ! open(999,file=trim(path(1))//'OH_FIELDS/jscalar_50N.txt',action='write',status='new')
  ! open(998,file=trim(path(1))//'OH_FIELDS/jscalar_50S.txt',action='write',status='new')


Matthias Langer's avatar
 
Matthias Langer committed
322
323
324
325
  ! Write basic information on the simulation to a file "header"
  ! and open files that are to be kept open throughout the simulation
  !******************************************************************

326
327
328
329
330
331
332
333
334
335
336
337
338
339
  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

340
  if (verbosity.gt.0) then
341
    print*,'call writeheader'
342
343
  endif

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

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

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

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

374
  if (verbosity.gt.0) then
375
    print*,'Initialize cloud-base mass fluxes for the convection scheme'
376
377
  endif

Matthias Langer's avatar
 
Matthias Langer committed
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
  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
  !********************************

394
  if (verbosity.gt.0) then
395
396
397
398
399
400
401
402
403
     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'
404
405
  endif

Matthias Langer's avatar
 
Matthias Langer committed
406
407
  call timemanager

408
409
410
411
412
413
414
415
! NIK 16.02.2005 
  write(*,*) '**********************************************'
  write(*,*) 'Total number of occurences of below-cloud scavenging', tot_blc_count
  write(*,*) 'Total number of occurences of in-cloud    scavenging', tot_inc_count
  write(*,*) '**********************************************'

  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
       &XPART MODEL RUN!'
Matthias Langer's avatar
 
Matthias Langer committed
416
417

end program flexpart