source: branches/flexpart91_hasod/src_parallel/netcdf_output_mod.f90 @ 10

Last change on this file since 10 was 10, checked in by hasod, 11 years ago

ADD: namelist input implemented for all common input files

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