Context Navigation

netcdf_output_mod.f90 @ 11d86e7

dev

Last change on this file since 11d86e7 was 11d86e7, checked in by Sabine <sabine.eckhardt@…>, 4 years ago
bugfix in the fluxoutput
Property mode set to `100644`
File size: 64.8 KB

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

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

Download in other formats:

Original Format