MANUAL/input/Remeshing.m

clear all
close all
clc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read in and reshape files
fid = fopen('pickup.OLD.data','r','b');
pickupold = fread(fid,inf,'float64');
fclose(fid);

fid = fopen('pload.OLD.jmd95z','r','b');
ploadold = fread(fid,inf,'real*8');
fclose(fid);

fid = fopen('shelftopo.OLDPIG.bin','r','b');
stpigold=fread(fid,inf,'real*8');
fclose(fid);

fid = fopen('shelftopo.OLDROUND.bin','r','b');
stroundold=fread(fid,inf,'real*8');
fclose(fid);

nx=1;
ny=200;
nz=100;

delz=10;

dz = delz*ones(1,nz);
zgp1 = [0,cumsum(dz)];
zc = .5*(zgp1(1:end-1)+zgp1(2:end));
zg = zgp1(1:end-1);
dz = diff(zgp1);

s3=nx*ny*nz;
s2=nx*ny;

olduvel=reshape(pickupold(1:s3),[nx ny nz]);
oldvvel=reshape(pickupold(1+(1*s3):2*s3),[nx ny nz]);
oldtheta=reshape(pickupold(1+2*s3:3*s3),[nx ny nz]);
oldsalt=reshape(pickupold(1+3*s3:4*s3),[nx ny nz]);
oldgunm1=reshape(pickupold(1+4*s3:5*s3),[nx ny nz]);
oldgvnm1=reshape(pickupold(1+5*s3:6*s3),[nx ny nz]);
oldetan=reshape(pickupold(1+6*s3:(6*s3)+(1*s2)),[nx ny]);
olddetahdt=reshape(pickupold(1+6*s3+(1*s2):(6*s3)+(2*s2)),[nx ny]);
oldetah=reshape(pickupold(1+6*s3+(2*s2):(6*s3)+(3*s2)),[nx ny]);

