FLEXPART_MPI.f90 16.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
!                                                                     *
! This file is part of FLEXPART.                                      *
!                                                                     *
! FLEXPART is free software: you can redistribute it and/or modify    *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or   *
! (at your option) any later version.                                 *
!                                                                     *
! FLEXPART is distributed in the hope that it will be useful,         *
! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
! GNU General Public License for more details.                        *
!                                                                     *
! You should have received a copy of the GNU General Public License   *
! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
!**********************************************************************

program flexpart

  !*****************************************************************************
  !                                                                            *
  !     This is the Lagrangian Particle Dispersion Model FLEXPART.             *
  !     The main program manages the reading of model run specifications, etc. *
  !     All actual computing is done within subroutine timemanager.            *
  !                                                                            *
  !     Author: A. Stohl                                                       *
  !                                                                            *
  !     18 May 1996                                                            *
  !                                                                            *
  !*****************************************************************************
35
36
37
38
39
40
41
42
  ! Changes:                                                                   *
  !   Unified ECMWF and GFS builds                                             *
  !   Marian Harustak, 12.5.2017                                               *
  !     - Added detection of metdata format using gributils routines           *
  !     - Distinguished calls to ecmwf/gfs gridcheck versions based on         *
  !       detected metdata format                                              *
  !     - Passed metdata format down to timemanager                            *
  !*****************************************************************************
43
44
45
46
47
48
49
50
51
52
53
54
55
56
  !                                                                            *
  ! Variables:                                                                 *
  !                                                                            *
  ! Constants:                                                                 *
  !                                                                            *
  !*****************************************************************************

  use point_mod
  use par_mod
  use com_mod
  use conv_mod
  use mpi_mod
  use netcdf_output_mod, only: writeheader_netcdf
  use random_mod, only: gasdev1
57
  use class_gribfile
58
59
60
61
62
63

  implicit none

  integer :: i,j,ix,jy,inest
  integer :: idummy = -320
  character(len=256) :: inline_options  !pathfile, flexversion, arg2
64
65
66
  integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN
  integer :: detectformat

67
68


69
70
  ! Initialize mpi
  !*********************
71
72
73
74
  call mpif_init

  if (mp_measure_time) call mpif_mtime('flexpart',0)

75
76
  ! Initialize arrays in com_mod 
  !*****************************
77
78
79

  if(.not.(lmpreader.and.lmp_use_reader)) call com_mod_allocate_part(maxpart_mpi)

