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

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since fe32dca was fe32dca, checked in by Espen Sollum ATMOS <eso@…>, 6 years ago

Renamed lnokernel, corrected default setting.

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