source: issm/oecreview/Archive/18296-19100/ISSM-18532-18533.diff@ 19102

Last change on this file since 19102 was 19102, checked in by Mathieu Morlighem, 10 years ago

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