source: flexpart.git/src/FLEXPART.f90 @ c8fc724

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since c8fc724 was c8fc724, checked in by Espen Sollum ATMOS <eso@…>, 7 years ago

Bugfix: using namelist species files would sometimes set incorrect below-cloud scavenging parameters when using more than 1 species. Update: per-species scavenging statistics are written out at end of simulation.

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