source: flexpart.git/src/FLEXPART.f90 @ 332fbbd

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file since 332fbbd was 332fbbd, checked in by Ignacio Pisso <ip@…>, 4 years ago

Revert "move license from headers to a different file"
Adding the full GPL license it was apparent that the FSF recommends to
have the statements as file headers

This reverts commit 3481cc180a6b19f71d9d2d04240a91d23f0a2800.

  • 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-11-12)'
83  verbosity=0
84
85  ! Read the pathnames where input/output files are stored
86  !*******************************************************
87
88  inline_options='none'
89  select case (iargc())
90  case (2)
91    call getarg(1,arg1)
92    pathfile=arg1
93    call getarg(2,arg2)
94    inline_options=arg2
95  case (1)
96    call getarg(1,arg1)
97    pathfile=arg1
98    if (arg1(1:1).eq.'-') then
99      write(pathfile,'(a11)') './pathnames'
100      inline_options=arg1
101    endif
102  case (0)
103    write(pathfile,'(a11)') './pathnames'
104  end select
105 
106  ! Print the GPL License statement
107  !*******************************************************
108  print*,'Welcome to FLEXPART ', trim(flexversion)
109  print*,'FLEXPART is free software released under the GNU General Public License.'
110 
111
112  ! Ingest inline options
113  !*******************************************************
114  if (inline_options(1:1).eq.'-') then
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
125    endif
126    !debug mode
127    iopt=index(inline_options,'d')
128    if (iopt.gt.0) then
129      debug_mode=.true.
130    endif
131    if (trim(inline_options).eq.'-i') then
132       print*, 'Info mode: provide detailed run specific information and stop'
133       verbosity=1
134       info_flag=1
135    endif
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
140    endif
141  endif
142           
143  if (verbosity.gt.0) then
144    print*, 'nxmax=',nxmax
145    print*, 'nymax=',nymax
146    print*, 'nzmax=',nzmax
147    print*,'nxshift=',nxshift
148  endif
149 
150  if (verbosity.gt.0) then
151    write(*,*) 'call readpaths'
152  endif
153  call readpaths
154 
155  if (verbosity.gt.1) then !show clock info
156     !print*,'length(4)',length(4)
157     !count=0,count_rate=1000
158     CALL SYSTEM_CLOCK(count_clock0, count_rate, count_max)
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
163  endif
164
165  ! Read the user specifications for the current model run
166  !*******************************************************
167
168  if (verbosity.gt.0) then
169    write(*,*) 'call readcommand'
170  endif
171  call readcommand
172  if (verbosity.gt.0) then
173    write(*,*) '    ldirect=', ldirect
174    write(*,*) '    ibdate,ibtime=',ibdate,ibtime
175    write(*,*) '    iedate,ietime=', iedate,ietime
176    if (verbosity.gt.1) then   
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     
180  endif
181
182  ! Initialize arrays in com_mod
183  !*****************************
184  call com_mod_allocate_part(maxpart)
185
186
187  ! Read the age classes to be used
188  !********************************
189  if (verbosity.gt.0) then
190    write(*,*) 'call readageclasses'
191  endif
192  call readageclasses
193
194  if (verbosity.gt.1) then   
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
197  endif     
198
199  ! Read, which wind fields are available within the modelling period
200  !******************************************************************
201
202  if (verbosity.gt.0) then
203    write(*,*) 'call readavailable'
204  endif 
205  call readavailable
206
207  ! Detect metdata format
208  !**********************
209  if (verbosity.gt.0) then
210    write(*,*) 'call detectformat'
211  endif
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'
221    stop
222  endif
223
224
225
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
234  ! Read the model grid specifications,
235  ! both for the mother domain and eventual nests
236  !**********************************************
237 
238  if (verbosity.gt.0) then
239     write(*,*) 'call gridcheck'
240  endif
241
242  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
243    call gridcheck_ecmwf
244  else
245    call gridcheck_gfs
246  end if
247
248  if (verbosity.gt.1) then   
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
251  endif     
252
253  if (verbosity.gt.0) then
254    write(*,*) 'call gridcheck_nests'
255  endif 
256  call gridcheck_nests
257
258  ! Read the output grid specifications
259  !************************************
260
261  if (verbosity.gt.0) then
262    write(*,*) 'call readoutgrid'
263  endif
264
265  call readoutgrid
266
267  if (nested_output.eq.1) then
268    call readoutgrid_nest
269    if (verbosity.gt.0) then
270      write(*,*) '# readoutgrid_nest'
271    endif
272  endif
273
274  ! Read the receptor points for which extra concentrations are to be calculated
275  !*****************************************************************************
276
277  if (verbosity.eq.1) then
278     print*,'call readreceptors'
279  endif
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
291  if (verbosity.gt.0) then
292    print*,'call readlanduse'
293  endif
294  call readlanduse
295
296  ! Assign fractional cover of landuse classes to each ECMWF grid point
297  !********************************************************************
298
299  if (verbosity.gt.0) then
300    print*,'call assignland'
301  endif
302  call assignland
303
304  ! Read the coordinates of the release locations
305  !**********************************************
306
307  if (verbosity.gt.0) then
308    print*,'call readreleases'
309  endif
310  call readreleases
311
312  ! Read and compute surface resistances to dry deposition of gases
313  !****************************************************************
314
315  if (verbosity.gt.0) then
316    print*,'call readdepo'
317  endif
318  call readdepo
319
320  ! Convert the release point coordinates from geografical to grid coordinates
321  !***************************************************************************
322
323  call coordtrafo 
324  if (verbosity.gt.0) then
325    print*,'call coordtrafo'
326  endif
327
328  ! Initialize all particles to non-existent
329  !*****************************************
330
331  if (verbosity.gt.0) then
332    print*,'Initialize all particles to non-existent'
333  endif
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
342    if (verbosity.gt.0) then
343      print*,'call readpartpositions'
344    endif
345    call readpartpositions
346  else
347    if (verbosity.gt.0) then
348      print*,'numpart=0, numparticlecount=0'
349    endif   
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
358  if (verbosity.gt.0) then
359    print*,'call outgrid_init'
360  endif
361  call outgrid_init
362  if (nested_output.eq.1) call outgrid_init_nest
363
364  ! Read the OH field
365  !******************
366
367  if (OHREA.eqv..TRUE.) then
368    if (verbosity.gt.0) then
369      print*,'call readOHfield'
370    endif
371    call readOHfield
372  endif
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
378#ifdef USE_NCF
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
392#endif
393
394  if (verbosity.gt.0) then
395    print*,'call writeheader'
396  endif
397
398  call writeheader
399  ! FLEXPART 9.2 ticket ?? write header in ASCII format
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
409    print*,'call openreceptors'
410  endif
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
417  if (verbosity.gt.0) then
418    print*,'discretize release times'
419  endif
420  do i=1,numpoint
421    ireleasestart(i)=nint(real(ireleasestart(i))/real(lsynctime))*lsynctime
422    ireleaseend(i)=nint(real(ireleaseend(i))/real(lsynctime))*lsynctime
423  end do
424
425  ! Initialize cloud-base mass fluxes for the convection scheme
426  !************************************************************
427
428  if (verbosity.gt.0) then
429    print*,'Initialize cloud-base mass fluxes for the convection scheme'
430  endif
431
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
445  ! Inform whether output kernel is used or not
446  !*********************************************
447  if (lroot) then
448    if (.not.lusekerneloutput) then
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
456  ! Calculate particle trajectories
457  !********************************
458
459  if (verbosity.gt.0) then
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'
469  endif
470
471  if (verbosity.gt.0) write (*,*) 'timemanager> call wetdepo'
472  call timemanager(metdata_format)
473 
474
475  if (verbosity.gt.0) then
476! NIK 16.02.2005
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
489 
490  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
491       &XPART MODEL RUN!'
492
493end program flexpart
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG