source: flexpart.git/src/FLEXPART.f90 @ c8fc724

devrelease-10univie
Last change on this file since c8fc724 was c8fc724, checked in by Espen Sollum ATMOS <eso@…>, 3 years ago

Bugfix: using namelist species files would sometimes set incorrect below-cloud scavenging parameters when using more than 1 species. Update: per-species scavenging statistics are written out at end of simulation.

  • Property mode set to 100644
File size: 13.5 KB
Line 
1!**********************************************************************
2! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
3! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
4! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
5!                                                                     *
6! This file is part of FLEXPART.                                      *
7!                                                                     *
8! FLEXPART is free software: you can redistribute it and/or modify    *
9! it under the terms of the GNU General Public License as published by*
10! the Free Software Foundation, either version 3 of the License, or   *
11! (at your option) any later version.                                 *
12!                                                                     *
13! FLEXPART is distributed in the hope that it will be useful,         *
14! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
15! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
16! GNU General Public License for more details.                        *
17!                                                                     *
18! You should have received a copy of the GNU General Public License   *
19! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
20!**********************************************************************
21
22program flexpart
23
24  !*****************************************************************************
25  !                                                                            *
26  !     This is the Lagrangian Particle Dispersion Model FLEXPART.             *
27  !     The main program manages the reading of model run specifications, etc. *
28  !     All actual computing is done within subroutine timemanager.            *
29  !                                                                            *
30  !     Author: A. Stohl                                                       *
31  !                                                                            *
32  !     18 May 1996                                                            *
33  !                                                                            *
34  !*****************************************************************************
35  !                                                                            *
36  ! Variables:                                                                 *
37  !                                                                            *
38  ! Constants:                                                                 *
39  !                                                                            *
40  !*****************************************************************************
41
42  use point_mod
43  use par_mod
44  use com_mod
45  use conv_mod
46  use netcdf_output_mod, only: writeheader_netcdf
47  use random_mod, only: gasdev1
48
49  implicit none
50
51  integer :: i,j,ix,jy,inest
52  integer :: idummy = -320
53  character(len=256) :: inline_options  !pathfile, flexversion, arg2
54
55
56  ! Initialize arrays in com_mod
57  !*****************************
58  call com_mod_allocate_part(maxpart)
59 
60
61  ! Generate a large number of random numbers
62  !******************************************
63
64  do i=1,maxrand-1,2
65    call gasdev1(idummy,rannumb(i),rannumb(i+1))
66  end do
67  call gasdev1(idummy,rannumb(maxrand),rannumb(maxrand-1))
68
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  verbosity=0
74
75  ! Read the pathnames where input/output files are stored
76  !*******************************************************
77
78  inline_options='none'
79  select case (iargc())
80  case (2)
81    call getarg(1,arg1)
82    pathfile=arg1
83    call getarg(2,arg2)
84    inline_options=arg2
85  case (1)
86    call getarg(1,arg1)
87    pathfile=arg1
88    if (arg1(1:1).eq.'-') then
89      write(pathfile,'(a11)') './pathnames'
90      inline_options=arg1
91    endif
92  case (0)
93    write(pathfile,'(a11)') './pathnames'
94  end select
95 
96  ! Print the GPL License statement
97  !*******************************************************
98  print*,'Welcome to FLEXPART ', trim(flexversion)
99  print*,'FLEXPART is free software released under the GNU General Public License.'
100 
101  if (inline_options(1:1).eq.'-') then
102    if (trim(inline_options).eq.'-v'.or.trim(inline_options).eq.'-v1') then
103       print*, 'Verbose mode 1: display detailed information during run'
104       verbosity=1
105    endif
106    if (trim(inline_options).eq.'-v2') then
107       print*, 'Verbose mode 2: display more detailed information during run'
108       verbosity=2
109    endif
110    if (trim(inline_options).eq.'-i') then
111       print*, 'Info mode: provide detailed run specific information and stop'
112       verbosity=1
113       info_flag=1
114    endif
115    if (trim(inline_options).eq.'-i2') then
116       print*, 'Info mode: provide more detailed run specific information and stop'
117       verbosity=2
118       info_flag=1
119    endif
120  endif
121           
122  if (verbosity.gt.0) then
123    write(*,*) 'call readpaths'
124  endif
125  call readpaths(pathfile)
126 
127  if (verbosity.gt.1) then !show clock info
128     !print*,'length(4)',length(4)
129     !count=0,count_rate=1000
130     CALL SYSTEM_CLOCK(count_clock0, count_rate, count_max)
131     !WRITE(*,*) 'SYSTEM_CLOCK',count, count_rate, count_max
132     !WRITE(*,*) 'SYSTEM_CLOCK, count_clock0', count_clock0
133     !WRITE(*,*) 'SYSTEM_CLOCK, count_rate', count_rate
134     !WRITE(*,*) 'SYSTEM_CLOCK, count_max', count_max
135  endif
136
137  ! Read the user specifications for the current model run
138  !*******************************************************
139
140  if (verbosity.gt.0) then
141    write(*,*) 'call readcommand'
142  endif
143  call readcommand
144  if (verbosity.gt.0) then
145    write(*,*) '    ldirect=', ldirect
146    write(*,*) '    ibdate,ibtime=',ibdate,ibtime
147    write(*,*) '    iedate,ietime=', iedate,ietime
148    if (verbosity.gt.1) then   
149      CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
150      write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
151    endif     
152  endif
153
154  ! Read the age classes to be used
155  !********************************
156  if (verbosity.gt.0) then
157    write(*,*) 'call readageclasses'
158  endif
159  call readageclasses
160
161  if (verbosity.gt.1) then   
162    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
163    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
164  endif     
165
166  ! Read, which wind fields are available within the modelling period
167  !******************************************************************
168
169  if (verbosity.gt.0) then
170    write(*,*) 'call readavailable'
171  endif 
172  call readavailable
173
174  ! If nested wind fields are used, allocate arrays
175  !************************************************
176
177  if (verbosity.gt.0) then
178    write(*,*) 'call com_mod_allocate_nests'
179  endif
180  call com_mod_allocate_nests
181
182  ! Read the model grid specifications,
183  ! both for the mother domain and eventual nests
184  !**********************************************
185 
186  if (verbosity.gt.0) then
187     write(*,*) 'call gridcheck'
188  endif
189
190  call gridcheck
191
192  if (verbosity.gt.1) then   
193    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
194    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
195  endif     
196
197  if (verbosity.gt.0) then
198    write(*,*) 'call gridcheck_nests'
199  endif 
200  call gridcheck_nests
201
202  ! Read the output grid specifications
203  !************************************
204
205  if (verbosity.gt.0) then
206    write(*,*) 'call readoutgrid'
207  endif
208
209  call readoutgrid
210
211  if (nested_output.eq.1) then
212    call readoutgrid_nest
213    if (verbosity.gt.0) then
214      write(*,*) '# readoutgrid_nest'
215    endif
216  endif
217
218  ! Read the receptor points for which extra concentrations are to be calculated
219  !*****************************************************************************
220
221  if (verbosity.eq.1) then
222     print*,'call readreceptors'
223  endif
224  call readreceptors
225
226  ! Read the physico-chemical species property table
227  !*************************************************
228  !SEC: now only needed SPECIES are read in readreleases.f
229  !call readspecies
230
231
232  ! Read the landuse inventory
233  !***************************
234
235  if (verbosity.gt.0) then
236    print*,'call readlanduse'
237  endif
238  call readlanduse
239
240  ! Assign fractional cover of landuse classes to each ECMWF grid point
241  !********************************************************************
242
243  if (verbosity.gt.0) then
244    print*,'call assignland'
245  endif
246  call assignland
247
248  ! Read the coordinates of the release locations
249  !**********************************************
250
251  if (verbosity.gt.0) then
252    print*,'call readreleases'
253  endif
254  call readreleases
255
256  ! Read and compute surface resistances to dry deposition of gases
257  !****************************************************************
258
259  if (verbosity.gt.0) then
260    print*,'call readdepo'
261  endif
262  call readdepo
263
264  ! Convert the release point coordinates from geografical to grid coordinates
265  !***************************************************************************
266
267  call coordtrafo 
268  if (verbosity.gt.0) then
269    print*,'call coordtrafo'
270  endif
271
272  ! Initialize all particles to non-existent
273  !*****************************************
274
275  if (verbosity.gt.0) then
276    print*,'Initialize all particles to non-existent'
277  endif
278  do j=1,maxpart
279    itra1(j)=-999999999
280  end do
281
282  ! For continuation of previous run, read in particle positions
283  !*************************************************************
284
285  if (ipin.eq.1) then
286    if (verbosity.gt.0) then
287      print*,'call readpartpositions'
288    endif
289    call readpartpositions
290  else
291    if (verbosity.gt.0) then
292      print*,'numpart=0, numparticlecount=0'
293    endif   
294    numpart=0
295    numparticlecount=0
296  endif
297
298  ! Calculate volume, surface area, etc., of all output grid cells
299  ! Allocate fluxes and OHfield if necessary
300  !***************************************************************
301
302  if (verbosity.gt.0) then
303    print*,'call outgrid_init'
304  endif
305  call outgrid_init
306  if (nested_output.eq.1) call outgrid_init_nest
307
308  ! Read the OH field
309  !******************
310
311  if (OHREA.eqv..TRUE.) then
312    if (verbosity.gt.0) then
313      print*,'call readOHfield'
314    endif
315    call readOHfield
316  endif
317
318  ! Write basic information on the simulation to a file "header"
319  ! and open files that are to be kept open throughout the simulation
320  !******************************************************************
321
322  if (lnetcdfout.eq.1) then
323    call writeheader_netcdf(lnest=.false.)
324  else
325    call writeheader
326  end if
327
328  if (nested_output.eq.1) then
329    if (lnetcdfout.eq.1) then
330      call writeheader_netcdf(lnest=.true.)
331    else
332      call writeheader_nest
333    endif
334  endif
335
336  if (verbosity.gt.0) then
337    print*,'call writeheader'
338  endif
339
340  call writeheader
341  ! FLEXPART 9.2 ticket ?? write header in ASCII format
342  call writeheader_txt
343  !if (nested_output.eq.1) call writeheader_nest
344  if (nested_output.eq.1.and.surf_only.ne.1) call writeheader_nest
345  if (nested_output.eq.1.and.surf_only.eq.1) call writeheader_nest_surf
346  if (nested_output.ne.1.and.surf_only.eq.1) call writeheader_surf
347
348  !open(unitdates,file=path(2)(1:length(2))//'dates')
349
350  if (verbosity.gt.0) then
351    print*,'call openreceptors'
352  endif
353  call openreceptors
354  if ((iout.eq.4).or.(iout.eq.5)) call openouttraj
355
356  ! Releases can only start and end at discrete times (multiples of lsynctime)
357  !***************************************************************************
358
359  if (verbosity.gt.0) then
360    print*,'discretize release times'
361  endif
362  do i=1,numpoint
363    ireleasestart(i)=nint(real(ireleasestart(i))/real(lsynctime))*lsynctime
364    ireleaseend(i)=nint(real(ireleaseend(i))/real(lsynctime))*lsynctime
365  end do
366
367  ! Initialize cloud-base mass fluxes for the convection scheme
368  !************************************************************
369
370  if (verbosity.gt.0) then
371    print*,'Initialize cloud-base mass fluxes for the convection scheme'
372  endif
373
374  do jy=0,nymin1
375    do ix=0,nxmin1
376      cbaseflux(ix,jy)=0.
377    end do
378  end do
379  do inest=1,numbnests
380    do jy=0,nyn(inest)-1
381      do ix=0,nxn(inest)-1
382        cbasefluxn(ix,jy,inest)=0.
383      end do
384    end do
385  end do
386
387  ! Inform whether output kernel is used or not
388  !*********************************************
389  if (lroot) then
390    if (lnokernel) then
391      write(*,*) "Concentrations are calculated without using kernel"
392    else
393      write(*,*) "Concentrations are calculated using kernel"
394    end if
395  end if
396
397
398  ! Calculate particle trajectories
399  !********************************
400
401  if (verbosity.gt.0) then
402     if (verbosity.gt.1) then   
403       CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
404       write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
405     endif
406     if (info_flag.eq.1) then
407       print*, 'info only mode (stop)'   
408       stop
409     endif
410     print*,'call timemanager'
411  endif
412
413  call timemanager
414
415! NIK 16.02.2005
416  do i=1,nspec
417    write(*,*) '**********************************************'
418    write(*,*) 'Scavenging statistics for species ', species(i), ':'
419    write(*,*) 'Total number of occurences of below-cloud scavenging', &
420         & tot_blc_count(i)
421    write(*,*) 'Total number of occurences of in-cloud    scavenging', &
422         & tot_inc_count(i)
423    write(*,*) '**********************************************'
424  end do
425 
426  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
427       &XPART MODEL RUN!'
428
429end program flexpart
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG