Context Navigation

concoutput.f90 @ 2eefa58

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

Last change on this file since 2eefa58 was 2eefa58, checked in by Espen Sollum ATMOS <eso@…>, 5 years ago
Added Ronas changes for inversion output
Property mode set to `100644`
File size: 26.6 KB

Line
1	!**********************************************************************
2	! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
3	! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
4	! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
5	! *
6	! This file is part of FLEXPART. *
7	! *
8	! FLEXPART is free software: you can redistribute it and/or modify *
9	! it under the terms of the GNU General Public License as published by*
10	! the Free Software Foundation, either version 3 of the License, or *
11	! (at your option) any later version. *
12	! *
13	! FLEXPART is distributed in the hope that it will be useful, *
14	! but WITHOUT ANY WARRANTY; without even the implied warranty of *
15	! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
16	! GNU General Public License for more details. *
17	! *
18	! You should have received a copy of the GNU General Public License *
19	! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
20	!**********************************************************************
21
22	subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
23	drygridtotalunc)
24	! i i o o
25	! o
26	!*****************************************************************************
27	! *
28	! Output of the concentration grid and the receptor concentrations. *
29	! *
30	! Author: A. Stohl *
31	! *
32	! 24 May 1995 *
33	! *
34	! 13 April 1999, Major update: if output size is smaller, dump output *
35	! in sparse matrix format; additional output of *
36	! uncertainty *
37	! *
38	! 05 April 2000, Major update: output of age classes; output for backward*
39	! runs is time spent in grid cell times total mass of *
40	! species. *
41	! *
42	! 17 February 2002, Appropriate dimensions for backward and forward runs *
43	! are now specified in file par_mod *
44	! *
45	! June 2006, write grid in sparse matrix with a single write command *
46	! in order to save disk space *
47	! *
48	! 2008 new sparse matrix format *
49	! *
50	!*****************************************************************************
51	! *
52	! Variables: *
53	! outnum number of samples *
54	! ncells number of cells with non-zero concentrations *
55	! sparse .true. if in sparse matrix format, else .false. *
56	! tot_mu 1 for forward, initial mass mixing ration for backw. runs *
57	! *
58	!*****************************************************************************
59
60	use unc_mod
61	use point_mod
62	use outg_mod
63	use par_mod
64	use com_mod
65	use mean_mod
66
67	implicit none
68
69	real(kind=dp) :: jul
70	integer :: itime,i,ix,jy,kz,ks,kp,l,iix,jjy,kzz,nage,jjjjmmdd,ihmmss
71	integer :: sp_count_i,sp_count_r
72	real :: sp_fact
73	real :: outnum,densityoutrecept(maxreceptor),xl,yl
74	! RLT
75	real :: densitydryrecept(maxreceptor)
76	real :: factor_dryrecept(maxreceptor)
77
78	!real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
79	! +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,
80	! + maxageclass)
81	!real wetgrid(0:numxgrid-1,0:numygrid-1,maxspec,maxpointspec_act,
82	! + maxageclass)
83	!real drygrid(0:numxgrid-1,0:numygrid-1,maxspec,
84	! + maxpointspec_act,maxageclass)
85	!real gridsigma(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,
86	! + maxpointspec_act,maxageclass),
87	! + drygridsigma(0:numxgrid-1,0:numygrid-1,maxspec,
88	! + maxpointspec_act,maxageclass),
89	! + wetgridsigma(0:numxgrid-1,0:numygrid-1,maxspec,
90	! + maxpointspec_act,maxageclass)
91	!real factor(0:numxgrid-1,0:numygrid-1,numzgrid)
92	!real sparse_dump_r(numxgridnumygridnumzgrid)
93	!integer sparse_dump_i(numxgridnumygridnumzgrid)
94
95	!real sparse_dump_u(numxgridnumygridnumzgrid)
96	real(dep_prec) :: auxgrid(nclassunc)
97	real(sp) :: gridtotal,gridsigmatotal,gridtotalunc
98	real(dep_prec) :: wetgridtotal,wetgridsigmatotal,wetgridtotalunc
99	real(dep_prec) :: drygridtotal,drygridsigmatotal,drygridtotalunc
100	real :: halfheight,dz,dz1,dz2,tot_mu(maxspec,maxpointspec_act)
101	real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
102	real,parameter :: weightair=28.97
103	logical :: sp_zer
104	logical,save :: init=.true.
105	character :: adate8,atime6
106	character(len=3) :: anspec
107	integer :: mind
108	! mind eso:added to ensure identical results between 2&3-fields versions
109	character(LEN=8),save :: file_stat='REPLACE'
110	logical :: ldates_file
111	logical :: lexist
112	integer :: ierr
113	character(LEN=100) :: dates_char
114
115	! Determine current calendar date, needed for the file name
116	!**********************************************************
117
118	jul=bdate+real(itime,kind=dp)/86400._dp
119	call caldate(jul,jjjjmmdd,ihmmss)
120	write(adate,'(i8.8)') jjjjmmdd
121	write(atime,'(i6.6)') ihmmss
122
123	! Overwrite existing dates file on first call, later append to it
124	! This fixes a bug where the dates file kept growing across multiple runs
125
126	! If 'dates' file exists in output directory, make a backup
127	inquire(file=path(2)(1:length(2))//'dates', exist=ldates_file)
128	if (ldates_file.and.init) then
129	open(unit=unitdates, file=path(2)(1:length(2))//'dates',form='formatted', &
130	&access='sequential', status='old', action='read', iostat=ierr)
131	open(unit=unittmp, file=path(2)(1:length(2))//'dates.bak', access='sequential', &
132	&status='replace', action='write', form='formatted', iostat=ierr)
133	do while (.true.)
134	read(unitdates, '(a)', iostat=ierr) dates_char
135	if (ierr.ne.0) exit
136	write(unit=unittmp, fmt='(a)', iostat=ierr, advance='yes') trim(dates_char)
137	end do
138	close(unit=unitdates)
139	close(unit=unittmp)
140	end if
141
142	open(unitdates,file=path(2)(1:length(2))//'dates', ACCESS='APPEND', STATUS=file_stat)
143	write(unitdates,'(a)') adate//atime
144	close(unitdates)
145
146	! Overwrite existing dates file on first call, later append to it
147	! This fixes a bug where the dates file kept growing across multiple runs
148	IF (init) THEN
149	file_stat='OLD'
150	init=.false.
151	END IF
152
153
154	! For forward simulations, output fields have dimension MAXSPEC,
155	! for backward simulations, output fields have dimension MAXPOINT.
156	! Thus, make loops either about nspec, or about numpoint
157	!*****************************************************************
158
159
160	if (ldirect.eq.1) then
161	do ks=1,nspec
162	do kp=1,maxpointspec_act
163	tot_mu(ks,kp)=1
164	end do
165	end do
166	else
167	do ks=1,nspec
168	do kp=1,maxpointspec_act
169	tot_mu(ks,kp)=xmass(kp,ks)
170	end do
171	end do
172	endif
173
174
175	!*******************************************************************
176	! Compute air density: sufficiently accurate to take it
177	! from coarse grid at some time
178	! Determine center altitude of output layer, and interpolate density
179	! data to that altitude
180	!*******************************************************************
181
182	mind=memind(2)
183	do kz=1,numzgrid
184	if (kz.eq.1) then
185	halfheight=outheight(1)/2.
186	else
187	halfheight=(outheight(kz)+outheight(kz-1))/2.
188	endif
189	do kzz=2,nz
190	if ((height(kzz-1).lt.halfheight).and. &
191	(height(kzz).gt.halfheight)) goto 46
192	end do
193	46 kzz=max(min(kzz,nz),2)
194	dz1=halfheight-height(kzz-1)
195	dz2=height(kzz)-halfheight
196	dz=dz1+dz2
197	do jy=0,numygrid-1
198	do ix=0,numxgrid-1
199	xl=outlon0+real(ix)*dxout
200	yl=outlat0+real(jy)*dyout
201	xl=(xl-xlon0)/dx
202	yl=(yl-ylat0)/dy !v9.1.1
203	iix=max(min(nint(xl),nxmin1),0)
204	jjy=max(min(nint(yl),nymin1),0)
205	! densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
206	! rho(iix,jjy,kzz-1,2)*dz2)/dz
207	densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
208	rho(iix,jjy,kzz-1,mind)*dz2)/dz
209	! RLT
210	densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,mind)*dz1+ &
211	rho_dry(iix,jjy,kzz-1,mind)*dz2)/dz
212	end do
213	end do
214	end do
215
216	do i=1,numreceptor
217	xl=xreceptor(i)
218	yl=yreceptor(i)
219	iix=max(min(nint(xl),nxmin1),0)
220	jjy=max(min(nint(yl),nymin1),0)
221	!densityoutrecept(i)=rho(iix,jjy,1,2)
222	densityoutrecept(i)=rho(iix,jjy,1,mind)
223	! RLT
224	densitydryrecept(i)=rho_dry(iix,jjy,1,mind)
225	end do
226
227	! RLT
228	! conversion factor for output relative to dry air
229	factor_drygrid=densityoutgrid/densitydrygrid
230	factor_dryrecept=densityoutrecept/densitydryrecept
231
232	! Output is different for forward and backward simulations
233	do kz=1,numzgrid
234	do jy=0,numygrid-1
235	do ix=0,numxgrid-1
236	if (ldirect.eq.1) then
237	factor3d(ix,jy,kz)=1.e12/volume(ix,jy,kz)/outnum
238	else
239	factor3d(ix,jy,kz)=real(abs(loutaver))/outnum
240	endif
241	end do
242	end do
243	end do
244
245	!*********************************************************************
246	! Determine the standard deviation of the mean concentration or mixing
247	! ratio (uncertainty of the output) and the dry and wet deposition
248	!*********************************************************************
249
250	gridtotal=0.
251	gridsigmatotal=0.
252	gridtotalunc=0.
253	wetgridtotal=0.
254	wetgridsigmatotal=0.
255	wetgridtotalunc=0.
256	drygridtotal=0.
257	drygridsigmatotal=0.
258	drygridtotalunc=0.
259
260	do ks=1,nspec
261
262	write(anspec,'(i3.3)') ks
263
264	if (DRYBKDEP.or.WETBKDEP) then !scavdep output
265	if (DRYBKDEP) &
266	open(unitoutgrid,file=path(2)(1:length(2))//'grid_drydep_'//adate// &
267	atime//'_'//anspec,form='unformatted')
268	if (WETBKDEP) &
269	open(unitoutgrid,file=path(2)(1:length(2))//'grid_wetdep_'//adate// &
270	atime//'_'//anspec,form='unformatted')
271	write(unitoutgrid) itime
272	else
273	if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
274	if (ldirect.eq.1) then
275	open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//adate// &
276	atime//'_'//anspec,form='unformatted')
277	else
278	open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// &
279	atime//'_'//anspec,form='unformatted')
280	endif
281	write(unitoutgrid) itime
282	endif
283	if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio
284	open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// &
285	atime//'_'//anspec,form='unformatted')
286	write(unitoutgridppt) itime
287	endif
288	endif ! if deposition output
289
290	do kp=1,maxpointspec_act
291	do nage=1,nageclass
292
293	do jy=0,numygrid-1
294	do ix=0,numxgrid-1
295
296	! WET DEPOSITION
297	if ((WETDEP).and.(ldirect.gt.0)) then
298	do l=1,nclassunc
299	auxgrid(l)=wetgridunc(ix,jy,ks,kp,l,nage)
300	end do
301	call mean(auxgrid,wetgrid(ix,jy), &
302	wetgridsigma(ix,jy),nclassunc)
303	! Multiply by number of classes to get total concentration
304	wetgrid(ix,jy)=wetgrid(ix,jy) &
305	*nclassunc
306	wetgridtotal=wetgridtotal+wetgrid(ix,jy)
307	! Calculate standard deviation of the mean
308	wetgridsigma(ix,jy)= &
309	wetgridsigma(ix,jy)* &
310	sqrt(real(nclassunc))
311	wetgridsigmatotal=wetgridsigmatotal+ &
312	wetgridsigma(ix,jy)
313	endif
314
315	! DRY DEPOSITION
316	if ((DRYDEP).and.(ldirect.gt.0)) then
317	do l=1,nclassunc
318	auxgrid(l)=drygridunc(ix,jy,ks,kp,l,nage)
319	end do
320	call mean(auxgrid,drygrid(ix,jy), &
321	drygridsigma(ix,jy),nclassunc)
322	! Multiply by number of classes to get total concentration
323	drygrid(ix,jy)=drygrid(ix,jy)* &
324	nclassunc
325	drygridtotal=drygridtotal+drygrid(ix,jy)
326	! Calculate standard deviation of the mean
327	drygridsigma(ix,jy)= &
328	drygridsigma(ix,jy)* &
329	sqrt(real(nclassunc))
330	drygridsigmatotal=drygridsigmatotal+ &
331	drygridsigma(ix,jy)
332	endif
333
334	! CONCENTRATION OR MIXING RATIO
335	do kz=1,numzgrid
336	do l=1,nclassunc
337	auxgrid(l)=gridunc(ix,jy,kz,ks,kp,l,nage)
338	end do
339	call mean(auxgrid,grid(ix,jy,kz), &
340	gridsigma(ix,jy,kz),nclassunc)
341	! Multiply by number of classes to get total concentration
342	grid(ix,jy,kz)= &
343	grid(ix,jy,kz)*nclassunc
344	gridtotal=gridtotal+grid(ix,jy,kz)
345	! Calculate standard deviation of the mean
346	gridsigma(ix,jy,kz)= &
347	gridsigma(ix,jy,kz)* &
348	sqrt(real(nclassunc))
349	gridsigmatotal=gridsigmatotal+ &
350	gridsigma(ix,jy,kz)
351	end do
352	end do
353	end do
354
355
356
357
358	!*******************************************************************
359	! Generate output: may be in concentration (ng/m3) or in mixing
360	! ratio (ppt) or both
361	! Output the position and the values alternated multiplied by
362	! 1 or -1, first line is number of values, number of positions
363	! For backward simulations, the unit is seconds, stored in grid_time
364	!*******************************************************************
365
366	! Concentration output
367	!*********************
368	if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
369
370	! Wet deposition
371	sp_count_i=0
372	sp_count_r=0
373	sp_fact=-1.
374	sp_zer=.true.
375	if ((ldirect.eq.1).and.(WETDEP)) then
376	do jy=0,numygrid-1
377	do ix=0,numxgrid-1
378	!oncentraion greater zero
379	if (wetgrid(ix,jy).gt.smallnum) then
380	if (sp_zer.eqv..true.) then ! first non zero value
381	sp_count_i=sp_count_i+1
382	sparse_dump_i(sp_count_i)=ix+jy*numxgrid
383	sp_zer=.false.
384	sp_fact=sp_fact*(-1.)
385	endif
386	sp_count_r=sp_count_r+1
387	sparse_dump_r(sp_count_r)= &
388	sp_fact1.e12wetgrid(ix,jy)/area(ix,jy)
389	! sparse_dump_u(sp_count_r)=
390	!+ 1.e12*wetgridsigma(ix,jy,ks,kp,nage)/area(ix,jy)
391	else ! concentration is zero
392	sp_zer=.true.
393	endif
394	end do
395	end do
396	else
397	sp_count_i=0
398	sp_count_r=0
399	endif
400	write(unitoutgrid) sp_count_i
401	write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
402	write(unitoutgrid) sp_count_r
403	write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
404	! write(unitoutgrid) sp_count_u
405	! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
406
407	! Dry deposition
408	sp_count_i=0
409	sp_count_r=0
410	sp_fact=-1.
411	sp_zer=.true.
412	if ((ldirect.eq.1).and.(DRYDEP)) then
413	do jy=0,numygrid-1
414	do ix=0,numxgrid-1
415	if (drygrid(ix,jy).gt.smallnum) then
416	if (sp_zer.eqv..true.) then ! first non zero value
417	sp_count_i=sp_count_i+1
418	sparse_dump_i(sp_count_i)=ix+jy*numxgrid
419	sp_zer=.false.
420	sp_fact=sp_fact*(-1.)
421	endif
422	sp_count_r=sp_count_r+1
423	sparse_dump_r(sp_count_r)= &
424	sp_fact* &
425	1.e12*drygrid(ix,jy)/area(ix,jy)
426	! sparse_dump_u(sp_count_r)=
427	!+ 1.e12*drygridsigma(ix,jy,ks,kp,nage)/area(ix,jy)
428	else ! concentration is zero
429	sp_zer=.true.
430	endif
431	end do
432	end do
433	else
434	sp_count_i=0
435	sp_count_r=0
436	endif
437	write(unitoutgrid) sp_count_i
438	write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
439	write(unitoutgrid) sp_count_r
440	write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
441	! write(,) sp_count_u
442	! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
443
444
445
446	! Concentrations
447	sp_count_i=0
448	sp_count_r=0
449	sp_fact=-1.
450	sp_zer=.true.
451	do kz=1,numzgrid
452	do jy=0,numygrid-1
453	do ix=0,numxgrid-1
454	if (grid(ix,jy,kz).gt.smallnum) then
455	if (sp_zer.eqv..true.) then ! first non zero value
456	sp_count_i=sp_count_i+1
457	sparse_dump_i(sp_count_i)= &
458	ix+jynumxgrid+kznumxgrid*numygrid
459	sp_zer=.false.
460	sp_fact=sp_fact*(-1.)
461	endif
462	sp_count_r=sp_count_r+1
463	if (lparticlecountoutput) then
464	sparse_dump_r(sp_count_r)= &
465	sp_fact* &
466	grid(ix,jy,kz)
467	else
468	sparse_dump_r(sp_count_r)= &
469	sp_fact* &
470	grid(ix,jy,kz)* &
471	factor3d(ix,jy,kz)/tot_mu(ks,kp)
472	end if
473
474
475	! if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)
476	! + write (,) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp
477	! sparse_dump_u(sp_count_r)=
478	!+ ,gridsigma(ix,jy,kz,ks,kp,nage)*
479	!+ factor(ix,jy,kz)/tot_mu(ks,kp)
480	else ! concentration is zero
481	sp_zer=.true.
482	endif
483	end do
484	end do
485	end do
486	write(unitoutgrid) sp_count_i
487	write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
488	write(unitoutgrid) sp_count_r
489	write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
490	! write(unitoutgrid) sp_count_u
491	! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
492
493
494
495	endif ! concentration output
496
497	! Mixing ratio output
498	!********************
499
500	if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio
501
502	! Wet deposition
503	sp_count_i=0
504	sp_count_r=0
505	sp_fact=-1.
506	sp_zer=.true.
507	if ((ldirect.eq.1).and.(WETDEP)) then
508	do jy=0,numygrid-1
509	do ix=0,numxgrid-1
510	if (wetgrid(ix,jy).gt.smallnum) then
511	if (sp_zer.eqv..true.) then ! first non zero value
512	sp_count_i=sp_count_i+1
513	sparse_dump_i(sp_count_i)= &
514	ix+jy*numxgrid
515	sp_zer=.false.
516	sp_fact=sp_fact*(-1.)
517	endif
518	sp_count_r=sp_count_r+1
519	sparse_dump_r(sp_count_r)= &
520	sp_fact* &
521	1.e12*wetgrid(ix,jy)/area(ix,jy)
522	! sparse_dump_u(sp_count_r)=
523	! + ,1.e12*wetgridsigma(ix,jy,ks,kp,nage)/area(ix,jy)
524	else ! concentration is zero
525	sp_zer=.true.
526	endif
527	end do
528	end do
529	else
530	sp_count_i=0
531	sp_count_r=0
532	endif
533	write(unitoutgridppt) sp_count_i
534	write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
535	write(unitoutgridppt) sp_count_r
536	write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
537	! write(unitoutgridppt) sp_count_u
538	! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
539
540
541	! Dry deposition
542	sp_count_i=0
543	sp_count_r=0
544	sp_fact=-1.
545	sp_zer=.true.
546	if ((ldirect.eq.1).and.(DRYDEP)) then
547	do jy=0,numygrid-1
548	do ix=0,numxgrid-1
549	if (drygrid(ix,jy).gt.smallnum) then
550	if (sp_zer.eqv..true.) then ! first non zero value
551	sp_count_i=sp_count_i+1
552	sparse_dump_i(sp_count_i)= &
553	ix+jy*numxgrid
554	sp_zer=.false.
555	sp_fact=sp_fact*(-1)
556	endif
557	sp_count_r=sp_count_r+1
558	sparse_dump_r(sp_count_r)= &
559	sp_fact* &
560	1.e12*drygrid(ix,jy)/area(ix,jy)
561	! sparse_dump_u(sp_count_r)=
562	! + ,1.e12*drygridsigma(ix,jy,ks,kp,nage)/area(ix,jy)
563	else ! concentration is zero
564	sp_zer=.true.
565	endif
566	end do
567	end do
568	else
569	sp_count_i=0
570	sp_count_r=0
571	endif
572	write(unitoutgridppt) sp_count_i
573	write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
574	write(unitoutgridppt) sp_count_r
575	write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
576	! write(unitoutgridppt) sp_count_u
577	! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
578
579
580	! Mixing ratios
581	sp_count_i=0
582	sp_count_r=0
583	sp_fact=-1.
584	sp_zer=.true.
585	do kz=1,numzgrid
586	do jy=0,numygrid-1
587	do ix=0,numxgrid-1
588	if (grid(ix,jy,kz).gt.smallnum) then
589	if (sp_zer.eqv..true.) then ! first non zero value
590	sp_count_i=sp_count_i+1
591	sparse_dump_i(sp_count_i)= &
592	ix+jynumxgrid+kznumxgrid*numygrid
593	sp_zer=.false.
594	sp_fact=sp_fact*(-1.)
595	endif
596	sp_count_r=sp_count_r+1
597	sparse_dump_r(sp_count_r)= &
598	sp_fact* &
599	1.e12*grid(ix,jy,kz) &
600	/volume(ix,jy,kz)/outnum* &
601	weightair/weightmolar(ks)/densityoutgrid(ix,jy,kz)
602	! sparse_dump_u(sp_count_r)=
603	!+ ,1.e12*gridsigma(ix,jy,kz,ks,kp,nage)/volume(ix,jy,kz)/
604	!+ outnum*weightair/weightmolar(ks)/
605	!+ densityoutgrid(ix,jy,kz)
606	else ! concentration is zero
607	sp_zer=.true.
608	endif
609	end do
610	end do
611	end do
612	write(unitoutgridppt) sp_count_i
613	write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
614	write(unitoutgridppt) sp_count_r
615	write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
616	! write(unitoutgridppt) sp_count_u
617	! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
618
619	endif ! output for ppt
620
621	end do
622	end do
623
624	close(unitoutgridppt)
625	close(unitoutgrid)
626
627	end do
628
629	! RLT Aug 2017
630	! Write out conversion factor for dry air
631	inquire(file=path(2)(1:length(2))//'factor_drygrid',exist=lexist)
632	if (lexist) then
633	! open and append
634	open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
635	status='old',action='write',access='append')
636	else
637	! create new
638	open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
639	status='new',action='write')
640	endif
641	sp_count_i=0
642	sp_count_r=0
643	sp_fact=-1.
644	sp_zer=.true.
645	do kz=1,numzgrid
646	do jy=0,numygrid-1
647	do ix=0,numxgrid-1
648	if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
649	if (sp_zer.eqv..true.) then ! first value not equal to one
650	sp_count_i=sp_count_i+1
651	sparse_dump_i(sp_count_i)= &
652	ix+jynumxgrid+kznumxgrid*numygrid
653	sp_zer=.false.
654	sp_fact=sp_fact*(-1.)
655	endif
656	sp_count_r=sp_count_r+1
657	sparse_dump_r(sp_count_r)= &
658	sp_fact*factor_drygrid(ix,jy,kz)
659	else ! factor is one
660	sp_zer=.true.
661	endif
662	end do
663	end do
664	end do
665	write(unitoutfactor) sp_count_i
666	write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
667	write(unitoutfactor) sp_count_r
668	write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
669	close(unitoutfactor)
670
671
672	if (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal
673	if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ &
674	wetgridtotal
675	if (drygridtotal.gt.0.) drygridtotalunc=drygridsigmatotal/ &
676	drygridtotal
677
678	! Dump of receptor concentrations
679
680	if (numreceptor.gt.0 .and. (iout.eq.2 .or. iout.eq.3) ) then
681	write(unitoutreceptppt) itime
682	do ks=1,nspec
683	write(unitoutreceptppt) (1.e12creceptor(i,ks)/outnum &
684	weightair/weightmolar(ks)/densityoutrecept(i),i=1,numreceptor)
685	end do
686	endif
687
688	! Dump of receptor concentrations
689
690	if (numreceptor.gt.0) then
691	write(unitoutrecept) itime
692	do ks=1,nspec
693	write(unitoutrecept) (1.e12*creceptor(i,ks)/outnum, &
694	i=1,numreceptor)
695	end do
696	endif
697
698	! RLT Aug 2017
699	! Write out conversion factor for dry air
700	if (numreceptor.gt.0) then
701	inquire(file=path(2)(1:length(2))//'factor_dryreceptor',exist=lexist)
702	if (lexist) then
703	! open and append
704	open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
705	status='old',action='write',access='append')
706	else
707	! create new
708	open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
709	status='new',action='write')
710	endif
711	write(unitoutfactor) itime
712	write(unitoutfactor) (factor_dryrecept(i),i=1,numreceptor)
713	close(unitoutfactor)
714	endif
715
716	! Reinitialization of grid
717	!*************************
718
719	! do ks=1,nspec
720	! do kp=1,maxpointspec_act
721	! do i=1,numreceptor
722	! creceptor(i,ks)=0.
723	! end do
724	! do jy=0,numygrid-1
725	! do ix=0,numxgrid-1
726	! do l=1,nclassunc
727	! do nage=1,nageclass
728	! do kz=1,numzgrid
729	! gridunc(ix,jy,kz,ks,kp,l,nage)=0.
730	! end do
731	! end do
732	! end do
733	! end do
734	! end do
735	! end do
736	! end do
737	creceptor(:,:)=0.
738	gridunc(:,:,:,:,:,:,:)=0.
739
740
741	end subroutine concoutput

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

Download in other formats:

Original Format