Context Navigation

netcdf_output_mod.f90 @ 4ad96c5

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

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

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

Download in other formats:

Original Format