source: flexpart.git/src/FLEXPART.f90

bugfixes+enhancements
Last change on this file was 8ad70c7, checked in by Pirmin Kaufmann <pirmin.kaufmann@…>, 16 months ago

Bumped bugfix version.

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