[MITgcm-devel] (not so) funny things happen in seaice_lsr and pickups

[Martin Losch]

Hi all,

I looks like I have found the problem that causes the "spontaneous"  
explosions in seaice_lsr (C-grid): Some of the metric terms were  
discretized in an, let's say, unfortunate way (actually, really wrong  
in a few places, all my fault of course). Also my fault: I never  
tested this part of the code properly (all the rigorous tests were  
carried out on cartesian grids or Dimitris CS510, where the metric  
terms are turned off).

It is very well possible (from my experience with an inconsistent  
discretization in the B-grid that caused the B-grid CS510 to crash),  
that these issues are responsible for Matt's problems, too. The  
crashes I observed are highly sensitive to the model trajectory and  
slight modifications already changed the time/timestep when the crash  
actually occurs. Using a fix like Matt's (use different thermodynamic  
code) just changes the solution, and the model may explode at a  
different time, that fortuitously is beyond Matt's integration period,  
so what works for Matt may not work otherwise.

Now my questions:
1. In order to code the metric terms correctly I need fields like  
tan(YG) and tan(YC) (in analogy to tanPhiAtU/V these would be called  
tanPhiAtZ and tanPhiAtC). On a lat/lon grid these fields should be  
identical to tanPhiAtC=tanPhiAtU, tanPhiAtZ=tanPhiAtV. However, I  
would prefer having these fields with the correct names. Do you agree?  
If you agree there's a second question:
2. Although these fields will only be used in pkg/seaice, I still need  
to initialize them in ini_grid (and subroutines) because of the  
rotated grid option (after ini_spherical_polar_grid, YC and YG are no  
longer the grid coordinate but the geographical coordinates). Should I  
define and set them with a CPP flag ALLOW_SEAICE/useSEAICE?
Alternatively on could derive tanPhiAtZ/C in seaice_init_fixed based  
on "averaging" or copying tanPhiAtU/V, but that method might break  
down, once we decide to bite the bullet and implement all metric terms  
also for the curvilinear grids.

How should I do it? the clean solution with 2 further, mostly unused  
2D grid-fields (defined in GRID.h), or the solution that is confined  
to pkg/seaice and include a hack for spherical grids and potentially  
more hacks for general curvilinear grids? While I prefer the first  
option, the second is perfectly OK, too.

Martin

PS. All my restart problems remain, and are really weird, but that's a  
different story.