Context Navigation

netcdf_output_mod.f90 @ ae43937

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

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

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

Download in other formats:

Original Format