Skip to content

GitLab

Commit 7e9777f5 authored by ronesy's avatar ronesy
Browse files

modified output format for inversion and included output of ratio of air_density to dry_air_density

parent 4c64400c
Hide whitespace changes
Inline Side-by-side
options/COMMAND
View file @ 7e9777f5
&COMMAND
LDIRECT= 1,
IBDATE= 20140907,
IBTIME= 060000,
IEDATE= 20140907,
IETIME= 120000,
LOUTSTEP= 3600,
LOUTAVER= 3600,
LOUTSAMPLE= 900,
ITSPLIT= 99999999,
LSYNCTIME= 900,
CTL= -5.0000000 ,
IFINE= 4,
IOUT= 9,
IPOUT= 0,
LSUBGRID= 1,
LCONVECTION= 1,
LAGESPECTRA= 0,
IPIN= 0,
IOUTPUTFOREACHRELEASE= 0,
IFLUX= 0,
MDOMAINFILL= 0,
IND_SOURCE= 1,
IND_RECEPTOR= 1,
MQUASILAG= 0,
NESTED_OUTPUT= 1,
LINIT_COND= 0,
SURF_ONLY= 0,
CBLFLAG= 0,
OHFIELDS_PATH= "../../flexin/",
LDIRECT= -1,
IBDATE= 20101222,
IBTIME= 000000,
IEDATE= 20110103,
IETIME= 000000,
LOUTSTEP= 86400,
LOUTAVER= 86400,
LOUTSAMPLE= 1800,
ITSPLIT= 999999999,
LSYNCTIME= 600,
CTL= -5.00,
IFINE= 4,
IOUT= 1,
IPOUT= 2,
LSUBGRID= 1,
LCONVECTION= 1,
LAGESPECTRA= 1,
IPIN= 0,
IOUTPUTFOREACHRELEASE= 1,
IFLUX= 0,
MDOMAINFILL= 0,
IND_SOURCE= 1,
IND_RECEPTOR= 2,
MQUASILAG= 0,
NESTED_OUTPUT= 0,
LINIT_COND= 2,
SURF_ONLY= 1,
LINVERSIONOUT= 1,
OHFIELDS_PATH= "/xnilu_wrk/rlt/NILU/FLEXDATA/",
/
options/SPECIES/SPECIES_029 0 → 100644
View file @ 7e9777f5
&SPECIES_PARAMS
PSPECIES="C-CO2 ",
PWEIGHTMOLAR= 12.0000000 ,
/
options/SPECIES/spec_overview
View file @ 7e9777f5
grep Tra SPEC*
SPECIES_001: PSPECIES=TRACER
SPECIES_002: PSPECIES=O3
SPECIES_003: PSPECIES=NO
SPECIES_004: PSPECIES=NO2
SPECIES_005: PSPECIES=HNO3
SPECIES_006: PSPECIES=HNO2
SPECIES_007: PSPECIES=H2O2
SPECIES_009: PSPECIES=HCHO
SPECIES_010: PSPECIES=PAN
SPECIES_011: PSPECIES=NH3
SPECIES_012: PSPECIES=SO4-aero
SPECIES_013: PSPECIES=NO3-aero
SPECIES_014: PSPECIES=I2-131
SPECIES_015: PSPECIES=I-131
SPECIES_016: PSPECIES=Cs-137
SPECIES_017: PSPECIES=Y-91
SPECIES_018: PSPECIES=Ru-106
SPECIES_019: PSPECIES=Kr-85
SPECIES_020: PSPECIES=Sr-90
SPECIES_021: PSPECIES=Xe-133
SPECIES_022: PSPECIES=CO
SPECIES_023: PSPECIES=SO2dryOH
SPECIES_024: PSPECIES=AIRTRACER
SPECIES_025: PSPECIES=AERO-TRACE
SPECIES_026: PSPECIES=CH4
SPECIES_027: PSPECIES=C2H6
SPECIES_028: PSPECIES=C3H8
SPECIES_029: PSPECIES=C-CO2
SPECIES_031: PSPECIES=PCB28
SPECIES_034: PSPECIES=G-HCH
src/advance.f90
View file @ 7e9777f5
......@@ -220,6 +220,11 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
endif
ixp=ix+1
jyp=jy+1
! fix (Espen 02.05.2017)
if(jyp>=nymax) then
write(*,*) 'WARNING: advance.f90 jyp>nymax'
jyp=jyp-1
endif
! Compute maximum mixing height around particle position
......
src/com_mod.f90
View file @ 7e9777f5
......@@ -120,6 +120,9 @@ module com_mod
integer :: lnetcdfout
! lnetcdfout 1 for netcdf grid output, 0 if not. Set in COMMAND (namelist input)
integer :: linversionout
! linversionout 1 for one grid_time output file for each release containing all timesteps
integer :: nageclass,lage(maxageclass)
! nageclass number of ageclasses for the age spectra calculation
......@@ -356,7 +359,9 @@ module com_mod
real :: clwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0 !liquid [kg/kg]
real :: ciwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0 !ice [kg/kg]
real :: clw(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0 !combined [m3/m3]
! RLT add pressure and dry air density
real :: prs(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
real :: rho_dry(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
real :: pv(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
real :: rho(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
real :: drhodz(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
......@@ -378,6 +383,7 @@ module com_mod
! uu,vv,ww [m/2] wind components in x,y and z direction
! uupol,vvpol [m/s] wind components in polar stereographic projection
! tt [K] temperature data
! prs air pressure
! qv specific humidity data
! pv (pvu) potential vorticity
! rho [kg/m3] air density
......
src/conccalc.f90
View file @ 7e9777f5
......@@ -100,6 +100,11 @@ subroutine conccalc(itime,weight)
jy=int(ytra1(i))
ixp=ix+1
jyp=jy+1
! fix (Espen 02.05.2017)
if(jyp>=nymax) then
write(*,*) 'WARNING: conccalc.f90 jyp>nymax'
jyp=jyp-1
endif
ddx=xtra1(i)-real(ix)
ddy=ytra1(i)-real(jy)
rddx=1.-ddx
......
src/concoutput.f90
View file @ 7e9777f5
......@@ -71,6 +71,9 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
integer :: sp_count_i,sp_count_r
real :: sp_fact
real :: outnum,densityoutrecept(maxreceptor),xl,yl
! RLT
real :: densitydryrecept(maxreceptor)
real :: factor_dryrecept(maxreceptor)
!real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
! +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,
......@@ -105,6 +108,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
! mind eso:added to ensure identical results between 2&3-fields versions
character(LEN=8),save :: file_stat='REPLACE'
logical :: ldates_file
logical :: lexist
integer :: ierr
character(LEN=100) :: dates_char
......@@ -200,6 +204,9 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
! rho(iix,jjy,kzz-1,2)*dz2)/dz
densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
rho(iix,jjy,kzz-1,mind)*dz2)/dz
! RLT
densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,mind)*dz1+ &
rho_dry(iix,jjy,kzz-1,mind)*dz2)/dz
end do
end do
end do
......@@ -211,8 +218,14 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
jjy=max(min(nint(yl),nymin1),0)
!densityoutrecept(i)=rho(iix,jjy,1,2)
densityoutrecept(i)=rho(iix,jjy,1,mind)
! RLT
densitydryrecept(i)=rho_dry(iix,jjy,1,mind)
end do
! RLT
! conversion factor for output relative to dry air
factor_drygrid=densityoutgrid/densitydrygrid
factor_dryrecept=densityoutrecept/densitydryrecept
! Output is different for forward and backward simulations
do kz=1,numzgrid
......@@ -351,7 +364,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
if ((ldirect.eq.1).and.(WETDEP)) then
do jy=0,numygrid-1
do ix=0,numxgrid-1
!oncentraion greater zero
!concentraion greater zero
if (wetgrid(ix,jy).gt.smallnum) then
if (sp_zer.eqv..true.) then ! first non zero value
sp_count_i=sp_count_i+1
......@@ -594,6 +607,49 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
end do
! RLT Aug 2017
! Write out conversion factor for dry air
inquire(file=path(2)(1:length(2))//'factor_drygrid',exist=lexist)
if (lexist) then
! open and append
open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
status='old',action='write',access='append')
else
! create new
open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
status='new',action='write')
endif
sp_count_i=0
sp_count_r=0
sp_fact=-1.
sp_zer=.true.
do kz=1,numzgrid
do jy=0,numygrid-1
do ix=0,numxgrid-1
if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
if (sp_zer.eqv..true.) then ! first value not equal to one
sp_count_i=sp_count_i+1
sparse_dump_i(sp_count_i)= &
ix+jy*numxgrid+kz*numxgrid*numygrid
sp_zer=.false.
sp_fact=sp_fact*(-1.)
endif
sp_count_r=sp_count_r+1
sparse_dump_r(sp_count_r)= &
sp_fact*factor_drygrid(ix,jy,kz)
else ! factor is one
sp_zer=.true.
endif
end do
end do
end do
write(unitoutfactor) sp_count_i
write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutfactor) sp_count_r
write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
close(unitoutfactor)
if (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal
if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ &
wetgridtotal
......@@ -620,7 +676,23 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
end do
endif
! RLT Aug 2017
! Write out conversion factor for dry air
if (numreceptor.gt.0) then
inquire(file=path(2)(1:length(2))//'factor_dryreceptor',exist=lexist)
if (lexist) then
! open and append
open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
status='old',action='write',access='append')
else
! create new
open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
status='new',action='write')
endif
write(unitoutfactor) itime
write(unitoutfactor) (factor_dryrecept(i),i=1,numreceptor)
close(unitoutfactor)
endif
! Reinitialization of grid
!*************************
......
src/concoutput_inversion.f90 0 → 100644
View file @ 7e9777f5
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! 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/>. *
!**********************************************************************
subroutine concoutput_inversion(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 file 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 *
!
! January 2017, Separate files by release but include all timesteps
! *
!*****************************************************************************
! *
! 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
use point_mod
use outg_mod
use par_mod
use com_mod
use mean_mod
implicit none
real(kind=dp) :: jul
integer :: itime,i,ix,jy,kz,ks,kp,l,iix,jjy,kzz,nage,jjjjmmdd,ihmmss
integer :: sp_count_i,sp_count_r
real :: sp_fact
real :: outnum,densityoutrecept(maxreceptor),xl,yl
! RLT
real :: densitydryrecept(maxreceptor)
real :: factor_dryrecept(maxreceptor)
!real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
! +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,
! + maxageclass)
!real wetgrid(0:numxgrid-1,0:numygrid-1,maxspec,maxpointspec_act,
! + maxageclass)
!real drygrid(0:numxgrid-1,0:numygrid-1,maxspec,
! + maxpointspec_act,maxageclass)
!real gridsigma(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,
! + maxpointspec_act,maxageclass),
! + drygridsigma(0:numxgrid-1,0:numygrid-1,maxspec,
! + maxpointspec_act,maxageclass),
! + wetgridsigma(0:numxgrid-1,0:numygrid-1,maxspec,
! + maxpointspec_act,maxageclass)
!real factor(0:numxgrid-1,0:numygrid-1,numzgrid)
!real sparse_dump_r(numxgrid*numygrid*numzgrid)
!integer sparse_dump_i(numxgrid*numygrid*numzgrid)
!real sparse_dump_u(numxgrid*numygrid*numzgrid)
real(dep_prec) :: auxgrid(nclassunc)
real(sp) :: gridtotal,gridsigmatotal,gridtotalunc
real(dep_prec) :: wetgridtotal,wetgridsigmatotal,wetgridtotalunc
real(dep_prec) :: drygridtotal,drygridsigmatotal,drygridtotalunc
real :: halfheight,dz,dz1,dz2,tot_mu(maxspec,maxpointspec_act)
real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
real,parameter :: weightair=28.97
logical :: sp_zer
character :: adate*8,atime*6
character(len=3) :: anspec
logical :: lexist
character :: areldate*8,areltime*6
if (verbosity.eq.1) then
print*,'inside concoutput_surf '
CALL SYSTEM_CLOCK(count_clock)
WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0
endif
! Determine current calendar date
!**********************************************************
jul=bdate+real(itime,kind=dp)/86400._dp
call caldate(jul,jjjjmmdd,ihmmss)
write(adate,'(i8.8)') jjjjmmdd
write(atime,'(i6.6)') ihmmss
!write(unitdates,'(a)') adate//atime
open(unitdates,file=path(2)(1:length(2))//'dates', ACCESS='APPEND')
write(unitdates,'(a)') adate//atime
close(unitdates)
print*, 'outnum:',outnum
print*, 'datetime:',adate//atime
! 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
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
if (verbosity.eq.1) then
print*,'concoutput_surf 2'
CALL SYSTEM_CLOCK(count_clock)
WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0
endif
!*******************************************************************
! Compute air density: sufficiently accurate to take it
! from coarse grid at some time
! 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)) goto 46
end do
46 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
xl=(xl-xlon0)/dx
yl=(yl-ylat0)/dy
iix=max(min(nint(xl),nxmin1),0)
jjy=max(min(nint(yl),nymin1),0)
densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
rho(iix,jjy,kzz-1,2)*dz2)/dz
! RLT
densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,2)*dz1+ &
rho_dry(iix,jjy,kzz-1,2)*dz2)/dz
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,2)
! RLT
densitydryrecept(i)=rho_dry(iix,jjy,1,2)
end do
! RLT
! conversion factor for output relative to dry air
factor_drygrid=densityoutgrid/densitydrygrid
factor_dryrecept=densityoutrecept/densitydryrecept
! Output is different for forward and backward simulations
do kz=1,numzgrid
do jy=0,numygrid-1
do ix=0,numxgrid-1
if (ldirect.eq.1) then
factor3d(ix,jy,kz)=1.e12/volume(ix,jy,kz)/outnum
else
factor3d(ix,jy,kz)=real(abs(loutaver))/outnum
endif
end do
end do
end do
!*********************************************************************
! Determine the standard deviation of the mean concentration or mixing
! ratio (uncertainty of the output) and the dry and wet deposition
!*********************************************************************
if (verbosity.eq.1) then
print*,'concoutput_surf 3 (sd)'
CALL SYSTEM_CLOCK(count_clock)
WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0
endif
gridtotal=0.
gridsigmatotal=0.
gridtotalunc=0.
wetgridtotal=0.
wetgridsigmatotal=0.
wetgridtotalunc=0.
drygridtotal=0.
drygridsigmatotal=0.
drygridtotalunc=0.
do ks=1,nspec
write(anspec,'(i3.3)') ks
do kp=1,maxpointspec_act
print*, 'itime = ',itime
print*, 'lage(1) = ',lage(1)
print*, 'ireleasestart(kp) = ',ireleasestart(kp)
print*, 'ireleaseend(kp) = ',ireleaseend(kp)
! check itime is within release and backward trajectory length
if (nageclass.eq.1) then
if ((itime.gt.ireleaseend(kp)).or.(itime.lt.(ireleasestart(kp)-lage(1)))) then
go to 10
endif
endif
! calculate date of release for filename
jul=bdate+real(ireleasestart(kp),kind=dp)/86400._dp ! this is the current day
call caldate(jul,jjjjmmdd,ihmmss)
write(areldate,'(i8.8)') jjjjmmdd
write(areltime,'(i6.6)') ihmmss
print*, areldate//areltime
! calculate date of field
jul=bdate+real(itime,kind=dp)/86400._dp
call caldate(jul,jjjjmmdd,ihmmss)
write(adate,'(i8.8)') jjjjmmdd
write(atime,'(i6.6)') ihmmss
print*, adate//atime
if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
if (ldirect.eq.1) then
! concentrations
inquire(file=path(2)(1:length(2))//'grid_conc_'//areldate// &
areltime//'_'//anspec,exist=lexist)
if(lexist) then
! open and append to existing file
open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//areldate// &
areltime//'_'//anspec,form='unformatted',status='old',action='write',access='append')
else
! open new file
open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//areldate// &
areltime//'_'//anspec,form='unformatted',status='new',action='write')
endif
else
! residence times
inquire(file=path(2)(1:length(2))//'grid_time_'//areldate// &
areltime//'_'//anspec,exist=lexist)
if(lexist) then
! open and append to existing file
open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//areldate// &
areltime//'_'//anspec,form='unformatted',status='old',action='write',access='append')
else
! open new file
open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//areldate// &
areltime//'_'//anspec,form='unformatted',status='new',action='write')
endif
endif
write(unitoutgrid) jjjjmmdd
write(unitoutgrid) ihmmss
endif
if ((iout.eq.2).or.(iout.eq.3)) then
! mixing ratio
inquire(file=path(2)(1:length(2))//'grid_pptv_'//areldate// &
areltime//'_'//anspec,exist=lexist)
if(lexist) then
! open and append to existing file
open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//areldate// &
areltime//'_'//anspec,form='unformatted',status='old',action='write',access='append')
else
! open new file
open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//areldate// &
areltime//'_'//anspec,form='unformatted',status='new',action='write')
endif
write(unitoutgridppt) jjjjmmdd
write(unitoutgridppt) ihmmss