Skip to content
Snippets Groups Projects
Commit cbd7de37 authored by Ignacio Pisso's avatar Ignacio Pisso
Browse files

from svn 

r31 | hasod | 2014-11-05 08:11:58 +0100 (Wed, 05 Nov 2014) | 2 lines
Changed paths:
   A /trunk/src/netcdf_output_mod.f90

ADD: netcdf module file
parent 37017c8f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/netcdf_output_mod.f90 0 → 100644
+ 1338
− 0
View file @ cbd7de37
!**********************************************************************
! Copyright 2013 *
! Dominik Brunner *
! *
! This file is part of FLEXPART-COSMO *
! *
! FLEXPART is free software: you can redistribute it and/or modify *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or *
! (at your option) any later version. *
! *
! FLEXPART is distributed in the hope that it will be useful, *
! but WITHOUT ANY WARRANTY; without even the implied warranty of *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
! GNU General Public License for more details. *
! *
! You should have received a copy of the GNU General Public License *
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
!*****************************************************************************
! *
! This module handles all gridded netcdf output for concentration or *
! residence time and wet and dry deposition output. *
! *
! - writeheader_netcdf generates file including all information previously *
! stored in separate header files *
! - concoutput_netcdf write concentration output and wet and dry deposition *
! *
! Author: D. Brunner *
! *
! 12 April 2013 *
! *
! HSO: 21 Oct 2014
! - added option to not writeout releases information by changing
! switch write_releases
! - additional updates for FLEXPART 9.x
!*****************************************************************************
#ifdef NETCDF_OUTPUT
module netcdf_output_mod
use netcdf
use point_mod, only: ireleasestart,ireleaseend,kindz,&
xpoint1,ypoint1,xpoint2,ypoint2,zpoint1,zpoint2,npart,xmass
use outg_mod, only: outheight,oroout,densityoutgrid,factor3d,volume,&
wetgrid,wetgridsigma,drygrid,drygridsigma,grid,gridsigma,&
area,arean,volumen, orooutn
use par_mod, only: dp, maxspec, maxreceptor, nclassunc,&
unitoutrecept,unitoutreceptppt, nxmax
use com_mod, only: path,length,ldirect,ibdate,ibtime,iedate,ietime, &
loutstep,loutaver,loutsample,outlon0,outlat0,&
numxgrid,numygrid,dxout,dyout,numzgrid, height, &
outlon0n,outlat0n,dxoutn,dyoutn,numxgridn,numygridn, &
nspec,maxpointspec_act,species,numpoint,&
dx,xlon0,dy,ylat0,compoint,method,lsubgrid,lconvection,&
ind_source,ind_receptor,nageclass,lage,&
drydep,wetdep,decay,weta,wetb, numbnests, &
weta_in,wetb_in,wetc_in,wetd_in, &
reldiff,henry,f0,density,dquer,dsigma,dryvel,&
weightmolar,ohreact,spec_ass,kao,vsetaver,&
! for concoutput_netcdf and concoutput_nest_netcdf
nxmin1,nymin1,nz,oro,oron,rho,rhon,&
memind,xresoln,yresoln,xrn, xln, yrn,yln,nxn,nyn,&
xreceptor,yreceptor,numreceptor,creceptor,iout, &
itsplit, lsynctime, ctl, ifine, lagespectra, ipin, &
ioutputforeachrelease, iflux, mdomainfill, mquasilag, &
nested_output, ipout, surf_only, linit_cond, &
flexversion
implicit none
! parameter for data compression (1-9, 9 = most aggressive)
integer, parameter :: deflate_level = 9
logical, parameter :: min_size = .false. ! if set true, redundant fields (topography) are not written to minimize file size
character(len=255), parameter :: institution = 'NILU'
integer :: tpointer
character(len=255) :: ncfname, ncfnamen
! netcdf dimension and variable IDs for main and nested output grid
integer, dimension(maxspec) :: specID,specIDppt, wdspecID,ddspecID
integer, dimension(maxspec) :: specIDn,specIDnppt, wdspecIDn,ddspecIDn
integer :: timeID, timeIDn
integer, dimension(6) :: dimids, dimidsn
integer, dimension(5) :: depdimids, depdimidsn
real,parameter :: eps=nxmax/3.e5
private:: writemetadata, output_units, nf90_err
! switch output of release point information on/off
logical, parameter :: write_releases = .true.
contains
!****************************************************************
! determine output units (see table 1 in Stohl et al., ACP 2005
!****************************************************************
subroutine output_units(units)
character(len=15), intent(out) :: units
if (ldirect.eq.1) then
! forward simulation
if (ind_source.eq.1) then
if (ind_receptor.eq.1) then
units = 'ng m-3' ! hes the kg in Tab1 is only indicating the units of the relase not the output
else
units = 'ng kg-1'
endif
else
if (ind_receptor.eq.1) then
units = 'ng m-3'
else
units = 'ng kg-1'
endif
endif
else
! backward simulation
if (ind_source.eq.1) then
if (ind_receptor.eq.1) then
units = 's'
else
units = 's m3 kg-1'
endif
else
if (ind_receptor.eq.1) then
units = 's kg m-3'
else
units = 's'
endif
endif
endif
end subroutine output_units
!****************************************************************
! write metadata to netCDF file
!****************************************************************
subroutine writemetadata(ncid,lnest)
integer, intent(in) :: ncid
logical, intent(in) :: lnest
integer :: status
character :: time*10,date*8,adate*8,atime*6
character(5) :: zone
character(255) :: login_name, host_name
! gather system information
call date_and_time(date,time,zone)
call getlog(login_name)
call hostnm(host_name)
! hes CF convention requires these attributes
call nf90_err(nf90_put_att(ncid, nf90_global, 'Conventions', 'CF-1.6'))
call nf90_err(nf90_put_att(ncid, nf90_global, 'title', 'FLEXPART model output'))
call nf90_err(nf90_put_att(ncid, nf90_global, 'institution', trim(institution)))
call nf90_err(nf90_put_att(ncid, nf90_global, 'source', trim(flexversion)//' model output'))
call nf90_err(nf90_put_att(ncid, nf90_global, 'history', date(1:4)//'-'//date(5:6)// &
'-'//date(7:8)//' '//time(1:2)//':'//time(3:4)//' '//zone//' created by '// &
trim(login_name)//' on '//trim(host_name)))
call nf90_err(nf90_put_att(ncid, nf90_global, 'references', &
'Stohl et al., Atmos. Chem. Phys., 2005, doi:10.5194/acp-5-2461-200'))
! attributes describing model run
!************************************************************************************
if (lnest) then
call nf90_err(nf90_put_att(ncid, nf90_global, 'outlon0', outlon0n))
call nf90_err(nf90_put_att(ncid, nf90_global, 'outlat0', outlat0n))
call nf90_err(nf90_put_att(ncid, nf90_global, 'dxout', dxoutn))
call nf90_err(nf90_put_att(ncid, nf90_global, 'dyout', dyoutn))
else
call nf90_err(nf90_put_att(ncid, nf90_global, 'outlon0', outlon0))
call nf90_err(nf90_put_att(ncid, nf90_global, 'outlat0', outlat0))
call nf90_err(nf90_put_att(ncid, nf90_global, 'dxout', dxout))
call nf90_err(nf90_put_att(ncid, nf90_global, 'dyout', dyout))
endif
! vertical levels stored in grid structure
! COMMAND file settings
call nf90_err(nf90_put_att(ncid, nf90_global, 'ldirect', ldirect))
write(adate,'(i8.8)') ibdate
write(atime,'(i6.6)') ibtime
call nf90_err(nf90_put_att(ncid, nf90_global, 'ibdate', adate))
call nf90_err(nf90_put_att(ncid, nf90_global, 'ibtime', atime))
write(adate,'(i8.8)') iedate
write(atime,'(i6.6)') ietime
call nf90_err(nf90_put_att(ncid, nf90_global, 'iedate', adate))
call nf90_err(nf90_put_att(ncid, nf90_global, 'ietime', atime))
call nf90_err(nf90_put_att(ncid, nf90_global, 'loutstep', loutstep))
call nf90_err(nf90_put_att(ncid, nf90_global, 'loutaver', loutaver))
call nf90_err(nf90_put_att(ncid, nf90_global, 'loutsample', loutsample))
call nf90_err(nf90_put_att(ncid, nf90_global, 'itsplit', itsplit))
call nf90_err(nf90_put_att(ncid, nf90_global, 'lsynctime', lsynctime))
call nf90_err(nf90_put_att(ncid, nf90_global, 'ctl', ctl))
call nf90_err(nf90_put_att(ncid, nf90_global, 'ifine', ifine))
call nf90_err(nf90_put_att(ncid, nf90_global, 'iout', iout))
call nf90_err(nf90_put_att(ncid, nf90_global, 'ipout', ipout))
call nf90_err(nf90_put_att(ncid, nf90_global, 'lsubgrid', lsubgrid))
call nf90_err(nf90_put_att(ncid, nf90_global, 'lconvection', lconvection))
call nf90_err(nf90_put_att(ncid, nf90_global, 'lagespectra', lagespectra))
call nf90_err(nf90_put_att(ncid, nf90_global, 'ipin', ipin))
call nf90_err(nf90_put_att(ncid, nf90_global, 'ioutputforeachrelease', ioutputforeachrelease))
call nf90_err(nf90_put_att(ncid, nf90_global, 'iflux', iflux))
call nf90_err(nf90_put_att(ncid, nf90_global, 'mdomainfill', mdomainfill))
call nf90_err(nf90_put_att(ncid, nf90_global, 'ind_source', ind_source))
call nf90_err(nf90_put_att(ncid, nf90_global, 'ind_receptor', ind_receptor))
call nf90_err(nf90_put_att(ncid, nf90_global, 'mquasilag', mquasilag))
call nf90_err(nf90_put_att(ncid, nf90_global, 'nested_output', nested_output))
call nf90_err(nf90_put_att(ncid, nf90_global, 'surf_only', surf_only))
call nf90_err(nf90_put_att(ncid, nf90_global, 'linit_cond', linit_cond))
end subroutine writemetadata
!****************************************************************
! netcdf error message handling
!****************************************************************
subroutine nf90_err(status)
integer, intent (in) :: status
if(status /= nf90_noerr) then
print *, trim(nf90_strerror(status))
stop 'Stopped'
end if
end subroutine nf90_err
!****************************************************************
! Create netcdf file and write header/metadata information
! lnest = .false. : Create main output file
! lnest = .true. : Create nested output file
!****************************************************************
subroutine writeheader_netcdf(lnest)
implicit none
logical, intent(in) :: lnest
integer :: ncid, sID, wdsID, ddsID
integer :: timeDimID, latDimID, lonDimID, levDimID
integer :: nspecDimID, npointDimID, nageclassDimID, ncharDimID, pointspecDimID
integer :: tID, lonID, latID, levID, poleID, lageID, oroID
integer :: rellng1ID, rellng2ID, rellat1ID, rellat2ID, relzz1ID, relzz2ID
integer :: relcomID, relkindzID, relstartID, relendID, relpartID, relxmassID
integer :: nnx, nny
integer, dimension(6) :: dIDs
integer, dimension(5) :: depdIDs
character(len=255) :: fname
character(len=15) :: units
character(len=10) :: fprefix
character(len=3) :: anspec
character :: adate*8,atime*6,timeunit*32
real, dimension(1000) :: coord
integer :: cache_size
integer, dimension(6) :: chunksizes
integer, dimension(5) :: dep_chunksizes
integer :: i,ix,jy
!************************
! Create netcdf file
!************************
if (ldirect.eq.1) then
write(adate,'(i8.8)') ibdate
write(atime,'(i6.6)') ibtime
fprefix = 'grid_conc_'
else
write(adate,'(i8.8)') iedate
write(atime,'(i6.6)') ietime
fprefix = 'grid_time_'
endif
if (lnest) then
fname = path(2)(1:length(2))//fprefix//adate//atime//'_nest.nc'
ncfnamen = fname
nnx = numxgridn
nny = numygridn
else
fname = path(2)(1:length(2))//fprefix//adate//atime//'.nc'
ncfname = fname
nnx = numxgrid
nny = numygrid
endif
cache_size = 16 * nnx * nny * numzgrid
! setting cache size in bytes. It is set to 4 times the largest data block that is written
! size_type x nx x ny x nz
! create file
call nf90_err(nf90_create(trim(fname), cmode = nf90_hdf5, ncid = ncid, &
cache_size = cache_size))
! create dimensions:
!*************************
! time
call nf90_err(nf90_def_dim(ncid, 'time', nf90_unlimited, timeDimID))
timeunit = 'seconds since '//adate(1:4)//'-'//adate(5:6)// &
'-'//adate(7:8)//' '//atime(1:2)//':'//atime(3:4)
! lon
call nf90_err(nf90_def_dim(ncid, 'longitude', nnx, lonDimID))
! lat
call nf90_err(nf90_def_dim(ncid, 'latitude', nny, latDimID))
! level
call nf90_err(nf90_def_dim(ncid, 'height', numzgrid, levDimID))
! number of species
call nf90_err(nf90_def_dim(ncid, 'numspec', nspec, nspecDimID))
! number of release points
call nf90_err(nf90_def_dim(ncid, 'pointspec', maxpointspec_act, pointspecDimID))
! number of age classes
call nf90_err(nf90_def_dim(ncid, 'nageclass', nageclass, nageclassDimID))
! dimension for release point characters
call nf90_err(nf90_def_dim(ncid, 'nchar', 45, ncharDimID))
! number of actual release points
call nf90_err(nf90_def_dim(ncid, 'numpoint', numpoint, npointDimID))
! create variables
!*************************
! time
call nf90_err(nf90_def_var(ncid, 'time', nf90_int, (/ timeDimID /), tID))
call nf90_err(nf90_put_att(ncid, tID, 'units', timeunit))
call nf90_err(nf90_put_att(ncid, tID, 'calendar', 'proleptic_gregorian'))
if (lnest) then
timeIDn = tID
else
timeID = tID
endif
! lon
call nf90_err(nf90_def_var(ncid, 'longitude', nf90_float, (/ lonDimID /), lonID))
call nf90_err(nf90_put_att(ncid, lonID, 'long_name', 'longitude in degree east'))
call nf90_err(nf90_put_att(ncid, lonID, 'axis', 'Lon'))
call nf90_err(nf90_put_att(ncid, lonID, 'units', 'degrees_east'))
call nf90_err(nf90_put_att(ncid, lonID, 'standard_name', 'grid_longitude'))
call nf90_err(nf90_put_att(ncid, lonID, 'description', 'grid cell centers'))
! lat
call nf90_err(nf90_def_var(ncid, 'latitude', nf90_float, (/ latDimID /), latID))
call nf90_err(nf90_put_att(ncid, latID, 'long_name', 'latitude in degree north'))
call nf90_err(nf90_put_att(ncid, latID, 'axis', 'Lat'))
call nf90_err(nf90_put_att(ncid, latID, 'units', 'degrees_north'))
call nf90_err(nf90_put_att(ncid, latID, 'standard_name', 'grid_latitude'))
call nf90_err(nf90_put_att(ncid, latID, 'description', 'grid cell centers'))
! height
call nf90_err(nf90_def_var(ncid, 'height', nf90_float, (/ levDimID /), levID))
! call nf90_err(nf90_put_att(ncid, levID, 'axis', 'Z'))
call nf90_err(nf90_put_att(ncid, levID, 'units', 'meters'))
call nf90_err(nf90_put_att(ncid, levID, 'positive', 'up'))
call nf90_err(nf90_put_att(ncid, levID, 'standard_name', 'height'))
call nf90_err(nf90_put_att(ncid, levID, 'long_name', 'height above ground'))
if (write_releases.eqv..true.) then
! release comment
call nf90_err(nf90_def_var(ncid, 'RELCOM', nf90_char, (/ ncharDimID,npointDimID /), &
relcomID))
call nf90_err(nf90_put_att(ncid, relcomID, 'long_name', 'release point name'))
! release longitude 1
call nf90_err(nf90_def_var(ncid, 'RELLNG1', nf90_float, (/ npointDimID /), rellng1ID))
call nf90_err(nf90_put_att(ncid, rellng1ID, 'units', 'degrees_east'))
call nf90_err(nf90_put_att(ncid, rellng1ID, 'long_name', &
'release longitude lower left corner'))
! release longitude 2
call nf90_err(nf90_def_var(ncid, 'RELLNG2', nf90_float, (/ npointDimID /), rellng2ID))
call nf90_err(nf90_put_att(ncid, rellng2ID, 'units', 'degrees_east'))
call nf90_err(nf90_put_att(ncid, rellng2ID, 'long_name', &
'release longitude upper right corner'))
! release latitude 1
call nf90_err(nf90_def_var(ncid, 'RELLAT1', nf90_float, (/ npointDimID /), rellat1ID))
call nf90_err(nf90_put_att(ncid, rellat1ID, 'units', 'degrees_north'))
call nf90_err(nf90_put_att(ncid, rellat1ID, 'long_name', &
'release latitude lower left corner'))
! release latitude 2
call nf90_err(nf90_def_var(ncid, 'RELLAT2', nf90_float, (/ npointDimID /), rellat2ID))
call nf90_err(nf90_put_att(ncid, rellat2ID, 'units', 'degrees_north'))
call nf90_err(nf90_put_att(ncid, rellat2ID, 'long_name', &
'release latitude upper right corner'))
! hes: if rotated_ll it would be convenient also to write the the release points in rotated_coordinates
! release height bottom
call nf90_err(nf90_def_var(ncid, 'RELZZ1', nf90_float, (/ npointDimID /), relzz1ID))
call nf90_err(nf90_put_att(ncid, relzz1ID, 'units', 'meters'))
call nf90_err(nf90_put_att(ncid, relzz1ID, 'long_name', 'release height bottom'))
! release height top
call nf90_err(nf90_def_var(ncid, 'RELZZ2', nf90_float, (/ npointDimID /), relzz2ID))
call nf90_err(nf90_put_att(ncid, relzz2ID, 'units', 'meters'))
call nf90_err(nf90_put_att(ncid, relzz2ID, 'long_name', 'release height top'))
! release kind
call nf90_err(nf90_def_var(ncid, 'RELKINDZ', nf90_int, (/ npointDimID /), relkindzID))
call nf90_err(nf90_put_att(ncid, relkindzID, 'long_name', 'release kind'))
! release start
call nf90_err(nf90_def_var(ncid, 'RELSTART', nf90_int, (/ npointDimID /), relstartID))
call nf90_err(nf90_put_att(ncid, relstartID, 'units', 'seconds'))
call nf90_err(nf90_put_att(ncid, relstartID, 'long_name', &
'release start relative to simulation start'))
! release end
call nf90_err(nf90_def_var(ncid, 'RELEND', nf90_int, (/ npointDimID /), relendID))
call nf90_err(nf90_put_att(ncid, relendID, 'units', 'seconds'))
call nf90_err(nf90_put_att(ncid, relendID, 'long_name', &
'release end relative to simulation start'))
! release particles
call nf90_err(nf90_def_var(ncid, 'RELPART', nf90_int, (/ npointDimID /), relpartID))
call nf90_err(nf90_put_att(ncid, relpartID, 'long_name', 'number of release particles'))
! release particle masses
call nf90_err(nf90_def_var(ncid, 'RELXMASS', nf90_float, (/ npointDimID, nspecDimID /), &
relxmassID))
call nf90_err(nf90_put_att(ncid, relxmassID, 'long_name', 'total release particle mass'))
end if
! age classes
call nf90_err(nf90_def_var(ncid, 'LAGE', nf90_int, (/ nageclassDimID /), lageID))
call nf90_err(nf90_put_att(ncid, lageID, 'units', 'seconds'))
call nf90_err(nf90_put_att(ncid, lageID, 'long_name', 'age class'))
! output orography
if (.not. min_size) then
call nf90_err(nf90_def_var(ncid, 'ORO', nf90_int, (/ lonDimID, latDimID /), oroID, &
deflate_level=deflate_level, chunksizes= (/ nnx, nny /)))
call nf90_err(nf90_put_att(ncid, oroID, 'standard_name', 'surface altitude'))
call nf90_err(nf90_put_att(ncid, oroID, 'long_name', 'outgrid surface altitude'))
call nf90_err(nf90_put_att(ncid, oroID, 'units', 'm'))
end if
! concentration output, wet and dry deposition variables (one per species)
call output_units(units)
dIDs = (/ londimid, latdimid, levdimid, timedimid, pointspecdimid, nageclassdimid /)
depdIDs = (/ londimid, latdimid, timedimid, pointspecdimid, nageclassdimid /)
if (lnest) then
dimidsn = dIDs
depdimidsn = depdIDs
else
dimids = dIDs
depdimids = depdIDs
endif
! set chunksizes according to largest written portion of data in an individual call to
! nf90_put_var
chunksizes = (/ nnx, nny, numzgrid, 1, 1, 1 /)
dep_chunksizes = (/ nnx, nny, 1, 1, 1 /)
do i = 1,nspec
write(anspec,'(i3.3)') i
! concentration output
if (iout.eq.1.or.iout.eq.3.or.iout.eq.5) then
call nf90_err(nf90_def_var(ncid,'spec'//anspec//'_mr', nf90_float, dIDs, sID , &
deflate_level = deflate_level, &
chunksizes = chunksizes ))
call nf90_err(nf90_put_att(ncid, sID, 'units', units))
call nf90_err(nf90_put_att(ncid, sID, 'long_name', species(i)))
call nf90_err(nf90_put_att(ncid, sID, 'decay', decay(i)))
call nf90_err(nf90_put_att(ncid, sID, 'weightmolar', weightmolar(i)))
call nf90_err(nf90_put_att(ncid, sID, 'ohreact', ohreact(i)))
call nf90_err(nf90_put_att(ncid, sID, 'kao', kao(i)))
call nf90_err(nf90_put_att(ncid, sID, 'vsetaver', vsetaver(i)))
call nf90_err(nf90_put_att(ncid, sID, 'spec_ass', spec_ass(i)))
if (lnest) then
specIDn(i) = sID
else
specID(i) = sID
endif
endif
! mixing ratio output
if (iout.eq.2.or.iout.eq.3) then
call nf90_err(nf90_def_var(ncid,'spec'//anspec//'_pptv', nf90_float, dIDs, sID , &
deflate_level = deflate_level, &
chunksizes = chunksizes ))
call nf90_err(nf90_put_att(ncid, sID, 'units', 'pptv'))
call nf90_err(nf90_put_att(ncid, sID, 'long_name', species(i)))
call nf90_err(nf90_put_att(ncid, sID, 'decay', decay(i)))
call nf90_err(nf90_put_att(ncid, sID, 'weightmolar', weightmolar(i)))
call nf90_err(nf90_put_att(ncid, sID, 'ohreact', ohreact(i)))
call nf90_err(nf90_put_att(ncid, sID, 'kao', kao(i)))
call nf90_err(nf90_put_att(ncid, sID, 'vsetaver', vsetaver(i)))
call nf90_err(nf90_put_att(ncid, sID, 'spec_ass', spec_ass(i)))
if (lnest) then
specIDnppt(i) = sID
else
specIDppt(i) = sID
endif
endif
! wet and dry deposition fields for forward runs
if (wetdep) then
call nf90_err(nf90_def_var(ncid,'WD_spec'//anspec, nf90_float, depdIDs, &
wdsID, deflate_level = deflate_level, &
chunksizes = dep_chunksizes))
call nf90_err(nf90_put_att(ncid, wdsID, 'units', '1e-12 kg m-2'))
call nf90_err(nf90_put_att(ncid, wdsID, 'weta', weta(i)))
call nf90_err(nf90_put_att(ncid, wdsID, 'wetb', wetb(i)))
call nf90_err(nf90_put_att(ncid, wdsID, 'weta_in', weta_in(i)))
call nf90_err(nf90_put_att(ncid, wdsID, 'wetb_in', wetb_in(i)))
call nf90_err(nf90_put_att(ncid, wdsID, 'wetc_in', wetc_in(i)))
call nf90_err(nf90_put_att(ncid, wdsID, 'wetd_in', wetd_in(i)))
call nf90_err(nf90_put_att(ncid, wdsID, 'dquer', dquer(i)))
call nf90_err(nf90_put_att(ncid, wdsID, 'henry', henry(i)))
if (lnest) then
wdspecIDn(i) = wdsID
else
wdspecID(i) = wdsID
endif
endif
if (drydep) then
call nf90_err(nf90_def_var(ncid,'DD_spec'//anspec, nf90_float, depdIDs, &
ddsID, deflate_level = deflate_level, &
chunksizes = dep_chunksizes))
call nf90_err(nf90_put_att(ncid, ddsID, 'units', '1e-12 kg m-2'))
call nf90_err(nf90_put_att(ncid, ddsID, 'dryvel', dryvel(i)))
call nf90_err(nf90_put_att(ncid, ddsID, 'reldiff', reldiff(i)))
call nf90_err(nf90_put_att(ncid, ddsID, 'henry', henry(i)))
call nf90_err(nf90_put_att(ncid, ddsID, 'f0', f0(i)))
call nf90_err(nf90_put_att(ncid, ddsID, 'dquer', dquer(i)))
call nf90_err(nf90_put_att(ncid, ddsID, 'density', density(i)))
call nf90_err(nf90_put_att(ncid, ddsID, 'dsigma', dsigma(i)))
if (lnest) then
ddspecIDn(i) = ddsID
else
ddspecID(i) = ddsID
endif
endif
end do
! global (metadata) attributes
!*******************************
call writemetadata(ncid,lnest)
! moves the file from define to data mode
call nf90_err(nf90_enddef(ncid))
! ! hes: inquire var definition
! do i = 1,nspec
! write(anspec,'(i3.3)') i
!
! ! concentration output
! if (iout.eq.1.or.iout.eq.3.or.iout.eq.5) then
! if (lnest) then
! sID = specIDn(i)
! else
! sID = specID(i)
! endif
! call nf90_err(nf90_inquire_variable(ncid, sID, chunksizes=inq_chunksizes))
! write(*,*) "Chunksizes for var "//anspec//": ", inq_chunksizes
! endif
! end do
! fill with data
!******************************
! longitudes (grid cell centers)
if (lnest) then
do i = 1,numxgridn
coord(i) = outlon0n + (i-0.5)*dxoutn
enddo
call nf90_err(nf90_put_var(ncid, lonID, coord(1:numxgridn)))
else
do i = 1,numxgrid
coord(i) = outlon0 + (i-0.5)*dxout
enddo
call nf90_err(nf90_put_var(ncid, lonID, coord(1:numxgrid)))
endif
! latitudes (grid cell centers)
if (lnest) then
do i = 1,numygridn
coord(i) = outlat0n + (i-0.5)*dyoutn
enddo
call nf90_err(nf90_put_var(ncid, latID, coord(1:numygridn)))
else
do i = 1,numygrid
coord(i) = outlat0 + (i-0.5)*dyout
enddo
call nf90_err(nf90_put_var(ncid, latID, coord(1:numygrid)))
endif
! levels
call nf90_err(nf90_put_var(ncid, levID, outheight(1:numzgrid)))
if (write_releases.eqv..true.) then
! release point information
do i = 1,numpoint
call nf90_err(nf90_put_var(ncid, relstartID, ireleasestart(i), (/i/)))
call nf90_err(nf90_put_var(ncid, relendID, ireleaseend(i), (/i/)))
call nf90_err(nf90_put_var(ncid, relkindzID, kindz(i), (/i/)))
call nf90_err(nf90_put_var(ncid, rellng1ID, xpoint1(i), (/i/)))
call nf90_err(nf90_put_var(ncid, rellng2ID, xpoint2(i), (/i/)))
call nf90_err(nf90_put_var(ncid, rellat1ID, ypoint1(i), (/i/)))
call nf90_err(nf90_put_var(ncid, rellat2ID, ypoint2(i), (/i/)))
call nf90_err(nf90_put_var(ncid, relzz1ID, zpoint1(i), (/i/)))
call nf90_err(nf90_put_var(ncid, relzz2ID, zpoint2(i), (/i/)))
call nf90_err(nf90_put_var(ncid, relpartID, npart(i), (/i/)))
if (i .le. 1000) then
call nf90_err(nf90_put_var(ncid, relcomID, compoint(i), (/1,i/), (/45,1/)))
else
call nf90_err(nf90_put_var(ncid, relcomID, 'NA', (/1,i/), (/45,1/)))
endif
call nf90_err(nf90_put_var(ncid, relxmassID, xmass(i,1:nspec), (/i,1/), (/1,nspec/)))
end do
end if
! age classes
call nf90_err(nf90_put_var(ncid, lageID, lage(1:nageclass)))
! orography
if (.not. min_size) then
if (lnest) then
call nf90_err(nf90_put_var(ncid, oroID, orooutn(0:(nnx-1), 0:(nny-1))))
else
call nf90_err(nf90_put_var(ncid, oroID, oroout(0:(nnx-1), 0:(nny-1))))
endif
end if
call nf90_err(nf90_close(ncid))
return
end subroutine writeheader_netcdf
subroutine concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
! i i o o
! o
!*****************************************************************************
! *
! Output of the concentration grid and the receptor concentrations. *
! *
! Author: A. Stohl *
! *
! 24 May 1995 *
! *
! 13 April 1999, Major update: if output size is smaller, dump output in *
! sparse matrix format; additional output of uncertainty *
! *
! 05 April 2000, Major update: output of age classes; output for backward*
! runs is time spent in grid cell times total mass of *
! species. *
! *
! 17 February 2002, Appropriate dimensions for backward and forward runs *
! are now specified in module par_mod *
! *
! June 2006, write grid in sparse matrix with a single write command *
! in order to save disk space *
! *
! 2008 new sparse matrix format *
! *
! February 2010, Dominik Brunner, Empa *
! Adapted for COSMO *
! Remark: calculation of density could be improved. *
! Currently, it is calculated for the lower left corner *
! of each output grid cell rather than for its center. *
! Furthermore, the average density could be calculated *
! from the difference in pressure at the top and bottom *
! of each cell rather than by interpolation. *
! *
! April 2013, Dominik Brunner, Empa *
! Adapted for netcdf output *
! *
!*****************************************************************************
! *
! Variables: *
! outnum number of samples *
! ncells number of cells with non-zero concentrations *
! sparse .true. if in sparse matrix format, else .false. *
! tot_mu 1 for forward, initial mass mixing ration for backw. runs *
! *
!*****************************************************************************
use unc_mod, only: gridunc,drygridunc,wetgridunc
implicit none
integer, intent(in) :: itime
real, intent(in) :: outnum
real, intent(out) :: gridtotalunc,wetgridtotalunc,drygridtotalunc
real :: densityoutrecept(maxreceptor)
integer :: ncid,kp,ks,kz,ix,jy,iix,jjy,kzz,kzzm1,ngrid
integer :: nage,i,l, jj
real :: tot_mu(maxspec,maxpointspec_act)
real :: halfheight,dz,dz1,dz2
real :: xl,yl,xlrot,ylrot,zagnd,zagndprev
real :: auxgrid(nclassunc),gridtotal,gridsigmatotal
real :: wetgridtotal,wetgridsigmatotal
real :: drygridtotal,drygridsigmatotal
real, parameter :: weightair=28.97
! open output file
call nf90_err(nf90_open(trim(ncfname), nf90_write, ncid))
! write time
tpointer = tpointer + 1
call nf90_err(nf90_put_var( ncid, timeID, itime, (/ tpointer /)))
! For forward simulations, output fields have dimension MAXSPEC,
! for backward simulations, output fields have dimension MAXPOINT.
! Thus, make loops either about nspec, or about numpoint
!*****************************************************************
if (ldirect.eq.1) then
do ks=1,nspec
do kp=1,maxpointspec_act
tot_mu(ks,kp)=1.0
end do
end do
else
do ks=1,nspec
do kp=1,maxpointspec_act
tot_mu(ks,kp)=xmass(kp,ks)
end do
end do
endif
!*******************************************************************
! Compute air density:
! brd134: we now take into account whether we are in the mother or in
! a nested domain (before only from mother domain)
! Determine center altitude of output layer, and interpolate density
! data to that altitude
!*******************************************************************
do kz=1,numzgrid
if (kz.eq.1) then
halfheight=outheight(1)/2.
else
halfheight=(outheight(kz)+outheight(kz-1))/2.
endif
do kzz=2,nz
if ((height(kzz-1).lt.halfheight).and. &
(height(kzz).gt.halfheight)) exit
end do
kzz=max(min(kzz,nz),2)
dz1=halfheight-height(kzz-1)
dz2=height(kzz)-halfheight
dz=dz1+dz2
do jy=0,numygrid-1
do ix=0,numxgrid-1
xl=outlon0+real(ix)*dxout
yl=outlat0+real(jy)*dyout
! grid index in mother domain
xl=(xl-xlon0)/dx
yl=(yl-ylat0)/dx
ngrid=0
do jj=numbnests,1,-1
if ( xl.gt.xln(jj)+eps .and. xl.lt.xrn(jj)-eps .and. &
yl.gt.yln(jj)+eps .and. yl.lt.yrn(jj)-eps ) then
ngrid=jj
exit
end if
end do
if (ngrid.eq.0) then
iix=max(min(nint(xl),nxmin1),0) ! if output grid cell is outside mother domain
jjy=max(min(nint(yl),nymin1),0)
densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,memind(2))*dz1+ &
rho(iix,jjy,kzz-1,memind(2))*dz2)/dz
else
xl=(xl-xln(ngrid))*xresoln(ngrid)
yl=(yl-yln(ngrid))*yresoln(ngrid)
iix=max(min(nint(xl),nxn(ngrid)-1),0)
jjy=max(min(nint(yl),nyn(ngrid)-1),0)
densityoutgrid(ix,jy,kz)=(rhon(iix,jjy,kzz,memind(2), ngrid)*dz1+ &
rhon(iix,jjy,kzz-1,memind(2), ngrid)*dz2)/dz
endif
end do
end do
end do
! brd134: for receptor points no option for nests yet to specify density
! and also altitude zreceptor not considered yet (needs revision)
do i=1,numreceptor
xl=xreceptor(i)
yl=yreceptor(i)
iix=max(min(nint(xl),nxmin1),0)
jjy=max(min(nint(yl),nymin1),0)
densityoutrecept(i)=rho(iix,jjy,1,memind(2))
end do
! Output is different for forward and backward simulations
if (ldirect.eq.1) then
do kz=1,numzgrid
do jy=0,numygrid-1
do ix=0,numxgrid-1
factor3d(ix,jy,kz)=1.e12/volume(ix,jy,kz)/outnum
end do
end do
end do
else
do kz=1,numzgrid
do jy=0,numygrid-1
do ix=0,numxgrid-1
factor3d(ix,jy,kz)=real(abs(loutaver))/outnum
end do
end do
end do
endif
!*********************************************************************
! Determine the standard deviation of the mean concentration or mixing
! ratio (uncertainty of the output) and the dry and wet deposition
!*********************************************************************
gridtotal=0.
gridsigmatotal=0.
gridtotalunc=0.
wetgridtotal=0.
wetgridsigmatotal=0.
wetgridtotalunc=0.
drygridtotal=0.
drygridsigmatotal=0.
drygridtotalunc=0.
do ks=1,nspec
do kp=1,maxpointspec_act
do nage=1,nageclass
do jy=0,numygrid-1
do ix=0,numxgrid-1
! WET DEPOSITION
if ((wetdep).and.(ldirect.gt.0)) then
do l=1,nclassunc
auxgrid(l)=wetgridunc(ix,jy,ks,kp,l,nage)
end do
call mean(auxgrid,wetgrid(ix,jy), &
wetgridsigma(ix,jy),nclassunc)
! Multiply by number of classes to get total concentration
wetgrid(ix,jy)=wetgrid(ix,jy)*real(nclassunc)
wetgridtotal=wetgridtotal+wetgrid(ix,jy)
! Calculate standard deviation of the mean
wetgridsigma(ix,jy)= &
wetgridsigma(ix,jy)* &
sqrt(real(nclassunc))
wetgridsigmatotal=wetgridsigmatotal+ &
wetgridsigma(ix,jy)
endif
! DRY DEPOSITION
if ((drydep).and.(ldirect.gt.0)) then
do l=1,nclassunc
auxgrid(l)=drygridunc(ix,jy,ks,kp,l,nage)
end do
call mean(auxgrid,drygrid(ix,jy), &
drygridsigma(ix,jy),nclassunc)
! Multiply by number of classes to get total concentration
drygrid(ix,jy)=drygrid(ix,jy)*real(nclassunc)
drygridtotal=drygridtotal+drygrid(ix,jy)
! Calculate standard deviation of the mean
drygridsigma(ix,jy)= &
drygridsigma(ix,jy)* &
sqrt(real(nclassunc))
drygridsigmatotal=drygridsigmatotal+ &
drygridsigma(ix,jy)
endif
! CONCENTRATION OR MIXING RATIO
do kz=1,numzgrid
do l=1,nclassunc
auxgrid(l)=gridunc(ix,jy,kz,ks,kp,l,nage)
end do
call mean(auxgrid,grid(ix,jy,kz), &
gridsigma(ix,jy,kz),nclassunc)
! Multiply by number of classes to get total concentration
grid(ix,jy,kz)= &
grid(ix,jy,kz)*real(nclassunc)
gridtotal=gridtotal+grid(ix,jy,kz)
! Calculate standard deviation of the mean
gridsigma(ix,jy,kz)= &
gridsigma(ix,jy,kz)* &
sqrt(real(nclassunc))
gridsigmatotal=gridsigmatotal+ &
gridsigma(ix,jy,kz)
end do
end do
end do
! print*,gridtotal,maxpointspec_act
!*******************************************************************
! Generate output: may be in concentration (ng/m3) or in mixing
! ratio (ppt) or both
! Output the position and the values alternated multiplied by
! 1 or -1, first line is number of values, number of positions
! For backward simulations, the unit is seconds, stored in grid_time
!*******************************************************************
! Concentration output
!*********************
if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
! Wet deposition
if ((ldirect.eq.1).and.(WETDEP)) then
call nf90_err(nf90_put_var(ncid,wdspecID(ks),1.e12*&
wetgrid(0:numxgrid-1,0:numygrid-1)/area(0:numxgrid-1,0:numygrid-1),&
(/ 1,1,tpointer,kp,nage /), (/ numxgrid,numygrid,1,1,1 /)))
endif
! Dry deposition
if ((ldirect.eq.1).and.(DRYDEP)) then
call nf90_err(nf90_put_var(ncid,ddspecID(ks),1.e12*&
drygrid(0:numxgrid-1,0:numygrid-1)/area(0:numxgrid-1,0:numygrid-1),&
(/ 1,1,tpointer,kp,nage /), (/ numxgrid,numygrid,1,1,1 /)))
endif
! Concentrations
call nf90_err(nf90_put_var(ncid,specID(ks),grid(0:numxgrid-1,0:numygrid-1,&
1:numzgrid)*factor3d(0:numxgrid-1,0:numygrid-1,1:numzgrid)/tot_mu(ks,kp),&
(/ 1,1,1,tpointer,kp,nage /), (/ numxgrid,numygrid,numzgrid,1,1,1 /) ))
endif ! concentration output
! Mixing ratio output
!********************
if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio
! Wet deposition
if ((ldirect.eq.1).and.(WETDEP)) then
call nf90_err(nf90_put_var(ncid,wdspecID(ks),1.e12*&
wetgrid(0:numxgrid-1,0:numygrid-1)/area(0:numxgrid-1,0:numygrid-1),&
(/ 1,1,tpointer,kp,nage /), (/ numxgrid,numygrid,1,1,1 /)))
endif
! Dry deposition
if ((ldirect.eq.1).and.(DRYDEP)) then
call nf90_err(nf90_put_var(ncid,ddspecID(ks),1.e12*&
drygrid(0:numxgrid-1,0:numygrid-1)/area(0:numxgrid-1,0:numygrid-1),&
(/ 1,1,tpointer,kp,nage /), (/ numxgrid,numygrid,1,1,1 /)))
endif
! Mixing ratios
call nf90_err(nf90_put_var(ncid,specIDppt(ks),weightair/weightmolar(ks)*&
grid(0:numxgrid-1,0:numygrid-1,1:numzgrid)*&
factor3d(0:numxgrid-1,0:numygrid-1,1:numzgrid)/&
densityoutgrid(0:numxgrid-1,0:numygrid-1,1:numzgrid),&
(/ 1,1,1,tpointer,kp,nage /), (/ numxgrid,numygrid,numzgrid,1,1,1 /)))
endif ! output for ppt
end do
end do
end do
! Close netCDF file
!**************************
call nf90_err(nf90_close(ncid))
if (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal
if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ &
wetgridtotal
if (drygridtotal.gt.0.) drygridtotalunc=drygridsigmatotal/ &
drygridtotal
! Dump of receptor concentrations
if (numreceptor.gt.0 .and. (iout.eq.2 .or. iout.eq.3) ) then
write(unitoutreceptppt) itime
do ks=1,nspec
write(unitoutreceptppt) (1.e12*creceptor(i,ks)/outnum* &
weightair/weightmolar(ks)/densityoutrecept(i),i=1,numreceptor)
end do
endif
! Dump of receptor concentrations
if (numreceptor.gt.0) then
write(unitoutrecept) itime
do ks=1,nspec
write(unitoutrecept) (1.e12*creceptor(i,ks)/outnum,i=1,numreceptor)
end do
endif
! Reinitialization of grid
!*************************
creceptor(1:numreceptor,1:nspec) = 0.
gridunc(:,:,:,1:nspec,:,:,1:nageclass) = 0.
end subroutine concoutput_netcdf
subroutine concoutput_surf_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
use unc_mod, only: gridunc,drygridunc,wetgridunc
implicit none
integer, intent(in) :: itime
real, intent(in) :: outnum
real, intent(out) :: gridtotalunc,wetgridtotalunc,drygridtotalunc
print*,'Netcdf output for surface only not yet implemented'
end subroutine concoutput_surf_netcdf
subroutine concoutput_nest_netcdf(itime,outnum)
! i i
!*****************************************************************************
! *
! Output of the concentration grid and the receptor concentrations. *
! *
! Author: A. Stohl *
! *
! 24 May 1995 *
! *
! 13 April 1999, Major update: if output size is smaller, dump output in *
! sparse matrix format; additional output of uncertainty *
! *
! 05 April 2000, Major update: output of age classes; output for backward*
! runs is time spent in grid cell times total mass of *
! species. *
! *
! 17 February 2002, Appropriate dimensions for backward and forward runs *
! are now specified in module par_mod *
! *
! June 2006, write grid in sparse matrix with a single write command *
! in order to save disk space *
! *
! 2008 new sparse matrix format *
! *
! 19 February 2010, Dominik Brunner, Empa: Adapted for COSMO *
! *
! April 2013, Dominik Brunner, Empa *
! Adapted for netcdf output *
! *
!*****************************************************************************
! *
! Variables: *
! itime current simulation time *
! outnum number of samples *
! *
!*****************************************************************************
use unc_mod, only: griduncn,drygriduncn,wetgriduncn
implicit none
integer, intent(in) :: itime
real, intent(in) :: outnum
real :: densityoutrecept(maxreceptor)
integer :: ncid,kp,ks,kz,ix,jy,iix,jjy,kzz,kzzm1,ngrid
integer :: nage,i,l, jj
real :: tot_mu(maxspec,maxpointspec_act)
real :: halfheight,dz,dz1,dz2
real :: xl,yl,xlrot,ylrot,zagnd,zagndprev
real :: auxgrid(nclassunc),gridtotal
real, parameter :: weightair=28.97
! open output file
call nf90_err(nf90_open(trim(ncfnamen), nf90_write, ncid))
! write time (do not increase time counter here, done in main output domain)
call nf90_err(nf90_put_var( ncid, timeID, itime, (/ tpointer /)))
! For forward simulations, output fields have dimension MAXSPEC,
! for backward simulations, output fields have dimension MAXPOINT.
! Thus, make loops either about nspec, or about numpoint
!*****************************************************************
if (ldirect.eq.1) then
do ks=1,nspec
do kp=1,maxpointspec_act
tot_mu(ks,kp)=1.0
end do
end do
else
do ks=1,nspec
do kp=1,maxpointspec_act
tot_mu(ks,kp)=xmass(kp,ks)
end do
end do
endif
!*******************************************************************
! Compute air density:
! brd134: we now take into account whether we are in the mother or in
! a nested domain (before only from mother domain)
! Determine center altitude of output layer, and interpolate density
! data to that altitude
!*******************************************************************
do kz=1,numzgrid
if (kz.eq.1) then
halfheight=outheight(1)/2.
else
halfheight=(outheight(kz)+outheight(kz-1))/2.
endif
do kzz=2,nz
if ((height(kzz-1).lt.halfheight).and. &
(height(kzz).gt.halfheight)) exit
end do
kzz=max(min(kzz,nz),2)
dz1=halfheight-height(kzz-1)
dz2=height(kzz)-halfheight
dz=dz1+dz2
do jy=0,numygridn-1
do ix=0,numxgridn-1
xl=outlon0n+real(ix)*dxoutn
yl=outlat0n+real(jy)*dyoutn
xl=(xl-xlon0)/dx
yl=(yl-ylat0)/dy
ngrid=0
do jj=numbnests,1,-1
if ( xl.gt.xln(jj)+eps .and. xl.lt.xrn(jj)-eps .and. &
yl.gt.yln(jj)+eps .and. yl.lt.yrn(jj)-eps ) then
ngrid=jj
exit
end if
end do
if (ngrid.eq.0) then
iix=max(min(nint(xl),nxmin1),0)
jjy=max(min(nint(yl),nymin1),0)
densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,memind(2))*dz1+ &
rho(iix,jjy,kzz-1,memind(2))*dz2)/dz
else
xl=(xl-xln(ngrid))*xresoln(ngrid)
yl=(yl-yln(ngrid))*yresoln(ngrid)
iix=max(min(nint(xl),nxn(ngrid)-1),0)
jjy=max(min(nint(yl),nyn(ngrid)-1),0)
densityoutgrid(ix,jy,kz)=(rhon(iix,jjy,kzz,memind(2), ngrid)*dz1+ &
rhon(iix,jjy,kzz-1,memind(2), ngrid)*dz2)/dz
endif
end do
end do
end do
do i=1,numreceptor
xl=xreceptor(i)
yl=yreceptor(i)
iix=max(min(nint(xl),nxmin1),0)
jjy=max(min(nint(yl),nymin1),0)
densityoutrecept(i)=rho(iix,jjy,1,memind(2))
end do
! Output is different for forward and backward simulations
if (ldirect.eq.1) then
do kz=1,numzgrid
do jy=0,numygridn-1
do ix=0,numxgridn-1
factor3d(ix,jy,kz)=1.e12/volumen(ix,jy,kz)/outnum
end do
end do
end do
else
do kz=1,numzgrid
do jy=0,numygridn-1
do ix=0,numxgridn-1
factor3d(ix,jy,kz)=real(abs(loutaver))/outnum
end do
end do
end do
endif
!*********************************************************************
! Determine the standard deviation of the mean concentration or mixing
! ratio (uncertainty of the output) and the dry and wet deposition
!*********************************************************************
gridtotal=0.
do ks=1,nspec
do kp=1,maxpointspec_act
do nage=1,nageclass
do jy=0,numygridn-1
do ix=0,numxgridn-1
! WET DEPOSITION
if ((WETDEP).and.(ldirect.gt.0)) then
do l=1,nclassunc
auxgrid(l)=wetgriduncn(ix,jy,ks,kp,l,nage)
end do
call mean(auxgrid,wetgrid(ix,jy), &
wetgridsigma(ix,jy),nclassunc)
! Multiply by number of classes to get total concentration
wetgrid(ix,jy)=wetgrid(ix,jy)*real(nclassunc)
! Calculate standard deviation of the mean
wetgridsigma(ix,jy)= &
wetgridsigma(ix,jy)* &
sqrt(real(nclassunc))
endif
! DRY DEPOSITION
if ((DRYDEP).and.(ldirect.gt.0)) then
do l=1,nclassunc
auxgrid(l)=drygriduncn(ix,jy,ks,kp,l,nage)
end do
call mean(auxgrid,drygrid(ix,jy), &
drygridsigma(ix,jy),nclassunc)
! Multiply by number of classes to get total concentration
drygrid(ix,jy)=drygrid(ix,jy)*real(nclassunc)
! Calculate standard deviation of the mean
drygridsigma(ix,jy)= &
drygridsigma(ix,jy)* &
sqrt(real(nclassunc))
endif
! CONCENTRATION OR MIXING RATIO
do kz=1,numzgrid
do l=1,nclassunc
auxgrid(l)=griduncn(ix,jy,kz,ks,kp,l,nage)
end do
call mean(auxgrid,grid(ix,jy,kz), &
gridsigma(ix,jy,kz),nclassunc)
! Multiply by number of classes to get total concentration
grid(ix,jy,kz)= &
grid(ix,jy,kz)*real(nclassunc)
gridtotal=gridtotal+grid(ix,jy,kz)
! Calculate standard deviation of the mean
gridsigma(ix,jy,kz)= &
gridsigma(ix,jy,kz)* &
sqrt(real(nclassunc))
end do
end do
end do
! print*,gridtotal,maxpointspec_act
!*******************************************************************
! Generate output: may be in concentration (ng/m3) or in mixing
! ratio (ppt) or both
! Output the position and the values alternated multiplied by
! 1 or -1, first line is number of values, number of positions
! For backward simulations, the unit is seconds, stored in grid_time
!*******************************************************************
! Concentration output
!*********************
if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
! Wet deposition
if ((ldirect.eq.1).and.(WETDEP)) then
call nf90_err(nf90_put_var(ncid,wdspecIDn(ks),1.e12*&
wetgrid(0:numxgridn-1,0:numygridn-1)/arean(0:numxgridn-1,0:numygridn-1),&
(/ 1,1,tpointer,kp,nage /), (/ numxgridn,numygridn,1,1,1 /)))
endif
! Dry deposition
if ((ldirect.eq.1).and.(DRYDEP)) then
call nf90_err(nf90_put_var(ncid,ddspecIDn(ks),1.e12*&
drygrid(0:numxgridn-1,0:numygridn-1)/arean(0:numxgridn-1,0:numygridn-1),&
(/ 1,1,tpointer,kp,nage /), (/ numxgridn,numygridn,1,1,1 /)))
endif
! Concentrations
call nf90_err(nf90_put_var(ncid,specIDn(ks),grid(0:numxgridn-1,0:numygridn-1,&
1:numzgrid)*factor3d(0:numxgridn-1,0:numygridn-1,1:numzgrid)/tot_mu(ks,kp),&
(/ 1,1,1,tpointer,kp,nage /), (/ numxgridn,numygridn,numzgrid,1,1,1 /)))
endif ! concentration output
! Mixing ratio output
!********************
if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio
! Wet deposition
if ((ldirect.eq.1).and.(WETDEP)) then
call nf90_err(nf90_put_var(ncid,wdspecIDn(ks),1.e12*&
wetgrid(0:numxgridn-1,0:numygridn-1)/arean(0:numxgridn-1,0:numygridn-1),&
(/ 1,1,tpointer,kp,nage /), (/ numxgridn,numygridn,1,1,1 /)))
endif
! Dry deposition
if ((ldirect.eq.1).and.(DRYDEP)) then
call nf90_err(nf90_put_var(ncid,ddspecIDn(ks),1.e12*&
drygrid(0:numxgridn-1,0:numygridn-1)/arean(0:numxgridn-1,0:numygridn-1),&
(/ 1,1,tpointer,kp,nage /), (/ numxgridn,numygridn,1,1,1 /)))
endif
! Mixing ratios
call nf90_err(nf90_put_var(ncid,specIDnppt(ks),weightair/weightmolar(ks)*&
grid(0:numxgridn-1,0:numygridn-1,1:numzgrid)*&
factor3d(0:numxgridn-1,0:numygridn-1,1:numzgrid)/&
densityoutgrid(0:numxgridn-1,0:numygridn-1,1:numzgrid),&
(/ 1,1,1,tpointer,kp,nage /), (/ numxgridn,numygridn,numzgrid,1,1,1 /)))
endif ! output for ppt
end do
end do
end do
! Close netCDF file
!**************************
call nf90_err(nf90_close(ncid))
! Reinitialization of grid
!*************************
creceptor(1:numreceptor,1:nspec) = 0.
griduncn(:,:,:,1:nspec,:,:,1:nageclass) = 0.
end subroutine concoutput_nest_netcdf
subroutine concoutput_surf_nest_netcdf(itime,outnum)
implicit none
integer, intent(in) :: itime
real, intent(in) :: outnum
print*,'Netcdf output for surface only not yet implemented'
end subroutine concoutput_surf_nest_netcdf
end module netcdf_output_mod
#endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment