[MITgcm-devel] adding ice dyn run time switches

[Jean-Michel Campin]

Hi Martin and Gael,

I thought this was on mitgcm-devel (but was not), and since it might
interest also others, I am switching to devel to add few comments.

Regarding this point:
> > I don't like these ****formula flags,  ...
> > I do see, that if there are more than two options,  ...
I think that it's OK when dealing with multiple choices (with one excluding 
the others), and a better documentation would avoid this problem of
not knowing which to pick.
It would be nice to have this ****formula or **select** things documented
in several places, in SEAICE_PARAMS.h, in seaice_summary.F (I tried to 
tweak the config_summary.F to report such flags, e.g., selectCoriMap:
(PID.TID 0000.0001) selectCoriMap = /* Coriolis Map options (0,1,2,3)*/
(PID.TID 0000.0001)                       2
(PID.TID 0000.0001)     0= f-Plane ; 1= Beta-Plane ; 2= Spherical ; 3= read from file
(PID.TID 0000.0001)     ;
but also where they appear in one of verification/*/input*/data.seaice
and ultimately in the manual (pkg/seaice section).

I have a question (to Martin) regarding seaice_summary.F :
I was tempted (and may have started) to not report parameters
that are irrevelant (e.g., tAlpha & sBeta when eosType <> LINEAR).
This would make the summary shorter and, may be, more useful
(assuming that all the relevant params are reported, which is 
 far to be the case right now, ice/snow/air density, air CP ...).
The first thhing I had in mind was to skip all the seaice-growth params
when we use thsice (or usePW79thermodynamics=F).
But I saw that you went a different path by adding a kind of
warning before printing the param value, like this one:
(PID.TID 0000.0001) MAX_HEFF has no effect because SEAICE_CAP_HEFF is undefined
(PID.TID 0000.0001) MAX_HEFF          = /* maximum ice thickness */
(PID.TID 0000.0001)                 1.000000000000000E+01
(PID.TID 0000.0001)     ;
So, which way should we follow ?

And to Gael: did not see much information in doc/tag-index 
related to your changes made on Friday (evening).

Cheers,
Jean-Michel

On Sat, Feb 11, 2012 at 12:09:30AM -0500, Gael Forget wrote:
> Hi Martin,
> Thanks for the replies. Have
> Not done anything about it yet.
> Will give more thought & email tomorrow,
> But I just wanted to let you know,
> The heavier clean up stuff is done.
> The rest should be transparent.
> Cheers,
> Gael
> 
> Sent from my iPhone
> 
> On Feb 10, 2012, at 3:30 AM, Martin Losch <Martin.Losch at awi.de> wrote:
> 
> > Hi Gael,
> > 
> > did not see any vectorization problems so far; I am not testing all flag combinations, but the ones I did did not cause a performance problem. I'll let you know if something problematic occurs. BTW, the ifort compiler has some vectorization capabilities, of if you want to you can test your code with that. A code-listing with compiler messages will give you some clue of the vectorization, try -vec_report5 or -vec_report[5] (see ifort man-page for more details).
> > 
> > I am not against adding more ice-strength formulations, on the contrary, I think we need this (although this relationship shouldn't be used to avoid too thick ice. I think that this "feature" has a thermodynamic cause). I though that Gunnar tried this. At least I talked to him about doing that. 
> > Can I ask you to include these flags around do loops, so like this:
> >      IF ( SEAICE_pressFormula==1 ) THEN
> >        DO j=1-Oly+1,sNy+Oly
> >         DO i=1-Olx+1,sNx+Olx
> >          PRESS0(i,j,bi,bj) = PSTAR*HEFF(i,j,bi,bj)
> >     &         *EXP(-20.0 _d 0*(ONE-AREA(i,j,bi,bj)))
> >         ENDDO
> >        ENDDO
> >    ELSEIF ( SEAICE_pressFormula==2 ) THEN
> >        DO j=1-Oly+1,sNy+Oly
> >         DO i=1-Olx+1,sNx+Olx
> >           PRESS0(I,J,bi,bj)=PSTAR
> >      &         *MAX(HEFF(i,j,bi,bj),HEFF(i,j,bi,bj)*HEFF(i,j,bi,bj)*
> >      &                              HEFF(i,j,bi,bj)/4. _d 0)
> >      &         *EXP(-20.0 _d 0*(SEAICE_area_max-AREA(i,j,bi,bj)))
> >         ENDDO
> >        ENDDO
> >     ELSE
> >       etc.
> >      ENDIF
> >        DO j=1-Oly+1,sNy+Oly
> >         DO i=1-Olx+1,sNx+Olx
> >          Zmax(I,J,bi,bj)   = SEAICE_zetaMaxFac*PRESS0(I,J,bi,bj)
> >          ZMIN(i,j,bi,bj)   = SEAICE_zetaMin
> >          PRESS0(i,j,bi,bj) = PRESS0(i,j,bi,bj)*HEFFM(i,j,bi,bj)
> >          TAUX(i,j,bi,bj)   = 0. _d 0
> >          TAUY(i,j,bi,bj)   = 0. _d 0
> > #ifdef SEAICE_ALLOW_FREEDRIFT
> >          uice_fd(i,j,bi,bj)= 0. _d 0
> >          vice_fd(i,j,bi,bj)= 0. _d 0
> > #endif
> >         ENDDO
> >        ENDDO
> > 
> > [Why are the loops not over the full range: I extended them in revision 1.4 to the present values to fix some problem, but they can probably be extened to cover the entire range, better for vectorization (o:, Could you test/fix that  while you are at it? If not, I'll do it later]
> > 
> > The same applies to your second idea, although I do feel that this is going to be more difficult. Again, please put these runtime flags around the do-loops (or use multipliers). I am a little concerned that the nearly un-intelligible code of e.g. seaice_lsr, will become even less intelligible if even more options are added so I plea that you do have this in mind when you make the changes. Do you plan to do this simulaneously for VP and EVP code? 
> > 
> > BTW, the stress tensor is zero when you see PRESS0=0 (and also ZMIN=0, but that's the default). If you want to keep the "press" contribution (I have not understood, when/why this can be necessary), then you'd have to set press0=0, zeta=0, eta=0, after calling seaice_calc_viscosities (or even within), in seaice_evp, it's probably enough to reset
> >             seaice_div    (I,J) = - pressC(I,J) * hEffM(I,J,bi,bj)
> >             seaice_tension(I,J) = 0. 
> >             seaice_shear  (I,J) = 0.
> > 
> > about the flag names:
> > 
> > seaiceDynPstar ( if false, then there is no stress either )
> > seaiceDynTilt  (I don't like "Tilt", but I guess this is the term that seaice modelers use)
> > seaiceDynRheo => seaiceDynStressTensor (if this is false, you still have "Rheo"logy)
> > The last one is problematic, because the pressure term (-P/2) is part of the constitutive VP-relation, and also part of the stress tensor. I find it a bit strange to separate these. Unless you generate a cavitating fluid. In that case I'd choose names according to this, e.g. seaiceDynPstar => seaiceDynStressTensor (if false there is no stress tensor) and
> > seaiceDynRheo => seaiceUseCavitatingFluid (implies all terms of the stress tensor = 0, except for the -P/2, which then makes also sense to the "uninitiated" user)
> > 
> > I don't like these ****formula flags, because it's not obvious from the name, what the actual value means. Recently I have "taught" a few people to use the MITgcm, and these flags are a nuisance when trying to explain a namelist to a newbie, because you always have to look up which integer you want to use, and then you have to remember, where this is described, etc....
> > I do see, that if there are more than two options, it's complicated to use something else, so I'll go along with it (grinding my teeth).
> > 
> > About your motivation: "last time I tried the tilt term I could not run for long (rstar stops)" 
> > What does that mean? With the surface elevation driving the ice, the model does not run? I've never seen that in my runs (and I do use the rstar code too). Are you sure that this isn't masking some other problem?
> > 
> > Martin
> > 
> > On Feb 9, 2012, at 11:39 PM, Gael Forget wrote:
> > 
> >> Hi Martin,
> >> 
> >> I would like to do two things in the seace dynamics, and would greatly 
> >> appreciate your opinion, or potential objection. Jean Michel 
> >> did not object. I believe they both are worth testing 
> >> and could be useful long-term At first I would certainly 
> >> add warnings to STDERR stating "potentially harmful", 
> >> "more testing needed" or something like that.  
> >> 
> >> The first thing is a SEAICE_pressFormula (tentative name) switch 
> >> that for a value of 1 would do the usual way and for 2 would be like
> >> <           PRESS0(I,J,bi,bj)=PSTAR
> >> <      &         *MAX(HEFF(i,j,bi,bj),HEFF(i,j,bi,bj)*HEFF(i,j,bi,bj)*
> >> <      &                              HEFF(i,j,bi,bj)/4. _d 0)
> >> <      &         *EXP(-20.0 _d 0*(SEAICE_area_max-AREA(i,j,bi,bj)))
> >> This ramps up the pressure term for very thick ice. I cannot run 
> >> my real fw setup with either this or the dreadful CAP_HEFF. 
> >> I know that JPL has played such games so, while I have not seen JPL 
> >> offer any code to MITgcm yet, SEAICE_pressFormula should make sense 
> >> long term and my version is a working placeholder I think. Makes sense?
> >> 
> >> The second thing are switches for the various terms in the dynamics. 
> >> Something like seaiceDynPstar, seaiceDynTilt, and seaiceDynRheo (tentative names) 
> >> that, when false, would respectively exclude the pstar term, the SSH tilt term, and 
> >> the rheology (stress tensor) terms. The last one I have not worked out yet.
> >> There are three motivations : 1) last time I tried the tilt term I could not run 
> >> for long (rstar stops) so I carry around and keep up to date an edited version
> >> of seaice_dynsolver.F, which does not make my life easier. 2) it would allow to 
> >> easily go from free drift, to cavitating fluid, to vp, and I would like to quantify 
> >> the relative importance of those changes 3) assuming that free drift or 
> >> cavitating fluid is a decent enough approximations in some cases (e.g. at coarse res, 
> >> or with a climatogical forcing) it would be useful for adjoint/estimation purposes.
> >> Makes sense?
> >> 
> >> In case you dont object, would you have better suggestions regarding the params names?
> >> 
> >> Side question : did last night's checkin deteriorate vectorization?
> >> 
> >> Cheers,
> >> Gael
> >> 
> >