Index: /issm/trunk-jpl/test/NightlyRun/test4003.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test4003.m (revision 23897) +++ /issm/trunk-jpl/test/NightlyRun/test4003.m (revision 23898) @@ -4,4 +4,5 @@ %Script control parameters steps=[1 2 3 4 5 6 7 8 9 10 11 12]; +steps=[ 7 8 9 10]; final_time=1/365; @@ -161,20 +162,17 @@ if perform(org,'CreateMesh'), - loaddata(org,'Parameters'); - loaddata(org,'Bathymetry'); - loaddata(org,'IceSheetGeometry'); - %create model: md=model(); %Grab lat,long from MITgcm: - lat=lat(:); - long=long(:); + long=readbin('run/XG.data',[3 200]); + long=[long long(:,end)]; long=[long; -105.0*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); - x=long; - y=lat; - + [x,y]=ll2xy(lat(:),long(:),-1); + + Nx=size(lat,1); Ny=size(lat,2) index=[]; % C D @@ -193,6 +191,6 @@ %fill mesh and model: md=meshconvert(md,index,x,y); - md.mesh.lat=lat; - md.mesh.long=long; + md.mesh.lat=lat(:); + md.mesh.long=long(:); savemodel(org,md); @@ -203,25 +201,23 @@ if perform(org,'MeshGeometry'), - loaddata(org,'Parameters'); loaddata(org,'CreateMesh'); - loaddata(org,'Bathymetry'); - loaddata(org,'IceSheetGeometry'); %transfer to vertices: - bathymetry=bathymetry(:); - iceshelf_mask=iceshelf_mask(:); - ice_mask=ice_mask(:); - thickness=thickness(:); - draft=draft(:); + 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=md.geometry.bed; - pos=find(iceshelf_mask); md.geometry.base(pos)=draft(pos); + 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); + pos=find((~ice_mask(:) & ice_mask(:)~=0) | thickness(:)==0); md.geometry.thickness(pos)=1; md.geometry.surface(pos)=(1-di)*md.geometry.thickness(pos); @@ -229,27 +225,8 @@ %level sets: - md.mask.groundedice_levelset=-ones(md.mesh.numberofvertices,1); - md.mask.ice_levelset=ones(md.mesh.numberofvertices,1); - - pos=find(ice_mask); md.mask.ice_levelset(pos)=-1; - pos=find(~iceshelf_mask & ice_mask); md.mask.groundedice_levelset(pos)=1; - - %identify edges of grounded ice: - groundedice_levelset=md.mask.groundedice_levelset; - for i=1:md.mesh.numberofelements, - m=groundedice_levelset(md.mesh.elements(i,:)); - if abs(sum(m))~=3, - pos=find(m==1); md.mask.groundedice_levelset(md.mesh.elements(i,pos))=0; - end - end - - %identify edges of ice: - ice_levelset=md.mask.ice_levelset; - for i=1:md.mesh.numberofelements, - m=ice_levelset(md.mesh.elements(i,:)); - if abs(sum(m))~=3, - pos=find(m==-1); md.mask.ice_levelset(md.mesh.elements(i,pos))=0; - end - end + md.mask.groundedice_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); @@ -259,10 +236,8 @@ if perform(org,'ParameterizeIce'), - loaddata(org,'Parameters'); - loaddata(org,'CreateMesh'); loaddata(org,'MeshGeometry'); %miscellaneous - md.miscellaneous.name='test4002'; + md.miscellaneous.name='test4003'; %initial velocity: @@ -281,5 +256,5 @@ 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 = [1*ones(md.mesh.numberofvertices,1); 1]; + md.smb.mass_balance = [0*ones(md.mesh.numberofvertices,1); 1]; %Flow law @@ -297,22 +272,19 @@ md.masstransport.spcthickness=NaN*ones(md.mesh.numberofvertices,1); - %deal with water: - pos=find(md.mask.ice_levelset>0); - md.stressbalance.spcvx(pos)=0; + %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; - md.stressbalance.spcvz(pos)=0; - md.masstransport.spcthickness(pos)=0; - - %get some flux at the ice divide: - pos=find(md.mesh.y==min(md.mesh.y)); - md.masstransport.spcthickness(pos)=md.geometry.thickness(pos); - md.stressbalance.spcvx(pos)=0; - md.stressbalance.spcvy(pos)=1500; %deal with boundaries, excluding icefront: - vertex_on_boundary=zeros(md.mesh.numberofvertices,1); - vertex_on_boundary(md.mesh.segments(:,1:2))=1; - pos=find(vertex_on_boundary & md.mask.groundedice_levelset<=0); - md.stressbalance.spcvx(pos)=0; + 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); @@ -335,5 +307,5 @@ %timestepping: md.timestepping.final_time=100; - md.timestepping.time_step=1; + md.timestepping.time_step=0.25; md.transient.isgroundingline=0; md.transient.isthermal=0; @@ -346,4 +318,6 @@ savemodel(org,md); + + plotmodel(md,'data',md.results.TransientSolution(end).Vel,'data',md.results.TransientSolution(end).Thickness) end % }}}