Issue with Prather scheme

[fmasson]

Hi everyone,

Context

I'm running nemo_3.6_stable, r6396, in ORCA1L75 using the configuration files of ShacoNemo?. My atmospheric forcing is a perturbed version of the Drakkar Forcing Set 5.2. Using a perturbed version of DFS5.2 probably increases the likelihood of the crash described below, but as you can read below the problem happens for realistic values of sea ice drift, so I think the problem is more structurally related to a NEMO issue. Before posting this message, I exchanged e-mails with Martin Vancoppenolle and Clément Rousset, who confirm that this problem can happen anyway.

Analysis

I'm compiling NEMO with the -fpe0 flag to make sure I catch floating point exceptions. Since a few months I've had a 'floating divide by zero' error every now and then. Instead of working around (changing time step, diffusivity parameters, tweaking restarts...) I decided to face the problem in depth and trace the error back to its origin. Hopefully this problem can also make ORCA1 coupled more stable and ORCA025 more stable too.

The problem is the following. The advection routines of LIM cause the problem. If the time step is odd, advection is first computed along direction y, then along direction x (if time step is even, it's the reverse). This step of advection is repeated over several sub-time steps (variable initad).

Excerpt of routine lim_trp, when jt is odd:

 DO jt = 1, initad
               CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
               
               CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &
                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
              

The calls to routines lim_adv_y and lim_adv_x imply the use of several arrays: the u- and v- velocities of the ice, the field to be advected, the area it occupies, and the first and second order moments of the advection scheme used (Prather's scheme).

By having a closer look at the routine lim_adv_y one can see that the array with area occupied by the field, named psm in the routine, can have negative values under certain circumstances:

Excerpt of routine lim_adv_y:

DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0
         DO ji = 1, jpi
            ...
            zalf         =  MAX( 0._wp, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj)
            ...
            zfm (ji,jj)  =  zalf  * psm(ji,jj)
            ...
            !  Readjust moments remaining in the box.
            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
            ...

In my case the scalar variable zalf is slightly above 1.0 while variables involved in its calculation are not particularly irrealistic:

In my case:

zalf = MAX(0.0, 2.81860230041890) * 10800 * 29917.5891308800 / 885576195.653601 = 1.0283...

As a consequence one entry of the array psm becomes negative.

The code then exits successfully the routine lim_adv_y, but when it enters the x-advection routinelim_adv_x, one entry of the array psm (that has not changed) is reset to 1e-20:

Excerpt of routine lim_adv_x

      psm (:,:)  = MAX( pcrh * e12t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )

(pcrh is a switch equal to 0 when advection is in the x- direction).

Later in the routine, we have a division by psm:

            zalf         =  MAX( 0._wp, put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji,jj)
            zalfq        =  zalf * zalf

But since one entry of psm is 1e-20, the value of zalf becomes 1e20 for one grid point; even worse, the value of zalfq (square of zalf) is 1e40.

The value advected becomes unrealistically high, the max of sea ice concentration quickly grows to 1e40 or bigger, and the model blows up eventually in the mechanical redistribution routine (limitd_me) because of the unrealistic high value of sea ice concentration.

So, the true problem is that the parameter zalf can have values larger than 1.0, while it shouldn't.

Fix

There are two possibilities:

• Upper-bound zalf by 1.0; but in that case the Prather scheme may not be conservative and leaks may appear.
• Run with more sub time-steps. I haven't tried that yet.

Best wishes,

François Massonnet

[vancop]

Thanks François, I think it is just a violation of CFL, because of your strong wind forcing. A solution would indeed be to implement more sub-time steps.
Best
Martin