View source for Flexible precision in GEOS-Chem

== Overview ==

=== Flexible precision was introduced in Fortran-90===

In most Fortran codes (including GEOS-Chem) you will see declarations such as:

 ! Integers
 INTEGER   :: I, J   ! 4-byte integer  
 INTEGER*4 :: K, L   ! 4-byte integer
 INTEGER*8 :: M, N   ! 8-byte integer
 
 ! Floating point
 REAL      :: A, B   ! 4-byte floating point
 REAL*4    :: C, D   ! 4-byte floating point
 REAL*8    :: E, F   ! 8-byte floating point 

etc.  Note that:

#On most compilers, <tt>INTEGER</tt> refers to a 4-byte integer.  You can make this default to an 8-byte integer by compiling with <tt>-i8</tt>.  In most circumstances it is OK to use 4-byte integers, unless you need to point to a memory location or are reading an 8-byte integer from a netCDF file.
#On most compilers, <tt>REAL</tt> refers to a 4-byte floating-point.  You can make this default to an 8-byte floating point by compiling with <tt>-r8</tt>.
#In some older Fortran codes, you will see the term <tt>DOUBLE PRECISION</tt>.  This is the same as <tt>REAL*8</tt> -- it is an 8-byte floating point.

Fortran 90 introduced a new precision concept.  You can replace these fixed data types with [https://www.nsc.liu.se/~boein/f77to90/c13.html declarations of arbitrary precision].  This is done with the <tt>SELECTED_REAL_KIND</tt> and <tt>SELECTED_INT_KIND</tt> functions, [https://www.nsc.liu.se/~boein/f77to90/a5.html#section6 which are described in more detail here].

--[[User:Bmy|Bob Y.]] 15:01, 7 November 2014 (EST)

=== Why are we implementing flexible precision in GEOS-Chem? ===

Long story short: we need to do this in order to interface GEOS-Chem into the NASA GEOS-5 GCM (and other GCM's) more efficiently.

Most of the floating-point variables in GEOS-Chem are declared as <tt>REAL*8</tt>.  We wrote the GEOS-Chem code in this manner, starting many years ago, when we were running on coarser-resolution grids (i.e. 4&deg; x 5&deg;).  Back then, memory was generally not an issue.

But many GCMs&mdash;including the GEOS-5 GCM&mdash;declare floating-point variables as <tt>REAL*4</tt>.  Because GCM's typically operate on very fine horizontal grids, conserving memory is of paramount concern.  

When we connect GEOS-Chem to the GEOS-5 GCM, for example, we have to copy <tt>REAL*4</tt> data from the GCM (such as the met fields, surface parameters, and other relevant quantities) into GEOS-Chem's <tt>REAL*8</tt> arrays.  This copying operation is very costly, as it requires extra memory and CPU cycles.  But if we can transform GEOS-Chem's <tt>REAL*8</tt> arrays into <tt>REAL*4</tt> arrays, then we could just point GEOS-Chem's arrays to the GCM arrays without having to do all of the extra operations associated with the copying process. 

'''Therefore, our goal is to recode GEOS-Chem so that you can select the floating-point precision that you want to use (either <tt>REAL*4</tt> or <tt>REAL*8</tt>) at compile time.'''  If you are going to connect GEOS-Chem to the GEOS-5 GCM, you can request all of the floating point variables to be declared as <tt>REAL*4</tt>, in order to match the variables in the GEOS-5 GCM.  But if you are using the "traditional" serial GEOS-Chem, you can request that the floating point varaibles be declared as <tt>REAL*8</tt>, for backwards compatibility with prior code.

--[[User:Bmy|Bob Y.]] 15:17, 7 November 2014 (EST)

== Methodology ==

We followed this procedure in order to implement flexible precision into GEOS-Chem:

'''1.''' We defined  module <tt>Headers/precision_mod.F</tt>.  This module defines a parameter that will be used to specify the precision of variables in other parts of the code.

 !------------------------------------------------------------------------------
 !                  GEOS-Chem Global Chemical Transport Model                  !
 !------------------------------------------------------------------------------
 !BOP
 !
 ! !MODULE: precision_mod.F
 !
 ! !DESCRIPTION: Module PRECISION\_MOD is used to change the precision of
 !  many variables throughout GEOS-Chem at compile-time.
 !\\
 !\\
 ! !INTERFACE:
 !
       MODULE PRECISION_MOD
 !
 ! !USES:
 ! 
       IMPLICIT NONE
       PRIVATE
 !
 ! !DEFINED PARAMETERS:
 ! 
       !=================================================================
       ! Set parameters for floating precision
       !
       ! FP will be set to either 4-byte or 8-byte precision at compile 
       ! time.  Most variables can now  declared with REAL(fp).
       !=================================================================
 #if defined( USE_REAL8 )
 
       ! Use 8-byte floating point precision when asked.
       INTEGER, PARAMETER, PUBLIC :: fp = KIND( REAL( 0.0, 8 ) )
 
 #else
 
       ! Use 4-byte floating point by default.
       INTEGER, PARAMETER, PUBLIC :: fp = KIND( REAL( 0.0, 4 ) )
 
 #endif
 
       !=================================================================
       ! Set parameters for fixed precision
       !
       ! Not all variables can be converted into the flexible precision.  
       ! Some may have to be still declared as either 4-byte or 8-byte 
       ! floating point.  Use these parameters for such variables.
       !=================================================================
 
       ! KIND parameter for 4-byte precision
       INTEGER, PARAMETER, PUBLIC :: f4 = KIND( REAL( 0.0, 4 ) )
       
       ! KIND parameter for 8-byte precision
       INTEGER, PARAMETER, PUBLIC :: f8 = KIND( REAL( 0.0, 8 ) )
 !
 ! !REMARKS:
 !  This module is designed to help avoid hard-coding precision.
 !
 ! !REVISION HISTORY:
 !  (1 ) Created. (myannetti, 11/04/14)
 !  23 Nov 2016 - R. Yantosca - Now rewrite KIND definitions to prevent 4-byte
 !                              and 8-byte variables from being elevated
 !                              when using -r8 (or equivalent flags)
 !EOP
 !-----------------------------------------------------------------------------
 !BOC
       END MODULE PRECISION_MOD
 !EOC

Instead of having to figure out the proper settings with <tt>SELECTED_REAL_KIND</tt>, we can just use the KIND command to return that for us.  

*<tt>KIND( REAL( 0.0, 8 ) )</tt> returns the proper "kind" value to define an 8-byte floating point
*<tt>KIND( REAL( 0.0, 4 ) )</tt> returns the proper "kind" value to define a 4-byte floating point

This value returned by the <tt>KIND</tt> function is saved in the constant named <tt>fp</tt>.  (NOTE: This stands for "flexible precision".  It does not have any connection to the [[GEOS-FP]] meteorology.)

Note that we have used an <tt>#if defined</tt> block to define the value of <tt>fp</tt>.  If we compile with <tt>-DUSE_REAL8</tt>, then the <tt>fp</tt> can be used to declare 8-byte floating-point variables.  Otherwise, <tt>fp</tt> can be used to declare 4-byte variables by default.

'''2.''' Add <tt>precision_mod.F</tt> to the dependencies listing in the Headers/Makefile.  Add this line:

 precision_mod.o: precision_mod.F

'''3.''' In the <tt>Makefille_header.mk</tt>, we added a new Makefile variable named <tt>PRECISION</tt>:

 # %%%%% Default to 8-byte precision unless specified otherwise %%%%%
 ifndef PRECISION
  PRECISION     :=8
 endif

This variable is set to 8 by default (because for now, we want to compile the "traditional" serial GEOS-Chem with <tt>REAL*8</tt> floating point precision, as it has always been compiled.  

<tt>PRECISION</tt> is used again further down in the <tt>Makefile_header.mk</tt> to add a C-preprocessor switch:

 # Add flexible precision declaration
 ifeq ($(PRECISION),8)
 USER_DEFS      += -DUSE_REAL8
 endif

The <tt>-DUSE_REAL8</tt> will define the <tt>USE_REAL8</tt> C-preprocessor switch, which in turn will automatically pick 8-byte floating point precision.


'''4.''' At the top of each GEOS-Chem module or routine (typically in the <tt>!USES:</tt> comment section), you can place a reference to <tt>precision_mod.F90</tt>.   For example, at the top of <tt>GeosCore/carbon_mod.F</tt>, you would add the line in RED:

 !------------------------------------------------------------------------------
 !                  GEOS-Chem Global Chemical Transport Model                  !
 !------------------------------------------------------------------------------
 !BOP
 !     
 ! !MODULE: carbon_mod
 !     
 ! !DESCRIPTION: Module CARBON\_MOD contains arrays and routines for performing
 !  a carbonaceous aerosol simulation.  Original code taken from Mian Chin's 
 !  GOCART model and modified accordingly. (rjp, bmy, 4/2/04, 6/30/10)
 !\\   
 !\\   
 ! !INTERFACE: 
 !
       MODULE CARBON_MOD
 !
 ! !USES:
 !
 !
       USE HCO_ERROR_MOD       ! For real precisions (hp)
       USE PRECISION_MOD       ! For GEOS-Chem precision (fp)
 
       IMPLICIT NONE
       PRIVATE

Note that HEMCO has its own precision parameters.  We'll leave those alone, because HEMCO ships as a separate package as well.


'''5.''' Look for all variables in the module that are declared as <tt>REAL*8</tt>.  Replace the <tt>REAL*8</tt> text with <tt>REAL(fp)</tt> instead.  So, in the above example, these lines:

       <span style="color:red">REAL*8</span>, ALLOCATABLE :: ANTH_BLKC(:,:,:)
       <span style="color:red">REAL*8</span>, ALLOCATABLE :: ANTH_ORGC(:,:,:)
       <span style="color:red">REAL*8</span>, ALLOCATABLE :: BIOB_BLKC(:,:,:)

would become 

       <span style="color:green">REAL(fp)</span>, ALLOCATABLE :: ANTH_BLKC(:,:,:)
       <span style="color:green">REAL(fp)</span>, ALLOCATABLE :: ANTH_ORGC(:,:,:)
       <span style="color:green">REAL(fp)</span>, ALLOCATABLE :: BIOB_BLKC(:,:,:)


'''6.''' IMPORTANT NOTE!  Any literal constants in scientific notation made with the Fortran <tt>d</tt> exponents have to be changed to <tt>e</tt>.  Also, the text <tt>_fp</tt> has to be appended to the exponent.  This tells Fortran that we are using a customized precision definition.  For example, the code:

       ! Molecules OH  per kg OH [molec/kg]
       REAL*8,  PARAMETER  :: XNUMOL_OH  = 6.022<span style="color:red">d23</span> / 17<span style="color:red">d-3</span>
       REAL*8,  PARAMETER  :: CM3PERM3   = 1.<span style="color:red">d6</span>
       REAL*8,  PARAMETER  :: TINY       = TINY(1.0)

would become instead:

       ! Molecules OH  per kg OH [molec/kg]
       REAL(fp), PARAMETER :: XNUMOL_OH  = 6.022<span style="color:green">e+23_fp</span> / 17<span style="color:green">e-3_fp</span>
       REAL(fp), PARAMETER :: CM3PERM3   = 1.<span style="color:green">e+6_fp</span>
       REAL(fp), PARAMETER :: TINY       = TINY(1.0<span style="color:green">_fp</span>)

'''7.''' We repeated the process in steps 4-6 for each GEOS-Chem source code file.  We typically modify a file or two at a time, and then run a difference test.  A difference test compares the code we are editing to a code with known behavior, such as the last accepted benchmarked version.

--[[User:Bmy|Bob Yantosca]] ([[User talk:Bmy|talk]]) 18:10, 23 November 2016 (UTC)

== Compiling GEOS-Chem ==

To compile GEOS-Chem for 8-byte floating-point precision, just use the same commands as you always would:

 make -j4 MET=geosfp GRID=4x5 TRACEBACK=y ...

(The ... denotes other compiler options, as described in [[Compiling GEOS-Chem]].)

But to compile with 4-byte floating-point precision, you must now use the <tt>PRECISION</tt> keyword:

 make -j4 MET=geosfp GRID=4x5 TRACEBACK=y <span style="color:green">PRECISION=4</span> ...

Eventually, <tt>PRECISION=4</tt> will be automatically set if you build GEOS-Chem with the <tt>hpc</tt> option.  This option compiles GEOS-Chem for use with the ESMF environment (such as is used in the GEOS-5 GCM).  Then you can just use this command:

 make -j4 MET=geosfp GRID=4x5 TRACEBACK=y ... hpc

--[[User:Bmy|Bob Y.]] 15:41, 7 November 2014 (EST)

== Validation ==

The following sections summarize our evaluation of GEOS-Chem with the <tt>PRECISION=4</tt> option:

=== Benchmarks ===

Bob Yantosca ([[GCST]]) performed a 1-month benchmark simulation with the <tt>PRECISION=4</tt> option.  The results are posted on 

Basic takeaways: The 1-month benchmark simulation with <tt>PRECISION=4</tt>:

# Used about 40% less memory than the [[GEOS-Chem v11-01 benchmark history#v11-01f|v11-01f 1-month benchmark simulation (compiled with <tt>PRECISION=8</tt>)]].
# Finished 20 minutes faster than the v11-01f 1-month benchmark simulation.
# Shows negligible differences w/r/t the v11-01f 1-month benchmark.

--[[User:Bmy|Bob Yantosca]] ([[User talk:Bmy|talk]]) 17:48, 22 February 2017 (UTC)

=== Mass conservation ===

The plots below illustrate the evolution of total CO2 tracer mass in <tt>geosfp_2x25_masscons</tt> simulations (met year 2013) using GEOS-Chem "Classic".  The <span style="color:red">'''RED LINE'''</span> indicates the initial total mass, for reference.

'''Figure 1: Using <tt>PRECISION=4</tt> and [[Boundary_layer_mixing#TURBDAY|TURBDAY]] (aka full) PBL mixing'''

When the <tt>PRECISION=4</tt> option is activated, we observe a consistent loss of mass starting after approximately one month.  Differences start appearing in the 3rd decimal place.

 TAU range    [hrs since 1/1/1985] :         245454.00000        254208.00000
 Initial mass [kg                ] : 2.87<span style="color:green">55168</span>0000000e+15
 Mass range   [kg                ] : 2.87<span style="color:red">48114</span>0000000e+15 2.87<span style="color:blue">56022</span>0000000e+15

[[Image:Mass_cons_prec4_1yr.png]]

'''Figure 2: Using <tt>PRECISION=8</tt> and [[Boundary_layer_mixing#TURBDAY|TURBDAY]] (aka full) PBL mixing''' 

With <tt>PRECISION=8</tt>, differences start to appear in the 12th decimal place.  This indicates numerical noise around a constant total tracer mass.

 TAU range    [hrs since 1/1/1985] :         245454.00000        254208.00000
 Initial mass [kg                ] : 2.87745967221<span style="color:green">413</span>e+15
 Mass range   [kg                ] : 2.87745967221<span style="color:red">277</span>e+15 2.87745967221<span style="color:blue">468</span>e+15

[[Image:Mass_cons_prec8_1yr.png]]

'''Figure 3: Using <tt>PRECISION=4</tt> and [[Boundary_layer_mixing#VDIFF|VDIFF]] (aka non-local) PBL mixing'''.  

Using <tt>PRECISION=4</tt> and the non-local boundary layer mixing scheme (VDIFF) results in a similar picture to Figure 1.  There is a consistent, significant mass loss with time.  Again, differences start to appear in the 3rd decimal place.

 TAU range    [hrs since 1/1/1985] :         245454.00000        254208.00000
 Initial mass [kg                ] : 2.87<span style="color:green">55171</span>0000000e+15
 Mass range   [kg                ] : 2.87<span style="color:red">44329</span>0000000e+15 2.87<span style="color:blue">55509</span>0000000e+15

[[Image:Mass_cons_prec4_1yr_vdiff.png]]

'''Figure 4: Using <tt>PRECISION=8</tt> and [[Boundary_layer_mixing#VDIFF|VDIFF]] (aka non-local) PBL mixing'''. 

Using <tt>PRECISION=8</tt> and the non-local boundary layer mixing scheme (VDIFF), we see mass conservation to the 8th decimal place.  This indicates numerical noise around a constant tracer mass.

 TAU range    [hrs since 1/1/1985] :         245454.00000        254208.00000
 Initial mass [kg                ] : 2.87745967<span style="color:green">264367</span>e+15
 Mass range   [kg                ] : 2.87745967<span style="color:red">257390</span>e+15 2.87745967<span style="color:blue">694125</span>e+15

[[Image:Mass_cons_prec8_1yr_vdiff.png]]

'''Summary'''

The plots above show that mass conservation is better achieved with the <tt>PRECISION=8</tt> option than with <tt>PRECISION=4</tt>.  

'''''[[User:Sebastian D. Eastham|Seb Eastham]] wrote:'''''

<blockquote>I would also advise that we tell users who want to do a simulation involving long-lived species (think stratosphere or CO2) that they might want to consider using <tt>PRECISION=8</tt>, which is the current GEOS-Chem default setting. One concern with the mass conservation error in <tt>PRECISION=4</tt> (or actually at any precision) is that I think it will be a stationary random variable&mdash;there's no guarantee that the mean of the error over time will tend towards zero. This isn't a problem when the error at each step is ~1e-10, but could be a bigger problem when it's 1e-5 because we have no guarantee that it won't accumulate in one direction out of random bad luck.</blockquote>

'''''Daniel Jacob replied:'''''

<blockquote>The performance gains from using <tt>PRECISION=4</tt> are most relevant for full-chemistry applications, not for CO2 applications.</blockquote>

'''''[[User:bmy|Bob Yantosca]] replied:'''''

<blockquote>We could...set <tt>PRECISION=4</tt> as the default option when you generate any of the full-chemistry run directories (Standard, Tropchem, UCX, SOA, SOA-SVPOA, MarinePOA, Aciduptake).   Then for other specialty simulations with long-lived species (CO2, CH4, maybe Hg, etc.), we could leave those as <tt>PRECISION=8</tt>.</blockquote>

'''''[[User:Sebastian D. Eastham|Seb Eastham]] replied:'''''

<blockquote>I would caution that there are still some species relevant to chemistry simulations which are long-lived and could be affected by mass balance – I’m thinking CFCs and N2O in particular. As long as users interested in these species are aware that there might be a random mass drift of ~0.01% per year, that should be fine.</blockquote>

--[[User:Bmy|Bob Yantosca]] ([[User talk:Bmy|talk]]) 20:43, 24 February 2017 (UTC)

=== Mass conservation implications for GCHP ===

[[User:Sebastian D. Eastham|Seb Eastham]] plotted the evolution of the total tracer mass versus time in [[GEOS-Chem HP|GCHP]] for several configurations, as shown below.

*The <span style="color:blue">BLUE</span> line represents a GCHP simulation using <tt>PRECISION=4</tt>.
*The <span style="color:orange">LIGHT ORANGE</span> line represents a GCHP simulation using the default <tt>PRECISION=8</tt>.

[[Image:Transport_error_v3_gchp.png]]

As you can see, the <span style="color:blue">GCHP simulation using the <tt>PRECISION=4</tt> option</span> causes a significant loss of mass, similar to the GEOS-Chem classic simulations from [[#Mass conservation|the preceding section]].  On the other hand, the <span style="color:orange">simulation using <tt>PRECISION=8</tt></span> results in nearly perfect mass conservation.

Seb writes:

<blockquote>A fairly major structural change is needed to ensure mass conservation (specifically, the [GCHP] internal state must be changed from <code>REAL*4</code> to <code>REAL*8</code>). This fix also invalidates "strong reproducibility"; there will be machine-precision-level differences in output when the code is run with different numbers of cores. This basically comes down to a choice; we can have two of the following:

#Mass conservation
#Strong reproducibility
#Speed

I thought 1 and 3 were most important but I figure that it's useful for us all to be aware.</blockquote>

--[[User:Bmy|Bob Yantosca]] ([[User talk:Bmy|talk]]) 15:17, 24 February 2017 (UTC)

== Previous issues that are now resolved ==

=== Update KIND parameters to facilitate interface with the Beijing Climate Center model ===

<span style="color:green">'''''This update was included in [[GEOS-Chem v11-01 benchmark history#v11-01k|v11-01k]] and approved on 19 Dec 2016'''''</span>

'''''Mike Long wrote:'''''

<blockquote>With the flexible precision implementation, GEOS-Chem is unable to compile when the Intel Fortran Compiler flag <tt>-r8</tt> is thrown in the compiler command. The <tt>-r8</tt> flag forces variables declared as <tt>REAL</tt> or <tt>KIND(0.0)</tt> to be double precision by default. I learned this because the Beijing Climate Model uses the <tt>-r8</tt> flag by default, and cannot compile without it.

The result is a multitude of errors that seem intractable. The solution I found was simple, but it may not work altogether -- i.e. it's implications haven't been tested yet. By setting the "f4" definition in precision_mod and "sp" definition in hco_error_mod from <tt>KIND( 0.0 )</tt> to <tt>KIND( REAL( 0.0, 4 ) )</tt> all variables declared as either f4 or sp are allowed to be single precision.</blockquote>

'''''[[User:Bmy|Bob Yantosca]] replied:'''''

<blockquote>I tested Mike's proposed fix on the [[Intel Fortran Compiler|ifort]], [[PGI Fortran compiler|pgfortran]], and [[GNU Fortran compiler|gfortran]] compilers and confirmed that it is the proper one.  We should use:</blockquote>

     KIND( REAL( <span style="color:green">0.0, 4</span> ) ) instead of KIND( <span style="color:red">0.0</span>   ) 
     KIND( REAL( <span style="color:green">0.0, 8</span> ) ) instead of KIND( <span style="color:red">0.0d0</span> ) 

<blockquote>in both <tt>Headers/precision_mod.F</tt> and <tt>HEMCO/Core/hco_error_mod.F90</tt>.  This will prevent variables that you want to keep as 4-byte reals (e.g. for netCDF I/O) from being promoted to 8-byte reals.

Also note, for gfortran we need to add the code in <span style="color:green">GREEN</span> to this variable setting in the <tt>Makefile_header.mk</tt>:</blockquote>

     R8                 := -fdefault-real-8 <span style="color:green">-fdouble-real-8</span>

<blockquote>This will prevent the default double precision size (which in GEOS-Chem should be 8 bytes) from inadvertently being elevated to 16 bytes.</blockquote>

--[[User:Bmy|Bob Yantosca]] ([[User talk:Bmy|talk]]) 17:33, 23 November 2016 (UTC)