source: issm/oecreview/Archive/23390-24306/ISSM-23877-23878.diff

Last change on this file was 24307, checked in by Mathieu Morlighem, 5 years ago

NEW: adding Archive/23390-24306
File size: 35.2 KB
RevLine 
[24307]1Index: ../trunk-jpl/test/NightlyRun/test4003.m
2===================================================================
3--- ../trunk-jpl/test/NightlyRun/test4003.m (nonexistent)
4+++ ../trunk-jpl/test/NightlyRun/test4003.m (revision 23878)
5@@ -0,0 +1,1058 @@
6+%Test Name: IceOceanCoupling Dan Goldberg'd setup
7+%ISSM/MITgcm coupled set-up
8+%
9+%Script control parameters
10+steps=[1 2 3 4 5 6 7 8 9 10 11 12];
11+final_time=1/365;
12+
13+%To download and recompile MITgcm from scratch:
14+rm -rf $ISSM_DIR/test/MITgcm/install
15+rm -rf $ISSM_DIR/test/MITgcm/build
16+rm $ISSM_DIR/test/MITgcm/code/SIZE.h
17+
18+%Organizer
19+mkdir Models
20+org=organizer('repository','Models/','prefix','IceOcean.','steps',steps);
21+
22+presentdirectory=pwd;
23+
24+% {{{ Parameters:
25+if perform(org,'Parameters'),
26+ Nx=3; %number of longitude cells
27+ Ny=200; %number of latitude cells
28+ Nz=90; %number of MITgcm vertical cells
29+ nPx=1; %number of MITgcm processes to use in x direction
30+ nPy=1; %number of MITgcm processes to use in y direction
31+ xgOrigin=1; %origin of longitude
32+ ygOrigin=1; %origin of latitude
33+ dLong=1; %longitude grid spacing
34+ dLat=1; %latitude grid spacing
35+ delZ=10; %thickness of vertical levels
36+ icefront_position_ratio=.75;
37+ ice_thickness=1000;
38+ rho_ice=917;
39+ rho_water=1028.14;
40+ di=rho_ice/rho_water;
41+
42+ % MITgcm initial and lateral boundary conditions
43+ iniSalt = 34.4; % initial salinity (PSU)
44+ iniTheta = -1.9; % initial potential temperature (deg C)
45+ obcSalt = 34.4; % open boundary salinity (PSU)
46+ obcTheta = 1.0; % open boundary potential temperature (deg C)
47+ mlDepth = 120.; % mixed layer depth (m)
48+ mlSalt = 33.4; % open boundary salinity (PSU)
49+ mlTheta = -1.9; % open boundary potential temperature (deg C)
50+ obcUvel = -0.1; % open boundary velocity (m/s)
51+
52+ MITgcmDeltaT=600; % MITgcm time step in seconds
53+ y2s=31536000; % year to seconds conversion, i.e., seconds per year
54+
55+ % start_time and time_step
56+ start_time=0; % in decimal years
57+ time_step=1/(365*24); % coupling interval in decimal years
58+ async_step_MITgcm_multiplier=1; % used to reduce run time for MItgcm
59+
60+ % bedrock/bathymetry
61+ hmax=1000;
62+ trough_depth=200;
63+ deltah=300;
64+ sea_level=1095;
65+
66+ % issm settings:
67+ numlayers=10;
68+
69+ savedata(org, Nx, Ny, nPx, nPy, Nz, dLong, dLat, delZ, xgOrigin, ...
70+ ygOrigin, icefront_position_ratio, ice_thickness, rho_ice, ...
71+ rho_water, di, hmax, trough_depth, deltah, sea_level, ...
72+ iniSalt, iniTheta, obcSalt, obcTheta, mlDepth, mlSalt, ...
73+ mlTheta, obcUvel, start_time, time_step, MITgcmDeltaT, y2s,...
74+ numlayers,async_step_MITgcm_multiplier);
75+end
76+% }}}
77+% {{{ Bathymetry:
78+if perform(org,'Bathymetry'),
79+
80+ loaddata(org,'Parameters');
81+ %create lat,long
82+ lat=(ygOrigin+dLat/2):dLat:(ygOrigin+Ny*dLat);
83+ long=(xgOrigin+dLong/2):dLong:(xgOrigin+Nx*dLong);
84+ [lat long]=meshgrid(lat,long);
85+
86+ longmin=min(long(:));
87+ longmax=max(long(:));
88+ latmin=min(lat(:));
89+ latmax=max(lat(:));
90+
91+ %create bedrock/bathymetry:
92+ bedrock=zeros(Nx,Ny);
93+ bedrock=hmax-deltah*tanh(pi*(2*(lat-latmin)./(latmax-latmin)-1))+ ...
94+ trough_depth*cos(2*pi*long./(longmax-longmin));
95+
96+ %save bathymetry file for MITgcm
97+ bathymetry=bedrock-sea_level;
98+ savedata(org,lat,long,bathymetry);
99+
100+end
101+% }}}
102+% {{{ IceSheetGeometry:
103+if perform(org,'IceSheetGeometry'),
104+
105+ loaddata(org,'Parameters');
106+ loaddata(org,'Bathymetry');
107+ latmin=min(lat(:));
108+ latmax=max(lat(:));
109+
110+ %put ice_thickness constant layer of ice over the bathymetry, unless it floats:
111+ s=size(bathymetry);
112+ thickness=ice_thickness*ones(s);
113+
114+ %figure out ice shelf:
115+ pos=find(-di*thickness>bathymetry);
116+ iceshelf_mask=zeros(s);
117+ iceshelf_mask(pos)=1;
118+
119+ ice_mask=ones(s);
120+ pos=find((lat-latmin)/(latmax-latmin)>(icefront_position_ratio));
121+ ice_mask(pos)=0;
122+ iceshelf_mask(pos)=0;
123+
124+ %compute draft of ice shelf:
125+ draft=bathymetry;
126+ pos=find(iceshelf_mask);
127+ draft(pos)=-di*thickness(pos);
128+ pos=find(~ice_mask);
129+ draft(pos)=0;
130+
131+ savedata(org,ice_mask,iceshelf_mask,draft,thickness);
132+end
133+% }}}
134+
135+%Configure MITgcm
136+% {{{ GetMITgcm:
137+if perform(org,'GetMITgcm'),
138+ system([pwd '/../MITgcm/get_mitgcm.sh']);
139+end
140+% }}}
141+% {{{ BuildMITgcm:
142+if perform(org,'BuildMITgcm'),
143+
144+ %load data:
145+ loaddata(org,'Parameters');
146+
147+ %specify computational grid in SIZE.h
148+ if ~exist('../MITgcm/code/SIZE.h')
149+ fidi=fopen('../MITgcm/code/SIZE.h.bak','r');
150+ fido=fopen('../MITgcm/code/SIZE.h','w');
151+ tline = fgetl(fidi);
152+ fprintf(fido,'%s\n',tline);
153+ while 1
154+ tline = fgetl(fidi);
155+ if ~ischar(tline), break, end
156+ %do the change here:
157+ if strcmpi(tline,' & sNx = 20,'),
158+ fprintf(fido,'%s%i%s\n',' & sNx = ',round(Nx/nPx),',');
159+ continue;
160+ end
161+ if strcmpi(tline,' & sNy = 20,'),
162+ fprintf(fido,'%s%i%s\n',' & sNy = ',round(Ny/nPy),',');
163+ continue;
164+ end
165+ if strcmpi(tline,' & nPx = 1,'),
166+ fprintf(fido,'%s%i%s\n',' & nPx = ',nPx,',');
167+ continue;
168+ end
169+ if strcmpi(tline,' & nPy = 2,'),
170+ fprintf(fido,'%s%i%s\n',' & nPy = ',nPy,',');
171+ continue;
172+ end
173+ fprintf(fido,'%s\n',tline);
174+ end
175+ %close files
176+ fclose(fidi);
177+ fclose(fido);
178+ end
179+
180+ system(['../MITgcm/build.sh generic ' pwd '/../MITgcm']);
181+end
182+% }}}
183+addpath(recursivepath([pwd '/../MITgcm']));
184+% {{{ RunUncoupledMITgcm:
185+if perform(org,'RunUncoupledMITgcm'),
186+
187+ %load data:
188+ loaddata(org,'Parameters');
189+ loaddata(org,'Bathymetry');
190+ loaddata(org,'IceSheetGeometry');
191+ endtime = round(MITgcmDeltaT * ...
192+ floor(time_step*y2s*async_step_MITgcm_multiplier/MITgcmDeltaT));
193+
194+ % {{{ prepare MITgcm
195+ % rename previous run directory and create new one
196+ if exist ('run.old')
197+ !\rm -rf run.old
198+ end
199+ if exist ('run')
200+ !\mv run run.old
201+ end
202+ !\mkdir run
203+ !\cp ../MITgcm/build/mitgcmuv run
204+ !\cp ../MITgcm/input/* run
205+ !\cp ../MITgcm/input/eedata_uncoupled run/eedata
206+
207+ %load data:
208+ loaddata(org,'Parameters');
209+
210+ % initial salinity
211+ S=iniSalt*ones(Nx,Ny,Nz);
212+ writebin('run/Salt.bin',S);
213+
214+ % initial temperature
215+ T=iniTheta*ones(Nx,Ny,Nz);
216+ writebin('run/Theta.bin',T);
217+
218+ % initial velocity
219+ Z=zeros(Nx,Ny,Nz);
220+ writebin('run/Uvel.bin',Z);
221+ writebin('run/Vvel.bin',Z);
222+
223+ % initial sea surface height
224+ Z=zeros(Nx,Ny);
225+ writebin('run/Etan.bin',Z);
226+
227+ % salinity boundary conditions
228+ S=obcSalt*ones(Ny,Nz);
229+ thk=delZ*ones(Nz,1);
230+ bot=cumsum(thk);
231+ ik=find(bot<=mlDepth);
232+ S(:,ik)=mlSalt;
233+ writebin('run/OBs.bin',S);
234+
235+ % temperature boundary conditions
236+ T=obcTheta*ones(Ny,Nz);
237+ T(:,ik)=mlTheta;
238+ writebin('run/OBt.bin',T);
239+
240+ % zonal velocity boundary conditions
241+ U=obcUvel*ones(Ny,Nz);
242+ writebin('run/OBu.bin',U);
243+
244+ % zero boundary conditions
245+ Z=zeros(Ny,Nz);
246+ writebin('run/zeros.bin',Z);
247+
248+ % build parameter file data.obcs
249+ fidi=fopen('../MITgcm/input/data.obcs','r');
250+ fido=fopen('run/data.obcs','w');
251+ tline = fgetl(fidi);
252+ fprintf(fido,'%s\n',tline);
253+ while 1
254+ tline = fgetl(fidi);
255+ if ~ischar(tline), break, end
256+ %do the change here:
257+ if strcmpi(tline,' OB_Iwest = 40*1,'),
258+ fprintf(fido,'%s%i%s\n',' OB_Iwest = ',Ny,'*1,');
259+ continue;
260+ end
261+ if strcmpi(tline,' OB_Ieast = 40*-1,'),
262+ fprintf(fido,'%s%i%s\n',' OB_Ieast = ',Ny,'*-1,');
263+ continue;
264+ end
265+ fprintf(fido,'%s\n',tline);
266+ end
267+ %close files
268+ fclose(fidi);
269+ fclose(fido);
270+
271+ %save bathymetry and bedrock in run directory
272+ writebin('run/bathymetry.bin',bathymetry);
273+ writebin('run/icetopo.bin',draft);
274+ % }}}
275+
276+ %start looping:
277+ for t=start_time:time_step:final_time,
278+ disp(['Year: ' num2str(t)])
279+ % {{{ generate MITgcm parameter file data
280+ fidi=fopen('../MITgcm/input/data','r');
281+ fido=fopen('run/data','w');
282+ tline = fgetl(fidi);
283+ fprintf(fido,'%s\n',tline);
284+ while 1
285+ tline = fgetl(fidi);
286+ if ~ischar(tline), break, end
287+ %do the change here:
288+ if strcmpi(tline,' xgOrigin = 0.0,'),
289+ fprintf(fido,'%s%i%s\n',' xgOrigin = ',xgOrigin,',');
290+ continue;
291+ end
292+ if strcmpi(tline,' ygOrigin = -80.0,'),
293+ fprintf(fido,'%s%i%s\n',' ygOrigin = ',ygOrigin,',');
294+ continue;
295+ end
296+ if strcmpi(tline,' delX = 20*0.25,'),
297+ fprintf(fido,'%s%i*%g%s\n',' delX = ',Nx,dLong,',');
298+ continue;
299+ end
300+ if strcmpi(tline,' delY = 20*0.25,'),
301+ fprintf(fido,'%s%i*%g%s\n',' delY = ',Ny,dLat,',');
302+ continue;
303+ end
304+ if strcmpi(tline,' delZ = 30*30.0,'),
305+ fprintf(fido,'%s%i*%g%s\n',' delZ = ',Nz,delZ,',');
306+ continue;
307+ end
308+ if strcmpi(tline,' endTime=2592000.,'),
309+ fprintf(fido,'%s%i%s\n',' endTime= ',endtime,',');
310+ continue;
311+ end
312+ if strcmpi(tline,' deltaT=1200.0,'),
313+ fprintf(fido,'%s%i%s\n',' deltaT= ',MITgcmDeltaT,',');
314+ continue;
315+ end
316+ if strcmpi(tline,' pChkptFreq=2592000.,'),
317+ fprintf(fido,'%s%i%s\n',' pChkptFreq= ',endtime,',');
318+ continue;
319+ end
320+ if strcmpi(tline,' taveFreq=2592000.,'),
321+ fprintf(fido,'%s%i%s\n',' taveFreq= ',endtime,',');
322+ continue;
323+ end
324+ if strcmpi(tline,' rhoConst=1030.,'),
325+ fprintf(fido,'%s%i%s\n',' rhoConst= ',rho_water,',');
326+ continue;
327+ end
328+ if strcmpi(tline,' rhoNil=1030.,'),
329+ fprintf(fido,'%s%i%s\n',' rhoNil= ',rho_water,',');
330+ continue;
331+ end
332+ fprintf(fido,'%s\n',tline);
333+ end
334+ %close files
335+ fclose(fidi);
336+ fclose(fido);
337+ % }}}
338+ % {{{ generate initial MITgcm conditions
339+
340+ ds=round(endtime/MITgcmDeltaT);
341+ if t>start_time
342+ % Read pickup file
343+ fnm=['run/pickup.' myint2str(ds,10) '.data'];
344+ U=readbin(fnm,[Nx Ny Nz],1,'real*8',0);
345+ V=readbin(fnm,[Nx Ny Nz],1,'real*8',1);
346+ T=readbin(fnm,[Nx Ny Nz],1,'real*8',2);
347+ S=readbin(fnm,[Nx Ny Nz],1,'real*8',3);
348+ E=readbin(fnm,[Nx Ny],1,'real*8',8*Nz);
349+ writebin('run/Salt.bin' ,S);
350+ writebin('run/Theta.bin',T);
351+ writebin('run/Uvel.bin' ,U);
352+ writebin('run/Vvel.bin' ,V);
353+ writebin('run/Etan.bin' ,E);
354+ end
355+
356+ % }}}
357+ % {{{ system call to run MITgcm
358+ cd run
359+ eval(['!mpirun -np ' int2str(nPx*nPy) ' ./mitgcmuv']);
360+ ts=round((t+time_step)*y2s/MITgcmDeltaT);
361+ eval(['!\mv STDERR.0000 STDERR_' myint2str(ts,10) '.data'])
362+ eval(['!\mv STDOUT.0000 STDOUT_' myint2str(ts,10) '.data'])
363+ eval(['!\cp hFacC.data hFacC_' myint2str(ts,10) '.data'])
364+ eval(['!\cp icetopo.bin icetopo_' myint2str(ts,10) '.data'])
365+ for fld={'S','T','U','V','Eta', ...
366+ 'SHICE_heatFluxtave','SHICE_fwFluxtave'}
367+ eval(['!\mv ' fld{1} '.' myint2str(ds,10) '.data ' ...
368+ fld{1} '_' myint2str(ts,10) '.data'])
369+ end
370+ cd ..
371+ % }}}
372+ end
373+end
374+% }}}
375+
376+%Configure ISSM
377+% {{{ CreateMesh:
378+if perform(org,'CreateMesh'),
379+
380+ loaddata(org,'Parameters');
381+ loaddata(org,'Bathymetry');
382+ loaddata(org,'IceSheetGeometry');
383+
384+ %create model:
385+ md=model();
386+
387+ %Grab lat,long from MITgcm:
388+ lat=lat(:);
389+ long=long(:);
390+
391+ %project lat,long:
392+ [x,y]=ll2xy(lat,long,-1);
393+
394+ index=[];
395+ % C D
396+ % A B
397+ for j=1:Ny-1,
398+ for i=1:Nx-1,
399+ A=(j-1)*Nx+i;
400+ B=(j-1)*Nx+i+1;
401+ C=j*Nx+i;
402+ D=j*Nx+i+1;
403+ index(end+1,:)=[A B C];
404+ index(end+1,:)=[C B D];
405+ end
406+ end
407+
408+ %fill mesh and model:
409+ md=meshconvert(md,index,x,y);
410+ md.mesh.lat=lat;
411+ md.mesh.long=long;
412+
413+ savemodel(org,md);
414+
415+end
416+% }}}
417+% {{{ MeshGeometry:
418+if perform(org,'MeshGeometry'),
419+
420+ loaddata(org,'Parameters');
421+ loaddata(org,'CreateMesh');
422+ loaddata(org,'Bathymetry');
423+ loaddata(org,'IceSheetGeometry');
424+
425+ %transfer to vertices:
426+ bathymetry=bathymetry(:);
427+ iceshelf_mask=iceshelf_mask(:);
428+ ice_mask=ice_mask(:);
429+ thickness=thickness(:);
430+ draft=draft(:);
431+
432+ %start filling some of the fields
433+ md.geometry.bed=bathymetry;
434+ md.geometry.thickness=thickness;
435+ md.geometry.base=md.geometry.bed;
436+ pos=find(iceshelf_mask); md.geometry.base(pos)=draft(pos);
437+ md.geometry.surface=md.geometry.base+md.geometry.thickness;
438+
439+ %nothing passes icefront:
440+ pos=find(~ice_mask);
441+ md.geometry.thickness(pos)=1;
442+ md.geometry.surface(pos)=(1-di)*md.geometry.thickness(pos);
443+ md.geometry.base(pos)=-di*md.geometry.thickness(pos);
444+
445+ %level sets:
446+ md.mask.groundedice_levelset=-ones(md.mesh.numberofvertices,1);
447+ md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
448+
449+ pos=find(ice_mask); md.mask.ice_levelset(pos)=-1;
450+ pos=find(~iceshelf_mask & ice_mask); md.mask.groundedice_levelset(pos)=1;
451+
452+ %identify edges of grounded ice:
453+ groundedice_levelset=md.mask.groundedice_levelset;
454+ for i=1:md.mesh.numberofelements,
455+ m=groundedice_levelset(md.mesh.elements(i,:));
456+ if abs(sum(m))~=3,
457+ pos=find(m==1); md.mask.groundedice_levelset(md.mesh.elements(i,pos))=0;
458+ end
459+ end
460+
461+ %identify edges of ice:
462+ ice_levelset=md.mask.ice_levelset;
463+ for i=1:md.mesh.numberofelements,
464+ m=ice_levelset(md.mesh.elements(i,:));
465+ if abs(sum(m))~=3,
466+ pos=find(m==-1); md.mask.ice_levelset(md.mesh.elements(i,pos))=0;
467+ end
468+ end
469+
470+ savemodel(org,md);
471+end
472+% }}}
473+% {{{ ParameterizeIce:
474+if perform(org,'ParameterizeIce'),
475+
476+ loaddata(org,'Parameters');
477+ loaddata(org,'CreateMesh');
478+ loaddata(org,'MeshGeometry');
479+
480+ %miscellaneous
481+ md.miscellaneous.name='test4002';
482+
483+ %initial velocity:
484+ md.initialization.vx=zeros(md.mesh.numberofvertices,1);
485+ md.initialization.vy=zeros(md.mesh.numberofvertices,1);
486+ md.initialization.vz=zeros(md.mesh.numberofvertices,1);
487+
488+ %friction:
489+ md.friction.coefficient=30*ones(md.mesh.numberofvertices,1);
490+ pos=find(md.mask.groundedice_levelset<=0);
491+ md.friction.coefficient(pos)=0;
492+ md.friction.p=ones(md.mesh.numberofelements,1);
493+ md.friction.q=ones(md.mesh.numberofelements,1);
494+
495+ %temperatures and surface mass balance:
496+ md.initialization.temperature=(273.15-20)*ones(md.mesh.numberofvertices,1);
497+ md.initialization.pressure=md.materials.rho_ice*md.constants.g*(md.geometry.surface-md.geometry.base);
498+ md.smb.mass_balance = [1*ones(md.mesh.numberofvertices,1); 1];
499+
500+ %Flow law
501+ md.materials.rheology_B=paterson(md.initialization.temperature);
502+ md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
503+ md.damage.D=zeros(md.mesh.numberofvertices,1);
504+ md.damage.spcdamage=NaN*ones(md.mesh.numberofvertices,1);
505+
506+ %the spcs going
507+ md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);
508+ md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
509+ md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);
510+ md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
511+ md.stressbalance.loadingforce=0*ones(md.mesh.numberofvertices,3);
512+ md.masstransport.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
513+
514+ %deal with water:
515+ pos=find(md.mask.ice_levelset>0);
516+ md.stressbalance.spcvx(pos)=0;
517+ md.stressbalance.spcvy(pos)=0;
518+ md.stressbalance.spcvz(pos)=0;
519+ md.masstransport.spcthickness(pos)=0;
520+
521+ %get some flux at the ice divide:
522+ pos=find(md.mesh.lat==min(md.mesh.lat));
523+ md.stressbalance.spcvy(pos)=200;
524+
525+ %deal with boundaries, excluding icefront:
526+ vertex_on_boundary=zeros(md.mesh.numberofvertices,1);
527+ vertex_on_boundary(md.mesh.segments(:,1:2))=1;
528+ pos=find(vertex_on_boundary & md.mask.groundedice_levelset<=0);
529+ md.stressbalance.spcvx(pos)=md.initialization.vx(pos);
530+ md.stressbalance.spcvy(pos)=md.initialization.vy(pos);
531+ md.stressbalance.spcvz(pos)=md.initialization.vz(pos);
532+ md.masstransport.spcthickness(pos)=md.geometry.thickness(pos);
533+
534+ md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
535+ md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
536+ md.thermal.spctemperature=[md.initialization.temperature; 1]; %impose observed temperature on surface
537+ md.basalforcings.geothermalflux=.064*ones(md.mesh.numberofvertices,1);
538+
539+ %flow equations:
540+ md=setflowequation(md,'SSA','all');
541+
542+ savemodel(org,md);
543+end
544+% }}}
545+% {{{ RunUncoupledISSM:
546+if perform(org,'RunUncoupledISSM'),
547+
548+ loaddata(org,'Parameters');
549+ loaddata(org,'ParameterizeIce');
550+
551+ %timestepping:
552+ md.timestepping.final_time=final_time;
553+ md.timestepping.time_step=time_step;
554+ md.transient.isgroundingline=1;
555+ md.transient.isthermal=0;
556+ md.groundingline.migration='SubelementMigration';
557+ md.groundingline.melt_interpolation='SubelementMelt2';
558+ md.groundingline.friction_interpolation='SubelementFriction2';
559+
560+ md.cluster=generic('name',oshostname(),'np',2);
561+ md=solve(md,'Transient');
562+
563+ savemodel(org,md);
564+end
565+% }}}
566+
567+%Run MITgcm/ISSM
568+% {{{ RunCoupledMITgcmISSM:
569+if perform(org,'RunCoupledMITgcmISSM'),
570+
571+ %load data:
572+ loaddata(org,'Parameters');
573+ loaddata(org,'ParameterizeIce');
574+ loaddata(org,'Bathymetry');
575+ loaddata(org,'IceSheetGeometry');
576+ endtime = round(MITgcmDeltaT * ...
577+ floor(time_step*y2s*async_step_MITgcm_multiplier/MITgcmDeltaT));
578+
579+ % {{{ prepare MITgcm
580+ % rename previous run directory and create new one
581+ if exist ('run.old')
582+ !\rm -rf run.old
583+ end
584+ if exist ('run')
585+ !\mv run run.old
586+ end
587+ !\mkdir run
588+ !\cp ../MITgcm/build/mitgcmuv run
589+ !\cp ../MITgcm/input/* run
590+ !\cp ../MITgcm/input/eedata_uncoupled run/eedata
591+
592+ % initial salinity
593+ S=iniSalt*ones(Nx,Ny,Nz);
594+ writebin('run/Salt.bin',S);
595+
596+ % initial temperature
597+ T=iniTheta*ones(Nx,Ny,Nz);
598+ writebin('run/Theta.bin',T);
599+
600+ % initial velocity
601+ Z=zeros(Nx,Ny,Nz);
602+ writebin('run/Uvel.bin',Z);
603+ writebin('run/Vvel.bin',Z);
604+
605+ % initial sea surface height
606+ Z=zeros(Nx,Ny);
607+ writebin('run/Etan.bin',Z);
608+
609+ % salinity boundary conditions
610+ S=obcSalt*ones(Ny,Nz);
611+ thk=delZ*ones(Nz,1);
612+ bot=cumsum(thk);
613+ ik=find(bot<=mlDepth);
614+ S(:,ik)=mlSalt;
615+ writebin('run/OBs.bin',S);
616+
617+ % temperature boundary conditions
618+ T=obcTheta*ones(Ny,Nz);
619+ T(:,ik)=mlTheta;
620+ writebin('run/OBt.bin',T);
621+
622+ % zonal velocity boundary conditions
623+ U=obcUvel*ones(Ny,Nz);
624+ writebin('run/OBu.bin',U);
625+
626+ % zero boundary conditions
627+ Z=zeros(Ny,Nz);
628+ writebin('run/zeros.bin',Z);
629+
630+ % build parameter file data.obcs
631+ fidi=fopen('../MITgcm/input/data.obcs','r');
632+ fido=fopen('run/data.obcs','w');
633+ tline = fgetl(fidi);
634+ fprintf(fido,'%s\n',tline);
635+ while 1
636+ tline = fgetl(fidi);
637+ if ~ischar(tline), break, end
638+ %do the change here:
639+ if strcmpi(tline,' OB_Iwest = 40*1,'),
640+ fprintf(fido,'%s%i%s\n',' OB_Iwest = ',Ny,'*1,');
641+ continue;
642+ end
643+ if strcmpi(tline,' OB_Ieast = 40*-1,'),
644+ fprintf(fido,'%s%i%s\n',' OB_Ieast = ',Ny,'*-1,');
645+ continue;
646+ end
647+ fprintf(fido,'%s\n',tline);
648+ end
649+ %close files
650+ fclose(fidi);
651+ fclose(fido);
652+
653+ %save bathymetry in MITgcm run directory
654+ writebin('run/bathymetry.bin',bathymetry);
655+ % }}}
656+
657+ % {{{ ISSM settings:
658+
659+ setenv('DYLD_LIBRARY_PATH', '/usr/local/gfortran/lib')
660+ %timestepping:
661+ md.timestepping.start_time=start_time;
662+ md.timestepping.final_time=final_time;
663+ md.timestepping.time_step=time_step;
664+ md.cluster=generic('name',oshostname(),'np',2);
665+ md.results.TransientSolution.Base=md.geometry.base;
666+ md.transient.isgroundingline=1;
667+ md.transient.isthermal=0;
668+ md.groundingline.migration='SubelementMigration';
669+ md.groundingline.melt_interpolation='SubelementMelt2';
670+ md.groundingline.friction_interpolation='SubelementFriction2';
671+
672+ % }}}
673+
674+ %start looping:
675+ results=md.results;
676+
677+ for t=start_time:time_step:final_time
678+ disp(['Year: ' num2str(t)])
679+
680+ %send draft from ISSM to MITgcm:
681+ draft=md.results.TransientSolution(end).Base;
682+ pos=find(md.mask.ice_levelset>0); draft(pos)=0;
683+ if t>start_time
684+ old_draft=readbin('run/icetopo.bin',[Nx,Ny]);
685+ end
686+ writebin('run/icetopo.bin',draft);
687+
688+ % {{{ generate MITgcm parameter file data
689+ fidi=fopen('../MITgcm/input/data','r');
690+ fido=fopen('run/data','w');
691+ tline = fgetl(fidi);
692+ fprintf(fido,'%s\n',tline);
693+ while 1
694+ tline = fgetl(fidi);
695+ if ~ischar(tline), break, end
696+ %do the change here:
697+ if strcmpi(tline,' xgOrigin = 0.0,'),
698+ fprintf(fido,'%s%i%s\n',' xgOrigin = ',xgOrigin,',');
699+ continue;
700+ end
701+ if strcmpi(tline,' ygOrigin = -80.0,'),
702+ fprintf(fido,'%s%i%s\n',' ygOrigin = ',ygOrigin,',');
703+ continue;
704+ end
705+ if strcmpi(tline,' delX = 20*0.25,'),
706+ fprintf(fido,'%s%i*%g%s\n',' delX = ',Nx,dLong,',');
707+ continue;
708+ end
709+ if strcmpi(tline,' delY = 20*0.25,'),
710+ fprintf(fido,'%s%i*%g%s\n',' delY = ',Ny,dLat,',');
711+ continue;
712+ end
713+ if strcmpi(tline,' delZ = 30*30.0,'),
714+ fprintf(fido,'%s%i*%g%s\n',' delZ = ',Nz,delZ,',');
715+ continue;
716+ end
717+ if strcmpi(tline,' endTime=2592000.,'),
718+ fprintf(fido,'%s%i%s\n',' endTime= ',endtime,',');
719+ continue;
720+ end
721+ if strcmpi(tline,' deltaT=1200.0,'),
722+ fprintf(fido,'%s%i%s\n',' deltaT= ',MITgcmDeltaT,',');
723+ continue;
724+ end
725+ if strcmpi(tline,' pChkptFreq=2592000.,'),
726+ fprintf(fido,'%s%i%s\n',' pChkptFreq= ',endtime,',');
727+ continue;
728+ end
729+ if strcmpi(tline,' taveFreq=2592000.,'),
730+ fprintf(fido,'%s%i%s\n',' taveFreq= ',endtime,',');
731+ continue;
732+ end
733+ if strcmpi(tline,' rhoConst=1030.,'),
734+ fprintf(fido,'%s%i%s\n',' rhoConst= ',rho_water,',');
735+ continue;
736+ end
737+ if strcmpi(tline,' rhoNil=1030.,'),
738+ fprintf(fido,'%s%i%s\n',' rhoNil= ',rho_water,',');
739+ continue;
740+ end
741+ fprintf(fido,'%s\n',tline);
742+ end
743+ %close files
744+ fclose(fidi);
745+ fclose(fido);
746+ % }}}
747+
748+ % {{{ generate initial MITgcm conditions
749+ ds=round(endtime/MITgcmDeltaT);
750+ if t>start_time
751+ % Read pickup file
752+ fnm=['run/pickup.' myint2str(ds,10) '.data'];
753+ U=readbin(fnm,[Nx Ny Nz],1,'real*8',0);
754+ V=readbin(fnm,[Nx Ny Nz],1,'real*8',1);
755+ T=readbin(fnm,[Nx Ny Nz],1,'real*8',2);
756+ S=readbin(fnm,[Nx Ny Nz],1,'real*8',3);
757+ E=readbin(fnm,[Nx Ny],1,'real*8',8*Nz);
758+
759+ % find indices of locations where ice shelf retreated
760+ h=readbin('run/hFacC.data',[Nx Ny Nz]);
761+ msk=sum(h,3);
762+ msk(find(msk))=1;
763+ [iw jw]=find(msk); % horizontal indices where there is water
764+ tmp=reshape(draft,[Nx,Ny])-old_draft;
765+ tmp(find(tmp<0))=0;
766+ [im jm]=find(tmp); % horizontal indices where there is melt
767+
768+ % Extrapolate T/S to locations where ice shelf retreated
769+ for i=1:length(im)
770+
771+ % first try vertical extrapolation
772+ in=find(h(im(i),jm(i),:));
773+ if length(in)>0;
774+ S(im(i),jm(i),1:min(in) ) = S(im(i),jm(i),min(in));
775+ T(im(i),jm(i),1:min(in) ) = T(im(i),jm(i),min(in));
776+ continue
777+ end
778+
779+ % if not succesful, use closest neighbor horizontal extrapolation
780+ [y c]=min((iw-im(i)).^2+(jw-jm(i)).^2);
781+ salt=squeeze(S(iw(c),jw(c),:)); % salinity profile of closest neighbor
782+ temp=squeeze(T(iw(c),jw(c),:)); % salinity profile of closest neighbor
783+ in=find(h(iw(c),jw(c),:));
784+ salt(1:min(in))=salt(min(in));
785+ temp(1:min(in))=temp(min(in));
786+ salt(max(in):end)=salt(max(in));
787+ temp(max(in):end)=temp(max(in));
788+ S(im(i),jm(i),:)=salt;
789+ T(im(i),jm(i),:)=temp;
790+ end
791+
792+ % Write initial conditions
793+ writebin('run/Salt.bin' ,S);
794+ writebin('run/Theta.bin',T);
795+ writebin('run/Uvel.bin' ,U);
796+ writebin('run/Vvel.bin' ,V);
797+ writebin('run/Etan.bin' ,E);
798+ end
799+ % }}}
800+
801+ % {{{ system call to run MITgcm
802+ cd run
803+ eval(['!mpiexec -np ' int2str(nPx*nPy) ' ./mitgcmuv']);
804+ ts=round((t+time_step)*y2s/MITgcmDeltaT);
805+ eval(['!\mv STDERR.0000 STDERR_' myint2str(ts,10) '.data'])
806+ eval(['!\mv STDOUT.0000 STDOUT_' myint2str(ts,10) '.data'])
807+ eval(['!\cp hFacC.data hFacC_' myint2str(ts,10) '.data'])
808+ eval(['!\cp icetopo.bin icetopo_' myint2str(ts,10) '.data'])
809+ for fld={'S','T','U','V','Eta', ...
810+ 'SHICE_heatFluxtave','SHICE_fwFluxtave'}
811+ eval(['!\mv ' fld{1} '.' myint2str(ds,10) '.data ' ...
812+ fld{1} '_' myint2str(ts,10) '.data'])
813+ end
814+ cd ..
815+ % }}}
816+
817+ %get melting rates from MITgcm
818+ %upward fresh water flux (kg/m^2/s):
819+ fnm=['run/SHICE_fwFluxtave_' myint2str(ts,10) '.data'];
820+ melting_rate=readbin(fnm,[Nx Ny]);
821+
822+ %send averaged melting rate to ISSM
823+ %downward fresh water flux (m/y):
824+ melting_rate=-melting_rate(:)*y2s/rho_ice;
825+ md.basalforcings.floatingice_melting_rate=melting_rate;
826+
827+ % {{{ run ISSM and recover results
828+
829+ md.timestepping.start_time=t;
830+ md.timestepping.final_time=t+time_step;;
831+ md=solve(md,'Transient');
832+
833+ base=md.results.TransientSolution(end).Base;
834+ thickness=md.results.TransientSolution(end).Thickness;
835+ md.geometry.base=base;
836+ md.geometry.thickness=thickness;
837+ md.geometry.surface=md.geometry.base+md.geometry.thickness;
838+ md.initialization.vx=md.results.TransientSolution(end).Vx;
839+ md.initialization.vy=md.results.TransientSolution(end).Vy;
840+ md.initialization.vel=md.results.TransientSolution(end).Vel;
841+ md.initialization.pressure=md.results.TransientSolution(end).Pressure;
842+ md.mask.groundedice_levelset=md.results.TransientSolution(end).MaskGroundediceLevelset;
843+ md.results.TransientSolution(end).FloatingiceMeltingRate=md.basalforcings.floatingice_melting_rate;
844+
845+ %save these results in the model, otherwise, they'll be wiped out
846+ results(end+1)=md.results;
847+
848+ % }}}
849+
850+
851+ end
852+
853+ md.results=results;
854+ savemodel(org,md);
855+end
856+% }}}
857+% {{{ RunCoupledMITgcmISSM2:
858+if perform(org,'RunCoupledMITgcmISSM2'),
859+
860+ loaddata(org,'Parameters');
861+ loaddata(org,'ParameterizeIce');
862+ loaddata(org,'Bathymetry');
863+ loaddata(org,'IceSheetGeometry');
864+ endtime = round(MITgcmDeltaT * floor(final_time*y2s/MITgcmDeltaT));
865+ outputtime = round(MITgcmDeltaT * floor(time_step*y2s/MITgcmDeltaT));
866+
867+ % {{{ prepare MITgcm
868+ % rename previous run directory and create new one
869+ if exist ('run.old')
870+ !\rm -rf run.old
871+ end
872+ if exist ('run')
873+ !\mv run run.old
874+ end
875+ !\mkdir run
876+ !\cp ../MITgcm/build/mitgcmuv run
877+ !\cp ../MITgcm/input/* run
878+
879+ % initial salinity
880+ S=iniSalt*ones(Nx,Ny,Nz);
881+ writebin('run/Salt.bin',S);
882+
883+ % initial temperature
884+ T=iniTheta*ones(Nx,Ny,Nz);
885+ writebin('run/Theta.bin',T);
886+
887+ % initial velocity
888+ Z=zeros(Nx,Ny,Nz);
889+ writebin('run/Uvel.bin',Z);
890+ writebin('run/Vvel.bin',Z);
891+
892+ % initial sea surface height
893+ Z=zeros(Nx,Ny);
894+ writebin('run/Etan.bin',Z);
895+
896+ % salinity boundary conditions
897+ S=obcSalt*ones(Ny,Nz);
898+ thk=delZ*ones(Nz,1);
899+ bot=cumsum(thk);
900+ ik=find(bot<=mlDepth);
901+ S(:,ik)=mlSalt;
902+ writebin('run/OBs.bin',S);
903+
904+ % temperature boundary conditions
905+ T=obcTheta*ones(Ny,Nz);
906+ T(:,ik)=mlTheta;
907+ writebin('run/OBt.bin',T);
908+
909+ % zonal velocity boundary conditions
910+ U=obcUvel*ones(Ny,Nz);
911+ writebin('run/OBu.bin',U);
912+
913+ % zero boundary conditions
914+ Z=zeros(Ny,Nz);
915+ writebin('run/zeros.bin',Z);
916+
917+ % build parameter file data.obcs
918+ fidi=fopen('../MITgcm/input/data.obcs','r');
919+ fido=fopen('run/data.obcs','w');
920+ tline = fgetl(fidi);
921+ fprintf(fido,'%s\n',tline);
922+ while 1
923+ tline = fgetl(fidi);
924+ if ~ischar(tline), break, end
925+ %do the change here:
926+ if strcmpi(tline,' OB_Iwest = 40*1,'),
927+ fprintf(fido,'%s%i%s\n',' OB_Iwest = ',Ny,'*1,');
928+ continue;
929+ end
930+ if strcmpi(tline,' OB_Ieast = 40*-1,'),
931+ fprintf(fido,'%s%i%s\n',' OB_Ieast = ',Ny,'*-1,');
932+ continue;
933+ end
934+ fprintf(fido,'%s\n',tline);
935+ end
936+ %close files
937+ fclose(fidi);
938+ fclose(fido);
939+
940+ %save bathymetry and bedrock in run directory
941+ writebin('run/bathymetry.bin',bathymetry);
942+ writebin('run/icetopo.bin',draft);
943+ % }}}
944+ % {{{ generate MITgcm parameter file data
945+ fidi=fopen('../MITgcm/input/data','r');
946+ fido=fopen('run/data','w');
947+ tline = fgetl(fidi);
948+ fprintf(fido,'%s\n',tline);
949+ while 1
950+ tline = fgetl(fidi);
951+ if ~ischar(tline), break, end
952+ %do the change here:
953+ if strcmpi(tline,' xgOrigin = 0.0,'),
954+ fprintf(fido,'%s%i%s\n',' xgOrigin = ',xgOrigin,',');
955+ continue;
956+ end
957+ if strcmpi(tline,' ygOrigin = -80.0,'),
958+ fprintf(fido,'%s%i%s\n',' ygOrigin = ',ygOrigin,',');
959+ continue;
960+ end
961+ if strcmpi(tline,' delX = 20*0.25,'),
962+ fprintf(fido,'%s%i*%g%s\n',' delX = ',Nx,dLong,',');
963+ continue;
964+ end
965+ if strcmpi(tline,' delY = 20*0.25,'),
966+ fprintf(fido,'%s%i*%g%s\n',' delY = ',Ny,dLat,',');
967+ continue;
968+ end
969+ if strcmpi(tline,' delZ = 30*30.0,'),
970+ fprintf(fido,'%s%i*%g%s\n',' delZ = ',Nz,delZ,',');
971+ continue;
972+ end
973+ if strcmpi(tline,' endTime=2592000.,'),
974+ fprintf(fido,'%s%i%s\n',' endTime= ',endtime,',');
975+ continue;
976+ end
977+ if strcmpi(tline,' deltaT=1200.0,'),
978+ fprintf(fido,'%s%i%s\n',' deltaT= ',MITgcmDeltaT,',');
979+ continue;
980+ end
981+ if strcmpi(tline,' pChkptFreq=2592000.,'),
982+ fprintf(fido,'%s%i%s\n',' pChkptFreq= ',endtime,',');
983+ continue;
984+ end
985+ if strcmpi(tline,' taveFreq=2592000.,'),
986+ fprintf(fido,'%s%i%s\n',' taveFreq= ',outputtime,',');
987+ continue;
988+ end
989+ if strcmpi(tline,' rhoConst=1030.,'),
990+ fprintf(fido,'%s%i%s\n',' rhoConst= ',rho_water,',');
991+ continue;
992+ end
993+ if strcmpi(tline,' rhoNil=1030.,'),
994+ fprintf(fido,'%s%i%s\n',' rhoNil= ',rho_water,',');
995+ continue;
996+ end
997+ fprintf(fido,'%s\n',tline);
998+ end
999+ %close files
1000+ fclose(fidi);
1001+ fclose(fido);
1002+ % }}}
1003+
1004+ md.transient.isoceancoupling=1;
1005+ md.transient.isgroundingline=1;
1006+ md.groundingline.migration='None';
1007+ md.groundingline.melt_interpolation='SubelementMelt2';
1008+ md.groundingline.friction_interpolation='SubelementFriction2';
1009+ md.timestepping.coupling_time=time_step;
1010+ md.timestepping.time_step=time_step;
1011+ md.timestepping.final_time=final_time-time_step;
1012+ md.cluster.npocean=nPx*nPy;
1013+ md.cluster.np=2;
1014+ md.cluster.executionpath=[pwd '/run'];
1015+ md.masstransport.requested_outputs={'default','BasalforcingsFloatingiceMeltingRate'};
1016+
1017+ md=solveiceocean(md,'Transient','runtimename',false);
1018+
1019+% %eval(['!mpiexec -np ' int2str(md.cluster.np) ' ' md.cluster.codepath '/issm_ocean.exe TransientSolution ' pwd ' ' md.miscellaneous.name ' ']);
1020+% eval(['!mpiexec -np ' int2str(md.cluster.np) ' ' md.cluster.codepath '/issm_ocean.exe TransientSolution ' pwd ' ' md.miscellaneous.name ' : -np ' int2str(nPx*nPy) ' ./mitgcmuv']);
1021+end
1022+% }}}
1023+
1024+%Fields and tolerances to track changes
1025+fnm=['run/SHICE_fwFluxtave.0000000006.data'];
1026+melting_rate_1=readbin(fnm,[Nx Ny]);
1027+fnm=['run/SHICE_fwFluxtave.0000000012.data'];
1028+melting_rate_2=readbin(fnm,[Nx Ny]);
1029+fnm=['run/SHICE_fwFluxtave.0000000018.data'];
1030+melting_rate_3=readbin(fnm,[Nx Ny]);
1031+fnm=['run/SHICE_fwFluxtave.0000000024.data'];
1032+melting_rate_4=readbin(fnm,[Nx Ny]);
1033+field_names ={'Base1','Melting1','Vx2','Vy2','Thickness2','Base2','MaskGroundediceLevelset2','FloatingiceMeltingRate2',...
1034+ 'Melting2','Vx3','Vy3','Thickness3','Base3','MaskGroundediceLevelset3','FloatingiceMeltingRate3',...
1035+ 'Melting3','Vx4','Vy4','Thickness4','Base4','MaskGroundediceLevelset4','FloatingiceMeltingRate4','Melting4'};
1036+field_tolerances={2e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,...
1037+ 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13,...
1038+ 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13 };
1039+field_values={...
1040+ (md.results.TransientSolution(1).Base),...
1041+ (melting_rate_1(:)),...
1042+ (md.results.TransientSolution(2).Vx),...
1043+ (md.results.TransientSolution(2).Vy),...
1044+ (md.results.TransientSolution(2).Thickness),...
1045+ (md.results.TransientSolution(2).Base),...
1046+ (md.results.TransientSolution(2).MaskGroundediceLevelset),...
1047+ (md.results.TransientSolution(2).BasalforcingsFloatingiceMeltingRate),...
1048+ (melting_rate_2(:)),...
1049+ (md.results.TransientSolution(3).Vx),...
1050+ (md.results.TransientSolution(3).Vy),...
1051+ (md.results.TransientSolution(3).Thickness),...
1052+ (md.results.TransientSolution(3).Base),...
1053+ (md.results.TransientSolution(3).MaskGroundediceLevelset),...
1054+ (md.results.TransientSolution(3).BasalforcingsFloatingiceMeltingRate),...
1055+ (melting_rate_3(:)),...
1056+ (md.results.TransientSolution(4).Vx),...
1057+ (md.results.TransientSolution(4).Vy),...
1058+ (md.results.TransientSolution(4).Thickness),...
1059+ (md.results.TransientSolution(4).Base),...
1060+ (md.results.TransientSolution(4).MaskGroundediceLevelset),...
1061+ (md.results.TransientSolution(4).BasalforcingsFloatingiceMeltingRate),...
1062+ (melting_rate_4(:)),...
1063+ };
Note: See TracBrowser for help on using the repository browser.