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

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file since e9e0f06 was e9e0f06, checked in by Sabine <sabine.eckhardt@…>, 5 years ago

Removed kao and ass_spec and obsolete lines in get_wetscav.f90

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