Context Navigation

readcommand.f90 @ 25

Last change on this file since 25 was 24, checked in by igpis, 10 years ago
version 9.2 beta. Changes from HH, AST, MC, NIK, IP. Changes in vert transform. New SPECIES input includes scavenging coefficients
File size: 22.8 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	subroutine readcommand
23
24	!*****************************************************************************
25	! *
26	! This routine reads the user specifications for the current model run. *
27	! *
28	! Author: A. Stohl *
29	! *
30	! 18 May 1996 *
31	! *
32	!*****************************************************************************
33	! *
34	! Variables: *
35	! bdate beginning date as Julian date *
36	! ctl factor by which time step must be smaller than *
37	! Lagrangian time scale *
38	! ibdate,ibtime beginnning date and time (YYYYMMDD, HHMISS) *
39	! ideltas [s] modelling period *
40	! iedate,ietime ending date and time (YYYYMMDD, HHMISS) *
41	! ifine reduction factor for vertical wind time step *
42	! outputforeachrel for forward runs it is possible either to create *
43	! one outputfield or several for each releasepoint *
44	! iflux switch to turn on (1)/off (0) flux calculations *
45	! iout 1 for conc. (residence time for backward runs) output,*
46	! 2 for mixing ratio output, 3 both, 4 for plume *
47	! trajectory output, 5 = options 1 and 4 *
48	! ipin 1 continue simulation with dumped particle data, 0 no *
49	! ipout 0 no particle dump, 1 every output time, 3 only at end*
50	! itsplit [s] time constant for particle splitting *
51	! loutaver [s] concentration output is an average over loutaver *
52	! seconds *
53	! loutsample [s] average is computed from samples taken every [s] *
54	! seconds *
55	! loutstep [s] time interval of concentration output *
56	! lsynctime [s] synchronisation time interval for all particles *
57	! lagespectra switch to turn on (1)/off (0) calculation of age *
58	! spectra *
59	! lconvection value of either 0 and 1 indicating mixing by *
60	! convection *
61	! = 0 .. no convection *
62	! + 1 .. parameterisation of mixing by subgrid-scale *
63	! convection = on *
64	! lsubgrid switch to turn on (1)/off (0) subgrid topography *
65	! parameterization *
66	! method method used to compute the particle pseudovelocities *
67	! mdomainfill 1 use domain-filling option, 0 not, 2 use strat. O3 *
68	! *
69	! Constants: *
70	! unitcommand unit connected to file COMMAND *
71	! *
72	!*****************************************************************************
73
74	use par_mod
75	use com_mod
76
77	implicit none
78
79	real(kind=dp) :: juldate
80	character(len=50) :: line
81	logical :: old
82	logical :: nml_COMMAND=.true. , nmlout=.true. !.false.
83	integer :: readerror
84
85	namelist /command/ &
86	ldirect, &
87	ibdate,ibtime, &
88	iedate,ietime, &
89	loutstep, &
90	loutaver, &
91	loutsample, &
92	itsplit, &
93	lsynctime, &
94	ctl, &
95	ifine, &
96	iout, &
97	ipout, &
98	lsubgrid, &
99	lconvection, &
100	lagespectra, &
101	ipin, &
102	ioutputforeachrelease, &
103	iflux, &
104	mdomainfill, &
105	ind_source, &
106	ind_receptor, &
107	mquasilag, &
108	nested_output, &
109	linit_cond, &
110	surf_only
111
112	! Presetting namelist command
113	ldirect=1
114	ibdate=20000101
115	ibtime=0
116	iedate=20000102
117	ietime=0
118	loutstep=10800
119	loutaver=10800
120	loutsample=900
121	itsplit=999999999
122	lsynctime=900
123	ctl=-5.0
124	ifine=4
125	iout=3
126	ipout=0
127	lsubgrid=1
128	lconvection=1
129	lagespectra=0
130	ipin=1
131	ioutputforeachrelease=1
132	iflux=1
133	mdomainfill=0
134	ind_source=1
135	ind_receptor=1
136	mquasilag=0
137	nested_output=0
138	linit_cond=0
139	surf_only=0
140
141	! Open the command file and read user options
142	! Namelist input first: try to read as namelist file
143	!**************************************************************************
144	open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
145	form='formatted',iostat=readerror)
146	! If fail, check if file does not exist
147	if (readerror.ne.0) then
148	print,'**ERROR: file COMMAND not found in '
149	print*, path(1)(1:length(1))//'COMMAND'
150	print*, 'Check your pathnames file.'
151	stop
152	endif
153
154	! print error code
155	!write(,) 'readcommand > readerror open=' , readerror
156	!probe first line
157	read (unitcommand,901) line
158	!write(,) 'index(line,COMMAND) =', index(line,'COMMAND')
159
160
161	!default is namelist input
162	! distinguish namelist from fixed text input
163	if (index(line,'COMMAND') .eq. 0) then
164	nml_COMMAND = .false.
165	!write(,) 'COMMAND file does not contain the string COMMAND in the first line'
166	endif
167	!write(,) 'readcommand > read as namelist? ' , nml_COMMAND
168	rewind(unitcommand)
169	read(unitcommand,command,iostat=readerror)
170
171	close(unitcommand)
172
173	!write(,) 'readcommand > readerror read=' , readerror
174	! If error in namelist format, try to open with old input code
175	! if (readerror.ne.0) then
176	! IP 21/5/2014 the previous line cause the old long format
177	! to be confused with namelist input
178
179	! use text input
180	if (nml_COMMAND .eqv. .false.) then
181
182	open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
183	err=999)
184
185	! Check the format of the COMMAND file (either in free format,
186	! or using formatted mask)
187	! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
188	!**************************************************************************
189
190	call skplin(9,unitcommand)
191	read (unitcommand,901) line
192	901 format (a)
193	if (index(line,'LDIRECT') .eq. 0) then
194	old = .false.
195	!write(,) 'readcommand old short'
196	else
197	old = .true.
198	!write(,) 'readcommand old long'
199	endif
200	rewind(unitcommand)
201
202
203	! Read parameters
204	!****************
205
206	call skplin(7,unitcommand)
207	if (old) call skplin(1,unitcommand)
208
209	read(unitcommand,*) ldirect
210	if (old) call skplin(3,unitcommand)
211	read(unitcommand,*) ibdate,ibtime
212	if (old) call skplin(3,unitcommand)
213	read(unitcommand,*) iedate,ietime
214	if (old) call skplin(3,unitcommand)
215	read(unitcommand,*) loutstep
216	if (old) call skplin(3,unitcommand)
217	read(unitcommand,*) loutaver
218	if (old) call skplin(3,unitcommand)
219	read(unitcommand,*) loutsample
220	if (old) call skplin(3,unitcommand)
221	read(unitcommand,*) itsplit
222	if (old) call skplin(3,unitcommand)
223	read(unitcommand,*) lsynctime
224	if (old) call skplin(3,unitcommand)
225	read(unitcommand,*) ctl
226	if (old) call skplin(3,unitcommand)
227	read(unitcommand,*) ifine
228	if (old) call skplin(3,unitcommand)
229	read(unitcommand,*) iout
230	if (old) call skplin(3,unitcommand)
231	read(unitcommand,*) ipout
232	if (old) call skplin(3,unitcommand)
233	read(unitcommand,*) lsubgrid
234	if (old) call skplin(3,unitcommand)
235	read(unitcommand,*) lconvection
236	if (old) call skplin(3,unitcommand)
237	read(unitcommand,*) lagespectra
238	if (old) call skplin(3,unitcommand)
239	read(unitcommand,*) ipin
240	if (old) call skplin(3,unitcommand)
241	read(unitcommand,*) ioutputforeachrelease
242	if (old) call skplin(3,unitcommand)
243	read(unitcommand,*) iflux
244	if (old) call skplin(3,unitcommand)
245	read(unitcommand,*) mdomainfill
246	if (old) call skplin(3,unitcommand)
247	read(unitcommand,*) ind_source
248	if (old) call skplin(3,unitcommand)
249	read(unitcommand,*) ind_receptor
250	if (old) call skplin(3,unitcommand)
251	read(unitcommand,*) mquasilag
252	if (old) call skplin(3,unitcommand)
253	read(unitcommand,*) nested_output
254	if (old) call skplin(3,unitcommand)
255	read(unitcommand,*) linit_cond
256	if (old) call skplin(3,unitcommand)
257	read(unitcommand,*) surf_only
258	close(unitcommand)
259
260	endif ! input format
261
262	! write command file in namelist format to output directory if requested
263	if (nmlout.eqv..true.) then
264	open(unitcommand,file=path(2)(1:length(2))//'COMMAND.namelist',err=1000)
265	write(unitcommand,nml=command)
266	close(unitcommand)
267	! open(unitheader,file=path(2)(1:length(2))//'header_nml',status='new',err=999)
268	! write(unitheader,NML=COMMAND)
269	!close(unitheader)
270	endif
271
272	ifine=max(ifine,1)
273
274	! Determine how Markov chain is formulated (for w or for w/sigw)
275	!***************************************************************
276
277	if (ctl.ge.0.1) then
278	turbswitch=.true.
279	else
280	turbswitch=.false.
281	ifine=1
282	endif
283	fine=1./real(ifine)
284	ctl=1./ctl
285
286	! Set the switches required for the various options for input/output units
287	!*************************************************************************
288	!AF Set the switches IND_REL and IND_SAMP for the release and sampling
289	!Af switches for the releasefile:
290	!Af IND_REL = 1 : xmass * rho
291	!Af IND_REL = 0 : xmass * 1
292
293	!Af switches for the conccalcfile:
294	!AF IND_SAMP = 0 : xmass * 1
295	!Af IND_SAMP = -1 : xmass / rho
296
297	!AF IND_SOURCE switches between different units for concentrations at the source
298	!Af NOTE that in backward simulations the release of computational particles
299	!Af takes place at the "receptor" and the sampling of p[articles at the "source".
300	!Af 1 = mass units
301	!Af 2 = mass mixing ratio units
302	!Af IND_RECEPTOR switches between different units for concentrations at the receptor
303	!Af 1 = mass units
304	!Af 2 = mass mixing ratio units
305
306	if ( ldirect .eq. 1 ) then ! FWD-Run
307	!Af set release-switch
308	if (ind_source .eq. 1 ) then !mass
309	ind_rel = 0
310	else ! mass mix
311	ind_rel = 1
312	endif
313	!Af set sampling switch
314	if (ind_receptor .eq. 1) then !mass
315	ind_samp = 0
316	else ! mass mix
317	ind_samp = -1
318	endif
319	elseif (ldirect .eq. -1 ) then !BWD-Run
320	!Af set sampling switch
321	if (ind_source .eq. 1 ) then !mass
322	ind_samp = -1
323	else ! mass mix
324	ind_samp = 0
325	endif
326	!Af set release-switch
327	if (ind_receptor .eq. 1) then !mass
328	ind_rel = 1
329	else ! mass mix
330	ind_rel = 0
331	endif
332	endif
333
334	!*************************************************************
335	! Check whether valid options have been chosen in file COMMAND
336	!*************************************************************
337
338	! Check options for initial condition output: Switch off for forward runs
339	!************************************************************************
340
341	if (ldirect.eq.1) linit_cond=0
342	if ((linit_cond.lt.0).or.(linit_cond.gt.2)) then
343	write(,) ' #### FLEXPART MODEL ERROR! INVALID OPTION #### '
344	write(,) ' #### FOR LINIT_COND IN FILE "COMMAND". #### '
345	stop
346	endif
347
348	! Check input dates
349	!******************
350
351	if (iedate.lt.ibdate) then
352	write(,) ' #### FLEXPART MODEL ERROR! BEGINNING DATE #### '
353	write(,) ' #### IS LARGER THAN ENDING DATE. CHANGE #### '
354	write(,) ' #### EITHER POINT 2 OR POINT 3 IN FILE #### '
355	write(,) ' #### "COMMAND". #### '
356	stop
357	else if (iedate.eq.ibdate) then
358	if (ietime.lt.ibtime) then
359	write(,) ' #### FLEXPART MODEL ERROR! BEGINNING TIME #### '
360	write(,) ' #### IS LARGER THAN ENDING TIME. CHANGE #### '
361	write(,) ' #### EITHER POINT 2 OR POINT 3 IN FILE #### '
362	write(,) ' #### "COMMAND". #### '
363	stop
364	endif
365	endif
366
367
368	! Determine kind of dispersion method
369	!************************************
370
371	if (ctl.gt.0.) then
372	method=1
373	mintime=minstep
374	else
375	method=0
376	mintime=lsynctime
377	endif
378
379	! Check whether a valid option for gridded model output has been chosen
380	!**********************************************************************
381
382	if ((iout.lt.1).or.(iout.gt.5)) then
383	write(,) ' #### FLEXPART MODEL ERROR! FILE COMMAND: #### '
384	write(,) ' #### IOUT MUST BE 1, 2, 3, 4, OR 5! #### '
385	stop
386	endif
387
388	!AF check consistency between units and volume mixing ratio
389	if ( ((iout.eq.2).or.(iout.eq.3)).and. &
390	(ind_source.gt.1 .or.ind_receptor.gt.1) ) then
391	write(,) ' #### FLEXPART MODEL ERROR! FILE COMMAND: #### '
392	write(,) ' #### VOLUME MIXING RATIO ONLY SUPPORTED #### '
393	write(,) ' #### FOR MASS UNITS (at the moment) #### '
394	stop
395	endif
396
397
398
399	! For quasilag output for each release is forbidden
400	!*****************************************************************************
401
402	if ((ioutputforeachrelease.eq.1).and.(mquasilag.eq.1)) then
403	write(,) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####'
404	write(,) '#### OUTPUTFOREACHRELEASE AND QUASILAGRANGIAN####'
405	write(,) '#### MODE IS NOT POSSIBLE ! ####'
406	stop
407	endif
408
409
410	! For quasilag backward is forbidden
411	!*****************************************************************************
412
413	if ((ldirect.lt.0).and.(mquasilag.eq.1)) then
414	write(,) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####'
415	write(,) '#### FOR BACKWARD RUNS, QUASILAGRANGIAN MODE ####'
416	write(,) '#### IS NOT POSSIBLE ! ####'
417	stop
418	endif
419
420
421	! For backward runs one releasefield for all releases makes no sense,
422	! For quasilag and domainfill ioutputforechrelease is forbidden
423	!*****************************************************************************
424
425	if ((ldirect.lt.0).and.(ioutputforeachrelease.eq.0)) then
426	write(,) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####'
427	write(,) '#### FOR BACKWARD RUNS, IOUTPUTFOREACHRLEASE ####'
428	write(,) '#### MUST BE SET TO ONE! ####'
429	stop
430	endif
431
432
433	! For backward runs one releasefield for all releases makes no sense,
434	! and is "forbidden"
435	!*****************************************************************************
436
437	if ((mdomainfill.eq.1).and.(ioutputforeachrelease.eq.1)) then
438	write(,) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####'
439	write(,) '#### FOR DOMAIN FILLING RUNS OUTPUT FOR ####'
440	write(,) '#### EACH RELEASE IS FORBIDDEN ! ####'
441	stop
442	endif
443
444
445	! For domain-filling trajectories, a plume centroid trajectory makes no sense,
446	! For backward runs, only residence time output (iout=1) or plume trajectories (iout=4),
447	! or both (iout=5) makes sense; other output options are "forbidden"
448	!*****************************************************************************
449
450	if (ldirect.lt.0) then
451	if ((iout.eq.2).or.(iout.eq.3)) then
452	write(,) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####'
453	write(,) '#### FOR BACKWARD RUNS, IOUT MUST BE 1,4,OR 5####'
454	stop
455	endif
456	endif
457
458
459	! For domain-filling trajectories, a plume centroid trajectory makes no sense,
460	! and is "forbidden"
461	!*****************************************************************************
462
463	if (mdomainfill.ge.1) then
464	if ((iout.eq.4).or.(iout.eq.5)) then
465	write(,) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####'
466	write(,) '#### FOR DOMAIN-FILLING TRAJECTORY OPTION, ####'
467	write(,) '#### IOUT MUST NOT BE SET TO 4 OR 5. ####'
468	stop
469	endif
470	endif
471
472
473
474	! Check whether a valid options for particle dump has been chosen
475	!****************************************************************
476
477	if ((ipout.ne.0).and.(ipout.ne.1).and.(ipout.ne.2)) then
478	write(,) ' #### FLEXPART MODEL ERROR! FILE COMMAND: #### '
479	write(,) ' #### IPOUT MUST BE 1, 2 OR 3! #### '
480	stop
481	endif
482
483	if(lsubgrid.ne.1.and.verbosity.eq.0) then
484	write(,) ' ---------------- '
485	write(,) ' INFORMATION: SUBGRIDSCALE TERRAIN EFFECT IS'
486	write(,) ' NOT PARAMETERIZED DURING THIS SIMULATION. '
487	write(,) ' ---------------- '
488	endif
489
490
491	! Check whether convection scheme is either turned on or off
492	!***********************************************************
493
494	if ((lconvection.ne.0).and.(lconvection.ne.1)) then
495	write(,) ' #### FLEXPART MODEL ERROR! FILE COMMAND: #### '
496	write(,) ' #### LCONVECTION MUST BE SET TO EITHER 1 OR 0#### '
497	stop
498	endif
499
500
501	! Check whether synchronisation interval is sufficiently short
502	!*************************************************************
503
504	if (lsynctime.gt.(idiffnorm/2)) then
505	write(,) ' #### FLEXPART MODEL ERROR! SYNCHRONISATION #### '
506	write(,) ' #### TIME IS TOO LONG. MAKE IT SHORTER. #### '
507	write(,) ' #### MINIMUM HAS TO BE: ', idiffnorm/2
508	stop
509	endif
510
511
512	! Check consistency of the intervals, sampling periods, etc., for model output
513	!*****************************************************************************
514
515	if (loutaver.eq.0) then
516	write(,) ' #### FLEXPART MODEL ERROR! TIME AVERAGE OF #### '
517	write(,) ' #### CONCENTRATION FIELD OUTPUT MUST NOT BE #### '
518	write(,) ' #### ZERO. #### '
519	write(,) ' #### CHANGE INPUT IN FILE COMMAND. #### '
520	stop
521	endif
522
523	if (loutaver.gt.loutstep) then
524	write(,) ' #### FLEXPART MODEL ERROR! TIME AVERAGE OF #### '
525	write(,) ' #### CONCENTRATION FIELD OUTPUT MUST NOT BE #### '
526	write(,) ' #### GREATER THAN INTERVAL OF OUTPUT. #### '
527	write(,) ' #### CHANGE INPUT IN FILE COMMAND. #### '
528	stop
529	endif
530
531	if (loutsample.gt.loutaver) then
532	write(,) ' #### FLEXPART MODEL ERROR! SAMPLING TIME OF #### '
533	write(,) ' #### CONCENTRATION FIELD OUTPUT MUST NOT BE #### '
534	write(,) ' #### GREATER THAN TIME AVERAGE OF OUTPUT. #### '
535	write(,) ' #### CHANGE INPUT IN FILE COMMAND. #### '
536	stop
537	endif
538
539	if (mod(loutaver,lsynctime).ne.0) then
540	write(,) ' #### FLEXPART MODEL ERROR! AVERAGING TIME OF #### '
541	write(,) ' #### CONCENTRATION FIELD MUST BE A MULTIPLE #### '
542	write(,) ' #### OF THE SYNCHRONISATION INTERVAL #### '
543	stop
544	endif
545
546	if ((loutaver/lsynctime).lt.2) then
547	write(,) ' #### FLEXPART MODEL ERROR! AVERAGING TIME OF #### '
548	write(,) ' #### CONCENTRATION FIELD MUST BE AT LEAST #### '
549	write(,) ' #### TWICE THE SYNCHRONISATION INTERVAL #### '
550	stop
551	endif
552
553	if (mod(loutstep,lsynctime).ne.0) then
554	write(,) ' #### FLEXPART MODEL ERROR! INTERVAL BETWEEN #### '
555	write(,) ' #### CONCENTRATION FIELDS MUST BE A MULTIPLE #### '
556	write(,) ' #### OF THE SYNCHRONISATION INTERVAL #### '
557	stop
558	endif
559
560	if ((loutstep/lsynctime).lt.2) then
561	write(,) ' #### FLEXPART MODEL ERROR! INTERVAL BETWEEN #### '
562	write(,) ' #### CONCENTRATION FIELDS MUST BE AT LEAST #### '
563	write(,) ' #### TWICE THE SYNCHRONISATION INTERVAL #### '
564	stop
565	endif
566
567	if (mod(loutsample,lsynctime).ne.0) then
568	write(,) ' #### FLEXPART MODEL ERROR! SAMPLING TIME OF #### '
569	write(,) ' #### CONCENTRATION FIELD MUST BE A MULTIPLE #### '
570	write(,) ' #### OF THE SYNCHRONISATION INTERVAL #### '
571	stop
572	endif
573
574	if (itsplit.lt.loutaver) then
575	write(,) ' #### FLEXPART MODEL ERROR! SPLITTING TIME FOR#### '
576	write(,) ' #### PARTICLES IS TOO SHORT. PLEASE INCREASE #### '
577	write(,) ' #### SPLITTING TIME CONSTANT. #### '
578	stop
579	endif
580
581	if ((mquasilag.eq.1).and.(iout.ge.4)) then
582	write(,) ' #### FLEXPART MODEL ERROR! CONFLICTING #### '
583	write(,) ' #### OPTIONS: IF MQUASILAG=1, PLUME #### '
584	write(,) ' #### TRAJECTORY OUTPUT IS IMPOSSIBLE. #### '
585	stop
586	endif
587
588	! Compute modeling time in seconds and beginning date in Julian date
589	!*******************************************************************
590
591	outstep=real(abs(loutstep))
592	if (ldirect.eq.1) then
593	bdate=juldate(ibdate,ibtime)
594	edate=juldate(iedate,ietime)
595	ideltas=nint((edate-bdate)*86400.)
596	else if (ldirect.eq.-1) then
597	loutaver=-1*loutaver
598	loutstep=-1*loutstep
599	loutsample=-1*loutsample
600	lsynctime=-1*lsynctime
601	bdate=juldate(iedate,ietime)
602	edate=juldate(ibdate,ibtime)
603	ideltas=nint((edate-bdate)*86400.)
604	else
605	write(,) ' #### FLEXPART MODEL ERROR! DIRECTION IN #### '
606	write(,) ' #### FILE "COMMAND" MUST BE EITHER -1 OR 1. #### '
607	stop
608	endif
609
610	return
611
612	999 write(,) ' #### FLEXPART MODEL ERROR! FILE "COMMAND" #### '
613	write(,) ' #### CANNOT BE OPENED IN THE DIRECTORY #### '
614	write(*,'(a)') path(1)(1:length(1))
615	stop
616
617	1000 write(,) ' #### FLEXPART MODEL ERROR! FILE "COMMAND" #### '
618	write(,) ' #### CANNOT BE OPENED IN THE DIRECTORY #### '
619	write(*,'(a)') path(2)(1:length(1))
620	stop
621	end subroutine readcommand

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

Download in other formats:

Original Format