source: flexpart.git/src/FLEXPART.f90 @ 3481cc1

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

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