Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

readinput.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 11 years ago
sources for flexwrf v3.1
File size: 85.0 KB

Line
1	!***********************************************************************
2	!* Copyright 2012,2013 *
3	!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine, *
4	!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
5	!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso, *
6	!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann, *
7	!* *
8	!* This file is part of FLEXPART WRF *
9	!* *
10	!* FLEXPART is free software: you can redistribute it 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 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. If not, see <http://www.gnu.org/licenses/>. *
22	!***********************************************************************
23	subroutine readinput
24	!*******************************************************************************
25	! *
26	! Reads the pathnames, where input/output files are expected to be. *
27	! The file pathnames must be available in the current working directory. *
28	! *
29	! Author: A. Stohl *
30	! *
31	! 1 February 1994 *
32	! *
33	!*******************************************************************************
34	! *
35	! Variables: *
36	! len(numpath) lengths of the path names *
37	! path(numpath) pathnames of input/output files *
38	! *
39	! Constants: *
40	! numpath number of pathnames to be read in *
41	! *
42	!*******************************************************************************
43	use par_mod
44	use com_mod
45	use outg_mod
46	use point_mod
47	! use xmass_mod
48
49	implicit none
50	! include 'includepar'
51	! include 'includecom'
52
53	integer :: i, j, numtable,numpoint2
54	integer :: hhh,mi,ss,pos_spec
55	character(len=50) :: line
56	logical :: old
57	integer :: idiff,ldat,ltim,wftime1(maxwf),numbwfn(maxnests),k
58	integer :: wftime1n(maxnests,maxwf),wftimen(maxnests,maxwf)
59	real(kind=dp) :: jul,juldate,beg,end,jul1,jul2
60	! double precision juldate,jul,beg,end,jul1,jul2
61	character(len=80) :: fname,wfname1(maxwf),wfname1n(maxnests,maxwf)
62	character(len=10) :: spec, wfspec1(maxwf),wfspec1n(maxnests,maxwf)
63	real :: outhelp,xr,xr1,yr,yr1
64	real :: xtmp, xtmp1, xtmp2, ytmp, ytmp1, ytmp2
65	! 10-mar-2006 rce - flexpart_wrf - eps should be a small dx/dy value in meters
66	! parameter(eps=1.e-4)
67	real,parameter :: eps=10.0
68	real :: x,y,xm,ym
69	character(len=16) :: receptor
70	integer :: numpartmax,id1,it1,id2,it2,idow,ihour
71	integer :: emitvar,stat,nparttot,nparttot2,outgriddef,outgriddefn
72	real :: vsh(ni),fracth(ni),schmih(ni),releaserate,releaserate2
73	character(len=3) :: aspecnumb
74	character(len=10) :: specname(maxtable)
75	real :: cun
76	logical :: spec_found
77	! Read the pathname information stored in unitpath
78	!*************************************************
79	nparttot=0
80	! print*,'path2'
81	! print*,'path',unitpath
82	! open(unitpath,file='flexwrf.input',status='old',err=801)
83	open(unitpath,file=inputname,status='old',err=801)
84
85	! jdf start read pathnames
86	!*************************************************
87	call skplin(1,unitpath)
88	do i=1,numpath
89	read(unitpath,'(a)',err=800) path(i)
90	length(i)=index(path(i),' ')-1
91	enddo
92
93	! Check whether any nested subdomains are to be used
94	!***************************************************
95
96	do i=1,maxnests
97	read(unitpath,'(a)') path(numpath+2*(i-1)+1)
98	read(unitpath,'(a)') path(numpath+2*(i-1)+2)
99	if (path(numpath+2*(i-1)+1)(1:5).eq.'=====') goto 30
100	length(numpath+2(i-1)+1)=index(path(numpath+2(i-1)+1),' ')-1
101	length(numpath+2(i-1)+2)=index(path(numpath+2(i-1)+2),' ')-1
102	enddo
103
104	if (path(numpath+2*(i-1)+1)(1:5).ne.'=====') then
105	write(,) ' #### FLEXPART MODEL ERROR! NUMBER OF PATHS #### '
106	write(,) ' #### IN FORMER PATHNAMES FILE GREATER THAN #### '
107	write(,) ' #### MAXNESTS IN PAR_MOD.F90, INCREASE THIS #### '
108	write(,) ' #### OR REMOVE EXTRA PATHS FROM PATHNAMES #### '
109	stop
110	endif
111
112	! Determine number of available nested domains
113	!*********************************************
114
115	30 numbnests=i-1
116
117	! jdf end read pathnames
118	!*************************************************
119	! jdf start read command
120	!*************************************************
121	!
122	!*******************************************************************************
123	! *
124	! Note: This is the FLEXPART_WRF version of subroutine readcommand. *
125	! *
126	! This routine reads the user specifications for the current model run. *
127	! *
128	! Author: A. Stohl *
129	! *
130	! 18 May 1996 *
131	! *
132	! Nov-Dec-2005, R. Easter - input turb_option, add_sfc_level, iouttype *
133	! *
134	!*******************************************************************************
135	! *
136	! Variables: *
137	! bdate beginning date as Julian date *
138	! ctl factor by which time step must be smaller than *
139	! Lagrangian time scale *
140	! edate ending date as Julian date *
141	! hhh hour *
142	! ibdate,ibtime beginnning date and time (YYYYMMDD, HHMISS) *
143	! ideltas [s] modelling period *
144	! iedate,ietime ending date and time (YYYYMMDD, HHMISS) *
145	! ifine reduction factor for vertical wind time step *
146	! iflux switch to turn on (1)/off (0) flux calculations *
147	! iout 1 for conc. (residence time for backward runs) output, *
148	! 2 for mixing ratio output, 3 both, 4 for plume *
149	! trajectory output, 5 = options 1 and 4 *
150	! ipin 1 continue simulation with dumped particle data, 0 no *
151	! ipout 0 no particle dump, 1 every output time, 3 only at end *
152	! itsplit [s] time constant for particle splitting *
153	! loutaver [s] concentration output is an average over loutaver seconds*
154	! loutsample [s] average is computed from samples taken every [s] seconds*
155	! loutstep [s] time interval of concentration output *
156	! lsynctime [s] synchronisation time interval for all particles *
157	! lagespectra switch to turn on (1)/off (0) calculation of age spectra*
158	! lconvection value of either 0 and 1 indicating mixing by convection *
159	! = 0 .. no convection *
160	! + 1 .. parameterisation of mixing by subgrid-scale *
161	! convection = on *
162	! lsubgrid switch to turn on (1)/off (0) subgrid topography *
163	! parameterization *
164	! method method used to compute the particle pseudovelocities *
165	! mdomainfill 1 use domain-filling option, 0 not, 2 use strat. O3 *
166	! mi minute *
167	! ss second *
168	! *
169	! Constants: *
170	! unitcommand unit connected to file COMMAND *
171	!
172	! 9-10-2007, W Wang *
173	! add turb_option_tke and turb_option_mytke
174	! 1 Oct, 2007
175	! add dt_conv input, time intervals to call convection, seconds
176	!*******************************************************************************
177
178
179	! Open the command file and read user options
180	!********************************************
181
182	read(unitcommand,*) ldirect
183	read(unitcommand,*) ibdate,ibtime
184	read(unitcommand,*) iedate,ietime
185	read(unitcommand,*) loutstep
186	read(unitcommand,*) loutaver
187	read(unitcommand,*) loutsample
188	read(unitcommand,*) itsplit
189	read(unitcommand,*) lsynctime
190	read(unitcommand,*) ctl
191	read(unitcommand,*) ifine
192	read(unitcommand,*) iout
193	read(unitcommand,*) ipout
194	read(unitcommand,*) lsubgrid
195	read(unitcommand,*) lconvection
196	read(unitcommand,*) dt_conv
197	read(unitcommand,*) lagespectra
198	read(unitcommand,*) ipin
199	read(unitcommand,*) iflux
200	read(unitcommand,*) ioutputforeachrelease
201	read(unitcommand,*) mdomainfill
202	read(unitcommand,*) ind_source
203	read(unitcommand,*) ind_receptor
204	read(unitcommand,*) nested_output
205	read(unitcommand,*) linit_cond
206	! if(nested_output.ge.1) then
207	! write(*,'(/a/a/a/)') &
208	! '* Nested grid output is not fully implemented *', &
209	! '* Set NESTED_OUTPUT=0 in COMMAND file *'
210	! stop
211	! endif
212	! FLEXPART_WRF - read turb_option, add_sfc_level
213	read(unitcommand,*) turb_option
214	read(unitcommand,*) cblflag ! added by mc for cbl option
215	! read(unitcommand,*) add_sfc_level
216	add_sfc_level=1
217	read(unitcommand,*) sfc_option
218	read(unitcommand,*) wind_option
219	read(unitcommand,*) time_option
220	read(unitcommand,*) outgrid_option
221	read(unitcommand,*) numpoint_option
222	read(unitcommand,*) iouttype
223	read(unitcommand,*) ncnumrec
224	read(unitcommand,*) option_verbose
225
226	ifine=max(ifine,1)
227
228	! Determine how Markov chain is formulated (for w or for w/sigw)
229	!***************************************************************
230
231	if (ctl.ge.0.1) then
232	turbswitch=.true.
233	else
234	turbswitch=.false.
235	ifine=1
236	endif
237	fine=1./real(ifine)
238	ctl=1./ctl
239
240	! Set the switches required for the various options for input/output units
241	!*************************************************************************
242	!AF Set the switches IND_REL and IND_SAMP for the release and sampling
243	!Af switches for the releasefile:
244	!Af IND_REL = 1 : xmass * rho
245	!Af IND_REL = 0 : xmass * 1
246
247	!Af switches for the conccalcfile:
248	!AF IND_SAMP = 0 : xmass * 1
249	!Af IND_SAMP = -1 : xmass / rho
250
251	!AF IND_SOURCE switches between different units for concentrations at the source
252	!Af NOTE that in backward simulations the release of computational particles
253	!Af takes place at the "receptor" and the sampling of p[articles at the "source".
254	!Af 1 = mass units
255	!Af 2 = mass mixing ratio units
256	!Af IND_RECEPTOR switches between different units for concentrations at the receptor
257	!Af 1 = mass units
258	!Af 2 = mass mixing ratio units
259
260	if ( ldirect .eq. 1 ) then ! FWD-Run
261	!Af set release-switch
262	if ( IND_SOURCE .eq. 1 ) then !mass
263	ind_rel = 0
264	else ! mass mix
265	ind_rel = 1
266	endif
267	!Af set sampling switch
268	if ( IND_RECEPTOR .eq. 1) then !mass
269	ind_samp = 0
270	else ! mass mix
271	ind_samp = -1
272	endif
273	elseif (ldirect .eq. -1 ) then !BWD-Run
274	!Af set sampling switch
275	if ( IND_SOURCE .eq. 1 ) then !mass
276	ind_samp = -1
277	else ! mass mix
278	ind_samp = 0
279	endif
280	!Af set release-switch
281	if ( IND_RECEPTOR .eq. 1) then !mass
282	ind_rel = 1
283	else ! mass mix
284	ind_rel = 0
285	endif
286	endif
287
288
289
290	!*************************************************************
291	! Check whether valid options have been chosen in file COMMAND
292	!*************************************************************
293
294	! Check input dates
295	!******************
296
297	if (iedate.lt.ibdate) then
298	write(,) ' #### FLEXPART MODEL ERROR! BEGINNING DATE #### '
299	write(,) ' #### IS LARGER THAN ENDING DATE. CHANGE #### '
300	write(,) ' #### EITHER POINT 2 OR POINT 3 #### '
301	stop
302	else if (iedate.eq.ibdate) then
303	if (ietime.lt.ibtime) then
304	write(,) ' #### FLEXPART MODEL ERROR! BEGINNING TIME #### '
305	write(,) ' #### IS LARGER THAN ENDING TIME. CHANGE #### '
306	write(,) ' #### EITHER POINT 2 OR POINT 3 #### '
307	stop
308	endif
309	endif
310
311
312	! Determine kind of dispersion method
313	!************************************
314
315	if (ctl.gt.0.) then
316	method=1
317	mintime=minstep
318	else
319	method=0
320	mintime=lsynctime
321	endif
322
323	! Check whether a valid option for gridded model output has been chosen
324	!**********************************************************************
325
326	! if ((iout.lt.0).or.(iout.gt.5)) then
327	! write(,) ' #### FLEXPART MODEL ERROR! FILE COMMAND: #### '
328	! write(,) ' #### IOUT MUST BE 0, 1, 2, 3, 4, OR 5! #### '
329	! stop
330	! endif
331
332	!AF check consistency between units and volume mixing ratio
333	if ( ((iout.eq.2).or.(iout.eq.3)).and. &
334	(ind_source.gt.1 .or.ind_receptor.gt.1) ) then
335	write(,) ' #### FLEXPART MODEL ERROR! : #### '
336	write(,) ' #### VOLUME MIXING RATIO ONLY SUPPORTED #### '
337	write(,) ' #### FOR MASS UNITS (at the moment) #### '
338	stop
339	endif
340
341	! For quasilag output for each release is forbidden
342	!*****************************************************************************
343
344	if ((ioutputforeachrelease.eq.1).and.(mquasilag.eq.1)) then
345	write(,) '#### FLEXPART MODEL ERROR! : ####'
346	write(,) '#### OUTPUTFOREACHRELEASE AND QUASILAGRANGIAN####'
347	write(,) '#### MODE IS NOT POSSIBLE ! ####'
348	stop
349	endif
350
351	! For quasilag backward is forbidden
352	!*****************************************************************************
353
354	if ((ldirect.lt.0).and.(mquasilag.eq.1)) then
355	write(,) '#### FLEXPART MODEL ERROR! : ####'
356	write(,) '#### FOR BACKWARD RUNS, QUASILAGRANGIAN MODE ####'
357	write(,) '#### IS NOT POSSIBLE ! ####'
358	stop
359	endif
360	! For backward runs one releasefield for all releases makes no sense,
361	! For quasilag and domainfill ioutputforechrelease is forbidden
362	!*****************************************************************************
363
364	! print*,'ioutput',ioutputforeachrelease,ldirect
365	if ((ldirect.lt.0).and.(ioutputforeachrelease.eq.0)) then
366	write(,) '#### FLEXPART MODEL ERROR! : ####'
367	write(,) '#### FOR BACKWARD RUNS, IOUTPUTFOREACHRLEASE ####'
368	write(,) '#### MUST BE SET TO ONE! ####'
369	stop
370	endif
371
372	! For backward runs one releasefield for all releases makes no sense,
373	! and is "forbidden"
374	!*****************************************************************************
375
376	if ((mdomainfill.eq.1).and.(ioutputforeachrelease.eq.1)) then
377	write(,) '#### FLEXPART MODEL ERROR! : ####'
378	write(,) '#### FOR DOMAIN FILLING RUNS OUTPUT FOR ####'
379	write(,) '#### EACH RELEASE IS FORBIDDEN ! ####'
380	stop
381	endif
382
383
384	! For domain-filling trajectories, a plume centroid trajectory makes no sense,
385	! For backward runs, only residence time output (iout=1) or plume trajectories (iout=4),
386	! or both (iout=5) makes sense; other output options are "forbidden"
387	!*****************************************************************************
388
389	if (ldirect.lt.0) then
390	if ((iout.eq.2).or.(iout.eq.3)) then
391	write(,) '#### FLEXPART MODEL ERROR! : ####'
392	write(,) '#### FOR BACKWARD RUNS, IOUT MUST BE 1,4,OR 5####'
393	stop
394	endif
395	endif
396	! For domain-filling trajectories, a plume centroid trajectory makes no sense,
397	! and is "forbidden"
398	!*****************************************************************************
399
400	if (mdomainfill.ge.1) then
401	if ((iout.eq.4).or.(iout.eq.5)) then
402	write(,) '#### FLEXPART MODEL ERROR! : ####'
403	write(,) '#### FOR DOMAIN-FILLING TRAJECTORY OPTION, ####'
404	write(,) '#### IOUT MUST NOT BE SET TO 4 OR 5. ####'
405	stop
406	endif
407	endif
408
409
410	! Check whether a valid options for particle dump has been chosen
411	!****************************************************************
412
413	if ((ipout.ne.0).and.(ipout.ne.1).and.(ipout.ne.2)) then
414	write(,) ' #### FLEXPART MODEL ERROR! : #### '
415	write(,) ' #### IPOUT MUST BE 1, 2 OR 3! #### '
416	stop
417	endif
418
419	if(lsubgrid.ne.1) then
420	write(,) ' ---------------- '
421	write(,) ' INFORMATION: SUBGRIDSCALE TERRAIN EFFECT IS'
422	write(,) ' NOT PARAMETERIZED DURING THIS SIMULATION. '
423	write(,) ' ---------------- '
424	endif
425
426	if ((ipout.ne.0).and.(ipout.ne.1).and.(ipout.ne.2)) then
427	write(,) ' #### FLEXPART MODEL ERROR! : #### '
428	write(,) ' #### IPOUT MUST BE 1, 2 OR 3! #### '
429	stop
430	endif
431
432
433	! Check whether convection scheme is either turned on or off
434	!***********************************************************
435
436	! if ((lconvection.ne.0).and.(lconvection.ne.1)) then
437	if ((lconvection.lt.0).or. (lconvection.gt.3)) then
438	write(,) ' #### FLEXPART MODEL ERROR! : #### '
439	! write(,) ' #### LCONVECTION MUST BE SET TO EITHER 1 OR 0#### '
440	write(,) ' #### LCONVECTION MUST BE SET TO 0 or 3 #### '
441
442	stop
443	endif
444
445
446	! Check whether synchronisation interval is sufficiently short
447	!*************************************************************
448
449	if (lsynctime.gt.(idiffnorm/2)) then
450	write(,) ' #### FLEXPART MODEL ERROR! SYNCHRONISATION #### '
451	write(,) ' #### TIME IS TOO LONG. MAKE IT SHORTER. #### '
452	stop
453	endif
454
455
456	! Check consistency of the intervals, sampling periods, etc., for model output
457	!*****************************************************************************
458
459	if (loutaver.eq.0) then
460	write(,) ' #### FLEXPART MODEL ERROR! TIME AVERAGE OF #### '
461	write(,) ' #### CONCENTRATION FIELD OUTPUT MUST NOT BE #### '
462	write(,) ' #### ZERO. #### '
463	write(,) ' #### CHANGE INPUT #### '
464	stop
465	endif
466
467	if (loutaver.gt.loutstep) then
468	write(,) ' #### FLEXPART MODEL ERROR! TIME AVERAGE OF #### '
469	write(,) ' #### CONCENTRATION FIELD OUTPUT MUST NOT BE #### '
470	write(,) ' #### GREATER THAN INTERVAL OF OUTPUT. #### '
471	write(,) ' #### CHANGE INPUT #### '
472	stop
473	endif
474
475	if (loutsample.gt.loutaver) then
476	write(,) ' #### FLEXPART MODEL ERROR! SAMPLING TIME OF #### '
477	write(,) ' #### CONCENTRATION FIELD OUTPUT MUST NOT BE #### '
478	write(,) ' #### GREATER THAN TIME AVERAGE OF OUTPUT. #### '
479	write(,) ' #### CHANGE INPUT #### '
480	stop
481	endif
482
483	if (mod(loutaver,lsynctime).ne.0) then
484	write(,) ' #### FLEXPART MODEL ERROR! AVERAGING TIME OF #### '
485	write(,) ' #### CONCENTRATION FIELD MUST BE A MULTIPLE #### '
486	write(,) ' #### OF THE SYNCHRONISATION INTERVAL #### '
487	stop
488	endif
489
490	if ((loutaver/lsynctime).lt.2) then
491	write(,) ' #### FLEXPART MODEL ERROR! AVERAGING TIME OF #### '
492	write(,) ' #### CONCENTRATION FIELD MUST BE AT LEAST #### '
493	write(,) ' #### TWICE THE SYNCHRONISATION INTERVAL #### '
494	stop
495	endif
496
497	if (mod(loutstep,lsynctime).ne.0) then
498	write(,) ' #### FLEXPART MODEL ERROR! INTERVAL BETWEEN #### '
499	write(,) ' #### CONCENTRATION FIELDS MUST BE A MULTIPLE #### '
500	write(,) ' #### OF THE SYNCHRONISATION INTERVAL #### '
501	stop
502	endif
503
504	if ((loutstep/lsynctime).lt.2) then
505	write(,) ' #### FLEXPART MODEL ERROR! INTERVAL BETWEEN #### '
506	write(,) ' #### CONCENTRATION FIELDS MUST BE AT LEAST #### '
507	write(,) ' #### TWICE THE SYNCHRONISATION INTERVAL #### '
508	stop
509	endif
510
511	if (mod(loutsample,lsynctime).ne.0) then
512	write(,) ' #### FLEXPART MODEL ERROR! SAMPLING TIME OF #### '
513	write(,) ' #### CONCENTRATION FIELD MUST BE A MULTIPLE #### '
514	write(,) ' #### OF THE SYNCHRONISATION INTERVAL #### '
515	stop
516	endif
517
518	if (itsplit.lt.loutaver) then
519	write(,) ' #### FLEXPART MODEL ERROR! SPLITTING TIME FOR#### '
520	write(,) ' #### PARTICLES IS TOO SHORT. PLEASE INCREASE #### '
521	write(,) ' #### SPLITTING TIME CONSTANT. #### '
522	stop
523	endif
524
525	if ((mquasilag.eq.1).and.(iout.ge.4)) then
526	write(,) ' #### FLEXPART MODEL ERROR! CONFLICTING #### '
527	write(,) ' #### OPTIONS: IF MQUASILAG=1, PLUME #### '
528	write(,) ' #### TRAJECTORY OUTPUT IS IMPOSSIBLE. #### '
529	stop
530	endif
531
532	if (ncnumrec.le.0) then
533	write(,) ' #### FLEXPART MODEL ERROR! NUMREC MUST #### '
534	write(,) ' #### BE GREATER THAN 0 #### '
535	stop
536	endif
537
538	! FLEXPART_WRF - check turb_option, add_sfc_level
539	if ((turb_option.ne.turb_option_none ) .and. &
540	(turb_option.ne.turb_option_diagnosed) .and. &
541	(turb_option.ne.turb_option_tke ) .and. &
542	(turb_option.ne.turb_option_mytke)) then
543	write(,) ' #### FLEXPART MODEL ERROR! #### '
544	write(,) ' #### TURB_OPTION MUST BE ONE OF: #### '
545	write(*,'(5x,5i5)') turb_option_none, turb_option_diagnosed, &
546	turb_option_tke,turb_option_mytke
547	write(,) ' #### --------------------------------- #### '
548	stop
549	endif
550
551	if ((add_sfc_level.ne.0) .and. (add_sfc_level.ne.1)) then
552	write(,) ' #### FLEXPART MODEL ERROR! #### '
553	write(,) ' #### ADD_SFC_LAYER MUST BE 0 or 1 #### '
554	stop
555	endif
556
557	if ((sfc_option.ne.sfc_option_diagnosed) .and. &
558	(sfc_option.ne.sfc_option_wrf )) then
559	write(,) ' #### FLEXPART MODEL ERROR! #### '
560	write(,) ' #### SFC_OPTION MUST BE ONE OF: #### '
561	write(*,'(5x,5i5)') sfc_option_diagnosed, sfc_option_wrf
562	write(,) ' #### --------------------------------- #### '
563	stop
564	endif
565
566	! iouttype -- convert negative values to 0; positive out of range values to 2
567	if (iouttype .lt. 0) iouttype = 0
568	if (iouttype .gt. 2) iouttype = 2
569
570	! Conversion of format HHHMISS to seconds
571	!****************************************
572
573	hhh=ideltas/10000
574	mi=(ideltas-10000*hhh)/100
575	ss=ideltas-10000hhh-100mi
576	ideltas=hhh3600+60mi+ss
577
578
579	! Compute modeling time in seconds and beginning date in Julian date
580	!*******************************************************************
581
582	outstep=real(abs(loutstep))
583	if (ldirect.eq.1) then
584	bdate=juldate(ibdate,ibtime)
585	edate=juldate(iedate,ietime)
586	ideltas=nint((edate-bdate)*86400.)
587	else if (ldirect.eq.-1) then
588	loutaver=-1*loutaver
589	loutstep=-1*loutstep
590	loutsample=-1*loutsample
591	lsynctime=-1*lsynctime
592	bdate=juldate(iedate,ietime)
593	edate=juldate(ibdate,ibtime)
594	ideltas=nint((edate-bdate)*86400.)
595	else
596	write(,) ' #### FLEXPART MODEL ERROR! DIRECTION IN #### '
597	write(,) ' #### INPUT FILE MUST BE EITHER -1 OR 1. #### '
598	stop
599	endif
600
601	! jdf end read command
602	!*************************************************
603	! jdf start read ageclasees
604	!*************************************************
605	!
606	!*******************************************************************************
607	! *
608	! This routine reads the age classes to be used for the current model run. *
609	! *
610	! Author: A. Stohl *
611	! *
612	! 20 March 2000 *
613	! *
614	!*******************************************************************************
615	! *
616	! Variables: *
617	! *
618	! Constants: *
619	! *
620	!*******************************************************************************
621
622	! If age spectra claculation is switched on,
623	! open the AGECLASSSES file and read user options
624	!************************************************
625
626	call skplin(1,unitpath)
627	read(unitpath,*) nageclass
628	read(unitpath,*) lage(1)
629	if (nageclass.gt.maxageclass) then
630	write(,) ' #### FLEXPART MODEL ERROR! NUMBER OF AGE #### '
631	write(,) ' #### CLASSES GREATER THAN MAXIMUM ALLOWED. #### '
632	write(,) ' #### CHANGE SETTINGS IN AGECLASSES OR #### '
633	write(,) ' #### RECOMPILE WITH LARGER MAXAGECLASS IN #### '
634	write(,) ' #### FILE PAR_MOD.F90 #### '
635	stop
636	endif
637
638	if (lage(1).le.0) then
639	write(,) ' #### FLEXPART MODEL ERROR! AGE OF FIRST #### '
640	write(,) ' #### CLASS MUST BE GREATER THAN ZERO. CHANGE #### '
641	write(,) ' #### SETTINGS IN AGECLASSES. #### '
642	stop
643	endif
644
645	! If age spectra calculation is switched off, set number of age classes
646	! to 1 and maximum age to a large number
647	!**********************************************************************
648
649
650	do i=2,nageclass
651	read(unitpath,*) lage(i)
652	! print*,'age',lage(i),i
653	if (lage(i).le.lage(i-1)) then
654	write(,) ' #### FLEXPART MODEL ERROR! AGE CLASSES #### '
655	write(,) ' #### MUST BE GIVEN IN TEMPORAL ORDER. #### '
656	write(,) ' #### CHANGE SETTINGS IN AGECLASSES. #### '
657	stop
658	endif
659	enddo
660	if (lagespectra.ne.1) then
661	nageclass=1
662	lage(nageclass)=999999999
663	endif
664
665	! jdf end read ageclasses
666	!*************************************************
667	! jdf start read available
668	!*************************************************
669	!
670	!*******************************************************************************
671	! *
672	! This routine reads the dates and times for which windfields are available. *
673	! *
674	! Authors: A. Stohl *
675	! *
676	! 6 February 1994 *
677	! 8 February 1999, Use of nested fields, A. Stohl *
678	! *
679	! 12 October 2005, R. Easter - *
680	! fname,wfname1,wfname1n changed from char10 to char80; *
681	! reads from unitavailab changed to free format *
682	! *
683	!*******************************************************************************
684	! *
685	! Variables: *
686	! bdate beginning date as Julian date *
687	! beg beginning date for windfields *
688	! end ending date for windfields *
689	! fname filename of wind field, help variable *
690	! ideltas [s] duration of modelling period *
691	! idiff time difference between 2 wind fields *
692	! idiffnorm normal time difference between 2 wind fields *
693	! idiffmax [s] maximum allowable time between 2 wind fields *
694	! jul julian date, help variable *
695	! numbwf actual number of wind fields *
696	! wfname(maxwf) file names of needed wind fields *
697	! wfspec(maxwf) file specifications of wind fields (e.g., if on disc) *
698	! wftime(maxwf) [s]times of wind fields relative to beginning time *
699	! wfname1,wfspec1,wftime1 = same as above, but only local (help variables) *
700	! *
701	! Constants: *
702	! maxwf maximum number of wind fields *
703	! unitavailab unit connected to file AVAILABLE *
704	! *
705	!*******************************************************************************
706
707	! Windfields are only used, if they are within the modelling period.
708	! However, 1 additional day at the beginning and at the end is used for
709	! interpolation. -> Compute beginning and ending date for the windfields.
710	!************************************************************************
711
712	if (ideltas.gt.0) then ! forward trajectories
713	beg=bdate-1.
714	end=bdate+dble(real(ideltas)/86400.)+dble(real(idiffmax)/ &
715	86400.)
716	else ! backward trajectories
717	beg=bdate+dble(real(ideltas)/86400.)-dble(real(idiffmax)/ &
718	86400.)
719	end=bdate+1.
720	endif
721
722	! Open the wind field availability file and read available wind fields
723	! within the modelling period.
724	!*********************************************************************
725
726	open(unitavailab,file=path(3)(1:length(3)),status='old', &
727	err=804)
728
729	do i=1,3
730	read(unitavailab,*)
731	enddo
732
733	numbwf=0
734	100 read(unitavailab,*,end=99) ldat,ltim,fname,spec
735	jul=juldate(ldat,ltim)
736	if ((jul.ge.beg).and.(jul.le.end)) then
737	numbwf=numbwf+1
738	if (numbwf.gt.maxwf) then ! check exceedance of dimension
739	write(,) 'Number of wind fields needed is too great.'
740	write(,) 'Reduce modelling period (file "COMMAND") or'
741	write(,) 'reduce number of wind fields (file "AVAILABLE").'
742	stop
743	endif
744
745	wfname1(numbwf)=fname
746	wfspec1(numbwf)=spec
747	wftime1(numbwf)=nint((jul-bdate)*86400.)
748	endif
749	goto 100 ! next wind field
750
751	99 continue
752
753	close(unitavailab)
754
755	! Open the wind field availability file and read available wind fields
756	! within the modelling period (nested grids)
757	!*********************************************************************
758
759	do k=1,numbnests
760	open(unitavailab,file=path(numpath+2*(k-1)+2) &
761	(1:length(numpath+2*(k-1)+2)),status='old',err=803)
762
763	do i=1,3
764	read(unitavailab,*)
765	enddo
766
767	numbwfn(k)=0
768	700 read(unitavailab,*,end=699) ldat,ltim,fname,spec
769	jul=juldate(ldat,ltim)
770	if ((jul.ge.beg).and.(jul.le.end)) then
771	numbwfn(k)=numbwfn(k)+1
772	if (numbwfn(k).gt.maxwf) then ! check exceedance of dimension
773	write(,) 'Number of nested wind fields is too great.'
774	write(,) 'Reduce modelling period (file "COMMAND") or'
775	write(,) 'reduce number of wind fields (file "AVAILABLE").'
776	stop
777	endif
778
779	wfname1n(k,numbwfn(k))=fname
780	wfspec1n(k,numbwfn(k))=spec
781	wftime1n(k,numbwfn(k))=nint((jul-bdate)*86400.)
782	endif
783	goto 700 ! next wind field
784
785	699 continue
786
787	enddo
788	close(unitavailab)
789
790
791	! Check wind field times of file AVAILABLE (expected to be in temporal order)
792	!****************************************************************************
793
794	if (numbwf.eq.0) then
795	write(,) ' #### FLEXPART MODEL ERROR! NO WIND FIELDS #### '
796	write(,) ' #### AVAILABLE FOR SELECTED TIME PERIOD. #### '
797	stop
798	endif
799
800	do i=2,numbwf
801	if (wftime1(i).le.wftime1(i-1)) then
802	write(,) 'FLEXPART ERROR: FILE AVAILABLE IS CORRUPT.'
803	write(,) 'THE WIND FIELDS ARE NOT IN TEMPORAL ORDER.'
804	write(,) 'PLEASE CHECK FIELD ',wfname1(i)
805	stop
806	endif
807	enddo
808
809	! Check wind field times of file AVAILABLE for the nested fields
810	! (expected to be in temporal order)
811	!***************************************************************
812
813	do k=1,numbnests
814	if (numbwfn(k).eq.0) then
815	write(,) '#### FLEXPART MODEL ERROR! NO WIND FIELDS ####'
816	write(,) '#### AVAILABLE FOR SELECTED TIME PERIOD. ####'
817	stop
818	endif
819
820	do i=2,numbwfn(k)
821	if (wftime1n(k,i).le.wftime1n(k,i-1)) then
822	write(,) 'FLEXTRA ERROR: FILE AVAILABLE IS CORRUPT. '
823	write(,) 'THE NESTED WIND FIELDS ARE NOT IN TEMPORAL ORDER.'
824	write(,) 'PLEASE CHECK FIELD ',wfname1n(k,i)
825	write(,) 'AT NESTING LEVEL ',k
826	stop
827	endif
828	enddo
829	enddo
830
831
832	! For backward trajectories, reverse the order of the windfields
833	!***************************************************************
834
835	do i=1,numbwf
836	wfdt(i)=-999999
837	enddo
838	if (ideltas.ge.0) then
839	do i=1,numbwf
840	wfname(i)=wfname1(i)
841	wfspec(i)=wfspec1(i)
842	wftime(i)=wftime1(i)
843	if(i.gt.1) wfdt(i)=wftime1(i)-wftime1(i-1)
844	enddo
845	do k=1,numbnests
846	do i=1,numbwfn(k)
847	wfnamen(k,i)=wfname1n(k,i)
848	wfspecn(k,i)=wfspec1n(k,i)
849	wftimen(k,i)=wftime1n(k,i)
850	enddo
851	enddo
852	else
853	do i=1,numbwf
854	wfname(numbwf-i+1)=wfname1(i)
855	wfspec(numbwf-i+1)=wfspec1(i)
856	wftime(numbwf-i+1)=wftime1(i)
857	! if(i.lt.numbwf) wfdt(numbwf-i+1)=wftime1(i+1)-wftime1(i)
858	if(i.gt.1) wfdt(numbwf-i+1)=wftime1(i)-wftime1(i-1)
859	enddo
860	do k=1,numbnests
861	do i=1,numbwfn(k)
862	wfnamen(k,numbwfn(k)-i+1)=wfname1n(k,i)
863	wfspecn(k,numbwfn(k)-i+1)=wfspec1n(k,i)
864	wftimen(k,numbwfn(k)-i+1)=wftime1n(k,i)
865	enddo
866	enddo
867	endif
868
869	! Check the time difference between the wind fields. If it is big,
870	! write a warning message. If it is too big, terminate the trajectory.
871	!*********************************************************************
872
873	do i=2,numbwf
874	idiff=abs(wftime(i)-wftime(i-1))
875	if (idiff.gt.idiffmax) then
876	write(,) 'FLEXPART WARNING: TIME DIFFERENCE BETWEEN TWO'
877	write(,) 'WIND FIELDS IS TOO BIG FOR TRANSPORT CALCULATION'
878	write(,) 'THEREFORE, TRAJECTORIES HAVE TO BE SKIPPED.'
879	else if (idiff.gt.idiffnorm) then
880	write(,) 'FLEXPART WARNING: TIME DIFFERENCE BETWEEN TWO'
881	write(,) 'WIND FIELDS IS BIG. THIS MAY CAUSE A DEGRADATION'
882	write(,) 'OF SIMULATION QUALITY.'
883	endif
884	enddo
885
886	do k=1,numbnests
887	if (numbwfn(k).ne.numbwf) then
888	write(,) 'FLEXTRA ERROR: THE AVAILABLE FILES FOR THE'
889	write(,) 'NESTED WIND FIELDS ARE NOT CONSISTENT WITH'
890	write(,) 'THE AVAILABLE FILE OF THE MOTHER DOMAIN. '
891	write(,) 'ERROR AT NEST LEVEL: ',k
892	stop
893	endif
894	do i=1,numbwf
895	if (wftimen(k,i).ne.wftime(i)) then
896	write(,) 'FLEXTRA ERROR: THE AVAILABLE FILES FOR THE'
897	write(,) 'NESTED WIND FIELDS ARE NOT CONSISTENT WITH'
898	write(,) 'THE AVAILABLE FILE OF THE MOTHER DOMAIN. '
899	write(,) 'ERROR AT NEST LEVEL: ',k
900	stop
901	endif
902	enddo
903	enddo
904
905	! Reset the times of the wind fields that are kept in memory to no time
906	!**********************************************************************
907
908	do i=1,2
909	memind(i)=i
910	memtime(i)=999999999
911	enddo
912
913	call gridcheck()
914	call gridcheck_nests()
915
916	! jdf end read available
917	!*************************************************
918	! jdf start read outgrid
919	!*************************************************
920	!
921	!*******************************************************************************
922	! *
923	! Note: This is the FLEXPART_WRF version of subroutine readoutgrid. *
924	! *
925	! This routine reads the user specifications for the output grid. *
926	! *
927	! Author: A. Stohl *
928	! *
929	! 4 June 1996 *
930	! *
931	! Dec 2005, R. Easter - *
932	! The output grid is defined by specifying its southwest and *
933	! northease corners, either in degrees-latlon or grid-meters *
934	! Changes to some read formats (wider fields). *
935	! Changed names of "lon0" & "lat0" variables *
936	! 10 Mar 2006, R. Easter - *
937	! Change eps from 1.0e-4 (degrees value) to 10.0 (meters value) *
938	! *
939	!*******************************************************************************
940	! *
941	! Variables: *
942	! dxout,dyout grid distance *
943	! numxgrid,numygrid,numzgrid grid dimensions *
944	! out_xm0,out_ym0 lower left corner of grid *
945	! outheight(maxzgrid) height levels of output grid [m] *
946	! *
947	! Constants: *
948	! unitpath unit connected to file OUTGRID *
949	! *
950	!*******************************************************************************
951
952	! Open the OUTGRID file and read output grid specifications
953	!**********************************************************
954
955	call skplin(1,unitpath)
956
957
958	! 1. Read horizontal grid specifications
959	!
960	! * NOTE *
961	! [xtmp1, ytmp1] are the coordinates of the southwest corner
962	! of the first (i.e., southwest or lower-left) output grid cell
963	! [xtmp2, ytmp2] are the coordinates of the northeast corner
964	! of the last (i.e,, northeast or upper-right) output grid cell
965	!****************************************
966
967	! read(unitpath,'(f15.8)') xtmp1
968	read(unitpath,*) xtmp1
969	! read(unitpath,'(f15.8)') ytmp1
970	read(unitpath,*) ytmp1
971	! read(unitpath,'(2x,i7)') numxgrid
972	read(unitpath,*) numxgrid
973	! read(unitpath,'(2x,i7)') numygrid
974	read(unitpath,*) numygrid
975	! read(unitpath,'(2x,i4)') outgriddef
976	read(unitpath,*) outgriddef
977	if (outgriddef.eq.1) then
978	! read(unitpath,'(f15.8)') xtmp2
979	read(unitpath,*) xtmp2
980	! read(unitpath,'(f15.8)') ytmp2
981	read(unitpath,*) ytmp2
982	else
983	! read(unitpath,'(f15.8)') dxout
984	! read(unitpath,'(f15.8)') dyout
985	read(unitpath,*) dxout
986	read(unitpath,*) dyout
987	xtmp2=dxout*real(numxgrid)+xtmp1
988	ytmp2=dyout*real(numygrid)+ytmp1
989	endif
990
991	if (option_verbose.eq. 1) then
992	write(*,'(/a)') 'readoutgrid diagnostics'
993	write(*,'(a,1p,2e18.10)') 'xtmp1, ytmp1 (in)', xtmp1, ytmp1
994	write(*,'(a,1p,2e18.10)') 'xtmp2, ytmp2 (in)', xtmp2, ytmp2
995	endif
996	if (outgrid_option .eq. 1) then
997	! In this case, the above inputs are the actual geographical lat/lon
998	! of the southwest & northeast corners of the output grid
999	! Need to convert from lat/lon to grid-meters
1000	outgrid_swlon = xtmp1
1001	outgrid_swlat = ytmp1
1002
1003	outgrid_nelon = xtmp2
1004	outgrid_nelat = ytmp2
1005	call ll_to_xymeter_wrf( outgrid_swlon, outgrid_swlat, &
1006	out_xm0, out_ym0 )
1007	call ll_to_xymeter_wrf( outgrid_nelon, outgrid_nelat, &
1008	xtmp, ytmp )
1009	! 10-mar-2006 rce
1010	! If out_xm0 is very close to mother grid sw corner (=xmet0), set it to that
1011	! If xtmp is very close to mother grid ne corner (=xmet0+nxmin1*dx), set it to that
1012	! Do similar for out_ym0 & ytmp
1013	if (abs(out_xm0-xmet0) .le. eps) out_xm0 = xmet0
1014	if (abs(out_ym0-ymet0) .le. eps) out_ym0 = ymet0
1015	xr1 = xmet0 + real(nxmin1)*dx
1016	yr1 = ymet0 + real(nymin1)*dy
1017	if (abs(xtmp-xr1) .le. eps) xtmp = xr1
1018	if (abs(ytmp-yr1) .le. eps) ytmp = yr1
1019	dxout = (xtmp - out_xm0)/numxgrid
1020	dyout = (ytmp - out_ym0)/numygrid
1021	if (outgrid_option.eq.1) then ! regular
1022	outlat0=outgrid_swlat
1023	outlon0=outgrid_swlon
1024	dyoutl=(outgrid_nelat-outgrid_swlat)/numygrid
1025	dxoutl= (outgrid_nelon-outgrid_swlon)/numxgrid
1026	endif
1027	else
1028	! In this case, the above inputs are in grid-meters
1029	! Need to convert from grid-meters to lat/lon
1030	out_xm0 = xtmp1
1031	out_ym0 = ytmp1
1032	dxout = (xtmp2 - xtmp1)/numxgrid
1033	dyout = (ytmp2 - ytmp1)/numygrid
1034	call xymeter_to_ll_wrf( xtmp1, ytmp1, &
1035	outgrid_swlon, outgrid_swlat )
1036	call xymeter_to_ll_wrf( xtmp2, ytmp2, &
1037	outgrid_nelon, outgrid_nelat )
1038	if (option_verbose.eq. 1) then
1039	write(*,'(f15.10,5x,a)') outgrid_swlon, 'outgrid_swlon'
1040	write(*,'(f15.10,5x,a)') outgrid_swlat, 'outgrid_swlat'
1041	write(*,'(f15.10,5x,a)') outgrid_nelon, 'outgrid_nelon'
1042	write(*,'(f15.10,5x,a)') outgrid_nelat, 'outgrid_nelat'
1043	endif
1044	end if
1045	if (option_verbose.eq. 1) then
1046	write(*,'(a,1p,2e18.10)') 'out_xm0, out_ym0 ', out_xm0, out_ym0
1047	write(*,'(a,1p,2e18.10)') 'dxout, dyout ', dxout, dyout
1048	endif
1049
1050	! Check validity of output grid (shall be within model domain)
1051	!*************************************************************
1052
1053	xr=out_xm0+real(numxgrid)*dxout
1054	yr=out_ym0+real(numygrid)*dyout
1055	xr1=xmet0+real(nxmin1+1)*dx
1056	yr1=ymet0+real(nymin1+1)*dy
1057	if (outgrid_option.lt.1) then
1058	! print*,'nx',nxmin1,nymin1
1059	if ((out_xm0+eps.lt.xmet0).or.(out_ym0+eps.lt.ymet0) &
1060	.or.(xr.gt.xr1+eps).or.(yr.gt.yr1+eps)) then
1061	write(,) ' #### FLEXPART MODEL ERROR! PART OF OUTPUT ####'
1062	write(,) ' #### GRID IS OUTSIDE MODEL DOMAIN. CHANGE ####'
1063	write(,) ' #### OUTGRID IN INPUT FILE ####'
1064	write(*,'(a)') path(1)(1:length(1))
1065	stop
1066	endif
1067	endif
1068	! if ((numxgrid.gt.maxxgrid).or.(numygrid.gt.maxygrid)) then
1069	! write(,) ' #### FLEXPART MODEL ERROR! DIMENSIONS OF ####'
1070	! write(,) ' #### OUTPUT GRID EXCEED MAXIMUM VALUES. ####'
1071	! write(,) ' #### CHANGE FILE $FLEXPART/options/OUTGRID. ####'
1072	! stop
1073	! endif
1074
1075	! 2. Vertical levels of output grid
1076	!**********************************
1077
1078	read(unitpath,*) numzgrid
1079	! if (numzgrid.gt.maxzgrid) then
1080	! write(,) ' #### FLEXPART MODEL ERROR! TOO MANY HEIGHT #### '
1081	! write(,) ' #### LEVELS ARE GIVEN FOR OUTPUT GRID. #### '
1082	! write(,) ' #### MAXIMUM NUMBER IS ',maxzgrid,' #### '
1083	! write(,) ' #### PLEASE MAKE CHANGES IN FILE OUTGRID. #### '
1084	! stop
1085	! endif
1086
1087	allocate(outheight(numzgrid) &
1088	,stat=stat)
1089	if (stat.ne.0) write(,)'ERROR: could not allocate outh'
1090	allocate(outheighthalf(numzgrid) &
1091	,stat=stat)
1092	if (stat.ne.0) write(,)'ERROR: could not allocate outh'
1093
1094	allocate(oroout(0:numxgrid-1,0:numygrid-1) &
1095	,stat=stat)
1096	if (stat.ne.0) write(,)'ERROR: could not allocate outh'
1097	allocate(area(0:numxgrid-1,0:numygrid-1) &
1098	,stat=stat)
1099	if (stat.ne.0) write(,)'ERROR: could not allocate outh'
1100	allocate(volume(0:numxgrid-1,0:numygrid-1,numzgrid) &
1101	,stat=stat)
1102	if (stat.ne.0) write(,)'ERROR: could not allocate outh'
1103	allocate(areaeast(0:numxgrid-1,0:numygrid-1,numzgrid) &
1104	,stat=stat)
1105	if (stat.ne.0) write(,)'ERROR: could not allocate outh'
1106	allocate(areanorth(0:numxgrid-1,0:numygrid-1,numzgrid) &
1107	,stat=stat)
1108	if (stat.ne.0) write(,)'ERROR: could not allocate outh'
1109
1110	do j=1,numzgrid
1111	read(unitpath,*) outhelp
1112	outheight(j)=outhelp
1113	enddo
1114
1115
1116	! Check whether vertical levels are specified in ascending order
1117	!***************************************************************
1118
1119	do j=2,numzgrid
1120	if (outheight(j).le.outheight(j-1)) then
1121	write(,) ' #### FLEXPART MODEL ERROR! YOUR SPECIFICATION#### '
1122	write(,) ' #### OF OUTPUT LEVELS IS CORRUPT AT LEVEL #### '
1123	write(,) ' #### ',j,' #### '
1124	write(,) ' #### PLEASE MAKE CHANGES IN OUTGRID. #### '
1125	endif
1126	enddo
1127
1128	! Determine the half levels, i.e. middle levels of the output grid
1129	!*****************************************************************
1130
1131	outheighthalf(1)=outheight(1)/2.
1132	do j=2,numzgrid
1133	outheighthalf(j)=(outheight(j-1)+outheight(j))/2.
1134	enddo
1135
1136	xoutshift=xmet0-out_xm0
1137	youtshift=ymet0-out_ym0
1138
1139	if(nested_output .eq. 1) then
1140	call skplin(1,unitpath)
1141	! read(unitpath,'(f15.8)') xtmp1
1142	! read(unitpath,'(f15.8)') ytmp1
1143	! read(unitpath,'(4x,i5)') numxgridn
1144	! read(unitpath,'(4x,i5)') numygridn
1145	! read(unitpath,'(2x,i4)') outgriddefn
1146	read(unitpath,*) xtmp1
1147	read(unitpath,*) ytmp1
1148	read(unitpath,*) numxgridn
1149	read(unitpath,*) numygridn
1150	read(unitpath,*) outgriddefn
1151	if (outgriddefn.eq.1) then
1152	! read(unitpath,'(f15.8)') xtmp2
1153	! read(unitpath,'(f15.8)') ytmp2
1154	read(unitpath,*) xtmp2
1155	read(unitpath,*) ytmp2
1156	else
1157	! read(unitpath,'(f15.8)') dxoutn
1158	! read(unitpath,'(f15.8)') dyoutn
1159	read(unitpath,*) dxoutn
1160	read(unitpath,*) dyoutn
1161	xtmp2=dxoutn*real(numxgridn)+xtmp1
1162	ytmp2=dyoutn*real(numygridn)+ytmp1
1163	endif
1164
1165	write(*,'(/a)') 'readoutgrid_nest diagnostics'
1166	write(*,'(a,1p,2e18.10)') 'xtmp1, ytmp1 (in)', xtmp1, ytmp1
1167	write(*,'(a,1p,2e18.10)') 'xtmp2, ytmp2 (in)', xtmp2, ytmp2
1168	if (outgrid_option .eq. 1) then
1169	! In this case, the above inputs are the actual geographical lat/lon
1170	! of the southwest & northeast corners of the output grid
1171	! Need to convert from lat/lon to grid-meters
1172	outgridn_swlon = xtmp1
1173	outgridn_swlat = ytmp1
1174	outgridn_nelon = xtmp2
1175	outgridn_nelat = ytmp2
1176	call ll_to_xymeter_wrf( outgridn_swlon, outgridn_swlat, &
1177	out_xm0n, out_ym0n )
1178	call ll_to_xymeter_wrf( outgridn_nelon, outgridn_nelat, &
1179	xtmp, ytmp )
1180	dxoutn = (xtmp - out_xm0n)/numxgridn
1181	dyoutn = (ytmp - out_ym0n)/numygridn
1182	if (outgrid_option.eq.1) then ! regular
1183	outlat0n=outgridn_swlat
1184	outlon0n=outgridn_swlon
1185	dyoutln=(outgridn_nelat-outgridn_swlat)/numygridn
1186	dxoutln= (outgridn_nelon-outgridn_swlon)/numxgridn
1187	endif
1188
1189	else
1190	! In this case, the above inputs are in grid-meters
1191	! Need to convert from grid-meters to lat/lon
1192	out_xm0n = xtmp1
1193	out_ym0n = ytmp1
1194	dxoutn = (xtmp2 - xtmp1)/numxgridn
1195	dyoutn = (ytmp2 - ytmp1)/numygridn
1196	call xymeter_to_ll_wrf( xtmp1, ytmp1, &
1197	outgridn_swlon, outgridn_swlat )
1198	call xymeter_to_ll_wrf( xtmp2, ytmp2, &
1199	outgridn_nelon, outgridn_nelat )
1200	write(*,'(f15.10,5x,a)') outgridn_swlon, 'outgridn_swlon'
1201	write(*,'(f15.10,5x,a)') outgridn_swlat, 'outgridn_swlat'
1202	write(*,'(f15.10,5x,a)') outgridn_nelon, 'outgridn_nelon'
1203	write(*,'(f15.10,5x,a)') outgridn_nelat, 'outgridn_nelat'
1204	endif
1205	write(*,'(a,1p,2e18.10)') 'out_xm0n, out_ym0n',out_xm0n,out_ym0n
1206	write(*,'(a,1p,2e18.10)') 'dxoutn, dyoutn ',dxoutn,dyoutn
1207
1208
1209	! Check validity of output grid (shall be within model domain)
1210	!*************************************************************
1211
1212	xr=out_xm0n+real(numxgridn)*dxoutn
1213	yr=out_ym0n+real(numygridn)*dyoutn
1214	xr1=xmet0+real(nxmin1+1)*dx
1215	yr1=ymet0+real(nymin1+1)*dy
1216	if ((out_xm0n+eps.lt.xmet0).or.(out_ym0n+eps.lt.ymet0) &
1217	.or.(xr.gt.xr1+eps).or.(yr.gt.yr1+eps)) then
1218	write(,) ' #### FLEXPART MODEL ERROR! PART OF OUTPUT ####'
1219	write(,) ' #### NEST IS OUTSIDE MODEL DOMAIN. CHANGE ####'
1220	write(,) ' #### OUTGRID ####'
1221	write(*,'(a)') path(1)(1:length(1))
1222	stop
1223	endif
1224	! if ((numxgridn.gt.maxxgridn).or.(numygridn.gt.maxygridn)) then
1225	! write(,) ' #### FLEXPART MODEL ERROR! DIMENSIONS OF ####'
1226	! write(,) ' #### OUTPUT NEST EXCEED MAXIMUM VALUES. ####'
1227	! write(,) ' #### CHANGE FILE $FLEXPART/options/OUTGRID. ####'
1228	! stop
1229	! endif
1230	xoutshiftn=xmet0-out_xm0n
1231	youtshiftn=ymet0-out_ym0n
1232	allocate(orooutn(0:numxgridn-1,0:numygridn-1) &
1233	,stat=stat)
1234	if (stat.ne.0) write(,)'ERROR: could not allocate orooutn'
1235	allocate(arean(0:numxgridn-1,0:numygridn-1) &
1236	,stat=stat)
1237	if (stat.ne.0) write(,)'ERROR: could not allocate arean'
1238	allocate(volumen(0:numxgridn-1,0:numygridn-1,numzgrid) &
1239	,stat=stat)
1240	if (stat.ne.0) write(,)'ERROR: could not allocate volumen'
1241
1242	endif
1243
1244	! jdf end read outgrid
1245	!*************************************************
1246	! jdf start read receptors
1247	!*************************************************
1248	!
1249	!*******************************************************************************
1250	! *
1251	! Note: This is the FLEXPART_WRF version of subroutine assignland. *
1252	! The computational grid is the WRF x-y grid rather than lat-lon. *
1253	! *
1254	! This routine reads the user specifications for the receptor points. *
1255	! *
1256	! Author: A. Stohl *
1257	! *
1258	! 1 August 1996 *
1259	! *
1260	! Oct 2005, R. Easter - change calc of receptorarea() *
1261	! Dec 2005, R. Easter - x/yrecptor values may be input either as *
1262	! degrees-latlon or grid-meters *
1263	! Changes to some read formats (wider fields). *
1264	! Changed names of "lon0" & "lat0" variables *
1265	! *
1266	!*******************************************************************************
1267	! *
1268	! Variables: *
1269	! receptorarea(maxreceptor) area of dxdy at location of receptor
1270	! receptorname(maxreceptor) names of receptors *
1271	! xreceptor,yreceptor coordinates of receptor points *
1272	! *
1273	! Constants: *
1274	! unitreceptor unit connected to file RECEPTORS *
1275	! *
1276	!*******************************************************************************
1277
1278	call skplin(1,unitpath)
1279	read(unitpath,*) numreceptor
1280	! For backward runs, do not allow receptor output. Thus, set number of receptors to zero
1281	!***************************************************************************************
1282
1283	if (ldirect.lt.0) then
1284	numreceptor=0
1285	endif
1286
1287	! Open the RECEPTORS file and read output grid specifications
1288	!************************************************************
1289	! Read the names and coordinates of the receptors
1290	!************************************************
1291
1292	if (numreceptor.gt.maxreceptor) then
1293	write(,) ' #### FLEXPART MODEL ERROR! TOO MANY RECEPTOR #### '
1294	write(,) ' #### POINTS ARE GIVEN. #### '
1295	write(,) ' #### MAXIMUM NUMBER IS ',maxreceptor,' #### '
1296	write(,) ' #### PLEASE MAKE CHANGES IN RECEPTORS #### '
1297	stop
1298	endif
1299	do j=1,numreceptor
1300	read(unitpath,'(4x,a16)') receptor
1301	read(unitpath,*) x
1302	read(unitpath,*) y
1303	receptorname(j)=receptor
1304
1305	write(*,'(/a,i5)') 'readreceptor diagnostics, j =', j
1306	write(*,'(a,1p,2e18.10)') 'x, y (in) ', x, y
1307	if (numpoint_option .eq. 1) then
1308	! In this case, the above inputs are the actual geographical lat/lon
1309	! Need to convert from lat/lon to grid-index coordinates
1310	receptor_lon(j) = x
1311	receptor_lat(j) = y
1312	call ll_to_xyindex_wrf( receptor_lon(j), receptor_lat(j), &
1313	xreceptor(j), yreceptor(j) )
1314	else
1315	! In this case, the above inputs are in grid-meters
1316	! Need to convert from grid-meters to grid-index coordinates, then to lat/lon
1317	xreceptor(j)=(x-xmet0)/dx
1318	yreceptor(j)=(y-ymet0)/dy
1319	call xyindex_to_ll_wrf( 0, xreceptor(j), yreceptor(j), &
1320	receptor_lon(j), receptor_lat(j) )
1321	write(*,'(f15.10,5x,a)') receptor_lon(j), 'receptor_lon'
1322	write(*,'(f15.10,5x,a)') receptor_lat(j), 'receptor_lat'
1323	end if
1324	write(*,'(a,1p,2e18.10)') 'x, yreceptor', &
1325	xreceptor(j), yreceptor(j)
1326
1327	! for FLEXPART_WRF, dx & dy are in meters,
1328	! receptorarea() appears to be area in m**2 of a mother grid cell
1329	! centers on x,y
1330	! xm=r_earthcos(ypi/180.)dx/180.pi
1331	! ym=r_earthdy/180.pi
1332	xm=dx
1333	ym=dy
1334	receptorarea(j)=xm*ym
1335	enddo
1336
1337	! jdf end read receptors
1338	!*************************************************
1339	! jdf start read species
1340	!*************************************************
1341	!
1342	!*******************************************************************************
1343	! *
1344	! This routine reads names and physical constants of chemical species/ *
1345	! radionuclides available with FLEXPART. *
1346	! *
1347	! Author: A. Stohl *
1348	! *
1349	! 11 July 1996 *
1350	! *
1351	!*******************************************************************************
1352	! *
1353	! Variables: *
1354	! decaytime(maxtable) half time for radiological decay *
1355	! specname(maxtable) names of chemical species, radionuclides *
1356	! wetscava, wetscavb Parameters for determining scavenging coefficient *
1357	! *
1358	! Constants: *
1359	! *
1360	!*******************************************************************************
1361
1362	! Open the SPECIES file and read species names and properties
1363	!************************************************************
1364
1365	read(unitpath,*)
1366	read(unitpath,*) numtable
1367	read(unitpath,*)
1368
1369	do i=1,numtable
1370	read(unitpath,21) specname(i),decaytime(i), &
1371	wetscava(i),wetscavb(i),drydiff(i),dryhenry(i),dryactiv(i), &
1372	partrho(i),partmean(i),partsig(i),dryvelo(i),weightmol(i), &
1373	ohreact(i),spec_ass(i),kao(i)
1374	! read(unitpath,21) specname(i),decay(i), &
1375	! weta(i),wetb(i),reldiff(i),henry(i),f0(i), &
1376	! density(i),dquer(i),dsigma(i),dryvel(i),weightmolar(i), &
1377	! ohreact(i),spec_ass(i),kao(i)
1378	pos_spec=i
1379
1380	if ((wetscava(i).gt.0.).and.(dryhenry(i).le.0.)) then
1381	if (partmean(i).le.0.) goto 996 ! no particle, no henry set
1382	endif
1383
1384
1385	if (spec_ass(i).gt.0) then
1386	spec_found=.FALSE.
1387	do j=1,pos_spec-1
1388	if (spec_ass(pos_spec).eq.specnum(j)) then
1389	spec_ass(pos_spec)=j
1390	spec_found=.TRUE.
1391	ASSSPEC=.TRUE.
1392	endif
1393	end do
1394	if (spec_found.eqv..FALSE.) then
1395	goto 997
1396	endif
1397	endif
1398
1399	if (partsig(i).eq.1.) partsig(i)=1.0001 ! avoid realing exception
1400	if (partsig(i).eq.0.) partsig(i)=1.0001 ! avoid realing exception
1401
1402	if ((drydiff(i).gt.0.).and.(partrho(i).gt.0.)) then
1403	write(,) '#### FLEXPART MODEL ERROR! SPECIES FORMAT ####'
1404	write(,) '#### IS CORRUPT. SPECIES CANNOT BE BOTH ####'
1405	write(,) '#### PARTICLE AND GAS. ####'
1406	stop
1407	endif
1408	enddo
1409
1410	21 format(4x,a10,f10.1,e11.1,f6.2,f7.1,e9.1,f5.1,e10.1,2e8.1,2f8.2, &
1411	e18.1,i18,f18.2)
1412
1413	! jdf end read species
1414	!*************************************************
1415	! jdf start read releases
1416	!*************************************************
1417	!
1418	!*******************************************************************************
1419	! *
1420	! This routine reads the release point specifications for the current *
1421	! model run. Several release points can be used at the same time. *
1422	! *
1423	! Note: This is the FLEXPART_WRF version of subroutine readreleases. *
1424	! The computational grid is the WRF x-y grid rather than lat-lon. *
1425	! *
1426	! Author: A. Stohl *
1427	! 18 May 1996 *
1428	! *
1429	! Update: 29 January 2001 *
1430	! Release altitude can be either in magl or masl *
1431	! *
1432	! Nov 2005, R. Easter - Do not adjust xpoint1 & 2 by +/-360 degrees *
1433	! Dec 2005, R. Easter - x/ypoint1/2 values may be input either as *
1434	! degrees-latlon or grid-meters *
1435	! *
1436	!*******************************************************************************
1437	! *
1438	! Variables: *
1439	! decay decay constant of species *
1440	! dquer [um] mean particle diameters *
1441	! dsigma e.g. dsigma=10 or dsigma=0.1 means that 68% of the mass *
1442	! are between 0.1dquer and 10dquer *
1443	! ireleasestart, ireleaseend [s] starting time and ending time of each release *
1444	! kindz 1: zpoint is in m agl, 2: zpoint is in m asl, 3: zpoint *
1445	! is in hPa *
1446	! npart number of particles to be released *
1447	! nspec number of species to be released *
1448	! density [kg/m3] density of the particles *
1449	! rm [s/m] Mesophyll resistance *
1450	! species name of species *
1451	! xmass total mass of each species *
1452	! xpoint1,ypoint1 geograf. coordinates of lower left corner of release area*
1453	! xpoint2,ypoint2 geograf. coordinates of upper right corner of release are*
1454	! weta, wetb parameters to determine the wet scavenging coefficient *
1455	! zpoint1,zpoint2 height range, over which release takes place *
1456	! *
1457	!*******************************************************************************
1458
1459	DEP=.false.
1460	DRYDEP=.false.
1461	WETDEP=.false.
1462	do i=1,maxspec
1463	DRYDEPSPEC(i)=.false.
1464	enddo
1465
1466	! Open the releases file and read user options
1467	!*********************************************
1468
1469	! Check the format of the RELEASES file (either in free format,
1470	! or using a formatted mask)
1471	! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
1472	!**************************************************************************
1473	call skplin(1,unitpath)
1474
1475	! Read the number of species and the link to the species information table
1476	! Assign species-specific parameters needed for physical processes
1477	!*************************************************************************
1478
1479	read(unitpath,*) nspec
1480	read(unitpath,*) emitvar
1481	if (nspec.gt.maxspec) then
1482	write(,) '#####################################################'
1483	write(,) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
1484	write(,) '#### ####'
1485	write(,) '#### ERROR - MAXIMUM NUMBER OF EMITTED SPECIES IS####'
1486	write(,) '#### TOO LARGE. PLEASE REDUCE NUMBER OF SPECIES. ####'
1487	write(,) '#####################################################'
1488	stop
1489	endif
1490	do i=1,nspec
1491	read(unitpath,*) link(i)
1492	species(i)=specname(link(i))
1493
1494	! For backward runs, only 1 species is allowed
1495	!*********************************************
1496
1497	if ((ldirect.lt.0).and.(nspec.gt.1)) then
1498	write(,) '#####################################################'
1499	write(,) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
1500	write(,) '#### FOR BACKWARD RUNS, ONLY 1 SPECIES IS ALLOWED####'
1501	write(,) '#####################################################'
1502	stop
1503	endif
1504
1505	! Molecular weight
1506	!*****************
1507
1508	weightmolar(i)=weightmol(link(i))
1509	if (((iout.eq.2).or.(iout.eq.3)).and. &
1510	(weightmolar(i).lt.0.)) then
1511	write(,) 'For mixing ratio output, valid molar weight'
1512	write(,) 'must be specified for all simulated species.'
1513	write(,) 'Check table SPECIES or choose concentration'
1514	write(,) 'output instead if molar weight is not known.'
1515	stop
1516	endif
1517
1518
1519	! Radioactive decay
1520	!******************
1521
1522	decay(i)=0.693147/decaytime(link(i)) !conversion half life to decay constant
1523
1524	! Wet deposition
1525	!***************
1526
1527	weta(i)=wetscava(link(i))
1528	wetb(i)=wetscavb(link(i))
1529
1530	! Dry deposition of gases
1531	!************************
1532
1533	reldiff(i)=drydiff(link(i)) ! Diffusivity rel. to H20
1534	henry(i)=dryhenry(link(i)) ! Henry constant
1535	f0(i)=dryactiv(link(i)) ! activity
1536	if (reldiff(i).gt.0.) &
1537	rm(i)=1./(henry(i)/3000.+100.*f0(i)) ! mesophyll resistance
1538
1539	! Dry deposition of particles
1540	!****************************
1541
1542	vsetaver(i)=0.
1543	density(i)=partrho(link(i)) ! Particle density
1544	dquer(i)=partmean(link(i))*1000000. ! Conversion m to um
1545	dsigma(i)=partsig(link(i))
1546	if (density(i).gt.0.) then ! Additional parameters
1547	! call part0(dquer(i),dsigma(i),density(i),fracth,schmih,vsh)
1548	call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh)
1549
1550	do j=1,ni
1551	fract(i,j)=fracth(j)
1552	schmi(i,j)=schmih(j)
1553	vset(i,j)=vsh(j)
1554	cunningham(i)=cunningham(i)+cun*fract(i,j)
1555	vsetaver(i)=vsetaver(i)-vset(i,j)*fract(i,j)
1556	enddo
1557	endif
1558
1559	! Dry deposition for constant deposition velocity
1560	!************************************************
1561
1562	dryvel(i)=dryvelo(link(i))*0.01 ! conversion to m/s
1563
1564	if (weta(i).gt.0.) WETDEP=.true.
1565	if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or. &
1566	(dryvel(i).gt.0.)) then
1567	DRYDEP=.true.
1568	DRYDEPSPEC(i)=.true.
1569	endif
1570
1571
1572	! Read in daily and day-of-week variation of emissions, if available
1573	!*******************************************************************
1574
1575	do j=1,24 ! initialize everything to no variation
1576	area_hour(i,j)=1.
1577	point_hour(i,j)=1.
1578	enddo
1579	do j=1,7
1580	area_dow(i,j)=1.
1581	point_dow(i,j)=1.
1582	enddo
1583
1584	! write(aspecnumb,'(i3.3)') link(i)
1585	! open(unitemissvar,file=path(1)(1:len(1))//'EMISSION_VARIATION_'
1586	! + //aspecnumb//'.dat',status='old',err=35)
1587	! read(unitemissvar,*)
1588	if(emitvar.eq.1) then
1589	do j=1,24 ! 24 hours, starting with 0-1 local time
1590	read(unitpath,*) ihour,area_hour(i,j),point_hour(i,j)
1591	enddo
1592	! read(unitemissvar,*)
1593	do j=1,7 ! 7 days of the week, starting with Monday
1594	read(unitpath,*) idow,area_dow(i,j),point_dow(i,j)
1595	enddo
1596	! close(unitemissvar)
1597	endif
1598
1599	!35 continue
1600	enddo
1601	if (WETDEP.or.DRYDEP) DEP=.true.
1602
1603
1604	! Read specifications for each release point
1605	!*******************************************
1606
1607	numpoint=0
1608	numpartmax=0
1609	releaserate=0.
1610	releaserate2=0.
1611	read(unitpath,*) numpoint
1612
1613	! numpoint2=numpoint+500
1614	numpoint2=numpoint+0
1615	allocate(ireleasestart(numpoint2) &
1616	,stat=stat)
1617	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
1618	allocate(ireleaseend(numpoint2) &
1619	,stat=stat)
1620	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
1621	allocate(xpoint1(numpoint2) &
1622	,stat=stat)
1623	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
1624	allocate(xpoint12(numpoint2) &
1625	,stat=stat)
1626	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
1627	allocate(xpoint22(numpoint2) &
1628	,stat=stat)
1629	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
1630	allocate(ypoint12(numpoint2) &
1631	,stat=stat)
1632	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
1633	allocate(ypoint22(numpoint2) &
1634	,stat=stat)
1635	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
1636	allocate(releases_swlon(numpoint2) &
1637	,stat=stat)
1638	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
1639	allocate(releases_swlat(numpoint2) &
1640	,stat=stat)
1641	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
1642	allocate(releases_nelon(numpoint2) &
1643	,stat=stat)
1644	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
1645	allocate(releases_nelat(numpoint2) &
1646	,stat=stat)
1647	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
1648	allocate(xpoint2(numpoint2) &
1649	,stat=stat)
1650	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
1651	allocate(ypoint1(numpoint2) &
1652	,stat=stat)
1653	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
1654	allocate(ypoint2(numpoint2) &
1655	,stat=stat)
1656	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
1657	allocate(zpoint1(numpoint2) &
1658	,stat=stat)
1659	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
1660	allocate(zpoint2(numpoint2) &
1661	,stat=stat)
1662	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
1663	allocate(kindz(numpoint2) &
1664	,stat=stat)
1665	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
1666	allocate(xmass(numpoint2,maxspec) &
1667	,stat=stat)
1668	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
1669	allocate(rho_rel(numpoint2) &
1670	,stat=stat)
1671	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
1672	allocate(npart(numpoint2) &
1673	,stat=stat)
1674	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
1675	! print*,'allocate xmassa',numpoint
1676	allocate(xmasssave(numpoint2) &
1677	,stat=stat)
1678	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
1679
1680	if (option_verbose.eq. 1) then
1681	write (,) ' Releasepoints allocated: ', numpoint
1682	endif
1683	do i=1,numpoint
1684	xmasssave(i)=0.
1685	end do
1686
1687
1688
1689	! if (numpoint.gt.maxpoint) goto 997
1690	do j=1,numpoint
1691	read(unitpath,*) id1,it1
1692	read(unitpath,*) id2,it2
1693	read(unitpath,*) xpoint1(j)
1694	read(unitpath,*) ypoint1(j)
1695	read(unitpath,*) xpoint2(j)
1696	read(unitpath,*) ypoint2(j)
1697	read(unitpath,*) kindz(j)
1698	read(unitpath,*) zpoint1(j)
1699	read(unitpath,*) zpoint2(j)
1700	read(unitpath,*) npart(j)
1701	do i=1,nspec
1702	read(unitpath,*) xmass(j,i)
1703	enddo
1704	nparttot=nparttot+npart(j)
1705	if (j.le.2000) then
1706	read(unitreleases,'(a40)',err=998) compoint(j)(1:40)
1707	else
1708	read(unitreleases,'(a40)',err=998) compoint(2001)(1:40)
1709	endif
1710
1711	! read(unitpath,'(a20)',err=998) compoint(j)(1:20)
1712
1713	if(option_verbose.ge.1) print*,'release location=',j
1714	! write(*,'(/a,i7)') 'readreleases diagnostics - numpoint = ', j
1715	! write(*,'(a,1p,2e18.10)') 'x, ypoint1 (in) ', &
1716	! xpoint1(j), ypoint1(j)
1717	! write(*,'(a,1p,2e18.10)') 'x, ypoint2 (in) ', &
1718	! xpoint2(j), ypoint2(j)
1719	! JB
1720	! In this case, the above inputs are the actual geographical lat/lon
1721	! of the southwest & northeast corners of the release area
1722	! Need to convert from lat/lon to grid-meters
1723	! if (outgrid_option.eq.1) then !regular
1724	if (numpoint_option.eq.1) then !regular
1725	xpoint12(j)=xpoint1(j)
1726	xpoint22(j)=xpoint2(j)
1727	ypoint12(j)=ypoint1(j)
1728	ypoint22(j)=ypoint2(j)
1729
1730
1731	releases_swlon(j) = xpoint1(j)
1732	releases_swlat(j) = ypoint1(j)
1733	releases_nelon(j) = xpoint2(j)
1734	releases_nelat(j) = ypoint2(j)
1735	call ll_to_xymeter_wrf( releases_swlon(j), releases_swlat(j), &
1736	xpoint1(j), ypoint1(j) )
1737	call ll_to_xymeter_wrf( releases_nelon(j), releases_nelat(j), &
1738	xpoint2(j), ypoint2(j) )
1739
1740	else
1741	! write(*,'(a,1p,2e18.10)') 'x, ypoint1 ',
1742	! & xpoint1(j), ypoint1(j)
1743	! write(*,'(a,1p,2e18.10)') 'x, ypoint2 ',
1744	! & xpoint2(j), ypoint2(j)
1745	! else
1746	! In this case, the above inputs are in grid-meters
1747	! Need to convert from grid-meters to lat/lon
1748	call xymeter_to_ll_wrf( xpoint1(j), ypoint1(j), &
1749	releases_swlon(j), releases_swlat(j) )
1750	call xymeter_to_ll_wrf( xpoint2(j), ypoint2(j), &
1751	releases_nelon(j), releases_nelat(j) )
1752	endif
1753	! write(*,'(f15.10,5x,a)') releases_swlon(j), 'releases_swlon'
1754	! write(*,'(f15.10,5x,a)') releases_swlat(j), 'releases_swlat'
1755	! write(*,'(f15.10,5x,a)') releases_nelon(j), 'releases_nelon'
1756	! write(*,'(f15.10,5x,a)') releases_nelat(j), 'releases_nelat'
1757	! end if
1758
1759	! If a release point contains no particles, stop and issue error message
1760	!***********************************************************************
1761
1762	if (npart(j).eq.0) then
1763	write(,) 'FLEXPART MODEL ERROR'
1764	write(,) 'RELEASES file is corrupt.'
1765	write(,) 'At least for one release point, there are zero'
1766	write(,) 'particles released. Make changes to RELEASES.'
1767	stop
1768	endif
1769
1770	! Check whether x coordinates of release point are within model domain
1771	!*********************************************************************
1772
1773	! FLEXPART_WRF - x & y coords are in meters, so the following lines
1774	! (which adjust longitude by +/-360 degrees) are not needed
1775	!
1776	! if (xpoint1(numpoint).lt.xlon0)
1777	! + xpoint1(numpoint)=xpoint1(numpoint)+360.
1778	! if (xpoint1(numpoint).gt.xlon0+(nxmin1)*dx)
1779	! + xpoint1(numpoint)=xpoint1(numpoint)-360.
1780	! if (xpoint2(numpoint).lt.xlon0)
1781	! + xpoint2(numpoint)=xpoint2(numpoint)+360.
1782	! if (xpoint2(numpoint).gt.xlon0+(nxmin1)*dx)
1783	! + xpoint2(numpoint)=xpoint2(numpoint)-360.
1784
1785	! Determine relative beginning and ending times of particle release
1786	!******************************************************************
1787
1788	jul1=juldate(id1,it1)
1789	jul2=juldate(id2,it2)
1790	if (jul1.gt.jul2) then
1791	write(,) 'FLEXPART MODEL ERROR'
1792	write(,) 'Release stops before it begins.'
1793	write(,) 'Make changes to file RELEASES.'
1794	stop
1795	endif
1796	if (mdomainfill.eq.0) then ! no domain filling
1797	if (ldirect.eq.1) then
1798	if ((jul1.lt.bdate).or.(jul2.gt.edate)) then
1799	write(,) 'FLEXPART MODEL ERROR'
1800	write(,) 'Release starts before simulation begins or ends'
1801	write(,) 'after simulation stops.'
1802	write(,) 'Make files COMMAND and RELEASES consistent.'
1803	stop
1804	endif
1805	ireleasestart(j)=int((jul1-bdate)*86400.)
1806	ireleaseend(j)=int((jul2-bdate)*86400.)
1807	else if (ldirect.eq.-1) then
1808	if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
1809	write(,) 'FLEXPART MODEL ERROR'
1810	write(,) 'Release starts before simulation begins or ends'
1811	write(,) 'after simulation stops.'
1812	write(,) 'Make files COMMAND and RELEASES consistent.'
1813	stop
1814	endif
1815	ireleasestart(j)=int((jul1-bdate)*86400.)
1816	ireleaseend(j)=int((jul2-bdate)*86400.)
1817	endif
1818	endif
1819
1820
1821	! Check, whether the total number of particles may exceed totally allowed
1822	! number of particles at some time during the simulation
1823	!************************************************************************
1824
1825	! Determine the release rate (particles per second) and total number
1826	! of particles released during the simulation
1827	!*******************************************************************
1828
1829	if (ireleasestart(j).ne.ireleaseend(j)) then
1830	releaserate=releaserate+real(npart(j))/ &
1831	real(ireleaseend(j)-ireleasestart(j))
1832	! print*,'release',ireleaseend(j),ireleasestart(j)
1833	else
1834	releaserate=99999999.
1835	endif
1836	if (ireleaseend(j)-ireleasestart(j).lt.lage(nageclass)) then
1837	releaserate2=releaserate2+real(npart(j))
1838	else
1839	releaserate2=releaserate2+real(npart(j))*real(lage(nageclass))/ &
1840	real(ireleaseend(j)-ireleasestart(j))
1841	endif
1842	numpartmax=numpartmax+npart(j)
1843
1844	if (ioutputforeachrelease.eq.1) then
1845	maxpointspec_act=numpoint
1846	else
1847	maxpointspec_act=1
1848	endif
1849
1850	enddo
1851
1852	maxpart=int(releaserate2*1.1)
1853	! print*,'maxpart',maxpart,releaserate2,nparttot
1854	! maxpart=40000000
1855
1856	! print*,'numpoint',numpoint
1857	! print*,'NPARTTOT',nparttot,lage(4),lage(5),lage(6),nageclass
1858	! print,'maxpart',maxpart,real(lage(nageclass)),releaserate,releaserate2,real(lage(nageclass))1.1*releaserate
1859	! nparttot2=nparttot+500000
1860	! nparttot2=nparttot
1861
1862	nparttot2=maxpart
1863	allocate(xmass1(nparttot2,nspec) &
1864	,stat=stat)
1865	! allocate(drydep1(nparttot2,nspec) &
1866	! ,stat=stat)
1867	if (stat.ne.0) write(,)'ERROR: could not allocate xmass1'
1868	allocate(itra1(nparttot2) &
1869	,stat=stat)
1870	if (stat.ne.0) write(,)'ERROR: could not allocate itra1'
1871	allocate(npoint(nparttot2) &
1872	,stat=stat)
1873	if (stat.ne.0) write(,)'ERROR: could not allocate npoint'
1874	allocate(nclass(nparttot2) &
1875	,stat=stat)
1876	if (stat.ne.0) write(,)'ERROR: could not allocate nclass'
1877	allocate(idt(nparttot2) &
1878	,stat=stat)
1879	if (stat.ne.0) write(,)'ERROR: could not allocate idt '
1880	allocate(itramem(nparttot2) &
1881	,stat=stat)
1882	if (stat.ne.0) write(,)'ERROR: could not allocate itrame'
1883	allocate(itrasplit(nparttot2) &
1884	,stat=stat)
1885	if (stat.ne.0) write(,)'ERROR: could not allocate itrasp'
1886	allocate(xtra1(nparttot2) &
1887	,stat=stat)
1888	if (stat.ne.0) write(,)'ERROR: could not allocate xtra1 '
1889	allocate(ytra1(nparttot2) &
1890	,stat=stat)
1891	if (stat.ne.0) write(,)'ERROR: could not allocate ytra1 '
1892	allocate(ztra1(nparttot2) &
1893	,stat=stat)
1894	if (stat.ne.0) write(,)'ERROR: could not allocate ztra1 '
1895
1896	do j=1,nparttot2
1897	itra1(j)=-999999999
1898	enddo
1899
1900
1901	if (releaserate.gt. &
1902	0.99*real(maxpart)/real(lage(nageclass))) then
1903	! if (numpartmax.gt.maxpart) then
1904	! write(,) '#####################################################'
1905	! write(,) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
1906	! write(,) '#### ####'
1907	! write(,) '####WARNING - TOTAL NUMBER OF PARTICLES SPECIFIED####'
1908	! write(,) '#### IN FILE "RELEASES" MAY AT SOME POINT DURING ####'
1909	! write(,) '#### THE SIMULATION EXCEED THE MAXIMUM ALLOWED ####'
1910	! write(,) '#### NUMBER (MAXPART).IF RELEASES DO NOT OVERLAP,####'
1911	! write(,) '#### FLEXPART CAN POSSIBLY COMPLETE SUCCESSFULLY.####'
1912	! write(,) '#### HOWEVER, FLEXPART MAY HAVE TO STOP ####'
1913	! write(,) '#### AT SOME TIME DURING THE SIMULATION. PLEASE ####'
1914	! write(,) '#### MAKE SURE THAT YOUR SETTINGS ARE CORRECT. ####'
1915	! write(,) '#####################################################'
1916	! write(,) 'Maximum release rate may be: ',releaserate, &
1917	! ' particles per second'
1918	! write(,) 'Maximum allowed release rate is: ', &
1919	! real(maxpart)/real(lage(nageclass)),' particles per second'
1920	! write(,) &
1921	! 'Total number of particles released during the simulation is: ', &
1922	! numpartmax
1923	! write(,) 'Maximum allowed number of particles is: ',maxpart
1924	! endif
1925	endif
1926
1927
1928	! Make a consistency check, whether the forward/backward switch is correctly set
1929	!*******************************************************************************
1930
1931	! if (ldirect.eq.1) then
1932	! if (maxpointspec.lt.nspec) then
1933	! write(,) '#####################################################'
1934	! write(,) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
1935	! write(,) '#### ####'
1936	! write(,) '#### ERROR - PARAMETER MAXPOINTSPEC IS NOT ####'
1937	! write(,) '#### CORRECTLY SET FOR A FORWARD SIMULATION. ####'
1938	! write(,) '#### CHANGE APPROPRIATELY IN FILE INCLUDEPAR. ####'
1939	! write(,) '#####################################################'
1940	! endif
1941	! else
1942	! if (maxpointspec.lt.numpoint) then
1943	! write(,) '#####################################################'
1944	! write(,) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
1945	! write(,) '#### ####'
1946	! write(,) '#### ERROR - PARAMETER MAXPOINTSPEC IS NOT ####'
1947	! write(,) '#### CORRECTLY SET FOR A BACKWARD SIMULATION. ####'
1948	! write(,) '#### CHANGE APPROPRIATELY IN FILE INCLUDEPAR. ####'
1949	! write(,) '#####################################################'
1950	! endif
1951	! endif
1952
1953	return
1954
1955	! jdf end read releases
1956	!*************************************************
1957
1958	997 write(,) '#####################################################'
1959	write(,) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
1960	write(,) '#### ####'
1961	write(,) '#### ERROR - NUMBER OF RELEASE POINTS SPECIFIED ####'
1962	write(,) '#### IN "RELEASES" EXCEEDS THE MAXIMUM ####'
1963	write(,) '#### ALLOWED NUMBER. ####'
1964	write(,) '#####################################################'
1965	stop
1966
1967
1968	998 write(,) '#####################################################'
1969	write(,) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
1970	write(,) '#### ####'
1971	write(,) '#### FATAL ERROR - "RELEASES" IS ####'
1972	write(,) '#### CORRUPT. PLEASE CHECK YOUR INPUTS FOR ####'
1973	write(,) '#### MISTAKES OR GET A NEW "RELEASES"- ####'
1974	write(,) '#####################################################'
1975	stop
1976
1977
1978	999 write(,) '#####################################################'
1979	write(,) ' FLEXPART MODEL SUBROUTINE READRELEASES: '
1980	write(,)
1981	write(,) 'FATAL ERROR - FILE CONTAINING PARTICLE RELEASE POINTS'
1982	write(,) 'POINTS IS NOT AVAILABLE OR YOU ARE NOT'
1983	write(,) 'PERMITTED FOR ANY ACCESS'
1984	write(,) '#####################################################'
1985	stop
1986
1987	800 write(,) ' #### TRAJECTORY MODEL ERROR! ERROR WHILE #### '
1988	write(,) ' #### READING FILE PATHNAMES. #### '
1989	stop
1990
1991	801 write(,) '#### TRAJECTORY MODEL ERROR! FILE '// inputname
1992	write(,) '#### CANNOT BE OPENED IN THE CURRENT WORKING #### '
1993	write(,) '#### DIRECTORY. #### '
1994	stop
1995
1996	803 write(,) ' #### FLEXPART MODEL ERROR! FILE #### '
1997	write(,'(a)') ' '//path(numpath+2(k-1)+2) &
1998	(1:length(numpath+2*(k-1)+2))
1999	write(,) ' #### CANNOT BE OPENED #### '
2000	stop
2001
2002	804 write(,) ' #### FLEXPART MODEL ERROR! FILE #### '
2003	write(*,'(a)') ' '//path(3)(1:length(3))
2004	write(,) ' #### CANNOT BE OPENED #### '
2005	stop
2006	996 write(,) '#####################################################'
2007	write(,) '#### FLEXPART MODEL ERROR! #### '
2008	write(,) '#### WET DEPOSITION SWITCHED ON, BUT NO HENRYS #### '
2009	write(,) '#### CONSTANT IS SET ####'
2010	write(,) '#### PLEASE MODIFY SPECIES DESCR. #### '
2011	write(,) '#####################################################'
2012	stop
2013
2014	end

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

Download in other formats:

Original Format