source: flexpart.git/src/netcdf_output_mod.f90 @ 6a678e3

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

Added option to use double precision for calculating the deposition fields

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