[MITgcm-support] Icefront returnin NaN based on topology.

Martin Losch Martin.Losch at awi.de
Thu Feb 20 03:45:27 EST 2014 

Hi Harry,

as Dimitris already said, melting and refreezing only affects the ocean. The ice doesn’t move, so you will not be able to simulate glacier retreat scenarios with these packages.

I think your nan’s are most likely the consequence of a CLF-violation. Have a look at the monitor output of advcfl_uvel etc. These numbers should increase to values above 1 a few timesteps before the first nan if I am right about this.

If you want to suppress the effect of melting/freezing on the ocean (for test purposes) you should set the exchange coefficients to zero. In data.shelfice, these are SHELFICEheat/saltTransCoeff and probably similar for data.icefront

M.

On Feb 20, 2014, at 7:14 AM, Harry A <gmtonys3 at hotmail.com> wrote:

> Hi Martin
> 
> Sorry, I thought you were just wanting to see the ins and outs of my setup.
> 
> I had previously had the pressure as a 2-d Grid a in the example case but then i changed it to 3-D as an experiment and it made the model produce results, so I kept it in after that, although it was probably not the cause of the fix.
> 
> I had the icefront and shelfice latent heats and heat capacities very large in there as I didn't want any of the ice wall to melt and hence reduce in size. Having removed those lines, so using the default values, the model now works and produces useful results. However, if I leave it going for more than a few hours simulation time it breaks down again and returns NaN answers, as the temperature seems to diverge quite rapidly: currently the first layer near the ice wall drops to about -1 C within 2 hours, which is what I'd hoped to avoid using those parameters.
> I've now experimented with setting the latent heat and heat capacity in both locations to be much higher, 9e10, and in this case the model does not diverge so badly.
> 
> I will try running it with smaller time steps to see if I can get it to stay on track with the default values for latent heat and heat capacity, as it may need some time to spin up.
> 
> So, i think my final question is just (which may be very simple and obvious), do the icefront and shelfice packages let the ice melt over time? i.e. if I start the initial shelfice with some depth in a certain location, will this depth change over time in the model, or does setting the ice up make the model keep it that way, so although it may remove energy from the water, it doesn't actually shrink as the model progresses...
> 
> Thank you so much for all of your assistance!
> 
> Harry
> 
> 
> 
> > From: Martin.Losch at awi.de
> > Date: Wed, 19 Feb 2014 11:09:00 +0100
> > To: mitgcm-support at mitgcm.org
> > Subject: Re: [MITgcm-support] Icefront returnin NaN based on topology.
> > 
> > Hi Harry,
> > 
> > the idea was, that I try to reproduce your problem. For that I would need the complete set up, so could you please provide your topography, or an example, where no explicit topography file is required.
> > 
> > I have some questions:
> > 1. What’s the rational behind (data.shelfic)
> > SHELFICElatentHeat = 9.0e99,
> > SHELFICEHeatCapacity_Cp=9.0e99,
> > ? This will most likely make the model explode
> > Same with ICEFRONTlatentHeat = 9.0e99, (data.icefront)
> > 2. file:
> > SHELFICEloadAnomalyFile = ‘potential.jmd95z'
> > does not contain the correct field. From your gendata.m is looks like you store a 3d field that’s something like an intergral expression of the hydrostatic pressure balance with constant density 998. The model will only read a 2D field so the first layer of your 3d field = .5*dz*9.81*998. That’s not what you want, because that can also break the model due to strong initial shocks (it depends on your ice shelf depth). YOu nee the anomaly of phyHyd which is the vertical integral of gravity*(rho-rho0) at (rather, close to) the bottom of the iceshelf. See isomip/input/gendata.m for an example, how to compute it so that it is completely consistent with the initial hydrograph and EOS (only required if you want to avoid any initial adjustment shocks); the first guess can actually be 0 (that’s definitely better than what you have right now).
> > 
> > I cannot tell you much about the icefront package, I have never used it, but I would “fix” the above things first (or give me a reproducable configuration).
> > 
> > Martin
> > 
> > 
> > On Feb 18, 2014, at 8:24 PM, Harry A <gmtonys3 at hotmail.com> wrote:
> > 
> > > Hi Martin
> > > 
> > > Attached is my gendata.m code which comes up with the topology, which is actually read from a CSV. The size of the simulation is small, 120 x 48 in the xy plane with grid size of 50m, and 10 intervals of 15.5 m in the z direction.
> > > 
> > > Also attached are the major data files. Other files such as the initial temperature and surface heat flux have been tested a lot and I believe are not the issue here. It is something in the icefront package which means it works OK when I turn it off, or if I turn it on and set the icelength input file to zero everywhere, but if I make this non-zero, as is in the gendata.m attached, it breaks down.
> > > 
> > > I've also not got any initial results to pull a pressure distribution from, so as my lake is very slowly moving, I have just used pressure = rho*g*h.
> > > In Size.h if I replace Nr everywhere with Nz then the model returns no output. (Or at least, none of the T.000000050.001.001.data files for T, U V W ....) but from looking at tutorial_Deep_convection I thought it would be ok to leave this as Nr seeing as it is done in that tutorial, which uses a similarly sized grid and Cartesian coordinates. 
> > > 
> > > Thank you so much for any assistance!
> > > 
> > > 
> > > Harry
> > > 
> > > Here is the contents of the gendata.m file attached:
> > > 
> > > close all
> > > ieee='b';
> > > accuracy='real*8';
> > > 
> > > nx=120;
> > > ny=48;
> > > nz=10;
> > > dz=15.5;
> > > dy=50;
> > > dx=dy;
> > > x=csvread('120by48topology.csv'); %This topology is between 0 and 150 meters, with zero's all around the outside so all boundaries are fixed.
> > > 
> > > %So making our topology...
> > > h=-1*x;
> > > h=2*round(h/2); %Making h multiples of 2 rather than rational numbers.
> > > h=h';
> > > %And the overall topology
> > > fid=fopen('topog.box','w',ieee); fwrite(fid,h,accuracy); fclose(fid);
> > > 
> > > %Now, making our shelfice input file
> > > icet= abs(h)>0;
> > > icetopo=icet*0.0; 
> > > icetopo(1:55,1:48)=1*h(1:55,1:48); %So now we make the shelf ice the same depth as the topology over the top 55 spaces in domain
> > > fid=fopen('icetopo.exp1','w','b'); fwrite(fid,icetopo,accuracy);fclose(fid);
> > > 
> > > 
> > > %so, now we generate our icefront.We just set it to be equal to the topology height at the front of the shelfice file, as is done in isomip
> > > iceheight=0*h;
> > > iceheight(56,1:48)=-1*h(55,1:48);%%%%%%%%% 
> > > fid=fopen('ifront.xuyun','w','b'); fwrite(fid,iceheight,accuracy); fclose(fid);
> > > 
> > > %Now making the icelength file
> > > icelength=iceheight*0.0+(dz/dy)*(iceheight~=0); %So, it will only have data in the same places that the iceheight does
> > > fid=fopen('ilength.xuyun','w','b'); fwrite(fid,icelength,accuracy); fclose(fid);
> > > 
> > > 
> > > 
> > > %Now making the potential distribution using simple pressure = rho*9.81*height
> > > potential=ones(nx,ny,nz);
> > > for k=1:nz
> > > potential(:,:,k)=9.81*dz*998*potential(:,:,k)*(k-0.5);
> > > end
> > > fid=fopen('potential.jmd95z','w','b'); fwrite(fid,potential,accuracy);fclose(fid);
> > > 
> > > 
> > > 
> > > 
> > > 
> > > 
> > > > From: Martin.Losch at awi.de
> > > > Date: Tue, 18 Feb 2014 13:17:27 +0100
> > > > To: mitgcm-support at mitgcm.org
> > > > Subject: Re: [MITgcm-support] Icefront returnin NaN based on topology.
> > > > 
> > > > Harry,
> > > > 
> > > > I think it would be useful to know about the details of the configuration. It would also help to have the illustrative example that you describe as a “code” and “input” directory (I am assuming it’s a very small domain), including the matlab (?) file to generate the different input files.
> > > > 
> > > > Martin
> > > > 
> > > > On Feb 16, 2014, at 2:41 AM, Harry A <gmtonys3 at hotmail.com> wrote:
> > > > 
> > > > > Hi Everyone
> > > > > 
> > > > > I'm working on a model including a glacier, so I've been implementing the Icefront package to add this functionality.
> > > > > However, when I am setting up the topology of my submerged ice, I've had NaN results for all outputted velocities and temperatures turning up at strange times.
> > > > > 
> > > > > I have a bottom topology over my whole domain with closed boundaries (Depth = 0 all around) To begin with, I just made a 3x3 grid in the middle of my domain ice, but depending on where I place it the model with either return all NaN or return good results. It seems even if I make this ice shallow (say 10 meters deep in 100m water) or as deep as the water at the certain location, or even deeper, I still get failure with some locations and not others.
> > > > > 
> > > > > To highlight the problem, one example is I can make a single grid point in the middle of my domain have ice 70.4952 meters deep (with the bathymetry file being exactly 70m deep in that location) and it will work, however if I make this 70.4953 or 70.4951 the model returns all NaNs.
> > > > > 
> > > > > Does anyone have any idea what could be causing these issues? Am I correct in saying the Ice depth file for Icefront ( ICEFRONTdepthFile) meant to be positive values ranging from 0 up to the depth of water at any given grid point ? And the ICEFRONTlengthFile should be the facial area of a grid element in the horizontal direction, so If I'm using dX=50, dY=50, dZ = 15, it should be 15*50 = 750?
> > > > > 
> > > > > Thanks for any assistance, I've been running through matlab scrips generating ice topologies for 6 hours now and there seems to be almost no rhyme or reason as to why some situations fail and some work just fine...
> > > > > 
> > > > > Harry
> > > > > _______________________________________________
> > > > > MITgcm-support mailing list
> > > > > MITgcm-support at mitgcm.org
> > > > > http://mitgcm.org/mailman/listinfo/mitgcm-support
> > > > 
> > > > 
> > > > _______________________________________________
> > > > MITgcm-support mailing list
> > > > MITgcm-support at mitgcm.org
> > > > http://mitgcm.org/mailman/listinfo/mitgcm-support
> > > <data><data.icefront><data.pkg><data.shelfice><CPP_OPTIONS.h><packages.conf><SIZE.h><gendata.m>_______________________________________________
> > > MITgcm-support mailing list
> > > MITgcm-support at mitgcm.org
> > > http://mitgcm.org/mailman/listinfo/mitgcm-support
> > 
> > 
> > _______________________________________________
> > MITgcm-support mailing list
> > MITgcm-support at mitgcm.org
> > http://mitgcm.org/mailman/listinfo/mitgcm-support
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-support

More information about the MITgcm-support mailing list