Perturbation to basic state variables

[Sougata18]

Hi! I want know how to perturb the basic state variables such as temperature, salinity etc. for a specific location in the global 1_deg model.
Thank you!

[dionhaefner]

There's no general recipe for this. You will have to find the index of the desired location (by looking it up in the coordinate arrays vs.xt, vs.yt), then change the variables at that locations (either in set_initial_conditions or in set_forcing, depending whether you want it to happen once or in every time step).

[Sougata18]

Sorry but I'm kinda new to this model, so can you just elaborate a little on the procedure to perturb the variable. For example, If I want to decrease salinity gradually (with a fixed amount per year) then what do I have to do ?

[dionhaefner]

The simplest way would probably be to add / subtract a constant term here:

https://github.com/team-ocean/veros/blob/master/veros/setups/global_1deg/global_1deg.py#L328

A negative term would indicate a constant freshwater influx at the surface which gradually reduces salinity.

[Sougata18]

For understanding the model code ,what can be a good start ?

[dionhaefner]

Besides what's in the docs, we unfortunately don't have any beginner's resources yet. If you have a specific question and share what you tried I am happy to help.

[Sougata18]

How to change the runtime of the global_1deg setup ?
I have tried settings.runlen=86400*365*3 to run it for 3 years but the output files does not contain any time step , i.e., Dimension of Time is 0.

[dionhaefner]

That sounds correct. Which commands did you use to create and run the setup? And can you share the output of the run?

[Sougata18]

I just changed the settings.runlen in the global_1deg.py setup from it's default settings.runlen=10*settings.dt_tracer to settings.runlen=86400*365*3 and ran the setup file.
files.zip

[dionhaefner]

Sorry, I meant the output to the terminal. Did it actually run for 3 years (which should take several days of real-time)?

[Sougata18]

I am actually using a HPC system.
I noticed that it contains a runtime error: solution diverged at iteration 37.
I just copy pasted the terminal output in a text file which is attached below. The error is something related to the sqrtsa = npx.sqrt(sa) in gsw.py in /veros/veros/core/density.

I have another question as you mentioned that the real run would take several days in real time; how much would it take if I want to do a 80-100 year simulation on a HPC system ?

Toutput.zip

[dionhaefner]

I noticed that it contains a runtime error: solution diverged at iteration 37.
I just copy pasted the terminal output in a text file which is attached below. The error is something related to the sqrtsa = npx.sqrt(sa) in gsw.py in /veros/veros/core/density.

So the solution became unstable and diverged. It seems to me like you also increased the time step? Like, the CFL diagnostic runs every model day, which should be at 48 time steps with the default step of 0.5 hours. But in the log I see its output after 24 time steps. That would explain why the setup diverged; the default values are close to the largest possible stable time steps.

I have another question as you mentioned that the real run would take several days in real time; how much would it take if I want to do a 80-100 year simulation on a HPC system ?

That depends how much hardware you want to throw at it :) On a single CPU node you will probably get something like 1 yr/day with JAX. By using multiple nodes or GPUs you can probably push that to something like 10 yr/day. But it all depends on your hardware and cluster architecture.

[Sougata18]

Thanks...changing the CFL diagnostic time step from 24 to 48 worked; also checked the default global_1deg code and the CFL diagnostic time step is 24 there.

[dionhaefner]

I meant the main model time step (dt_tracer and dt_mom). The CFL diagnostic has no influence on the model, I just used it to deduce that you must have changed the time step from the default of 1800.

[Sougata18]

Okay!!
yes I also changed that before. But in the docs I guess I have read that the dt_tracer can be anything >= dt_mom. But anyways got it now.

[Sougata18]

If I'm doing 100 years of simulation, what would be the approx spin up time for global_1deg.py ?

[dionhaefner]

Depends on what part of the ocean you need spun up. 100 years sounds about right for equilibration of the MOC, but it's always safer to check (plot some mean variable over time and make sure it's converged). Or ask your advisor :)

[Sougata18]

I have been trying to decrease the salinity in the North Atlantic using the 4 degree model. But the method seems to be not working. Did I miss something ?

I made a mask of 1s and 0s stored in the custom variable vs.freshwater

salt_mask_raw=xr.open_dataset('freshwater_mask.nc')
salt_mask_data = salt_mask_raw.freshwater_mask
longitude=salt_mask_data.longitude
latitude=salt_mask_data.latitude

salt_mask_data_interpolate = veros.tools.interpolate((latitude,longitude),salt_mask_data,(vs.xt,vs.yt), kind='nearest')
 vs.freshwater = vs.freshwater.at[:, :].set(salt_mask_data_interpolate)

Then used the following code to decrease salinity in the mased region -

ampl = 1.0  # Decrease rate of 10 PSU per year
t_start = 4055.56  # Start time of the salinity decrease
t_duration = 10.0  # Duration of the salinity decrease in years
t_sim = vs.time / (360.*86400.)  # Current simulation time in years

cond = jnp.less_equal(t_sim-t_start, t_duration)
decrease_years = jnp.where(cond, jnp.minimum(t_sim - t_start, t_duration), 0.0)
forc_val = jnp.where(cond, ampl * decrease_years, 0.0)
vs.forc_salt_surface = jnp.where(cond, vs.forc_salt_surface - (vs.freshwater * forc_val), vs.forc_salt_surface)

[dionhaefner]

Define not working?

Code looks fine to me for the most part, but please format it with triple backticks next time ```

The part that stands out to me is this:

cond = jnp.less_equal(t_sim-t_start, t_duration)
decrease_years = jnp.where(cond, jnp.minimum(t_sim - t_start, t_duration), 0.0)
forc_val = jnp.where(cond, ampl * decrease_years, 0.0)
vs.forc_salt_surface = jnp.where(cond, vs.forc_salt_surface - (vs.freshwater * forc_val), vs.forc_salt_surface)

That seems complicated for a simple ramp. I would implement this like so:

forc_val = ampl * jnp.clip((t_sim - t_start) / t_duration, 0, 1)
vs.forc_salt_surface = vs.forc_salt_surface - (vs.freshwater * forc_val)

If this still doesn't work I'd recommend you run your setup with NumPy on your local machine and plot some of the variables here (forc_salt_surface, freshwater).

[Sougata18]

So the code didn't work.
But to plot force_salt_surface and freshwater it has to be added to the digontics before spin up, right ?

[dionhaefner]

You can use any Python code when running with NumPy:

def set_forcing(self, ...):
    ...
    import matplotlib.pyplot as plt
    plt.imshow(vs.force_salt_surface)
    plt.show()