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

univie
Last change on this file since c0884a8 was c0884a8, checked in by pesei <petra seibert at univie ac at>, 6 years ago

replace CTBTO code for checking type of GRIB

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