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

Wed Feb 19 05:09:00 EST 2014

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