Context Navigation

concoutput_inversion.f90 @ e048999

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

Last change on this file since e048999 was e048999, checked in by Ignacio Pisso <ip@…>, 5 years ago
replace the nonstandard call flush () with the F2003 standard flush () in concoutput_inversion.f90
Property mode set to `100644`
File size: 31.1 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_inversion(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	! January 2017, Separate files by release but include all timesteps
51	! *
52	!*****************************************************************************
53	! *
54	! Variables: *
55	! outnum number of samples *
56	! ncells number of cells with non-zero concentrations *
57	! sparse .true. if in sparse matrix format, else .false. *
58	! tot_mu 1 for forward, initial mass mixing ration for backw. runs *
59	! *
60	!*****************************************************************************
61
62	use unc_mod
63	use point_mod
64	use outg_mod
65	use par_mod
66	use com_mod
67	use mean_mod
68
69	implicit none
70
71	real(kind=dp) :: jul
72	integer :: itime,i,ix,jy,kz,ks,kp,l,iix,jjy,kzz,nage,jjjjmmdd,ihmmss
73	integer :: sp_count_i,sp_count_r
74	real :: sp_fact
75	real :: outnum,densityoutrecept(maxreceptor),xl,yl
76	! RLT
77	real :: densitydryrecept(maxreceptor)
78	real :: factor_dryrecept(maxreceptor)
79
80	!real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
81	! +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,
82	! + maxageclass)
83	!real wetgrid(0:numxgrid-1,0:numygrid-1,maxspec,maxpointspec_act,
84	! + maxageclass)
85	!real drygrid(0:numxgrid-1,0:numygrid-1,maxspec,
86	! + maxpointspec_act,maxageclass)
87	!real gridsigma(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,
88	! + maxpointspec_act,maxageclass),
89	! + drygridsigma(0:numxgrid-1,0:numygrid-1,maxspec,
90	! + maxpointspec_act,maxageclass),
91	! + wetgridsigma(0:numxgrid-1,0:numygrid-1,maxspec,
92	! + maxpointspec_act,maxageclass)
93	!real factor(0:numxgrid-1,0:numygrid-1,numzgrid)
94	!real sparse_dump_r(numxgridnumygridnumzgrid)
95	!integer sparse_dump_i(numxgridnumygridnumzgrid)
96	!real sparse_dump_u(numxgridnumygridnumzgrid)
97
98	real(dep_prec) :: auxgrid(nclassunc)
99	real(sp) :: gridtotal,gridsigmatotal,gridtotalunc
100	real(dep_prec) :: wetgridtotal,wetgridsigmatotal,wetgridtotalunc
101	real(dep_prec) :: drygridtotal,drygridsigmatotal,drygridtotalunc
102	real :: halfheight,dz,dz1,dz2,tot_mu(maxspec,maxpointspec_act)
103	real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
104	real,parameter :: weightair=28.97
105	logical :: sp_zer
106	character :: adate8,atime6
107	character(len=3) :: anspec
108	logical :: lexist
109	character :: areldate8,areltime6
110	logical,save :: lstart=.true.
111	logical,save,allocatable,dimension(:) :: lstartrel
112	integer :: ierr
113	character(LEN=100) :: dates_char
114	integer, parameter :: unitrelnames=654
115
116
117	if(lstart) then
118	allocate(lstartrel(maxpointspec_act))
119	lstartrel(:)=.true.
120	endif
121	print*, 'lstartrel = ',lstartrel
122
123	if (verbosity.eq.1) then
124	print*,'inside concoutput_inversion '
125	CALL SYSTEM_CLOCK(count_clock)
126	WRITE(,) 'SYSTEM_CLOCK',count_clock - count_clock0
127	endif
128
129	! Determine current calendar date
130	!**********************************************************
131
132	jul=bdate+real(itime,kind=dp)/86400._dp
133	call caldate(jul,jjjjmmdd,ihmmss)
134	write(adate,'(i8.8)') jjjjmmdd
135	write(atime,'(i6.6)') ihmmss
136	! write(unitdates,'(a)') adate//atime
137
138
139	! Prepare output files for dates
140	!**********************************************************
141
142	! Overwrite existing dates file on first call, later append to it
143	! If 'dates' file exists in output directory, copy to new file dates.old
144	inquire(file=path(2)(1:length(2))//'dates', exist=lexist)
145	if (lexist.and.lstart) then
146	! copy contents of existing dates file to dates.old
147	print*, 'warning: replacing old dates file'
148	open(unit=unitdates, file=path(2)(1:length(2))//'dates',form='formatted', &
149	&access='sequential', status='old', action='read', iostat=ierr)
150	open(unit=unittmp, file=path(2)(1:length(2))//'dates.old', access='sequential', &
151	&status='replace', action='write', form='formatted', iostat=ierr)
152	do while (.true.)
153	read(unitdates, '(a)', iostat=ierr) dates_char
154	if (ierr.ne.0) exit
155	write(unit=unittmp, fmt='(a)', iostat=ierr, advance='yes') trim(dates_char)
156	end do
157	close(unit=unitdates)
158	close(unit=unittmp)
159	! create new dates file
160	open(unit=unitdates, file=path(2)(1:length(2))//'dates',form='formatted', &
161	&access='sequential', status='replace', iostat=ierr)
162	close(unit=unitdates)
163	endif
164
165	open(unitdates,file=path(2)(1:length(2))//'dates', ACCESS='APPEND')
166	write(unitdates,'(a)') adate//atime
167	close(unitdates)
168
169	!CGZ: Make a filename with names of releases
170	if (lstart) then
171	open(unit=unitrelnames, file=path(2)(1:length(2))//'releases_out',form='formatted', &
172	&access='sequential', status='replace', iostat=ierr)
173	close(unitrelnames)
174	endif
175
176	print*, 'after creating dates files: lstart = ',lstart
177	! print*, 'outnum:',outnum
178	! print*, 'datetime:',adate//atime
179
180
181	! For forward simulations, output fields have dimension MAXSPEC,
182	! for backward simulations, output fields have dimension MAXPOINT.
183	! Thus, make loops either about nspec, or about numpoint
184	!*****************************************************************
185
186
187	if (ldirect.eq.1) then
188	do ks=1,nspec
189	do kp=1,maxpointspec_act
190	tot_mu(ks,kp)=1
191	end do
192	end do
193	else
194	do ks=1,nspec
195	do kp=1,maxpointspec_act
196	tot_mu(ks,kp)=xmass(kp,ks)
197	end do
198	end do
199	endif
200
201
202	if (verbosity.eq.1) then
203	print*,'concoutput_inversion 2'
204	CALL SYSTEM_CLOCK(count_clock)
205	WRITE(,) 'SYSTEM_CLOCK',count_clock - count_clock0
206	endif
207
208	!*******************************************************************
209	! Compute air density: sufficiently accurate to take it
210	! from coarse grid at some time
211	! Determine center altitude of output layer, and interpolate density
212	! data to that altitude
213	!*******************************************************************
214
215	do kz=1,numzgrid
216	if (kz.eq.1) then
217	halfheight=outheight(1)/2.
218	else
219	halfheight=(outheight(kz)+outheight(kz-1))/2.
220	endif
221	do kzz=2,nz
222	if ((height(kzz-1).lt.halfheight).and. &
223	(height(kzz).gt.halfheight)) goto 46
224	end do
225	46 kzz=max(min(kzz,nz),2)
226	dz1=halfheight-height(kzz-1)
227	dz2=height(kzz)-halfheight
228	dz=dz1+dz2
229	do jy=0,numygrid-1
230	do ix=0,numxgrid-1
231	xl=outlon0+real(ix)*dxout
232	yl=outlat0+real(jy)*dyout
233	xl=(xl-xlon0)/dx
234	yl=(yl-ylat0)/dy
235	iix=max(min(nint(xl),nxmin1),0)
236	jjy=max(min(nint(yl),nymin1),0)
237	densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
238	rho(iix,jjy,kzz-1,2)*dz2)/dz
239	! RLT
240	densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,2)*dz1+ &
241	rho_dry(iix,jjy,kzz-1,2)*dz2)/dz
242	end do
243	end do
244	end do
245
246	do i=1,numreceptor
247	xl=xreceptor(i)
248	yl=yreceptor(i)
249	iix=max(min(nint(xl),nxmin1),0)
250	jjy=max(min(nint(yl),nymin1),0)
251	densityoutrecept(i)=rho(iix,jjy,1,2)
252	! RLT
253	densitydryrecept(i)=rho_dry(iix,jjy,1,2)
254	end do
255
256	! RLT
257	! conversion factor for output relative to dry air
258	factor_drygrid=densityoutgrid/densitydrygrid
259	factor_dryrecept=densityoutrecept/densitydryrecept
260
261	! Output is different for forward and backward simulations
262	do kz=1,numzgrid
263	do jy=0,numygrid-1
264	do ix=0,numxgrid-1
265	if (ldirect.eq.1) then
266	factor3d(ix,jy,kz)=1.e12/volume(ix,jy,kz)/outnum
267	else
268	factor3d(ix,jy,kz)=real(abs(loutaver))/outnum
269	endif
270	end do
271	end do
272	end do
273
274	!*********************************************************************
275	! Determine the standard deviation of the mean concentration or mixing
276	! ratio (uncertainty of the output) and the dry and wet deposition
277	!*********************************************************************
278
279	if (verbosity.eq.1) then
280	print*,'concoutput_inversion 3 (sd)'
281	CALL SYSTEM_CLOCK(count_clock)
282	WRITE(,) 'SYSTEM_CLOCK',count_clock - count_clock0
283	endif
284	gridtotal=0.
285	gridsigmatotal=0.
286	gridtotalunc=0.
287	wetgridtotal=0.
288	wetgridsigmatotal=0.
289	wetgridtotalunc=0.
290	drygridtotal=0.
291	drygridsigmatotal=0.
292	drygridtotalunc=0.
293
294	do ks=1,nspec
295
296	write(anspec,'(i3.3)') ks
297
298	! loop over releases
299	do kp=1,maxpointspec_act
300
301	print*, 'itime = ',itime
302	!print*, 'lage(1) = ',lage(1)
303	print*, 'ireleasestart(kp) = ',ireleasestart(kp)
304	print*, 'ireleaseend(kp) = ',ireleaseend(kp)
305
306	! check itime is within release and backward trajectory length
307	if (nageclass.eq.1) then
308	if ((itime.gt.ireleaseend(kp)).or.(itime.lt.(ireleasestart(kp)-lage(1)))) then
309	go to 10
310	endif
311	endif
312
313	! calculate date of release for filename
314	jul=bdate+real(ireleasestart(kp),kind=dp)/86400._dp ! this is the current day
315	call caldate(jul,jjjjmmdd,ihmmss)
316	write(areldate,'(i8.8)') jjjjmmdd
317	write(areltime,'(i6.6)') ihmmss
318	print*, 'areldate/areltime = ',areldate//areltime
319
320	! calculate date of field
321	jul=bdate+real(itime,kind=dp)/86400._dp
322	call caldate(jul,jjjjmmdd,ihmmss)
323	write(adate,'(i8.8)') jjjjmmdd
324	write(atime,'(i6.6)') ihmmss
325	! print*, adate//atime
326
327	if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
328	if (ldirect.eq.1) then
329	! concentrations
330	inquire(file=path(2)(1:length(2))//'grid_conc_'//areldate// &
331	areltime//'_'//anspec,exist=lexist)
332	if(lexist.and..not.lstartrel(kp)) then
333	! open and append to existing file
334	open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//areldate// &
335	areltime//'_'//anspec,form='unformatted',status='old',action='write',access='append')
336	else
337	! open new file
338	open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//areldate// &
339	areltime//'_'//anspec,form='unformatted',status='replace',action='write')
340	endif
341	else
342	! residence times
343	inquire(file=path(2)(1:length(2))//'grid_time_'//areldate// &
344	areltime//'_'//anspec,exist=lexist)
345	if(lexist.and..not.lstartrel(kp)) then
346	! open and append to existing file
347	open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//areldate// &
348	areltime//'_'//anspec,form='unformatted',status='old',action='write',access='append')
349	else
350	! open new file
351	open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//areldate// &
352	areltime//'_'//anspec,form='unformatted',status='replace',action='write')
353	! add part of the filename to a file (similar to dates) for easier post-processing
354	open(unit=unitrelnames, file=path(2)(1:length(2))//'releases_out',form='formatted', &
355	& access='APPEND', iostat=ierr)
356	write(unitrelnames,'(a)') areldate//areltime//'_'//anspec
357	!call flush(unit=unitrelnames)
358	flush(unit=unitrelnames)
359	close(unitrelnames)
360	endif
361	endif
362	write(unitoutgrid) jjjjmmdd
363	write(unitoutgrid) ihmmss
364	endif
365
366	if ((iout.eq.2).or.(iout.eq.3)) then
367	! mixing ratio
368	inquire(file=path(2)(1:length(2))//'grid_pptv_'//areldate// &
369	areltime//'_'//anspec,exist=lexist)
370	if(lexist.and..not.lstartrel(kp)) then
371	! open and append to existing file
372	open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//areldate// &
373	areltime//'_'//anspec,form='unformatted',status='old',action='write',access='append')
374	else
375	! open new file
376	open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//areldate// &
377	areltime//'_'//anspec,form='unformatted',status='replace',action='write')
378	endif
379	write(unitoutgridppt) jjjjmmdd
380	write(unitoutgridppt) ihmmss
381	endif
382
383	lstartrel(kp)=.false.
384
385	do nage=1,nageclass
386
387	do jy=0,numygrid-1
388	do ix=0,numxgrid-1
389
390	!! WET DEPOSITION
391	! if ((WETDEP).and.(ldirect.gt.0)) then
392	! do l=1,nclassunc
393	! auxgrid(l)=wetgridunc(ix,jy,ks,kp,l,nage)
394	! end do
395	! call mean(auxgrid,wetgrid(ix,jy), &
396	! wetgridsigma(ix,jy),nclassunc)
397	!! Multiply by number of classes to get total concentration
398	! wetgrid(ix,jy)=wetgrid(ix,jy) &
399	! *nclassunc
400	! wetgridtotal=wetgridtotal+wetgrid(ix,jy)
401	!! Calculate standard deviation of the mean
402	! wetgridsigma(ix,jy)= &
403	! wetgridsigma(ix,jy)* &
404	! sqrt(real(nclassunc))
405	! wetgridsigmatotal=wetgridsigmatotal+ &
406	! wetgridsigma(ix,jy)
407	! endif
408
409	!! DRY DEPOSITION
410	! if ((DRYDEP).and.(ldirect.gt.0)) then
411	! do l=1,nclassunc
412	! auxgrid(l)=drygridunc(ix,jy,ks,kp,l,nage)
413	! end do
414	! call mean(auxgrid,drygrid(ix,jy), &
415	! drygridsigma(ix,jy),nclassunc)
416	!! Multiply by number of classes to get total concentration
417	! drygrid(ix,jy)=drygrid(ix,jy)* &
418	! nclassunc
419	! drygridtotal=drygridtotal+drygrid(ix,jy)
420	!! Calculate standard deviation of the mean
421	! drygridsigma(ix,jy)= &
422	! drygridsigma(ix,jy)* &
423	! sqrt(real(nclassunc))
424	!125 drygridsigmatotal=drygridsigmatotal+ &
425	! drygridsigma(ix,jy)
426	! endif
427
428	! CONCENTRATION OR MIXING RATIO
429	do kz=1,numzgrid
430	do l=1,nclassunc
431	auxgrid(l)=gridunc(ix,jy,kz,ks,kp,l,nage)
432	end do
433	call mean(auxgrid,grid(ix,jy,kz), &
434	gridsigma(ix,jy,kz),nclassunc)
435	! Multiply by number of classes to get total concentration
436	grid(ix,jy,kz)= &
437	grid(ix,jy,kz)*nclassunc
438	gridtotal=gridtotal+grid(ix,jy,kz)
439	! Calculate standard deviation of the mean
440	gridsigma(ix,jy,kz)= &
441	gridsigma(ix,jy,kz)* &
442	sqrt(real(nclassunc))
443	gridsigmatotal=gridsigmatotal+ &
444	gridsigma(ix,jy,kz)
445	end do
446	end do
447	end do
448
449
450	!*******************************************************************
451	! Generate output: may be in concentration (ng/m3) or in mixing
452	! ratio (ppt) or both
453	! Output the position and the values alternated multiplied by
454	! 1 or -1, first line is number of values, number of positions
455	! For backward simulations, the unit is seconds, stored in grid_time
456	!*******************************************************************
457
458	if (verbosity.eq.1) then
459	print*,'concoutput_inversion 4 (output)'
460	CALL SYSTEM_CLOCK(count_clock)
461	WRITE(,) 'SYSTEM_CLOCK',count_clock - count_clock0
462	endif
463
464	! Concentration output
465	!*********************
466
467	if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
468
469	if (verbosity.eq.1) then
470	print*,'concoutput_inversion (Wet deposition)'
471	CALL SYSTEM_CLOCK(count_clock)
472	WRITE(,) 'SYSTEM_CLOCK',count_clock - count_clock0
473	endif
474
475	! Wet deposition
476	! sp_count_i=0
477	! sp_count_r=0
478	! sp_fact=-1.
479	! sp_zer=.true.
480	! if ((ldirect.eq.1).and.(WETDEP)) then
481	! do jy=0,numygrid-1
482	! do ix=0,numxgrid-1
483	!! concentration greater zero
484	! if (wetgrid(ix,jy).gt.smallnum) then
485	! if (sp_zer.eqv..true.) then ! first non zero value
486	! sp_count_i=sp_count_i+1
487	! sparse_dump_i(sp_count_i)=ix+jy*numxgrid
488	! sp_zer=.false.
489	! sp_fact=sp_fact*(-1.)
490	! endif
491	! sp_count_r=sp_count_r+1
492	! sparse_dump_r(sp_count_r)= &
493	! sp_fact1.e12wetgrid(ix,jy)/area(ix,jy)
494	! sparse_dump_u(sp_count_r)= &
495	! 1.e12*wetgridsigma(ix,jy)/area(ix,jy)
496	! else ! concentration is zero
497	! sp_zer=.true.
498	! endif
499	! end do
500	! end do
501	! else
502	! sp_count_i=0
503	! sp_count_r=0
504	! endif
505	! write(unitoutgrid) sp_count_i
506	! write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
507	! write(unitoutgrid) sp_count_r
508	! write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
509	!! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
510	!
511	! if (verbosity.eq.1) then
512	! print*,'concoutput_surf (Dry deposition)'
513	! CALL SYSTEM_CLOCK(count_clock)
514	! WRITE(,) 'SYSTEM_CLOCK',count_clock - count_clock0
515	! endif
516	!! Dry deposition
517	! sp_count_i=0
518	! sp_count_r=0
519	! sp_fact=-1.
520	! sp_zer=.true.
521	! if ((ldirect.eq.1).and.(DRYDEP)) then
522	! do jy=0,numygrid-1
523	! do ix=0,numxgrid-1
524	! if (drygrid(ix,jy).gt.smallnum) then
525	! if (sp_zer.eqv..true.) then ! first non zero value
526	! sp_count_i=sp_count_i+1
527	! sparse_dump_i(sp_count_i)=ix+jy*numxgrid
528	! sp_zer=.false.
529	! sp_fact=sp_fact*(-1.)
530	! endif
531	! sp_count_r=sp_count_r+1
532	! sparse_dump_r(sp_count_r)= &
533	! sp_fact* &
534	! 1.e12*drygrid(ix,jy)/area(ix,jy)
535	! sparse_dump_u(sp_count_r)= &
536	! 1.e12*drygridsigma(ix,jy)/area(ix,jy)
537	! else ! concentration is zero
538	! sp_zer=.true.
539	! endif
540	! end do
541	! end do
542	! else
543	! sp_count_i=0
544	! sp_count_r=0
545	! endif
546	! write(unitoutgrid) sp_count_i
547	! write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
548	! write(unitoutgrid) sp_count_r
549	! write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
550	!! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
551
552	if (verbosity.eq.1) then
553	print*,'concoutput_inversion (Concentrations)'
554	CALL SYSTEM_CLOCK(count_clock)
555	WRITE(,) 'SYSTEM_CLOCK',count_clock - count_clock0
556	endif
557
558	! Concentrations
559
560	! surf_only write only 1st layer
561
562	sp_count_i=0
563	sp_count_r=0
564	sp_fact=-1.
565	sp_zer=.true.
566	do kz=1,1
567	do jy=0,numygrid-1
568	do ix=0,numxgrid-1
569	if (grid(ix,jy,kz).gt.smallnum) then
570	if (sp_zer.eqv..true.) then ! first non zero value
571	sp_count_i=sp_count_i+1
572	sparse_dump_i(sp_count_i)= &
573	ix+jynumxgrid+kznumxgrid*numygrid
574	sp_zer=.false.
575	sp_fact=sp_fact*(-1.)
576	endif
577	sp_count_r=sp_count_r+1
578	sparse_dump_r(sp_count_r)= &
579	sp_fact* &
580	grid(ix,jy,kz)* &
581	factor3d(ix,jy,kz)/tot_mu(ks,kp)
582	! if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)
583	! + write (,) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp
584	sparse_dump_u(sp_count_r)= &
585	gridsigma(ix,jy,kz)* &
586	factor3d(ix,jy,kz)/tot_mu(ks,kp)
587
588	else ! concentration is zero
589	sp_zer=.true.
590	endif
591	end do
592	end do
593	end do
594	write(unitoutgrid) sp_count_i
595	write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
596	write(unitoutgrid) sp_count_r
597	write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
598	! write(unitoutgrid) sp_count_r
599	! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
600
601	endif ! concentration output
602
603	! Mixing ratio output
604	!********************
605
606	if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio
607	!
608	!! Wet deposition
609	! sp_count_i=0
610	! sp_count_r=0
611	! sp_fact=-1.
612	! sp_zer=.true.
613	! if ((ldirect.eq.1).and.(WETDEP)) then
614	! do jy=0,numygrid-1
615	! do ix=0,numxgrid-1
616	! if (wetgrid(ix,jy).gt.smallnum) then
617	! if (sp_zer.eqv..true.) then ! first non zero value
618	! sp_count_i=sp_count_i+1
619	! sparse_dump_i(sp_count_i)= &
620	! ix+jy*numxgrid
621	! sp_zer=.false.
622	! sp_fact=sp_fact*(-1.)
623	! endif
624	! sp_count_r=sp_count_r+1
625	! sparse_dump_r(sp_count_r)= &
626	! sp_fact* &
627	! 1.e12*wetgrid(ix,jy)/area(ix,jy)
628	! sparse_dump_u(sp_count_r)= &
629	! 1.e12*wetgridsigma(ix,jy)/area(ix,jy)
630	! else ! concentration is zero
631	! sp_zer=.true.
632	! endif
633	! end do
634	! end do
635	! else
636	! sp_count_i=0
637	! sp_count_r=0
638	! endif
639	! write(unitoutgridppt) sp_count_i
640	! write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
641	! write(unitoutgridppt) sp_count_r
642	! write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
643	!! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
644	!
645
646	! Dry deposition
647	! sp_count_i=0
648	! sp_count_r=0
649	! sp_fact=-1.
650	! sp_zer=.true.
651	! if ((ldirect.eq.1).and.(DRYDEP)) then
652	! do jy=0,numygrid-1
653	! do ix=0,numxgrid-1
654	! if (drygrid(ix,jy).gt.smallnum) then
655	! if (sp_zer.eqv..true.) then ! first non zero value
656	! sp_count_i=sp_count_i+1
657	! sparse_dump_i(sp_count_i)= &
658	! ix+jy*numxgrid
659	! sp_zer=.false.
660	! sp_fact=sp_fact*(-1)
661	! endif
662	! sp_count_r=sp_count_r+1
663	! sparse_dump_r(sp_count_r)= &
664	! sp_fact* &
665	! 1.e12*drygrid(ix,jy)/area(ix,jy)
666	! sparse_dump_u(sp_count_r)= &
667	! 1.e12*drygridsigma(ix,jy)/area(ix,jy)
668	! else ! concentration is zero
669	! sp_zer=.true.
670	! endif
671	! end do
672	! end do
673	! else
674	! sp_count_i=0
675	! sp_count_r=0
676	! endif
677	! write(unitoutgridppt) sp_count_i
678	! write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
679	! write(unitoutgridppt) sp_count_r
680	! write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
681	!! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
682	!
683
684	! Mixing ratios
685
686	! surf_only write only 1st layer
687
688	sp_count_i=0
689	sp_count_r=0
690	sp_fact=-1.
691	sp_zer=.true.
692	do kz=1,1
693	do jy=0,numygrid-1
694	do ix=0,numxgrid-1
695	if (grid(ix,jy,kz).gt.smallnum) then
696	if (sp_zer.eqv..true.) then ! first non zero value
697	sp_count_i=sp_count_i+1
698	sparse_dump_i(sp_count_i)= &
699	ix+jynumxgrid+kznumxgrid*numygrid
700	sp_zer=.false.
701	sp_fact=sp_fact*(-1.)
702	endif
703	sp_count_r=sp_count_r+1
704	sparse_dump_r(sp_count_r)= &
705	sp_fact* &
706	1.e12*grid(ix,jy,kz) &
707	/volume(ix,jy,kz)/outnum* &
708	weightair/weightmolar(ks)/densityoutgrid(ix,jy,kz)
709	sparse_dump_u(sp_count_r)= &
710	1.e12*gridsigma(ix,jy,kz)/volume(ix,jy,kz)/ &
711	outnum*weightair/weightmolar(ks)/ &
712	densityoutgrid(ix,jy,kz)
713	else ! concentration is zero
714	sp_zer=.true.
715	endif
716	end do
717	end do
718	end do
719	write(unitoutgridppt) sp_count_i
720	write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
721	write(unitoutgridppt) sp_count_r
722	write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
723	! write(unitoutgridppt) sp_count_r
724	! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
725
726	endif ! output for ppt
727
728	end do ! nageclass
729
730	close(unitoutgridppt)
731	close(unitoutgrid)
732
733	! itime is outside range
734	10 continue
735
736	end do ! maxpointspec_act
737
738	end do ! nspec
739
740	! RLT Aug 2017
741	! Write out conversion factor for dry air
742	inquire(file=path(2)(1:length(2))//'factor_drygrid',exist=lexist)
743	if (lexist.and..not.lstart) then
744	! open and append
745	open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
746	status='old',action='write',access='append')
747	else
748	! create new
749	open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
750	status='replace',action='write')
751	endif
752	sp_count_i=0
753	sp_count_r=0
754	sp_fact=-1.
755	sp_zer=.true.
756	do kz=1,1
757	do jy=0,numygrid-1
758	do ix=0,numxgrid-1
759	if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
760	if (sp_zer.eqv..true.) then ! first value not equal to one
761	sp_count_i=sp_count_i+1
762	sparse_dump_i(sp_count_i)= &
763	ix+jynumxgrid+kznumxgrid*numygrid
764	sp_zer=.false.
765	sp_fact=sp_fact*(-1.)
766	endif
767	sp_count_r=sp_count_r+1
768	sparse_dump_r(sp_count_r)= &
769	sp_fact*factor_drygrid(ix,jy,kz)
770	else ! factor is one
771	sp_zer=.true.
772	endif
773	end do
774	end do
775	end do
776	write(unitoutfactor) sp_count_i
777	write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
778	write(unitoutfactor) sp_count_r
779	write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
780	close(unitoutfactor)
781
782
783	if (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal
784	! if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ &
785	! wetgridtotal
786	! if (drygridtotal.gt.0.) drygridtotalunc=drygridsigmatotal/ &
787	! drygridtotal
788
789	! Dump of receptor concentrations
790
791	if (numreceptor.gt.0 .and. (iout.eq.2 .or. iout.eq.3) ) then
792	write(unitoutreceptppt) itime
793	do ks=1,nspec
794	write(unitoutreceptppt) (1.e12creceptor(i,ks)/outnum &
795	weightair/weightmolar(ks)/densityoutrecept(i),i=1,numreceptor)
796	end do
797	endif
798
799	! Dump of receptor concentrations
800
801	if (numreceptor.gt.0) then
802	write(unitoutrecept) itime
803	do ks=1,nspec
804	write(unitoutrecept) (1.e12*creceptor(i,ks)/outnum, &
805	i=1,numreceptor)
806	end do
807	endif
808
809	! RLT Aug 2017
810	! Write out conversion factor for dry air
811	if (numreceptor.gt.0) then
812	inquire(file=path(2)(1:length(2))//'factor_dryreceptor',exist=lexist)
813	if (lexist.and..not.lstart) then
814	! open and append
815	open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
816	status='old',action='write',access='append')
817	else
818	! create new
819	open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
820	status='replace',action='write')
821	endif
822	write(unitoutfactor) itime
823	write(unitoutfactor) (factor_dryrecept(i),i=1,numreceptor)
824	close(unitoutfactor)
825	endif
826
827	! reset lstart
828	if (lstart) then
829	lstart=.false.
830	endif
831	print*, 'after writing output files: lstart = ',lstart
832
833
834	! Reinitialization of grid
835	!*************************
836
837	do ks=1,nspec
838	do kp=1,maxpointspec_act
839	do i=1,numreceptor
840	creceptor(i,ks)=0.
841	end do
842	do jy=0,numygrid-1
843	do ix=0,numxgrid-1
844	do l=1,nclassunc
845	do nage=1,nageclass
846	do kz=1,numzgrid
847	gridunc(ix,jy,kz,ks,kp,l,nage)=0.
848	end do
849	end do
850	end do
851	end do
852	end do
853	end do
854	end do
855
856	end subroutine concoutput_inversion

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

Download in other formats:

Original Format