Additions for netcdf output

src/FLEXPART.f90

@@ -54,6 +54,7 @@ program flexpart_NorESM
   use par_mod
   use com_mod
   use conv_mod
+  use netcdf_output_mod, only: writeheader_netcdf
 
   implicit none
 
@@ -190,6 +191,14 @@ program flexpart_NorESM
   ! and open files that are to be kept open throughout the simulation
   !******************************************************************
 
+#ifdef USE_NCF
+  if (lnetcdfout.eq.1) then 
+    call writeheader_netcdf(lnest=.false.)
+  else 
+    call writeheader
+  end if
+#endif
+  
   call writeheader
   if (nested_output.eq.1) call writeheader_nest
   open(unitdates,file=path(2)(1:length(2))//'dates')

src/com_mod.f90

@@ -51,7 +51,7 @@ module com_mod
   ! bdate                   beginning date of simulation (julian date)
   ! edate                   ending date of simulation (julian date)
 
-
+  
   integer :: ldirect,ideltas
 
   ! ldirect                 1 for forward, -1 for backward simulation
@@ -68,6 +68,11 @@ module com_mod
   ! method                  indicator which dispersion method is to be used
   ! outstep = real(abs(loutstep))
 
+  integer :: lnetcdfout=0
+  integer, parameter :: mpi_mode = 0
+  ! lnetcdfout   1 for netcdf grid output, 0 if not.
+  ! mpi_mode   0 hardcoded (added for compability with flexpart-mpi routines)
+
   real :: ctl,fine
   integer :: ifine,iout,ipout,ipin,iflux,mdomainfill
   integer :: mquasilag,nested_output,ind_source,ind_receptor

src/makefile_nilu

+SHELL = /bin/bash
+
+#  NETCDF OUTPUT
+#    To add support for output in netCDF format, append `ncf=yes` to the
+#    `make` command
+
+### Enable netCDF output?
+ifeq ($(ncf), yes)
+	NCOPT = -DUSE_NCF -lnetcdff	
+else
+	NCOPT = -UUSE_NCF
+endif
+
+F90       = gfortran
+INCPATH  = /usr/include/
+LIBPATH = /usr/lib/
+LNK = -o
+CMPL = -c
+LIBS = -lnetcdf -lnetcdff $(NCOPT)
+FFLAGS = -O0 -g -m64 -cpp -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -ffree-form \
+         -fbacktrace $(NCOPT)
+LDFLAGS  = $(FFLAGS) -L$(LIBPATH) -I$(INCPATH) $(LIBS)
+MAIN = flexpartnoresm
+#
+
+
+
+OBJECTS_NCF = netcdf_output_mod.o
+
+MODOBJS = \
+par_mod.o               com_mod.o \
+conv_mod.o              hanna_mod.o \
+interpol_mod.o          cmapf_mod.o \
+unc_mod.o               oh_mod.o \
+xmass_mod.o             flux_mod.o \
+point_mod.o             outg_mod.o \
+noresm_variables.o netcdf_output_mod.o \
+mean_mod.o
+
+
+OBJECTS = advance.o \
+allocatedumarray.o \
+allocatedumarray_parts.o \
+assignland.o \
+boundcond_domainfill.o \
+calcfluxes.o \
+calcmatrix.o \
+calcpar.o \
+calcpar_nests.o \
+calcpv.o \
+calcpv_nests.o \
+caldate_noleap.o \
+centerofmass.o \
+check_variable.o \
+clustering.o \
+conccalc.o \
+concoutput.o \
+concoutput_nest.o \
+convect43c.o \
+convmix.o \
+coordtrafo.o \
+dewpoint.o \
+distance2.o \
+distance.o \
+drydepokernel.o \
+drydepokernel_nest.o \
+dynamic_viscosity.o \
+erf.o \
+ew.o \
+FLEXPART.o \
+fluxoutput.o \
+getfields_with_dpdt.o \
+getrb.o \
+getrc.o \
+get_settling.o \
+getvdep.o \
+getvdep_nests.o \
+gridcheck.o \
+hanna1.o \
+hanna.o \
+hanna_short.o \
+init_domainfill.o \
+initial_cond_calc.o \
+initial_cond_output.o \
+initialize.o \
+interpol_all.o \
+interpol_all_nests.o \
+interpol_misslev.o \
+interpol_misslev_nests.o \
+interpol_rain.o \
+interpol_rain_nests.o \
+interpol_vdep.o \
+interpol_vdep_nests.o \
+interpol_wind.o \
+interpol_wind_nests.o \
+interpol_wind_short.o \
+interpol_wind_short_nests.o \
+juldate_noleap.o \
+mean.o \
+netcdf_output_mod.o \
+obukhov.o \
+ohreaction.o \
+openouttraj.o \
+openreceptors.o \
+outgrid_init.o \
+outgrid_init_nest.o \
+part0.o \
+partdep.o \
+partoutput.o \
+partoutput_short.o \
+pbl_profile.o \
+plumetraj.o \
+psih.o \
+psim.o \
+qvsat.o \
+raerod.o \
+random.o \
+readageclasses.o \
+readavailable.o \
+readcommand.o \
+read_delta_ps_intime.o \
+readdepo.o \
+readlanduse.o \
+readOHfield.o \
+readoutgrid.o \
+readoutgrid_nest.o \
+readpartpositions.o \
+readpaths.o \
+readreceptors.o \
+readreleases.o \
+readspecies.o \
+readwind.o \
+redist.o \
+releaseparticles.o \
+richardson.o \
+scalev.o \
+shift_field_0.o \
+shift_field.o \
+skplin.o \
+sort2.o \
+timemanager.o \
+transform_omega_etadot.o \
+verttransform.o \
+verttransform_nests.o \
+verttransform_omega.o \
+wetdepo.o \
+wetdepokernel.o \
+wetdepokernel_nest.o \
+windalign.o \
+writeheader.o \
+writeheader_nest.o \
+zenithangle.o
+
+$(MAIN): $(MODOBJS) $(OBJECTS)
+	$(F90) $(LNK) $(MAIN) *.o $(LIBS)
+
+$(OBJECTS): $(MODOBJS)
+
+%.o: %.f90
+	$(F90) $(LDFLAGS) $(CMPL) $<
+
+clean:
+	rm *.o *.mod
+

src/mean_mod.f90

+! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
+! SPDX-License-Identifier: GPL-3.0-or-later
+
+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_dss
+    module procedure mean_mixed_dsd
+
+  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
+
+    integer,intent(in) :: number
+    real(sp), intent(in) :: x_sp(number)
+    real(sp), intent(out) ::xm,xs
+    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
+
+    integer,intent(in) :: number
+    real(dp), intent(in) :: x_dp(number)
+    real(dp), intent(out) ::xm,xs
+    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_dss(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 in, sp out, sp out)              *
+!*****************************************************************************
+!                                                                            *
+! 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
+
+    integer,intent(in) :: number
+    real(dp), intent(in) :: x_dp(number)
+    real(sp), intent(out) ::xm,xs
+    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+real(x_dp(i),kind=sp)
+      xq=xq+real(x_dp(i),kind=sp)*real(x_dp(i),kind=sp)
+    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_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=real(xl,kind=sp)/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

src/par_mod.f90

@@ -44,11 +44,29 @@ module par_mod
   integer, parameter :: yearorigin=-4800 !look in juldate for this: comment by mc
 
   !****************************************************************
-  ! Parameter defining KIND parameter for "double precision"
+  ! Parameters defining KIND parameter for double/single precision
   !****************************************************************
 
   integer,parameter :: dp=selected_real_kind(P=15)
+  integer,parameter :: sp=selected_real_kind(6)
 
+  !****************************************************************
+  ! dep_prec sets the precision for deposition calculations (sp or 
+  ! dp). sp is default, dp can be used for increased precision.
+  !****************************************************************
+
+  integer,parameter :: dep_prec=sp
+
+  !****************************************************************
+  ! Set to F to disable use of kernel for concentrations/deposition
+  !****************************************************************
+
+  logical, parameter :: lusekerneloutput=.true.
+
+  !*********************************************************************
+  ! Set to T to change output units to number of particles per grid cell
+  !*********************************************************************
+  logical, parameter :: lparticlecountoutput=.false.
 
   !***********************************************************
   ! Number of directories/files used for FLEXPART input/output
@@ -65,6 +83,9 @@ module par_mod
 
   real,parameter :: pi=3.14159265, r_earth=6.371e6, r_air=287.05, ga=9.81
   real,parameter :: cpa=1004.6, kappa=0.286, pi180=pi/180., vonkarman=0.4
+  ! additional constants RLT Aug-2017
+  real,parameter :: rgas=8.31447 
+  real,parameter :: r_water=461.495
 
   ! pi                      number "pi"
   ! pi180                   pi/180.
@@ -74,6 +95,8 @@ module par_mod
   ! cpa                     specific heat for dry air
   ! kappa                   exponent of formula for potential temperature
   ! vonkarman               von Karman constant
+  ! rgas                    universal gas constant [J/mol/K]
+  ! r_water                 specific gas constant for water vapor [J/kg/K]
 
   real,parameter :: karman=0.40, href=15., convke=2.0
   real,parameter :: hmixmin=100., hmixmax=4500., turbmesoscale=0.16
@@ -205,9 +228,11 @@ module par_mod
   integer,parameter :: maxpart=20000000
   integer,parameter :: maxspec=4
 
+  real,parameter :: minmass=0.0001
 
   ! maxpart                 Maximum number of particles
   ! maxspec                 Maximum number of chemical species per release
+  ! minmass                 Terminate particles carrying less mass
 
   ! maxpoint is also set dynamically during runtime
   ! maxpoint                Maximum number of release locations
@@ -266,5 +291,8 @@ module par_mod
   integer,parameter :: unitdates=94, unitheader=90, unitshortpart=95
   integer,parameter :: unitboundcond=89
   integer,parameter :: unitdiagnostic1=72,unitdiagnostic2=73,unitdiagnostic3=74 !number from 70 to 79 reserved for diagnostic
+  integer,parameter :: unittmp=101
+! RLT
+  integer,parameter :: unitoutfactor=102
 
 end module par_mod

src/timemanager.f90

@@ -91,6 +91,7 @@ subroutine timemanager
   use oh_mod
   use par_mod
   use com_mod
+  use netcdf_output_mod
 
   implicit none
 
@@ -298,9 +299,25 @@ subroutine timemanager
 
       if ((itime.eq.loutend).and.(outnum.gt.0.)) then
         if ((iout.le.3.).or.(iout.eq.5)) then
-          call concoutput(itime,outnum,gridtotalunc, &
-               wetgridtotalunc,drygridtotalunc)
-          if (nested_output.eq.1) call concoutput_nest(itime,outnum)
+          if (lnetcdfout.eq.1) then 
+#ifdef USE_NCF
+            call concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
+#endif
+          else
+            call concoutput(itime,outnum,gridtotalunc, &
+                 wetgridtotalunc,drygridtotalunc)
+          end if
+          if (nested_output.eq.1) then
+            if (lnetcdfout.eq.1) then
+#ifdef USE_NCF
+              call concoutput_nest_netcdf(itime,outnum)
+#endif
+            else
+              call concoutput_nest(itime,outnum)              
+            endif
+          end if
+            
+
           outnum=0.
         endif
         if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)

