Context Navigation

concoutput_inversion.f90 @ 10bfff9

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

Last change on this file since 10bfff9 was 10bfff9, checked in by ronesy <rlt@…>, 5 years ago
removed flush statement from concoutput_inversion
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	close(unitrelnames)
358	endif
359	endif
360	write(unitoutgrid) jjjjmmdd
361	write(unitoutgrid) ihmmss
362	endif
363
364	if ((iout.eq.2).or.(iout.eq.3)) then
365	! mixing ratio
366	inquire(file=path(2)(1:length(2))//'grid_pptv_'//areldate// &
367	areltime//'_'//anspec,exist=lexist)
368	if(lexist.and..not.lstartrel(kp)) then
369	! open and append to existing file
370	open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//areldate// &
371	areltime//'_'//anspec,form='unformatted',status='old',action='write',access='append')
372	else
373	! open new file
374	open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//areldate// &
375	areltime//'_'//anspec,form='unformatted',status='replace',action='write')
376	endif
377	write(unitoutgridppt) jjjjmmdd
378	write(unitoutgridppt) ihmmss
379	endif
380
381	lstartrel(kp)=.false.
382
383	do nage=1,nageclass
384
385	do jy=0,numygrid-1
386	do ix=0,numxgrid-1
387
388	!! WET DEPOSITION
389	! if ((WETDEP).and.(ldirect.gt.0)) then
390	! do l=1,nclassunc
391	! auxgrid(l)=wetgridunc(ix,jy,ks,kp,l,nage)
392	! end do
393	! call mean(auxgrid,wetgrid(ix,jy), &
394	! wetgridsigma(ix,jy),nclassunc)
395	!! Multiply by number of classes to get total concentration
396	! wetgrid(ix,jy)=wetgrid(ix,jy) &
397	! *nclassunc
398	! wetgridtotal=wetgridtotal+wetgrid(ix,jy)
399	!! Calculate standard deviation of the mean
400	! wetgridsigma(ix,jy)= &
401	! wetgridsigma(ix,jy)* &
402	! sqrt(real(nclassunc))
403	! wetgridsigmatotal=wetgridsigmatotal+ &
404	! wetgridsigma(ix,jy)
405	! endif
406
407	!! DRY DEPOSITION
408	! if ((DRYDEP).and.(ldirect.gt.0)) then
409	! do l=1,nclassunc
410	! auxgrid(l)=drygridunc(ix,jy,ks,kp,l,nage)
411	! end do
412	! call mean(auxgrid,drygrid(ix,jy), &
413	! drygridsigma(ix,jy),nclassunc)
414	!! Multiply by number of classes to get total concentration
415	! drygrid(ix,jy)=drygrid(ix,jy)* &
416	! nclassunc
417	! drygridtotal=drygridtotal+drygrid(ix,jy)
418	!! Calculate standard deviation of the mean
419	! drygridsigma(ix,jy)= &
420	! drygridsigma(ix,jy)* &
421	! sqrt(real(nclassunc))
422	!125 drygridsigmatotal=drygridsigmatotal+ &
423	! drygridsigma(ix,jy)
424	! endif
425
426	! CONCENTRATION OR MIXING RATIO
427	do kz=1,numzgrid
428	do l=1,nclassunc
429	auxgrid(l)=gridunc(ix,jy,kz,ks,kp,l,nage)
430	end do
431	call mean(auxgrid,grid(ix,jy,kz), &
432	gridsigma(ix,jy,kz),nclassunc)
433	! Multiply by number of classes to get total concentration
434	grid(ix,jy,kz)= &
435	grid(ix,jy,kz)*nclassunc
436	gridtotal=gridtotal+grid(ix,jy,kz)
437	! Calculate standard deviation of the mean
438	gridsigma(ix,jy,kz)= &
439	gridsigma(ix,jy,kz)* &
440	sqrt(real(nclassunc))
441	gridsigmatotal=gridsigmatotal+ &
442	gridsigma(ix,jy,kz)
443	end do
444	end do
445	end do
446
447
448	!*******************************************************************
449	! Generate output: may be in concentration (ng/m3) or in mixing
450	! ratio (ppt) or both
451	! Output the position and the values alternated multiplied by
452	! 1 or -1, first line is number of values, number of positions
453	! For backward simulations, the unit is seconds, stored in grid_time
454	!*******************************************************************
455
456	if (verbosity.eq.1) then
457	print*,'concoutput_inversion 4 (output)'
458	CALL SYSTEM_CLOCK(count_clock)
459	WRITE(,) 'SYSTEM_CLOCK',count_clock - count_clock0
460	endif
461
462	! Concentration output
463	!*********************
464
465	if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
466
467	if (verbosity.eq.1) then
468	print*,'concoutput_inversion (Wet deposition)'
469	CALL SYSTEM_CLOCK(count_clock)
470	WRITE(,) 'SYSTEM_CLOCK',count_clock - count_clock0
471	endif
472
473	! Wet deposition
474	! sp_count_i=0
475	! sp_count_r=0
476	! sp_fact=-1.
477	! sp_zer=.true.
478	! if ((ldirect.eq.1).and.(WETDEP)) then
479	! do jy=0,numygrid-1
480	! do ix=0,numxgrid-1
481	!! concentration greater zero
482	! if (wetgrid(ix,jy).gt.smallnum) then
483	! if (sp_zer.eqv..true.) then ! first non zero value
484	! sp_count_i=sp_count_i+1
485	! sparse_dump_i(sp_count_i)=ix+jy*numxgrid
486	! sp_zer=.false.
487	! sp_fact=sp_fact*(-1.)
488	! endif
489	! sp_count_r=sp_count_r+1
490	! sparse_dump_r(sp_count_r)= &
491	! sp_fact1.e12wetgrid(ix,jy)/area(ix,jy)
492	! sparse_dump_u(sp_count_r)= &
493	! 1.e12*wetgridsigma(ix,jy)/area(ix,jy)
494	! else ! concentration is zero
495	! sp_zer=.true.
496	! endif
497	! end do
498	! end do
499	! else
500	! sp_count_i=0
501	! sp_count_r=0
502	! endif
503	! write(unitoutgrid) sp_count_i
504	! write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
505	! write(unitoutgrid) sp_count_r
506	! write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
507	!! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
508	!
509	! if (verbosity.eq.1) then
510	! print*,'concoutput_surf (Dry deposition)'
511	! CALL SYSTEM_CLOCK(count_clock)
512	! WRITE(,) 'SYSTEM_CLOCK',count_clock - count_clock0
513	! endif
514	!! Dry deposition
515	! sp_count_i=0
516	! sp_count_r=0
517	! sp_fact=-1.
518	! sp_zer=.true.
519	! if ((ldirect.eq.1).and.(DRYDEP)) then
520	! do jy=0,numygrid-1
521	! do ix=0,numxgrid-1
522	! if (drygrid(ix,jy).gt.smallnum) then
523	! if (sp_zer.eqv..true.) then ! first non zero value
524	! sp_count_i=sp_count_i+1
525	! sparse_dump_i(sp_count_i)=ix+jy*numxgrid
526	! sp_zer=.false.
527	! sp_fact=sp_fact*(-1.)
528	! endif
529	! sp_count_r=sp_count_r+1
530	! sparse_dump_r(sp_count_r)= &
531	! sp_fact* &
532	! 1.e12*drygrid(ix,jy)/area(ix,jy)
533	! sparse_dump_u(sp_count_r)= &
534	! 1.e12*drygridsigma(ix,jy)/area(ix,jy)
535	! else ! concentration is zero
536	! sp_zer=.true.
537	! endif
538	! end do
539	! end do
540	! else
541	! sp_count_i=0
542	! sp_count_r=0
543	! endif
544	! write(unitoutgrid) sp_count_i
545	! write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
546	! write(unitoutgrid) sp_count_r
547	! write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
548	!! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
549
550	if (verbosity.eq.1) then
551	print*,'concoutput_inversion (Concentrations)'
552	CALL SYSTEM_CLOCK(count_clock)
553	WRITE(,) 'SYSTEM_CLOCK',count_clock - count_clock0
554	endif
555
556	! Concentrations
557
558	! surf_only write only 1st layer
559
560	sp_count_i=0
561	sp_count_r=0
562	sp_fact=-1.
563	sp_zer=.true.
564	do kz=1,1
565	do jy=0,numygrid-1
566	do ix=0,numxgrid-1
567	if (grid(ix,jy,kz).gt.smallnum) then
568	if (sp_zer.eqv..true.) then ! first non zero value
569	sp_count_i=sp_count_i+1
570	sparse_dump_i(sp_count_i)= &
571	ix+jynumxgrid+kznumxgrid*numygrid
572	sp_zer=.false.
573	sp_fact=sp_fact*(-1.)
574	endif
575	sp_count_r=sp_count_r+1
576	sparse_dump_r(sp_count_r)= &
577	sp_fact* &
578	grid(ix,jy,kz)* &
579	factor3d(ix,jy,kz)/tot_mu(ks,kp)
580	! if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)
581	! + write (,) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp
582	sparse_dump_u(sp_count_r)= &
583	gridsigma(ix,jy,kz)* &
584	factor3d(ix,jy,kz)/tot_mu(ks,kp)
585
586	else ! concentration is zero
587	sp_zer=.true.
588	endif
589	end do
590	end do
591	end do
592	write(unitoutgrid) sp_count_i
593	write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
594	write(unitoutgrid) sp_count_r
595	write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
596	! write(unitoutgrid) sp_count_r
597	! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
598
599	endif ! concentration output
600
601	! Mixing ratio output
602	!********************
603
604	if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio
605	!
606	!! Wet deposition
607	! sp_count_i=0
608	! sp_count_r=0
609	! sp_fact=-1.
610	! sp_zer=.true.
611	! if ((ldirect.eq.1).and.(WETDEP)) then
612	! do jy=0,numygrid-1
613	! do ix=0,numxgrid-1
614	! if (wetgrid(ix,jy).gt.smallnum) then
615	! if (sp_zer.eqv..true.) then ! first non zero value
616	! sp_count_i=sp_count_i+1
617	! sparse_dump_i(sp_count_i)= &
618	! ix+jy*numxgrid
619	! sp_zer=.false.
620	! sp_fact=sp_fact*(-1.)
621	! endif
622	! sp_count_r=sp_count_r+1
623	! sparse_dump_r(sp_count_r)= &
624	! sp_fact* &
625	! 1.e12*wetgrid(ix,jy)/area(ix,jy)
626	! sparse_dump_u(sp_count_r)= &
627	! 1.e12*wetgridsigma(ix,jy)/area(ix,jy)
628	! else ! concentration is zero
629	! sp_zer=.true.
630	! endif
631	! end do
632	! end do
633	! else
634	! sp_count_i=0
635	! sp_count_r=0
636	! endif
637	! write(unitoutgridppt) sp_count_i
638	! write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
639	! write(unitoutgridppt) sp_count_r
640	! write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
641	!! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
642	!
643
644	! Dry deposition
645	! sp_count_i=0
646	! sp_count_r=0
647	! sp_fact=-1.
648	! sp_zer=.true.
649	! if ((ldirect.eq.1).and.(DRYDEP)) then
650	! do jy=0,numygrid-1
651	! do ix=0,numxgrid-1
652	! if (drygrid(ix,jy).gt.smallnum) then
653	! if (sp_zer.eqv..true.) then ! first non zero value
654	! sp_count_i=sp_count_i+1
655	! sparse_dump_i(sp_count_i)= &
656	! ix+jy*numxgrid
657	! sp_zer=.false.
658	! sp_fact=sp_fact*(-1)
659	! endif
660	! sp_count_r=sp_count_r+1
661	! sparse_dump_r(sp_count_r)= &
662	! sp_fact* &
663	! 1.e12*drygrid(ix,jy)/area(ix,jy)
664	! sparse_dump_u(sp_count_r)= &
665	! 1.e12*drygridsigma(ix,jy)/area(ix,jy)
666	! else ! concentration is zero
667	! sp_zer=.true.
668	! endif
669	! end do
670	! end do
671	! else
672	! sp_count_i=0
673	! sp_count_r=0
674	! endif
675	! write(unitoutgridppt) sp_count_i
676	! write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
677	! write(unitoutgridppt) sp_count_r
678	! write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
679	!! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
680	!
681
682	! Mixing ratios
683
684	! surf_only write only 1st layer
685
686	sp_count_i=0
687	sp_count_r=0
688	sp_fact=-1.
689	sp_zer=.true.
690	do kz=1,1
691	do jy=0,numygrid-1
692	do ix=0,numxgrid-1
693	if (grid(ix,jy,kz).gt.smallnum) then
694	if (sp_zer.eqv..true.) then ! first non zero value
695	sp_count_i=sp_count_i+1
696	sparse_dump_i(sp_count_i)= &
697	ix+jynumxgrid+kznumxgrid*numygrid
698	sp_zer=.false.
699	sp_fact=sp_fact*(-1.)
700	endif
701	sp_count_r=sp_count_r+1
702	sparse_dump_r(sp_count_r)= &
703	sp_fact* &
704	1.e12*grid(ix,jy,kz) &
705	/volume(ix,jy,kz)/outnum* &
706	weightair/weightmolar(ks)/densityoutgrid(ix,jy,kz)
707	sparse_dump_u(sp_count_r)= &
708	1.e12*gridsigma(ix,jy,kz)/volume(ix,jy,kz)/ &
709	outnum*weightair/weightmolar(ks)/ &
710	densityoutgrid(ix,jy,kz)
711	else ! concentration is zero
712	sp_zer=.true.
713	endif
714	end do
715	end do
716	end do
717	write(unitoutgridppt) sp_count_i
718	write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
719	write(unitoutgridppt) sp_count_r
720	write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
721	! write(unitoutgridppt) sp_count_r
722	! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
723
724	endif ! output for ppt
725
726	end do ! nageclass
727
728	close(unitoutgridppt)
729	close(unitoutgrid)
730
731	! itime is outside range
732	10 continue
733
734	end do ! maxpointspec_act
735
736	end do ! nspec
737
738	! RLT Aug 2017
739	! Write out conversion factor for dry air
740	inquire(file=path(2)(1:length(2))//'factor_drygrid',exist=lexist)
741	if (lexist.and..not.lstart) then
742	! open and append
743	open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
744	status='old',action='write',access='append')
745	else
746	! create new
747	open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
748	status='replace',action='write')
749	endif
750	sp_count_i=0
751	sp_count_r=0
752	sp_fact=-1.
753	sp_zer=.true.
754	do kz=1,1
755	do jy=0,numygrid-1
756	do ix=0,numxgrid-1
757	if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
758	if (sp_zer.eqv..true.) then ! first value not equal to one
759	sp_count_i=sp_count_i+1
760	sparse_dump_i(sp_count_i)= &
761	ix+jynumxgrid+kznumxgrid*numygrid
762	sp_zer=.false.
763	sp_fact=sp_fact*(-1.)
764	endif
765	sp_count_r=sp_count_r+1
766	sparse_dump_r(sp_count_r)= &
767	sp_fact*factor_drygrid(ix,jy,kz)
768	else ! factor is one
769	sp_zer=.true.
770	endif
771	end do
772	end do
773	end do
774	write(unitoutfactor) sp_count_i
775	write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
776	write(unitoutfactor) sp_count_r
777	write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
778	close(unitoutfactor)
779
780
781	if (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal
782	! if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ &
783	! wetgridtotal
784	! if (drygridtotal.gt.0.) drygridtotalunc=drygridsigmatotal/ &
785	! drygridtotal
786
787	! Dump of receptor concentrations
788
789	if (numreceptor.gt.0 .and. (iout.eq.2 .or. iout.eq.3) ) then
790	write(unitoutreceptppt) itime
791	do ks=1,nspec
792	write(unitoutreceptppt) (1.e12creceptor(i,ks)/outnum &
793	weightair/weightmolar(ks)/densityoutrecept(i),i=1,numreceptor)
794	end do
795	endif
796
797	! Dump of receptor concentrations
798
799	if (numreceptor.gt.0) then
800	write(unitoutrecept) itime
801	do ks=1,nspec
802	write(unitoutrecept) (1.e12*creceptor(i,ks)/outnum, &
803	i=1,numreceptor)
804	end do
805	endif
806
807	! RLT Aug 2017
808	! Write out conversion factor for dry air
809	if (numreceptor.gt.0) then
810	inquire(file=path(2)(1:length(2))//'factor_dryreceptor',exist=lexist)
811	if (lexist.and..not.lstart) then
812	! open and append
813	open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
814	status='old',action='write',access='append')
815	else
816	! create new
817	open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
818	status='replace',action='write')
819	endif
820	write(unitoutfactor) itime
821	write(unitoutfactor) (factor_dryrecept(i),i=1,numreceptor)
822	close(unitoutfactor)
823	endif
824
825	! reset lstart
826	if (lstart) then
827	lstart=.false.
828	endif
829	print*, 'after writing output files: lstart = ',lstart
830
831
832	! Reinitialization of grid
833	!*************************
834
835	do ks=1,nspec
836	do kp=1,maxpointspec_act
837	do i=1,numreceptor
838	creceptor(i,ks)=0.
839	end do
840	do jy=0,numygrid-1
841	do ix=0,numxgrid-1
842	do l=1,nclassunc
843	do nage=1,nageclass
844	do kz=1,numzgrid
845	gridunc(ix,jy,kz,ks,kp,l,nage)=0.
846	end do
847	end do
848	end do
849	end do
850	end do
851	end do
852	end do
853
854	end subroutine concoutput_inversion

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

Download in other formats:

Original Format