Context Navigation

netcdf_output_mod.f90

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

Last change on this file was 23547f3, checked in by Ignacio Pisso <ip@…>, 4 years ago
remove tabs from files of the form src/*.f90
Property mode set to `100644`
File size: 56.3 KB

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

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

Download in other formats:

Original Format