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

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