Context Navigation

netcdf_output_mod.f90 @ cbad0f1

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

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

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

Download in other formats:

Original Format