Context Navigation

readcommand.f90 @ 5f2e8f6

flexpart-noresm

Last change on this file since 5f2e8f6 was 5f2e8f6, checked in by Ignacio Pisso <Ignacio.Pisso@…>, 8 years ago
new flexpart noresm code in src
Property mode set to `100755`
File size: 21.3 KB

Line
1	!**********************************************************************
2	! Copyright 2016 *
3	! Andreas Stohl, Massimo Cassiani, Petra Seibert, A. Frank, *
4	! Gerhard Wotawa, Caroline Forster, Sabine Eckhardt, John Burkhart, *
5	! Harald Sodemann, Ignacio Pisso *
6	! *
7	! This file is part of FLEXPART-NorESM *
8	! *
9	! FLEXPART-NorESM is free software: you can redistribute it *
10	! and/or modify *
11	! it under the terms of the GNU General Public License as published by*
12	! the Free Software Foundation, either version 3 of the License, or *
13	! (at your option) any later version. *
14	! *
15	! FLEXPART-NorESM is distributed in the hope that it will be useful, *
16	! but WITHOUT ANY WARRANTY; without even the implied warranty of *
17	! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
18	! GNU General Public License for more details. *
19	! *
20	! You should have received a copy of the GNU General Public License *
21	! along with FLEXPART-NorESM. *
22	! If not, see <http://www.gnu.org/licenses/>. *
23	!**********************************************************************
24
25	subroutine readcommand
26
27	!*****************************************************************************
28	! *
29	! This routine reads the user specifications for the current model run. *
30	! *
31	! Author: A. Stohl *
32	! *
33	! 18 May 1996 *
34	! *
35	! modified by M. Cassiani, 2016 *
36	! added flag for alterntive methods to compute vertical velocity *
37	! *
38	!*****************************************************************************
39	! *
40	! Variables: *
41	! bdate beginning date as Julian date *
42	! ctl factor by which time step must be smaller than *
43	! Lagrangian time scale *
44	! ibdate,ibtime beginnning date and time (YYYYMMDD, HHMISS) *
45	! ideltas [s] modelling period *
46	! iedate,ietime ending date and time (YYYYMMDD, HHMISS) *
47	! ifine reduction factor for vertical wind time step *
48	! outputforeachrel for forward runs it is possible either to create *
49	! one outputfield or several for each releasepoint *
50	! iflux switch to turn on (1)/off (0) flux calculations *
51	! iout 1 for conc. (residence time for backward runs) output,*
52	! 2 for mixing ratio output, 3 both, 4 for plume *
53	! trajectory output, 5 = options 1 and 4 *
54	! ipin 1 continue simulation with dumped particle data, 0 no *
55	! ipout 0 no particle dump, 1 every output time, 3 only at end*
56	! itsplit [s] time constant for particle splitting *
57	! loutaver [s] concentration output is an average over loutaver *
58	! seconds *
59	! loutsample [s] average is computed from samples taken every [s] *
60	! seconds *
61	! loutstep [s] time interval of concentration output *
62	! lsynctime [s] synchronisation time interval for all particles *
63	! lagespectra switch to turn on (1)/off (0) calculation of age *
64	! spectra *
65	! lconvection value of either 0 and 1 indicating mixing by *
66	! convection *
67	! = 0 .. no convection *
68	! + 1 .. parameterisation of mixing by subgrid-scale *
69	! convection = on *
70	! lsubgrid switch to turn on (1)/off (0) subgrid topography *
71	! parameterization *
72	! method method used to compute the particle pseudovelocities *
73	! mdomainfill 1 use domain-filling option, 0 not, 2 use strat. O3 *
74	! *
75	! Constants: *
76	! unitcommand unit connected to file COMMAND *
77	! *
78	!*****************************************************************************
79
80	use par_mod
81	use com_mod
82
83	implicit none
84
85	real(kind=dp) :: juldate
86	character(len=50) :: line
87	logical :: old
88
89
90	! Open the command file and read user options
91	!********************************************
92
93
94	open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
95	err=999)
96
97	! Check the format of the COMMAND file (either in free format,
98	! or using formatted mask)
99	! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
100	!**************************************************************************
101
102	call skplin(9,unitcommand)
103	read (unitcommand,901) line
104	901 format (a)
105	if (index(line,'LDIRECT') .eq. 0) then
106	old = .false.
107	else
108	old = .true.
109	endif
110	rewind(unitcommand)
111
112	! Read parameters
113	!****************
114
115	call skplin(7,unitcommand)
116	if (old) call skplin(1,unitcommand)
117
118	read(unitcommand,*) ldirect
119	if (old) call skplin(3,unitcommand)
120	read(unitcommand,*) ibdate,ibtime
121	if (old) call skplin(3,unitcommand)
122	read(unitcommand,*) iedate,ietime
123	if (old) call skplin(3,unitcommand)
124	read(unitcommand,*) loutstep
125	if (old) call skplin(3,unitcommand)
126	read(unitcommand,*) loutaver
127	if (old) call skplin(3,unitcommand)
128	read(unitcommand,*) loutsample
129	if (old) call skplin(3,unitcommand)
130	read(unitcommand,*) itsplit
131	if (old) call skplin(3,unitcommand)
132	read(unitcommand,*) lsynctime
133	if (old) call skplin(3,unitcommand)
134	read(unitcommand,*) ctl
135	if (old) call skplin(3,unitcommand)
136	read(unitcommand,*) ifine
137	if (old) call skplin(3,unitcommand)
138	read(unitcommand,*) iout
139	if (old) call skplin(3,unitcommand)
140	read(unitcommand,*) ipout
141	if (old) call skplin(3,unitcommand)
142	read(unitcommand,*) lsubgrid
143	if (old) call skplin(3,unitcommand)
144	read(unitcommand,*) lconvection
145	if (old) call skplin(3,unitcommand)
146	read(unitcommand,*) lagespectra
147	if (old) call skplin(3,unitcommand)
148	read(unitcommand,*) ipin
149	if (old) call skplin(3,unitcommand)
150	read(unitcommand,*) ioutputforeachrelease
151	if (old) call skplin(3,unitcommand)
152	read(unitcommand,*) iflux
153	if (old) call skplin(3,unitcommand)
154	read(unitcommand,*) mdomainfill
155	if (old) call skplin(3,unitcommand)
156	read(unitcommand,*) ind_source
157	if (old) call skplin(3,unitcommand)
158	read(unitcommand,*) ind_receptor
159	if (old) call skplin(3,unitcommand)
160	read(unitcommand,*) mquasilag
161	if (old) call skplin(3,unitcommand)
162	read(unitcommand,*) nested_output
163	if (old) call skplin(3,unitcommand)
164	read(unitcommand,*) linit_cond
165	if (old) call skplin(3,unitcommand) !added by mc
166	read(unitcommand,*) method_w
167	close(unitcommand)
168
169
170	close(unitcommand)
171
172	ifine=max(ifine,1)
173
174
175	! Determine how Markov chain is formulated (for w or for w/sigw)
176	!***************************************************************
177
178	if (ctl.ge.0.1) then
179	turbswitch=.true.
180	else
181	turbswitch=.false.
182	ifine=1
183	endif
184	fine=1./real(ifine)
185	ctl=1./ctl
186
187	! Set the switches required for the various options for input/output units
188	!*************************************************************************
189	!AF Set the switches IND_REL and IND_SAMP for the release and sampling
190	!Af switches for the releasefile:
191	!Af IND_REL = 1 : xmass * rho
192	!Af IND_REL = 0 : xmass * 1
193
194	!Af switches for the conccalcfile:
195	!AF IND_SAMP = 0 : xmass * 1
196	!Af IND_SAMP = -1 : xmass / rho
197
198	!AF IND_SOURCE switches between different units for concentrations at the source
199	!Af NOTE that in backward simulations the release of computational particles
200	!Af takes place at the "receptor" and the sampling of p[articles at the "source".
201	!Af 1 = mass units
202	!Af 2 = mass mixing ratio units
203	!Af IND_RECEPTOR switches between different units for concentrations at the receptor
204	!Af 1 = mass units
205	!Af 2 = mass mixing ratio units
206
207	if ( ldirect .eq. 1 ) then ! FWD-Run
208	!Af set release-switch
209	if (ind_source .eq. 1 ) then !mass
210	ind_rel = 0
211	else ! mass mix
212	ind_rel = 1
213	endif
214	!Af set sampling switch
215	if (ind_receptor .eq. 1) then !mass
216	ind_samp = 0
217	else ! mass mix
218	ind_samp = -1
219	endif
220	elseif (ldirect .eq. -1 ) then !BWD-Run
221	!Af set sampling switch
222	if (ind_source .eq. 1 ) then !mass
223	ind_samp = -1
224	else ! mass mix
225	ind_samp = 0
226	endif
227	!Af set release-switch
228	if (ind_receptor .eq. 1) then !mass
229	ind_rel = 1
230	else ! mass mix
231	ind_rel = 0
232	endif
233	endif
234
235	!*************************************************************
236	! Check whether valid options have been chosen in file COMMAND
237	!*************************************************************
238
239	! Check options for initial condition output: Switch off for forward runs
240	!************************************************************************
241
242	if (ldirect.eq.1) linit_cond=0
243	if ((linit_cond.lt.0).or.(linit_cond.gt.2)) then
244	write(,) ' #### FLEXPART MODEL ERROR! INVALID OPTION #### '
245	write(,) ' #### FOR LINIT_COND IN FILE "COMMAND". #### '
246	stop
247	endif
248
249	! Check input dates
250	!******************
251
252	if (iedate.lt.ibdate) then
253	write(,) ' #### FLEXPART MODEL ERROR! BEGINNING DATE #### '
254	write(,) ' #### IS LARGER THAN ENDING DATE. CHANGE #### '
255	write(,) ' #### EITHER POINT 2 OR POINT 3 IN FILE #### '
256	write(,) ' #### "COMMAND". #### '
257	stop
258	else if (iedate.eq.ibdate) then
259	if (ietime.lt.ibtime) then
260	write(,) ' #### FLEXPART MODEL ERROR! BEGINNING TIME #### '
261	write(,) ' #### IS LARGER THAN ENDING TIME. CHANGE #### '
262	write(,) ' #### EITHER POINT 2 OR POINT 3 IN FILE #### '
263	write(,) ' #### "COMMAND". #### '
264	stop
265	endif
266	else if (bdate.lt.yearorigin) then !modified by mc
267	write(,) ' #### FLEXPART MODEL ERROR! BEGINNING TIME #### ' !modified by mc
268	write(,) ' #### must be larger than yearzero: Difference in time ####' !modified by mc
269	write(,) ' #### computed using no-leap calendar from an ####' !modified by mc
270	write(,) ' #### arbitrary starting yearorigin #### ' !modified by mc
271	write(,) ' #### "COMMAND". #### ' !modified by mc
272	stop
273	endif
274
275	! Determine kind of dispersion method
276	!************************************
277
278	if (ctl.gt.0.) then
279	method=1
280	mintime=minstep
281	else
282	method=0
283	mintime=lsynctime
284	endif
285
286	! Check whether a valid option for gridded model output has been chosen
287	!**********************************************************************
288
289	if ((iout.lt.1).or.(iout.gt.5)) then
290	write(,) ' #### FLEXPART MODEL ERROR! FILE COMMAND: #### '
291	write(,) ' #### IOUT MUST BE 1, 2, 3, 4, OR 5! #### '
292	stop
293	endif
294
295	!AF check consistency between units and volume mixing ratio
296	if ( ((iout.eq.2).or.(iout.eq.3)).and. &
297	(ind_source.gt.1 .or.ind_receptor.gt.1) ) then
298	write(,) ' #### FLEXPART MODEL ERROR! FILE COMMAND: #### '
299	write(,) ' #### VOLUME MIXING RATIO ONLY SUPPORTED #### '
300	write(,) ' #### FOR MASS UNITS (at the moment) #### '
301	stop
302	endif
303
304
305
306	! For quasilag output for each release is forbidden
307	!*****************************************************************************
308
309	if ((ioutputforeachrelease.eq.1).and.(mquasilag.eq.1)) then
310	write(,) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####'
311	write(,) '#### OUTPUTFOREACHRELEASE AND QUASILAGRANGIAN####'
312	write(,) '#### MODE IS NOT POSSIBLE ! ####'
313	stop
314	endif
315
316
317	! For quasilag backward is forbidden
318	!*****************************************************************************
319
320	if ((ldirect.lt.0).and.(mquasilag.eq.1)) then
321	write(,) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####'
322	write(,) '#### FOR BACKWARD RUNS, QUASILAGRANGIAN MODE ####'
323	write(,) '#### IS NOT POSSIBLE ! ####'
324	stop
325	endif
326
327
328	! For backward runs one releasefield for all releases makes no sense,
329	! For quasilag and domainfill ioutputforechrelease is forbidden
330	!*****************************************************************************
331
332	if ((ldirect.lt.0).and.(ioutputforeachrelease.eq.0)) then
333	write(,) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####'
334	write(,) '#### FOR BACKWARD RUNS, IOUTPUTFOREACHRLEASE ####'
335	write(,) '#### MUST BE SET TO ONE! ####'
336	stop
337	endif
338
339
340	! For backward runs one releasefield for all releases makes no sense,
341	! and is "forbidden"
342	!*****************************************************************************
343
344	if ((mdomainfill.eq.1).and.(ioutputforeachrelease.eq.1)) then
345	write(,) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####'
346	write(,) '#### FOR DOMAIN FILLING RUNS OUTPUT FOR ####'
347	write(,) '#### EACH RELEASE IS FORBIDDEN ! ####'
348	stop
349	endif
350
351
352	! For domain-filling trajectories, a plume centroid trajectory makes no sense,
353	! For backward runs, only residence time output (iout=1) or plume trajectories (iout=4),
354	! or both (iout=5) makes sense; other output options are "forbidden"
355	!*****************************************************************************
356
357	if (ldirect.lt.0) then
358	if ((iout.eq.2).or.(iout.eq.3)) then
359	write(,) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####'
360	write(,) '#### FOR BACKWARD RUNS, IOUT MUST BE 1,4,OR 5####'
361	stop
362	endif
363	endif
364
365
366	! For domain-filling trajectories, a plume centroid trajectory makes no sense,
367	! and is "forbidden"
368	!*****************************************************************************
369
370	if (mdomainfill.ge.1) then
371	if ((iout.eq.4).or.(iout.eq.5)) then
372	write(,) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####'
373	write(,) '#### FOR DOMAIN-FILLING TRAJECTORY OPTION, ####'
374	write(,) '#### IOUT MUST NOT BE SET TO 4 OR 5. ####'
375	stop
376	endif
377	endif
378
379
380
381	! Check whether a valid options for particle dump has been chosen
382	!****************************************************************
383
384	if ((ipout.ne.0).and.(ipout.ne.1).and.(ipout.ne.2)) then
385	write(,) ' #### FLEXPART MODEL ERROR! FILE COMMAND: #### '
386	write(,) ' #### IPOUT MUST BE 1, 2 OR 3! #### '
387	stop
388	endif
389
390	if(lsubgrid.ne.1) then
391	write(,) ' ---------------- '
392	write(,) ' INFORMATION: SUBGRIDSCALE TERRAIN EFFECT IS'
393	write(,) ' NOT PARAMETERIZED DURING THIS SIMULATION. '
394	write(,) ' ---------------- '
395	endif
396
397
398	! Check whether convection scheme is either turned on or off
399	!***********************************************************
400
401	if ((lconvection.ne.0).and.(lconvection.ne.1)) then
402	write(,) ' #### FLEXPART MODEL ERROR! FILE COMMAND: #### '
403	write(,) ' #### LCONVECTION MUST BE SET TO EITHER 1 OR 0#### '
404	stop
405	endif
406
407
408	! Check whether synchronisation interval is sufficiently short
409	!*************************************************************
410
411	if (lsynctime.gt.(idiffnorm/2)) then
412	write(,) ' #### FLEXPART MODEL ERROR! SYNCHRONISATION #### '
413	write(,) ' #### TIME IS TOO LONG. MAKE IT SHORTER. #### '
414	write(,) ' #### MINIMUM HAS TO BE: ', idiffnorm/2
415	stop
416	endif
417
418
419	! Check consistency of the intervals, sampling periods, etc., for model output
420	!*****************************************************************************
421
422	if (loutaver.eq.0) then
423	write(,) ' #### FLEXPART MODEL ERROR! TIME AVERAGE OF #### '
424	write(,) ' #### CONCENTRATION FIELD OUTPUT MUST NOT BE #### '
425	write(,) ' #### ZERO. #### '
426	write(,) ' #### CHANGE INPUT IN FILE COMMAND. #### '
427	stop
428	endif
429
430	if (loutaver.gt.loutstep) then
431	write(,) ' #### FLEXPART MODEL ERROR! TIME AVERAGE OF #### '
432	write(,) ' #### CONCENTRATION FIELD OUTPUT MUST NOT BE #### '
433	write(,) ' #### GREATER THAN INTERVAL OF OUTPUT. #### '
434	write(,) ' #### CHANGE INPUT IN FILE COMMAND. #### '
435	stop
436	endif
437
438	if (loutsample.gt.loutaver) then
439	write(,) ' #### FLEXPART MODEL ERROR! SAMPLING TIME OF #### '
440	write(,) ' #### CONCENTRATION FIELD OUTPUT MUST NOT BE #### '
441	write(,) ' #### GREATER THAN TIME AVERAGE OF OUTPUT. #### '
442	write(,) ' #### CHANGE INPUT IN FILE COMMAND. #### '
443	stop
444	endif
445
446	if (mod(loutaver,lsynctime).ne.0) then
447	write(,) ' #### FLEXPART MODEL ERROR! AVERAGING TIME OF #### '
448	write(,) ' #### CONCENTRATION FIELD MUST BE A MULTIPLE #### '
449	write(,) ' #### OF THE SYNCHRONISATION INTERVAL #### '
450	stop
451	endif
452
453	if ((loutaver/lsynctime).lt.2) then
454	write(,) ' #### FLEXPART MODEL ERROR! AVERAGING TIME OF #### '
455	write(,) ' #### CONCENTRATION FIELD MUST BE AT LEAST #### '
456	write(,) ' #### TWICE THE SYNCHRONISATION INTERVAL #### '
457	stop
458	endif
459
460	if (mod(loutstep,lsynctime).ne.0) then
461	write(,) ' #### FLEXPART MODEL ERROR! INTERVAL BETWEEN #### '
462	write(,) ' #### CONCENTRATION FIELDS MUST BE A MULTIPLE #### '
463	write(,) ' #### OF THE SYNCHRONISATION INTERVAL #### '
464	stop
465	endif
466
467	if ((loutstep/lsynctime).lt.2) then
468	write(,) ' #### FLEXPART MODEL ERROR! INTERVAL BETWEEN #### '
469	write(,) ' #### CONCENTRATION FIELDS MUST BE AT LEAST #### '
470	write(,) ' #### TWICE THE SYNCHRONISATION INTERVAL #### '
471	stop
472	endif
473
474	if (mod(loutsample,lsynctime).ne.0) then
475	write(,) ' #### FLEXPART MODEL ERROR! SAMPLING TIME OF #### '
476	write(,) ' #### CONCENTRATION FIELD MUST BE A MULTIPLE #### '
477	write(,) ' #### OF THE SYNCHRONISATION INTERVAL #### '
478	stop
479	endif
480
481	if (itsplit.lt.loutaver) then
482	write(,) ' #### FLEXPART MODEL ERROR! SPLITTING TIME FOR#### '
483	write(,) ' #### PARTICLES IS TOO SHORT. PLEASE INCREASE #### '
484	write(,) ' #### SPLITTING TIME CONSTANT. #### '
485	stop
486	endif
487
488	if ((mquasilag.eq.1).and.(iout.ge.4)) then
489	write(,) ' #### FLEXPART MODEL ERROR! CONFLICTING #### '
490	write(,) ' #### OPTIONS: IF MQUASILAG=1, PLUME #### '
491	write(,) ' #### TRAJECTORY OUTPUT IS IMPOSSIBLE. #### '
492	stop
493	endif
494
495	! Compute modeling time in seconds and beginning date in Julian date
496	!*******************************************************************
497
498	outstep=real(abs(loutstep))
499	if (ldirect.eq.1) then
500	bdate=juldate(ibdate,ibtime)
501	edate=juldate(iedate,ietime)
502	ideltas=nint((edate-bdate)*86400.)
503	else if (ldirect.eq.-1) then
504	loutaver=-1*loutaver
505	loutstep=-1*loutstep
506	loutsample=-1*loutsample
507	lsynctime=-1*lsynctime
508	bdate=juldate(iedate,ietime)
509	edate=juldate(ibdate,ibtime)
510	ideltas=nint((edate-bdate)*86400.)
511	else
512	write(,) ' #### FLEXPART MODEL ERROR! DIRECTION IN #### '
513	write(,) ' #### FILE "COMMAND" MUST BE EITHER -1 OR 1. #### '
514	stop
515	endif
516
517	return
518
519	999 write(,) ' #### FLEXPART MODEL ERROR! FILE "COMMAND" #### '
520	write(,) ' #### CANNOT BE OPENED IN THE DIRECTORY #### '
521	write(*,'(a)') path(1)(1:length(1))
522	stop
523
524	end subroutine readcommand

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

Download in other formats:

Original Format