Interpolation scheme (in time) used for JRA55 forcing

[GMSaran]

Hello,

Thank you for the great code!

I have a small doubt on how the atmospheric forcing/fluxes on the ocean is computed from the 3 hourly data of JRA55. I run global simulations using ClimaOcean.

An example case: I simulate a day, which results in the atmospheric velocity data having 9 time instances (assuming the instances are at 0,3,6,12....24 hours). I modified the atmospheric velocity as follows:

for j = 1:4
atmosphere.velocities.u.data[:,:,:,j] .= 0.0
end

for j = 5
atmosphere.velocities.u.data[:,:,:,j] .= 5.0
end

for j = 6:9
atmosphere.velocities.u.data[:,:,:,j] .= 0.0
end

If I simulate using the above modification, the maximum ocean horizontal velocity is around 10^(-7) up to 3 hours (which makes sense with no atmospheric forcing, velocity is low). But after 3 hours maximum ocean horizontal velocity starts ramping up. By the end of 6 hours, it is 7.7*10^(-3).

If I set atmosphere.velocities.u.data[:,:,:,:] .= 0.0 (everything is zero at all times), then after six hours horizontal velocity stays at 10^(-7).

So I was wondering how the fluxes are computed for each time step using the 9 original data slices.

Any answer would be much appreciated!

[glwagner]

From what you have reported, it sounds like there may be a bug. However, I cannot investigate further without seeing your code. Can you provide the entire code you are running so that the issue can be reproduced? If you can reduce your code to a case that is as simple and minimal as possible (for example, a single column model), that will help a lot.

Also just some misc nodes:

• by convention we use n as the time-index (j is the latitudinal / y-index). Using n will reduce confusion and speed up solving this problem!
• To set an entire time slice from a FieldTimeSeries, I suggest using

set!(atmosphere.velocities.u[n], 5)

this will work with functions and arrays as well as scalars in place of 5. I also think it is easier to read.

