/[MITgcm]/MITgcm_contrib/shelfice_remeshing/input/gendata_vero.m
ViewVC logotype

Contents of /MITgcm_contrib/shelfice_remeshing/input/gendata_vero.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.2 - (show annotations) (download)
Thu Sep 10 14:41:58 2015 UTC (9 years, 10 months ago) by dgoldberg
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +0 -0 lines
FILE REMOVED
*** empty log message ***

1 %Verion of gendata.m modified by Vero
2 %This is a matlab script that generates the input data
3
4
5 % the configuation approximately the ISOMIP experiment no. 1
6 % require matlab functions for equation of state
7
8
9 % Dimensions of grid
10 nx=1;
11 ny=200;
12 nz=100;
13 delz = 10;
14
15 hfacMin = 0.2;
16
17 dlat = 0.125/16; dy=dlat;
18 dlon = 0.125/4; dx=dlon;
19
20 %eos = 'linear';
21 eos = 'jmd95z';
22 % eos = 'mdjwf';
23
24 acc = 'real*8';
25
26 long = [-105.5:dlon:-105.5];
27 lonc = long+dlon/2;
28 latg = [-75.4457:dlat:-73.8809-dlat];
29 latc = latg+dlat/2;
30 size(latc)
31
32 % Nominal depth of model (meters)
33 H = -1000; %water depth in the ice shelf cavity
34 Hmin = -900; % deepest point of cavern
35 Hmax = -300; % shallowest point of cavern
36 dHdy = (Hmax-Hmin)/(max(latc)-min(latc)); %Slope of ice shelf: if denominator = nx, shelf will cover the whole domain
37
38 bathy = ones(nx,ny)*H; %For flat bathymetry: bathy = ones(nx,ny)*H;
39 %bathy(1,:) = 0;
40 %bathy(2,:) = 0;
41 %bathy(nx,:) = 0;
42 %bathy(nx-1,:) = 0;
43 bathy(:,1) = 0;
44 %bathy(:,ny) = 0;
45 fid=fopen('bathymetry.pig.bin','w','b'); fwrite(fid,bathy,acc);fclose(fid);
46
47
48 dz = delz*ones(1,nz);
49 zgp1 = [0,cumsum(dz)];
50 zc = .5*(zgp1(1:end-1)+zgp1(2:end));
51 zg = zgp1(1:end-1);
52 dz = diff(zgp1);
53 sprintf('delZ = %d * %7.6g,',nz,dz)
54
55
56 T_sfc = -1.9;
57 T_bottom = 2;
58 del_T = (T_bottom - T_sfc)/(59*delz);
59
60 for iz = 1:nz;
61
62
63 tref(iz) = T_sfc + del_T*((iz-30)*delz);
64 if iz<=30;
65 tref(iz)=-1.9;
66 end
67 if iz>=90
68 tref(iz) =2;
69 end
70 end
71
72 S_sfc = 34.2;
73 S_bottom = 34.7;
74 del_S = (S_bottom - S_sfc)/(59*delz);
75
76 for iz = 1:nz;
77
78
79 sref(iz) = S_sfc + del_S*((iz-30)*delz);
80 if iz<=30;
81 sref(iz)=34.2;
82 end
83 if iz>=90
84 sref(iz) =34.7;
85 end
86 end
87
88 % Gravity
89 gravity=9.81;
90 rhoConst = 1030;
91 % compute potential field underneath ice shelf
92 talpha = 2e-4;
93 sbeta = 7.4e-4;
94 % tref = -1.9*ones(nz,1);
95 t = tref;
96 % sref = 34.4*ones(nz,1);
97 s = sref;
98 gravity = 9.81;
99 k=1;
100 dzm = abs([zg(1)-zc(1) .5*diff(zc)]);
101 dzp = abs([.5*diff(zc) zc(end)-zg(end)]);
102 p = abs(zc)*gravity*rhoConst*1e-4;
103 dp = p;
104 kp = 0;
105
106
107
108 while rms(dp) > 1e-13
109 phiHydF(k) = 0;
110 p0 = p;
111 kp = kp+1
112 for k = 1:nz
113 switch eos
114 case 'linear'
115 drho = rhoConst*(1-talpha*(t(k)-tref(k))+sbeta*(s(k)-sref(k)))-rhoConst;
116 case 'jmd95z'
117 drho = densjmd95(s(k),t(k),p(k))-rhoConst;
118 case 'mdjwf'
119 drho = densmdjwf(s(k),t(k),p(k))-rhoConst;
120 otherwise
121 error(sprintf('unknown EOS: %s',eos))
122 end
123 phiHydC(k) = phiHydF(k) + dzm(k)*gravity*drho/rhoConst;
124 phiHydF(k+1) = phiHydC(k) + dzp(k)*gravity*drho/rhoConst;
125 end
126 switch eos
127 case 'mdjwf'
128 p = (gravity*rhoConst*abs(zc) + phiHydC*rhoConst)/gravity/rhoConst;
129 end
130 dp = p-p0;
131 end
132
133
134 %Modify icetopo (shape of ice shelf cavity)
135
136 %icetopo = ones(ny,1)*min(Hmin + 2*dHdy*(lonc(nx)-long),Hmax);
137 B=min(Hmin + 2*dHdy*(latc(ny)-latg),Hmax);
138 B=fliplr(B);
139
140 icetopo = ones(nx,1)*B;
141 %icetopo = icetopo';
142 icetopo(:,101:end)=0;
143
144
145
146
147 % adjust topo so that no hfac is smaller than hfacMin
148
149 for ix=1:nx
150 for iy=1:ny
151 k=max(find(abs(zg)<abs(icetopo(ix,iy))));
152
153 if(~isempty(k))
154
155
156 hfacTemp = (icetopo(ix,iy) - (-zg(k+1)))/delz;
157
158 if (hfacTemp < hfacMin)
159 if (hfacTemp < hfacMin/2)
160 hfacTemp = 0;
161 else
162 hfacTemp = hfacMin;
163 end
164 end
165
166 else
167
168 hfacTemp = 0;
169
170 end
171
172 icetopo(ix,iy) = icetopo(ix,iy) + hfacTemp;
173
174 end
175 end
176
177 % phi anomaly (relative to hydrostatic with rho_const) at icetopo
178
179 phi0surf = zeros(nx,ny);
180
181 for ix=1:nx
182 for iy=1:ny
183 k=max(find(abs(zg)<abs(icetopo(ix,iy))));
184 if isempty(k)
185 k=0;
186 end
187 if k>0
188
189 dr = -zg(k) - icetopo(ix,iy);
190
191 if (dr>=delz/2)
192 phi0surf(ix,iy) = phiHydF(k) + (delz-dr) * (phiHydC(k)-phiHydF(k))/(delz/2);
193 else
194 phi0surf(ix,iy) = phiHydC(k) + (delz/2-dr) * (phiHydF(k+1)-phiHydC(k))/(delz/2);
195 end
196
197 end
198 end
199 end
200
201 mass = phi0surf / gravity - rhoConst * icetopo;
202
203 fid = fopen('shelftopo.pig.bin','w','b'); fwrite(fid,icetopo,'real*8'); fclose(fid);
204 fid = fopen('shelficemassinit.bin','w','b'); fwrite(fid,mass,'real*8'); fclose(fid);
205 fid = fopen('pload.pig.jmd95z','w','b'); fwrite(fid,phi0surf,'real*8'); fclose(fid);
206
207 etainit = zeros(size(phi0surf));
208
209 % new topography: icetopo rounded to the nearest k * deltaZ
210 % eta_init set to make difference
211
212 icetopo2 = icetopo;
213
214 for ix=1:nx
215 for iy=1:ny
216 k=max(find(abs(zg)<abs(icetopo2(ix,iy))));
217 if isempty(k)
218 k=0;
219 else
220
221 dr = 1-(-zg(k) - icetopo2(ix,iy))/delz;
222 if (dr > .5)
223 % bring Ro_surf *up* to closest grid face & make etainit negative
224 % to compensate
225 icetopo2(ix,iy) = -zg(k);
226 etainit(ix,iy) = (dr-1)*delz;
227 else
228 % bring Ro_surf *down* to closest grid face & make etainit pos
229 % to compensate
230 icetopo2(ix,iy) = -zg(k+1);
231 etainit(ix,iy) = (dr)*delz;
232 end
233
234 end
235 end
236 end
237
238 fid = fopen('shelftopo.round.bin','w','b'); fwrite(fid,icetopo2,'real*8'); fclose(fid);
239 fid = fopen('etainit.round.bin','w','b'); fwrite(fid,etainit,'real*8'); fclose(fid);

  ViewVC Help
Powered by ViewVC 1.1.22