Index: /issm/trunk-jpl/test/MITgcm/code_4003/cpl_issm.F =================================================================== --- /issm/trunk-jpl/test/MITgcm/code_4003/cpl_issm.F (revision 26718) +++ /issm/trunk-jpl/test/MITgcm/code_4003/cpl_issm.F (revision 26719) @@ -67,5 +67,9 @@ & MPI_COMM_MODEL,mpiRC) _END_MASTER( myThid ) -C print*, 'Ocean received CouplingTime: ', CouplingTime +#ifdef ALLOW_DEBUG + IF ( debugLevel .GE. debLevD ) THEN + print*, 'Ocean received CouplingTime: ', CouplingTime + ENDIF +#endif C I2. MITgcm sends grid size (NX and NY) @@ -135,5 +139,9 @@ call MPI_Recv(IceModelTime,1,MPI_DOUBLE,0,10001001, & toissmcomm,mpistatus,mpiRC) -C print*, 'Ocean received IceModelTime: ', IceModelTime +#ifdef ALLOW_DEBUG + IF ( debugLevel .GE. debLevD ) THEN + print*, 'Ocean received IceModelTime: ', IceModelTime + ENDIF +#endif _END_MASTER( myThid ) ENDIF @@ -176,6 +184,9 @@ ENDIF CALL BAR2( myThid ) -C print*,'Done Sending shelficeFreshWaterFlux array.' - +#ifdef ALLOW_DEBUG + IF ( debugLevel .GE. debLevD ) THEN + print*,'Done Sending shelficeFreshWaterFlux array.' + ENDIF +#endif ENDIF C End recurring step C3. Index: /issm/trunk-jpl/test/MITgcm/input_4003/data =================================================================== --- /issm/trunk-jpl/test/MITgcm/input_4003/data (revision 26718) +++ /issm/trunk-jpl/test/MITgcm/input_4003/data (revision 26719) @@ -46,5 +46,5 @@ #- readBinaryPrec=64, - writeBinaryPrec=64, +#writeBinaryPrec=64, useSingleCpuIO=.TRUE., #globalFiles=.TRUE., Index: /issm/trunk-jpl/test/NightlyRun/test4003.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test4003.m (revision 26718) +++ /issm/trunk-jpl/test/NightlyRun/test4003.m (revision 26719) @@ -3,7 +3,5 @@ % %Script control parameters -steps=1:12; -steps=1:6; -final_time=1/365; +steps=1:11; %To download and recompile MITgcm from scratch: @@ -20,53 +18,29 @@ % {{{ 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=8; %number of MITgcm processes to use in y direction - xgOrigin=0; %origin of longitude - ygOrigin=-75.5; %origin of latitude - dLong=.125; %longitude grid spacing - dLat=dLong/16; %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=450; % 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 - - % 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); + 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=8; % number of MITgcm processes to use in y direction + xgOrigin=0; % origin of longitude + ygOrigin=-75.5; % origin of latitude + dLong=.125; % longitude grid spacing + dLat=dLong/16; % latitude grid spacing + delZ=10; % thickness of vertical levels (m) + gravity= 9.81; % gravity (m^2/s) + rho_ice=917; + rho_water=1030; + di=rho_ice/rho_water; + prec = 'real*8'; % precision of MITgcm input binary files + + % bathymetry and ice sheet geometry + H = -900; % water depth in the ice shelf cavity + Hmin = -600; % deepest point of cavern + Hmax = -300; % shallowest point of cavern + jEnd = Ny*3/4; % where ice-shelf ends + j2 = jEnd+1; + + savedata(org, Nx, Ny, nPx, nPy, Nz, dLong, dLat, delZ, xgOrigin, ygOrigin, ... + rho_ice, rho_water, di, H, Hmin, Hmax, jEnd, j2, gravity, prec); end % }}} @@ -75,12 +49,4 @@ loaddata(org,'Parameters'); - - v0 = 2e3; - h0 = 800; - - hfacMin = 0.2; - - eos = 'jmd95z'; - prec = 'real*8'; %create lat,lon @@ -89,27 +55,8 @@ long = xgOrigin+[0:Nx-1]*dLong; lonc = long+dLong/2; + [lat lon]=meshgrid(latc,lonc); zC=-delZ*([1:Nz]-0.5); zF=-delZ*[0:Nz]; - - lat=(ygOrigin+dLat/2):dLat:(ygOrigin+Ny*dLat); - lon=(xgOrigin+dLong/2):dLong:(xgOrigin+Nx*dLong); - [lat lon]=meshgrid(lat,lon); - - longmin=min(lon(:)); - longmax=max(lon(:)); - latmin=min(lat(:)); - latmax=max(lat(:)); - -% Gravity -gravity= 9.81; -rhoConst= 1030; - -% Nominal depth of model (meters) -H = -900; %water depth in the ice shelf cavity -Hmin = -600; % deepest point of cavern -Hmax = -300; % shallowest point of cavern -jEnd = Ny*3/4; % where ice-shelf ends -dHdy = (Hmax-Hmin)/dLat/(jEnd-2); %Slope of ice shelf - + %create bathymetry: bathymetry = ones(Nx,Ny)*H; @@ -119,71 +66,4 @@ savedata(org,lat,lon,bathymetry); -j2=jEnd+1; -hIce=bathymetry; -for i=1:Nx - hIce(i,:)=Hmin+dHdy*[-1:Ny-2]*dLat; - hIce(i,1)=0; hIce(i,j2:Ny)=0; -end - -var=([1:Ny]-2)/(jEnd-2); -dMdt_fy=bathymetry; -for i=1:Nx - dMdt_fy(i,:)=-cos(pi*var); - dMdt_fy(i,1)=0; dMdt_fy(i,j2:Ny)=0; -end - -regMsk=ones(Nx,Ny); -regMsk(:,1)=0; regMsk(:,j2:Ny)=2; - -%- rate of change due to ice-stream dynamics -rateDyn=rhoConst*0.1/3600; sfx='r02'; - -dMdt=rateDyn*dMdt_fy; - -dz = delZ*ones(1,Nz); -zgp1 = [0,cumsum(dz)]; -zc = .5*(zgp1(1:end-1)+zgp1(2:end)); -zg = zgp1(1:end-1); -dz = diff(zgp1); - -T_sfc = -1.9; -T_bot = 2; -del_T = (T_bot - T_sfc)/(59*delZ); -tref=0*dz; -for k = 1:Nz; - tref(k) = T_sfc + del_T*((k-20)*delZ); - tref(k)= max(T_sfc,min(tref(k),T_bot)); -end - -S_sfc = 34.2; -S_bot = 34.7; -del_S = (S_bot - S_sfc)/(59*delZ); -sref=0*dz; -for k = 1:Nz; - sref(k) = S_sfc + del_S*((k-20)*delZ); - sref(k)= max(S_sfc,min(sref(k),S_bot)); -end -pEOS=-rhoConst*gravity*zC; % in Pa -pEOS=pEOS*1.e-4; % in dBar -rhoAn=densjmd95(sref,tref,pEOS); -rhoAn=rhoAn-rhoConst; - -pF=-rhoConst*gravity*zF*1.e-4; % in dBar -rhoUp=densjmd95(sref,tref,pF(2:end)); -rhoDw=densjmd95(sref,tref,pF(1:Nz)); -dRho=rhoUp(1:Nz-1)-rhoDw(2:Nz); -NSq=-gravity*dRho/delZ/rhoConst; - -mnV=min(NSq); MxV=max(NSq); Avr=mean(NSq); - -zax=[1:Nz]; - -v1=2.5e-2; -var=1+Nz-2*zax; var=var/(Nz-1); -vobc=v1*var; - -rhoAvr=rhoConst-1.345; -mIce=-rhoAvr*hIce; - end % }}} @@ -191,31 +71,24 @@ 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); + loaddata(org,'Parameters'); + loaddata(org,'Bathymetry'); + latmin=min(lat(:)); + latmax=max(lat(:)); + + dHdy = (Hmax-Hmin)/dLat/(jEnd-2); %Slope of ice shelf + draft=bathymetry; + for i=1:Nx + draft(i,:)=Hmin+dHdy*[-1:Ny-2]*dLat; + draft(i,1)=0; draft(i,j2:Ny)=0; + end + + ice_mask=ones(Nx,Ny); + ice_mask(:,j2:Ny)=0; + iceshelf_mask=ice_mask; + iceshelf_mask(:,1)=0; + thickness=abs(draft)/di; + thickness(:,1)=thickness(:,2)+draft(:,2); + + savedata(org,ice_mask,iceshelf_mask,draft,thickness); end % }}} @@ -235,4 +108,9 @@ % {{{ RunUncoupledMITgcm: if perform(org,'RunUncoupledMITgcm'), + + loaddata(org,'Parameters'); + loaddata(org,'Bathymetry'); + loaddata(org,'IceSheetGeometry'); + % rename previous run directory and create new one if exist ('run.old') @@ -246,27 +124,53 @@ !\cp ../MITgcm/input_4003/* run !\cp ../MITgcm/input_4003/eedata_uncoupled run/eedata - cd run - - namF='bathy_flat.bin'; + + % {{{ Construct MITgcm binary input files + namF='run/bathy_flat.bin'; fid=fopen(namF,'w','b'); fwrite(fid,bathymetry,prec);fclose(fid); - namF='shelficeTopo.Lin.bin'; - fid=fopen(namF,'w','b'); fwrite(fid,hIce,prec);fclose(fid); - - namF='under_Ice_mask.bin'; - fid=fopen(namF,'w','b'); fwrite(fid,regMsk,prec);fclose(fid); - - namF=sprintf('%s.%s.%s','shelfice_dMdt',sfx,'bin'); + namF='run/shelficeTopo.Lin.bin'; + fid=fopen(namF,'w','b'); fwrite(fid,draft,prec);fclose(fid); + + var=([1:Ny]-2)/(jEnd-2); + dMdt_fy=zeros(Nx,Ny); + for i=1:Nx + dMdt_fy(i,:)=-cos(pi*var); + dMdt_fy(i,1)=0; dMdt_fy(i,j2:Ny)=0; + end + %- rate of change due to ice-stream dynamics + rateDyn=rho_water*0.1/3600; sfx='r02'; + dMdt=rateDyn*dMdt_fy; + namF=sprintf('%s.%s.%s','run/shelfice_dMdt',sfx,'bin'); fid=fopen(namF,'w','b'); fwrite(fid,dMdt,prec);fclose(fid); - namF='temp_obc.bin'; + T_sfc = -1.9; + T_bot = 2; + del_T = (T_bot - T_sfc)/(59*delZ); + tref=zeros(1,Nz); + for k = 1:Nz; + tref(k) = T_sfc + del_T*((k-20)*delZ); + tref(k)= max(T_sfc,min(tref(k),T_bot)); + end + namF='run/temp_obc.bin'; tref=[tref; tref; tref]; fid=fopen(namF,'w','b'); fwrite(fid,tref,prec);fclose(fid); - namF='salt_obc.bin'; + S_sfc = 34.2; + S_bot = 34.7; + del_S = (S_bot - S_sfc)/(59*delZ); + sref=zeros(1,Nz); + for k = 1:Nz; + sref(k) = S_sfc + del_S*((k-20)*delZ); + sref(k)= max(S_sfc,min(sref(k),S_bot)); + end + namF='run/salt_obc.bin'; sref=[sref; sref; sref]; fid=fopen(namF,'w','b'); fwrite(fid,sref,prec);fclose(fid); - namF='vVel_obc.bin'; + zax=[1:Nz]; + v1=2.5e-2; + var=1+Nz-2*zax; var=var/(Nz-1); + vobc=v1*var; + namF='run/vVel_obc.bin'; vobc=[vobc; vobc; vobc]; fid=fopen(namF,'w','b'); fwrite(fid,vobc,prec);fclose(fid); @@ -276,5 +180,5 @@ var(i,j,:)=tref(1,:); end, end - namF='temp_ini.bin'; + namF='run/temp_ini.bin'; fid=fopen(namF,'w','b'); fwrite(fid,var,prec);fclose(fid); @@ -282,10 +186,14 @@ var(i,j,:)=sref(1,:); end, end - namF='salt_ini.bin'; + namF='run/salt_ini.bin'; fid=fopen(namF,'w','b'); fwrite(fid,var,prec);fclose(fid); - namF='shelficeMass.Lin.bin'; + rhoAvr=rho_water-1.345; + mIce=-rhoAvr*draft; + namF='run/shelficeMass.Lin.bin'; fid=fopen(namF,'w','b'); fwrite(fid,mIce,prec);fclose(fid); - + % }}} + + cd run eval(['!mpirun -np ' int2str(nPx*nPy) ' ./mitgcmuv']); cd .. @@ -297,37 +205,37 @@ if perform(org,'CreateMesh'), - %create model: - md=model(); - - %Grab lat,long from MITgcm: - long=readbin('run/XG.data',[Nx Ny]); - long=[long long(:,end)]; long=[long; -105.1250*ones(1,size(long,2))]; - lat=readbin('run/YG.data',[Nx Ny]); - 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); + loaddata(org,'Parameters'); + + %create model: + md=model(); + + %Grab lat,long from MITgcm: + long=readbin('run/XG.data',[Nx Ny]); + long=[long long(:,end)]; long=[long; 2*long(Nx,:)-long(Nx-1,:)]; + lat=readbin('run/YG.data',[Nx Ny]); + lat=[lat 2*lat(:,Ny)-lat(:,Ny-1)]; lat=[lat; lat(end,:)]; + + %project lat,long: + [x,y]=ll2xy(lat(:),long(:),-1); + index=[]; + % C D + % A B + for j=1:Ny, + for i=1:Nx, + A=(j-1)*(Nx+1)+i; + B=(j-1)*(Nx+1)+i+1; + C=j*(Nx+1)+i; + D=j*(Nx+1)+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 @@ -336,34 +244,35 @@ 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); + loaddata(org,'Parameters'); + loaddata(org,'CreateMesh'); + loaddata(org,'Bathymetry'); + loaddata(org,'IceSheetGeometry'); + + %transfer to vertices: + bathymetry=[bathymetry bathymetry(:,end)]; bathymetry=[bathymetry(1,:); bathymetry]; + ice_mask=[ice_mask ice_mask(:,end)]; ice_mask=[ice_mask(1,:); ice_mask]; + iceshelf_mask=[iceshelf_mask iceshelf_mask(:,end)]; iceshelf_mask=[iceshelf_mask(1,:); iceshelf_mask]; + 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 % }}} @@ -463,340 +372,120 @@ 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/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 -% }}} + loaddata(org,'Parameters'); + loaddata(org,'Bathymetry'); + loaddata(org,'IceSheetGeometry'); + loaddata(org,'ParameterizeIce'); + md=loadmodel(org,'RunUncoupledISSM'); + + % {{{ 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/input_4003/* run + !\cp ../MITgcm/input_4003/eedata_uncoupled run/eedata + + % {{{ Construct MITgcm binary input files + namF='run/bathy_flat.bin'; + fid=fopen(namF,'w','b'); fwrite(fid,bathymetry,prec);fclose(fid); + + namF='run/shelficeTopo.Lin.bin'; + fid=fopen(namF,'w','b'); fwrite(fid,draft,prec);fclose(fid); + + var=([1:Ny]-2)/(jEnd-2); + dMdt_fy=zeros(Nx,Ny); + for i=1:Nx + dMdt_fy(i,:)=-cos(pi*var); + dMdt_fy(i,1)=0; dMdt_fy(i,j2:Ny)=0; + end + %- rate of change due to ice-stream dynamics + rateDyn=rho_water*0.1/3600; sfx='r02'; + dMdt=rateDyn*dMdt_fy; + namF=sprintf('%s.%s.%s','run/shelfice_dMdt',sfx,'bin'); + fid=fopen(namF,'w','b'); fwrite(fid,dMdt,prec);fclose(fid); + + T_sfc = -1.9; + T_bot = 2; + del_T = (T_bot - T_sfc)/(59*delZ); + tref=zeros(1,Nz); + for k = 1:Nz; + tref(k) = T_sfc + del_T*((k-20)*delZ); + tref(k)= max(T_sfc,min(tref(k),T_bot)); + end + namF='run/temp_obc.bin'; + tref=[tref; tref; tref]; + fid=fopen(namF,'w','b'); fwrite(fid,tref,prec);fclose(fid); + + S_sfc = 34.2; + S_bot = 34.7; + del_S = (S_bot - S_sfc)/(59*delZ); + sref=zeros(1,Nz); + for k = 1:Nz; + sref(k) = S_sfc + del_S*((k-20)*delZ); + sref(k)= max(S_sfc,min(sref(k),S_bot)); + end + namF='run/salt_obc.bin'; + sref=[sref; sref; sref]; + fid=fopen(namF,'w','b'); fwrite(fid,sref,prec);fclose(fid); + + zax=[1:Nz]; + v1=2.5e-2; + var=1+Nz-2*zax; var=var/(Nz-1); + vobc=v1*var; + namF='run/vVel_obc.bin'; + vobc=[vobc; vobc; vobc]; + fid=fopen(namF,'w','b'); fwrite(fid,vobc,prec);fclose(fid); + + var=zeros(Nx,Ny,Nz); + for i=1:Nx, for j=1:Ny + var(i,j,:)=tref(1,:); + end, end + namF='run/temp_ini.bin'; + fid=fopen(namF,'w','b'); fwrite(fid,var,prec);fclose(fid); + + for i=1:Nx, for j=1:Ny + var(i,j,:)=sref(1,:); + end, end + namF='run/salt_ini.bin'; + fid=fopen(namF,'w','b'); fwrite(fid,var,prec);fclose(fid); + + rhoAvr=rho_water-1.345; + mIce=-rhoAvr*draft; + namF='run/shelficeMass.Lin.bin'; + fid=fopen(namF,'w','b'); fwrite(fid,mIce,prec);fclose(fid); + % }}} + % }}} + + md.timestepping.coupling_time=1/24/365; % 1 hour in decimal years + md.timestepping.time_step=1/24/365; % 1 hour in decimal years + md.timestepping.final_time=23/24/365; % hour 23 in decimal years + 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 +% }}}