Context Navigation

netcdf_output_mod.f90 @ 0ecc1fe

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

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

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

Download in other formats:

Original Format