%Test Name: IceOceanCoupling Dan Goldberg'd setup
%ISSM/MITgcm coupled set-up
%
%Script control parameters
steps=1:11;
final_time=1/365;

%To download and recompile MITgcm from scratch:
!rm -rf ${ISSM_DIR}/test/MITgcm/install_dngoldberg
!rm -rf ${ISSM_DIR}/test/MITgcm/build/*
!rm -rf Models

%Organizer
!mkdir Models
org=organizer('repository','Models','prefix','IceOcean.','steps',steps);

presentdirectory=pwd;

% {{{ Parameters:
if perform(org,'Parameters'),
	Nx=3; %number of longitude cells
	Ny=200; %number of latitude cells
	Nz=90; %number of MITgcm vertical cells
	nPx=1; %number of MITgcm processes to use in x direction
	nPy=1; %number of MITgcm processes to use in y direction
	xgOrigin=1; %origin of longitude
	ygOrigin=1; %origin of latitude
	dLong=1; %longitude grid spacing
	dLat=1; %latitude grid spacing
	delZ=10; %thickness of vertical levels
	icefront_position_ratio=.75;
	ice_thickness=1000;
	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 and time_step
	start_time=0; % in decimal years
	time_step=1/(365*24); % coupling interval in decimal years
	async_step_MITgcm_multiplier=1; % used to reduce run time for MItgcm

	% bedrock/bathymetry
	hmax=1000;
	trough_depth=200;
	deltah=300;
	sea_level=1095;

	% issm settings:
	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,...
		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_dngoldberg.sh']);
end
% }}}
% {{{ BuildMITgcm:
if perform(org,'BuildMITgcm'),
    system(['../MITgcm/build_4003.sh generic ' pwd '/../MITgcm']);
end
% }}}
addpath(recursivepath([pwd '/../MITgcm']));
% {{{ RunUncoupledMITgcm:
if perform(org,'RunUncoupledMITgcm'),
    % 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/install_dngoldberg/verification/shelfice_remeshing/input/* run
    !\cp ../MITgcm/input_4003/* run
    !\cp ../MITgcm/input_4003/data_uncoupled run/data
    !\cp ../MITgcm/input_4003/eedata_uncoupled run/eedata
    cd run
    eval(['!mpirun -np ' int2str(nPx*nPy) ' ./mitgcmuv']);
	 cd ..
end
% }}}

%Configure ISSM
% {{{ CreateMesh:
if perform(org,'CreateMesh'),

	%create model:
	md=model();

	%Grab lat,long from MITgcm:
	long=readbin('run/XG.data',[3 200]);
	long=[long long(:,end)]; long=[long; -105.1250*ones(1,size(long,2))];
	lat=readbin('run/YG.data',[3 200]);
	lat=[lat -73.8832*ones(size(lat,1),1)]; lat=[lat; lat(end,:)];

	%project lat,long:
	[x,y]=ll2xy(lat(:),long(:),-1);

	Nx=size(lat,1); Ny=size(lat,2)
	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,'CreateMesh');

	%transfer to vertices:
	bathymetry=readbin('run/bathy.box',[3 200],1,'real*8');
	bathymetry=[bathymetry bathymetry(:,end)]; bathymetry=[bathymetry(1,:); bathymetry];
	iceshelf_mask=-1*ones(size(bathymetry));
	ice_mask=readbin('run/hmask3.box',[3 200],1,'real*8');
	ice_mask=[ice_mask ice_mask(:,end)]; ice_mask=[ice_mask(1,:); ice_mask];
	thickness=readbin('run/h0.bin',[3 200],1,'real*8');
	thickness=[thickness thickness(:,end)]; thickness=[thickness; thickness(end,:)];

	%start filling some of the fields
	md.geometry.bed=bathymetry(:);
	md.geometry.thickness=thickness(:);
	md.geometry.base=-917/1028*md.geometry.thickness;
	md.geometry.surface=md.geometry.base+md.geometry.thickness;

	%nothing passes icefront:
	pos=find((~ice_mask(:) & ice_mask(:)~=0) | thickness(:)==0);
	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.ocean_levelset=iceshelf_mask(:);
	md.mask.ice_levelset=-ones(md.mesh.numberofvertices,1);

	pos=find(~ice_mask(:) & thickness(:)==0); md.mask.ice_levelset(pos)=1;

	savemodel(org,md);
end
% }}}
% {{{ ParameterizeIce:
if perform(org,'ParameterizeIce'),

	loaddata(org,'MeshGeometry');

	%miscellaneous
	md.miscellaneous.name='test4003';

	%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=0*ones(md.mesh.numberofvertices,1);
	pos=find(md.mask.ocean_levelset>0);
	md.friction.coefficient(pos)=5;
	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-22)*ones(md.mesh.numberofvertices,1);
	md.initialization.pressure=md.materials.rho_ice*md.constants.g*(md.geometry.surface-md.geometry.base);
	md.smb.mass_balance = 0*ones(md.mesh.numberofvertices,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);

	%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);

	%get some flux at the ice divide:
	pos=find(md.mesh.lat==min(md.mesh.lat));
	md.masstransport.spcthickness(pos)=md.geometry.thickness(pos);
	md.stressbalance.spcvx(pos)=-800;
	md.stressbalance.spcvy(pos)=0;

	%deal with boundaries, excluding icefront:
	pos=find(md.mesh.long==min(md.mesh.long) | md.mesh.long==max(md.mesh.long));
	md.stressbalance.spcvy(pos)=0;

	point1=find(md.mesh.y==min(md.mesh.y)); point2=find(md.mesh.x==max(md.mesh.x));
	costheta=(md.mesh.x(point2)-md.mesh.x(point1))/sqrt((md.mesh.x(point2)-md.mesh.x(point1)).^2+(md.mesh.y(point2)-md.mesh.y(point1)).^2);
	sintheta=(md.mesh.y(point2)-md.mesh.y(point1))/sqrt((md.mesh.x(point2)-md.mesh.x(point1)).^2+(md.mesh.y(point2)-md.mesh.y(point1)).^2);
	md.stressbalance.referential(:,1:3)=repmat([costheta,sintheta,0],md.mesh.numberofvertices,1);
	md.stressbalance.referential(:,4:6)=repmat([-sintheta,costheta,0],md.mesh.numberofvertices,1);

	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
% }}}
% {{{ RunUncoupledISSM:
if perform(org,'RunUncoupledISSM'),

	loaddata(org,'Parameters');
	loaddata(org,'ParameterizeIce');

	%timestepping:
	md.timestepping.final_time=100;
	md.timestepping.time_step=0.5;
	md.transient.isgroundingline=0;
	md.transient.isthermal=0;
	md.groundingline.migration='SubelementMigration';
	md.groundingline.melt_interpolation='NoMeltOnPartiallyFloating';
	md.groundingline.friction_interpolation='SubelementFriction2';
	md.masstransport.stabilization=0;

	md.cluster=generic('name',oshostname(),'np',2);
	md=solve(md,'Transient');

	savemodel(org,md);

	plotmodel(md,'data',md.results.TransientSolution(end).Vel,'data',md.results.TransientSolution(end).Thickness)
end
% }}}

%Run MITgcm/ISSM
% {{{ RunCoupledMITgcmISSM:
if perform(org,'RunCoupledMITgcmISSM'),

	loaddata(org,'Parameters');
	loaddata(org,'ParameterizeIce');

	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=0;
		  md.masstransport.requested_outputs={'default','BasalforcingsFloatingiceMeltingRate'};

        % }}}
        % {{{ 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/install_dngoldberg/verification/shelfice_remeshing/input/* run
        !\cp ../MITgcm/input_4003/* run
        % }}}
		  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.transient.requested_outputs={'default','MaskOceanLevelset'};

		  md=solveiceocean(md,'Transient','runtimename',false);

%	%eval(['!mpiexec -np ' int2str(md.cluster.np) ' ' md.cluster.codepath '/issm_ocean.exe TransientSolution ' pwd ' ' md.miscellaneous.name ' ']);
%	eval(['!mpiexec -np ' int2str(md.cluster.np) ' ' md.cluster.codepath '/issm_ocean.exe TransientSolution ' pwd ' ' md.miscellaneous.name ' : -np ' int2str(nPx*nPy) ' ./mitgcmuv']);
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.ocean_levelset=md.results.TransientSolution(end).MaskOceanLevelset;
		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
% }}}

%Fields and tolerances to track changes
fnm=['run/SHICE_fwFlux.0000011400.data'];
melting_rate_1=readbin(fnm,[Nx Ny]);
fnm=['run/SHICE_fwFlux.0000011450.data'];
melting_rate_2=readbin(fnm,[Nx Ny]);
fnm=['run/SHICE_fwFlux.0000011500.data'];
melting_rate_3=readbin(fnm,[Nx Ny]);
fnm=['run/SHICE_fwFlux.0000011550.data'];
melting_rate_4=readbin(fnm,[Nx Ny]);
field_names     ={'Base1','Melting1','Vx2','Vy2','Thickness2','Base2','MaskOceanLevelset2','FloatingiceMeltingRate2',...
	'Melting2','Vx3','Vy3','Thickness3','Base3','MaskOceanLevelset3','FloatingiceMeltingRate3',...
	'Melting3','Vx4','Vy4','Thickness4','Base4','MaskOceanLevelset4','FloatingiceMeltingRate4','Melting4'};
field_tolerances={2e-13,1e-13,7e-13,7e-13,1e-13,1e-13,1e-13,1e-12,...
	1e-13, 1e-12, 1e-12, 1e-13, 1e-13, 1e-13, 1e-12,...
	1e-13, 1e-12, 1e-12, 1e-13, 1e-13, 1e-13, 1e-12, 1e-13 };
field_values={...
	(md.results.TransientSolution(1).Base),...
	(melting_rate_1(:)),...
	(md.results.TransientSolution(2).Vx),...
	(md.results.TransientSolution(2).Vy),...
	(md.results.TransientSolution(2).Thickness),...
	(md.results.TransientSolution(2).Base),...
	(md.results.TransientSolution(2).MaskOceanLevelset),...
	(md.results.TransientSolution(2).BasalforcingsFloatingiceMeltingRate),...
	(melting_rate_2(:)),...
	(md.results.TransientSolution(3).Vx),...
	(md.results.TransientSolution(3).Vy),...
	(md.results.TransientSolution(3).Thickness),...
	(md.results.TransientSolution(3).Base),...
	(md.results.TransientSolution(3).MaskOceanLevelset),...
	(md.results.TransientSolution(3).BasalforcingsFloatingiceMeltingRate),...
	(melting_rate_3(:)),...
	(md.results.TransientSolution(4).Vx),...
	(md.results.TransientSolution(4).Vy),...
	(md.results.TransientSolution(4).Thickness),...
	(md.results.TransientSolution(4).Base),...
	(md.results.TransientSolution(4).MaskOceanLevelset),...
	(md.results.TransientSolution(4).BasalforcingsFloatingiceMeltingRate),...
	(melting_rate_4(:)),...
	};
