Context Navigation

concoutput_mpi.f90 @ f3054ea

GFS_025dev

Last change on this file since f3054ea was f3054ea, checked in by Espen Sollum <eso@…>, 4 years ago
Changed from grib_api to eccodes. MPI: implemented linversionout=1; fixed calculation of grid_initial fields.
Property mode set to `100644`
File size: 26.7 KB

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

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

Download in other formats:

Original Format