source: flexpart.git/src/FLEXPART.f90

10.4.1_pesei
Last change on this file was 92fab65, checked in by Ignacio Pisso <ip@…>, 8 months ago

add SPDX-License-Identifier to all .f90 files

  • Property mode set to 100644
File size: 14.0 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
[07c3e71]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
[adead08]64  flexversion='Version '//trim(flexversion_major)//'.4 (2019-11-12)'
[8a65cb0]65  verbosity=0
66
[f13406c]67  ! Read the pathnames where input/output files are stored
68  !*******************************************************
69
[8a65cb0]70  inline_options='none'
[f13406c]71  select case (iargc())
[8a65cb0]72  case (2)
[f13406c]73    call getarg(1,arg1)
74    pathfile=arg1
75    call getarg(2,arg2)
76    inline_options=arg2
[8a65cb0]77  case (1)
[f13406c]78    call getarg(1,arg1)
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.'
[414a5e5]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
114       print*, 'Info mode: provide detailed run specific information and stop'
[b4d29ce]115       verbosity=1
116       info_flag=1
117    endif
[8a65cb0]118    if (trim(inline_options).eq.'-i2') then
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
[8a65cb0]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
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  !**********************************************
[f13406c]219 
220  if (verbosity.gt.0) then
[8a65cb0]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
[8a65cb0]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
270  ! Read the landuse inventory
271  !***************************
272
[f13406c]273  if (verbosity.gt.0) then
[8a65cb0]274    print*,'call readlanduse'
[f13406c]275  endif
[e200b7a]276  call readlanduse
277
278  ! Assign fractional cover of landuse classes to each ECMWF grid point
279  !********************************************************************
280
[f13406c]281  if (verbosity.gt.0) then
[8a65cb0]282    print*,'call assignland'
[f13406c]283  endif
[e200b7a]284  call assignland
285
286  ! Read the coordinates of the release locations
287  !**********************************************
288
[f13406c]289  if (verbosity.gt.0) then
[8a65cb0]290    print*,'call readreleases'
[f13406c]291  endif
[e200b7a]292  call readreleases
293
294  ! Read and compute surface resistances to dry deposition of gases
295  !****************************************************************
296
[f13406c]297  if (verbosity.gt.0) then
[8a65cb0]298    print*,'call readdepo'
[f13406c]299  endif
[e200b7a]300  call readdepo
301
302  ! Convert the release point coordinates from geografical to grid coordinates
303  !***************************************************************************
304
[f13406c]305  call coordtrafo 
306  if (verbosity.gt.0) then
[8a65cb0]307    print*,'call coordtrafo'
[f13406c]308  endif
[e200b7a]309
310  ! Initialize all particles to non-existent
311  !*****************************************
312
[f13406c]313  if (verbosity.gt.0) then
[8a65cb0]314    print*,'Initialize all particles to non-existent'
[f13406c]315  endif
[e200b7a]316  do j=1,maxpart
317    itra1(j)=-999999999
318  end do
319
320  ! For continuation of previous run, read in particle positions
321  !*************************************************************
322
323  if (ipin.eq.1) then
[f13406c]324    if (verbosity.gt.0) then
[8a65cb0]325      print*,'call readpartpositions'
[f13406c]326    endif
[e200b7a]327    call readpartpositions
328  else
[f13406c]329    if (verbosity.gt.0) then
[8a65cb0]330      print*,'numpart=0, numparticlecount=0'
[f13406c]331    endif   
[e200b7a]332    numpart=0
333    numparticlecount=0
334  endif
335
336  ! Calculate volume, surface area, etc., of all output grid cells
337  ! Allocate fluxes and OHfield if necessary
338  !***************************************************************
339
[f13406c]340  if (verbosity.gt.0) then
[8a65cb0]341    print*,'call outgrid_init'
[f13406c]342  endif
[e200b7a]343  call outgrid_init
344  if (nested_output.eq.1) call outgrid_init_nest
345
346  ! Read the OH field
347  !******************
348
[f13406c]349  if (OHREA.eqv..TRUE.) then
350    if (verbosity.gt.0) then
[8a65cb0]351      print*,'call readOHfield'
[f13406c]352    endif
[b4d29ce]353    call readOHfield
[f13406c]354  endif
[e200b7a]355
356  ! Write basic information on the simulation to a file "header"
357  ! and open files that are to be kept open throughout the simulation
358  !******************************************************************
359
[a9cf4b1]360#ifdef USE_NCF
[8a65cb0]361  if (lnetcdfout.eq.1) then
362    call writeheader_netcdf(lnest=.false.)
363  else
364    call writeheader
365  end if
366
367  if (nested_output.eq.1) then
368    if (lnetcdfout.eq.1) then
369      call writeheader_netcdf(lnest=.true.)
370    else
371      call writeheader_nest
372    endif
373  endif
[a9cf4b1]374#endif
[8a65cb0]375
[f13406c]376  if (verbosity.gt.0) then
[8a65cb0]377    print*,'call writeheader'
[f13406c]378  endif
379
[e200b7a]380  call writeheader
[8a65cb0]381  ! FLEXPART 9.2 ticket ?? write header in ASCII format
[f13406c]382  call writeheader_txt
383  !if (nested_output.eq.1) call writeheader_nest
384  if (nested_output.eq.1.and.surf_only.ne.1) call writeheader_nest
385  if (nested_output.eq.1.and.surf_only.eq.1) call writeheader_nest_surf
386  if (nested_output.ne.1.and.surf_only.eq.1) call writeheader_surf
387
388  !open(unitdates,file=path(2)(1:length(2))//'dates')
389
390  if (verbosity.gt.0) then
[8a65cb0]391    print*,'call openreceptors'
[f13406c]392  endif
[e200b7a]393  call openreceptors
394  if ((iout.eq.4).or.(iout.eq.5)) call openouttraj
395
396  ! Releases can only start and end at discrete times (multiples of lsynctime)
397  !***************************************************************************
398
[f13406c]399  if (verbosity.gt.0) then
[8a65cb0]400    print*,'discretize release times'
[f13406c]401  endif
[e200b7a]402  do i=1,numpoint
[b4d29ce]403    ireleasestart(i)=nint(real(ireleasestart(i))/real(lsynctime))*lsynctime
404    ireleaseend(i)=nint(real(ireleaseend(i))/real(lsynctime))*lsynctime
[e200b7a]405  end do
406
407  ! Initialize cloud-base mass fluxes for the convection scheme
408  !************************************************************
409
[f13406c]410  if (verbosity.gt.0) then
[8a65cb0]411    print*,'Initialize cloud-base mass fluxes for the convection scheme'
[f13406c]412  endif
413
[e200b7a]414  do jy=0,nymin1
415    do ix=0,nxmin1
416      cbaseflux(ix,jy)=0.
417    end do
418  end do
419  do inest=1,numbnests
420    do jy=0,nyn(inest)-1
421      do ix=0,nxn(inest)-1
422        cbasefluxn(ix,jy,inest)=0.
423      end do
424    end do
425  end do
426
[4c64400]427  ! Inform whether output kernel is used or not
428  !*********************************************
429  if (lroot) then
[fe32dca]430    if (.not.lusekerneloutput) then
[4c64400]431      write(*,*) "Concentrations are calculated without using kernel"
432    else
433      write(*,*) "Concentrations are calculated using kernel"
434    end if
435  end if
436
437
[e200b7a]438  ! Calculate particle trajectories
439  !********************************
440
[f13406c]441  if (verbosity.gt.0) then
[8a65cb0]442     if (verbosity.gt.1) then   
443       CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
444       write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
445     endif
446     if (info_flag.eq.1) then
447       print*, 'info only mode (stop)'   
448       stop
449     endif
450     print*,'call timemanager'
[f13406c]451  endif
452
[2eefa58]453  if (verbosity.gt.0) write (*,*) 'timemanager> call wetdepo'
[6ecb30a]454  call timemanager(metdata_format)
[2eefa58]455 
[e200b7a]456
[e9e0f06]457  if (verbosity.gt.0) then
[8a65cb0]458! NIK 16.02.2005
[e9e0f06]459    do i=1,nspec
460      if (tot_inc_count(i).gt.0) then
461         write(*,*) '**********************************************'
462         write(*,*) 'Scavenging statistics for species ', species(i), ':'
463         write(*,*) 'Total number of occurences of below-cloud scavenging', &
464           & tot_blc_count(i)
465         write(*,*) 'Total number of occurences of in-cloud    scavenging', &
466           & tot_inc_count(i)
467         write(*,*) '**********************************************'
468      endif
469    end do
470  endif
[c8fc724]471 
[8a65cb0]472  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
473       &XPART MODEL RUN!'
[e200b7a]474
475end program flexpart
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG