source: flexpart.git/src/netcdf_output_mod.f90 @ 38b7917

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

Completed handling of nested wind fields with cloud water (for wet deposition).

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