Context Navigation

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

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

Download in other formats:

Original Format