Add skin wrapper+tests

.gitignore

@@ -7,3 +7,6 @@ _version.py
 .eggs/
 build/
 __pycache__
+mod_aerobulk_wrap.*.so
+source/fortran/mod_aerobulk_wrap-f2pywrappers2.f90
+source/fortran/mod_aerobulk_wrapmodule.c

source/aerobulk/__init__.py

@@ -6,4 +6,4 @@
     __version__ = "unknown"
     pass
 
-from .flux import flux_noskin
+from .flux import flux_noskin, flux_skin

source/aerobulk/flux.py

@@ -1,9 +1,143 @@
 import aerobulk.aerobulk.mod_aerobulk_wrap as aero
 
+VALID_ALGOS = ["coare3p0", "coare3p6", "ecmwf", "ncar", "andreas"]
+VALID_ALGOS_SKIN = ["coare3p0", "coare3p6", "ecmwf"]
+
+
+def _check_algo(algo, valids):
+    if algo not in valids:
+        raise ValueError(f"Algorithm {algo} not valid. Choose from {valids}.")
+
 
 def flux_noskin(
     sst, t_zt, hum_zt, u_zu, v_zu, slp=101000.0, algo="coare3p0", zt=10, zu=2, niter=1
 ):
-    return aero.mod_aerobulk_wrapper.aerobulk_model_noskin(
+    """Python wrapper for aerobulk without skin correction.
+    !ATTENTION If input not provided in correct units, will crash.
+
+    Parameters
+    ----------
+    sst : array_like
+        Bulk sea surface temperature [K]
+    t_zt : array_like
+        Absolute air temperature at height zt [K]
+    hum_zt : array_like
+        air humidity at zt, can be given as:
+            - specific humidity                   [kg/kg]
+            - dew-point temperature                   [K]
+            - relative humidity                       [%]
+            => type should normally be recognized based on value range
+    u_zu : array_like
+        zonal wind speed at zu                    [m/s]
+    v_zu : array_like
+        meridional wind speed at zu               [m/s]
+    slp : array_like, optional
+        mean sea-level pressure                   [Pa] ~101000 Pa,
+        by default 101000.0
+    algo : str, optional
+        Algorithm, can be one of: "coare3p0", "coare3p6", "ecmwf", "ncar", "andreas",
+        by default "coare3p0"
+    zt : int, optional
+        height for temperature and spec. hum. of air           [m],
+        by default 10
+    zu : int, optional
+        height for wind (10m = traditional anemometric height  [m],
+        by default 2
+    niter : int, optional
+        Number of iteration steps used in the algorithm,
+        by default 1
+
+    Returns
+    -------
+    ql : array_like
+        Latent heat flux                                     [W/m^2]
+    qh : array_like
+        Sensible heat flux                                   [W/m^2]
+    taux : array_like
+        zonal wind stress                                    [N/m^2]
+    tauy : array_like
+        meridional wind stress                                    [N/m^2]
+    evap : array_like
+        evaporation         [mm/s] aka [kg/m^2/s] (usually <0, as ocean loses water!)
+    """
+    _check_algo(algo, VALID_ALGOS)
+    ql, qh, taux, tauy, evap = aero.mod_aerobulk_wrapper.aerobulk_model_noskin(
         algo, zt, zu, sst, t_zt, hum_zt, u_zu, v_zu, slp, niter
     )
+    return ql, qh, taux, tauy, evap
+
+
+def flux_skin(
+    sst,
+    t_zt,
+    hum_zt,
+    u_zu,
+    v_zu,
+    rad_sw,
+    rad_lw,
+    slp=101000.0,
+    algo="coare3p0",
+    zt=10,
+    zu=2,
+    niter=1,
+):
+    """Python wrapper for aerobulk with skin correction.
+    !ATTENTION If input not provided in correct units, will crash.
+
+    Parameters
+    ----------
+    sst : array_like
+        Bulk sea surface temperature [K]
+    t_zt : array_like
+        Absolute air temperature at height zt [K]
+    hum_zt : array_like
+        air humidity at zt, can be given as:
+            - specific humidity                   [kg/kg]
+            - dew-point temperature                   [K]
+            - relative humidity                       [%]
+            => type should normally be recognized based on value range
+    u_zu : array_like
+        zonal wind speed at zu                    [m/s]
+    v_zu : array_like
+        meridional wind speed at zu               [m/s]
+    rad_sw : array_like
+        downwelling shortwave radiation at the surface (>0)   [W/m^2]
+    rad_lw : array_like
+        rad_lw : downwelling longwave radiation at the surface  (>0)   [W/m^2]
+    slp : array_like, optional
+        mean sea-level pressure                   [Pa] ~101000 Pa,
+        by default 101000.0
+    algo : str, optional
+        Algorithm, can be one of: "coare3p0", "coare3p6", "ecmwf",
+        by default "coare3p0"
+    zt : int, optional
+        height for temperature and spec. hum. of air           [m],
+        by default 10
+    zu : int, optional
+        height for wind (10m = traditional anemometric height  [m],
+        by default 2
+    niter : int, optional
+        Number of iteration steps used in the algorithm,
+        by default 1
+
+    Returns
+    -------
+    ql : array_like
+        Latent heat flux                                     [W/m^2]
+    qh : array_like
+        Sensible heat flux                                   [W/m^2]
+    taux : array_like
+        zonal wind stress                                    [N/m^2]
+    tauy : array_like
+        meridional wind stress                                    [N/m^2]
+    t_s : array_like
+        skin temperature    [K]    (only when l_use_skin_schemes=TRUE)
+    evap : array_like
+        evaporation         [mm/s] aka [kg/m^2/s] (usually <0, as ocean loses water!)
+    """
+    _check_algo(algo, VALID_ALGOS_SKIN)
+
+    ql, qh, taux, tauy, t_s, evap = aero.mod_aerobulk_wrapper.aerobulk_model_skin(
+        algo, zt, zu, sst, t_zt, hum_zt, u_zu, v_zu, slp, rad_sw, rad_lw, niter
+    )
+    return ql, qh, taux, tauy, t_s, evap

source/fortran/mod_aerobulk_wrap.f90

@@ -8,7 +8,7 @@ MODULE mod_aerobulk_wrapper
 
   PRIVATE
 
-  PUBLIC :: AEROBULK_MODEL_NOSKIN
+  PUBLIC :: AEROBULK_MODEL_NOSKIN, AEROBULK_MODEL_SKIN
 
 CONTAINS
 
@@ -45,4 +45,38 @@ SUBROUTINE AEROBULK_MODEL_NOSKIN( Ni, Nj, Nt,          &
 
   END SUBROUTINE AEROBULK_MODEL_NOSKIN
 
+  SUBROUTINE AEROBULK_MODEL_SKIN( Ni, Nj, Nt,          &
+     &                       calgo, zt, zu, sst, t_zt,   &
+     &                       hum_zt, U_zu, V_zu, slp,    &
+     &                       rad_sw, rad_lw,             &
+     &                       Niter,                      &
+     &                       QL, QH, Tau_x, Tau_y, T_s, Evap)
+
+    INTEGER,                  INTENT(in)  :: Ni, Nj, Nt
+    CHARACTER(len=*),         INTENT(in)  :: calgo
+    REAL(8),                 INTENT(in)  :: zt, zu
+    REAL(8), DIMENSION(Ni,Nj,Nt), INTENT(in)  :: sst, t_zt, hum_zt, U_zu, V_zu, slp, rad_sw, rad_lw
+    INTEGER,                  INTENT(in) :: Niter
+    REAL(8), DIMENSION(Ni,Nj,Nt), INTENT(out) :: QL, QH, Tau_x, Tau_y, T_s, Evap
+
+    INTEGER n
+
+    nb_iter = Niter
+
+    !! initialize based on first timestep
+    CALL AEROBULK_INIT(Nt, calgo, sst(:, :, 1), t_zt(:, :, 1), &
+        &              hum_zt(:, :, 1), U_zu(:, :, 1), V_zu(:, :, 1), slp(:, :, 1), &
+        &              l_use_skin=.TRUE., prsw=rad_sw(:, :, 1), prlw=rad_lw(:, :, 1))
+
+    DO n = 1, Nt
+       CALL AEROBULK_COMPUTE(n, calgo, zt, zu, sst(:, :, n), t_zt(:, :, n), &
+          &                  hum_zt(:, :, n), U_zu(:, :, n), V_zu(:, :, n), slp(:, :, n),  &
+          &                  QL(:, :, n), QH(:, :, n), Tau_x(:, :, n), Tau_y(:, :, n),     &
+          &                  rad_sw=rad_sw(:, :, n), rad_lw=rad_lw(:, :, n), T_s=T_s(:, :, n), Evp=Evap(:, :, n))
+    END DO
+
+    CALL AEROBULK_BYE()
+
+  END SUBROUTINE AEROBULK_MODEL_SKIN
+
 END MODULE mod_aerobulk_wrapper

source/fortran/mod_aerobulk_wrap.pyf

@@ -27,6 +27,29 @@ python module mod_aerobulk_wrap ! in
                 real(kind=8) dimension(ni,nj,nt),intent(out),depend(nj,ni,nt) :: tau_y
                 real(kind=8) dimension(ni,nj,nt),intent(out),depend(nj,ni,nt) :: evap
             end subroutine aerobulk_model_noskin
+            subroutine aerobulk_model_skin(ni,nj,nt,calgo,zt,zu,sst,t_zt,hum_zt,u_zu,v_zu,slp,rad_sw,rad_lw,niter,ql,qh,tau_x,tau_y,t_s,evap) ! in :mod_aerobulk_wrap:mod_aerobulk_wrap.f90:mod_aerobulk_wrapper
+                integer, optional,intent(in),check(shape(sst, 0) == ni),depend(sst) :: ni=shape(sst, 0)
+                integer, optional,intent(in),check(shape(sst, 1) == nj),depend(sst) :: nj=shape(sst, 1)
+                integer, optional,intent(in),check(shape(sst, 2) == nt),depend(sst) :: nt=shape(sst, 2)
+                character*(*) intent(in) :: calgo
+                real(kind=8) intent(in) :: zt
+                real(kind=8) intent(in) :: zu
+                real(kind=8) dimension(ni,nj,nt),intent(in) :: sst
+                real(kind=8) dimension(ni,nj,nt),intent(in),depend(nj,ni,nt) :: t_zt
+                real(kind=8) dimension(ni,nj,nt),intent(in),depend(nj,ni,nt) :: hum_zt
+                real(kind=8) dimension(ni,nj,nt),intent(in),depend(nj,ni,nt) :: u_zu
+                real(kind=8) dimension(ni,nj,nt),intent(in),depend(nj,ni,nt) :: v_zu
+                real(kind=8) dimension(ni,nj,nt),intent(in),depend(nj,ni,nt) :: slp
+                real(kind=8) dimension(ni,nj,nt),intent(in),depend(nj,ni,nt) :: rad_sw
+                real(kind=8) dimension(ni,nj,nt),intent(in),depend(nj,ni,nt) :: rad_lw
+                integer intent(in) :: niter
+                real(kind=8) dimension(ni,nj,nt),intent(out),depend(nj,ni,nt) :: ql
+                real(kind=8) dimension(ni,nj,nt),intent(out),depend(nj,ni,nt) :: qh
+                real(kind=8) dimension(ni,nj,nt),intent(out),depend(nj,ni,nt) :: tau_x
+                real(kind=8) dimension(ni,nj,nt),intent(out),depend(nj,ni,nt) :: tau_y
+                real(kind=8) dimension(ni,nj,nt),intent(out),depend(nj,ni,nt) :: t_s
+                real(kind=8) dimension(ni,nj,nt),intent(out),depend(nj,ni,nt) :: evap
+            end subroutine aerobulk_model_skin
         end module mod_aerobulk_wrapper
     end interface
 end python module mod_aerobulk_wrap

test/test.py

@@ -1,8 +1,15 @@
 import numpy as np
-from aerobulk import flux_noskin
+import pytest
+from aerobulk import flux_noskin, flux_skin
+from aerobulk.flux import VALID_ALGOS, VALID_ALGOS_SKIN
 
 
-def test_smoke_test():
+@pytest.mark.parametrize(
+    "algo",
+    VALID_ALGOS
+    + [pytest.param("wrong", marks=pytest.mark.xfail(strict=True, raises=ValueError))],
+)
+def test_smoke_noskin(algo):
     # I am very confused about shape.
     # We want Nt to be the *slowest varying dimension*
     # For C-order arrays, would give (Nt, Ny, Nx)
@@ -17,4 +24,32 @@ def test_smoke_test():
     u_zu = np.full(shape, 1.0, order=order)
     v_zu = np.full(shape, -1.0, order=order)
     slp = np.full(shape, 101000.0, order=order)
-    flux_noskin(sst, t_zt, hum_zt, u_zu, v_zu, slp)
+    flux_noskin(sst, t_zt, hum_zt, u_zu, v_zu, slp, algo=algo)
+
+
+@pytest.mark.parametrize(
+    "algo",
+    VALID_ALGOS_SKIN
+    + [
+        pytest.param("wrong", marks=pytest.mark.xfail(strict=True, raises=ValueError)),
+        pytest.param("ncar", marks=pytest.mark.xfail(strict=True, raises=ValueError)),
+    ],
+)
+def test_smoke_skin(algo):
+    # I am very confused about shape.
+    # We want Nt to be the *slowest varying dimension*
+    # For C-order arrays, would give (Nt, Ny, Nx)
+    # For F-order arrays, would give (Nx, Ny, Nt)
+    # I expected that the shape passed from Python would reverse when it hits Fortran,
+    # but this did not happen.
+    shape = (200, 200, 10)
+    order = "F"
+    sst = np.full(shape, 290.0, order=order)
+    t_zt = np.full(shape, 280.0, order=order)
+    hum_zt = np.full(shape, 0.001, order=order)
+    u_zu = np.full(shape, 1.0, order=order)
+    v_zu = np.full(shape, -1.0, order=order)
+    rad_sw = np.full(shape, 0.0, order=order)
+    rad_lw = np.full(shape, 350.0, order=order)
+    slp = np.full(shape, 101000.0, order=order)
+    flux_skin(sst, t_zt, hum_zt, u_zu, v_zu, rad_sw, rad_lw, slp, algo=algo)

test/test_fortran.py

@@ -1,6 +1,6 @@
 import numpy as np
 import pytest
-from aerobulk import flux_noskin
+from aerobulk import flux_noskin, flux_skin
 
 
 @pytest.mark.parametrize(
@@ -28,7 +28,7 @@
         ),
     ],
 )
-def test_fortran_code(algo, expected):
+def test_fortran_noskin(algo, expected):
     # Compare python wrapper output with fortran results from (...)
 
     # inputs are the same for all test cases in the fortran example
@@ -49,3 +49,79 @@ def test_fortran_code(algo, expected):
     np.testing.assert_allclose(evap, expected["evap"])
     np.testing.assert_allclose(taux, expected["taux"])
     np.testing.assert_allclose(tauy, expected["tauy"])
+
+
+@pytest.mark.parametrize(
+    "algo, expected",
+    [
+        (
+            "coare3p0",
+            {
+                "ql": -81.3846741,
+                "qh": -15.1545086,
+                "evap": -2.87061906,
+                "taux": 3.57834995e-02,
+                "tauy": 0.0,
+                "t_s": 21.7219677,
+            },
+        ),
+        (
+            "coare3p6",
+            {
+                "ql": -83.0788422,
+                "qh": -15.3865499,
+                "evap": -2.93033028,
+                "taux": 3.21817845e-02,
+                "tauy": 0.0,
+                "t_s": 21.7057972,
+            },
+        ),
+        (
+            "ecmwf",
+            {
+                "ql": -80.2958984,
+                "qh": -14.3822346,
+                "evap": -2.83224440,
+                "taux": 3.84389125e-02,
+                "tauy": 0.0,
+                "t_s": 21.7325401,
+            },
+        ),
+    ],
+)
+def test_fortran_skin(algo, expected):
+    # Compare python wrapper output with fortran results from (...)
+
+    # inputs are the same for all test cases in the fortran example
+    rt0 = 273.15  # conversion to Kelvin
+    sst = rt0 + 22  # in Kelvin
+    t_zt = rt0 + 20
+    hum_zt = 0.012  # kg/kg
+    u_zu = 5  # m/s
+    v_zu = 0
+    slp = 101000.0  # Pa
+    rad_sw = 0
+    rad_lw = 350
+    ql, qh, taux, tauy, t_s, evap = flux_skin(
+        sst,
+        t_zt,
+        hum_zt,
+        u_zu,
+        v_zu,
+        rad_sw,
+        rad_lw,
+        slp,
+        algo=algo,
+        zt=2,
+        zu=10,
+        niter=10,
+    )
+    evap = evap * 3600 * 24  # convert to mm/day
+    t_s = t_s - rt0
+
+    np.testing.assert_allclose(ql, expected["ql"])
+    np.testing.assert_allclose(qh, expected["qh"])
+    np.testing.assert_allclose(evap, expected["evap"])
+    np.testing.assert_allclose(taux, expected["taux"])
+    np.testing.assert_allclose(tauy, expected["tauy"])
+    np.testing.assert_allclose(t_s, expected["t_s"])