add unit test functionality using the Fortan package manager

INSTALL.md

@@ -0,0 +1,43 @@
+# How to build 
+
+Mocsy has several build systems. 
+
+## CMake (preferred, default) 
+
+Ensure you have a functioning Fortran compiler and a CMake of version at least 3.22, then simply do:
+```
+mkdir build
+cd build
+cmake ..
+make -j 
+```
+
+To specifiy where to install mocsy use: `cmake -DCMAKE_INSTALL_PREFIX=$MOCSY_PREFIX ../` and use `make -j install`
+
+## Make 
+
+Simply ensure you have a Fortran compiler and use `make libmocsy.a` 
+
+
+## FPM 
+
+Download and install the [Fortran Package Manager](https://fpm.fortran-lang.org/). You can install the latest version by doing:
+
+```
+git clone https://github.com/fortran-lang/fpm
+cd fpm
+./install.sh
+```
+
+This will put the `fpm` executable in your `~/.local/bin` if you don't have this in your `$PATH` simply do: `export PATH=$PATH:$HOME/.local/bin`
+
+To test the mocsy build use: `fpm test --profile release`
+
+To install mocsy to a defined location use: `fpm install --prefix $MOCSY_PREFIX --profile release`
+
+*NOTE*: The FPM install only supports using double precision for mocsy. 
+
+
+
+
+

fpm.toml

@@ -0,0 +1,38 @@
+name = "mocsy"
+version = "0.1.0"
+license = "MIT"
+author = "ACCESS-NRI"
+maintainer = "jorgegalvez1694@gmail.com"
+copyright = "Copyright 2024, MOCSY?"
+[preprocess]
+[preprocess.cpp]
+suffixes = ["F90", "f90"]
+macros = ["USE_PRECISION=2"]
+
+[library]
+source-dir = "src"
+
+[build]
+auto-executables = false 
+auto-tests = false
+auto-examples = false 
+module-naming = false
+[[test]]
+name = "mocsy_tests"
+source-dir = "test-suite"
+main = "main_tests.f90"
+[install]
+library = true
+
+[fortran]
+implicit-typing = true
+implicit-external = true
+source-form = "free"
+
+[dependencies]
+test-drive.git = "https://github.com/fortran-lang/test-drive"
+[extra.fortitude.check]
+select = ["C", "E", "S"]
+ignore = ["C003"]
+line-length = 132
+per-file-ignores = { "pic_blas_interfaces.F90" = ["C071"] }

src/mocsy_constants.F90

@@ -86,12 +86,12 @@ SUBROUTINE constants(K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa,  &
   !     K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa, K1p, K2p, K3p, Ksi
   !     St, Ft, Bt
 
+
 #if USE_PRECISION == 2
 #   define SGLE(x)    (x)
 #else
 #   define SGLE(x)    REAL(x)
 #endif
-
 ! Input variables
   !>     number of records
 !f2py intent(hide) n
@@ -288,7 +288,7 @@ SUBROUTINE constants(K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa,  &
 !          (see Dickson et al., Best Practices Guide, 2007, Chap. 5, p. 7, including footnote)
         tempot68 = (tempot - 0.0002_rx) / 0.99975_rx
 !       b) Compute "in-situ Temperature" from "Potential Temperature" (both on IPTS 68)
-        tempis68 = sw_temp(sal(i), tempot68, p, SGLE(0.0D0) )
+        tempis68 = sw_temp(sal(i), tempot68, p, SGLE(0.0_rx) )
 !       c) Convert the in-situ temp on older IPTS 68 scale to modern scale (ITS 90)
         tempis = 0.99975_rx*tempis68 + 0.0002_rx
 !       Note: parts (a) and (c) above are tiny corrections;
@@ -435,7 +435,6 @@ SUBROUTINE constants(K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa,  &
                 (148.0248d0 + 137.1942d0*sqrts + 1.62142d0*s) +   &
                 (-24.4344d0 - 25.085d0*sqrts - 0.2474d0*s) *      &
                 dlogtk + 0.053105d0*sqrts*tk)
-
 !       K1p = [H][H2PO4]/[H3PO4]
 !       (seawater scale)
 !       DOE(1994) eq 7.2.20 with footnote using data from Millero (1974)
@@ -496,8 +495,8 @@ SUBROUTINE constants(K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa,  &
 !       Dickson (1990, J. chem. Thermodynamics 22, 113)
         Ks_0p = EXP(-4276.1d0*invtk + 141.328d0 - 23.093d0*dlogtk          &
                 + (-13856.d0*invtk + 324.57d0 - 47.986d0*dlogtk) * sqrtis  &
-                + (35474.d0*invtk - 771.54 + 114.723d0*dlogtk) * is      &
-                - 2698.d0*invtk*is**1.5 + 1776.d0*invtk*is2              &
+                + (35474.d0*invtk - 771.54d0 + 114.723d0*dlogtk) * is      &
+                - 2698.d0*invtk*is**1.5d0 + 1776.d0*invtk*is2              &
                 + LOG(1.0d0 - 0.001005d0*s))
 
 !       Kf = [H][F]/[HF]
@@ -516,6 +515,7 @@ SUBROUTINE constants(K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa,  &
            STOP
         ENDIF
 
+
 !       Kspc (calcite) - apparent solubility product of calcite
 !       (no scale)
 !       Kspc = [Ca2+] [CO32-] when soln is in equilibrium w/ calcite
@@ -547,7 +547,7 @@ SUBROUTINE constants(K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa,  &
            deltav(ipc)  =  a0(ipc) + a1(ipc) *t + a2(ipc) *t*t
            deltak(ipc)   = (b0(ipc)  + b1(ipc) *t + b2(ipc) *t*t)
            lnkpok0(ipc)  = (-(deltav(ipc)) &
-                +(0.5d0*deltak(ipc) * prb) &
+                +(0.5_r8*deltak(ipc) * prb) &
                 )                         * prb/(R*tk)
         END DO
 

src/mocsy_derivauto.F90

@@ -339,8 +339,9 @@ SUBROUTINE derivauto(ph_deriv, pco2_deriv, fco2_deriv, co2_deriv, hco3_deriv, co
              .OR.  dic(i) > 1e+3 &
              .OR.  sil(i) > 1e+3 &
              .OR. phos(i) > 1e+3) THEN
-           PRINT *, 'i, icount, tempot, sal,    alk,    dic,    sil,    phos =', &
-                     i, icount, tempot, sal(i), alk(i), dic(i), sil(i), phos(i)
+            print *, " WARNING, unreasonable input variables !"
+            PRINT *, 'i, icount, tempot, sal,    alk,    dic,    sil,    phos =', &
+                      i, icount, tempot, sal(i), alk(i), dic(i), sil(i), phos(i)
         ENDIF
 !       Zero out any negative salinity, phosphate, silica, dic, and alk
         IF (sal(i) < 0.0) THEN

src/mocsy_derivnum.F90

@@ -349,6 +349,7 @@ SUBROUTINE derivnum (dh_dx, dpco2_dx, dfco2_dx, dco2_dx, dhco3_dx,
     adepth(1) = depth(i)
     alat(1)   = lat(i)
 
+
     IF (deriv_K) THEN
         ! Choose value of absolute perturbation
         abs_delta = K_values(var_index) * 1.d-3   ! 0.1 percent of Kx value
@@ -400,6 +401,11 @@ SUBROUTINE derivnum (dh_dx, dpco2_dx, dfco2_dx, dco2_dx, dhco3_dx,
 
         ! Determine two slightly different values of selected input value
         abs_delta = input_value * rel_delta_x
+        if (abs_delta < 1.0e-12_rx) then
+          ! this prevents a division by 0!!!
+          print *, "Warning: Input value too small for reliable numerical derivative, cycling"
+          cycle
+        end if
         ainput1(1) = input_value - abs_delta
         ainput2(1) = input_value + abs_delta
         ! Compute total absolue delta
@@ -483,7 +489,6 @@ SUBROUTINE derivnum (dh_dx, dpco2_dx, dfco2_dx, dco2_dx, dhco3_dx,
     h(1,2) = 10**(-ph(1,2))
 
     ! Compute derivatives by method of centered difference
-
     dh_dx(i)      = (h(1,2)    - h(1,1))    / dx
     dpco2_dx(i)   = (pco2(1,2) - pco2(1,1)) / dx
     dfco2_dx(i)   = (fco2(1,2) - fco2(1,1)) / dx

src/mocsy_helpers.f90

@@ -0,0 +1,33 @@
+!! General helpers
+module mocsy_test_helpers
+  !! simple reusable helpers for random things
+   use mocsy_singledouble, only: rx
+   implicit none
+
+   private
+
+   real(rx), parameter :: tol_dp = 1.0e-12_rx
+   public :: is_equal
+
+   interface is_equal
+      module procedure is_equal_dp
+   end interface is_equal
+
+contains
+
+   elemental function is_equal_dp(a, b, tol) result(res)
+      real(rx), intent(in) :: a, b
+      real(rx), intent(in), optional :: tol
+      real(rx) :: working_tol
+      logical :: res
+
+      if(present(tol)) then 
+        working_tol = tol
+      else 
+        working_tol = tol_dp
+      end if
+
+      res = abs(a - b) < working_tol
+   end function is_equal_dp
+
+end module mocsy_test_helpers

src/mocsy_vars.F90

@@ -346,8 +346,10 @@ SUBROUTINE vars(ph, pco2, fco2, co2, hco3, co3, OmegaA, OmegaC, BetaD, rhoSW, p,
                .OR.  dic(i) > 1e+3_rx &
                .OR.  sil(i) > 1e+3_rx &
                .OR. phos(i) > 1e+3_rx) THEN
-             PRINT *, 'i, icount, tempot, sal,    alk,    dic,    sil,    phos =', &
-                       i, icount, tempot, sal(i), alk(i), dic(i), sil(i), phos(i)
+
+            print *, " WARNING, unreasonable input variables !"
+              PRINT *, 'i, icount, tempot, sal,    alk,    dic,    sil,    phos =', &
+                        i, icount, tempot, sal(i), alk(i), dic(i), sil(i), phos(i)
           ENDIF
         ENDIF
 !       Zero out any negative salinity, phosphate, silica, dic, and alk
@@ -706,9 +708,9 @@ SUBROUTINE vars_pertK(ph, pco2, fco2, co2, hco3, co3, OmegaA, OmegaC,     &
 !       This is the case for most models and some data
 !       a) Convert the pot. temp on today's "ITS 90" scale to older IPTS 68 scale
 !          (see Dickson et al., Best Practices Guide, 2007, Chap. 5, p. 7, including footnote)
-        tempot68 = (tempot - 0.0002) / 0.99975
+        tempot68 = (tempot - 0.0002_rx) / 0.99975_rx
 !       b) Compute "in-situ Temperature" from "Potential Temperature" (both on IPTS 68)
-        tempis68 = sw_temp(sal(i), SGLE(tempot68), p, SGLE(0.d0) )
+        tempis68 = sw_temp(sal(i), SGLE(tempot68), p, SGLE(0.0_rx) )
 !       c) Convert the in-situ temp on older IPTS 68 scale to modern scale (ITS 90)
         tempis90 = 0.99975*tempis68 + 0.0002
 !       Note: parts (a) and (c) above are tiny corrections;
@@ -740,8 +742,9 @@ SUBROUTINE vars_pertK(ph, pco2, fco2, co2, hco3, co3, OmegaA, OmegaC,     &
              .OR.  dic(i) > 1e+3 &
              .OR.  sil(i) > 1e+3 &
              .OR. phos(i) > 1e+3) THEN
-           PRINT *, 'i, icount, tempot, sal,    alk,    dic,    sil,    phos =', &
-                     i, icount, tempot, sal(i), alk(i), dic(i), sil(i), phos(i)
+            print *, " WARNING, unreasonable input variables !"
+            PRINT *, 'i, icount, tempot, sal,    alk,    dic,    sil,    phos =', &
+                      i, icount, tempot, sal(i), alk(i), dic(i), sil(i), phos(i)
         ENDIF
 !       Zero out any negative salinity, phosphate, silica, dic, and alk
         IF (sal(i) < 0.0) THEN

test-suite/main_tests.f90

@@ -0,0 +1,56 @@
+program pic_tester
+   use, intrinsic :: iso_fortran_env, only: error_unit
+   use testdrive, only: run_testsuite, new_testsuite, testsuite_type, &
+   & select_suite, run_selected, get_argument
+   use test_mocsy, only: collect_mocsy_suite
+   implicit none
+   integer :: stat, is
+   !integer(default_int), parameter :: ntest_suites = 9
+    !! number of tests, this number needs to be modified and equal to the number of files we have with unit tests
+   character(len=:), allocatable :: suite_name, test_name
+   type(testsuite_type), allocatable :: testsuites(:)
+   character(len=*), parameter :: fmt = '("#", *(1x, a))'
+
+   stat = 0
+   !allocate (testsuites(ntest_suites))
+   testsuites = [ &
+                new_testsuite("mocsy", collect_mocsy_suite) &
+                ]
+   ! here you add another test suite to the array
+
+   call get_argument(1, suite_name)
+   call get_argument(2, test_name)
+
+   if (allocated(suite_name)) then
+      is = select_suite(testsuites, suite_name)
+      if (is > 0 .and. is <= size(testsuites)) then
+         if (allocated(test_name)) then
+            write (error_unit, fmt) "Suite:", testsuites(is)%name
+            call run_selected(testsuites(is)%collect, test_name, error_unit, stat)
+            if (stat < 0) then
+               error stop 1
+            end if
+         else
+            write (error_unit, fmt) "Testing:", testsuites(is)%name
+            call run_testsuite(testsuites(is)%collect, error_unit, stat)
+         end if
+      else
+         write (error_unit, fmt) "Available testsuites"
+         do is = 1, size(testsuites)
+            write (error_unit, fmt) "-", testsuites(is)%name
+         end do
+         error stop 1
+      end if
+   else
+      do is = 1, size(testsuites)
+         write (error_unit, fmt) "Testing all:", testsuites(is)%name
+         call run_testsuite(testsuites(is)%collect, error_unit, stat)
+      end do
+   end if
+
+   if (stat > 0) then
+      write (error_unit, "(i0, 1x, a)") stat, "test(s) failed!"
+      error stop 1
+   end if
+
+end program pic_tester