Context Navigation

netcdf_output_mod.f90 @ e2132cf

dev

Last change on this file since e2132cf was e2132cf, checked in by Ignacio Pisso <ip@…>, 3 years ago
fix permissions
Property mode set to `100644`
File size: 68.6 KB

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

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

Download in other formats:

Original Format