MOM6 makes use of Fortran2003 and later extensions but only when supported by all available compilers. We try to avoid the use of very modern Fortran constructs that can limit portability. To help keep the code readily understandable, this page makes recommendations about how to use or not use the various features of the modern Fortran language.
Code style is typically a personal choice, but when styles clash it can lead to discord. These standards have been adopted in an attempt to promote harmony and clarity.
if test with no then, which could be 2 space indent as though the then and endif were present. I.e.,
if (test) &
var = expression
and
if (test) then
var = expression
endif
should have the same indenting on var = expression so that there would be no need to change the indenting if another statement were added to the if block, or if the latter expression were converted to the former.
if(a==0) is legal fortran but bad style. Use if (a==0) instead.if (a == 0) is even better, since == is a language token.a = b is strongly preferred over a=b
do i=is,ie is acceptablecall fn(arg1, arg_name=a) is strongly preferred over call fn(arg1, arg_name = a)call fn(arg1, arg2, arg3)), but usually not after the comma separating array indices for multidimensional arrays (e.g., A(i,j,k) = B(i,j,k)).Some compilers handle very long lines of code gracefully, but MOM6 needs to adhere to the Fortran standard, which is 132 characters for code, after any macro expansion. MOM6 does use macros for some memory declarations, so we need to build in some added space in setting MOM6 guidelines:
del_rho_int (delta rho at interface).delta_rho).
do and if constructs should be terminated with the combined enddo and endif statements, respectively.end token (e.g. end program [label])
program, module, type, subroutine, function, interface, selecti,j,k are used for cell-center, layer-center references, e.g. h(i,j,k), T(i+1,j,k).I,J are used for staggered, cell-edge references, e.g. u(I,j,k), v(i,J,k), q(I,J,k), u(I-1,j,k). We use a north-east staggering convention so the I means i+1/2 and I-1 means i-1/2.K is used for the interface (between layer) references, e.g. del_t(i,j,K) = T(i,j,K+1) - T(i,j,K). The vertical staggering is such that interface K=1 is above layer k=1 so that K means k-1/2 and K+1 means k+1/2., only modifierimplicit none ; private
implicit nonetv%S = 0. is forbidden.tv%S(:,:,:) = 0.S_tmp(:,:,:) = tv%S(:,:,:)tmp(:,:) = 1.0 / G%areaT(:,:) might have zeros in the halo region.call post_data(id_AT, G%areaT(:,:)*tv%(T(:,:,1)) is wrong because it can use uninitialized data in halos.intent(in), intent(out), intent(inout).Divisions are prone to NaNs and relatively expensive. An optimizing compiler will often rearrange math which makes debugging divisions by zero harder to catch.
Q(i,j) * G%IareaT(i,j) instead of Q(i,j) / G%areaT(i,j).B / C * D which is ambiguous to humans (not the compiler)
( B * D ) / CA / ( A + B / C)
( A * C ) / ( A * C + B)Floating point operations are sensitive to the order of operations (associativity), which can not generally be guaranteed due to compiler serialization and optimization.
Addition operations must be done in pairs. When more than one addition is required, the order should be specified using parentheses.
z = a + b + cz = (a + b) + cIdeally, the order of operation should be chosen to give the best accuracy. For example, if a = 1. b = -1. and c = 1.e-20, then the order should be chosen to preserve the residual value.
a + (b + c) == 0.(a + b) + c == 1.e-20Not only does this impact reproducibility, but the second choice is more accurate and avoids a potential division by zero.
All operations should be ordered, but no particular ordering is enforced. Contributors are encouraged to consider the most accurate order of operations. In some cases the order of sums can be chosen to give expressions that yield identical answers if the underlying horizontal coordinate is rotated by 90 or 180 degrees, which would define the preferred order of operations.
sum() intrinsicWe avoid the Fortran sum() intrinsic since the result is dependent on the order of operations within the summation. Using explicit loops allows us to define the order of summation. So
a = sum(b(:))
should be
a = 0.
do k = 1, nz
a = a + b(k)
enddo
The prod() and matmul() intrinsics should also not be used.
Floating point operations across MPI ranks are volatile, since the order can change depending on the state of the network. Functions such as MPI_Reduce will not generally be reproducible when used for floating point arithmetic.
When performing summations over MPI ranks, use the reproducing_sum function.
use MOM_coms, only: reproducing_sum
...
sum = reproducing_sum(array(:,:))
Multiplication is also non-associative and thus not reproducible, but the impact is typically small. Results may depend on the order of operations, most often in the least significant bit of the fractional component.
In single precision, if a = b = 1 + 2**-23 and c = 1.5, then the following calculations differ:
(a*b)*c = 1.50000036a*(b*c) = 1.50000048Parentheses in multiplication operations are currently not enforced, but contributors should consider using them when applicable.
Use of transcendental functions, such as trigonometric functions, non-integer powers, and logarithms, are often implementation-dependent and should be avoided when possible.
The exponent operator, a**b should be used sparingly, since compilers will often internally replace it with pow(a, b), which is often computed as a transcendental function, exp(b * log(a)). Even small integral power, such as a**3, have been known to be replaced with pow(a, 3). To maximize reproducibility, integral powers should be explicitly computed, e.g. a3 = a * a * a.
Square roots (a**0.5) should always use the sqrt() intrinsic. An IEEE-754 compliant sqrt function must be exactly rounded.
Cube roots (a**(1./3.)) should be avoided, the MOM6 intrinsic cuberoot should be used. This is not exactly rounded, but it is reproducible.
Every module follows this pattern:
!> Brief module description
module MOM_module_name
! This file is part of MOM6. See LICENSE.md for the license.
use MOM_some_module, only : specific_symbol
use MOM_other, only : other_thing
implicit none ; private
#include <MOM_memory.h>
public :: exported_routine_1, exported_routine_2
!> Control structure for this module
type, public :: module_CS ; private
real :: param !< Description [units]
integer :: id_diag = -1 !< Diagnostic ID for some_field
end type module_CS
contains
!> Initialize the module, read parameters, register diagnostics
subroutine module_init(Time, G, GV, US, param_file, diag, CS)
! ... call log_version, get_param, register_diag_field ...
end subroutine module_init
!> Deallocate module memory
subroutine module_end(CS)
! ... cleanup ...
end subroutine module_end
!> \namespace MOM_module_name
!! Extended description, references, and equations
end module MOM_module_name
Key rules:
_init and _end subroutines for lifecycle management! Local variables comment separates dummy arguments from local declarationsallocatable over pointer for control structure membersMOM_something.F90 contains module MOM_something)module_CS (e.g., energetic_PBL_CS), always public but opaque (with private contents)id_diag_name, initialized to -1I (e.g., IdxCu = 1/dxCu, IareaT = 1/areaT)G (ocean_grid_type), GV (verticalGrid_type), US (unit_scale_type)Array dimensions use preprocessor macros from MOM_memory.h:
SZI_(G), SZJ_(G), SZK_(GV) for explicit-shape cell-center arraysSZIB_(G), SZJB_(G), SZKB_(GV) for explicit-shape face/edge-point arraysNIMEM_, NJMEM_, NKMEM_ for allocatable arraysMOM6 uses a dimensional annotation system for every real variable. Units are documented in square brackets at the end of comments, using a two-part notation:
[rescaled_dimensions ~> MKS_equivalent]
| Symbol | Physical Dimension | MKS Unit |
|---|---|---|
| Z | Vertical depth/distance | m |
| L | Horizontal length | m |
| T | Time | s |
| H | Layer thickness | m (Boussinesq) or kg m-2 |
| R | Density | kg m-3 |
| Q | Enthalpy | J kg-1 |
| C | Temperature | degC |
| S | Salinity | ppt |
| A | Arbitrary/generic units | a |
real :: velocity !< Horizontal velocity [L T-1 ~> m s-1]
real :: pressure !< Hydrostatic pressure [R L2 T-2 ~> Pa]
real :: thickness !< Layer thickness [H ~> m or kg m-2]
real :: diffusivity !< Vertical diffusivity [Z2 T-1 ~> m2 s-1]
real :: slope !< Isopycnal slope [Z L-1 ~> nondim]
real :: efficiency !< Mixing efficiency [nondim]
real :: field !< A field in arbitrary units [A]
real :: Z_to_m !< Scaling factor [m Z-1 ~> 1]
[brackets] at the end of its commentH L2 not L2 H)[H ~> m or kg m-2] when units differ by approximation[T2 Z-1 ~> s2 m-1], not [H Z T-1 / H Z2 T-3 = T2 Z-1 ~> s2 m-1]m-1, s-2, kg-3 (no slashes like 1/m)[nondim][A] or [A ~> a], never [arbitrary][target source-1 ~> 1], e.g., [Z m-1 ~> 1]!> for documentation comments on the following entity!< for inline documentation on the preceding entity (same line)!! for multi-line continuation (no blank lines between)!>@{ and !>@} for grouping related declarations!> header describing purpose!< or !> including units!< including units\f$ ... \f$ (inline) or \f[ ... \f] (display)end module using \namespace