Commit 4c64400c authored by Espen Sollum's avatar Espen Sollum
Browse files

Bugfix for double precision dry deposition calculation. Added (hardcoded)...

Bugfix for double precision dry deposition calculation. Added (hardcoded) option to not use output kernel. Version/date string updated.
parent 16b61a5e
......@@ -66,9 +66,10 @@ program flexpart
end do
call gasdev1(idummy,rannumb(maxrand),rannumb(maxrand-1))
! FLEXPART version string
flexversion_major = '10' ! Major version number, also used for species file names
flexversion='Version '//trim(flexversion_major)//'.0beta (2015-05-01)'
flexversion='Version '//trim(flexversion_major)//'.1beta (2016-11-02)'
verbosity=0
! Read the pathnames where input/output files are stored
......@@ -383,6 +384,17 @@ program flexpart
end do
end do
! Inform whether output kernel is used or not
!*********************************************
if (lroot) then
if (lnokernel) then
write(*,*) "Concentrations are calculated without using kernel"
else
write(*,*) "Concentrations are calculated using kernel"
end if
end if
! Calculate particle trajectories
!********************************
......
......@@ -62,7 +62,9 @@ program flexpart
! Initialize arrays in com_mod
!*****************************
if (.not.lmpreader) call com_mod_allocate_part(maxpart_mpi)
if(.not.(lmpreader.and.lmp_use_reader)) call com_mod_allocate_part(maxpart_mpi)
! Generate a large number of random numbers
!******************************************
......@@ -78,7 +80,7 @@ program flexpart
! FLEXPART version string
flexversion_major = '10' ! Major version number, also used for species file names
! flexversion='Ver. 10 Beta MPI (2015-05-01)'
flexversion='Ver. '//trim(flexversion_major)//' Beta MPI (2015-05-01)'
flexversion='Ver. '//trim(flexversion_major)//'.1beta MPI (2016-11-02)'
verbosity=0
! Read the pathnames where input/output files are stored
......@@ -305,7 +307,7 @@ program flexpart
print*,'Initialize all particles to non-existent'
endif
if (.not.lmpreader) then
if (.not.(lmpreader.and.lmp_use_reader)) then
do j=1, size(itra1) ! maxpart_mpi
itra1(j)=-999999999
end do
......@@ -319,7 +321,7 @@ program flexpart
print*,'call readpartpositions'
endif
! readwind process skips this step
if (lmp_use_reader.and..not.lmpreader) call readpartpositions
if (.not.(lmpreader.and.lmp_use_reader)) call readpartpositions
else
if (verbosity.gt.0 .and. lroot) then
print*,'numpart=0, numparticlecount=0'
......@@ -426,6 +428,16 @@ program flexpart
end do
end do
! Inform whether output kernel is used or not
!*********************************************
if (lroot) then
if (lnokernel) then
write(*,*) "Concentrations are calculated without using kernel"
else
write(*,*) "Concentrations are calculated using kernel"
end if
end if
! Calculate particle trajectories
!********************************
......
......@@ -131,9 +131,9 @@ module com_mod
! gdomainfill .T., if domain-filling is global, .F. if not
!ZHG SEP 2015 wheather or not to read clouds from GRIB
logical :: readclouds
logical :: readclouds=.false.
!ESO DEC 2015 whether or not both clwc and ciwc are present (if so they are summed)
logical :: sumclouds
logical :: sumclouds=.false.
logical,dimension(maxnests) :: readclouds_nest, sumclouds_nest
......@@ -742,7 +742,7 @@ module com_mod
logical, parameter :: nmlout=.true.
! These variables are used to avoid having separate versions of
! files in cases where differences with MPI version is minor (eso)
! files in cases where differences with MPI version are minor (eso)
!*****************************************************************
integer :: mpi_mode=0 ! .gt. 0 if running MPI version
logical :: lroot=.true. ! true if serial version, or if MPI .and. root process
......
......@@ -180,7 +180,7 @@ subroutine conccalc(itime,weight)
! If a particle is close to the domain boundary, do not use the kernel either.
!*****************************************************************************
if ((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
if (lnokernel.or.(itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
(xl.gt.real(numxgrid-1)-0.5).or. &
(yl.gt.real(numygrid-1)-0.5)) then ! no kernel, direct attribution to grid cell
if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
......
......@@ -193,7 +193,7 @@ subroutine conccalc(itime,weight)
! If a particle is close to the domain boundary, do not use the kernel either.
!*****************************************************************************
if ((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
if (lnokernel.or.(itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
(xl.gt.real(numxgrid-1)-0.5).or. &
(yl.gt.real(numygrid-1)-0.5)) then ! no kernel, direct attribution to grid cell
if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
......
......@@ -39,6 +39,11 @@ subroutine drydepokernel(nunc,deposit,x,y,nage,kp)
! deposit amount (kg) to be deposited *
! *
!*****************************************************************************
! Changes:
! eso 10/2016: Added option to disregard kernel
!
!*****************************************************************************
use unc_mod
use par_mod
......@@ -46,7 +51,8 @@ subroutine drydepokernel(nunc,deposit,x,y,nage,kp)
implicit none
real :: x,y,deposit(maxspec),ddx,ddy,xl,yl,wx,wy,w
real(dep_prec), dimension(maxspec) :: deposit
real :: x,y,ddx,ddy,xl,yl,wx,wy,w
integer :: ix,jy,ixp,jyp,ks,nunc,nage,kp
......@@ -73,20 +79,35 @@ subroutine drydepokernel(nunc,deposit,x,y,nage,kp)
wy=0.5+ddy
endif
! If no kernel is used, direct attribution to grid cell
!******************************************************
if (lnokernel) then
do ks=1,nspec
if ((abs(deposit(ks)).gt.0).and.DRYDEPSPEC(ks)) then
if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
(jy.le.numygrid-1)) then
drygridunc(ix,jy,ks,kp,nunc,nage)= &
drygridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)
end if
end if
end do
else ! use kernel
! Determine mass fractions for four grid points
!**********************************************
do ks=1,nspec
do ks=1,nspec
if ((abs(deposit(ks)).gt.0).and.DRYDEPSPEC(ks)) then
if ((abs(deposit(ks)).gt.0).and.DRYDEPSPEC(ks)) then
if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
(jy.le.numygrid-1)) then
w=wx*wy
drygridunc(ix,jy,ks,kp,nunc,nage)= &
drygridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w
continue
endif
w=wx*wy
drygridunc(ix,jy,ks,kp,nunc,nage)= &
drygridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w
continue
endif
if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgrid-1).and. &
(jyp.le.numygrid-1)) then
......@@ -111,6 +132,7 @@ subroutine drydepokernel(nunc,deposit,x,y,nage,kp)
endif
end do
end do
end if
end subroutine drydepokernel
......@@ -49,7 +49,8 @@ subroutine drydepokernel_nest(nunc,deposit,x,y,nage,kp)
implicit none
real :: x,y,deposit(maxspec),ddx,ddy,xl,yl,wx,wy,w
real(dep_prec), dimension(maxspec) :: deposit
real :: x,y,ddx,ddy,xl,yl,wx,wy,w
integer :: ix,jy,ixp,jyp,ks,kp,nunc,nage
......
......@@ -41,7 +41,8 @@ module wind_mod
! Maximum dimensions of the input mother grids
!*********************************************
integer,parameter :: nxmax=721,nymax=361,nuvzmax=64,nwzmax=64,nzmax=64
! integer,parameter :: nxmax=721,nymax=361,nuvzmax=64,nwzmax=64,nzmax=64
integer,parameter :: nxmax=361,nymax=181,nuvzmax=64,nwzmax=64,nzmax=64
integer,parameter :: nxshift=0 ! for GFS or FNL
!*********************************************
......
......@@ -420,18 +420,19 @@ subroutine gridcheck
! Output of grid info
!********************
write(*,*)
write(*,*)
write(*,'(a,2i7)') 'Vertical levels in NCEP data: ', &
nuvz,nwz
write(*,*)
write(*,'(a)') 'Mother domain:'
write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Longitude range: ', &
xlon0,' to ',xlon0+(nx-1)*dx,' Grid distance: ',dx
write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Latitude range : ', &
ylat0,' to ',ylat0+(ny-1)*dy,' Grid distance: ',dy
write(*,*)
if (lroot) then
write(*,*)
write(*,*)
write(*,'(a,2i7)') 'Vertical levels in NCEP data: ', &
nuvz,nwz
write(*,*)
write(*,'(a)') 'Mother domain:'
write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Longitude range: ', &
xlon0,' to ',xlon0+(nx-1)*dx,' Grid distance: ',dx
write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Latitude range : ', &
ylat0,' to ',ylat0+(ny-1)*dy,' Grid distance: ',dy
write(*,*)
end if
! CALCULATE VERTICAL DISCRETIZATION OF ECMWF MODEL
! PARAMETER akm,bkm DESCRIBE THE HYBRID "ETA" COORDINATE SYSTEM
......
......@@ -418,7 +418,7 @@ shift_field_0.o: par_mod.o
timemanager.o: com_mod.o flux_mod.o netcdf_output_mod.o oh_mod.o outg_mod.o \
par_mod.o point_mod.o unc_mod.o xmass_mod.o
timemanager_mpi.o: com_mod.o flux_mod.o mpi_mod.o oh_mod.o outg_mod.o \
par_mod.o point_mod.o unc_mod.o xmass_mod.o
par_mod.o point_mod.o unc_mod.o xmass_mod.o netcdf_output_mod.o
unc_mod.o: par_mod.o
verttransform.o: cmapf_mod.o com_mod.o par_mod.o
verttransform_gfs.o: cmapf_mod.o com_mod.o par_mod.o
......
......@@ -28,7 +28,9 @@ module mean_mod
interface mean
module procedure mean_sp
module procedure mean_dp
module procedure mean_mixed_prec
module procedure mean_mixed_dss
module procedure mean_mixed_dsd
end interface mean
contains
......@@ -141,7 +143,7 @@ contains
end subroutine mean_dp
subroutine mean_mixed_prec(x_dp,xm,xs,number)
subroutine mean_mixed_dss(x_dp,xm,xs,number)
!*****************************************************************************
! *
......@@ -149,7 +151,7 @@ contains
! *
! AUTHOR: Andreas Stohl, 25 January 1994 *
! *
! Mixed precision version ESO 2016 (dp input, sp output) *
! Mixed precision version ESO 2016 (dp in, sp out, sp out) *
!*****************************************************************************
! *
! Variables: *
......@@ -191,5 +193,59 @@ contains
xs=sqrt(xaux/real(number-1,kind=sp))
endif
end subroutine mean_mixed_prec
end subroutine mean_mixed_dss
subroutine mean_mixed_dsd(x_dp,xm,xs_dp,number)
!*****************************************************************************
! *
! This subroutine calculates mean and standard deviation of a given element.*
! *
! AUTHOR: Andreas Stohl, 25 January 1994 *
! *
! Mixed precision version ESO 2016 (dp in, sp out, dp out) *
!*****************************************************************************
! *
! Variables: *
! x_dp(number) field of input data *
! xm mean *
! xs_dp standard deviation *
! number number of elements of field x_dp *
! *
! Constants: *
! eps tiny number *
! *
!*****************************************************************************
use par_mod, only: sp,dp
implicit none
integer,intent(in) :: number
real(dp), intent(in) :: x_dp(number)
real(sp), intent(out) ::xm
real(dp), intent(out) ::xs_dp
real(dp) :: xl,xq,xaux
real(dp),parameter :: eps=1.0e-30_dp
integer :: i
xl=0._dp
xq=0._dp
do i=1,number
xl=xl+x_dp(i)
xq=xq+x_dp(i)*x_dp(i)
end do
xm=xl/real(number,kind=sp)
xaux=xq-xl*xl/real(number,kind=dp)
if (xaux.lt.eps) then
xs_dp=0._dp
else
xs_dp=sqrt(xaux/real(number-1,kind=dp))
endif
end subroutine mean_mixed_dsd
end module mean_mod
......@@ -330,7 +330,7 @@ contains
! Set maxpart per process
! eso 08/2016: Increase maxpart per process, in case of unbalanced distribution
maxpart_mpi=int(mp_maxpart_factor*maxpart/mp_partgroup_np)
maxpart_mpi=int(mp_maxpart_factor*real(maxpart)/real(mp_partgroup_np))
! Do not allocate particle data arrays for readwind process
if (lmpreader.and.lmp_use_reader) then
......@@ -831,8 +831,8 @@ contains
! Increase numpart, if necessary
numpart=max(numpart,ipart)
deallocate(itra1_tmp,npoint_tmp,nclass_tmp,idt_tmp,itramem_tmp,itrasplit_tmp,&
& xtra1_tmp,ytra1_tmp,ztra1_tmp,xmass1_tmp)
deallocate(itra1_tmp,npoint_tmp,nclass_tmp,idt_tmp,itramem_tmp,&
&itrasplit_tmp,xtra1_tmp,ytra1_tmp,ztra1_tmp,xmass1_tmp)
if (mp_dev_mode) then
write(*,*) "numpart after: ", numpart
......@@ -2339,7 +2339,9 @@ contains
! implicit synchronisation between all processes takes place here
!
! TODO
! take into account nested wind fields
! NB: take into account nested wind fields by using a separate
! subroutine mpif_gf_request_nest (and separate request arrays for the
! nested variables)
!
!*******************************************************************************
! use com_mod,only: readclouds
......
......@@ -18,14 +18,15 @@
! 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 *
! - writeheader_netcdf generates file including all information previously *
! stored in separate header files *
! - concoutput_netcdf write concentration output and wet and dry deposition *
! - concoutput_netcdf write concentration output and wet and dry deposition *
! *
! Author: D. Brunner *
! *
......@@ -892,12 +893,12 @@ subroutine concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridto
gridtotal=0.
gridsigmatotal=0.
gridtotalunc=0.
wetgridtotal=0.
wetgridsigmatotal=0.
wetgridtotalunc=0.
drygridtotal=0.
drygridsigmatotal=0.
drygridtotalunc=0.
wetgridtotal=0._dep_prec
wetgridsigmatotal=0._dep_prec
wetgridtotalunc=0._dep_prec
drygridtotal=0._dep_prec
drygridsigmatotal=0._dep_prec
drygridtotalunc=0._dep_prec
do ks=1,nspec
......@@ -921,7 +922,7 @@ subroutine concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridto
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,kind=dep_prec)
wetgrid(ix,jy)=wetgrid(ix,jy)*real(nclassunc,kind=sp)
wetgridtotal=wetgridtotal+wetgrid(ix,jy)
! Calculate standard deviation of the mean
wetgridsigma(ix,jy)= &
......@@ -945,12 +946,12 @@ subroutine concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridto
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)
drygrid(ix,jy)=drygrid(ix,jy)*real(nclassunc,kind=sp)
drygridtotal=drygridtotal+drygrid(ix,jy)
! Calculate standard deviation of the mean
drygridsigma(ix,jy)= &
drygridsigma(ix,jy)* &
sqrt(real(nclassunc))
sqrt(real(nclassunc, kind=dep_prec))
drygridsigmatotal=drygridsigmatotal+ &
drygridsigma(ix,jy)
endif
......@@ -1053,8 +1054,8 @@ subroutine concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridto
if (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal
if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ &
wetgridtotal
if (drygridtotal.gt.0.) drygridtotalunc=drygridsigmatotal/ &
drygridtotal
if (drygridtotal.gt.0.) drygridtotalunc=real(drygridsigmatotal/ &
drygridtotal, kind=dep_prec)
! Dump of receptor concentrations
......@@ -1297,7 +1298,7 @@ subroutine concoutput_nest_netcdf(itime,outnum)
! Calculate standard deviation of the mean
wetgridsigma(ix,jy)= &
wetgridsigma(ix,jy)* &
sqrt(real(nclassunc))
sqrt(real(nclassunc,kind=dep_prec))
endif
! DRY DEPOSITION
......@@ -1318,7 +1319,7 @@ subroutine concoutput_nest_netcdf(itime,outnum)
! Calculate standard deviation of the mean
drygridsigma(ix,jy)= &
drygridsigma(ix,jy)* &
sqrt(real(nclassunc))
sqrt(real(nclassunc,kind=dep_prec))
endif
! CONCENTRATION OR MIXING RATIO
......
......@@ -21,7 +21,7 @@
module outg_mod
use par_mod, only: dep_prec
use par_mod, only: dep_prec, sp
implicit none
......@@ -38,8 +38,8 @@ module outg_mod
real,allocatable, dimension (:,:,:) :: densityoutgrid
real,allocatable, dimension (:,:,:) :: factor3d
real,allocatable, dimension (:,:,:) :: grid
real(dep_prec),allocatable, dimension (:,:) :: wetgrid
real(dep_prec),allocatable, dimension (:,:) :: drygrid
real(sp),allocatable, dimension (:,:) :: wetgrid
real(sp),allocatable, dimension (:,:) :: drygrid
real,allocatable, dimension (:,:,:) :: gridsigma
real(dep_prec),allocatable, dimension (:,:) :: drygridsigma
real(dep_prec),allocatable, dimension (:,:) :: wetgridsigma
......
......@@ -58,6 +58,12 @@ module par_mod
integer,parameter :: dep_prec=sp
!****************************************************************
! Set to T to disable use of kernel for concentrations/deposition
!****************************************************************
logical, parameter :: lnokernel=.false.
!***********************************************************
! Number of directories/files used for FLEXPART input/output
!***********************************************************
......
......@@ -235,12 +235,12 @@ subroutine readavailable
do i=2,numbwf
idiff=abs(wftime(i)-wftime(i-1))
if (idiff.gt.idiffmax) then
if (idiff.gt.idiffmax.and.lroot) then
write(*,*) 'FLEXPART WARNING: TIME DIFFERENCE BETWEEN TWO'
write(*,*) 'WIND FIELDS IS TOO BIG FOR TRANSPORT CALCULATION.&
&'
write(*,*) 'THEREFORE, TRAJECTORIES HAVE TO BE SKIPPED.'
else if (idiff.gt.idiffnorm) then
else if (idiff.gt.idiffnorm.and.lroot) then
write(*,*) 'FLEXPART WARNING: TIME DIFFERENCE BETWEEN TWO'
write(*,*) 'WIND FIELDS IS BIG. THIS MAY CAUSE A DEGRADATION'
write(*,*) 'OF SIMULATION QUALITY.'
......
......@@ -169,10 +169,12 @@ subroutine readcommand
901 format (a)
if (index(line,'LDIRECT') .eq. 0) then
old = .false.
write(*,*) 'COMMAND in old short format, please update to namelist format'
if (lroot) write(*,*) 'COMMAND in old short format, &
&please update to namelist format'
else
old = .true.
write(*,*) 'COMMAND in old long format, please update to namelist format'
if (lroot) write(*,*) 'COMMAND in old long format, &
&please update to namelist format'
endif
rewind(unitcommand)
......
......@@ -217,6 +217,8 @@ subroutine readreleases
rewind(unitreleases)
if (nspec.gt.maxspec) goto 994
! allocate arrays of matching size for number of species (namelist output)
deallocate(mass)
allocate(mass(nspec),stat=stat)
......@@ -573,6 +575,7 @@ subroutine readreleases
endif
endif
! Determine the release rate (particles per second) and total number
! of particles released during the simulation
!*******************************************************************
......@@ -632,6 +635,10 @@ subroutine readreleases
endif
endif
if (lroot) then
write(*,FMT='(A,ES14.7)') ' Total mass released:', sum(xmass(1:numpoint,1:nspec))
end if
return
994 write(*,*) '#####################################################'
......
......@@ -25,7 +25,7 @@ module unc_mod
implicit none
real,allocatable ,dimension (:,:,:,:,:,:,:) :: gridunc
real,allocatable, dimension (:,:,:,:,:,:,:) :: gridunc
real,allocatable, dimension (:,:,:,:,:,:,:) :: griduncn
real(dep_prec),allocatable, dimension (:,:,:,:,:,:) :: drygridunc
real(dep_prec),allocatable, dimension (:,:,:,:,:,:) :: drygriduncn
......
......@@ -39,6 +39,10 @@ subroutine wetdepokernel(nunc,deposit,x,y,nage,kp)
! deposit amount (kg) to be deposited *
! *
!*****************************************************************************
! Changes:
! eso 10/2016: Added option to disregard kernel
!
!*****************************************************************************
use unc_mod
use par_mod
......@@ -72,7 +76,19 @@ subroutine wetdepokernel(nunc,deposit,x,y,nage,kp)
wy=0.5+ddy
endif
! If no kernel is used, direct attribution to grid cell
!******************************************************
if (lnokernel) then
do ks=1,nspec
if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
(jy.le.numygrid-1)) then
wetgridunc(ix,jy,ks,kp,nunc,nage)= &
wetgridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)
end if
end do
else ! use kernel
! Determine mass fractions for four grid points
!**********************************************
......@@ -106,5 +122,6 @@ subroutine wetdepokernel(nunc,deposit,x,y,nage,kp)
wetgridunc(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w
endif
end do