AUTO/input/rwpu.m

clear all
close all
clc

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

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

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

fid = fopen('shelftopo.OLDROUND.bin','r','b');
shelftopoold=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]);

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

threshold=5;
[etax,etay] =find(oldetan>threshold);

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

for i = 1:size(etax,1)
    etaz(i)=((shelftopoold(etax(i),etay(i)))/10*-1)+1;
end
     
OLDX=etax;
OLDY=etay;
OLDZ=etaz;
NEWZ=etaz-1;

shelftoponew=shelftopoold;
pignew=pigold;
shelftoponew(etax,etay)=shelftopoold(etax,etay)+delz;
pignew(etax,etay)=pigold(etax,etay)+delz;

newetan=oldetan;
newetan(etax,etay)=oldetan(etax,etay)-delz;

newdetahdt=olddetahdt;
%newetah=oldetah;
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.2050;
S_bottom = 34.6967;
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(pignew(ix,iy))));
    if isempty(k)
    k=0;
    end
    if k>0
    
      dr = -zg(k) - pignew(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;


for i=1:size(OLDX)
    newsalt(OLDX(i),OLDY(i),NEWZ(i))=oldsalt(OLDX(i),OLDY(i),OLDZ(i));
    newtheta(OLDX(i),OLDY(i),NEWZ(i))=oldtheta(OLDX(i),OLDY(i),OLDZ(i));
  
    newuvel(OLDX(i),OLDY(i),NEWZ(i))=olduvel(OLDX(i),OLDY(i),OLDZ(i));
    %newuvel(OLDX(i)+1,OLDY(i),NEWZ(i))=olduvel(OLDX(i)+1,OLDY(i),OLDZ(i));
    newvvel(OLDX(i),OLDY(i),NEWZ(i))=oldvvel(OLDX(i),OLDY(i),OLDZ(i));    
    newvvel(OLDX(i),OLDY(i)+1,NEWZ(i))=oldvvel(OLDX(i),OLDY(i)+1,OLDZ(i)); 
    
    newgunm1(OLDX(i),OLDY(i),NEWZ(i))=0;
    %newgunm1(OLDX(i)+1,OLDY(i),NEWZ(i))=0;
    newgvnm1(OLDX(i),OLDY(i),NEWZ(i))=0;
    newgvnm1(OLDX(i),OLDY(i)+1,NEWZ(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';
  
   shelftoponew=reshape(shelftoponew,[nx*ny 1]);
   ploadnew=reshape(ploadnew,[nx*ny 1]);


 fid = fopen('pickup.0000025920.data','w','b'); fwrite(fid,pickupnew,'real*8'); 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,shelftoponew,'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			pickupold = fread(fid,inf,'real*8');
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			pigold=fread(fid,inf,'real*8');
17			fclose(fid);
18
19			fid = fopen('shelftopo.OLDROUND.bin','r','b');
20			shelftopoold=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			shelftopoold=reshape(shelftopoold,[nx ny]);
49			ploadold=reshape(ploadold,[nx ny]);
50			pigold=reshape(pigold,[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			threshold=5;
56			[etax,etay] =find(oldetan>threshold);
57
58			%Change single, not all cells >threshold
59			%etax=etax(3);
60			%etay=etay(3);
61
62			for i = 1:size(etax,1)
63			etaz(i)=((shelftopoold(etax(i),etay(i)))/10*-1)+1;
64			end
65
66			OLDX=etax;
67			OLDY=etay;
68			OLDZ=etaz;
69			NEWZ=etaz-1;
70
71			shelftoponew=shelftopoold;
72			pignew=pigold;
73			shelftoponew(etax,etay)=shelftopoold(etax,etay)+delz;
74			pignew(etax,etay)=pigold(etax,etay)+delz;
75
76			newetan=oldetan;
77			newetan(etax,etay)=oldetan(etax,etay)-delz;
78
79			newdetahdt=olddetahdt;
80			%newetah=oldetah;
81			newetah=newetan;
82
83			ploadnew=ploadold;
84			rhoConst = 1030;
85			talpha = 2e-4;
86			sbeta = 7.4e-4;
87			gravity = 9.81;
88			k=1;
89			dzm = abs([zg(1)-zc(1) .5*diff(zc)]);
90			dzp = abs([.5*diff(zc) zc(end)-zg(end)]);
91			p = abs(zc)gravityrhoConst*1e-4;
92			dp = p;
93			kp = 0;
94			eos = 'jmd95z';
95			T_sfc = -1.9;
96			T_bottom = 2;
97			del_T = (T_bottom - T_sfc)/(59*delz);
98
99			for iz = 1:nz;
100
101
102			tref(iz) = T_sfc + del_T((iz-30)delz);
103			if iz<=30;
104			tref(iz)=-1.9;
105			end
106			if iz>=90
107			tref(iz) =2;
108			end
109			end
110
111			S_sfc = 34.2050;
112			S_bottom = 34.6967;
113			del_S = (S_bottom - S_sfc)/(59*delz);
114
115			for iz = 1:nz;
116
117
118			sref(iz) = S_sfc + del_S((iz-30)delz);
119			if iz<=30;
120			sref(iz)=34.2;
121			end
122			if iz>=90
123			sref(iz) =34.7;
124			end
125			end
126
127			t = tref;
128
129			s = sref;
130
131
132			while rms(dp) > 1e-13
133			phiHydF(k) = 0;
134			p0 = p;
135			kp = kp+1;
136			for k = 1:nz
137			switch eos
138			case 'linear'
139			drho = rhoConst(1-talpha(t(k)-tref(k))+sbeta*(s(k)-sref(k)))-rhoConst;
140			case 'jmd95z'
141			drho = densjmd95(s(k),t(k),p(k))-rhoConst;
142			case 'mdjwf'
143			drho = densmdjwf(s(k),t(k),p(k))-rhoConst;
144			otherwise
145			error(sprintf('unknown EOS: %s',eos))
146			end
147			phiHydC(k) = phiHydF(k) + dzm(k)gravitydrho/rhoConst;
148			phiHydF(k+1) = phiHydC(k) + dzp(k)gravitydrho/rhoConst;
149			end
150			switch eos
151			case 'mdjwf'
152			p = (gravityrhoConstabs(zc) + phiHydC*rhoConst)/gravity/rhoConst;
153			end
154			dp = p-p0;
155			end
156
157
158			for ix=1:nx
159			for iy=1:ny
160			k=max(find(abs(zg)<abs(pignew(ix,iy))));
161			if isempty(k)
162			k=0;
163			end
164			if k>0
165
166			dr = -zg(k) - pignew(ix,iy);
167
168			if (dr>=delz/2)
169			ploadnew(ix,iy) = phiHydF(k) + (delz-dr) * (phiHydC(k)-phiHydF(k))/(delz/2);
170			else
171			ploadnew(ix,iy) = phiHydC(k) + (delz/2-dr) * (phiHydF(k+1)-phiHydC(k))/(delz/2);
172			end
173
174			end
175			end
176			end
177
178
179
180
181			newsalt=oldsalt;
182			newtheta=oldtheta;
183			newuvel=olduvel;
184			newvvel=oldvvel;
185			newgunm1=oldgunm1;
186			newgvnm1=oldgvnm1;
187
188
189
190			for i=1:size(OLDX)
191			newsalt(OLDX(i),OLDY(i),NEWZ(i))=oldsalt(OLDX(i),OLDY(i),OLDZ(i));
192			newtheta(OLDX(i),OLDY(i),NEWZ(i))=oldtheta(OLDX(i),OLDY(i),OLDZ(i));
193
194			newuvel(OLDX(i),OLDY(i),NEWZ(i))=olduvel(OLDX(i),OLDY(i),OLDZ(i));
195			%newuvel(OLDX(i)+1,OLDY(i),NEWZ(i))=olduvel(OLDX(i)+1,OLDY(i),OLDZ(i));
196			newvvel(OLDX(i),OLDY(i),NEWZ(i))=oldvvel(OLDX(i),OLDY(i),OLDZ(i));
197			newvvel(OLDX(i),OLDY(i)+1,NEWZ(i))=oldvvel(OLDX(i),OLDY(i)+1,OLDZ(i));
198
199			newgunm1(OLDX(i),OLDY(i),NEWZ(i))=0;
200			%newgunm1(OLDX(i)+1,OLDY(i),NEWZ(i))=0;
201			newgvnm1(OLDX(i),OLDY(i),NEWZ(i))=0;
202			newgvnm1(OLDX(i),OLDY(i)+1,NEWZ(i))=0;
203
204			end
205
206
207
208			% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
209			% % %Write out pickup file
210			uvel=reshape(newuvel,[1 nxnynz]);
211			vvel=reshape(newvvel,[1 nxnynz]);
212			theta=reshape(newtheta,[1 nxnynz]);
213			salt=reshape(newsalt,[1 nxnynz]);
214			gunm1=reshape(newgunm1,[1 nxnynz]);
215			gvnm1=reshape(newgvnm1,[1 nxnynz]);
216			etan=reshape(newetan,[1 nx*ny]);
217			detahdt=reshape(newdetahdt,[1 nx*ny]);
218			etah=reshape(newetah,[1 nx*ny]);
219
220			pickupnew=[uvel,vvel,theta,salt,gunm1,gvnm1,etan,detahdt,etah];
221			pickupnew=pickupnew';
222
223			shelftoponew=reshape(shelftoponew,[nx*ny 1]);
224			ploadnew=reshape(ploadnew,[nx*ny 1]);
225
226
227			fid = fopen('pickup.0000025920.data','w','b'); fwrite(fid,pickupnew,'real*8'); fclose(fid);
228			fid = fopen('pload.pig.jmd95z','w','b'); fwrite(fid,ploadnew,'real*8'); fclose(fid)
229			fid = fopen('shelftopo.round.bin','w','b'); fwrite(fid,shelftoponew,'real*8'); fclose(fid)
230
231
232
233			%
234
235