coupled

Author: simone-silvestri

ext/ClimaOceanPythonCallExt/veros_ocean_simulation.jl

@@ -1,5 +1,6 @@
 using CondaPkg
 
+using Oceananigans.Grids: topology
 using ClimaOcean.OceanSeaIceModels: reference_density, heat_capacity, SeaIceSimulation
 
 import Oceananigans.Fields: set!
@@ -43,10 +44,10 @@ function remove_outputs(setup::Symbol)
 end
 
 const CCField2D = Field{<:Center, <:Center, <:Nothing}
-const FCField2D = Field{<:Face, <:Center,   <:Nothing}
+const FCField2D = Field{<:Face,   <:Center, <:Nothing}
 const CFField2D = Field{<:Center, <:Face,   <:Nothing}
 
-function set!(field::CCField2D, pyarray::Py, k=pyconvert(Int, pyarray.shape[2]))
+function set!(field::Field2D, pyarray::Py, k=pyconvert(Int, pyarray.shape[2]))
     array = PyArray(pyarray)
     Nx, Ny, Nz = size(array)
     set!(field, view(array, 3:Nx-2, 3:Ny-2, k, 1))
@@ -56,7 +57,9 @@ end
 function set!(field::FCField2D, pyarray::Py, k=pyconvert(Int, pyarray.shape[2]))
     array = PyArray(pyarray)
     Nx, Ny, Nz = size(array)
-    set!(field, view(array, 3:Nx-2, 3:Ny-2, k, 1))
+    TX, TY, _  = topology(field.grid)
+    i_indices  = TX == Periodic ? UnitRange(3, Nx-2) : UnitRange(2, Nx-2)
+    set!(field, view(array, i_indices, 3:Ny-2, k, 1))
     return field
 end
 

ext/ClimaOceanPythonCallExt/veros_state_exchanger.jl

@@ -83,11 +83,13 @@ end
 @inline get_radiative_forcing(ocean::VerosOceanSimulation) = nothing
 
 function fill_up_net_fluxes!(ocean::VerosOceanSimulation, net_ocean_fluxes)
-    nx = pyconvert(Int, ocean.setup.state.settings.nx)
-    ny = pyconvert(Int, ocean.setup.state.settings.ny)
-    t1 = parent(net_ocean_fluxes.u)[:, 1:44, 1]
-    t2 = parent(net_ocean_fluxes.v)[:, 1:44, 1]
+    nx = pyconvert(Int, ocean.setup.state.settings.nx) + 4
+    ny = pyconvert(Int, ocean.setup.state.settings.ny) + 4
+    t1 = view(parent(net_ocean_fluxes.u), 1:nx, 1:ny, 1)
+    t2 = view(parent(net_ocean_fluxes.v), 1:nx, 1:ny, 1)
 
+    # TODO: Remove this when they fix veros to
+    # be able to force with an array instead of Climatology
     ta = zeros(size(t1)..., 12)
     tb = zeros(size(t2)..., 12)
     for t in 1:12
@@ -97,5 +99,16 @@ function fill_up_net_fluxes!(ocean::VerosOceanSimulation, net_ocean_fluxes)
 
     veros_set!(ocean, "taux", ta)
     veros_set!(ocean, "tauy", tb)
+
+    # TODO: Add heat flux and salinity flux when they
+    # fix veros to be able to force with prescribed boundary
+    # conditions rather than restoring
+
+    # t1 = view(parent(net_ocean_fluxes.T), 1:nx, 2:ny+1, 1)
+    # t2 = view(parent(net_ocean_fluxes.S), 1:nx, 2:ny+1, 1)
+
+    # veros_set!(ocean, "temp_flux", t1)
+    # veros_set!(ocean, "salt_flux", t2)
+
     return nothing
 end