Changeset 23988

Timestamp:

06/05/19 10:08:12 (6 years ago)

Author:

seroussi

Message:

NEW: working on coupling test4003

File:

• issm/trunk-jpl/test/NightlyRun/test4003.m (modified) (4 diffs)

issm/trunk-jpl/test/NightlyRun/test4003.m ¶

@@ -3 +3 @@
 %
 %Script control parameters
-steps=[1 2 3 4 5 6 7 8 9 10 11 12];
+steps=[1 2 3 4 5 6 7 8 9 10 11];
 final_time=1/365;
 
@@ -326 +326 @@
 % {{{ RunCoupledMITgcmISSM: 
 if perform(org,'RunCoupledMITgcmISSM'),
-
-        %load data: 
-        loaddata(org,'Parameters');
-        loaddata(org,'ParameterizeIce');
-        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
-        !\cp ../MITgcm/input/eedata_uncoupled run/eedata
-
-        % 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.cluster=generic('name',oshostname(),'np',2);
-        md.results.TransientSolution.Base=md.geometry.base;
-        md.transient.isgroundingline=1;
-        md.transient.isthermal=0;
-        md.groundingline.migration='SubelementMigration';
-        md.groundingline.melt_interpolation='NoMeltOnPartiallyFloating';
-        md.groundingline.friction_interpolation='SubelementFriction2';
-
-        % }}}
-
-        %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 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(['!mpiexec -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;
-                 md.basalforcings.floatingice_melting_rate=melting_rate;
-
-            % {{{ run ISSM and recover results 
-
-            md.timestepping.start_time=t;
-            md.timestepping.final_time=t+time_step;;
-                md=solve(md,'Transient');
-
-                base=md.results.TransientSolution(end).Base;
-                thickness=md.results.TransientSolution(end).Thickness;
-                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
-% }}}
-% {{{ RunCoupledMITgcmISSM2: 
-if perform(org,'RunCoupledMITgcmISSM2'),
         
         loaddata(org,'Parameters');
         loaddata(org,'ParameterizeIce');
-        loaddata(org,'Bathymetry');
-        loaddata(org,'IceSheetGeometry');
+
+        md=loadmodel(org,'RunUncoupledISSM');
+
                 endtime = round(MITgcmDeltaT * floor(final_time*y2s/MITgcmDeltaT));
                 outputtime = round(MITgcmDeltaT * floor(time_step*y2s/MITgcmDeltaT));
 
+        % {{{ prepare ISSM: start from the steady-state
+
+                  md.geometry.base=md.results.TransientSolution(end).Base;
+                  md.geometry.surface=md.results.TransientSolution(end).Surface;
+                  md.geometry.thickness=md.results.TransientSolution(end).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.transient.isoceancoupling=1;
+                  md.transient.isgroundingline=1;
+                  md.masstransport.requested_outputs={'default','BasalforcingsFloatingiceMeltingRate'};
+
+        % }}}
         % {{{ prepare MITgcm 
         % rename previous run directory and create new one
@@ -760 +486 @@
         % }}}
                 
-        md.transient.isoceancoupling=1;
-        md.transient.isgroundingline=1;
-        md.groundingline.migration='None';
-        md.groundingline.melt_interpolation='SubelementMelt2';
-        md.groundingline.friction_interpolation='SubelementFriction2';
-        md.timestepping.coupling_time=time_step;
-   md.timestepping.time_step=time_step;
-        md.timestepping.final_time=final_time-time_step;
-        md.cluster.npocean=nPx*nPy;
-        md.cluster.np=2;
-        md.cluster.executionpath=[pwd '/run'];
-        md.masstransport.requested_outputs={'default','BasalforcingsFloatingiceMeltingRate'};
-
-        md=solveiceocean(md,'Transient','runtimename',false);
+                  md.timestepping.coupling_time=time_step;
+                  md.timestepping.time_step=time_step;
+                  md.timestepping.final_time=final_time-time_step;
+                  md.cluster.npocean=nPx*nPy;
+                  md.cluster.np=2;
+                  md.cluster.executionpath=[pwd '/run'];
+
+                  md=solveiceocean(md,'Transient','runtimename',false);
 
 %       %eval(['!mpiexec -np ' int2str(md.cluster.np) ' ' md.cluster.codepath '/issm_ocean.exe TransientSolution ' pwd ' ' md.miscellaneous.name ' ']);
@@ -779 +499 @@
 end
 % }}}
+% {{{ RunCoupledMITgcmISSMMatlabCoupling: 
+if perform(org,'RunCoupledMITgcmISSMMatlabCoupling'),
+
+        %load data: 
+        loaddata(org,'Parameters');
+        loaddata(org,'ParameterizeIce');
+        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
+        !\cp ../MITgcm/input/eedata_uncoupled run/eedata
+
+        % 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.cluster=generic('name',oshostname(),'np',2);
+        md.results.TransientSolution.Base=md.geometry.base;
+        md.transient.isgroundingline=1;
+        md.transient.isthermal=0;
+        md.groundingline.migration='SubelementMigration';
+        md.groundingline.melt_interpolation='NoMeltOnPartiallyFloating';
+        md.groundingline.friction_interpolation='SubelementFriction2';
+
+        % }}}
+
+        %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 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(['!mpiexec -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;
+                 md.basalforcings.floatingice_melting_rate=melting_rate;
+
+            % {{{ run ISSM and recover results 
+
+            md.timestepping.start_time=t;
+            md.timestepping.final_time=t+time_step;;
+                md=solve(md,'Transient');
+
+                base=md.results.TransientSolution(end).Base;
+                thickness=md.results.TransientSolution(end).Thickness;
+                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
+% }}}