Context Navigation

netcdf_output_mod.f90 @ b423ff6

NetCDF

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

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

Download in other formats:

Original Format