stroundold=reshape(stroundold,[nx ny]);
ploadold=reshape(ploadold,[nx ny]);
stpigold=reshape(stpigold,[nx ny]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Find where eta indicates a change in ice shelf topo above threshold.
%Change eta, iceshelf topo and r0surfto reflect change in mesh

splitthreshold=5;
mergethreshold=-5.6;


[splitx,splity] =find(oldetan>splitthreshold & stroundold<0);
[mergex,mergey] =find(oldetan<mergethreshold & stroundold<0);

%[splitx,splity] =find(oldetan>splitthreshold);
%[mergex,mergey] =find(oldetan>mergethreshold);

%Change single, not all cells >threshold
%etax=etax(3);
%etay=etay(3);


for i = 1:length(splitx)
    splitz(i)=((stroundold(splitx(i),splity(i)))/10*-1)+1;
end

for i = 1:length(mergex)
    mergez(i)=((stroundold(mergex(i),mergey(i)))/10*-1)+1;
end


OSPLITX=splitx;
OSPLITY=splity;
OSPLITZ=splitz;
NSPLITZ=splitz-1;

OMERGEX=mergex;
OMERGEY=mergey;
OMERGEZ=mergez;
NMERGEZ=mergez+1;

stroundnew=stroundold;
stpignew=stpigold;
stroundnew(splitx,splity)=stroundold(splitx,splity)+delz;
stpignew(splitx,splity)=stpigold(splitx,splity)+delz;
stroundnew(mergex,mergey)=stroundold(mergex,mergey)-delz;
stpignew(mergex,mergey)=stpigold(mergex,mergey)-delz;

newetan=oldetan;
newetan(splitx,splity)=oldetan(splitx,splity)-delz;
newetan(mergex,mergey)=oldetan(mergex,mergey)+delz;

newdetahdt=olddetahdt;
%%%%%%%%%%%%%%%%%%%%%%
%Etah=etan for now
newetah=newetan;

ploadnew=ploadold;
rhoConst = 1030;
talpha = 2e-4;
sbeta  = 7.4e-4;
gravity = 9.81;
k=1;
dzm = abs([zg(1)-zc(1) .5*diff(zc)]);
dzp = abs([.5*diff(zc) zc(end)-zg(end)]);
p = abs(zc)*gravity*rhoConst*1e-4;
dp = p;
kp = 0;
eos = 'jmd95z';
T_sfc = -1.9;
T_bottom = 2;
del_T = (T_bottom - T_sfc)/(59*delz);

for iz = 1:nz;
    tref(iz) = T_sfc + del_T*((iz-30)*delz);
    if iz<=30;
        tref(iz)=-1.9;
    end
    if iz>=90
        tref(iz) =2;
    end
end

S_sfc = 34.2;
S_bottom = 34.7;
del_S = (S_bottom - S_sfc)/(59*delz);

for iz = 1:nz;
    
    
    sref(iz) = S_sfc + del_S*((iz-30)*delz);
    if iz<=30;
        sref(iz)=34.2;
    end
    if iz>=90
        sref(iz) =34.7;
    end
end

t    = tref;

s    = sref;


while rms(dp) > 1e-13
    phiHydF(k) = 0;
    p0 = p;
    kp = kp+1;
    for k = 1:nz
        switch eos
            case 'linear'
                drho = rhoConst*(1-talpha*(t(k)-tref(k))+sbeta*(s(k)-sref(k)))-rhoConst;
            case 'jmd95z'
                drho = densjmd95(s(k),t(k),p(k))-rhoConst;
            case 'mdjwf'
                drho = densmdjwf(s(k),t(k),p(k))-rhoConst;
            otherwise
                error(sprintf('unknown EOS: %s',eos))
        end
        phiHydC(k)   = phiHydF(k) + dzm(k)*gravity*drho/rhoConst;
        phiHydF(k+1) = phiHydC(k) + dzp(k)*gravity*drho/rhoConst;
    end
    switch eos
        case 'mdjwf'
            p = (gravity*rhoConst*abs(zc) + phiHydC*rhoConst)/gravity/rhoConst;
    end
    dp = p-p0;
end


for ix=1:nx
    for iy=1:ny
        k=max(find(abs(zg)<abs(stpignew(ix,iy))));
        if isempty(k)
            k=0;
        end
        if k>0
            
            dr = -zg(k) - stpignew(ix,iy);
            
            if (dr>=delz/2)
                ploadnew(ix,iy) = phiHydF(k) + (delz-dr) * (phiHydC(k)-phiHydF(k))/(delz/2);
            else
                ploadnew(ix,iy) = phiHydC(k) + (delz/2-dr) * (phiHydF(k+1)-phiHydC(k))/(delz/2);
            end
            
        end
    end
end


newsalt=oldsalt;
newtheta=oldtheta;
newuvel=olduvel;
newvvel=oldvvel;
newgunm1=oldgunm1;
newgvnm1=oldgvnm1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Split cells, new cells take values of the cell they are split from

for i=1:length(OSPLITX)
    newsalt(OSPLITX(i),OSPLITY(i),NSPLITZ(i))=oldsalt(OSPLITX(i),OSPLITY(i),OSPLITZ(i));
    newtheta(OSPLITX(i),OSPLITY(i),NSPLITZ(i))=oldtheta(OSPLITX(i),OSPLITY(i),OSPLITZ(i));
    
    newuvel(OSPLITX(i),OSPLITY(i),NSPLITZ(i))=olduvel(OSPLITX(i),OSPLITY(i),OSPLITZ(i));
    %newuvel(OLDX(i)+1,OLDY(i),NEWZ(i))=olduvel(OLDX(i)+1,OLDY(i),OLDZ(i));
    newvvel(OSPLITX(i),OSPLITY(i),NSPLITZ(i))=oldvvel(OSPLITX(i),OSPLITY(i),OSPLITZ(i));
    newvvel(OSPLITX(i),OSPLITY(i)+1,NSPLITZ(i))=oldvvel(OSPLITX(i),OSPLITY(i)+1,OSPLITZ(i));
    
    newgunm1(OSPLITX(i),OSPLITY(i),NSPLITZ(i))=0;
    %newgunm1(OLDX(i)+1,OLDY(i),NEWZ(i))=0;
    newgvnm1(OSPLITX(i),OSPLITY(i),NSPLITZ(i))=0;
    newgvnm1(OSPLITX(i),OSPLITY(i)+1,NSPLITZ(i))=0;
    
end

for i=1:length(OMERGEX)
    newsalt(OMERGEX(i),OMERGEY(i),NMERGEZ(i))=((oldsalt(OMERGEX(i),OMERGEY(i),OMERGEZ(i))*(dz+oldetan(OMERGEX(i),OMERGEY(i))))...
        +(dz*oldsalt(OMERGEX(i),OMERGEY(i),NMERGEZ(i))))/(2*dz+oldetan(OMERGEX(i),OMERGEY(i)));
    newtheta(OMERGEX(i),OMERGEY(i),NMERGEZ(i))=((oldtheta(OMERGEX(i),OMERGEY(i),OMERGEZ(i))*(dz+oldetan(OMERGEX(i),OMERGEY(i))))...
        +(dz*oldtheta(OMERGEX(i),OMERGEY(i),NMERGEZ(i))))/(2*dz+oldetan(OMERGEX(i),OMERGEY(i)));
    
     newsalt(OMERGEX(i),OMERGEY(i),OMERGEZ(i))=0;
    newtheta(OMERGEX(i),OMERGEY(i),OMERGEZ(i))=0;
    
    newuvel(OMERGEX(i),OMERGEY(i),NMERGEZ(i))=((olduvel(OMERGEX(i),OMERGEY(i),OMERGEZ(i))*(dz+oldetan(OMERGEX(i),OMERGEY(i))))...
        +(dz*olduvel(OMERGEX(i),OMERGEY(i),NMERGEZ(i))))/(2*dz+oldetan(OMERGEX(i),OMERGEY(i)));
    %newuvel(OMERGEX(i)+1,OMERGEY(i),NMERGEZ(i))=((olduvel(OMERGEX(i)+1,OMERGEY(i),OMERGEZ(i))*(dz+oldetan(OMERGEX(i)+1,OMERGEY(i))))...
    %    +(dz*olduvel(OMERGEX(i)+1,OMERGEY(i),NMERGEZ(i))))/(2*dz+oldetan(OMERGEX(i)+1,OMERGEY(i)));
    newvvel(OMERGEX(i),OMERGEY(i),NMERGEZ(i))=((oldvvel(OMERGEX(i),OMERGEY(i),OMERGEZ(i))*(dz+oldetan(OMERGEX(i),OMERGEY(i))))...
        +(dz*oldvvel(OMERGEX(i),OMERGEY(i),NMERGEZ(i))))/(2*dz+oldetan(OMERGEX(i),OMERGEY(i)));
    newvvel(OMERGEX(i),OMERGEY(i)+1,NMERGEZ(i))=((oldvvel(OMERGEX(i),OMERGEY(i)+1,OMERGEZ(i))*(dz+oldetan(OMERGEX(i),OMERGEY(i))))...
        +(dz*oldvvel(OMERGEX(i),OMERGEY(i)+1,NMERGEZ(i))))/(2*dz+oldetan(OMERGEX(i),OMERGEY(i)));
    
    newgunm1(OMERGEX(i),OMERGEY(i),NMERGEZ(i))=0;
    %newgunm1(OMERGEX(i)+1,OMERGEY(i),NMERGEZ(i))=0;
    newgvnm1(OMERGEX(i),OMERGEY(i),NMERGEZ(i))=0;
    newgvnm1(OMERGEX(i),OMERGEY(i)+1,NMERGEZ(i))=0;
    
end


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % %Write out pickup file
uvel=reshape(newuvel,[1 nx*ny*nz]);
vvel=reshape(newvvel,[1 nx*ny*nz]);
theta=reshape(newtheta,[1 nx*ny*nz]);
salt=reshape(newsalt,[1 nx*ny*nz]);
gunm1=reshape(newgunm1,[1 nx*ny*nz]);
gvnm1=reshape(newgvnm1,[1 nx*ny*nz]);
etan=reshape(newetan,[1 nx*ny]);
detahdt=reshape(newdetahdt,[1 nx*ny]);
etah=reshape(newetah,[1 nx*ny]);

pickupnew=[uvel,vvel,theta,salt,gunm1,gvnm1,etan,detahdt,etah];
pickupnew=pickupnew';

stroundnew=reshape(stroundnew,[nx*ny 1]);
ploadnew=reshape(ploadnew,[nx*ny 1]);


fid = fopen('pickup.0000004320.data','w','b'); fwrite(fid,pickupnew,'float64'); fclose(fid);
fid = fopen('pload.pig.jmd95z','w','b'); fwrite(fid,ploadnew,'real*8'); fclose(fid)
fid = fopen('shelftopo.round.bin','w','b'); fwrite(fid,stroundnew,'real*8'); fclose(fid)


1	dgoldberg	1.1	clear all
2			close all
3			clc
4
5			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6			% Read in and reshape files
7			fid = fopen('pickup.OLD.data','r','b');
8	dgoldberg	1.2	pickupold = fread(fid,inf,'float64');
9	dgoldberg	1.1	fclose(fid);
10
11			fid = fopen('pload.OLD.jmd95z','r','b');
12			ploadold = fread(fid,inf,'real*8');
13			fclose(fid);
14
15			fid = fopen('shelftopo.OLDPIG.bin','r','b');
16			stpigold=fread(fid,inf,'real*8');
17			fclose(fid);
18
19			fid = fopen('shelftopo.OLDROUND.bin','r','b');
20			stroundold=fread(fid,inf,'real*8');
21			fclose(fid);
22
23			nx=1;
24			ny=200;
25			nz=100;
26
27			delz=10;
28
29			dz = delz*ones(1,nz);
30			zgp1 = [0,cumsum(dz)];
31			zc = .5*(zgp1(1:end-1)+zgp1(2:end));
32			zg = zgp1(1:end-1);
33			dz = diff(zgp1);
34
35			s3=nxnynz;
36			s2=nx*ny;
37
38			olduvel=reshape(pickupold(1:s3),[nx ny nz]);
39			oldvvel=reshape(pickupold(1+(1s3):2s3),[nx ny nz]);
40			oldtheta=reshape(pickupold(1+2s3:3s3),[nx ny nz]);
41			oldsalt=reshape(pickupold(1+3s3:4s3),[nx ny nz]);
42			oldgunm1=reshape(pickupold(1+4s3:5s3),[nx ny nz]);
43			oldgvnm1=reshape(pickupold(1+5s3:6s3),[nx ny nz]);
44			oldetan=reshape(pickupold(1+6s3:(6s3)+(1*s2)),[nx ny]);
45			olddetahdt=reshape(pickupold(1+6s3+(1s2):(6s3)+(2s2)),[nx ny]);
46			oldetah=reshape(pickupold(1+6s3+(2s2):(6s3)+(3s2)),[nx ny]);
47
48			stroundold=reshape(stroundold,[nx ny]);
49			ploadold=reshape(ploadold,[nx ny]);
50			stpigold=reshape(stpigold,[nx ny]);
51			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
52			%Find where eta indicates a change in ice shelf topo above threshold.
53			%Change eta, iceshelf topo and r0surfto reflect change in mesh
54
55			splitthreshold=5;
56	dgoldberg	1.2	mergethreshold=-5.6;
57	dgoldberg	1.1
58
59			[splitx,splity] =find(oldetan>splitthreshold & stroundold<0);
60			[mergex,mergey] =find(oldetan<mergethreshold & stroundold<0);
61
62			%[splitx,splity] =find(oldetan>splitthreshold);
63			%[mergex,mergey] =find(oldetan>mergethreshold);
64
65			%Change single, not all cells >threshold
66			%etax=etax(3);
67			%etay=etay(3);
68
69
70
71			for i = 1:length(splitx)
72			splitz(i)=((stroundold(splitx(i),splity(i)))/10*-1)+1;
73			end
74
75			for i = 1:length(mergex)
76			mergez(i)=((stroundold(mergex(i),mergey(i)))/10*-1)+1;
77			end
78
79
80			OSPLITX=splitx;
81			OSPLITY=splity;
82			OSPLITZ=splitz;
83			NSPLITZ=splitz-1;
84
85			OMERGEX=mergex;
86			OMERGEY=mergey;
87			OMERGEZ=mergez;
88			NMERGEZ=mergez+1;
89
90			stroundnew=stroundold;
91			stpignew=stpigold;
92			stroundnew(splitx,splity)=stroundold(splitx,splity)+delz;
93			stpignew(splitx,splity)=stpigold(splitx,splity)+delz;
94			stroundnew(mergex,mergey)=stroundold(mergex,mergey)-delz;
95			stpignew(mergex,mergey)=stpigold(mergex,mergey)-delz;
96
97			newetan=oldetan;
98			newetan(splitx,splity)=oldetan(splitx,splity)-delz;
99			newetan(mergex,mergey)=oldetan(mergex,mergey)+delz;
100
101			newdetahdt=olddetahdt;
102			%%%%%%%%%%%%%%%%%%%%%%
103			%Etah=etan for now
104			newetah=newetan;
105
106			ploadnew=ploadold;
107			rhoConst = 1030;
108			talpha = 2e-4;
109			sbeta = 7.4e-4;
110			gravity = 9.81;
111			k=1;
112			dzm = abs([zg(1)-zc(1) .5*diff(zc)]);
113			dzp = abs([.5*diff(zc) zc(end)-zg(end)]);
114			p = abs(zc)gravityrhoConst*1e-4;
115			dp = p;
116			kp = 0;
117			eos = 'jmd95z';
118			T_sfc = -1.9;
119			T_bottom = 2;
120			del_T = (T_bottom - T_sfc)/(59*delz);
121
122			for iz = 1:nz;
123			tref(iz) = T_sfc + del_T((iz-30)delz);
124			if iz<=30;
125			tref(iz)=-1.9;
126			end
127			if iz>=90
128			tref(iz) =2;
129			end
130			end
131
132			S_sfc = 34.2;
133			S_bottom = 34.7;
134			del_S = (S_bottom - S_sfc)/(59*delz);
135
136			for iz = 1:nz;
137
138
139			sref(iz) = S_sfc + del_S((iz-30)delz);
140			if iz<=30;
141			sref(iz)=34.2;
142			end
143			if iz>=90
144			sref(iz) =34.7;
145			end
146			end
147
148			t = tref;
149
150			s = sref;
151
152
153			while rms(dp) > 1e-13
154			phiHydF(k) = 0;
155			p0 = p;
156			kp = kp+1;
157			for k = 1:nz
158			switch eos
159			case 'linear'
160			drho = rhoConst(1-talpha(t(k)-tref(k))+sbeta*(s(k)-sref(k)))-rhoConst;
161			case 'jmd95z'
162			drho = densjmd95(s(k),t(k),p(k))-rhoConst;
163			case 'mdjwf'
164			drho = densmdjwf(s(k),t(k),p(k))-rhoConst;
165			otherwise
166			error(sprintf('unknown EOS: %s',eos))
167			end
168			phiHydC(k) = phiHydF(k) + dzm(k)gravitydrho/rhoConst;
169			phiHydF(k+1) = phiHydC(k) + dzp(k)gravitydrho/rhoConst;
170			end
171			switch eos
172			case 'mdjwf'
173			p = (gravityrhoConstabs(zc) + phiHydC*rhoConst)/gravity/rhoConst;
174			end
175			dp = p-p0;
176			end
177
178
179			for ix=1:nx
180			for iy=1:ny
181			k=max(find(abs(zg)<abs(stpignew(ix,iy))));
182			if isempty(k)
183			k=0;
184			end
185			if k>0
186
187			dr = -zg(k) - stpignew(ix,iy);
188
189			if (dr>=delz/2)
190			ploadnew(ix,iy) = phiHydF(k) + (delz-dr) * (phiHydC(k)-phiHydF(k))/(delz/2);
191			else
192			ploadnew(ix,iy) = phiHydC(k) + (delz/2-dr) * (phiHydF(k+1)-phiHydC(k))/(delz/2);
193			end
194
195			end
196			end
197			end
198
199
200
201
202			newsalt=oldsalt;
203			newtheta=oldtheta;
204			newuvel=olduvel;
205			newvvel=oldvvel;
206			newgunm1=oldgunm1;
207			newgvnm1=oldgvnm1;
208
209			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
210			% Split cells, new cells take values of the cell they are split from
211
212			for i=1:length(OSPLITX)
213			newsalt(OSPLITX(i),OSPLITY(i),NSPLITZ(i))=oldsalt(OSPLITX(i),OSPLITY(i),OSPLITZ(i));
214			newtheta(OSPLITX(i),OSPLITY(i),NSPLITZ(i))=oldtheta(OSPLITX(i),OSPLITY(i),OSPLITZ(i));
215
216			newuvel(OSPLITX(i),OSPLITY(i),NSPLITZ(i))=olduvel(OSPLITX(i),OSPLITY(i),OSPLITZ(i));
217			%newuvel(OLDX(i)+1,OLDY(i),NEWZ(i))=olduvel(OLDX(i)+1,OLDY(i),OLDZ(i));
218			newvvel(OSPLITX(i),OSPLITY(i),NSPLITZ(i))=oldvvel(OSPLITX(i),OSPLITY(i),OSPLITZ(i));
219			newvvel(OSPLITX(i),OSPLITY(i)+1,NSPLITZ(i))=oldvvel(OSPLITX(i),OSPLITY(i)+1,OSPLITZ(i));
220
221			newgunm1(OSPLITX(i),OSPLITY(i),NSPLITZ(i))=0;
222			%newgunm1(OLDX(i)+1,OLDY(i),NEWZ(i))=0;
223			newgvnm1(OSPLITX(i),OSPLITY(i),NSPLITZ(i))=0;
224			newgvnm1(OSPLITX(i),OSPLITY(i)+1,NSPLITZ(i))=0;
225
226			end
227
228			for i=1:length(OMERGEX)
229			newsalt(OMERGEX(i),OMERGEY(i),NMERGEZ(i))=((oldsalt(OMERGEX(i),OMERGEY(i),OMERGEZ(i))*(dz+oldetan(OMERGEX(i),OMERGEY(i))))...
230			+(dzoldsalt(OMERGEX(i),OMERGEY(i),NMERGEZ(i))))/(2dz+oldetan(OMERGEX(i),OMERGEY(i)));
231			newtheta(OMERGEX(i),OMERGEY(i),NMERGEZ(i))=((oldtheta(OMERGEX(i),OMERGEY(i),OMERGEZ(i))*(dz+oldetan(OMERGEX(i),OMERGEY(i))))...
232			+(dzoldtheta(OMERGEX(i),OMERGEY(i),NMERGEZ(i))))/(2dz+oldetan(OMERGEX(i),OMERGEY(i)));
233
234	dgoldberg	1.2	newsalt(OMERGEX(i),OMERGEY(i),OMERGEZ(i))=0;
235			newtheta(OMERGEX(i),OMERGEY(i),OMERGEZ(i))=0;
236
237	dgoldberg	1.1	newuvel(OMERGEX(i),OMERGEY(i),NMERGEZ(i))=((olduvel(OMERGEX(i),OMERGEY(i),OMERGEZ(i))*(dz+oldetan(OMERGEX(i),OMERGEY(i))))...
238			+(dzolduvel(OMERGEX(i),OMERGEY(i),NMERGEZ(i))))/(2dz+oldetan(OMERGEX(i),OMERGEY(i)));
239			%newuvel(OMERGEX(i)+1,OMERGEY(i),NMERGEZ(i))=((olduvel(OMERGEX(i)+1,OMERGEY(i),OMERGEZ(i))*(dz+oldetan(OMERGEX(i)+1,OMERGEY(i))))...
240			% +(dzolduvel(OMERGEX(i)+1,OMERGEY(i),NMERGEZ(i))))/(2dz+oldetan(OMERGEX(i)+1,OMERGEY(i)));
241			newvvel(OMERGEX(i),OMERGEY(i),NMERGEZ(i))=((oldvvel(OMERGEX(i),OMERGEY(i),OMERGEZ(i))*(dz+oldetan(OMERGEX(i),OMERGEY(i))))...
242			+(dzoldvvel(OMERGEX(i),OMERGEY(i),NMERGEZ(i))))/(2dz+oldetan(OMERGEX(i),OMERGEY(i)));
243			newvvel(OMERGEX(i),OMERGEY(i)+1,NMERGEZ(i))=((oldvvel(OMERGEX(i),OMERGEY(i)+1,OMERGEZ(i))*(dz+oldetan(OMERGEX(i),OMERGEY(i))))...
244			+(dzoldvvel(OMERGEX(i),OMERGEY(i)+1,NMERGEZ(i))))/(2dz+oldetan(OMERGEX(i),OMERGEY(i)));
245
246			newgunm1(OMERGEX(i),OMERGEY(i),NMERGEZ(i))=0;
247			%newgunm1(OMERGEX(i)+1,OMERGEY(i),NMERGEZ(i))=0;
248			newgvnm1(OMERGEX(i),OMERGEY(i),NMERGEZ(i))=0;
249			newgvnm1(OMERGEX(i),OMERGEY(i)+1,NMERGEZ(i))=0;
250
251			end
252
253
254
255			% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
256			% % %Write out pickup file
257			uvel=reshape(newuvel,[1 nxnynz]);
258			vvel=reshape(newvvel,[1 nxnynz]);
259			theta=reshape(newtheta,[1 nxnynz]);
260			salt=reshape(newsalt,[1 nxnynz]);
261			gunm1=reshape(newgunm1,[1 nxnynz]);
262			gvnm1=reshape(newgvnm1,[1 nxnynz]);
263			etan=reshape(newetan,[1 nx*ny]);
264			detahdt=reshape(newdetahdt,[1 nx*ny]);
265			etah=reshape(newetah,[1 nx*ny]);
266
267			pickupnew=[uvel,vvel,theta,salt,gunm1,gvnm1,etan,detahdt,etah];
268			pickupnew=pickupnew';
269
270			stroundnew=reshape(stroundnew,[nx*ny 1]);
271			ploadnew=reshape(ploadnew,[nx*ny 1]);
272
273
274	dgoldberg	1.2	fid = fopen('pickup.0000004320.data','w','b'); fwrite(fid,pickupnew,'float64'); fclose(fid);
275	dgoldberg	1.1	fid = fopen('pload.pig.jmd95z','w','b'); fwrite(fid,ploadnew,'real*8'); fclose(fid)
276			fid = fopen('shelftopo.round.bin','w','b'); fwrite(fid,stroundnew,'real*8'); fclose(fid)
277
278
279