source: trunk/src/netcdf_output_mod.f90

Last change on this file was 31, checked in by hasod, 9 years ago

ADD: netcdf module file

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