When writing code into github issues/discussions, please use triple backticks (```) to format your code. Long code will be hard to read without the formatting.

[GMSaran]

using ClimaOcean
using Oceananigans
using Oceananigans.Units
using Dates
using Printf
using Statistics
using NCDatasets

# We start by constructing an underlying TripolarGrid at ~1 degree resolution,

arch = CPU()
Nx = 360
Ny = 180
Nz = 5

####################################################################################################################################

depth = 4000meters

underlying_grid = TripolarGrid(arch; size = (Nx, Ny, Nz), halo = (5, 5, 4), z=(-4000,0) )

bottom_height = regrid_bathymetry(underlying_grid;
                                  minimum_depth = 100,
                                  interpolation_passes = 5,
                                  major_basins = 1)

grid = ImmersedBoundaryGrid(underlying_grid, PartialCellBottom(bottom_height; minimum_fractional_cell_height=0.2);
                            active_cells_map=true)

####################################################################################################################################

# ### Closures

eddy_closure = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=2e3, κ_symmetric=2e3)
vertical_mixing = ClimaOcean.OceanSimulations.default_ocean_closure()
horizontal_viscosity = HorizontalScalarDiffusivity(ν=4000)

####################################################################################################################################

# ### Ocean simulation

free_surface       = SplitExplicitFreeSurface(grid; substeps=70)
momentum_advection = WENOVectorInvariant(order=5)
tracer_advection   = WENO(order=5)

ocean = ocean_simulation(grid; momentum_advection, tracer_advection, free_surface, 
                         closure=(eddy_closure, horizontal_viscosity, vertical_mixing))

@info "We've built an ocean simulation with model:"
@show ocean.model

####################################################################################################################################

### Atmospheric forcing

atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(40), start_date = DateTime(1990, 1, 30), 
                                       end_date = DateTime(1990, 1, 31),
                                       include_rivers_and_icebergs = false)

####################################################################################################################################

    ctemp = 0

    atmosphere.velocities.v.data[:,:,:,:]    .= 0.0
    atmosphere.velocities.u.data[:,:,:,:]    .= 0.0
   
    for n = 1:4
      atmosphere.velocities.u.data[:,1:90,:,n] .= 0.0
    end 

    for n =  5
      atmosphere.velocities.u.data[:,1:90,:,n] .= 5.0
    end 

    for n = 6:9
      atmosphere.velocities.u.data[:,1:90,:,n] .= 0.0
    end 

    atmosphere.pressure.data[:,:,:,:] .= 100000

    atmosphere.tracers.T.data[:,:,:,:] .= ctemp + 273.15

    atmosphere.freshwater_flux.rain.data[:,:,:,:] .= 0.0
    atmosphere.freshwater_flux.snow.data[:,:,:,:] .= 0.0

    atmosphere.downwelling_radiation.shortwave.data[:,:,:,:] .= 0.0
    atmosphere.downwelling_radiation.longwave.data[:,:,:,:] .= 0.0

########################################################################################################################
 
#### Coupled simulation

atmosphere_ocean_flux_formulation =   CoefficientBasedFluxes(drag_coefficient = 1e-2,
                                      heat_transfer_coefficient = 0,
                                      vapor_flux_coefficient = 0)

interfaces = ClimaOcean.OceanSeaIceModels.ComponentInterfaces(atmosphere, ocean; 
                                                                    atmosphere_ocean_fluxes=atmosphere_ocean_flux_formulation,                                     
                                                                    ocean_heat_capacity = 1e13)

coupled_model = OceanSeaIceModel(ocean;
                 atmosphere, interfaces)

########################################################################################################################

#### Ready to run

simulation = Simulation(coupled_model; Δt=5minutes, stop_time=1440minutes)

wall_time = Ref(time_ns())

function progress(sim)

    AA      = sim.model.ocean.model.velocities.u.data

    BB      = sim.model.interfaces.atmosphere_ocean_interface.fluxes.x_momentum.data 

    amax    = maximum(AA)

    step_time = 1e-9 * (time_ns() - wall_time[])

    @info @sprintf("Time: %s, n: %d, Δt: %s, max|u|: (%.2e)\n",
                   prettytime(sim), iteration(sim), prettytime(sim.Δt),
                    amax)

     wall_time[] = time_ns()

     return nothing
end

add_callback!(simulation, progress, IterationInterval(1))

run!(simulation)

 

 

[glwagner]

Yes the formatting is good! I think it might be helpful to reduce the code to a single column case. There is an example of a single column simulation, and we can do the same thing in that code (eg set the velocity to 0 everywhere, and then set the velocity to some value for a single time point). It will be easier to disentangle issues with a single column model and much faster to run. Also, I don't think all of the model physics that are added are really necessary. We can simply use the default ocean_simulation(grid).

One question about your code...

you wrote

    ctemp = 0

    atmosphere.velocities.v.data[:,:,:,:]    .= 0.0
    atmosphere.velocities.u.data[:,:,:,:]    .= 0.0
   
    for n = 1:4
      atmosphere.velocities.u.data[:,1:90,:,n] .= 0.0
    end 

    for n =  5
      atmosphere.velocities.u.data[:,1:90,:,n] .= 5.0
    end 

    for n = 6:9
      atmosphere.velocities.u.data[:,1:90,:,n] .= 0.0
    end 

    atmosphere.pressure.data[:,:,:,:] .= 100000

    atmosphere.tracers.T.data[:,:,:,:] .= ctemp + 273.15

    atmosphere.freshwater_flux.rain.data[:,:,:,:] .= 0.0
    atmosphere.freshwater_flux.snow.data[:,:,:,:] .= 0.0

    atmosphere.downwelling_radiation.shortwave.data[:,:,:,:] .= 0.0
    atmosphere.downwelling_radiation.longwave.data[:,:,:,:] .= 0.0

This modulates the velocity field between 1:90, but not everywhere (1:180).

I think to investigate the problem we want, I suggest using this simpler more compact form, and also to change the atmosphere velocity everywhere (not just 1:90):

parent(atmosphere.velocities.u) .= 0
parent(atmosphere.velocities.v) .= 0
parent(atmosphere.tracers.T) .= 20 + 273.15
parent(atmosphere.freshwater_flux.rain) .= 0
parent(atmosphere.freshwater_flux.snow) .= 0
parent(atmosphere.downwelling_radiation.shortwave) .= 0
parent(atmosphere.downwelling_radiation.longwave) .= 0

# Set one time point to 10 m/s:
set!(atmosphere.velocities.u[5], 10)

Note that indexing into the field time series directly vs indexing into .data is the same.

[GMSaran]

Hello,

Sorry for the late reply and thank you so much for the answer! As you suggested, I tried the same type of forcing in an example code named "single_column_os_papa_simulation.jl" provided for the CLIMAOCEAN package.

Short summary: I modified the atmospheric forcing such that it is non-zero only at the 5th time index (hence at 12 hrs of simulation time). I monitored the simulation, and the velocity field was zero up to 9 hrs of simulation time and slowly started to evolve in time after that. This makes sense if interpolation is used. A specific non-zero forcing has an influence 3 hrs/1 time index (frequency of the JRA55 forcing files) before and after as well.

However, the original code I posted earlier has the same error and I cant figure out why? The error being the atmospheric forcing starts to have an influence nine hours/ 3 time indexes before, which is quite different from the single column case.

Any advice would be much appreciated, cheers. Please let me know if some other info will be useful, Im happy to provide it.

[glwagner]

The issue I pointed out refers to the fact that your code does not overwrite all of the atmospheric state at all time points, but only part of it at certain indices. Therefore it is expected that there is non-zero atmospheric forcing at all time points somewhere in the domain?

[GMSaran]

Apologies I forgot to mention that I modified the original code with the small snippet (given below) you gave. The results were same as the original code I initially was simulating.

parent(atmosphere.velocities.u) .= 0
parent(atmosphere.velocities.v) .= 0
parent(atmosphere.tracers.T) .= 20 + 273.15
parent(atmosphere.freshwater_flux.rain) .= 0
parent(atmosphere.freshwater_flux.snow) .= 0
parent(atmosphere.downwelling_radiation.shortwave) .= 0
parent(atmosphere.downwelling_radiation.longwave) .= 0
 
set!(atmosphere.velocities.u[5], 10)

For the above modifications, the velocity field starts evolving after 3 hrs of simulation time. Also If I comment out the final line (set!(atmosphere.velocities.u[5], 10)), then maximum value of velocity field stays at 10^(-7) m/s (so basically not evolving).

If I use

set!(atmosphere.velocities.u[6], 10)

the velocity field starts evolving at 6 hrs of simulation time.

cheers!

[glwagner]

Ok, so the summary is that you are finding that your modifications work for a single column simulation, but not for a non-single-column simulation? Is that right?

Can you try changing your setup so that you start from an "empty" atmospheric state (rather than JRA55)? This may avoid bugs associated with "forgetting" to zero some of the fields.