80
81
82
83
84
85
86
87
88
89
90
91
92

  ! Generate a large number of random numbers
  !******************************************

  ! eso: Different seed for each MPI process
  idummy=idummy+mp_seed

  do i=1,maxrand-1,2
    call gasdev1(idummy,rannumb(i),rannumb(i+1))
  end do
  call gasdev1(idummy,rannumb(maxrand),rannumb(maxrand-1))

  ! 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
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
  verbosity=0

  ! Read the pathnames where input/output files are stored
  !*******************************************************

  inline_options='none'
  select case (iargc())
  case (2)
    call getarg(1,arg1)
    pathfile=arg1
    call getarg(2,arg2)
    inline_options=arg2
  case (1)
    call getarg(1,arg1)
    pathfile=arg1
    if (arg1(1:1).eq.'-') then
      write(pathfile,'(a11)') './pathnames'
      inline_options=arg1 
    endif
  case (0)
    write(pathfile,'(a11)') './pathnames'
  end select
  
  if (lroot) then
  ! Print the GPL License statement
  !*******************************************************
    print*,'Welcome to FLEXPART ', trim(flexversion)
    print*,'FLEXPART is free software released under the GNU General Public License.'
  endif
  
  if (inline_options(1:1).eq.'-') then
    if (trim(inline_options).eq.'-v'.or.trim(inline_options).eq.'-v1') then
      print*, 'Verbose mode 1: display detailed information during run'
      verbosity=1
    endif
    if (trim(inline_options).eq.'-v2') then
      print*, 'Verbose mode 2: display more detailed information during run'
      verbosity=2
    endif
    if (trim(inline_options).eq.'-i') then
      print*, 'Info mode: provide detailed run specific information and stop'
      verbosity=1
      info_flag=1
    endif
    if (trim(inline_options).eq.'-i2') then
      print*, 'Info mode: provide more detailed run specific information and stop'
      verbosity=2
      info_flag=1
    endif
  endif
           
  if (verbosity.gt.0) then
    write(*,*) 'call readpaths'
  endif 
  call readpaths(pathfile)
 
  if (verbosity.gt.1) then !show clock info 
     !print*,'length(4)',length(4)
     !count=0,count_rate=1000
    call system_clock(count_clock0, count_rate, count_max)
     !WRITE(*,*) 'SYSTEM_CLOCK',count, count_rate, count_max
     !WRITE(*,*) 'SYSTEM_CLOCK, count_clock0', count_clock0
     !WRITE(*,*) 'SYSTEM_CLOCK, count_rate', count_rate
     !WRITE(*,*) 'SYSTEM_CLOCK, count_max', count_max
  endif

  ! Read the user specifications for the current model run
  !*******************************************************

  if (verbosity.gt.0) then
    write(*,*) 'call readcommand'
  endif
  call readcommand
  if (verbosity.gt.0 .and. lroot) then
    write(*,*) '    ldirect=', ldirect
    write(*,*) '    ibdate,ibtime=',ibdate,ibtime
    write(*,*) '    iedate,ietime=', iedate,ietime
    if (verbosity.gt.1) then   
      CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
      write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
175
    endif
176
177
  endif

178
179
180
181
  ! Exit if trying to run backwards
  if (ldirect.le.0) then
    write(*,FMT='(80("#"))')
    write(*,*) '#### FLEXPART_MPI> ERROR: ', &
182
183
         & 'MPI version not (yet) working with backward runs. '
    write(*,*) '#### Use the serial version instead.'
184
    write(*,FMT='(80("#"))')
185
186
    ! call mpif_finalize
    ! stop
187
188
189
  end if


190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209

! Read the age classes to be used
!********************************
  if (verbosity.gt.0 .and. lroot) then
    write(*,*) 'call readageclasses'
  endif
  call readageclasses

  if (verbosity.gt.1 .and. lroot) then   
    call system_clock(count_clock, count_rate, count_max)
    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
  endif     

  ! Read, which wind fields are available within the modelling period
  !******************************************************************

  if (verbosity.gt.0 .and. lroot) then
    write(*,*) 'call readavailable'
  endif  
  call readavailable
210

