FLEXPART_MPI.f90 14.8 KB
Newer Older
1 2
! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
! SPDX-License-Identifier: GPL-3.0-or-later
3

4 5 6 7 8 9 10 11 12 13 14 15 16
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                                                            *
  !                                                                            *
  !*****************************************************************************
17 18 19 20 21 22 23 24
  ! 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                            *
  !*****************************************************************************
25 26 27 28 29 30 31 32 33 34 35 36 37
  !                                                                            *
  ! Variables:                                                                 *
  !                                                                            *
  ! Constants:                                                                 *
  !                                                                            *
  !*****************************************************************************

  use point_mod
  use par_mod
  use com_mod
  use conv_mod
  use mpi_mod
  use random_mod, only: gasdev1
38
  use class_gribfile
39

40 41 42 43
#ifdef USE_NCF
  use netcdf_output_mod, only: writeheader_netcdf
#endif

44 45 46 47 48
  implicit none

  integer :: i,j,ix,jy,inest
  integer :: idummy = -320
  character(len=256) :: inline_options  !pathfile, flexversion, arg2
49 50
  integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN
  integer :: detectformat
51 52
  integer(selected_int_kind(16)), dimension(maxspec) :: tot_b=0, &
       & tot_i=0
53 54


55 56
  ! Initialize mpi
  !*********************
57 58 59 60
  call mpif_init

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

61
  
62 63 64 65 66 67 68 69 70 71 72 73
  ! 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
74
  flexversion_major = '10' ! Major version number, also used for species file names
Espen Sollum's avatar
Espen Sollum committed
75
  flexversion='Ver. '//trim(flexversion_major)//'.2beta MPI (2017-08-01)'
76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 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
  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 
130
  call readpaths
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
 
  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
156
    endif
157 158
  endif

159 160 161 162 163
  ! Initialize arrays in com_mod 
  !*****************************

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

164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183

! 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
184

185 186 187 188 189 190
  ! Detect metdata format
  !**********************

  metdata_format = detectformat()

  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
191
    if (lroot) print *,'ECMWF metdata detected'
192
  elseif (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
193
    if (lroot) print *,'NCEP metdata detected'
194
  else
195
    if (lroot) print *,'Unknown metdata format'
196
    stop
197 198 199 200
  endif



201 202 203 204 205 206 207
  ! 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
208 209 210 211 212 213 214 215 216

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

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

217 218 219 220 221 222
  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
    call gridcheck_ecmwf
  else
    call gridcheck_gfs
  end if

223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 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

  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

316
  if (.not.(lmpreader.and.lmp_use_reader)) then
317 318 319 320
    do j=1, size(itra1) ! maxpart_mpi
      itra1(j)=-999999999
    end do
  end if
321 322 323 324 325 326 327 328

  ! 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
329
    ! readwind process skips this step
330
    if (.not.(lmpreader.and.lmp_use_reader)) call readpartpositions
331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364
  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
  !******************************************************************

365
  if (mp_measure_time) call mpif_mtime('iotime',0)
366

367 368 369 370 371 372 373 374 375 376 377 378 379 380
  if (lroot) then ! MPI: this part root process only
#ifdef USE_NCF
    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
381
    endif
382
#endif
383 384 385 386 387 388 389 390

    if (verbosity.gt.0) then
      print*,'call writeheader'
    endif

    call writeheader
! FLEXPART 9.2 ticket ?? write header in ASCII format 
    call writeheader_txt
391

392 393 394 395 396
    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) 

397
  if (mp_measure_time) call mpif_mtime('iotime',1)
398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435

  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

436 437 438 439
  ! Inform whether output kernel is used or not
  !*********************************************

  if (lroot) then
440
    if (.not.lusekerneloutput) then
441 442 443 444 445
      write(*,*) "Concentrations are calculated without using kernel"
    else
      write(*,*) "Concentrations are calculated using kernel"
    end if
  end if
446 447 448 449 450 451 452 453 454 455 456

! 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
457 458 459 460 461 462
  if (info_flag.eq.1) then
    print*, 'info only mode (stop)'    
    call mpif_finalize
    stop
  endif

463

464
  call timemanager(metdata_format)
465

466

467
! NIK 16.02.2005 
468 469
  if (mp_partgroup_pid.ge.0) then ! Skip for readwind process 
    call MPI_Reduce(tot_blc_count, tot_b, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
Espen Sollum's avatar
Espen Sollum committed
470
         & mp_comm_used, mp_ierr)
471
    call MPI_Reduce(tot_inc_count, tot_i, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
Espen Sollum's avatar
Espen Sollum committed
472 473
         & mp_comm_used, mp_ierr)
  end if
474
  
475 476

  if (lroot) then
477 478 479 480
    do i=1,nspec
      write(*,*) '**********************************************'
      write(*,*) 'Scavenging statistics for species ', species(i), ':'
      write(*,*) 'Total number of occurences of below-cloud scavenging', &
481 482
           & tot_b(i)
!           & tot_blc_count(i)
483
      write(*,*) 'Total number of occurences of in-cloud    scavenging', &
484 485
           & tot_i(i)
!           & tot_inc_count(i)
486 487
      write(*,*) '**********************************************'
    end do
Espen Sollum's avatar
Espen Sollum committed
488

489 490 491 492 493 494 495 496 497
    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