src/unc_mod.f90

@@ -21,15 +21,24 @@
 
 module unc_mod
 
+  use par_mod, only:dep_prec
+
   implicit none
 
   real,allocatable, dimension (:,:,:,:,:,:,:) :: gridunc
+  real,allocatable, dimension (:,:,:,:,:,:,:) :: gridunc0
   real,allocatable, dimension (:,:,:,:,:,:,:) :: griduncn
   real,allocatable, dimension (:,:,:,:,:,:) :: drygridunc
   real,allocatable, dimension (:,:,:,:,:,:) :: drygriduncn
   real,allocatable, dimension (:,:,:,:,:,:) :: wetgridunc
   real,allocatable, dimension (:,:,:,:,:,:) :: wetgriduncn
 
+  ! For sum of individual contributions, used for the MPI version
+  real(dep_prec),allocatable, dimension (:,:,:,:,:,:) :: drygridunc0
+  real(dep_prec),allocatable, dimension (:,:,:,:,:,:) :: drygriduncn0
+  real(dep_prec),allocatable, dimension (:,:,:,:,:,:) :: wetgridunc0
+  real(dep_prec),allocatable, dimension (:,:,:,:,:,:) :: wetgriduncn0
+
   real,allocatable, dimension (:,:,:,:,:) :: init_cond
 
 end module unc_mod