How about setting up a vertical temperature penetration effect in the ocean and solve the NaN problem？

[wpan081]

Hi all! I have the following question for you all, thanks a lot for reading! I've read discussion #3218 and #2832, but it seems that my problem is still not solved.
For the setting of temperature in the general arithmetic example, we usually only set the sea table to FieldBoundaryConditions, e.g.:

@inline surface_temperature_flux(x,y,t,p) = - (p.Q⁰ * max(cos(2π * (t / 86400 - 0.5)), 0) + p.Qʰ) / (p.ρₒ * p.cᴾ) #diurnal cycle
Q⁰ = 500  # W m⁻², maximum solar radiation
Qʰ = -159.15  # W m⁻², surface heat flux
ρₒ = 1026 # kg m⁻³, average density at the surface of the world ocean
cᴾ = 3991 # J K⁻¹ kg⁻¹, typical heat capacity for seawater
#We impose a temperature gradient `dTdz` both initially and at the bottom of the domain, 
#culminating in the boundary conditions on temperature,
dTdz = 0.01 # K m⁻¹
T_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(surface_temperature_flux, parameters=(Q⁰=Q⁰, Qʰ=Qʰ, ρₒ=ρₒ, cᴾ=cᴾ)),
                                bottom = GradientBoundaryCondition(dTdz))

But in reality, the vertical change of temperature in the ocean conforms to the following equation 1: I_0*(A * exp(-z / g1) + (1.0 - A) * exp(-z / g2)). It has to do with the propagation of light. Then I use the forcing term to reflect this role. The code is:

T0(z) = T₀ + dTdz * z
gaussian_mask = GaussianMask{:z}(center=-grid.Lz, width=grid.Lz/10)
T_sponge = Relaxation(rate = 1/hour,
                      target = LinearTarget{:z}(intercept=T₀, gradient=dTdz),
                      mask = gaussian_mask)
#Derive equation 1
@inline temperature_forcing(x, y, z, t, p) = -(p.Q⁰ * max(cos(2π * (t / 86400 - 0.5)), 0) + p.Qʰ) / (p.ρₒ * p.cᴾ)* ((-p.A / p.g1) * exp(-z / p.g1) + (-(1.0 - p.A)/ p.g2) * exp(-z / p.g2))
#@inline temperature_forcing(x, y, z, t, p) = -(p.Q⁰ * max(cos(2π * (t / 86400 - 0.5)), 0) + p.Qʰ) / (p.ρₒ * p.cᴾ)*(p.A * exp(-z / p.g1) + (1.0 - p.A) * exp(-z / p.g2))
T_forcing = Forcing(temperature_forcing,
                    parameters = (Q⁰=Q⁰, Qʰ=Qʰ, A=0.7, g1=0.4, g2=8.0, ρₒ=ρₒ, cᴾ=cᴾ))
model = NonhydrostaticModel(; grid,
                advection = WENOFifthOrder(),
              timestepper = :RungeKutta3,
                  tracers = (:T, :S),
                 coriolis = FPlane(f=1e-4),
                 buoyancy = buoyancy,
                  closure = AnisotropicMinimumDissipation(),
                  stokes_drift = UniformStokesDrift(∂z_uˢ=∂z_uˢ),
      boundary_conditions = (u=u_bcs, T=T_bcs, S=S_bcs),
                  forcing = (u=u_sponge, v=v_sponge, w=w_sponge, T=(T_forcing, T_sponge), S=S_sponge))

And the TimeStepWizard set as:
wizard = TimeStepWizard(cfl=0.5, max_change=1.1, max_Δt=0.5minute)
But now the problem that will be encountered is as follows:

[ Info: Initializing simulation...
Iteration: 0000, time: 0 seconds, Δt: 11 seconds, max(|w|) = 8.2e-06 ms⁻¹, wall time: 0 seconds
[ Info:     ... simulation initialization complete (57.660 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (57.797 seconds).
[ Info: time = 1671.874849202, iteration = 100: NaN found in field u. Stopping simulation.
Iteration: 0100, time: 27.865 minutes, Δt: 28.531 seconds, max(|w|) = 0.0e+00 ms⁻¹, wall time: 2.588 minutes

the simulation cannot converge. What should I do? Thanks a lot!

[glwagner]

The depth dependence should probably be implemented with exp(z), not exp(-z). I can't see the code you use to create your grid, but generally z increases upward (think of "z" as "height", z=0 at the surface and z is negative within the ocean).

One way to check this is to plot your forcing function temperature_forcing as a function of z to inspect that it behaves as expected.

[wpan081]

Thanks a lot! It's amazing! I removed the minus sign in front of the z and now the problem is solved, thanks again!