source: flexpart.git/src/FLEXPART.f90

10.4.1_pesei
Last change on this file was 92fab65, checked in by Ignacio Pisso <ip@…>, 8 months ago

add SPDX-License-Identifier to all .f90 files

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