source: branches/petra/src/FLEXPART.f90 @ 36

Last change on this file since 36 was 36, checked in by pesei, 7 years ago

Implement switch for incremental deposition, see ticket:113 and many small changes, see changelog.txt

File size: 12.6 KB
Line 
1!**********************************************************************
2! Copyright 1998-2015                                                 *
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
47#ifdef NETCDF_OUTPUT
48  use netcdf_output_mod, only: writeheader_netcdf
49#endif
50
51
52  implicit none
53
54  integer :: i,j,ix,jy,inest
55  integer :: idummy = -320
56  character(len=256) :: inline_options  !pathfile, flexversion, arg2
57
58
59  ! Generate a large number of random numbers
60  !******************************************
61
62  do i=1,maxrand-1,2
63    call gasdev1(idummy,rannumb(i),rannumb(i+1))
64  end do
65  call gasdev1(idummy,rannumb(maxrand),rannumb(maxrand-1))
66
67  ! FLEXPART version string
68  flexversion='Version 9.2_beta.r34 (2015-02-12)'
69  verbosity=0
70
71  ! Read the pathnames where input/output files are stored
72  !*******************************************************
73
74  inline_options='none'
75  select case (iargc())
76  case (2)
77    call getarg(1,arg1)
78    pathfile=arg1
79    call getarg(2,arg2)
80    inline_options=arg2
81  case (1)
82    call getarg(1,arg1)
83    pathfile=arg1
84    if (arg1(1:1).eq.'-') then
85      write(pathfile,'(a11)') './pathnames'
86      inline_options=arg1
87    endif
88  case (0)
89    write(pathfile,'(a11)') './pathnames'
90  end select
91 
92  ! Print the GPL License statement
93  !*******************************************************
94  print*,'Welcome to FLEXPART ', trim(flexversion)
95  print*,'FLEXPART is free software released under the GNU General Public License.'
96 
97  if (inline_options(1:1).eq.'-') then
98    if (trim(inline_options).eq.'-v'.or.trim(inline_options).eq.'-v1') then
99       print*, 'Verbose mode 1: display detailed information during run'
100       verbosity=1
101    endif
102    if (trim(inline_options).eq.'-v2') then
103       print*, 'Verbose mode 2: display more detailed information during run'
104       verbosity=2
105    endif
106    if (trim(inline_options).eq.'-i') then
107       print*, 'Info mode: provide detailed run specific information and stop'
108       verbosity=1
109       info_flag=1
110    endif
111    if (trim(inline_options).eq.'-i2') then
112       print*, 'Info mode: provide more detailed run specific information and stop'
113       verbosity=2
114       info_flag=1
115    endif
116  endif
117           
118  if (verbosity.gt.0) then
119    write(*,*) 'call readpaths'
120  endif
121  call readpaths
122 
123  if (verbosity.gt.1) then !show clock info
124     !print*,'length(4)',length(4)
125     !count=0,count_rate=1000
126     CALL SYSTEM_CLOCK(count_clock0, count_rate, count_max)
127     !WRITE(*,*) 'SYSTEM_CLOCK',count, count_rate, count_max
128     !WRITE(*,*) 'SYSTEM_CLOCK, count_clock0', count_clock0
129     !WRITE(*,*) 'SYSTEM_CLOCK, count_rate', count_rate
130     !WRITE(*,*) 'SYSTEM_CLOCK, count_max', count_max
131  endif
132
133  ! Read the user specifications for the current model run
134  !*******************************************************
135
136  if (verbosity.gt.0) then
137    write(*,*) 'call readcommand'
138  endif
139  call readcommand
140  if (verbosity.gt.0) then
141    write(*,*) '    ldirect=', ldirect
142    write(*,*) '    ibdate,ibtime=',ibdate,ibtime
143    write(*,*) '    iedate,ietime=', iedate,ietime
144    if (verbosity.gt.1) then   
145      CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
146      write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
147    endif     
148  endif
149
150  ! Read the age classes to be used
151  !********************************
152  if (verbosity.gt.0) then
153    write(*,*) 'call readageclasses'
154  endif
155  call readageclasses
156
157  if (verbosity.gt.1) then   
158    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
159    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
160  endif     
161
162  ! Read, which wind fields are available within the modelling period
163  !******************************************************************
164
165  if (verbosity.gt.0) then
166    write(*,*) 'call readavailable'
167  endif 
168  call readavailable
169
170  ! Read the model grid specifications,
171  ! both for the mother domain and eventual nests
172  !**********************************************
173 
174  if (verbosity.gt.0) then
175     write(*,*) 'call gridcheck'
176  endif
177
178  call gridcheck
179
180  if (verbosity.gt.1) then   
181    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
182    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
183  endif     
184
185  if (verbosity.gt.0) then
186    write(*,*) 'call gridcheck_nests'
187  endif 
188  call gridcheck_nests
189
190  ! Read the output grid specifications
191  !************************************
192
193  if (verbosity.gt.0) then
194    write(*,*) 'call readoutgrid'
195  endif
196
197  call readoutgrid
198
199  if (nested_output.eq.1) then
200    call readoutgrid_nest
201    if (verbosity.gt.0) then
202      write(*,*) '# readoutgrid_nest'
203    endif
204  endif
205
206  ! Read the receptor points for which extra concentrations are to be calculated
207  !*****************************************************************************
208
209  if (verbosity.eq.1) then
210     print*,'call readreceptors'
211  endif
212  call readreceptors
213
214  ! Read the physico-chemical species property table
215  !*************************************************
216  !SEC: now only needed SPECIES are read in readreleases.f
217  !call readspecies
218
219
220  ! Read the landuse inventory
221  !***************************
222
223  if (verbosity.gt.0) then
224    print*,'call readlanduse'
225  endif
226  call readlanduse
227
228  ! Assign fractional cover of landuse classes to each ECMWF grid point
229  !********************************************************************
230
231  if (verbosity.gt.0) then
232    print*,'call assignland'
233  endif
234  call assignland
235
236  ! Read the coordinates of the release locations
237  !**********************************************
238
239  if (verbosity.gt.0) then
240    print*,'call readreleases'
241  endif
242  call readreleases
243
244  ! Read and compute surface resistances to dry deposition of gases
245  !****************************************************************
246
247  if (verbosity.gt.0) then
248    print*,'call readdepo'
249  endif
250  call readdepo
251
252  ! Convert the release point coordinates from geografical to grid coordinates
253  !***************************************************************************
254
255  call coordtrafo 
256  if (verbosity.gt.0) then
257    print*,'call coordtrafo'
258  endif
259
260  ! Initialize all particles to non-existent
261  !*****************************************
262
263  if (verbosity.gt.0) then
264    print*,'Initialize all particles to non-existent'
265  endif
266  do j=1,maxpart
267    itra1(j)=-999999999
268  end do
269
270  ! For continuation of previous run, read in particle positions
271  !*************************************************************
272
273  if (ipin.eq.1) then
274    if (verbosity.gt.0) then
275      print*,'call readpartpositions'
276    endif
277    call readpartpositions
278  else
279    if (verbosity.gt.0) then
280      print*,'numpart=0, numparticlecount=0'
281    endif   
282    numpart=0
283    numparticlecount=0
284  endif
285
286  ! Calculate volume, surface area, etc., of all output grid cells
287  ! Allocate fluxes and OHfield if necessary
288  !***************************************************************
289
290  if (verbosity.gt.0) then
291    print*,'call outgrid_init'
292  endif
293  call outgrid_init
294  if (nested_output.eq.1) call outgrid_init_nest
295
296  ! Read the OH field
297  !******************
298
299  if (OHREA.eqv..TRUE.) then
300    if (verbosity.gt.0) then
301      print*,'call readOHfield'
302    endif
303    call readOHfield
304  endif
305
306  ! Write basic information on the simulation to a file "header"
307  ! and open files that are to be kept open throughout the simulation
308  !******************************************************************
309
310  if (lnetcdfout.eq.1) then
311#ifdef NETCDF_OUTPUT
312    call writeheader_netcdf(lnest = .false.)
313#endif
314  else
315    call writeheader
316  end if
317
318  if (nested_output.eq.1) then
319    if (lnetcdfout.eq.1) then
320#ifdef NETCDF_OUTPUT
321      call writeheader_netcdf(lnest = .true.)
322#endif
323    else
324      call writeheader_nest
325    endif
326  endif
327
328  if (verbosity.gt.0) then
329    print*,'call writeheader'
330  endif
331
332  call writeheader
333  ! FLEXPART 9.2 ticket ?? write header in ASCII format
334  call writeheader_txt
335  !if (nested_output.eq.1) call writeheader_nest
336  if (nested_output.eq.1.and.surf_only.ne.1) call writeheader_nest
337  if (nested_output.eq.1.and.surf_only.eq.1) call writeheader_nest_surf
338  if (nested_output.ne.1.and.surf_only.eq.1) call writeheader_surf
339
340  if (lnetcdfout.ne.1) then
341    open(unitdates,file=path(2)(1:length(2))//'dates')
342  end if
343
344  if (verbosity.gt.0) then
345    print*,'call openreceptors'
346  endif
347  call openreceptors
348  if ((iout.eq.4).or.(iout.eq.5)) call openouttraj
349
350  ! Releases can only start and end at discrete times (multiples of lsynctime)
351  !***************************************************************************
352
353  if (verbosity.gt.0) then
354    print*,'discretize release times'
355  endif
356  do i=1,numpoint
357    ireleasestart(i)=nint(real(ireleasestart(i))/real(lsynctime))*lsynctime
358    ireleaseend(i)=nint(real(ireleaseend(i))/real(lsynctime))*lsynctime
359  end do
360
361  ! Initialize cloud-base mass fluxes for the convection scheme
362  !************************************************************
363
364  if (verbosity.gt.0) then
365    print*,'Initialize cloud-base mass fluxes for the convection scheme'
366  endif
367
368  do jy=0,nymin1
369    do ix=0,nxmin1
370      cbaseflux(ix,jy)=0.
371    end do
372  end do
373  do inest=1,numbnests
374    do jy=0,nyn(inest)-1
375      do ix=0,nxn(inest)-1
376        cbasefluxn(ix,jy,inest)=0.
377      end do
378    end do
379  end do
380
381  ! Calculate particle trajectories
382  !********************************
383
384  if (verbosity.gt.0) then
385     if (verbosity.gt.1) then   
386       call system_clock(count_clock, count_rate, count_max)
387       write(*,*) 'System clock',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
388     endif
389     if (info_flag.eq.1) then
390       print*, 'Info only mode (stop)'   
391       stop
392     endif
393     print*,'call timemanager'
394  endif
395
396  call timemanager
397
398  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLEXPART MODEL RUN!'
399
400  ! output wall time
401  if (verbosity .gt. 0) then
402    call system_clock(count_clock,count_rate)
403    tins=(count_clock - count_clock0)/real(count_rate)
404    print*,'Wall time ',tins,'s, ',tins/60,'min, ',tins/3600,'h.'
405  endif
406
407end program flexpart
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG