[MITgcm-devel] seaice dynamics in ocean 2-D

[Martin Losch]

Hi Jean-Michel,

that’s annoying. I have changed the dynamics solver(s) of pkg/seaice a lot since then, see attached file changes_to_seaice.txt, where I compiled all (documented) changes to pkg/seaice since ckpt65j. Below are the ones that I think may have affected your experiement, with comments; the most likely candidate is this: checkpoint65r (2015/12/21), seaice_lsr.F r1.119. I hope this helps. 

Martin

*** this will definitely change your results at the truncation level ***
checkpoint65m (2015/06/15)
o pkg/seaice:
  - seaice_calc_viscosity.F (r1.30): change computation of deltaC**2 to ensure
    positiveness, modify a few comments and improve variable names,
  - the changed computation of deltaC affects most seaice-related
    experiments at the truncation level.
    seaice_itd.thermo is affected by far the most, with
    only 7 digits of agreement remaining; update experiments.

*** this is unlikely a problem, but you never know, especially with the new parameter that may lead to a different default ***
checkpoint65n (2015/07/29)
o pkg/seaice:
  - Rename SEAICE_VECTORIZE_LSR_ZEBRA to SEAICE_LSR_ZEBRA and define it
    in lab_sea forward experiment (changes results -> updated)
o pkg/seaice:
    add new parameter SEAICE_deltaMin that is used ***only***
    for regularizing Delta (and nothing else, like the parameter
    which is also used for all sorts of things SEAICE_EPS).
    Defaults to SEAICE_EPS for backward compatibility only.

*** this should not affect your results, but it’s quite a big change to many subroutines (seaice_calc_lhs.F, seaice_dynsolver.F, seaice_evp.F, seaice_lsr.F, seaice_preconditioner.F) ****
checkpoint65o (2015/09/14)
o pkg/seaice: introduce new parameter SEAICEscaleSurfStress (default = .FALSE.)
  - if TRUE scale ice-ocean and ice-atmosphere stress acting on sea ice
    by AREA according to Connelly et al. (2004), JPO.
  - For EVP and the LSR solver, the implicit part of the drag term is
    the only non-zero term in the denominator when there is no ice (which
    prevented running the model with zero ice-ocean drag). If
    SEAICEscaleSurfStress = .TRUE., this results in a division by zero
    (or zero main diagonals BU/BV) which need to be caught. A practical
    consequence is that for open water, the momentum equation reduce
    to trivially 0 = 0 + 0 (for EVP). For LSR, BU/BV are reset to 1, if
    they happen to be zero, often leading to a non-zero solution over
    open water.

*** since you are using SEAICE_OLx/y=1, this is actually a very likely candidate, before that all LSR solutions with SEAICE_OLx/y>0 where clearly wrong (stripes in solutions along tile edges), JFNK solutions also improved, but not significantly, because LSR is only used as a preconditioner and the Newton step fixes this “automatically"***
checkpoint65r (2015/12/21)
o pkg/seaice: (seaice_lsr.F r1.119)
    fix bug in tridiagonal solvers for SEAICE_OLx/y>0, this affects the JFNK
    solver (offline.dyn_jfnk) and global_ocean.cs32x15.seaice -> updated

*** another candidate, but unlikely to have caused the large differences ***
checkpoint65t (2016/02/10)
o pkg/seaice/seaice_lsr.F: replace 1./SEAICE_deltaTdyn with recip_deltaT,
  affects some experiments at truncation level:
  - global_ocean.cs32x15.icedyn (11 digits),
  - global_ocean.cs32x15.seaice (11 digits),
  - seaice_itd (12 digits), seaice_itd.lipscomb07 (12 digits) -> update
  - also update global_ocean.cs32x15 adjoint and tangent-linear results


-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: changes_to_seaice.txt
URL: <http://mitgcm.org/pipermail/mitgcm-devel/attachments/20160629/59801116/attachment.txt>
-------------- next part --------------

> On 29 Jun 2016, at 01:34, Jean-Michel Campin <jmc at mit.edu> wrote:
> 
> Hi Martin,
> 
> I am trying to reproduce results from ~1.5 year ago (close to checkpoint65j), 
> with a 2-D (y-z) southern ocean channel, and am am getting very different seaice 
> velocity.
> I am using pkg/thsice thermodyn with pkg/seaice dynamics, deltaT=1.h
> and my data.seaice is attached.
> 
> I can reproduce the old results when I use:
> - full code from checkpoint65j (Feb 25, 2015) 
> - or just pkg/seaice from checkpoint65j into latest code
> 
> I made a short (10.iter) comparison between c65j (d12a) and latest code (d12c),
> starting from a pickup of the old run. The largest difference are on vIce as 
> seem on this plot (x-axis = latitude, magenta and red are d12a results
> at iter 1 and 9 respectively; green, yellow, ceilan and blue are d12c results
> at iter 1 3 5 and 9).
> 
> Do you have any suggestion of what could cause these large difference ?
> 
> Cheers,
> Jean-Michel
> <data.seaice><fig_uIce_d12.jpg>_______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel