source: flexpart.git/src/FLEXPART_MPI.f90 @ 3481cc1

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file since 3481cc1 was 3481cc1, checked in by Ignacio Pisso <ip@…>, 4 years ago

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