Commit 6a678e3f authored by Espen Sollum's avatar Espen Sollum
Browse files

Added option to use double precision for calculating the deposition fields

parent fddc6ec7
......@@ -62,6 +62,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
use outg_mod
use par_mod
use com_mod
use mean_mod
implicit none
......@@ -89,9 +90,10 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
!integer sparse_dump_i(numxgrid*numygrid*numzgrid)
!real sparse_dump_u(numxgrid*numygrid*numzgrid)
real :: auxgrid(nclassunc),gridtotal,gridsigmatotal,gridtotalunc
real :: wetgridtotal,wetgridsigmatotal,wetgridtotalunc
real :: drygridtotal,drygridsigmatotal,drygridtotalunc
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
......
......@@ -67,6 +67,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
use par_mod
use com_mod
use mpi_mod
use mean_mod
implicit none
......@@ -94,9 +95,10 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
!integer sparse_dump_i(numxgrid*numygrid*numzgrid)
!real sparse_dump_u(numxgrid*numygrid*numzgrid)
real :: auxgrid(nclassunc),gridtotal,gridsigmatotal,gridtotalunc
real :: wetgridtotal,wetgridsigmatotal,wetgridtotalunc
real :: drygridtotal,drygridsigmatotal,drygridtotalunc
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
......
......@@ -60,6 +60,7 @@ subroutine concoutput_nest(itime,outnum)
use outg_mod
use par_mod
use com_mod
use mean_mod
implicit none
......@@ -87,7 +88,7 @@ subroutine concoutput_nest(itime,outnum)
!integer sparse_dump_i(numxgrid*numygrid*numzgrid)
!real sparse_dump_u(numxgrid*numygrid*numzgrid)
real :: auxgrid(nclassunc)
real(dep_prec) :: auxgrid(nclassunc)
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
......
......@@ -64,6 +64,7 @@ subroutine concoutput_nest(itime,outnum)
use par_mod
use com_mod
use mpi_mod
use mean_mod
implicit none
......@@ -91,7 +92,7 @@ subroutine concoutput_nest(itime,outnum)
!integer sparse_dump_i(numxgrid*numygrid*numzgrid)
!real sparse_dump_u(numxgrid*numygrid*numzgrid)
real :: auxgrid(nclassunc)
real(dep_prec) :: auxgrid(nclassunc)
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
......
......@@ -62,6 +62,7 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
use outg_mod
use par_mod
use com_mod
use mean_mod
implicit none
......@@ -89,9 +90,10 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
!integer sparse_dump_i(numxgrid*numygrid*numzgrid)
!real sparse_dump_u(numxgrid*numygrid*numzgrid)
real :: auxgrid(nclassunc),gridtotal,gridsigmatotal,gridtotalunc
real :: wetgridtotal,wetgridsigmatotal,wetgridtotalunc
real :: drygridtotal,drygridsigmatotal,drygridtotalunc
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
......
......@@ -68,6 +68,7 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
use par_mod
use com_mod
use mpi_mod
use mean_mod
implicit none
......@@ -95,9 +96,10 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
!integer sparse_dump_i(numxgrid*numygrid*numzgrid)
!real sparse_dump_u(numxgrid*numygrid*numzgrid)
real :: auxgrid(nclassunc),gridtotal,gridsigmatotal,gridtotalunc
real :: wetgridtotal,wetgridsigmatotal,wetgridtotalunc
real :: drygridtotal,drygridsigmatotal,drygridtotalunc
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
......
......@@ -60,6 +60,7 @@ subroutine concoutput_surf_nest(itime,outnum)
use outg_mod
use par_mod
use com_mod
use mean_mod
implicit none
......@@ -87,7 +88,7 @@ subroutine concoutput_surf_nest(itime,outnum)
!integer sparse_dump_i(numxgrid*numygrid*numzgrid)
!real sparse_dump_u(numxgrid*numygrid*numzgrid)
real :: auxgrid(nclassunc)
real(dep_prec) :: auxgrid(nclassunc)
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
......
......@@ -63,6 +63,7 @@ subroutine concoutput_surf_nest(itime,outnum)
use par_mod
use com_mod
use mpi_mod
use mean_mod
implicit none
......@@ -90,7 +91,7 @@ subroutine concoutput_surf_nest(itime,outnum)
!integer sparse_dump_i(numxgrid*numygrid*numzgrid)
!real sparse_dump_u(numxgrid*numygrid*numzgrid)
real :: auxgrid(nclassunc)
real(dep_prec) :: auxgrid(nclassunc)
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
......
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
......@@ -81,7 +81,7 @@ interpol_mod.o cmapf_mod.o \
unc_mod.o oh_mod.o \
xmass_mod.o flux_mod.o \
point_mod.o outg_mod.o \
random_mod.o
mean_mod.o random_mod.o
MPI_MODOBJS = \
mpi_mod.o
......@@ -166,7 +166,7 @@ pbl_profile.o readOHfield.o\
juldate.o \
interpol_vdep.o interpol_rain.o \
hanna.o wetdepokernel.o \
mean.o wetdepo.o \
wetdepo.o \
hanna_short.o windalign.o \
hanna1.o \
gridcheck_nests.o \
......@@ -292,18 +292,18 @@ cmapf_mod.o: par_mod.o
com_mod.o: par_mod.o
conccalc.o: com_mod.o outg_mod.o par_mod.o unc_mod.o
conccalc_mpi.o: com_mod.o mpi_mod.o outg_mod.o par_mod.o unc_mod.o
concoutput.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o
concoutput.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
concoutput_mpi.o: com_mod.o mpi_mod.o outg_mod.o par_mod.o point_mod.o \
unc_mod.o
concoutput_nest.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o
unc_mod.o mean_mod.o
concoutput_nest.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
concoutput_nest_mpi.o: com_mod.o mpi_mod.o outg_mod.o par_mod.o point_mod.o \
unc_mod.o
concoutput_surf.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o
unc_mod.o mean_mod.o
concoutput_surf.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
concoutput_surf_mpi.o: com_mod.o mpi_mod.o outg_mod.o par_mod.o point_mod.o \
unc_mod.o
concoutput_surf_nest.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o
unc_mod.o mean_mod.o
concoutput_surf_nest.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
concoutput_surf_nest_mpi.o: com_mod.o mpi_mod.o outg_mod.o par_mod.o \
point_mod.o unc_mod.o
point_mod.o unc_mod.o mean_mod.o
conv_mod.o: par_mod.o
convect43c.o: conv_mod.o par_mod.o
convmix.o: com_mod.o conv_mod.o flux_mod.o par_mod.o
......@@ -358,13 +358,15 @@ interpol_wind_nests.o: com_mod.o interpol_mod.o par_mod.o
interpol_wind_short.o: com_mod.o interpol_mod.o par_mod.o
interpol_wind_short_nests.o: com_mod.o interpol_mod.o par_mod.o
juldate.o: par_mod.o
mean_mod.o: par_mod.o
mpi_mod.o: com_mod.o par_mod.o unc_mod.o
netcdf_output_mod.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o
netcdf_output_mod.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
obukhov.o: par_mod.o
obukhov_gfs.o: par_mod.o
ohreaction.o: com_mod.o oh_mod.o par_mod.o
openouttraj.o: com_mod.o par_mod.o point_mod.o
openreceptors.o: com_mod.o par_mod.o
outg_mod.o: par_mod.o
outgrid_init.o: com_mod.o flux_mod.o oh_mod.o outg_mod.o par_mod.o unc_mod.o
outgrid_init_nest.o: com_mod.o outg_mod.o par_mod.o unc_mod.o
par_mod.o : $(wind_mod)
......@@ -375,7 +377,7 @@ partoutput_mpi.o: com_mod.o mpi_mod.o par_mod.o
partoutput_short.o: com_mod.o par_mod.o
partoutput_short_mpi.o: com_mod.o mpi_mod.o par_mod.o
pbl_profile.o: par_mod.o
plumetraj.o: com_mod.o par_mod.o point_mod.o
plumetraj.o: com_mod.o par_mod.o point_mod.o mean_mod.o
psih.o: par_mod.o
psim.o: par_mod.o
raerod.o: par_mod.o
......
!**********************************************************************
! 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 mean(x,xm,xs,number)
!*****************************************************************************
! *
! This subroutine calculates mean and standard deviation of a given element.*
! *
! AUTHOR: Andreas Stohl, 25 January 1994 *
! *
!*****************************************************************************
! *
! Variables: *
! x(number) field of input data *
! xm mean *
! xs standard deviation *
! number number of elements of field x *
! *
! Constants: *
! eps tiny number *
! *
!*****************************************************************************
implicit none
integer :: number,i
real :: x(number),xm,xs,xl,xq,xaux
real,parameter :: eps=1.0e-30
xl=0.
xq=0.
do i=1,number
xl=xl+x(i)
xq=xq+x(i)*x(i)
end do
xm=xl/real(number)
xaux=xq-xl*xl/real(number)
if (xaux.lt.eps) then
xs=0.
else
xs=sqrt(xaux/real(number-1))
endif
end subroutine mean
!**********************************************************************
! 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/>. *
!**********************************************************************
module mean_mod
public
! Interface to select single or double precision version of the 'mean'
! function from type of input arguments ("function overloading")
!********************************************************************
interface mean
module procedure mean_sp
module procedure mean_dp
module procedure mean_mixed_prec
end interface mean
contains
subroutine mean_sp(x_sp,xm,xs,number)
!*****************************************************************************
! *
! This subroutine calculates mean and standard deviation of a given element.*
! *
! AUTHOR: Andreas Stohl, 25 January 1994 *
! *
! Single precision version ESO 2016 *
!*****************************************************************************
! *
! Variables: *
! x_sp(number) field of input data *
! xm mean *
! xs standard deviation *
! number number of elements of field x_sp *
! *
! Constants: *
! eps tiny number *
! *
!*****************************************************************************
use par_mod, only: sp
implicit none
! integer :: number,i
! real(sp) :: x_sp(number),xm,xs,xl,xq,xaux
! real(sp),parameter :: eps=1.0e-30
real(sp), intent(in) :: x_sp(number)
real(sp), intent(out) ::xm,xs
integer,intent(in) :: number
real(sp) :: xl,xq,xaux
real(sp),parameter :: eps=1.0e-30
integer :: i
xl=0.
xq=0.
do i=1,number
xl=xl+x_sp(i)
xq=xq+x_sp(i)*x_sp(i)
end do
xm=xl/real(number,kind=sp)
xaux=xq-xl*xl/real(number,kind=sp)
if (xaux.lt.eps) then
xs=0.
else
xs=sqrt(xaux/real(number-1,kind=sp))
endif
end subroutine mean_sp
subroutine mean_dp(x_dp,xm,xs,number)
!*****************************************************************************
! *
! This subroutine calculates mean and standard deviation of a given element.*
! *
! AUTHOR: Andreas Stohl, 25 January 1994 *
! *
! Double precision version ESO 2016 *
!*****************************************************************************
! *
! Variables: *
! x_dp(number) field of input data *
! xm mean *
! xs standard deviation *
! number number of elements of field x_dp *
! *
! Constants: *
! eps tiny number *
! *
!*****************************************************************************
use par_mod, only: dp
implicit none
real(dp), intent(in) :: x_dp(number)
real(dp), intent(out) ::xm,xs
integer,intent(in) :: number
real(dp) :: xl,xq,xaux
real(dp),parameter :: eps=1.0e-30
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=dp)
xaux=xq-xl*xl/real(number,kind=dp)
if (xaux.lt.eps) then
xs=0._dp
else
xs=sqrt(xaux/real(number-1,kind=dp))
endif
end subroutine mean_dp
subroutine mean_mixed_prec(x_dp,xm,xs,number)
!*****************************************************************************
! *
! This subroutine calculates mean and standard deviation of a given element.*
! *
! AUTHOR: Andreas Stohl, 25 January 1994 *
! *
! Mixed precision version ESO 2016 (dp input, sp output) *
!*****************************************************************************
! *
! Variables: *
! x_dp(number) field of input data *
! xm mean *
! xs standard deviation *
! number number of elements of field x_dp *
! *
! Constants: *
! eps tiny number *
! *
!*****************************************************************************
use par_mod, only: sp,dp
implicit none
real(dp), intent(in) :: x_dp(number)
real(sp), intent(out) ::xm,xs
integer,intent(in) :: number
real(sp) :: xl,xq,xaux
real(sp),parameter :: eps=1.0e-30
integer :: i
xl=0._sp
xq=0._sp
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=sp)
if (xaux.lt.eps) then
xs=0._sp
else
xs=sqrt(xaux/real(number-1,kind=sp))
endif
end subroutine mean_mixed_prec
end module mean_mod
......@@ -52,6 +52,7 @@ module mpi_mod
! loops over particles. Will be all processes *
! unless a dedicated process runs getfields/readwind *
! lmp_sync If .false., use asynchronous MPI *
! mp_cp Real precision to use for deposition fields *
! *
! *
! *
......@@ -83,7 +84,7 @@ module mpi_mod
integer :: mp_seed=0
integer, parameter :: mp_sp=MPI_REAL4, mp_dp=MPI_REAL8
integer, parameter :: mp_pp=mp_sp
integer :: mp_cp
integer, parameter :: id_root=0 ! master process
! MPI tags/requests for send/receive operation
......@@ -162,14 +163,14 @@ contains
! mpi_mode default 0, set to 2/3 if running MPI version
! mp_np number of running processes, decided at run-time
!***********************************************************************
use par_mod, only: maxpart, numwfmem
use par_mod, only: maxpart, numwfmem, dep_prec
use com_mod, only: mpi_mode
implicit none
integer :: i,j,s,addmaxpart=0
! each process gets an ID (mp_pid) in the range 0,..,mp_np-1
! Each process gets an ID (mp_pid) in the range 0,..,mp_np-1
call MPI_INIT(mp_ierr)
if (mp_ierr /= 0) goto 100
call MPI_COMM_RANK(MPI_COMM_WORLD, mp_pid, mp_ierr)
......@@ -178,7 +179,7 @@ contains
if (mp_ierr /= 0) goto 100
! this variable is used to handle subroutines common to parallel/serial version
! Variable mpi_mode is used to handle subroutines common to parallel/serial version
if (lmp_sync) then
mpi_mode=2 ! hold 2 windfields in memory
else
......@@ -189,6 +190,19 @@ contains
lroot = .false.
end if
! Set MPI precision to use for transferring deposition fields
!************************************************************
if (dep_prec==dp) then
mp_cp = MPI_REAL8
if (lroot) write(*,*) 'Using double precision for deposition fields'
else if (dep_prec==sp) then
mp_cp = MPI_REAL4
if (lroot) write(*,*) 'Using single precision for deposition fields'
else
write(*,*) 'ERROR: something went wrong setting MPI real precision'
stop
end if
! Check for sensible combination of parameters
!*********************************************
......@@ -1670,13 +1684,13 @@ contains
end if
if ((WETDEP).and.(ldirect.gt.0)) then
call MPI_Reduce(wetgridunc, wetgridunc0, grid_size2d, mp_sp, MPI_SUM, id_root, &
call MPI_Reduce(wetgridunc, wetgridunc0, grid_size2d, mp_cp, MPI_SUM, id_root, &
& mp_comm_used, mp_ierr)
if (mp_ierr /= 0) goto 600
end if
if ((DRYDEP).and.(ldirect.gt.0)) then
call MPI_Reduce(drygridunc, drygridunc0, grid_size2d, mp_sp, MPI_SUM, id_root, &
call MPI_Reduce(drygridunc, drygridunc0, grid_size2d, mp_cp, MPI_SUM, id_root, &
& mp_comm_used, mp_ierr)
if (mp_ierr /= 0) goto 600
end if
......@@ -1746,13 +1760,13 @@ contains
end if
if ((WETDEP).and.(ldirect.gt.0)) then
call MPI_Reduce(wetgriduncn, wetgriduncn0, grid_size2d, mp_sp, MPI_SUM, id_root, &
call MPI_Reduce(wetgriduncn, wetgriduncn0, grid_size2d, mp_cp, MPI_SUM, id_root, &
& mp_comm_used, mp_ierr)
if (mp_ierr /= 0) goto 600
end if
if ((DRYDEP).and.(ldirect.gt.0)) then
call MPI_Reduce(drygriduncn, drygriduncn0, grid_size2d, mp_sp, MPI_SUM, id_root, &
call MPI_Reduce(drygriduncn, drygriduncn0, grid_size2d, mp_cp, MPI_SUM, id_root, &
& mp_comm_used, mp_ierr)
if (mp_ierr /= 0) goto 600
end if
......
......@@ -35,6 +35,10 @@
! - added option to not writeout releases information by changing
! switch write_releases
! - additional updates for FLEXPART 9.x
!
! ESO 2016
! - Deposition fields can be calculated in double precision, see variable
! 'dep_prec' in par_mod
!*****************************************************************************
......@@ -47,7 +51,7 @@ module netcdf_output_mod
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,&
use par_mod, only: dep_prec, sp, dp, maxspec, maxreceptor, nclassunc,&
unitoutrecept,unitoutreceptppt, nxmax,unittmp
use com_mod, only: path,length,ldirect,ibdate,ibtime,iedate,ietime, &
loutstep,loutaver,loutsample,outlon0,outlat0,&
......@@ -70,6 +74,8 @@ module netcdf_output_mod
nested_output, ipout, surf_only, linit_cond, &
flexversion,mpi_mode
use mean_mod
implicit none
! include 'netcdf.inc'
......@@ -722,16 +728,21 @@ subroutine concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridto