Add CMake build system and Pkgconf pc file

CMakeLists.txt

@@ -0,0 +1,111 @@
+# https://cmake.org/cmake/help/v3.6/module/FindPkgConfig.html
+# 3.11: add_library() doesn't require dummy source
+# 3.12: target_link_libraries() supports OBJECT targets
+# 3.20: cmake_path()
+cmake_minimum_required(VERSION 3.20)
+project(libgtracers
+  VERSION 2024.08.1
+  DESCRIPTION "libgtracers contains a collection of tracers and associated code for use with the MOM and GOLD ocean models."
+  LANGUAGES Fortran
+)
+
+include(GNUInstallDirs)
+
+if(NOT CMAKE_BUILD_TYPE MATCHES "^(Debug|Release|RelWithDebInfo|MinSizeRel)$")
+  message(STATUS "Setting build type to 'Relwithdebinfo' as none was specified.")
+  set(CMAKE_BUILD_TYPE
+      "Relwithdebinfo"
+      CACHE STRING "Choose the type of build." FORCE)
+  set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "Debug" "Release" "MinSizeRel" "RelWithDebInfo")
+endif()
+
+if(NOT CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel|GNU)$")
+  message(
+    WARNING "Compiler not officially supported: ${CMAKE_Fortran_COMPILER_ID}")
+endif()
+
+cmake_path(APPEND CMAKE_Fortran_MODULE_DIRECTORY ${CMAKE_BINARY_DIR} "include")
+
+################################################################################
+# Fortran
+################################################################################
+
+if(CMAKE_Fortran_COMPILER_ID MATCHES "GNU")
+
+  # Copied from MOM5/bin/mkmf.template.nci.gfortran
+  set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fcray-pointer -fdefault-real-8 -ffree-line-length-none -fno-range-check -Waliasing -Wampersand -Warray-bounds -Wcharacter-truncation -Wconversion -Wline-truncation -Wintrinsics-std -Wsurprising -Wno-tabs -Wunderflow -Wunused-parameter -Wintrinsic-shadow -Wno-align-commons")
+  set(CMAKE_Fortran_FLAGS_RELEASE "-O2")
+  set(CMAKE_Fortran_FLAGS_RELWITHDEBINFO "-g")
+  set(CMAKE_Fortran_FLAGS_DEBUG "-O0 -g -W -fbounds-check")
+
+elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "Intel")
+
+  # Copied from MOM5/bin/mkmf.template.nci
+  set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fno-alias -safe-cray-ptr -fpe0 -ftz -assume byterecl -i4 -r8 -traceback -nowarn -check noarg_temp_created -assume nobuffered_io -convert big_endian -grecord-gcc-switches -align all")
+  set(CMAKE_Fortran_FLAGS_RELEASE "-g3 -O2 -xCORE-AVX2 -debug all -check none")
+  set(CMAKE_Fortran_FLAGS_RELWITHDEBINFO "-g3 -O2 -xCORE-AVX2 -debug all -check none")
+  set(CMAKE_Fortran_FLAGS_DEBUG "-g3 -O0 -debug all -check -check noarg_temp_created -check nopointer -warn -warn noerrors -ftrapuv")
+
+elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "IntelLLVM")
+
+  # Copied from MOM5/bin/mkmf.template.nci
+  set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fno-alias -safe-cray-ptr -fpe0 -ftz -assume byterecl -i4 -r8 -traceback -nowarn -check noarg_temp_created -assume nobuffered_io -convert big_endian -grecord-gcc-switches -align all")
+  set(CMAKE_Fortran_FLAGS_RELEASE "-g3 -O2 -xCORE-AVX2 -debug all -check none")
+  set(CMAKE_Fortran_FLAGS_RELWITHDEBINFO "-g3 -O2 -xCORE-AVX2 -debug all -check none")
+  set(CMAKE_Fortran_FLAGS_DEBUG "-g3 -O0 -debug all -check -check noarg_temp_created -check nopointer -warn -warn noerrors -ftrapuv")
+
+else()
+  message(WARNING "Fortran compiler with ID ${CMAKE_Fortran_COMPILER_ID} will be used with CMake default options")
+endif()
+
+################################################################################
+
+find_package(MPI REQUIRED COMPONENTS Fortran)
+find_package(FMS COMPONENTS R8 REQUIRED)
+find_package(NetCDF REQUIRED COMPONENTS Fortran)
+find_package(PkgConfig REQUIRED)
+pkg_check_modules(MOCSY REQUIRED IMPORTED_TARGET "mocsy")
+
+add_library(gtracers STATIC)
+
+target_sources(gtracers PRIVATE
+  generic_tracers/FMS_coupler_util.F90
+  generic_tracers/FMS_ocmip2_co2calc.F90
+  generic_tracers/generic_abiotic.F90
+  generic_tracers/generic_age.F90
+  generic_tracers/generic_argon.F90
+  generic_tracers/generic_BLING.F90
+  generic_tracers/generic_blres.F90
+  generic_tracers/generic_CFC.F90
+  generic_tracers/generic_COBALT.F90
+  generic_tracers/generic_ERGOM.F90
+  generic_tracers/generic_miniBLING.F90
+  generic_tracers/generic_SF6.F90
+  generic_tracers/generic_TOPAZ.F90
+  generic_tracers/generic_tracer.F90
+  generic_tracers/generic_tracer_utils.F90
+  generic_tracers/generic_WOMBATlite.F90
+)
+
+target_link_libraries(gtracers PUBLIC
+			FMS::fms_r8
+			NetCDF::NetCDF_Fortran
+			PkgConfig::MOCSY)
+
+install(TARGETS gtracers
+        EXPORT GFDLGTracersTargets)
+
+install(DIRECTORY ${CMAKE_Fortran_MODULE_DIRECTORY}/ TYPE INCLUDE)
+
+install(EXPORT GFDLGTracersTargets
+  FILE GFDLGTracersTargets.cmake
+  NAMESPACE GFDLGTracers::
+  DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/GFDLGTracers
+)
+
+configure_file(${CMAKE_CURRENT_SOURCE_DIR}/gfdl-generic-tracers.pc.in
+               ${CMAKE_CURRENT_BINARY_DIR}/gfdl-generic-tracers.pc @ONLY)
+
+install(FILES ${CMAKE_CURRENT_BINARY_DIR}/gfdl-generic-tracers.pc
+              DESTINATION ${CMAKE_INSTALL_LIBDIR}/pkgconfig
+              COMPONENT utilities)

generic_tracers/FMS_ocmip2_co2calc.F90

@@ -209,7 +209,7 @@ subroutine FMS_ocmip2_co2calc(dope_vec, mask,                      &
 !       local variables
 !
 integer :: isc, iec, jsc, jec
-integer :: i,j
+integer :: i,j,ipc
 real :: alpha_internal
 real :: bt
 real :: co2star_internal
@@ -229,17 +229,32 @@ subroutine FMS_ocmip2_co2calc(dope_vec, mask,                      &
 real :: ks
 real :: ksi
 real :: kw
+real :: Kspa
+real :: Kspc
+real :: calcium
 real :: log100
 real :: s2
 real :: scl
 real :: sqrtis
 real :: sqrts
+real :: s15
 real :: st
 real :: tk
 real :: tk100
 real :: tk1002
 real :: logf_of_s
+real :: prb
 real :: salinity
+real :: R
+real :: total2free_0p, total2free, ks_0p, kf_0p, free2SWS, total2SWS, SWS2total, free2SWS_0p, total2SWS_0p
+real, dimension(12) :: deltav, deltak, lnkpok0
+real, dimension(12) :: a0, a1, a2, b0, b1, b2
+DATA a0 /-25.5,    -15.82,      -29.48,    -20.02,     -18.03,   -9.78,     -48.76,    -45.96,    -14.51,    -23.12,    -26.57,    -29.48/
+DATA a1 /0.1271,   -0.0219,     0.1622,    0.1119,     0.0466,   -0.0090,   0.5304,    0.5304,    0.1211,    0.1758,    0.2020,    0.1622/
+DATA a2 /0.0,      0.0,         -2.608e-3, -1.409e-3,  0.316e-3, -0.942e-3, 0.0,       0.0,       -0.321e-3, -2.647e-3, -3.042e-3, -2.6080e-3/
+DATA b0 /-3.08e-3, 1.13e-3,     -2.84e-3,  -5.13e-3,   -4.53e-3, -3.91e-3,  -11.76e-3, -11.76e-3, -2.67e-3,  -5.15e-3,  -4.08e-3,  -2.84e-3/
+DATA b1 /0.0877e-3, -0.1475e-3, 0.0,        0.0794e-3, 0.09e-3,  0.054e-3,  0.3692e-3, 0.3692e-3, 0.0427e-3, 0.09e-3,   0.0714e-3, 0.0/
+DATA b2 /12*0.0/
 
 character(len=10) :: co2_calc_method
 
@@ -254,8 +269,7 @@ subroutine FMS_ocmip2_co2calc(dope_vec, mask,                      &
   co2_calc_method = 'ocmip2'
 end if
 
-
-
+R = 83.14472
 
 ! Set the loop indices.
   isc = dope_vec%isc ; iec = dope_vec%iec
@@ -364,8 +378,10 @@ subroutine FMS_ocmip2_co2calc(dope_vec, mask,                      &
         sqrtis    = sqrt(max(0.0,is))
         s2        = s_in(i,j) * s_in(i,j)
         sqrts     = sqrt(max(0.0,s_in(i,j)))
+        s15       = sqrts**3
         scl       = s_in(i,j) / 1.80655
         logf_of_s = log(1.0 - 0.001005 * s_in(i,j))
+        prb       = zt(i,j) / 10.0
 !
 ! k0 from Weiss 1974
 !
@@ -478,6 +494,85 @@ subroutine FMS_ocmip2_co2calc(dope_vec, mask,                      &
 !
         ft = 0.000067 / 18.9984 * scl
 !
+! Omega saturation states
+!
+        if (present(omega_arag)) then
+          Kspa = 10.0**(-171.945 - 0.077993*tk + 2903.293/tk + 71.595*LOG10(tk) &
+                      + (-0.068393 + 0.0017276*tk + 88.135/tk)*sqrts &
+                      - 0.10018*s_in(i,j) + 0.0059415*s15 )
+        endif
+        if (present(omega_calc)) then
+          Kspc = 10.0**(-171.9065 - 0.077993*tk + 2839.319/tk + 71.595*LOG10(tk) &
+                      + (-0.77712 + 0.0028426*tk + 178.34/tk)*sqrts &
+                      - 0.07711*s_in(i,j) + 0.0041249*s15 )
+        endif
+
+      !!! ------------------------------------------- !!! 
+      !!! Pressure corrections to the above constants !!! (Pearse J. Buchanan, July 2024)
+      !!! ------------------------------------------- !!!      (copied from mocsy)
+      ! This must be done to ensure that the co2 sys    !
+      ! calculations are accurate beneath the surface   !
+
+      do ipc=1,12
+        deltav(ipc)  = a0(ipc) + a1(ipc)*t_in(i,j) + a2(ipc)*t_in(i,j)*t_in(i,j)
+        deltak(ipc)  = b0(ipc) + b1(ipc)*t_in(i,j) + b2(ipc)*t_in(i,j)*t_in(i,j)
+        lnkpok0(ipc) = ( -deltav(ipc) + (0.5*deltak(ipc) * prb) )*prb/(R*tk)
+      enddo
+
+      ! Conversion factor total -> free scale at pressure zero
+      total2free_0p  = 1.0/(1.0 + st/ks)   ! Kfree = Ktotal*total2free
+      ks_0p = ks*1
+      ks = ks * EXP(lnkpok0(5))
+      ! Conversion factor total -> free scale
+      total2free     = 1.0/(1.0 + st/ks)   ! Kfree = Ktotal*total2free
+
+      kf_0p = kf * total2free_0p
+      kf    = kf_0p * EXP(lnkpok0(6))
+      kf    = kf / total2free
+
+      ! Convert between seawater and total hydrogen (pH) scales
+      free2SWS  = 1.0 + st/ks + ft/(kf*total2free)  ! using Kf on free scale
+      total2SWS = total2free * free2SWS             ! KSWS = Ktotal*total2SWS
+      SWS2total = 1.0 / total2SWS
+      ! Conversion at pressure zero
+      free2SWS_0p  = 1.0 + st/ks_0p + ft/(kf_0p)    ! using Kf on free scale
+      total2SWS_0p = total2free_0p * free2SWS_0p    ! KSWS = Ktotal*total2SWS
+
+      ! Convert from Total to Seawater scale before pressure correction
+      ! Must change to SEAWATER scale: Kb
+      kb = kb*total2SWS_0p
+
+      ! Already on SEAWATER scale: K1p, K2p, K3p, Kb, Ksi, Kw
+
+      ! Other contants (keep on another scale):
+      !    - K0         (independent of pH scale, already pressure corrected)
+      !    - Ks         (already on Free scale;   already pressure corrected)
+      !    - Kf         (already on Total scale;  already pressure corrected)
+      !    - Kspc, Kspa (independent of pH scale; pressure-corrected below)
+
+      ! Perform actual pressure correction (on seawater scale)
+      k1 = k1*EXP(lnkpok0(1))
+      k2 = k2*EXP(lnkpok0(2))
+      kb = kb*EXP(lnkpok0(3))
+      kw = kw*EXP(lnkpok0(4))
+      Kspc = Kspc*EXP(lnkpok0(7))
+      Kspa = Kspa*EXP(lnkpok0(8))
+      k1p = k1p*EXP(lnkpok0(9))
+      k2p = k2p*EXP(lnkpok0(10))
+      k3p = k3p*EXP(lnkpok0(11))
+      ksi = ksi*EXP(lnkpok0(12))
+
+      ! Convert back to original total scale:
+      k1  = k1 *SWS2total
+      k2  = k2 *SWS2total
+      k1p = k1p*SWS2total
+      k2p = k2p*SWS2total
+      k3p = k3p*SWS2total
+      kb  = kb *SWS2total
+      ksi = ksi*SWS2total
+      kw  = kw *SWS2total
+
+!
 !***********************************************************************
 !
 ! Calculate [H+] total when DIC and TA are known at T, S and 1 atm.
@@ -515,6 +610,12 @@ subroutine FMS_ocmip2_co2calc(dope_vec, mask,                      &
                          k1 * htotal(i,j) + k1 * k2)
         if (present(co2star)) co2star(i,j) = co2star_internal
         if (present(co3_ion)) co3_ion(i,j) = co2star_internal * k1 * k2 / htotal2
+        if (present(omega_arag) .or. present(omega_calc)) then
+          co3_ion(i,j) = co2star_internal * k1 * k2 / htotal2
+          calcium = (0.02128/40.078)*s_in(i,j)/1.80655
+          if (present(omega_arag)) omega_arag(i,j) = (calcium * co3_ion(i,j)) / Kspa
+          if (present(omega_calc)) omega_calc(i,j) = (calcium * co3_ion(i,j)) / Kspc
+        endif
 !
 ! Weiss & Price (1980, Mar. Chem., 8, 347-359; Eq 13 with table 6
 !                values)
@@ -548,7 +649,12 @@ subroutine FMS_ocmip2_co2calc(dope_vec, mask,                      &
       if (present(pco2surf)) then  !{
         pCO2surf(i,j) = 0.0
       endif  !}
-
+      if (present(omega_arag)) then
+        omega_arag(i,j) = 0.0
+      endif
+      if (present(omega_calc)) then
+        omega_calc(i,j) = 0.0
+      endif
 
     endif  !}mask