Context Navigation

nc_output_mod.f90 @ 9

Last change on this file since 9 was 9, checked in by hasod, 11 years ago
ADD: netcdf output module from Stephan Henne
File size: 51.7 KB

Line
1	!**********************************************************************
2	! Copyright 2013 *
3	! Dominik Brunner *
4	! *
5	! This file is part of FLEXPART-COSMO *
6	! *
7	! FLEXPART is free software: you can redistribute it and/or modify *
8	! it under the terms of the GNU General Public License as published by*
9	! the Free Software Foundation, either version 3 of the License, or *
10	! (at your option) any later version. *
11	! *
12	! FLEXPART is distributed in the hope that it will be useful, *
13	! but WITHOUT ANY WARRANTY; without even the implied warranty of *
14	! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15	! GNU General Public License for more details. *
16	! *
17	! You should have received a copy of the GNU General Public License *
18	! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
19	!**********************************************************************
20
21	!*****************************************************************************
22	! *
23	! This module handles all gridded netcdf output for concentration or *
24	! residence time and wet and dry deposition output. *
25	! *
26	! - writeheader_ncdf generates file including all information previously *
27	! stored in separate header files *
28	! - concoutput_ncdf write concentration output and wet and dry deposition *
29	! *
30	! Author: D. Brunner *
31	! *
32	! 12 April 2013 *
33	! *
34	!*****************************************************************************
35
36	#ifdef NETCDF_OUTPUT
37
38	module nc_output_mod
39
40	use netcdf
41
42	use point_mod, only: ireleasestart,ireleaseend,kindz,&
43	xpoint1,ypoint1,xpoint2,ypoint2,zpoint1,zpoint2,npart,xmass
44	use outg_mod, only: outheight,oroout,densityoutgrid,factor3d,volume,&
45	wetgrid,wetgridsigma,drygrid,drygridsigma,grid,gridsigma,&
46	area,arean,volumen, orooutn
47	use par_mod, only: dp, maxspec, maxreceptor, nclassunc,&
48	unitoutrecept,unitoutreceptppt, nxmax
49	use com_mod, only: path,length,ldirect,ibdate,ibtime,iedate,ietime, &
50	loutstep,loutaver,loutsample,outlon0,outlat0,&
51	numxgrid,numygrid,dxout,dyout,numzgrid, height, &
52	outlon0n,outlat0n,dxoutn,dyoutn,numxgridn,numygridn, &
53	nspec,maxpointspec_act,species,numpoint,&
54	dx,xlon0,dy,ylat0,compoint,method,lsubgrid,lconvection,&
55	ind_source,ind_receptor,nageclass,lage,&
56	drydep,wetdep,decay,weta,wetb, numbnests, &
57	reldiff,henry,f0,density,dquer,dsigma,dryvel,&
58	weightmolar,ohreact,spec_ass,kao,vsetaver,&
59	! for concoutput_ncdf and concoutput_nest_ncdf
60	nxmin1,nymin1,nz,oro,oron,rho,rhon,&
61	memind,xresoln,yresoln,xrn, xln, yrn,yln,nxn,nyn,&
62	xreceptor,yreceptor,numreceptor,creceptor,iout, &
63	itsplit, lsynctime, ctl, ifine, lagespectra, ipin, &
64	ioutputforeachrelease, iflux, mdomainfill, mquasilag, &
65	nested_output, ipout
66
67	implicit none
68
69	! parameter for data compression (1-9, 9 = most aggressive)
70	integer, parameter :: deflate_level = 9
71	logical, parameter :: min_size = .false. ! if set true, redundant fields (topography) are not written to minimize file size
72	character(len=255), parameter :: institution = 'Empa, Duebendorf, Switzerland'
73
74	integer :: tpointer
75	character(len=255) :: ncfname, ncfnamen
76
77	! netcdf dimension and variable IDs for main and nested output grid
78	integer, dimension(maxspec) :: specID,specIDppt, wdspecID,ddspecID
79	integer, dimension(maxspec) :: specIDn,specIDnppt, wdspecIDn,ddspecIDn
80	integer :: timeID, timeIDn
81	integer, dimension(6) :: dimids, dimidsn
82	integer, dimension(5) :: depdimids, depdimidsn
83	real,parameter :: eps=nxmax/3.e5
84
85	private:: writemetadata, output_units, nf90_err
86
87	contains
88
89	!****************************************************************
90	! determine output units (see table 1 in Stohl et al., ACP 2005
91	!****************************************************************
92	subroutine output_units(units)
93	character(len=15), intent(out) :: units
94	if (ldirect.eq.1) then
95	! forward simulation
96	if (ind_source.eq.1) then
97	if (ind_receptor.eq.1) then
98	units = 'ng m-3' ! hes the kg in Tab1 is only indicating the units of the relase not the output
99	else
100	units = 'ng kg-1'
101	endif
102	else
103	if (ind_receptor.eq.1) then
104	units = 'ng m-3'
105	else
106	units = 'ng kg-1'
107	endif
108	endif
109	else
110	! backward simulation
111	if (ind_source.eq.1) then
112	if (ind_receptor.eq.1) then
113	units = 's'
114	else
115	units = 's m3 kg-1'
116	endif
117	else
118	if (ind_receptor.eq.1) then
119	units = 's kg m-3'
120	else
121	units = 's'
122	endif
123	endif
124	endif
125	end subroutine output_units
126
127
128	!****************************************************************
129	! write metadata to netCDF file
130	!****************************************************************
131	subroutine writemetadata(ncid,lnest)
132
133	integer, intent(in) :: ncid
134	logical, intent(in) :: lnest
135	integer :: status
136	character :: time10,date8,adate8,atime6
137	character(5) :: zone
138	character(255) :: login_name, host_name
139
140	! gather system information
141	call date_and_time(date,time,zone)
142	call getlog(login_name)
143	call hostnm(host_name)
144
145	! hes CF convention requires these attributes
146	call nf90_err(nf90_put_att(ncid, nf90_global, 'Conventions', 'CF-1.6'))
147	call nf90_err(nf90_put_att(ncid, nf90_global, 'title', 'FLEXPART model output'))
148	! call nf90_err(nf90_put_att(ncid, nf90_global, 'institution', trim(institution)))
149	! call nf90_err(nf90_put_att(ncid, nf90_global, 'source', trim(version_str)//' model output'))
150	call nf90_err(nf90_put_att(ncid, nf90_global, 'history', date(1:4)//'-'//date(5:6)// &
151	'-'//date(7:8)//' '//time(1:2)//':'//time(3:4)//' '//zone//' created by '// &
152	trim(login_name)//' on '//trim(host_name)))
153	call nf90_err(nf90_put_att(ncid, nf90_global, 'references', &
154	'Stohl et al., Atmos. Chem. Phys., 2005, doi:10.5194/acp-5-2461-200'))
155
156	! attributes describing model run
157	!************************************************************************************
158
159	if (lnest) then
160	call nf90_err(nf90_put_att(ncid, nf90_global, 'outlon0', outlon0n))
161	call nf90_err(nf90_put_att(ncid, nf90_global, 'outlat0', outlat0n))
162	call nf90_err(nf90_put_att(ncid, nf90_global, 'dxout', dxoutn))
163	call nf90_err(nf90_put_att(ncid, nf90_global, 'dyout', dyoutn))
164	else
165	call nf90_err(nf90_put_att(ncid, nf90_global, 'outlon0', outlon0))
166	call nf90_err(nf90_put_att(ncid, nf90_global, 'outlat0', outlat0))
167	call nf90_err(nf90_put_att(ncid, nf90_global, 'dxout', dxout))
168	call nf90_err(nf90_put_att(ncid, nf90_global, 'dyout', dyout))
169	endif
170	! vertical levels stored in grid structure
171
172	! COMMAND file settings
173	call nf90_err(nf90_put_att(ncid, nf90_global, 'ldirect', ldirect))
174	write(adate,'(i8.8)') ibdate
175	write(atime,'(i6.6)') ibtime
176	call nf90_err(nf90_put_att(ncid, nf90_global, 'ibdate', adate))
177	call nf90_err(nf90_put_att(ncid, nf90_global, 'ibtime', atime))
178	write(adate,'(i8.8)') iedate
179	write(atime,'(i6.6)') ietime
180	call nf90_err(nf90_put_att(ncid, nf90_global, 'iedate', adate))
181	call nf90_err(nf90_put_att(ncid, nf90_global, 'ietime', atime))
182	call nf90_err(nf90_put_att(ncid, nf90_global, 'loutstep', loutstep))
183	call nf90_err(nf90_put_att(ncid, nf90_global, 'loutaver', loutaver))
184	call nf90_err(nf90_put_att(ncid, nf90_global, 'loutsample', loutsample))
185	call nf90_err(nf90_put_att(ncid, nf90_global, 'itsplit', itsplit))
186	call nf90_err(nf90_put_att(ncid, nf90_global, 'lsynctime', lsynctime))
187	call nf90_err(nf90_put_att(ncid, nf90_global, 'ctl', ctl))
188	call nf90_err(nf90_put_att(ncid, nf90_global, 'ifine', ifine))
189	call nf90_err(nf90_put_att(ncid, nf90_global, 'iout', iout))
190	call nf90_err(nf90_put_att(ncid, nf90_global, 'ipout', ipout))
191	call nf90_err(nf90_put_att(ncid, nf90_global, 'lsubgrid', lsubgrid))
192	call nf90_err(nf90_put_att(ncid, nf90_global, 'lconvection', lconvection))
193	call nf90_err(nf90_put_att(ncid, nf90_global, 'lagespectra', lagespectra))
194	call nf90_err(nf90_put_att(ncid, nf90_global, 'ipin', ipin))
195	call nf90_err(nf90_put_att(ncid, nf90_global, 'ioutputforeachrelease', ioutputforeachrelease))
196	call nf90_err(nf90_put_att(ncid, nf90_global, 'iflux', iflux))
197	call nf90_err(nf90_put_att(ncid, nf90_global, 'mdomainfill', mdomainfill))
198	call nf90_err(nf90_put_att(ncid, nf90_global, 'ind_source', ind_source))
199	call nf90_err(nf90_put_att(ncid, nf90_global, 'ind_receptor', ind_receptor))
200	call nf90_err(nf90_put_att(ncid, nf90_global, 'mquasilag', mquasilag))
201	call nf90_err(nf90_put_att(ncid, nf90_global, 'nested_output', nested_output))
202
203	end subroutine writemetadata
204
205
206	!****************************************************************
207	! netcdf error message handling
208	!****************************************************************
209	subroutine nf90_err(status)
210	integer, intent (in) :: status
211	if(status /= nf90_noerr) then
212	print *, trim(nf90_strerror(status))
213	stop 'Stopped'
214	end if
215	end subroutine nf90_err
216
217
218	!****************************************************************
219	! Create netcdf file and write header/metadata information
220	! lnest = .false. : Create main output file
221	! lnest = .true. : Create nested output file
222	!****************************************************************
223	subroutine writeheader_ncdf(lnest)
224
225	implicit none
226
227	logical, intent(in) :: lnest
228
229	integer :: ncid, sID, wdsID, ddsID
230	integer :: timeDimID, latDimID, lonDimID, levDimID
231	integer :: nspecDimID, npointDimID, nageclassDimID, ncharDimID
232	integer :: tID, lonID, latID, levID, poleID, lageID, oroID
233	integer :: rellng1ID, rellng2ID, rellat1ID, rellat2ID, relzz1ID, relzz2ID
234	integer :: relcomID, relkindzID, relstartID, relendID, relpartID, relxmassID
235	integer :: nnx, nny
236	integer, dimension(6) :: dIDs
237	integer, dimension(5) :: depdIDs
238	character(len=255) :: fname
239	character(len=15) :: units
240	character(len=10) :: fprefix
241	character(len=3) :: anspec
242	character :: adate8,atime6,timeunit*32
243	real, dimension(1000) :: coord
244
245	integer :: cache_size
246	integer, dimension(6) :: chunksizes
247	integer, dimension(5) :: dep_chunksizes
248
249	integer :: i,ix,jy
250
251	!************************
252	! Create netcdf file
253	!************************
254
255	if (ldirect.eq.1) then
256	write(adate,'(i8.8)') ibdate
257	write(atime,'(i6.6)') ibtime
258	fprefix = 'grid_conc_'
259	else
260	write(adate,'(i8.8)') iedate
261	write(atime,'(i6.6)') ietime
262	fprefix = 'grid_time_'
263	endif
264
265	if (lnest) then
266	fname = path(2)(1:length(2))//fprefix//adate//atime//'_nest.nc'
267	ncfnamen = fname
268	nnx = numxgridn
269	nny = numygridn
270	else
271	fname = path(2)(1:length(2))//fprefix//adate//atime//'.nc'
272	ncfname = fname
273	nnx = numxgrid
274	nny = numygrid
275	endif
276
277	cache_size = 16 * nnx * nny * numzgrid
278
279	! setting cache size in bytes. It is set to 4 times the largest data block that is written
280	! size_type x nx x ny x nz
281	! create file
282
283	call nf90_err(nf90_create(trim(fname), cmode = nf90_hdf5, ncid = ncid, &
284	cache_size = cache_size))
285
286	! create dimensions:
287	!*************************
288	! time
289	call nf90_err(nf90_def_dim(ncid, 'time', nf90_unlimited, timeDimID))
290	timeunit = 'seconds since '//adate(1:4)//'-'//adate(5:6)// &
291	'-'//adate(7:8)//' '//atime(1:2)//':'//atime(3:4)
292
293	! lon
294	call nf90_err(nf90_def_dim(ncid, 'lon', nnx, lonDimID))
295	! lat
296	call nf90_err(nf90_def_dim(ncid, 'lat', nny, latDimID))
297
298	! level
299	call nf90_err(nf90_def_dim(ncid, 'level', numzgrid, levDimID))
300	! number of species
301	call nf90_err(nf90_def_dim(ncid, 'numspec', nspec, nspecDimID))
302	! number of release points
303	call nf90_err(nf90_def_dim(ncid, 'numpoint', numpoint, npointDimID))
304	! number of age classes
305	call nf90_err(nf90_def_dim(ncid, 'nageclass', nageclass, nageclassDimID))
306	! dimension for release point characters
307	call nf90_err(nf90_def_dim(ncid, 'nchar', 45, ncharDimID))
308
309	! create variables
310	!*************************
311
312	! time
313	call nf90_err(nf90_def_var(ncid, 'time', nf90_int, (/ timeDimID /), tID))
314	call nf90_err(nf90_put_att(ncid, tID, 'units', timeunit))
315	call nf90_err(nf90_put_att(ncid, tID, 'calendar', 'proleptic_gregorian'))
316	if (lnest) then
317	timeIDn = tID
318	else
319	timeID = tID
320	endif
321
322	! lon
323	call nf90_err(nf90_def_var(ncid, 'lon', nf90_float, (/ lonDimID /), lonID))
324	call nf90_err(nf90_put_att(ncid, lonID, 'long_name', 'longitude in degree east'))
325
326	call nf90_err(nf90_put_att(ncid, lonID, 'axis', 'X'))
327	call nf90_err(nf90_put_att(ncid, lonID, 'units', 'degrees'))
328	call nf90_err(nf90_put_att(ncid, lonID, 'standard_name', 'grid_longitude'))
329	call nf90_err(nf90_put_att(ncid, lonID, 'description', 'grid cell centers'))
330	! lat
331	call nf90_err(nf90_def_var(ncid, 'lat', nf90_float, (/ latDimID /), latID))
332	call nf90_err(nf90_put_att(ncid, latID, 'long_name', 'latitude in degree north'))
333
334	call nf90_err(nf90_put_att(ncid, latID, 'axis', 'Y'))
335	call nf90_err(nf90_put_att(ncid, latID, 'units', 'degrees'))
336	call nf90_err(nf90_put_att(ncid, latID, 'standard_name', 'grid_latitude'))
337	call nf90_err(nf90_put_att(ncid, latID, 'description', 'grid cell centers'))
338
339
340	! level
341	call nf90_err(nf90_def_var(ncid, 'level', nf90_float, (/ levDimID /), levID))
342	call nf90_err(nf90_put_att(ncid, levID, 'axis', 'Z'))
343	call nf90_err(nf90_put_att(ncid, levID, 'units', 'meters'))
344	call nf90_err(nf90_put_att(ncid, levID, 'positive', 'up'))
345	call nf90_err(nf90_put_att(ncid, levID, 'standard_name', 'height'))
346	call nf90_err(nf90_put_att(ncid, levID, 'long_name', 'height above ground'))
347
348	! release comment
349	call nf90_err(nf90_def_var(ncid, 'RELCOM', nf90_char, (/ ncharDimID,npointDimID /), &
350	relcomID))
351	call nf90_err(nf90_put_att(ncid, relcomID, 'long_name', 'release point name'))
352	! release longitude 1
353	call nf90_err(nf90_def_var(ncid, 'RELLNG1', nf90_float, (/ npointDimID /), rellng1ID))
354	call nf90_err(nf90_put_att(ncid, rellng1ID, 'units', 'degrees_east'))
355	call nf90_err(nf90_put_att(ncid, rellng1ID, 'long_name', &
356	'release longitude lower left corner'))
357	! release longitude 2
358	call nf90_err(nf90_def_var(ncid, 'RELLNG2', nf90_float, (/ npointDimID /), rellng2ID))
359	call nf90_err(nf90_put_att(ncid, rellng2ID, 'units', 'degrees_east'))
360	call nf90_err(nf90_put_att(ncid, rellng2ID, 'long_name', &
361	'release longitude upper right corner'))
362	! release latitude 1
363	call nf90_err(nf90_def_var(ncid, 'RELLAT1', nf90_float, (/ npointDimID /), rellat1ID))
364	call nf90_err(nf90_put_att(ncid, rellat1ID, 'units', 'degrees_north'))
365	call nf90_err(nf90_put_att(ncid, rellat1ID, 'long_name', &
366	'release latitude lower left corner'))
367	! release latitude 2
368	call nf90_err(nf90_def_var(ncid, 'RELLAT2', nf90_float, (/ npointDimID /), rellat2ID))
369	call nf90_err(nf90_put_att(ncid, rellat2ID, 'units', 'degrees_north'))
370	call nf90_err(nf90_put_att(ncid, rellat2ID, 'long_name', &
371	'release latitude upper right corner'))
372
373	! hes: if rotated_ll it would be convenient also to write the the release points in rotated_coordinates
374
375	! release height bottom
376	call nf90_err(nf90_def_var(ncid, 'RELZZ1', nf90_float, (/ npointDimID /), relzz1ID))
377	call nf90_err(nf90_put_att(ncid, relzz1ID, 'units', 'meters'))
378	call nf90_err(nf90_put_att(ncid, relzz1ID, 'long_name', 'release height bottom'))
379	! release height top
380	call nf90_err(nf90_def_var(ncid, 'RELZZ2', nf90_float, (/ npointDimID /), relzz2ID))
381	call nf90_err(nf90_put_att(ncid, relzz2ID, 'units', 'meters'))
382	call nf90_err(nf90_put_att(ncid, relzz2ID, 'long_name', 'release height top'))
383	! release kind
384	call nf90_err(nf90_def_var(ncid, 'RELKINDZ', nf90_int, (/ npointDimID /), relkindzID))
385	call nf90_err(nf90_put_att(ncid, relkindzID, 'long_name', 'release kind'))
386	! release start
387	call nf90_err(nf90_def_var(ncid, 'RELSTART', nf90_int, (/ npointDimID /), relstartID))
388	call nf90_err(nf90_put_att(ncid, relstartID, 'units', 'seconds'))
389	call nf90_err(nf90_put_att(ncid, relstartID, 'long_name', &
390	'release start relative to simulation start'))
391	! release end
392	call nf90_err(nf90_def_var(ncid, 'RELEND', nf90_int, (/ npointDimID /), relendID))
393	call nf90_err(nf90_put_att(ncid, relendID, 'units', 'seconds'))
394	call nf90_err(nf90_put_att(ncid, relendID, 'long_name', &
395	'release end relative to simulation start'))
396	! release particles
397	call nf90_err(nf90_def_var(ncid, 'RELPART', nf90_int, (/ npointDimID /), relpartID))
398	call nf90_err(nf90_put_att(ncid, relpartID, 'long_name', 'number of release particles'))
399	! release particle masses
400	call nf90_err(nf90_def_var(ncid, 'RELXMASS', nf90_float, (/ npointDimID, nspecDimID /), &
401	relxmassID))
402	call nf90_err(nf90_put_att(ncid, relxmassID, 'long_name', 'total release particle mass'))
403
404	! age classes
405	call nf90_err(nf90_def_var(ncid, 'LAGE', nf90_int, (/ nageclassDimID /), lageID))
406	call nf90_err(nf90_put_att(ncid, lageID, 'units', 'seconds'))
407	call nf90_err(nf90_put_att(ncid, lageID, 'long_name', 'age class'))
408
409	! output orography
410	if (.not. min_size) then
411	call nf90_err(nf90_def_var(ncid, 'ORO', nf90_int, (/ lonDimID, latDimID /), oroID, &
412	deflate_level=deflate_level, chunksizes= (/ nnx, nny /)))
413	call nf90_err(nf90_put_att(ncid, oroID, 'standard_name', 'surface altitude'))
414	call nf90_err(nf90_put_att(ncid, oroID, 'long_name', 'outgrid surface altitude'))
415	call nf90_err(nf90_put_att(ncid, oroID, 'units', 'm'))
416	end if
417
418
419
420	! concentration output, wet and dry deposition variables (one per species)
421	call output_units(units)
422
423	dIDs = (/ londimid, latdimid, levdimid, timedimid, npointdimid, nageclassdimid /)
424	depdIDs = (/ londimid, latdimid, timedimid, npointdimid, nageclassdimid /)
425	if (lnest) then
426	dimidsn = dIDs
427	depdimidsn = depdIDs
428	else
429	dimids = dIDs
430	depdimids = depdIDs
431	endif
432
433	! set chunksizes according to largest written portion of data in an individual call to
434	! nf90_put_var
435	chunksizes = (/ nnx, nny, numzgrid, 1, 1, 1 /)
436	dep_chunksizes = (/ nnx, nny, 1, 1, 1 /)
437
438	do i = 1,nspec
439	write(anspec,'(i3.3)') i
440
441	! concentration output
442	if (iout.eq.1.or.iout.eq.3.or.iout.eq.5) then
443	call nf90_err(nf90_def_var(ncid,'spec'//anspec, nf90_float, dIDs, sID , &
444	deflate_level = deflate_level, &
445	chunksizes = chunksizes ))
446	call nf90_err(nf90_put_att(ncid, sID, 'units', units))
447	call nf90_err(nf90_put_att(ncid, sID, 'long_name', species(i)))
448	call nf90_err(nf90_put_att(ncid, sID, 'decay', decay(i)))
449	call nf90_err(nf90_put_att(ncid, sID, 'weightmolar', weightmolar(i)))
450	call nf90_err(nf90_put_att(ncid, sID, 'ohreact', ohreact(i)))
451	!call nf90_err(nf90_put_att(ncid, sID, 'kao', kao(i)) ! what is this??
452	call nf90_err(nf90_put_att(ncid, sID, 'vsetaver', vsetaver(i)))
453	call nf90_err(nf90_put_att(ncid, sID, 'spec_ass', spec_ass(i)))
454
455	if (lnest) then
456	specIDn(i) = sID
457	else
458	specID(i) = sID
459	endif
460	endif
461
462	! mixing ratio output
463	if (iout.eq.2.or.iout.eq.3) then
464	call nf90_err(nf90_def_var(ncid,'spec'//anspec//'_pptv', nf90_float, dIDs, sID , &
465	deflate_level = deflate_level, &
466	chunksizes = chunksizes ))
467	call nf90_err(nf90_put_att(ncid, sID, 'units', 'pptv'))
468	call nf90_err(nf90_put_att(ncid, sID, 'long_name', species(i)))
469	call nf90_err(nf90_put_att(ncid, sID, 'decay', decay(i)))
470	call nf90_err(nf90_put_att(ncid, sID, 'weightmolar', weightmolar(i)))
471	call nf90_err(nf90_put_att(ncid, sID, 'ohreact', ohreact(i)))
472	!call nf90_err(nf90_put_att(ncid, sID, 'kao', kao(i)) ! what is this??
473	call nf90_err(nf90_put_att(ncid, sID, 'vsetaver', vsetaver(i)))
474	call nf90_err(nf90_put_att(ncid, sID, 'spec_ass', spec_ass(i)))
475
476	if (lnest) then
477	specIDnppt(i) = sID
478	else
479	specIDppt(i) = sID
480	endif
481	endif
482
483	! wet and dry deposition fields for forward runs
484	if (wetdep) then
485	call nf90_err(nf90_def_var(ncid,'WD_spec'//anspec, nf90_float, depdIDs, &
486	wdsID, deflate_level = deflate_level, &
487	chunksizes = dep_chunksizes))
488	call nf90_err(nf90_put_att(ncid, wdsID, 'units', '1e-12 kg m-2'))
489	call nf90_err(nf90_put_att(ncid, wdsID, 'weta', weta(i)))
490	call nf90_err(nf90_put_att(ncid, wdsID, 'wetb', wetb(i)))
491	call nf90_err(nf90_put_att(ncid, wdsID, 'dquer', dquer(i)))
492	call nf90_err(nf90_put_att(ncid, wdsID, 'henry', henry(i)))
493	if (lnest) then
494	wdspecIDn(i) = wdsID
495	else
496	wdspecID(i) = wdsID
497	endif
498	endif
499	if (drydep) then
500	call nf90_err(nf90_def_var(ncid,'DD_spec'//anspec, nf90_float, depdIDs, &
501	ddsID, deflate_level = deflate_level, &
502	chunksizes = dep_chunksizes))
503	call nf90_err(nf90_put_att(ncid, ddsID, 'units', '1e-12 kg m-2'))
504	call nf90_err(nf90_put_att(ncid, ddsID, 'dryvel', dryvel(i)))
505	call nf90_err(nf90_put_att(ncid, ddsID, 'reldiff', reldiff(i)))
506	call nf90_err(nf90_put_att(ncid, ddsID, 'henry', henry(i)))
507	call nf90_err(nf90_put_att(ncid, ddsID, 'f0', f0(i)))
508	call nf90_err(nf90_put_att(ncid, ddsID, 'dquer', dquer(i)))
509	call nf90_err(nf90_put_att(ncid, ddsID, 'density', density(i)))
510	call nf90_err(nf90_put_att(ncid, ddsID, 'dsigma', dsigma(i)))
511	if (lnest) then
512	ddspecIDn(i) = ddsID
513	else
514	ddspecID(i) = ddsID
515	endif
516	endif
517	end do
518
519
520	! global (metadata) attributes
521	!*******************************
522	call writemetadata(ncid,lnest)
523
524
525	! moves the file from define to data mode
526	call nf90_err(nf90_enddef(ncid))
527
528	! ! hes: inquire var definition
529	! do i = 1,nspec
530	! write(anspec,'(i3.3)') i
531	!
532	! ! concentration output
533	! if (iout.eq.1.or.iout.eq.3.or.iout.eq.5) then
534	! if (lnest) then
535	! sID = specIDn(i)
536	! else
537	! sID = specID(i)
538	! endif
539	! call nf90_err(nf90_inquire_variable(ncid, sID, chunksizes=inq_chunksizes))
540	! write(,) "Chunksizes for var "//anspec//": ", inq_chunksizes
541	! endif
542	! end do
543
544
545	! fill with data
546	!******************************
547	! longitudes (grid cell centers)
548	if (lnest) then
549	do i = 1,numxgridn
550	coord(i) = outlon0n + (i-0.5)*dxoutn
551	enddo
552	call nf90_err(nf90_put_var(ncid, lonID, coord(1:numxgridn)))
553	else
554	do i = 1,numxgrid
555	coord(i) = outlon0 + (i-0.5)*dxout
556	enddo
557	call nf90_err(nf90_put_var(ncid, lonID, coord(1:numxgrid)))
558	endif
559	! latitudes (grid cell centers)
560	if (lnest) then
561	do i = 1,numygridn
562	coord(i) = outlat0n + (i-0.5)*dyoutn
563	enddo
564	call nf90_err(nf90_put_var(ncid, latID, coord(1:numygridn)))
565	else
566	do i = 1,numygrid
567	coord(i) = outlat0 + (i-0.5)*dyout
568	enddo
569	call nf90_err(nf90_put_var(ncid, latID, coord(1:numygrid)))
570	endif
571	! levels
572	call nf90_err(nf90_put_var(ncid, levID, outheight(1:numzgrid)))
573
574	! release point information
575	do i = 1,numpoint
576	call nf90_err(nf90_put_var(ncid, relstartID, ireleasestart(i), (/i/)))
577	call nf90_err(nf90_put_var(ncid, relendID, ireleaseend(i), (/i/)))
578	call nf90_err(nf90_put_var(ncid, relkindzID, kindz(i), (/i/)))
579	call nf90_err(nf90_put_var(ncid, rellng1ID, xpoint1(i), (/i/)))
580	call nf90_err(nf90_put_var(ncid, rellng2ID, xpoint2(i), (/i/)))
581	call nf90_err(nf90_put_var(ncid, rellat1ID, ypoint1(i), (/i/)))
582	call nf90_err(nf90_put_var(ncid, rellat2ID, ypoint2(i), (/i/)))
583	call nf90_err(nf90_put_var(ncid, relzz1ID, zpoint1(i), (/i/)))
584	call nf90_err(nf90_put_var(ncid, relzz2ID, zpoint2(i), (/i/)))
585	call nf90_err(nf90_put_var(ncid, relpartID, npart(i), (/i/)))
586	if (i .le. 1000) then
587	call nf90_err(nf90_put_var(ncid, relcomID, compoint(i), (/1,i/), (/45,1/)))
588	else
589	call nf90_err(nf90_put_var(ncid, relcomID, 'NA', (/1,i/), (/45,1/)))
590	endif
591	call nf90_err(nf90_put_var(ncid, relxmassID, xmass(i,1:nspec), (/i,1/), (/1,nspec/)))
592	end do
593
594	! age classes
595	call nf90_err(nf90_put_var(ncid, lageID, lage(1:nageclass)))
596
597	! orography
598	if (.not. min_size) then
599	if (lnest) then
600	call nf90_err(nf90_put_var(ncid, oroID, orooutn(0:(nnx-1), 0:(nny-1))))
601	else
602	call nf90_err(nf90_put_var(ncid, oroID, oroout(0:(nnx-1), 0:(nny-1))))
603	endif
604	end if
605
606	call nf90_err(nf90_close(ncid))
607
608	return
609
610	end subroutine writeheader_ncdf
611
612
613	subroutine concoutput_ncdf(itime,outnum,gridtotalunc,wetgridtotalunc, &
614	drygridtotalunc)
615	! i i o o
616	! o
617	!*****************************************************************************
618	! *
619	! Output of the concentration grid and the receptor concentrations. *
620	! *
621	! Author: A. Stohl *
622	! *
623	! 24 May 1995 *
624	! *
625	! 13 April 1999, Major update: if output size is smaller, dump output in *
626	! sparse matrix format; additional output of uncertainty *
627	! *
628	! 05 April 2000, Major update: output of age classes; output for backward*
629	! runs is time spent in grid cell times total mass of *
630	! species. *
631	! *
632	! 17 February 2002, Appropriate dimensions for backward and forward runs *
633	! are now specified in module par_mod *
634	! *
635	! June 2006, write grid in sparse matrix with a single write command *
636	! in order to save disk space *
637	! *
638	! 2008 new sparse matrix format *
639	! *
640	! February 2010, Dominik Brunner, Empa *
641	! Adapted for COSMO *
642	! Remark: calculation of density could be improved. *
643	! Currently, it is calculated for the lower left corner *
644	! of each output grid cell rather than for its center. *
645	! Furthermore, the average density could be calculated *
646	! from the difference in pressure at the top and bottom *
647	! of each cell rather than by interpolation. *
648	! *
649	! April 2013, Dominik Brunner, Empa *
650	! Adapted for netcdf output *
651	! *
652	!*****************************************************************************
653	! *
654	! Variables: *
655	! outnum number of samples *
656	! ncells number of cells with non-zero concentrations *
657	! sparse .true. if in sparse matrix format, else .false. *
658	! tot_mu 1 for forward, initial mass mixing ration for backw. runs *
659	! *
660	!*****************************************************************************
661
662	use unc_mod, only: gridunc,drygridunc,wetgridunc
663
664
665	implicit none
666
667	integer, intent(in) :: itime
668	real, intent(in) :: outnum
669	real, intent(out) :: gridtotalunc,wetgridtotalunc,drygridtotalunc
670	real :: densityoutrecept(maxreceptor)
671	integer :: ncid,kp,ks,kz,ix,jy,iix,jjy,kzz,kzzm1,ngrid
672	integer :: nage,i,l, jj
673	real :: tot_mu(maxspec,maxpointspec_act)
674	real :: halfheight,dz,dz1,dz2
675	real :: xl,yl,xlrot,ylrot,zagnd,zagndprev
676	real :: auxgrid(nclassunc),gridtotal,gridsigmatotal
677	real :: wetgridtotal,wetgridsigmatotal
678	real :: drygridtotal,drygridsigmatotal
679
680	real, parameter :: weightair=28.97
681
682
683	! open output file
684	call nf90_err(nf90_open(trim(ncfname), nf90_write, ncid))
685
686	! write time
687	tpointer = tpointer + 1
688	call nf90_err(nf90_put_var( ncid, timeID, itime, (/ tpointer /)))
689
690	! For forward simulations, output fields have dimension MAXSPEC,
691	! for backward simulations, output fields have dimension MAXPOINT.
692	! Thus, make loops either about nspec, or about numpoint
693	!*****************************************************************
694
695	if (ldirect.eq.1) then
696	do ks=1,nspec
697	do kp=1,maxpointspec_act
698	tot_mu(ks,kp)=1.0
699	end do
700	end do
701	else
702	do ks=1,nspec
703	do kp=1,maxpointspec_act
704	tot_mu(ks,kp)=xmass(kp,ks)
705	end do
706	end do
707	endif
708
709
710	!*******************************************************************
711	! Compute air density:
712	! brd134: we now take into account whether we are in the mother or in
713	! a nested domain (before only from mother domain)
714	! Determine center altitude of output layer, and interpolate density
715	! data to that altitude
716	!*******************************************************************
717
718	do kz=1,numzgrid
719	if (kz.eq.1) then
720	halfheight=outheight(1)/2.
721	else
722	halfheight=(outheight(kz)+outheight(kz-1))/2.
723	endif
724	do kzz=2,nz
725	if ((height(kzz-1).lt.halfheight).and. &
726	(height(kzz).gt.halfheight)) exit
727	end do
728	kzz=max(min(kzz,nz),2)
729	dz1=halfheight-height(kzz-1)
730	dz2=height(kzz)-halfheight
731	dz=dz1+dz2
732
733	do jy=0,numygrid-1
734	do ix=0,numxgrid-1
735	xl=outlon0+real(ix)*dxout
736	yl=outlat0+real(jy)*dyout
737	! grid index in mother domain
738	xl=(xl-xlon0)/dx
739	yl=(yl-ylat0)/dx
740
741	ngrid=0
742	do jj=numbnests,1,-1
743	if ( xl.gt.xln(jj)+eps .and. xl.lt.xrn(jj)-eps .and. &
744	yl.gt.yln(jj)+eps .and. yl.lt.yrn(jj)-eps ) then
745	ngrid=jj
746	exit
747	end if
748	end do
749
750	if (ngrid.eq.0) then
751	iix=max(min(nint(xl),nxmin1),0) ! if output grid cell is outside mother domain
752	jjy=max(min(nint(yl),nymin1),0)
753
754	densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,memind(2))*dz1+ &
755	rho(iix,jjy,kzz-1,memind(2))*dz2)/dz
756	else
757	xl=(xl-xln(ngrid))*xresoln(ngrid)
758	yl=(yl-yln(ngrid))*yresoln(ngrid)
759	iix=max(min(nint(xl),nxn(ngrid)-1),0)
760	jjy=max(min(nint(yl),nyn(ngrid)-1),0)
761
762	densityoutgrid(ix,jy,kz)=(rhon(iix,jjy,kzz,memind(2), ngrid)*dz1+ &
763	rhon(iix,jjy,kzz-1,memind(2), ngrid)*dz2)/dz
764	endif
765	end do
766	end do
767	end do
768
769	! brd134: for receptor points no option for nests yet to specify density
770	! and also altitude zreceptor not considered yet (needs revision)
771	do i=1,numreceptor
772	xl=xreceptor(i)
773	yl=yreceptor(i)
774	iix=max(min(nint(xl),nxmin1),0)
775	jjy=max(min(nint(yl),nymin1),0)
776	densityoutrecept(i)=rho(iix,jjy,1,memind(2))
777	end do
778
779	! Output is different for forward and backward simulations
780	if (ldirect.eq.1) then
781	do kz=1,numzgrid
782	do jy=0,numygrid-1
783	do ix=0,numxgrid-1
784	factor3d(ix,jy,kz)=1.e12/volume(ix,jy,kz)/outnum
785	end do
786	end do
787	end do
788	else
789	do kz=1,numzgrid
790	do jy=0,numygrid-1
791	do ix=0,numxgrid-1
792	factor3d(ix,jy,kz)=real(abs(loutaver))/outnum
793	end do
794	end do
795	end do
796	endif
797
798	!*********************************************************************
799	! Determine the standard deviation of the mean concentration or mixing
800	! ratio (uncertainty of the output) and the dry and wet deposition
801	!*********************************************************************
802
803	gridtotal=0.
804	gridsigmatotal=0.
805	gridtotalunc=0.
806	wetgridtotal=0.
807	wetgridsigmatotal=0.
808	wetgridtotalunc=0.
809	drygridtotal=0.
810	drygridsigmatotal=0.
811	drygridtotalunc=0.
812
813	do ks=1,nspec
814
815	do kp=1,maxpointspec_act
816	do nage=1,nageclass
817
818	do jy=0,numygrid-1
819	do ix=0,numxgrid-1
820
821	! WET DEPOSITION
822	if ((wetdep).and.(ldirect.gt.0)) then
823	do l=1,nclassunc
824	auxgrid(l)=wetgridunc(ix,jy,ks,kp,l,nage)
825	end do
826	call mean(auxgrid,wetgrid(ix,jy), &
827	wetgridsigma(ix,jy),nclassunc)
828	! Multiply by number of classes to get total concentration
829	wetgrid(ix,jy)=wetgrid(ix,jy)*real(nclassunc)
830	wetgridtotal=wetgridtotal+wetgrid(ix,jy)
831	! Calculate standard deviation of the mean
832	wetgridsigma(ix,jy)= &
833	wetgridsigma(ix,jy)* &
834	sqrt(real(nclassunc))
835	wetgridsigmatotal=wetgridsigmatotal+ &
836	wetgridsigma(ix,jy)
837	endif
838
839	! DRY DEPOSITION
840	if ((drydep).and.(ldirect.gt.0)) then
841	do l=1,nclassunc
842	auxgrid(l)=drygridunc(ix,jy,ks,kp,l,nage)
843	end do
844	call mean(auxgrid,drygrid(ix,jy), &
845	drygridsigma(ix,jy),nclassunc)
846	! Multiply by number of classes to get total concentration
847	drygrid(ix,jy)=drygrid(ix,jy)*real(nclassunc)
848	drygridtotal=drygridtotal+drygrid(ix,jy)
849	! Calculate standard deviation of the mean
850	drygridsigma(ix,jy)= &
851	drygridsigma(ix,jy)* &
852	sqrt(real(nclassunc))
853	drygridsigmatotal=drygridsigmatotal+ &
854	drygridsigma(ix,jy)
855	endif
856
857	! CONCENTRATION OR MIXING RATIO
858	do kz=1,numzgrid
859	do l=1,nclassunc
860	auxgrid(l)=gridunc(ix,jy,kz,ks,kp,l,nage)
861	end do
862	call mean(auxgrid,grid(ix,jy,kz), &
863	gridsigma(ix,jy,kz),nclassunc)
864	! Multiply by number of classes to get total concentration
865	grid(ix,jy,kz)= &
866	grid(ix,jy,kz)*real(nclassunc)
867	gridtotal=gridtotal+grid(ix,jy,kz)
868	! Calculate standard deviation of the mean
869	gridsigma(ix,jy,kz)= &
870	gridsigma(ix,jy,kz)* &
871	sqrt(real(nclassunc))
872	gridsigmatotal=gridsigmatotal+ &
873	gridsigma(ix,jy,kz)
874	end do
875	end do
876	end do
877
878
879	!*******************************************************************
880	! Generate output: may be in concentration (ng/m3) or in mixing
881	! ratio (ppt) or both
882	! Output the position and the values alternated multiplied by
883	! 1 or -1, first line is number of values, number of positions
884	! For backward simulations, the unit is seconds, stored in grid_time
885	!*******************************************************************
886
887	! Concentration output
888	!*********************
889	if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
890
891	! Wet deposition
892	if ((ldirect.eq.1).and.(WETDEP)) then
893	call nf90_err(nf90_put_var(ncid,wdspecID(ks),1.e12*&
894	wetgrid(0:numxgrid-1,0:numygrid-1)/area(0:numxgrid-1,0:numygrid-1),&
895	(/ 1,1,tpointer,kp,nage /), (/ numxgrid,numygrid,1,1,1 /)))
896	endif
897
898	! Dry deposition
899	if ((ldirect.eq.1).and.(DRYDEP)) then
900	call nf90_err(nf90_put_var(ncid,ddspecID(ks),1.e12*&
901	drygrid(0:numxgrid-1,0:numygrid-1)/area(0:numxgrid-1,0:numygrid-1),&
902	(/ 1,1,tpointer,kp,nage /), (/ numxgrid,numygrid,1,1,1 /)))
903	endif
904
905	! Concentrations
906	call nf90_err(nf90_put_var(ncid,specID(ks),grid(0:numxgrid-1,0:numygrid-1,&
907	1:numzgrid)*factor3d(0:numxgrid-1,0:numygrid-1,1:numzgrid)/tot_mu(ks,kp),&
908	(/ 1,1,1,tpointer,kp,nage /), (/ numxgrid,numygrid,numzgrid,1,1,1 /) ))
909
910	endif ! concentration output
911
912	! Mixing ratio output
913	!********************
914
915	if ((iout.eq.2)) then ! mixing ratio
916
917	! Wet deposition
918	if ((ldirect.eq.1).and.(WETDEP)) then
919	call nf90_err(nf90_put_var(ncid,wdspecID(ks),1.e12*&
920	wetgrid(0:numxgrid-1,0:numygrid-1)/area(0:numxgrid-1,0:numygrid-1),&
921	(/ 1,1,tpointer,kp,nage /), (/ numxgrid,numygrid,1,1,1 /)))
922	endif
923
924	! Dry deposition
925	if ((ldirect.eq.1).and.(DRYDEP)) then
926	call nf90_err(nf90_put_var(ncid,ddspecID(ks),1.e12*&
927	drygrid(0:numxgrid-1,0:numygrid-1)/area(0:numxgrid-1,0:numygrid-1),&
928	(/ 1,1,tpointer,kp,nage /), (/ numxgrid,numygrid,1,1,1 /)))
929	endif
930
931	! Mixing ratios
932	call nf90_err(nf90_put_var(ncid,specIDppt(ks),weightair/weightmolar(ks)*&
933	grid(0:numxgrid-1,0:numygrid-1,1:numzgrid)*&
934	factor3d(0:numxgrid-1,0:numygrid-1,1:numzgrid)/&
935	densityoutgrid(0:numxgrid-1,0:numygrid-1,1:numzgrid),&
936	(/ 1,1,1,tpointer,kp,nage /), (/ numxgrid,numygrid,numzgrid,1,1,1 /)))
937
938	endif ! output for ppt
939
940	end do
941	end do
942
943	end do
944
945	! Close netCDF file
946	!**************************
947	call nf90_err(nf90_close(ncid))
948
949
950	if (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal
951	if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ &
952	wetgridtotal
953	if (drygridtotal.gt.0.) drygridtotalunc=drygridsigmatotal/ &
954	drygridtotal
955
956	! Dump of receptor concentrations
957
958	if (numreceptor.gt.0 .and. (iout.eq.2 .or. iout.eq.3) ) then
959	write(unitoutreceptppt) itime
960	do ks=1,nspec
961	write(unitoutreceptppt) (1.e12creceptor(i,ks)/outnum &
962	weightair/weightmolar(ks)/densityoutrecept(i),i=1,numreceptor)
963	end do
964	endif
965
966	! Dump of receptor concentrations
967
968	if (numreceptor.gt.0) then
969	write(unitoutrecept) itime
970	do ks=1,nspec
971	write(unitoutrecept) (1.e12*creceptor(i,ks)/outnum,i=1,numreceptor)
972	end do
973	endif
974
975
976	! Reinitialization of grid
977	!*************************
978
979	creceptor(1:numreceptor,1:nspec) = 0.
980	gridunc(:,:,:,1:nspec,:,:,1:nageclass) = 0.
981
982
983	end subroutine concoutput_ncdf
984
985
986	subroutine concoutput_nest_ncdf(itime,outnum)
987	! i i
988	!*****************************************************************************
989	! *
990	! Output of the concentration grid and the receptor concentrations. *
991	! *
992	! Author: A. Stohl *
993	! *
994	! 24 May 1995 *
995	! *
996	! 13 April 1999, Major update: if output size is smaller, dump output in *
997	! sparse matrix format; additional output of uncertainty *
998	! *
999	! 05 April 2000, Major update: output of age classes; output for backward*
1000	! runs is time spent in grid cell times total mass of *
1001	! species. *
1002	! *
1003	! 17 February 2002, Appropriate dimensions for backward and forward runs *
1004	! are now specified in module par_mod *
1005	! *
1006	! June 2006, write grid in sparse matrix with a single write command *
1007	! in order to save disk space *
1008	! *
1009	! 2008 new sparse matrix format *
1010	! *
1011	! 19 February 2010, Dominik Brunner, Empa: Adapted for COSMO *
1012	! *
1013	! April 2013, Dominik Brunner, Empa *
1014	! Adapted for netcdf output *
1015	! *
1016	!*****************************************************************************
1017	! *
1018	! Variables: *
1019	! itime current simulation time *
1020	! outnum number of samples *
1021	! *
1022	!*****************************************************************************
1023
1024	use unc_mod, only: griduncn,drygriduncn,wetgriduncn
1025
1026	implicit none
1027
1028	integer, intent(in) :: itime
1029	real, intent(in) :: outnum
1030	real :: densityoutrecept(maxreceptor)
1031	integer :: ncid,kp,ks,kz,ix,jy,iix,jjy,kzz,kzzm1,ngrid
1032	integer :: nage,i,l, jj
1033	real :: tot_mu(maxspec,maxpointspec_act)
1034	real :: halfheight,dz,dz1,dz2
1035	real :: xl,yl,xlrot,ylrot,zagnd,zagndprev
1036	real :: auxgrid(nclassunc)
1037	real, parameter :: weightair=28.97
1038
1039	! open output file
1040	call nf90_err(nf90_open(trim(ncfnamen), nf90_write, ncid))
1041
1042	! write time (do not increase time counter here, done in main output domain)
1043	call nf90_err(nf90_put_var( ncid, timeID, itime, (/ tpointer /)))
1044
1045	! For forward simulations, output fields have dimension MAXSPEC,
1046	! for backward simulations, output fields have dimension MAXPOINT.
1047	! Thus, make loops either about nspec, or about numpoint
1048	!*****************************************************************
1049
1050	if (ldirect.eq.1) then
1051	do ks=1,nspec
1052	do kp=1,maxpointspec_act
1053	tot_mu(ks,kp)=1.0
1054	end do
1055	end do
1056	else
1057	do ks=1,nspec
1058	do kp=1,maxpointspec_act
1059	tot_mu(ks,kp)=xmass(kp,ks)
1060	end do
1061	end do
1062	endif
1063
1064
1065	!*******************************************************************
1066	! Compute air density:
1067	! brd134: we now take into account whether we are in the mother or in
1068	! a nested domain (before only from mother domain)
1069	! Determine center altitude of output layer, and interpolate density
1070	! data to that altitude
1071	!*******************************************************************
1072
1073	do kz=1,numzgrid
1074	if (kz.eq.1) then
1075	halfheight=outheight(1)/2.
1076	else
1077	halfheight=(outheight(kz)+outheight(kz-1))/2.
1078	endif
1079	do kzz=2,nz
1080	if ((height(kzz-1).lt.halfheight).and. &
1081	(height(kzz).gt.halfheight)) exit
1082	end do
1083	kzz=max(min(kzz,nz),2)
1084	dz1=halfheight-height(kzz-1)
1085	dz2=height(kzz)-halfheight
1086	dz=dz1+dz2
1087
1088	do jy=0,numygridn-1
1089	do ix=0,numxgridn-1
1090	xl=outlon0n+real(ix)*dxoutn
1091	yl=outlat0n+real(jy)*dyoutn
1092	xl=(xl-xlon0)/dx
1093	yl=(yl-ylat0)/dy
1094
1095	ngrid=0
1096	do jj=numbnests,1,-1
1097	if ( xl.gt.xln(jj)+eps .and. xl.lt.xrn(jj)-eps .and. &
1098	yl.gt.yln(jj)+eps .and. yl.lt.yrn(jj)-eps ) then
1099	ngrid=jj
1100	exit
1101	end if
1102	end do
1103
1104	if (ngrid.eq.0) then
1105	iix=max(min(nint(xl),nxmin1),0)
1106	jjy=max(min(nint(yl),nymin1),0)
1107
1108	densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,memind(2))*dz1+ &
1109	rho(iix,jjy,kzz-1,memind(2))*dz2)/dz
1110	else
1111	xl=(xl-xln(ngrid))*xresoln(ngrid)
1112	yl=(yl-yln(ngrid))*yresoln(ngrid)
1113	iix=max(min(nint(xl),nxn(ngrid)-1),0)
1114	jjy=max(min(nint(yl),nyn(ngrid)-1),0)
1115	densityoutgrid(ix,jy,kz)=(rhon(iix,jjy,kzz,memind(2), ngrid)*dz1+ &
1116	rhon(iix,jjy,kzz-1,memind(2), ngrid)*dz2)/dz
1117	endif
1118
1119	end do
1120	end do
1121	end do
1122
1123	do i=1,numreceptor
1124	xl=xreceptor(i)
1125	yl=yreceptor(i)
1126	iix=max(min(nint(xl),nxmin1),0)
1127	jjy=max(min(nint(yl),nymin1),0)
1128	densityoutrecept(i)=rho(iix,jjy,1,memind(2))
1129	end do
1130
1131	! Output is different for forward and backward simulations
1132	if (ldirect.eq.1) then
1133	do kz=1,numzgrid
1134	do jy=0,numygridn-1
1135	do ix=0,numxgridn-1
1136	factor3d(ix,jy,kz)=1.e12/volumen(ix,jy,kz)/outnum
1137	end do
1138	end do
1139	end do
1140	else
1141	do kz=1,numzgrid
1142	do jy=0,numygridn-1
1143	do ix=0,numxgridn-1
1144	factor3d(ix,jy,kz)=real(abs(loutaver))/outnum
1145	end do
1146	end do
1147	end do
1148	endif
1149
1150	!*********************************************************************
1151	! Determine the standard deviation of the mean concentration or mixing
1152	! ratio (uncertainty of the output) and the dry and wet deposition
1153	!*********************************************************************
1154
1155	do ks=1,nspec
1156
1157	do kp=1,maxpointspec_act
1158	do nage=1,nageclass
1159
1160	do jy=0,numygridn-1
1161	do ix=0,numxgridn-1
1162	! WET DEPOSITION
1163	if ((WETDEP).and.(ldirect.gt.0)) then
1164	do l=1,nclassunc
1165	auxgrid(l)=wetgriduncn(ix,jy,ks,kp,l,nage)
1166	end do
1167	call mean(auxgrid,wetgrid(ix,jy), &
1168	wetgridsigma(ix,jy),nclassunc)
1169	! Multiply by number of classes to get total concentration
1170	wetgrid(ix,jy)=wetgrid(ix,jy)*real(nclassunc)
1171	! Calculate standard deviation of the mean
1172	wetgridsigma(ix,jy)= &
1173	wetgridsigma(ix,jy)* &
1174	sqrt(real(nclassunc))
1175	endif
1176
1177	! DRY DEPOSITION
1178	if ((DRYDEP).and.(ldirect.gt.0)) then
1179	do l=1,nclassunc
1180	auxgrid(l)=drygriduncn(ix,jy,ks,kp,l,nage)
1181	end do
1182	call mean(auxgrid,drygrid(ix,jy), &
1183	drygridsigma(ix,jy),nclassunc)
1184	! Multiply by number of classes to get total concentration
1185	drygrid(ix,jy)=drygrid(ix,jy)*real(nclassunc)
1186	! Calculate standard deviation of the mean
1187	drygridsigma(ix,jy)= &
1188	drygridsigma(ix,jy)* &
1189	sqrt(real(nclassunc))
1190	endif
1191
1192	! CONCENTRATION OR MIXING RATIO
1193	do kz=1,numzgrid
1194	do l=1,nclassunc
1195	auxgrid(l)=griduncn(ix,jy,kz,ks,kp,l,nage)
1196	end do
1197	call mean(auxgrid,grid(ix,jy,kz), &
1198	gridsigma(ix,jy,kz),nclassunc)
1199	! Multiply by number of classes to get total concentration
1200	grid(ix,jy,kz)= &
1201	grid(ix,jy,kz)*real(nclassunc)
1202	! Calculate standard deviation of the mean
1203	gridsigma(ix,jy,kz)= &
1204	gridsigma(ix,jy,kz)* &
1205	sqrt(real(nclassunc))
1206	end do
1207	end do
1208	end do
1209
1210
1211	!*******************************************************************
1212	! Generate output: may be in concentration (ng/m3) or in mixing
1213	! ratio (ppt) or both
1214	! Output the position and the values alternated multiplied by
1215	! 1 or -1, first line is number of values, number of positions
1216	! For backward simulations, the unit is seconds, stored in grid_time
1217	!*******************************************************************
1218
1219	! Concentration output
1220	!*********************
1221	if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
1222
1223	! Wet deposition
1224	if ((ldirect.eq.1).and.(WETDEP)) then
1225	call nf90_err(nf90_put_var(ncid,wdspecIDn(ks),1.e12*&
1226	wetgrid(0:numxgridn-1,0:numygridn-1)/arean(0:numxgridn-1,0:numygridn-1),&
1227	(/ 1,1,tpointer,kp,nage /), (/ numxgridn,numygridn,1,1,1 /)))
1228	endif
1229
1230	! Dry deposition
1231	if ((ldirect.eq.1).and.(DRYDEP)) then
1232	call nf90_err(nf90_put_var(ncid,ddspecIDn(ks),1.e12*&
1233	drygrid(0:numxgridn-1,0:numygridn-1)/arean(0:numxgridn-1,0:numygridn-1),&
1234	(/ 1,1,tpointer,kp,nage /), (/ numxgridn,numygridn,1,1,1 /)))
1235	endif
1236
1237	! Concentrations
1238	call nf90_err(nf90_put_var(ncid,specIDn(ks),grid(0:numxgridn-1,0:numygridn-1,&
1239	1:numzgrid)*factor3d(0:numxgridn-1,0:numygridn-1,1:numzgrid)/tot_mu(ks,kp),&
1240	(/ 1,1,1,tpointer,kp,nage /), (/ numxgridn,numygridn,numzgrid,1,1,1 /)))
1241
1242	endif ! concentration output
1243
1244	! Mixing ratio output
1245	!********************
1246
1247	if ((iout.eq.2)) then ! mixing ratio
1248
1249	! Wet deposition
1250	if ((ldirect.eq.1).and.(WETDEP)) then
1251	call nf90_err(nf90_put_var(ncid,wdspecIDn(ks),1.e12*&
1252	wetgrid(0:numxgridn-1,0:numygridn-1)/arean(0:numxgridn-1,0:numygridn-1),&
1253	(/ 1,1,tpointer,kp,nage /), (/ numxgridn,numygridn,1,1,1 /)))
1254	endif
1255
1256	! Dry deposition
1257	if ((ldirect.eq.1).and.(DRYDEP)) then
1258	call nf90_err(nf90_put_var(ncid,ddspecIDn(ks),1.e12*&
1259	drygrid(0:numxgridn-1,0:numygridn-1)/arean(0:numxgridn-1,0:numygridn-1),&
1260	(/ 1,1,tpointer,kp,nage /), (/ numxgridn,numygridn,1,1,1 /)))
1261	endif
1262
1263	! Mixing ratios
1264	call nf90_err(nf90_put_var(ncid,specIDnppt(ks),weightair/weightmolar(ks)*&
1265	grid(0:numxgridn-1,0:numygridn-1,1:numzgrid)*&
1266	factor3d(0:numxgridn-1,0:numygridn-1,1:numzgrid)/&
1267	densityoutgrid(0:numxgridn-1,0:numygridn-1,1:numzgrid),&
1268	(/ 1,1,1,tpointer,kp,nage /), (/ numxgridn,numygridn,numzgrid,1,1,1 /)))
1269
1270	endif ! output for ppt
1271
1272	end do
1273	end do
1274
1275	end do
1276
1277	! Close netCDF file
1278	!**************************
1279	call nf90_err(nf90_close(ncid))
1280
1281
1282	! Reinitialization of grid
1283	!*************************
1284
1285	creceptor(1:numreceptor,1:nspec) = 0.
1286	griduncn(:,:,:,1:nspec,:,:,1:nageclass) = 0.
1287
1288
1289	end subroutine concoutput_nest_ncdf
1290
1291
1292	end module nc_output_mod
1293
1294
1295	#endif

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

Download in other formats:

Original Format