211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
  ! Detect metdata format
  !**********************

  metdata_format = detectformat()

  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
    print *,'ECMWF metdata detected'
  elseif (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
    print *,'NCEP metdata detected'
  else
    print *,'Unknown metdata format'
    return
  endif



227
228
229
230
231
232
233
  ! If nested wind fields are used, allocate arrays
  !************************************************

  if (verbosity.gt.0 .and. lroot) then
    write(*,*) 'call com_mod_allocate_nests'
  endif
  call com_mod_allocate_nests
234
235
236
237
238
239
240
241
242

! Read the model grid specifications,
! both for the mother domain and eventual nests
!**********************************************

  if (verbosity.gt.0 .and. lroot) then
    write(*,*) 'call gridcheck'
  endif

243
244
245
246
247
248
  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
    call gridcheck_ecmwf
  else
    call gridcheck_gfs
  end if

249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341

  if (verbosity.gt.1 .and. lroot) then   
    call system_clock(count_clock, count_rate, count_max)
    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
  endif      


  if (verbosity.gt.0 .and. lroot) then
    write(*,*) 'call gridcheck_nests'
  endif  
  call gridcheck_nests


! Read the output grid specifications
!************************************

  if (verbosity.gt.0 .and. lroot) then
    WRITE(*,*) 'call readoutgrid'
  endif

  call readoutgrid

  if (nested_output.eq.1) then
    call readoutgrid_nest
    if (verbosity.gt.0.and.lroot) then
      WRITE(*,*) '# readoutgrid_nest'
    endif
  endif

  ! Read the receptor points for which extra concentrations are to be calculated
  !*****************************************************************************

  if (verbosity.eq.1 .and. lroot) then
    print*,'call readreceptors'
  endif
  call readreceptors

  ! Read the physico-chemical species property table
  !*************************************************
  !SEC: now only needed SPECIES are read in readreleases.f
  !call readspecies


  ! Read the landuse inventory
  !***************************

  if (verbosity.gt.0 .and. lroot) then
    print*,'call readlanduse'
  endif
  call readlanduse


! Assign fractional cover of landuse classes to each ECMWF grid point
!********************************************************************

  if (verbosity.gt.0 .and. lroot) then
    print*,'call assignland'
  endif
  call assignland

  ! Read the coordinates of the release locations
  !**********************************************

  if (verbosity.gt.0 .and. lroot) then
    print*,'call readreleases'
  endif
  call readreleases


! Read and compute surface resistances to dry deposition of gases
!****************************************************************

  if (verbosity.gt.0 .and. lroot) then
    print*,'call readdepo'
  endif
  call readdepo

  ! Convert the release point coordinates from geografical to grid coordinates
  !***************************************************************************

  call coordtrafo  
  if (verbosity.gt.0 .and. lroot) then
    print*,'call coordtrafo'
  endif


  ! Initialize all particles to non-existent
  !*****************************************

  if (verbosity.gt.0 .and. lroot) then
    print*,'Initialize all particles to non-existent'
  endif

342
  if (.not.(lmpreader.and.lmp_use_reader)) then
343
344
345
346
    do j=1, size(itra1) ! maxpart_mpi
      itra1(j)=-999999999
    end do
  end if
347
348
349
350
351
352
353
354

  ! For continuation of previous run, read in particle positions
  !*************************************************************

  if (ipin.eq.1) then
    if (verbosity.gt.0 .and. lroot) then
      print*,'call readpartpositions'
    endif
355
    ! readwind process skips this step
356
    if (.not.(lmpreader.and.lmp_use_reader)) call readpartpositions
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
  else
    if (verbosity.gt.0 .and. lroot) then
      print*,'numpart=0, numparticlecount=0'
    endif
    numpart=0
    numparticlecount=0
  endif


  ! Calculate volume, surface area, etc., of all output grid cells
  ! Allocate fluxes and OHfield if necessary
  !***************************************************************


  if (verbosity.gt.0 .and. lroot) then
    print*,'call outgrid_init'
  endif
  call outgrid_init
  if (nested_output.eq.1) call outgrid_init_nest

  ! Read the OH field
  !******************

  if (OHREA.eqv..true.) then
    if (verbosity.gt.0 .and. lroot) then
      print*,'call readOHfield'
    endif
    call readOHfield
  endif

  ! Write basic information on the simulation to a file "header"
  ! and open files that are to be kept open throughout the simulation
  !******************************************************************

391
392
393
  if (mp_measure_time) call mpif_mtime('iotime',0)
  if (lroot) then ! MPI: this part root process only

394
395
396
397
398
399
400
401
402
403
404
405
406
407
  if (lnetcdfout.eq.1) then 
    call writeheader_netcdf(lnest=.false.)
  else 
    call writeheader
  end if

  if (nested_output.eq.1) then
    if (lnetcdfout.eq.1) then
      call writeheader_netcdf(lnest=.true.)
    else
      call writeheader_nest
    endif
  endif

408
!
409
410
411
412
413
414
415
416
417
418
419
420
421
    if (verbosity.gt.0) then
      print*,'call writeheader'
    endif

    call writeheader
! FLEXPART 9.2 ticket ?? write header in ASCII format 
    call writeheader_txt
!if (nested_output.eq.1) call writeheader_nest
    if (nested_output.eq.1.and.surf_only.ne.1) call writeheader_nest
    if (nested_output.eq.1.and.surf_only.eq.1) call writeheader_nest_surf
    if (nested_output.ne.1.and.surf_only.eq.1) call writeheader_surf
  end if ! (mpif_pid == 0) 

422
  if (mp_measure_time) call mpif_mtime('iotime',0)
423

424
  !open(unitdates,file=path(2)(1:length(2))//'dates')
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462

  if (verbosity.gt.0 .and. lroot) then
    print*,'call openreceptors'
  endif
  call openreceptors
  if ((iout.eq.4).or.(iout.eq.5)) call openouttraj

  ! Releases can only start and end at discrete times (multiples of lsynctime)
  !***************************************************************************

  if (verbosity.gt.0 .and. lroot) then
    print*,'discretize release times'
  endif
  do i=1,numpoint
    ireleasestart(i)=nint(real(ireleasestart(i))/real(lsynctime))*lsynctime
    ireleaseend(i)=nint(real(ireleaseend(i))/real(lsynctime))*lsynctime
  end do

  ! Initialize cloud-base mass fluxes for the convection scheme
  !************************************************************

  if (verbosity.gt.0 .and. lroot) then
    print*,'Initialize cloud-base mass fluxes for the convection scheme'
  endif

  do jy=0,nymin1
    do ix=0,nxmin1
      cbaseflux(ix,jy)=0.
    end do
  end do
  do inest=1,numbnests
    do jy=0,nyn(inest)-1
      do ix=0,nxn(inest)-1
        cbasefluxn(ix,jy,inest)=0.
      end do
    end do
  end do

463
464
465
466
467
468
469
470
471
472
  ! Inform whether output kernel is used or not
  !*********************************************

  if (lroot) then
    if (lnokernel) then
      write(*,*) "Concentrations are calculated without using kernel"
    else
      write(*,*) "Concentrations are calculated using kernel"
    end if
  end if
473
474
475
476
477
478
479
480
481
482
483

! Calculate particle trajectories
!********************************

  if (verbosity.gt.0.and. lroot) then
    if (verbosity.gt.1) then   
      CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
      WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
    endif
    print*,'call timemanager'
  endif
484
485
486
487
488
489
  if (info_flag.eq.1) then
    print*, 'info only mode (stop)'    
    call mpif_finalize
    stop
  endif

490

491
  call timemanager(metdata_format)
492

493

494
! NIK 16.02.2005 
Espen Sollum's avatar
Espen Sollum committed
495
  if (lroot) then
496
    call MPI_Reduce(MPI_IN_PLACE, tot_blc_count, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
Espen Sollum's avatar
Espen Sollum committed
497
         & mp_comm_used, mp_ierr)
498
    call MPI_Reduce(MPI_IN_PLACE, tot_inc_count, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
Espen Sollum's avatar
Espen Sollum committed
499
500
         & mp_comm_used, mp_ierr)
  else
Espen Sollum's avatar
Espen Sollum committed
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, &
Espen Sollum's avatar
Espen Sollum committed
503
           & mp_comm_used, mp_ierr)
504
      call MPI_Reduce(tot_inc_count, 0, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
Espen Sollum's avatar
Espen Sollum committed
505
506
           & mp_comm_used, mp_ierr)
    end if
Espen Sollum's avatar
Espen Sollum committed
507
  end if
508
509

  if (lroot) then
510
511
512
513
514
515
516
517
518
    do i=1,nspec
      write(*,*) '**********************************************'
      write(*,*) 'Scavenging statistics for species ', species(i), ':'
      write(*,*) 'Total number of occurences of below-cloud scavenging', &
           & tot_blc_count(i)
      write(*,*) 'Total number of occurences of in-cloud    scavenging', &
           & tot_inc_count(i)
      write(*,*) '**********************************************'
    end do
Espen Sollum's avatar
Espen Sollum committed
519

520
521
522
523
524
525
526
527
528
    write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
         &XPART MODEL RUN!'
  end if

  if (mp_measure_time) call mpif_mtime('flexpart',1)

  call mpif_finalize

end program flexpart