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

../trunk-jpl/test/NightlyRun/test4001.m

@@ +1 @@
+%ISSM/MITgcm coupled set-up
+%
+%Script control parameters
+steps=1:15;
+final_time=1;
+
+%Organizer
+mkdir Models
+org=organizer('repository','Models/','prefix','IceOcean.2013#','steps',steps);
+
+presentdirectory=pwd; 
+addpath([pwd '/../MITgcm/tools']);
+
+% {{{ Parameters: 
+if perform(org,'Parameters'),
+        Nx=20; %number of longitude cells
+        Ny=40; %number of latitude cells
+        Nz=30; %number of MITgcm vertical cells
+        nPx=2; %number of MITgcm processes to use in x direction
+        nPy=4; %number of MITgcm processes to use in y direction
+        xgOrigin=0; %origin of longitude
+        ygOrigin=-80; %origin of latitude
+        dLong=.25; %longitude grid spacing
+        dLat=.05; %latitude grid spacing
+        delZ=30; %thickness of vertical levels
+        icefront_position_ratio=.75; 
+        ice_thickness=100;
+        rho_ice=917;
+        rho_water=1028.14;
+        di=rho_ice/rho_water;
+        
+        % MITgcm initial and lateral boundary conditions
+        iniSalt  = 34.4; % initial salinity (PSU)
+        iniTheta = -1.9; % initial potential temperature (deg C)
+        obcSalt  = 34.4; % open boundary salinity (PSU)
+        obcTheta =  1.0; % open boundary potential temperature (deg C)
+        mlDepth  = 120.; % mixed layer depth (m)
+        mlSalt   = 33.4; % open boundary salinity (PSU)
+        mlTheta  = -1.9; % open boundary potential temperature (deg C)
+        obcUvel  = -0.1; % open boundary velocity (m/s)        
+
+        MITgcmDeltaT=600; % MITgcm time step in seconds
+        y2s=31536000; % year to seconds conversion, i.e., seconds per year
+
+        % start_time, final_time, and time_step
+        start_time=0; % in decimal years
+        time_step=1/12; % coupling interval in decimal years
+        async_step_MITgcm_multiplier=1/30; % used to reduce run time for MItgcm
+        
+        % bedrock/bathymetry
+        hmax=1000;
+        trough_depth=200;
+        deltah=300;
+        sea_level=1095;
+
+        % issm settings:
+        higherorder=0;
+        numlayers=10;
+
+        savedata(org, Nx, Ny, nPx, nPy, Nz, dLong, dLat, delZ, xgOrigin, ...
+                 ygOrigin, icefront_position_ratio, ice_thickness, rho_ice, ...
+                 rho_water, di, hmax, trough_depth, deltah, sea_level, ...
+                 iniSalt, iniTheta, obcSalt, obcTheta, mlDepth, mlSalt, ...
+                 mlTheta, obcUvel, start_time, time_step, MITgcmDeltaT, y2s,...
+                 higherorder,numlayers,async_step_MITgcm_multiplier);
+end
+% }}}
+% {{{ Bathymetry: 
+if perform(org,'Bathymetry'),
+
+    loaddata(org,'Parameters');
+    %create lat,long
+    lat=(ygOrigin+dLat/2):dLat:(ygOrigin+Ny*dLat);
+    long=(xgOrigin+dLong/2):dLong:(xgOrigin+Nx*dLong);
+    [lat long]=meshgrid(lat,long);
+
+    longmin=min(long(:));
+    longmax=max(long(:));
+    latmin=min(lat(:));
+    latmax=max(lat(:));
+
+    %create bedrock/bathymetry:
+    bedrock=zeros(Nx,Ny);
+    bedrock=hmax-deltah*tanh(pi*(2*(lat-latmin)./(latmax-latmin)-1))+ ...
+            trough_depth*cos(2*pi*long./(longmax-longmin));
+
+    %save bathymetry file for MITgcm
+    bathymetry=bedrock-sea_level;
+    savedata(org,lat,long,bathymetry);
+
+end
+% }}}
+% {{{ IceSheetGeometry: 
+if perform(org,'IceSheetGeometry'),
+        
+        loaddata(org,'Parameters');
+        loaddata(org,'Bathymetry');
+        latmin=min(lat(:));
+        latmax=max(lat(:));
+
+        %put ice_thickness constant layer of ice over the bathymetry, unless it floats: 
+        s=size(bathymetry);
+        thickness=ice_thickness*ones(s);
+
+        %figure out ice shelf: 
+        pos=find(-di*thickness>bathymetry);
+        iceshelf_mask=zeros(s);
+        iceshelf_mask(pos)=1;
+
+        ice_mask=ones(s);
+        pos=find((lat-latmin)/(latmax-latmin)>(icefront_position_ratio));
+        ice_mask(pos)=0;
+        iceshelf_mask(pos)=0;
+
+        %compute draft of ice shelf: 
+        draft=bathymetry;
+        pos=find(iceshelf_mask);
+        draft(pos)=-di*thickness(pos);
+        pos=find(~ice_mask); 
+        draft(pos)=0;
+
+        savedata(org,ice_mask,iceshelf_mask,draft,thickness);
+end
+% }}}
+
+%Configure MITgcm
+% {{{ GetMITgcm: 
+if perform(org,'GetMITgcm'),
+  %system([pwd '/../MITgcm/get_mitgcm.sh']);
+end
+% }}}
+% {{{ BuildMITgcm: 
+if perform(org,'BuildMITgcm'),
+
+    %load data: 
+    loaddata(org,'Parameters');
+
+    %specify computational grid in SIZE.h
+    fidi=fopen('../MITgcm/code/SIZE.h.bak','r');
+    fido=fopen('../MITgcm/code/SIZE.h','w');
+    tline = fgetl(fidi);
+    fprintf(fido,'%s\n',tline);
+    while 1
+        tline = fgetl(fidi);
+        if ~ischar(tline), break, end
+        %do the change here: 
+        if strcmpi(tline,'     &           sNx =  20,'),
+            fprintf(fido,'%s%i%s\n','     &           sNx =  ',round(Nx/nPx),',');
+            continue;
+        end
+        if strcmpi(tline,'     &           sNy =  20,'),
+            fprintf(fido,'%s%i%s\n','     &           sNy =  ',round(Ny/nPy),',');
+            continue;
+        end
+        if strcmpi(tline,'     &           nPx =   1,'),
+            fprintf(fido,'%s%i%s\n','     &           nPx = ',nPx,',');
+            continue;
+        end
+        if strcmpi(tline,'     &           nPy =   2,'),
+            fprintf(fido,'%s%i%s\n','     &           nPy = ',nPy,',');
+            continue;
+        end
+        fprintf(fido,'%s\n',tline);
+    end
+    %close  files
+    fclose(fidi);
+    fclose(fido);
+
+        system('./build_mitgcm.sh generic');
+end
+% }}}
+% {{{ RunUncoupledMITgcm: 
+if perform(org,'RunUncoupledMITgcm'),
+
+    %load data: 
+    loaddata(org,'Parameters');
+    loaddata(org,'Bathymetry');
+    loaddata(org,'IceSheetGeometry');
+    endtime = round(MITgcmDeltaT * ...
+     floor(time_step*y2s*async_step_MITgcm_multiplier/MITgcmDeltaT));
+
+    % {{{ prepare MITgcm 
+    % rename previous run directory and create new one
+    if exist ('run.old')
+        !\rm -rf run.old
+    end
+    if exist ('run')
+        !\mv run run.old
+    end
+    !\mkdir run
+    !\cp ../MITgcm/build/mitgcmuv run
+    !\cp ../MITgcm/input/* run
+    
+    %load data: 
+    loaddata(org,'Parameters');
+
+    % initial salinity
+    S=iniSalt*ones(Nx,Ny,Nz);
+    writebin('run/Salt.bin',S);
+    
+    % initial temperature
+    T=iniTheta*ones(Nx,Ny,Nz);
+    writebin('run/Theta.bin',T);
+    
+    % initial velocity
+    Z=zeros(Nx,Ny,Nz);
+    writebin('run/Uvel.bin',Z);
+    writebin('run/Vvel.bin',Z);
+    
+    % initial sea surface height
+    Z=zeros(Nx,Ny);
+    writebin('run/Etan.bin',Z);
+
+    % salinity boundary conditions
+    S=obcSalt*ones(Ny,Nz);
+    thk=delZ*ones(Nz,1);
+    bot=cumsum(thk);
+    ik=find(bot<=mlDepth);
+    S(:,ik)=mlSalt;
+    writebin('run/OBs.bin',S);
+
+    % temperature boundary conditions
+    T=obcTheta*ones(Ny,Nz);
+    T(:,ik)=mlTheta;
+    writebin('run/OBt.bin',T);
+
+    % zonal velocity boundary conditions
+    U=obcUvel*ones(Ny,Nz);
+    writebin('run/OBu.bin',U);
+    
+    % zero boundary conditions
+    Z=zeros(Ny,Nz);
+    writebin('run/zeros.bin',Z);
+
+    % build parameter file data.obcs
+    fidi=fopen('../MITgcm/input/data.obcs','r');
+    fido=fopen('run/data.obcs','w');
+    tline = fgetl(fidi);
+    fprintf(fido,'%s\n',tline);
+    while 1
+        tline = fgetl(fidi);
+        if ~ischar(tline), break, end
+        %do the change here: 
+        if strcmpi(tline,' OB_Iwest = 40*1,'),
+            fprintf(fido,'%s%i%s\n',' OB_Iwest = ',Ny,'*1,');
+            continue;
+        end
+        if strcmpi(tline,' OB_Ieast = 40*-1,'),
+            fprintf(fido,'%s%i%s\n',' OB_Ieast = ',Ny,'*-1,');
+            continue;
+        end
+        fprintf(fido,'%s\n',tline);
+    end
+    %close  files
+    fclose(fidi);
+    fclose(fido);
+
+    %save bathymetry and bedrock in run directory
+    writebin('run/bathymetry.bin',bathymetry);
+    writebin('run/icetopo.bin',draft);
+    % }}}
+
+    %start looping:  
+    for t=start_time:time_step:final_time,
+        disp(['Year: ' num2str(t)])
+        % {{{ generate MITgcm parameter file data 
+        fidi=fopen('../MITgcm/input/data','r');
+        fido=fopen('run/data','w');
+        tline = fgetl(fidi);
+        fprintf(fido,'%s\n',tline);
+        while 1
+            tline = fgetl(fidi);
+            if ~ischar(tline), break, end
+            %do the change here: 
+            if strcmpi(tline,' xgOrigin = 0.0,'),
+                fprintf(fido,'%s%i%s\n',' xgOrigin = ',xgOrigin,',');
+                continue;
+            end
+            if strcmpi(tline,' ygOrigin = -80.0,'),
+                fprintf(fido,'%s%i%s\n',' ygOrigin = ',ygOrigin,',');
+                continue;
+            end
+            if strcmpi(tline,' delX = 20*0.25,'),
+                fprintf(fido,'%s%i*%g%s\n',' delX = ',Nx,dLong,',');
+                continue;
+            end
+            if strcmpi(tline,' delY = 20*0.25,'),
+                fprintf(fido,'%s%i*%g%s\n',' delY = ',Ny,dLat,',');
+                continue;
+            end
+            if strcmpi(tline,' delZ = 30*30.0,'),
+                fprintf(fido,'%s%i*%g%s\n',' delZ = ',Nz,delZ,',');
+                continue;
+            end
+            if strcmpi(tline,' endTime=2592000.,'),
+                fprintf(fido,'%s%i%s\n',' endTime= ',endtime,',');
+                continue;
+            end
+            if strcmpi(tline,' deltaT=1200.0,'),
+                fprintf(fido,'%s%i%s\n',' deltaT= ',MITgcmDeltaT,',');
+                continue;
+            end
+            if strcmpi(tline,' pChkptFreq=2592000.,'),
+                fprintf(fido,'%s%i%s\n',' pChkptFreq= ',endtime,',');
+                continue;
+            end
+            if strcmpi(tline,' taveFreq=2592000.,'),
+                fprintf(fido,'%s%i%s\n',' taveFreq= ',endtime,',');
+                continue;
+            end
+            if strcmpi(tline,' rhoConst=1030.,'),
+                fprintf(fido,'%s%i%s\n',' rhoConst= ',rho_water,',');
+                continue;
+            end
+            if strcmpi(tline,' rhoNil=1030.,'),
+                fprintf(fido,'%s%i%s\n',' rhoNil= ',rho_water,',');
+                continue;
+            end
+            fprintf(fido,'%s\n',tline);
+        end
+        %close  files
+        fclose(fidi);
+        fclose(fido);
+        % }}}
+        % {{{ generate initial MITgcm conditions
+
+        ds=round(endtime/MITgcmDeltaT);
+        if t>start_time
+            % Read pickup file
+            fnm=['run/pickup.' myint2str(ds,10) '.data'];
+            U=readbin(fnm,[Nx Ny Nz],1,'real*8',0);
+            V=readbin(fnm,[Nx Ny Nz],1,'real*8',1);
+            T=readbin(fnm,[Nx Ny Nz],1,'real*8',2);
+            S=readbin(fnm,[Nx Ny Nz],1,'real*8',3);
+            E=readbin(fnm,[Nx Ny],1,'real*8',8*Nz);
+            writebin('run/Salt.bin' ,S);
+            writebin('run/Theta.bin',T);
+            writebin('run/Uvel.bin' ,U);
+            writebin('run/Vvel.bin' ,V);
+            writebin('run/Etan.bin' ,E);
+        end
+
+        % }}}
+        % {{{ system call to run MITgcm
+        cd run
+        eval(['!mpirun -np ' int2str(nPx*nPy) ' ./mitgcmuv']);
+        ts=round((t+time_step)*y2s/MITgcmDeltaT);
+        eval(['!\mv STDERR.0000 STDERR_' myint2str(ts,10) '.data'])
+        eval(['!\mv STDOUT.0000 STDOUT_' myint2str(ts,10) '.data'])
+        eval(['!\cp hFacC.data hFacC_' myint2str(ts,10) '.data'])
+        eval(['!\cp icetopo.bin icetopo_' myint2str(ts,10) '.data'])
+        for fld={'S','T','U','V','Eta', ...
+                 'SHICE_heatFluxtave','SHICE_fwFluxtave'}
+            eval(['!\mv ' fld{1} '.' myint2str(ds,10) '.data ' ...
+                  fld{1} '_' myint2str(ts,10) '.data'])
+        end
+        cd ..
+        % }}}
+    end
+end
+% }}}
+
+%Configure ISSM
+% {{{ CreateMesh: 
+if perform(org,'CreateMesh'),
+        
+        loaddata(org,'Parameters');
+        loaddata(org,'Bathymetry');
+        loaddata(org,'IceSheetGeometry');
+
+        %create model: 
+        md=model();
+        
+        %Grab lat,long from MITgcm: 
+        lat=lat(:);
+        long=long(:);
+
+        %project lat,long: 
+        [x,y]=ll2xy(lat,long,-1);
+
+        index=[];
+
+        %  C  D
+        %  A  B 
+        for j=1:Ny-1,
+                for i=1:Nx-1, 
+                        A=(j-1)*Nx+i;
+                        B=(j-1)*Nx+i+1;
+                        C=j*Nx+i;
+                        D=j*Nx+i+1;
+                        index(end+1,:)=[A B C];
+                        index(end+1,:)=[C B D];
+                end
+        end
+
+        %fill mesh and model: 
+        md=meshconvert(md,index,x,y);
+        md.mesh.lat=lat;
+        md.mesh.long=long;
+
+        savemodel(org,md);
+
+end
+% }}}
+% {{{ MeshGeometry: 
+if perform(org,'MeshGeometry'),
+        
+        loaddata(org,'Parameters');
+        loaddata(org,'CreateMesh');
+        loaddata(org,'Bathymetry');
+        loaddata(org,'IceSheetGeometry');
+
+        %transfer to vertices: 
+        bathymetry=bathymetry(:);
+        iceshelf_mask=iceshelf_mask(:);
+        ice_mask=ice_mask(:);
+        thickness=thickness(:);
+        draft=draft(:);
+
+        %start filling some of the fields 
+        md.geometry.bed=bathymetry;
+        md.geometry.thickness=thickness;
+        md.geometry.base=md.geometry.bed; 
+        pos=find(iceshelf_mask); md.geometry.base(pos)=draft(pos);
+        md.geometry.surface=md.geometry.base+md.geometry.thickness;
+
+        %nothing passes icefront: 
+        pos=find(~ice_mask);
+        md.geometry.thickness(pos)=1;
+        md.geometry.surface(pos)=(1-di)*md.geometry.thickness(pos);
+        md.geometry.base(pos)=-di*md.geometry.thickness(pos);
+
+        %level sets: 
+        md.mask.groundedice_levelset=-ones(md.mesh.numberofvertices,1);
+        md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
+
+        pos=find(ice_mask); md.mask.ice_levelset(pos)=-1; 
+        pos=find(~iceshelf_mask & ice_mask); md.mask.groundedice_levelset(pos)=1;
+
+        %identify edges of grounded ice: 
+        groundedice_levelset=md.mask.groundedice_levelset;
+        for i=1:md.mesh.numberofelements,
+                m=groundedice_levelset(md.mesh.elements(i,:));
+                if abs(sum(m))~=3,
+                        pos=find(m==1); md.mask.groundedice_levelset(md.mesh.elements(i,pos))=0;
+                end
+        end
+
+        %identify edges of ice: 
+        ice_levelset=md.mask.ice_levelset;
+        for i=1:md.mesh.numberofelements,
+                m=ice_levelset(md.mesh.elements(i,:));
+                if abs(sum(m))~=3,
+                        pos=find(m==-1); md.mask.ice_levelset(md.mesh.elements(i,pos))=0;
+                end
+        end
+
+        savemodel(org,md);
+end
+% }}}
+% {{{ ParameterizeIce: 
+if perform(org,'ParameterizeIce'),
+        
+        loaddata(org,'Parameters');
+        loaddata(org,'CreateMesh');
+        loaddata(org,'MeshGeometry');
+
+        %miscellaneous
+        md.miscellaneous.name='IceOcean';
+
+        %initial velocity: 
+        md.initialization.vx=zeros(md.mesh.numberofvertices,1);
+        md.initialization.vy=zeros(md.mesh.numberofvertices,1);
+        md.initialization.vz=zeros(md.mesh.numberofvertices,1);
+
+        %friction: 
+        md.friction.coefficient=30*ones(md.mesh.numberofvertices,1);
+        pos=find(md.mask.groundedice_levelset<=0);
+        md.friction.coefficient(pos)=0;
+        
+        md.friction.p=ones(md.mesh.numberofelements,1);
+        md.friction.q=ones(md.mesh.numberofelements,1);
+
+        %temperatures and surface mass balance:
+        md.initialization.temperature=(273.15-20)*ones(md.mesh.numberofvertices,1);
+        md.initialization.pressure=md.materials.rho_ice*md.constants.g*(md.geometry.surface-md.geometry.base);
+        md.surfaceforcings.mass_balance = [1*ones(md.mesh.numberofvertices,1); 1];
+
+        %Flow law 
+        md.materials.rheology_B=paterson(md.initialization.temperature);
+        md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+        md.damage.D=zeros(md.mesh.numberofvertices,1);
+        md.damage.spcdamage=NaN*ones(md.mesh.numberofvertices,1);
+        
+        % {{{ Deal with boundary conditions: we have the level sets right, just a matter of getting
+
+        %the spcs going.  
+        md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);
+        md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
+        md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);
+        md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
+        md.stressbalance.loadingforce=0*ones(md.mesh.numberofvertices,3);
+        md.masstransport.spcthickness=NaN*ones(md.mesh.numberofvertices,1); 
+
+        %deal with water: 
+        pos=find(md.mask.ice_levelset>0); 
+        md.stressbalance.spcvx(pos)=0;
+        md.stressbalance.spcvy(pos)=0;
+        md.stressbalance.spcvz(pos)=0;
+        md.masstransport.spcthickness(pos)=0;
+
+        %get some flux at the ice divide: 
+        pos=find(md.mesh.lat==min(md.mesh.lat));
+        md.stressbalance.spcvy(pos)=200;
+
+        %deal with boundaries, excluding icefront: 
+        vertex_on_boundary=zeros(md.mesh.numberofvertices,1);
+        vertex_on_boundary(md.mesh.segments(:,1:2))=1;
+        pos=find(vertex_on_boundary & md.mask.groundedice_levelset<=0);
+        md.stressbalance.spcvx(pos)=md.initialization.vx(pos);
+        md.stressbalance.spcvy(pos)=md.initialization.vy(pos);
+        md.stressbalance.spcvz(pos)=md.initialization.vz(pos);
+        md.masstransport.spcthickness(pos)=md.geometry.thickness(pos);
+        % }}}
+
+        md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+        md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
+        md.thermal.spctemperature=[md.initialization.temperature; 1]; %impose observed temperature on surface
+        md.basalforcings.geothermalflux=.064*ones(md.mesh.numberofvertices,1);
+
+        %flow equations: 
+        md=setflowequation(md,'SSA','all');
+
+        savemodel(org,md);
+end
+% }}}
+% {{{ Extrude: 
+if perform(org,'Extrude'),
+        
+        loaddata(org,'Parameters');
+        loaddata(org,'ParameterizeIce');
+
+        if higherorder,
+                md=extrude(md,numlayers,3);
+                md=setflowequation(md,'HO','all');
+
+                %water needs to be spc'd: 
+                pos=find(md.mask.ice_levelset>0); 
+                md.stressbalance.spcvx(pos)=0;
+                md.stressbalance.spcvy(pos)=0;
+                md.stressbalance.spcvz(pos)=0;
+                md.masstransport.spcthickness(pos)=0;
+        end
+
+        savemodel(org,md);
+end
+% }}}
+% {{{ RunUncoupledISSM: 
+if perform(org,'RunUncoupledISSM'),
+        
+        loaddata(org,'Parameters');
+        loaddata(org,'Extrude');
+
+        %timestepping: 
+        md.timestepping.final_time=final_time;
+        md.timestepping.time_step=time_step;
+        md.timestepping.time_adapt=0;
+        md.transient.isgroundingline=1;
+        md.transient.isthermal=0;
+        md.groundingline.migration='SubelementMigration2';
+
+        md.cluster=generic('name',oshostname(),'np',2);
+        md=solve(md,TransientSolutionEnum);
+
+        savemodel(org,md);
+end
+% }}}
+
+%Run MITgcm/ISSM
+% {{{ RunCoupledMITgcmISSM: 
+if perform(org,'RunCoupledMITgcmISSM'),
+
+        %load data: 
+        loaddata(org,'Parameters');
+        loaddata(org,'Extrude');
+        loaddata(org,'Bathymetry');
+        loaddata(org,'IceSheetGeometry');
+        endtime = round(MITgcmDeltaT * ...
+         floor(time_step*y2s*async_step_MITgcm_multiplier/MITgcmDeltaT));
+
+        % {{{ prepare MITgcm 
+        % rename previous run directory and create new one
+        if exist ('run.old')
+            !\rm -rf run.old
+        end
+        if exist ('run')
+            !\mv run run.old
+        end
+        !\mkdir run
+        !\cp ../MITgcm/build/mitgcmuv run
+        !\cp ../MITgcm/input/* run
+
+        %load data: 
+        loaddata(org,'Parameters');
+
+        % initial salinity
+        S=iniSalt*ones(Nx,Ny,Nz);
+        writebin('run/Salt.bin',S);
+
+        % initial temperature
+        T=iniTheta*ones(Nx,Ny,Nz);
+        writebin('run/Theta.bin',T);
+
+        % initial velocity
+        Z=zeros(Nx,Ny,Nz);
+        writebin('run/Uvel.bin',Z);
+        writebin('run/Vvel.bin',Z);
+
+        % initial sea surface height
+        Z=zeros(Nx,Ny);
+        writebin('run/Etan.bin',Z);
+
+        % salinity boundary conditions
+        S=obcSalt*ones(Ny,Nz);
+        thk=delZ*ones(Nz,1);
+        bot=cumsum(thk);
+        ik=find(bot<=mlDepth);
+        S(:,ik)=mlSalt;
+        writebin('run/OBs.bin',S);
+
+        % temperature boundary conditions
+        T=obcTheta*ones(Ny,Nz);
+        T(:,ik)=mlTheta;
+        writebin('run/OBt.bin',T);
+
+        % zonal velocity boundary conditions
+        U=obcUvel*ones(Ny,Nz);
+        writebin('run/OBu.bin',U);
+
+        % zero boundary conditions
+        Z=zeros(Ny,Nz);
+        writebin('run/zeros.bin',Z);
+
+        % build parameter file data.obcs
+        fidi=fopen('../MITgcm/input/data.obcs','r');
+        fido=fopen('run/data.obcs','w');
+        tline = fgetl(fidi);
+        fprintf(fido,'%s\n',tline);
+        while 1
+            tline = fgetl(fidi);
+            if ~ischar(tline), break, end
+            %do the change here: 
+            if strcmpi(tline,' OB_Iwest = 40*1,'),
+                fprintf(fido,'%s%i%s\n',' OB_Iwest = ',Ny,'*1,');
+                continue;
+            end
+            if strcmpi(tline,' OB_Ieast = 40*-1,'),
+                fprintf(fido,'%s%i%s\n',' OB_Ieast = ',Ny,'*-1,');
+                continue;
+            end
+            fprintf(fido,'%s\n',tline);
+        end
+        %close  files
+        fclose(fidi);
+        fclose(fido);
+
+        %save bathymetry in MITgcm run directory
+        writebin('run/bathymetry.bin',bathymetry);
+        % }}}
+            
+        % {{{ ISSM settings:
+
+        setenv('DYLD_LIBRARY_PATH', '/usr/local/gfortran/lib') 
+        %timestepping: 
+        md.timestepping.start_time=start_time;
+        md.timestepping.final_time=final_time;
+        md.timestepping.time_step=time_step;
+        md.timestepping.time_adapt=0;
+        md.cluster=generic('name',oshostname(),'np',2);
+        md.results.TransientSolution.Base=md.geometry.base;
+        md.transient.isgroundingline=1;
+        md.transient.isthermal=0;
+        md.groundingline.migration='SubelementMigration2';
+
+        % }}}
+
+        %start looping:
+        results=md.results;
+
+        for t=start_time:time_step:final_time
+            disp(['Year: ' num2str(t)])
+
+            %send draft from ISSM to MITgcm:
+            draft=md.results.TransientSolution(end).Base;
+            pos=find(md.mask.ice_levelset>0); draft(pos)=0;
+            if md.mesh.dimension==3, 
+                %collapse onto bottom layer: 
+                draft=project2d(md,draft,1);
+            end
+            if t>start_time
+                old_draft=readbin('run/icetopo.bin',[Nx,Ny]);
+            end
+            writebin('run/icetopo.bin',draft);
+
+            % {{{ generate MITgcm parameter file data 
+            fidi=fopen('../MITgcm/input/data','r');
+            fido=fopen('run/data','w');
+            tline = fgetl(fidi);
+            fprintf(fido,'%s\n',tline);
+            while 1
+                tline = fgetl(fidi);
+                if ~ischar(tline), break, end
+                %do the change here: 
+                if strcmpi(tline,' xgOrigin = 0.0,'),
+                    fprintf(fido,'%s%i%s\n',' xgOrigin = ',xgOrigin,',');
+                    continue;
+                end
+                if strcmpi(tline,' ygOrigin = -80.0,'),
+                    fprintf(fido,'%s%i%s\n',' ygOrigin = ',ygOrigin,',');
+                    continue;
+                end
+                if strcmpi(tline,' delX = 20*0.25,'),
+                    fprintf(fido,'%s%i*%g%s\n',' delX = ',Nx,dLong,',');
+                    continue;
+                end
+                if strcmpi(tline,' delY = 20*0.25,'),
+                    fprintf(fido,'%s%i*%g%s\n',' delY = ',Ny,dLat,',');
+                    continue;
+                end
+                if strcmpi(tline,' delZ = 30*30.0,'),
+                    fprintf(fido,'%s%i*%g%s\n',' delZ = ',Nz,delZ,',');
+                    continue;
+                end
+                if strcmpi(tline,' endTime=2592000.,'),
+                    fprintf(fido,'%s%i%s\n',' endTime= ',endtime,',');
+                    continue;
+                end
+                if strcmpi(tline,' deltaT=1200.0,'),
+                    fprintf(fido,'%s%i%s\n',' deltaT= ',MITgcmDeltaT,',');
+                    continue;
+                end
+                if strcmpi(tline,' pChkptFreq=2592000.,'),
+                    fprintf(fido,'%s%i%s\n',' pChkptFreq= ',endtime,',');
+                    continue;
+                end
+                if strcmpi(tline,' taveFreq=2592000.,'),
+                    fprintf(fido,'%s%i%s\n',' taveFreq= ',endtime,',');
+                    continue;
+                end
+                if strcmpi(tline,' rhoConst=1030.,'),
+                    fprintf(fido,'%s%i%s\n',' rhoConst= ',rho_water,',');
+                    continue;
+                end
+                if strcmpi(tline,' rhoNil=1030.,'),
+                    fprintf(fido,'%s%i%s\n',' rhoNil= ',rho_water,',');
+                    continue;
+                end
+                fprintf(fido,'%s\n',tline);
+            end
+            %close  files
+            fclose(fidi);
+            fclose(fido);
+            % }}}
+
+            % {{{ generate initial MITgcm conditions
+            ds=round(endtime/MITgcmDeltaT);
+            if t>start_time
+                % Read pickup file
+                fnm=['run/pickup.' myint2str(ds,10) '.data'];
+                U=readbin(fnm,[Nx Ny Nz],1,'real*8',0);
+                V=readbin(fnm,[Nx Ny Nz],1,'real*8',1);
+                T=readbin(fnm,[Nx Ny Nz],1,'real*8',2);
+                S=readbin(fnm,[Nx Ny Nz],1,'real*8',3);
+                E=readbin(fnm,[Nx Ny],1,'real*8',8*Nz);
+
+                % find indices of locations where ice shelf retreated
+                h=readbin('run/hFacC.data',[Nx Ny Nz]);
+                msk=sum(h,3);
+                msk(find(msk))=1;
+                [iw jw]=find(msk); % horizontal indices where there is water
+                tmp=reshape(draft,[Nx,Ny])-old_draft;
+                tmp(find(tmp<0))=0;
+                [im jm]=find(tmp); % horizontal indices where there is melt
+                
+                % Extrapolate T/S to locations where ice shelf retreated
+                for i=1:length(im)
+
+                    % first try vertical extrapolation
+                    in=find(h(im(i),jm(i),:));
+                    if length(in)>0;
+                        S(im(i),jm(i),1:min(in)  ) = S(im(i),jm(i),min(in));
+                        T(im(i),jm(i),1:min(in)  ) = T(im(i),jm(i),min(in));
+                        continue
+                    end
+
+                    % if not succesful, use closest neighbor horizontal extrapolation
+                    [y c]=min((iw-im(i)).^2+(jw-jm(i)).^2);
+                    salt=squeeze(S(iw(c),jw(c),:)); % salinity profile of closest neighbor
+                    temp=squeeze(T(iw(c),jw(c),:)); % salinity profile of closest neighbor
+                    in=find(h(iw(c),jw(c),:));
+                    salt(1:min(in))=salt(min(in));
+                    temp(1:min(in))=temp(min(in));
+                    salt(max(in):end)=salt(max(in));
+                    temp(max(in):end)=temp(max(in));
+                    S(im(i),jm(i),:)=salt;
+                    T(im(i),jm(i),:)=temp;
+                end
+
+                % Write initial conditions
+                writebin('run/Salt.bin' ,S);
+                writebin('run/Theta.bin',T);
+                writebin('run/Uvel.bin' ,U);
+                writebin('run/Vvel.bin' ,V);
+                writebin('run/Etan.bin' ,E);
+            end
+            % }}}
+
+            % {{{ system call to run MITgcm 
+            cd run
+            eval(['!mpirun -np ' int2str(nPx*nPy) ' ./mitgcmuv']);
+            ts=round((t+time_step)*y2s/MITgcmDeltaT);
+            eval(['!\mv STDERR.0000 STDERR_' myint2str(ts,10) '.data'])
+            eval(['!\mv STDOUT.0000 STDOUT_' myint2str(ts,10) '.data'])
+            eval(['!\cp hFacC.data hFacC_' myint2str(ts,10) '.data'])
+            eval(['!\cp icetopo.bin icetopo_' myint2str(ts,10) '.data'])
+            for fld={'S','T','U','V','Eta', ...
+                     'SHICE_heatFluxtave','SHICE_fwFluxtave'}
+                eval(['!\mv ' fld{1} '.' myint2str(ds,10) '.data ' ...
+                      fld{1} '_' myint2str(ts,10) '.data'])
+            end
+            cd ..
+            % }}}
+
+            %get melting rates from MITgcm
+            %upward fresh water flux (kg/m^2/s):
+            fnm=['run/SHICE_fwFluxtave_' myint2str(ts,10) '.data'];
+            melting_rate=readbin(fnm,[Nx Ny]);
+
+            %send averaged melting rate to ISSM
+            %downward fresh water flux (m/y):
+            melting_rate=-melting_rate(:)*y2s/rho_ice;
+                if md.mesh.dimension==3,
+                        md.basalforcings.floatingice_melting_rate=project3d(md,'vector',melting_rate,'type','node');
+                else
+                        md.basalforcings.floatingice_melting_rate=melting_rate;
+                end
+
+            % {{{ run ISSM and recover results 
+
+            md.timestepping.start_time=t;
+            md.timestepping.final_time=t+time_step;;
+                md=solve(md,TransientSolutionEnum);
+
+                base=md.results.TransientSolution(end).Base;
+                thickness=md.results.TransientSolution(end).Thickness;
+                if md.mesh.dimension==3,
+                        md.mesh.z=base+thickness./md.geometry.thickness.*(md.mesh.z-md.geometry.bed);
+                        md.initialization.vz=md.results.TransientSolution(end).Vz;
+                end
+                md.geometry.base=base;
+                md.geometry.thickness=thickness;
+                md.geometry.surface=md.geometry.base+md.geometry.thickness;
+                md.initialization.vx=md.results.TransientSolution(end).Vx;
+                md.initialization.vy=md.results.TransientSolution(end).Vy;
+                md.initialization.vel=md.results.TransientSolution(end).Vel;
+                md.initialization.pressure=md.results.TransientSolution(end).Pressure;
+                md.mask.groundedice_levelset=md.results.TransientSolution(end).MaskGroundediceLevelset;
+                md.results.TransientSolution(end).FloatingiceMeltingRate=md.basalforcings.floatingice_melting_rate;
+                
+                %save these results in the model, otherwise, they'll be wiped out
+                results(end+1)=md.results;
+
+                % }}}
+
+        end
+
+        md.results=results;
+        savemodel(org,md);
+end
+% }}}