source: flexpart.git/src/FLEXPART_MPI.f90 @ 6ecb30a

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

Merged changes from CTBTO project

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