source: flexpart.git/src/FLEXPART_MPI.f90

bugfixes+enhancements
Last change on this file was dba4221, checked in by Pirmin Kaufmann <pirmin.kaufmann@…>, 18 months ago

Bugfixes:

  • options/SPECIES/SPECIES_009: corrected wrong number format, replaced comma by decimal point
  • options/SPECIES/SPECIES_028: corrected wrong number format, moved sign of exponent to after the E
  • options/SPECIES/specoverview.f90: added namelist parameters that appear in SPECIES files but were missing here
  • src/FLEXPART.f90: replaced compiler-specific command line argument routines by standard Fortran intrinsic routines
  • src/FLEXPART_MPI.f90: ditto
  • src/gridcheck_ecmwf.f90: corrected handling of vertical levels when input files do not contain uppermost layers
  • src/gridcheck_nests.f90: ditto
  • src/readwind_ecmwf.f90: corrected handling of vertical levels when input files do not contain uppermost layers
  • readwind_ecmwf_mpi.f90: ditto

Code enhancements:

  • options/OUTGRID: added comments describing contents
  • options/SPECIES/SPECIES_*: aligned comments
  • options/SPECIES/specoverview.f90: removed commented lines, rectified lines indenting
  • src/FLEXPART.f90: rectified lines indenting, updated date in version string
  • src/FLEXPART_MPI.f90: ditto, and realigned code with src/FLEXPART.f90
  • src/gridcheck_*.f90: added code to write out name of file before it is opened (helps a lot when an input file causes troubles)
  • src/par_mod.f90: added comment explaining relevance of nuvzmax for GRIB input
  • src/readreleases.f90: write out warning if too few particles are used to randomize release
  • src/readspecies.f90: write out name of SPECIES file before it is read
  • src/readwind_*.f90: write out name of input file before opening it
  • src/writeheader_txt.f90: removed wrong comment
  • Property mode set to 100644
File size: 15.5 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  use mpi_mod
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  integer(selected_int_kind(16)), dimension(maxspec) :: tot_b=0, &
52       & tot_i=0
53
54
55  ! Initialize mpi
56  !*********************
57  call mpif_init
58
59  if (mp_measure_time) call mpif_mtime('flexpart',0)
60
61 
62  ! Generate a large number of random numbers
63  !******************************************
64
65  ! eso: Different seed for each MPI process
66  idummy=idummy+mp_seed
67
68  do i=1,maxrand-1,2
69    call gasdev1(idummy,rannumb(i),rannumb(i+1))
70  end do
71  call gasdev1(idummy,rannumb(maxrand),rannumb(maxrand-1))
72
73
74  ! FLEXPART version string
75  flexversion_major = '10' ! Major version number, also used for species file names
76  flexversion='Version '//trim(flexversion_major)//'.4.1 MPI (2022-11-09)'
77  verbosity=0
78
79  ! Read the pathnames where input/output files are stored
80  !*******************************************************
81
82  inline_options='none'
83  select case (command_argument_count())  ! Portable standard Fortran intrinsic procedure
84  case (2)
85    call get_command_argument(1,arg1)     ! Portable standard Fortran intrinsic procedure
86    pathfile=arg1
87    call get_command_argument(2,arg2)     ! Portable standard Fortran intrinsic procedure
88    inline_options=arg2
89  case (1)
90    call get_command_argument(1,arg1)     ! Portable standard Fortran intrinsic procedure
91    pathfile=arg1
92    if (arg1(1:1).eq.'-') then
93      write(pathfile,'(a11)') './pathnames'
94      inline_options=arg1
95    endif
96  case (0)
97    write(pathfile,'(a11)') './pathnames'
98  end select
99 
100  ! Print the GPL License statement
101  !*******************************************************
102  if (lroot) then
103    print*,'Welcome to FLEXPART ', trim(flexversion)
104    print*,'FLEXPART is free software released under the GNU General Public License.'
105  endif
106
107
108  ! Ingest inline options
109  !*******************************************************
110  if (inline_options(1:1).eq.'-') then
111    print*,'inline_options:',inline_options
112    !verbose mode
113    iopt=index(inline_options,'v')
114    if (iopt.gt.0) then
115      verbosity=1
116      !print*, iopt, inline_options(iopt+1:iopt+1)
117      if  (trim(inline_options(iopt+1:iopt+1)).eq.'2') then
118        print*, 'Verbose mode 2: display more detailed information during run'
119        verbosity=2
120      endif
121    endif
122    !debug mode
123    iopt=index(inline_options,'d')
124    if (iopt.gt.0) then
125      debug_mode=.true.
126    endif
127    if (trim(inline_options).eq.'-i') then
128      print*, 'Info mode: provide detailed run specific information and stop'
129      verbosity=1
130      info_flag=1
131    endif
132    if (trim(inline_options).eq.'-i2') then
133      print*, 'Info mode: provide more detailed run specific information and stop'
134      verbosity=2
135      info_flag=1
136    endif
137  endif
138           
139  if (verbosity.gt.0) then
140    print*, 'nxmax=',nxmax
141    print*, 'nymax=',nymax
142    print*, 'nzmax=',nzmax
143    print*,'nxshift=',nxshift
144  endif
145 
146  if (verbosity.gt.0) then
147    write(*,*) 'call readpaths'
148  endif
149  call readpaths
150 
151  if (verbosity.gt.1) then !show clock info
152     !print*,'length(4)',length(4)
153     !count=0,count_rate=1000
154     call system_clock(count_clock0, count_rate, count_max)
155     !WRITE(*,*) 'SYSTEM_CLOCK',count, count_rate, count_max
156     !WRITE(*,*) 'SYSTEM_CLOCK, count_clock0', count_clock0
157     !WRITE(*,*) 'SYSTEM_CLOCK, count_rate', count_rate
158     !WRITE(*,*) 'SYSTEM_CLOCK, count_max', count_max
159  endif
160
161  ! Read the user specifications for the current model run
162  !*******************************************************
163
164  if (verbosity.gt.0) then
165    write(*,*) 'call readcommand'
166  endif
167  call readcommand
168  if (verbosity.gt.0 .and. lroot) then
169    write(*,*) '    ldirect=', ldirect
170    write(*,*) '    ibdate,ibtime=',ibdate,ibtime
171    write(*,*) '    iedate,ietime=', iedate,ietime
172    if (verbosity.gt.1) then   
173      CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
174      write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
175    endif
176  endif
177
178  ! Initialize arrays in com_mod
179  !*****************************
180  if(.not.(lmpreader.and.lmp_use_reader)) call com_mod_allocate_part(maxpart_mpi)
181
182
183  ! Read the age classes to be used
184  !********************************
185  if (verbosity.gt.0 .and. lroot) then
186    write(*,*) 'call readageclasses'
187  endif
188  call readageclasses
189
190  if (verbosity.gt.1 .and. lroot) then   
191    call system_clock(count_clock, count_rate, count_max)
192    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
193  endif     
194
195  ! Read, which wind fields are available within the modelling period
196  !******************************************************************
197
198  if (verbosity.gt.0 .and. lroot) then
199    write(*,*) 'call readavailable'
200  endif 
201  call readavailable
202
203  ! Detect metdata format
204  !**********************
205  if (verbosity.gt.0 .and. lroot) then
206    write(*,*) 'call detectformat'
207  endif
208
209  metdata_format = detectformat()
210
211  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
212    if (lroot) print *,'ECMWF metdata detected'
213  elseif (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
214    if (lroot) print *,'NCEP metdata detected'
215  else
216    if (lroot) print *,'Unknown metdata format'
217    stop
218  endif
219
220
221
222  ! If nested wind fields are used, allocate arrays
223  !************************************************
224
225  if (verbosity.gt.0 .and. lroot) then
226    write(*,*) 'call com_mod_allocate_nests'
227  endif
228  call com_mod_allocate_nests
229
230  ! Read the model grid specifications,
231  ! both for the mother domain and eventual nests
232  !**********************************************
233
234  if (verbosity.gt.0 .and. lroot) then
235    write(*,*) 'call gridcheck'
236  endif
237
238  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
239    call gridcheck_ecmwf
240  else
241    call gridcheck_gfs
242  end if
243
244  if (verbosity.gt.1 .and. lroot) then   
245    call system_clock(count_clock, count_rate, count_max)
246    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
247  endif     
248
249  if (verbosity.gt.0 .and. lroot) then
250    write(*,*) 'call gridcheck_nests'
251  endif 
252  call gridcheck_nests
253
254  ! Read the output grid specifications
255  !************************************
256
257  if (verbosity.gt.0 .and. lroot) then
258    write(*,*) 'call readoutgrid'
259  endif
260
261  call readoutgrid
262
263  if (nested_output.eq.1) then
264    call readoutgrid_nest
265    if (verbosity.gt.0.and.lroot) then
266      write(*,*) '# readoutgrid_nest'
267    endif
268  endif
269
270  ! Read the receptor points for which extra concentrations are to be calculated
271  !*****************************************************************************
272
273  if (verbosity.eq.1 .and. lroot) then
274    print*,'call readreceptors'
275  endif
276  call readreceptors
277
278  ! Read the physico-chemical species property table
279  !*************************************************
280  !SEC: now only needed SPECIES are read in readreleases.f
281  !call readspecies
282
283  ! Read the landuse inventory
284  !***************************
285
286  if (verbosity.gt.0 .and. lroot) then
287    print*,'call readlanduse'
288  endif
289  call readlanduse
290
291  ! Assign fractional cover of landuse classes to each ECMWF grid point
292  !********************************************************************
293
294  if (verbosity.gt.0 .and. lroot) then
295    print*,'call assignland'
296  endif
297  call assignland
298
299  ! Read the coordinates of the release locations
300  !**********************************************
301
302  if (verbosity.gt.0 .and. lroot) then
303    print*,'call readreleases'
304  endif
305  call readreleases
306
307  ! Read and compute surface resistances to dry deposition of gases
308  !****************************************************************
309
310  if (verbosity.gt.0 .and. lroot) then
311    print*,'call readdepo'
312  endif
313  call readdepo
314
315  ! Convert the release point coordinates from geografical to grid coordinates
316  !***************************************************************************
317
318  call coordtrafo 
319  if (verbosity.gt.0 .and. lroot) then
320    print*,'call coordtrafo'
321  endif
322
323  ! Initialize all particles to non-existent
324  !*****************************************
325
326  if (verbosity.gt.0 .and. lroot) then
327    print*,'Initialize all particles to non-existent'
328  endif
329
330  if (.not.(lmpreader.and.lmp_use_reader)) then
331    do j=1, size(itra1) ! maxpart_mpi
332      itra1(j)=-999999999
333    end do
334  end if
335
336  ! For continuation of previous run, read in particle positions
337  !*************************************************************
338
339  if (ipin.eq.1) then
340    if (verbosity.gt.0 .and. lroot) then
341      print*,'call readpartpositions'
342    endif
343    ! readwind process skips this step
344    if (.not.(lmpreader.and.lmp_use_reader)) call readpartpositions
345  else
346    if (verbosity.gt.0 .and. lroot) then
347      print*,'numpart=0, numparticlecount=0'
348    endif
349    numpart=0
350    numparticlecount=0
351  endif
352
353  ! Calculate volume, surface area, etc., of all output grid cells
354  ! Allocate fluxes and OHfield if necessary
355  !***************************************************************
356
357  if (verbosity.gt.0 .and. lroot) then
358    print*,'call outgrid_init'
359  endif
360  call outgrid_init
361  if (nested_output.eq.1) call outgrid_init_nest
362
363  ! Read the OH field
364  !******************
365
366  if (OHREA.eqv..true.) then
367    if (verbosity.gt.0 .and. lroot) then
368      print*,'call readOHfield'
369    endif
370    call readOHfield
371  endif
372
373  ! Write basic information on the simulation to a file "header"
374  ! and open files that are to be kept open throughout the simulation
375  !******************************************************************
376
377  if (mp_measure_time) call mpif_mtime('iotime',0)
378
379  if (lroot) then ! MPI: this part root process only
380#ifdef USE_NCF
381    if (lnetcdfout.eq.1) then
382      call writeheader_netcdf(lnest=.false.)
383    else
384      call writeheader
385    end if
386   
387    if (nested_output.eq.1) then
388      if (lnetcdfout.eq.1) then
389        call writeheader_netcdf(lnest=.true.)
390      else
391        call writeheader_nest
392      endif
393    endif
394#endif
395
396    if (verbosity.gt.0) then
397      print*,'call writeheader'
398    endif
399
400    call writeheader
401    ! FLEXPART 9.2 ticket ?? write header in ASCII format
402    call writeheader_txt
403
404    if (nested_output.eq.1.and.surf_only.ne.1) call writeheader_nest
405    if (nested_output.eq.1.and.surf_only.eq.1) call writeheader_nest_surf
406    if (nested_output.ne.1.and.surf_only.eq.1) call writeheader_surf
407  end if ! (mpif_pid == 0)
408
409  if (mp_measure_time) call mpif_mtime('iotime',1)
410
411  if (verbosity.gt.0 .and. lroot) then
412    print*,'call openreceptors'
413  endif
414  call openreceptors
415  if ((iout.eq.4).or.(iout.eq.5)) call openouttraj
416
417  ! Releases can only start and end at discrete times (multiples of lsynctime)
418  !***************************************************************************
419
420  if (verbosity.gt.0 .and. lroot) then
421    print*,'discretize release times'
422  endif
423  do i=1,numpoint
424    ireleasestart(i)=nint(real(ireleasestart(i))/real(lsynctime))*lsynctime
425    ireleaseend(i)=nint(real(ireleaseend(i))/real(lsynctime))*lsynctime
426  end do
427
428  ! Initialize cloud-base mass fluxes for the convection scheme
429  !************************************************************
430
431  if (verbosity.gt.0 .and. lroot) then
432    print*,'Initialize cloud-base mass fluxes for the convection scheme'
433  endif
434
435  do jy=0,nymin1
436    do ix=0,nxmin1
437      cbaseflux(ix,jy)=0.
438    end do
439  end do
440  do inest=1,numbnests
441    do jy=0,nyn(inest)-1
442      do ix=0,nxn(inest)-1
443        cbasefluxn(ix,jy,inest)=0.
444      end do
445    end do
446  end do
447
448  ! Inform whether output kernel is used or not
449  !*********************************************
450
451  if (lroot) then
452    if (.not.lusekerneloutput) then
453      write(*,*) "Concentrations are calculated without using kernel"
454    else
455      write(*,*) "Concentrations are calculated using kernel"
456    end if
457  end if
458
459  ! Calculate particle trajectories
460  !********************************
461
462  if (verbosity.gt.0.and. lroot) then
463    if (verbosity.gt.1) then   
464      CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
465      WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
466    endif
467    print*,'call timemanager'
468  endif
469  if (info_flag.eq.1) then
470    print*, 'info only mode (stop)'   
471    call mpif_finalize
472    stop
473  endif
474
475
476  call timemanager(metdata_format)
477
478
479! NIK 16.02.2005
480  if (mp_partgroup_pid.ge.0) then ! Skip for readwind process
481    call MPI_Reduce(tot_blc_count, tot_b, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
482         & mp_comm_used, mp_ierr)
483    call MPI_Reduce(tot_inc_count, tot_i, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
484         & mp_comm_used, mp_ierr)
485  end if
486 
487
488  if (lroot) then
489    do i=1,nspec
490      write(*,*) '**********************************************'
491      write(*,*) 'Scavenging statistics for species ', species(i), ':'
492      write(*,*) 'Total number of occurences of below-cloud scavenging', &
493           & tot_b(i)
494!           & tot_blc_count(i)
495      write(*,*) 'Total number of occurences of in-cloud    scavenging', &
496           & tot_i(i)
497!           & tot_inc_count(i)
498      write(*,*) '**********************************************'
499    end do
500
501    write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
502         &XPART MODEL RUN!'
503  end if
504
505  if (mp_measure_time) call mpif_mtime('flexpart',1)
506
507  call mpif_finalize
508
509end program flexpart
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG