Commit 1070b4c7 authored by Espen Sollum's avatar Espen Sollum
Browse files

Merge branch 'espen' into dev

parents 8c0ae9ba a8035213
...@@ -24,7 +24,7 @@ module com_mod ...@@ -24,7 +24,7 @@ module com_mod
! Variables defining where FLEXPART input/output files are stored ! Variables defining where FLEXPART input/output files are stored
!**************************************************************** !****************************************************************
character :: path(numpath+2*maxnests)*120 character :: path(numpath+2*maxnests)*200
integer :: length(numpath+2*maxnests) integer :: length(numpath+2*maxnests)
character(len=256) :: pathfile, flexversion, flexversion_major, arg1, arg2 character(len=256) :: pathfile, flexversion, flexversion_major, arg1, arg2
character(len=256) :: ohfields_path character(len=256) :: ohfields_path
......
...@@ -58,6 +58,9 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, & ...@@ -58,6 +58,9 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
integer :: sp_count_i,sp_count_r integer :: sp_count_i,sp_count_r
real :: sp_fact real :: sp_fact
real :: outnum,densityoutrecept(maxreceptor),xl,yl real :: outnum,densityoutrecept(maxreceptor),xl,yl
! RLT
real :: densitydryrecept(maxreceptor)
real :: factor_dryrecept(maxreceptor)
!real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid), !real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
! +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act, ! +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,
...@@ -92,6 +95,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, & ...@@ -92,6 +95,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
! mind eso: added to ensure identical results between 2&3-fields versions ! mind eso: added to ensure identical results between 2&3-fields versions
character(LEN=8),save :: file_stat='REPLACE' character(LEN=8),save :: file_stat='REPLACE'
logical :: ldates_file logical :: ldates_file
logical :: lexist
integer :: ierr integer :: ierr
character(LEN=100) :: dates_char character(LEN=100) :: dates_char
...@@ -193,6 +197,9 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, & ...@@ -193,6 +197,9 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
! rho(iix,jjy,kzz-1,2)*dz2)/dz ! rho(iix,jjy,kzz-1,2)*dz2)/dz
densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ & densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
rho(iix,jjy,kzz-1,mind)*dz2)/dz 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 end do
end do end do
...@@ -204,8 +211,14 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, & ...@@ -204,8 +211,14 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
jjy=max(min(nint(yl),nymin1),0) jjy=max(min(nint(yl),nymin1),0)
!densityoutrecept(i)=rho(iix,jjy,1,2) !densityoutrecept(i)=rho(iix,jjy,1,2)
densityoutrecept(i)=rho(iix,jjy,1,mind) densityoutrecept(i)=rho(iix,jjy,1,mind)
! RLT
densitydryrecept(i)=rho_dry(iix,jjy,1,mind)
end do 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 ! Output is different for forward and backward simulations
do kz=1,numzgrid do kz=1,numzgrid
...@@ -258,11 +271,9 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, & ...@@ -258,11 +271,9 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
endif endif
write(unitoutgrid) itime write(unitoutgrid) itime
endif endif
if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio
open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// & open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// &
atime//'_'//anspec,form='unformatted') atime//'_'//anspec,form='unformatted')
write(unitoutgridppt) itime write(unitoutgridppt) itime
endif endif
endif ! if deposition output endif ! if deposition output
...@@ -606,6 +617,49 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, & ...@@ -606,6 +617,49 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
end do 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 (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal
if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ & if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ &
wetgridtotal wetgridtotal
...@@ -632,7 +686,23 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, & ...@@ -632,7 +686,23 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
end do end do
endif 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 ! Reinitialization of grid
!************************* !*************************
......
...@@ -587,7 +587,6 @@ subroutine concoutput_nest(itime,outnum) ...@@ -587,7 +587,6 @@ subroutine concoutput_nest(itime,outnum)
write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r) write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
close(unitoutfactor) close(unitoutfactor)
! Reinitialization of grid ! Reinitialization of grid
!************************* !*************************
......
...@@ -39,7 +39,6 @@ subroutine concoutput_nest(itime,outnum) ...@@ -39,7 +39,6 @@ subroutine concoutput_nest(itime,outnum)
! * ! *
!***************************************************************************** !*****************************************************************************
use unc_mod use unc_mod
use point_mod use point_mod
use outg_mod use outg_mod
...@@ -55,6 +54,9 @@ subroutine concoutput_nest(itime,outnum) ...@@ -55,6 +54,9 @@ subroutine concoutput_nest(itime,outnum)
integer :: sp_count_i,sp_count_r integer :: sp_count_i,sp_count_r
real :: sp_fact real :: sp_fact
real :: outnum,densityoutrecept(maxreceptor),xl,yl real :: outnum,densityoutrecept(maxreceptor),xl,yl
! RLT
real :: densitydryrecept(maxreceptor)
real :: factor_dryrecept(maxreceptor)
!real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid), !real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
! +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act, ! +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,
...@@ -81,6 +83,7 @@ subroutine concoutput_nest(itime,outnum) ...@@ -81,6 +83,7 @@ subroutine concoutput_nest(itime,outnum)
logical :: sp_zer logical :: sp_zer
character :: adate*8,atime*6 character :: adate*8,atime*6
character(len=3) :: anspec character(len=3) :: anspec
logical :: lexist
integer :: mind integer :: mind
! mind eso:added to ensure identical results between 2&3-fields versions ! mind eso:added to ensure identical results between 2&3-fields versions
...@@ -157,6 +160,9 @@ subroutine concoutput_nest(itime,outnum) ...@@ -157,6 +160,9 @@ subroutine concoutput_nest(itime,outnum)
! rho(iix,jjy,kzz-1,2)*dz2)/dz ! rho(iix,jjy,kzz-1,2)*dz2)/dz
densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ & densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
rho(iix,jjy,kzz-1,mind)*dz2)/dz 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 end do
end do end do
...@@ -168,8 +174,14 @@ subroutine concoutput_nest(itime,outnum) ...@@ -168,8 +174,14 @@ subroutine concoutput_nest(itime,outnum)
jjy=max(min(nint(yl),nymin1),0) jjy=max(min(nint(yl),nymin1),0)
!densityoutrecept(i)=rho(iix,jjy,1,2) !densityoutrecept(i)=rho(iix,jjy,1,2)
densityoutrecept(i)=rho(iix,jjy,1,mind) densityoutrecept(i)=rho(iix,jjy,1,mind)
! RLT
densitydryrecept(i)=rho_dry(iix,jjy,1,mind)
end do 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 ! Output is different for forward and backward simulations
do kz=1,numzgrid do kz=1,numzgrid
...@@ -192,6 +204,16 @@ subroutine concoutput_nest(itime,outnum) ...@@ -192,6 +204,16 @@ subroutine concoutput_nest(itime,outnum)
do ks=1,nspec do ks=1,nspec
write(anspec,'(i3.3)') ks write(anspec,'(i3.3)') ks
if (DRYBKDEP.or.WETBKDEP) then !scavdep output
if (DRYBKDEP) &
open(unitoutgrid,file=path(2)(1:length(2))//'grid_drydep_nest_'//adate// &
atime//'_'//anspec,form='unformatted')
if (WETBKDEP) &
open(unitoutgrid,file=path(2)(1:length(2))//'grid_wetdep_nest_'//adate// &
atime//'_'//anspec,form='unformatted')
write(unitoutgrid) itime
else
if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
if (ldirect.eq.1) then if (ldirect.eq.1) then
open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_nest_' & open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_nest_' &
...@@ -204,6 +226,7 @@ subroutine concoutput_nest(itime,outnum) ...@@ -204,6 +226,7 @@ subroutine concoutput_nest(itime,outnum)
endif endif
write(unitoutgrid) itime write(unitoutgrid) itime
endif endif
endif
if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio
open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_nest_' & open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_nest_' &
...@@ -534,30 +557,72 @@ subroutine concoutput_nest(itime,outnum) ...@@ -534,30 +557,72 @@ subroutine concoutput_nest(itime,outnum)
end do end do
! RLT Aug 2017
! Write out conversion factor for dry air
! Reinitialization of grid inquire(file=path(2)(1:length(2))//'factor_drygrid_nest',exist=lexist)
!************************* if (lexist) then
! open and append
do ks=1,nspec open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
do kp=1,maxpointspec_act status='old',action='write',access='append')
do i=1,numreceptor else
creceptor(i,ks)=0. ! create new
end do open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
do jy=0,numygridn-1 status='new',action='write')
do ix=0,numxgridn-1 endif
do l=1,nclassunc sp_count_i=0
do nage=1,nageclass sp_count_r=0
do kz=1,numzgrid sp_fact=-1.
griduncn(ix,jy,kz,ks,kp,l,nage)=0. sp_zer=.true.
end do do kz=1,numzgrid
end do do jy=0,numygridn-1
end do do ix=0,numxgridn-1
end do 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*numxgridn+kz*numxgridn*numygridn
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 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)
! Reinitialization of grid
!*************************
! do ks=1,nspec
! do kp=1,maxpointspec_act
! do i=1,numreceptor
! creceptor(i,ks)=0.
! end do
! do jy=0,numygridn-1
! do ix=0,numxgridn-1
! do l=1,nclassunc
! do nage=1,nageclass
! do kz=1,numzgrid
! griduncn(ix,jy,kz,ks,kp,l,nage)=0.
! end do
! end do
! end do
! end do
! end do
! end do
! end do
creceptor(:,:)=0.
griduncn(:,:,:,:,:,:,:)=0.
if (mp_measure_time) call mpif_mtime('iotime',1) if (mp_measure_time) call mpif_mtime('iotime',1)
! if (mp_measure_time) then ! if (mp_measure_time) then
! call cpu_time(mp_root_time_end) ! call cpu_time(mp_root_time_end)
......
...@@ -632,23 +632,23 @@ subroutine concoutput_surf_nest(itime,outnum) ...@@ -632,23 +632,23 @@ subroutine concoutput_surf_nest(itime,outnum)
!************************* !*************************
do ks=1,nspec do ks=1,nspec
do kp=1,maxpointspec_act do kp=1,maxpointspec_act
do i=1,numreceptor do i=1,numreceptor
creceptor(i,ks)=0. creceptor(i,ks)=0.
end do end do
do jy=0,numygridn-1 do jy=0,numygridn-1
do ix=0,numxgridn-1 do ix=0,numxgridn-1
do l=1,nclassunc do l=1,nclassunc
do nage=1,nageclass do nage=1,nageclass
do kz=1,numzgrid do kz=1,numzgrid
griduncn(ix,jy,kz,ks,kp,l,nage)=0. griduncn(ix,jy,kz,ks,kp,l,nage)=0.
end do
end do end do
end do end do
end do end do
end do end do
end do end do
end do end do
end do
if (mp_measure_time) call mpif_mtime('iotime',1) if (mp_measure_time) call mpif_mtime('iotime',1)
! if (mp_measure_time) then ! if (mp_measure_time) then
......
...@@ -77,7 +77,7 @@ CONTAINS ...@@ -77,7 +77,7 @@ CONTAINS
! fills in attributes that can be accessed through methods in this ! fills in attributes that can be accessed through methods in this
! module ! module
USE grib_api use grib_api
IMPLICIT NONE IMPLICIT NONE
......
...@@ -54,7 +54,7 @@ subroutine gridcheck_ecmwf ...@@ -54,7 +54,7 @@ subroutine gridcheck_ecmwf
! * ! *
!********************************************************************** !**********************************************************************
use grib_api use eccodes
use par_mod use par_mod
use com_mod use com_mod
use conv_mod use conv_mod
......
...@@ -55,7 +55,7 @@ subroutine gridcheck_gfs ...@@ -55,7 +55,7 @@ subroutine gridcheck_gfs
! * ! *
!********************************************************************** !**********************************************************************
use grib_api use eccodes
use par_mod use par_mod
use com_mod use com_mod
use conv_mod use conv_mod
......
...@@ -18,7 +18,7 @@ subroutine gridcheck_nests ...@@ -18,7 +18,7 @@ subroutine gridcheck_nests
! CHANGE: 03/12/2008, Harald Sodemann, change to f90 grib_api * ! CHANGE: 03/12/2008, Harald Sodemann, change to f90 grib_api *
!***************************************************************************** !*****************************************************************************
use grib_api use eccodes
use par_mod use par_mod
use com_mod use com_mod
......
!**********************************************************************
! 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 initial_cond_output_inversion(itime)
! i
!*****************************************************************************
! *
! Output of the initial condition sensitivity field. *
! *
! 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 *
! *
!*****************************************************************************
! *
! Variables: *
! ncells number of cells with non-zero concentrations *
! sparse .true. if in sparse matrix format, else .false. *
! *
!*****************************************************************************
use unc_mod
use point_mod
use outg_mod
use par_mod
use com_mod
implicit none
integer :: itime,i,ix,jy,kz,ks,kp,sp_count_i,sp_count_r
integer :: jjjjmmdd, ihmmss
real(kind=dp) :: jul
real :: sp_fact,fact_recept
real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
logical :: sp_zer,lexist
logical,save :: listart=.true.
logical,save,allocatable,dimension(:) :: listartrel
character :: adate*8,atime*6
character :: areldate*8,areltime*6
character(len=3) :: anspec
if(listart) then
allocate(listartrel(maxpointspec_act))
listartrel(:)=.true.
endif
print*, 'listartrel = ',listartrel
!*********************************************************************
! Determine the standard deviation of the mean concentration or mixing
! ratio (uncertainty of the output) and the dry and wet deposition
!*********************************************************************
do ks=1,nspec
write(anspec,'(i3.3)') ks
do kp=1,maxpointspec_act
! calculate date of release
jul=bdate+real(ireleasestart(kp),kind=dp)/86400._dp ! this is the current day
call caldate(jul,jjjjmmdd,ihmmss)