Context Navigation

FLEXPART.f90 @ 07c3e71

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format