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	clear all
2	close all
3	clc
4
5	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6	% Read in and reshape files
7	fid = fopen('pickup.OLD.data','r','b');
8	pickupold = fread(fid,inf,'float64');
9	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	mergethreshold=-5.6;
57
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	newsalt(OMERGEX(i),OMERGEY(i),OMERGEZ(i))=0;
235	newtheta(OMERGEX(i),OMERGEY(i),OMERGEZ(i))=0;
236
237	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	fid = fopen('pickup.0000004320.data','w','b'); fwrite(fid,pickupnew,'float64'); fclose(fid);
275	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