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

GFS_025dev
Last change on this file since a756649 was a756649, checked in by Espen Sollum <eso@…>, 4 years ago

Minor modifications in dev branch

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