source: flexpart.git/src/FLEXPART_MPI.f90 @ 20963b1

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

Backwards deposition for the MPI version

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