Context Navigation

netcdf_output_mod.f90 @ 8a65cb0

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

Last change on this file since 8a65cb0 was 8a65cb0, checked in by Espen Sollum ATMOS <espen@…>, 9 years ago
Added code, makefile for dev branch
Property mode set to `100755`
File size: 54.9 KB

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

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

Download in other formats:

Original Format