Context Navigation

netcdf_output_mod.f90 @ d8107c2

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

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

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

Download in other formats:

Original Format