source: flexpart.git/src/netcdf_output_mod.f90 @ 4c64400

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

Bugfix for double precision dry deposition calculation. Added (hardcoded) option to not use output kernel. Version/date string updated.

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