source: flexpart.git/src/FLEXPART_MPI.f90 @ 5184a7c

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

Enable settling with multiple species if from separate releases

  • 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. '//trim(flexversion_major)//'.1beta MPI (2016-11-02)'
83  verbosity=0
84
85  ! Read the pathnames where input/output files are stored
86  !*******************************************************
87
88  inline_options='none'
89  select case (iargc())
90  case (2)
91    call getarg(1,arg1)
92    pathfile=arg1
93    call getarg(2,arg2)
94    inline_options=arg2
95  case (1)
96    call getarg(1,arg1)
97    pathfile=arg1
98    if (arg1(1:1).eq.'-') then
99      write(pathfile,'(a11)') './pathnames'
100      inline_options=arg1
101    endif
102  case (0)
103    write(pathfile,'(a11)') './pathnames'
104  end select
105 
106  if (lroot) then
107  ! Print the GPL License statement
108  !*******************************************************
109    print*,'Welcome to FLEXPART ', trim(flexversion)
110    print*,'FLEXPART is free software released under the GNU General Public License.'
111  endif
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 .and. lroot) 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  ! Exit if trying to run backwards
167  if (ldirect.le.0) then
168    write(*,FMT='(80("#"))')
169    write(*,*) '#### FLEXPART_MPI> ERROR: ', &
170         & 'MPI version not (yet) working with backward runs. '
171    write(*,*) '#### Use the serial version instead.'
172    write(*,FMT='(80("#"))')
173    ! call mpif_finalize
174    ! stop
175  end if
176
177
178
179! Read the age classes to be used
180!********************************
181  if (verbosity.gt.0 .and. lroot) then
182    write(*,*) 'call readageclasses'
183  endif
184  call readageclasses
185
186  if (verbosity.gt.1 .and. lroot) then   
187    call system_clock(count_clock, count_rate, count_max)
188    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
189  endif     
190
191  ! Read, which wind fields are available within the modelling period
192  !******************************************************************
193
194  if (verbosity.gt.0 .and. lroot) then
195    write(*,*) 'call readavailable'
196  endif 
197  call readavailable
198
199  ! If nested wind fields are used, allocate arrays
200  !************************************************
201
202  if (verbosity.gt.0 .and. lroot) then
203    write(*,*) 'call com_mod_allocate_nests'
204  endif
205  call com_mod_allocate_nests
206
207! Read the model grid specifications,
208! both for the mother domain and eventual nests
209!**********************************************
210
211  if (verbosity.gt.0 .and. lroot) then
212    write(*,*) 'call gridcheck'
213  endif
214
215  call gridcheck
216
217  if (verbosity.gt.1 .and. lroot) then   
218    call system_clock(count_clock, count_rate, count_max)
219    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
220  endif     
221
222
223  if (verbosity.gt.0 .and. lroot) then
224    write(*,*) 'call gridcheck_nests'
225  endif 
226  call gridcheck_nests
227
228
229! Read the output grid specifications
230!************************************
231
232  if (verbosity.gt.0 .and. lroot) then
233    WRITE(*,*) 'call readoutgrid'
234  endif
235
236  call readoutgrid
237
238  if (nested_output.eq.1) then
239    call readoutgrid_nest
240    if (verbosity.gt.0.and.lroot) then
241      WRITE(*,*) '# readoutgrid_nest'
242    endif
243  endif
244
245  ! Read the receptor points for which extra concentrations are to be calculated
246  !*****************************************************************************
247
248  if (verbosity.eq.1 .and. lroot) then
249    print*,'call readreceptors'
250  endif
251  call readreceptors
252
253  ! Read the physico-chemical species property table
254  !*************************************************
255  !SEC: now only needed SPECIES are read in readreleases.f
256  !call readspecies
257
258
259  ! Read the landuse inventory
260  !***************************
261
262  if (verbosity.gt.0 .and. lroot) then
263    print*,'call readlanduse'
264  endif
265  call readlanduse
266
267
268! Assign fractional cover of landuse classes to each ECMWF grid point
269!********************************************************************
270
271  if (verbosity.gt.0 .and. lroot) then
272    print*,'call assignland'
273  endif
274  call assignland
275
276  ! Read the coordinates of the release locations
277  !**********************************************
278
279  if (verbosity.gt.0 .and. lroot) then
280    print*,'call readreleases'
281  endif
282  call readreleases
283
284
285! Read and compute surface resistances to dry deposition of gases
286!****************************************************************
287
288  if (verbosity.gt.0 .and. lroot) then
289    print*,'call readdepo'
290  endif
291  call readdepo
292
293  ! Convert the release point coordinates from geografical to grid coordinates
294  !***************************************************************************
295
296  call coordtrafo 
297  if (verbosity.gt.0 .and. lroot) then
298    print*,'call coordtrafo'
299  endif
300
301
302  ! Initialize all particles to non-existent
303  !*****************************************
304
305  if (verbosity.gt.0 .and. lroot) then
306    print*,'Initialize all particles to non-existent'
307  endif
308
309  if (.not.(lmpreader.and.lmp_use_reader)) then
310    do j=1, size(itra1) ! maxpart_mpi
311      itra1(j)=-999999999
312    end do
313  end if
314
315  ! For continuation of previous run, read in particle positions
316  !*************************************************************
317
318  if (ipin.eq.1) then
319    if (verbosity.gt.0 .and. lroot) then
320      print*,'call readpartpositions'
321    endif
322    ! readwind process skips this step
323    if (.not.(lmpreader.and.lmp_use_reader)) call readpartpositions
324  else
325    if (verbosity.gt.0 .and. lroot) then
326      print*,'numpart=0, numparticlecount=0'
327    endif
328    numpart=0
329    numparticlecount=0
330  endif
331
332
333  ! Calculate volume, surface area, etc., of all output grid cells
334  ! Allocate fluxes and OHfield if necessary
335  !***************************************************************
336
337
338  if (verbosity.gt.0 .and. lroot) then
339    print*,'call outgrid_init'
340  endif
341  call outgrid_init
342  if (nested_output.eq.1) call outgrid_init_nest
343
344  ! Read the OH field
345  !******************
346
347  if (OHREA.eqv..true.) then
348    if (verbosity.gt.0 .and. lroot) then
349      print*,'call readOHfield'
350    endif
351    call readOHfield
352  endif
353
354  ! Write basic information on the simulation to a file "header"
355  ! and open files that are to be kept open throughout the simulation
356  !******************************************************************
357
358  if (mp_measure_time) call mpif_mtime('iotime',0)
359  if (lroot) then ! MPI: this part root process only
360
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
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  end if ! (mpif_pid == 0)
388
389  if (mp_measure_time) call mpif_mtime('iotime',0)
390
391  !open(unitdates,file=path(2)(1:length(2))//'dates')
392
393  if (verbosity.gt.0 .and. lroot) then
394    print*,'call openreceptors'
395  endif
396  call openreceptors
397  if ((iout.eq.4).or.(iout.eq.5)) call openouttraj
398
399  ! Releases can only start and end at discrete times (multiples of lsynctime)
400  !***************************************************************************
401
402  if (verbosity.gt.0 .and. lroot) then
403    print*,'discretize release times'
404  endif
405  do i=1,numpoint
406    ireleasestart(i)=nint(real(ireleasestart(i))/real(lsynctime))*lsynctime
407    ireleaseend(i)=nint(real(ireleaseend(i))/real(lsynctime))*lsynctime
408  end do
409
410  ! Initialize cloud-base mass fluxes for the convection scheme
411  !************************************************************
412
413  if (verbosity.gt.0 .and. lroot) then
414    print*,'Initialize cloud-base mass fluxes for the convection scheme'
415  endif
416
417  do jy=0,nymin1
418    do ix=0,nxmin1
419      cbaseflux(ix,jy)=0.
420    end do
421  end do
422  do inest=1,numbnests
423    do jy=0,nyn(inest)-1
424      do ix=0,nxn(inest)-1
425        cbasefluxn(ix,jy,inest)=0.
426      end do
427    end do
428  end do
429
430  ! Inform whether output kernel is used or not
431  !*********************************************
432
433  if (lroot) then
434    if (lnokernel) then
435      write(*,*) "Concentrations are calculated without using kernel"
436    else
437      write(*,*) "Concentrations are calculated using kernel"
438    end if
439  end if
440
441! Calculate particle trajectories
442!********************************
443
444  if (verbosity.gt.0.and. lroot) then
445    if (verbosity.gt.1) then   
446      CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
447      WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
448    endif
449    print*,'call timemanager'
450  endif
451  if (info_flag.eq.1) then
452    print*, 'info only mode (stop)'   
453    call mpif_finalize
454    stop
455  endif
456
457
458  call timemanager
459
460
461! NIK 16.02.2005
462  if (lroot) then
463    call MPI_Reduce(MPI_IN_PLACE, tot_blc_count, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
464         & mp_comm_used, mp_ierr)
465    call MPI_Reduce(MPI_IN_PLACE, tot_inc_count, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
466         & mp_comm_used, mp_ierr)
467  else
468    if (mp_partgroup_pid.ge.0) then ! Skip for readwind process
469      call MPI_Reduce(tot_blc_count, 0, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
470           & mp_comm_used, mp_ierr)
471      call MPI_Reduce(tot_inc_count, 0, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
472           & mp_comm_used, mp_ierr)
473    end if
474  end if
475
476  if (lroot) then
477    do i=1,nspec
478      write(*,*) '**********************************************'
479      write(*,*) 'Scavenging statistics for species ', species(i), ':'
480      write(*,*) 'Total number of occurences of below-cloud scavenging', &
481           & tot_blc_count(i)
482      write(*,*) 'Total number of occurences of in-cloud    scavenging', &
483           & tot_inc_count(i)
484      write(*,*) '**********************************************'
485    end do
486
487    write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
488         &XPART MODEL RUN!'
489  end if
490
491  if (mp_measure_time) call mpif_mtime('flexpart',1)
492
493  call mpif_finalize
494
495end program flexpart
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG