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

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

update flexversion to 10.3

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