Index: /issm/trunk-jpl/test/NightlyRun/test4003.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test4003.m (revision 23987) +++ /issm/trunk-jpl/test/NightlyRun/test4003.m (revision 23988) @@ -3,5 +3,5 @@ % %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,301 +326,27 @@ % {{{ 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,18 +486,12 @@ % }}} - 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,2 +499,291 @@ 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 +% }}}