source: flexpart.git/src/netcdf_output_mod.f90 @ f28aa0a

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since f28aa0a was f28aa0a, checked in by Espen Sollum ATMOS <eso@…>, 8 years ago

netCDF output; added (hardcoded) option for writing grid cell volume/area

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