source: flexpart.git/src/FLEXPART.f90

release-10
Last change on this file was adead08, checked in by Ignacio Pisso <ip@…>, 5 weeks ago

Merge branch 'master' into release-10 in order to update the master
branch to reflect recent changes.
The updated master branch is described in the GMD paper by Pisso et al.
(2019)

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