Index: /issm/trunk-jpl/jenkins/ross-debian_linux-iceocean =================================================================== --- /issm/trunk-jpl/jenkins/ross-debian_linux-iceocean (revision 24455) +++ /issm/trunk-jpl/jenkins/ross-debian_linux-iceocean (revision 24456) @@ -10,4 +10,5 @@ --prefix=$ISSM_DIR \ --disable-static \ + --with-ocean \ --with-matlab-dir=$MATLAB_PATH \ --with-fortran-lib="-L/usr/lib/gcc/x86_64-linux-gnu/8 -lgfortran" \ @@ -44,5 +45,5 @@ mpich install-3.3.sh lapack install-3.8-linux.sh - petsc install-3.11-linux.sh + petsc install-3.7-linux.sh triangle install-linux64.sh chaco install.sh Index: /issm/trunk-jpl/test/NightlyRun/test4001.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test4001.m (revision 24455) +++ /issm/trunk-jpl/test/NightlyRun/test4001.m (revision 24456) @@ -6,68 +6,69 @@ final_time=1; -!rm -rf $ISSM_DIR/test/MITgcm/install -!rm -rf $ISSM_DIR/test/MITgcm/build/* -!rm $ISSM_DIR/test/MITgcm/code/SIZE.h +!rm -rf ${ISSM_DIR}/test/MITgcm/install +!rm -rf ${ISSM_DIR}/test/MITgcm/build/* +!rm -f ${ISSM_DIR}/test/MITgcm/code/SIZE.h +!rm -rf Models %Organizer -mkdir Models -org=organizer('repository','Models/','prefix','IceOcean.','steps',steps); - -presentdirectory=pwd; - -% {{{ Parameters: +!mkdir Models +org=organizer('repository','Models','prefix','IceOcean.','steps',steps); + +presentdirectory=pwd; + +% {{{ Parameters: if perform(org,'Parameters'), - Nx=20; %number of longitude cells - Ny=40; %number of latitude cells - Nz=30; %number of MITgcm vertical cells - nPx=2; %number of MITgcm processes to use in x direction - nPy=4; %number of MITgcm processes to use in y direction - xgOrigin=0; %origin of longitude - ygOrigin=-80; %origin of latitude - dLong=.25; %longitude grid spacing - dLat=.05; %latitude grid spacing - delZ=30; %thickness of vertical levels - icefront_position_ratio=.75; - ice_thickness=100; - 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, final_time, and time_step - start_time=0; % in decimal years - time_step=1/12; % coupling interval in decimal years - async_step_MITgcm_multiplier=1/30; % 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: + Nx=20; %number of longitude cells + Ny=40; %number of latitude cells + Nz=30; %number of MITgcm vertical cells + nPx=2; %number of MITgcm processes to use in x direction + nPy=4; %number of MITgcm processes to use in y direction + xgOrigin=0; %origin of longitude + ygOrigin=-80; %origin of latitude + dLong=.25; %longitude grid spacing + dLat=.05; %latitude grid spacing + delZ=30; %thickness of vertical levels + icefront_position_ratio=.75; + ice_thickness=100; + 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, final_time, and time_step + start_time=0; % in decimal years + time_step=1/12; % coupling interval in decimal years + async_step_MITgcm_multiplier=1/30; % 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'), @@ -94,47 +95,47 @@ end % }}} -% {{{ IceSheetGeometry: +% {{{ 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); + + 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: +% {{{ GetMITgcm: if perform(org,'GetMITgcm'), system([pwd '/../MITgcm/get_mitgcm.sh']); end % }}} -% {{{ BuildMITgcm: +% {{{ BuildMITgcm: if perform(org,'BuildMITgcm'), - %load data: + %load data: loaddata(org,'Parameters'); @@ -148,5 +149,5 @@ tline = fgetl(fidi); if ~ischar(tline), break, end - %do the change here: + %do the change here: if strcmpi(tline,' & sNx = 20,'), fprintf(fido,'%s%i%s\n',' & sNx = ',round(Nx/nPx),','); @@ -176,15 +177,15 @@ % }}} addpath(recursivepath([pwd '/../MITgcm'])); -% {{{ RunUncoupledMITgcm: +% {{{ RunUncoupledMITgcm: if perform(org,'RunUncoupledMITgcm'), - %load data: + %load data: loaddata(org,'Parameters'); loaddata(org,'Bathymetry'); loaddata(org,'IceSheetGeometry'); - endtime = round(MITgcmDeltaT * ... - floor(time_step*y2s*async_step_MITgcm_multiplier/MITgcmDeltaT)); - - % {{{ prepare MITgcm + 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') @@ -198,6 +199,6 @@ !\cp ../MITgcm/input/* run !\cp ../MITgcm/input/eedata_uncoupled run/eedata - - %load data: + + %load data: loaddata(org,'Parameters'); @@ -205,14 +206,14 @@ 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); @@ -235,5 +236,5 @@ U=obcUvel*ones(Ny,Nz); writebin('run/OBu.bin',U); - + % zero boundary conditions Z=zeros(Ny,Nz); @@ -248,5 +249,5 @@ tline = fgetl(fidi); if ~ischar(tline), break, end - %do the change here: + %do the change here: if strcmpi(tline,' OB_Iwest = 40*1,'), fprintf(fido,'%s%i%s\n',' OB_Iwest = ',Ny,'*1,'); @@ -268,8 +269,8 @@ % }}} - %start looping: + %start looping: for t=start_time:time_step:final_time, disp(['Year: ' num2str(t)]) - % {{{ generate MITgcm parameter file data + % {{{ generate MITgcm parameter file data fidi=fopen('../MITgcm/input/data','r'); fido=fopen('run/data','w'); @@ -279,5 +280,5 @@ tline = fgetl(fidi); if ~ischar(tline), break, end - %do the change here: + %do the change here: if strcmpi(tline,' xgOrigin = 0.0,'), fprintf(fido,'%s%i%s\n',' xgOrigin = ',xgOrigin,','); @@ -369,207 +370,207 @@ %Configure ISSM -% {{{ CreateMesh: +% {{{ CreateMesh: 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(:); - - %project lat,long: - [x,y]=ll2xy(lat,long,-1); - - 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: + + loaddata(org,'Parameters'); + loaddata(org,'Bathymetry'); + loaddata(org,'IceSheetGeometry'); + + %create model: + md=model(); + + %Grab lat,long from MITgcm: + lat=lat(:); + long=long(:); + + %project lat,long: + [x,y]=ll2xy(lat,long,-1); + + 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,'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(:); - - %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.surface=md.geometry.base+md.geometry.thickness; - - %nothing passes icefront: - pos=find(~ice_mask); - 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.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 - - savemodel(org,md); -end -% }}} -% {{{ ParameterizeIce: + + 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(:); + + %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.surface=md.geometry.base+md.geometry.thickness; + + %nothing passes icefront: + pos=find(~ice_mask); + 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.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 + + savemodel(org,md); +end +% }}} +% {{{ ParameterizeIce: if perform(org,'ParameterizeIce'), - - loaddata(org,'Parameters'); - loaddata(org,'CreateMesh'); - loaddata(org,'MeshGeometry'); - - %miscellaneous - md.miscellaneous.name='test4001'; - - %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=30*ones(md.mesh.numberofvertices,1); - pos=find(md.mask.groundedice_levelset<=0); - md.friction.coefficient(pos)=0; - 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-20)*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]; - - %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); - - %deal with water: - pos=find(md.mask.ice_levelset>0); - md.stressbalance.spcvx(pos)=0; - 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.lat==min(md.mesh.lat)); - md.stressbalance.spcvy(pos)=200; - - %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)=md.initialization.vx(pos); - md.stressbalance.spcvy(pos)=md.initialization.vy(pos); - md.stressbalance.spcvz(pos)=md.initialization.vz(pos); - md.masstransport.spcthickness(pos)=md.geometry.thickness(pos); - - 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: + + loaddata(org,'Parameters'); + loaddata(org,'CreateMesh'); + loaddata(org,'MeshGeometry'); + + %miscellaneous + md.miscellaneous.name='test4001'; + + %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=30*ones(md.mesh.numberofvertices,1); + pos=find(md.mask.groundedice_levelset<=0); + md.friction.coefficient(pos)=0; + 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-20)*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]; + + %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); + + %deal with water: + pos=find(md.mask.ice_levelset>0); + md.stressbalance.spcvx(pos)=0; + 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.lat==min(md.mesh.lat)); + md.stressbalance.spcvy(pos)=200; + + %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)=md.initialization.vx(pos); + md.stressbalance.spcvy(pos)=md.initialization.vy(pos); + md.stressbalance.spcvz(pos)=md.initialization.vz(pos); + md.masstransport.spcthickness(pos)=md.geometry.thickness(pos); + + 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=final_time; - md.timestepping.time_step=time_step; - md.transient.isgroundingline=1; - md.transient.isthermal=0; - md.groundingline.migration='SubelementMigration'; - md.groundingline.melt_interpolation='SubelementMelt2'; - md.groundingline.friction_interpolation='SubelementFriction2'; - - md.cluster=generic('name',oshostname(),'np',2); - md=solve(md,'Transient'); - - savemodel(org,md); + + loaddata(org,'Parameters'); + loaddata(org,'ParameterizeIce'); + + %timestepping: + md.timestepping.final_time=final_time; + md.timestepping.time_step=time_step; + md.transient.isgroundingline=1; + md.transient.isthermal=0; + md.groundingline.migration='SubelementMigration'; + md.groundingline.melt_interpolation='SubelementMelt2'; + md.groundingline.friction_interpolation='SubelementFriction2'; + + md.cluster=generic('name',oshostname(),'np',2); + md=solve(md,'Transient'); + + savemodel(org,md); end % }}} %Run MITgcm/ISSM -% {{{ RunCoupledMITgcmISSM: +% {{{ RunCoupledMITgcmISSM: if perform(org,'RunCoupledMITgcmISSM'), - %load data: - loaddata(org,'Parameters'); - loaddata(org,'ParameterizeIce'); - loaddata(org,'Bathymetry'); - loaddata(org,'IceSheetGeometry'); + %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 + % {{{ prepare MITgcm % rename previous run directory and create new one if exist ('run.old') @@ -584,5 +585,5 @@ !\cp ../MITgcm/input/eedata_uncoupled run/eedata - %load data: + %load data: loaddata(org,'Parameters'); @@ -633,5 +634,5 @@ tline = fgetl(fidi); if ~ischar(tline), break, end - %do the change here: + %do the change here: if strcmpi(tline,' OB_Iwest = 40*1,'), fprintf(fido,'%s%i%s\n',' OB_Iwest = ',Ny,'*1,'); @@ -651,31 +652,31 @@ 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='SubelementMelt2'; - md.groundingline.friction_interpolation='SubelementFriction2'; - - % }}} - - %start looping: - results=md.results; - - for t=start_time:time_step:final_time + + % {{{ 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='SubelementMelt2'; + 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; + %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]); @@ -683,49 +684,49 @@ 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); + % {{{ 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 + %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,' 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,','); @@ -736,12 +737,12 @@ continue; end - fprintf(fido,'%s\n',tline); - end - %close files - fclose(fidi); - fclose(fido); - % }}} - - % {{{ generate initial MITgcm conditions + fprintf(fido,'%s\n',tline); + end + %close files + fclose(fidi); + fclose(fido); + % }}} + + % {{{ generate initial MITgcm conditions ds=round(endtime/MITgcmDeltaT); if t>start_time @@ -762,5 +763,5 @@ 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) @@ -796,5 +797,5 @@ % }}} - % {{{ system call to run MITgcm + % {{{ system call to run MITgcm cd run eval(['!mpirun -np ' int2str(nPx*nPy) ' ./mitgcmuv']); @@ -812,41 +813,41 @@ % }}} - %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); + %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 % }}} @@ -862,32 +863,32 @@ melting_rate_4=readbin(fnm,[Nx Ny]); field_names ={'Base1','Melting1','Vx2','Vy2','Thickness2','Base2','MaskGroundediceLevelset2','FloatingiceMeltingRate2',... - 'Melting2','Vx3','Vy3','Thickness3','Base3','MaskGroundediceLevelset3','FloatingiceMeltingRate3',... - 'Melting3','Vx4','Vy4','Thickness4','Base4','MaskGroundediceLevelset4','FloatingiceMeltingRate4','Melting4'}; + 'Melting2','Vx3','Vy3','Thickness3','Base3','MaskGroundediceLevelset3','FloatingiceMeltingRate3',... + 'Melting3','Vx4','Vy4','Thickness4','Base4','MaskGroundediceLevelset4','FloatingiceMeltingRate4','Melting4'}; field_tolerances={2e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,... - 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13,... - 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13 }; + 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13,... + 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13 }; field_values={... - (md.results(1).TransientSolution(end).Base),... - (melting_rate_1(:)),... - (md.results(2).TransientSolution(end).Vx),... - (md.results(2).TransientSolution(end).Vy),... - (md.results(2).TransientSolution(end).Thickness),... - (md.results(2).TransientSolution(end).Base),... - (md.results(2).TransientSolution(end).MaskGroundediceLevelset),... - (md.results(2).TransientSolution(end).FloatingiceMeltingRate),... - (melting_rate_2(:)),... - (md.results(3).TransientSolution(end).Vx),... - (md.results(3).TransientSolution(end).Vy),... - (md.results(3).TransientSolution(end).Thickness),... - (md.results(3).TransientSolution(end).Base),... - (md.results(3).TransientSolution(end).MaskGroundediceLevelset),... - (md.results(3).TransientSolution(end).FloatingiceMeltingRate),... - (melting_rate_3(:)),... - (md.results(4).TransientSolution(end).Vx),... - (md.results(4).TransientSolution(end).Vy),... - (md.results(4).TransientSolution(end).Thickness),... - (md.results(4).TransientSolution(end).Base),... - (md.results(4).TransientSolution(end).MaskGroundediceLevelset),... - (md.results(4).TransientSolution(end).FloatingiceMeltingRate),... - (melting_rate_4(:)),... - }; + (md.results(1).TransientSolution(end).Base),... + (melting_rate_1(:)),... + (md.results(2).TransientSolution(end).Vx),... + (md.results(2).TransientSolution(end).Vy),... + (md.results(2).TransientSolution(end).Thickness),... + (md.results(2).TransientSolution(end).Base),... + (md.results(2).TransientSolution(end).MaskGroundediceLevelset),... + (md.results(2).TransientSolution(end).FloatingiceMeltingRate),... + (melting_rate_2(:)),... + (md.results(3).TransientSolution(end).Vx),... + (md.results(3).TransientSolution(end).Vy),... + (md.results(3).TransientSolution(end).Thickness),... + (md.results(3).TransientSolution(end).Base),... + (md.results(3).TransientSolution(end).MaskGroundediceLevelset),... + (md.results(3).TransientSolution(end).FloatingiceMeltingRate),... + (melting_rate_3(:)),... + (md.results(4).TransientSolution(end).Vx),... + (md.results(4).TransientSolution(end).Vy),... + (md.results(4).TransientSolution(end).Thickness),... + (md.results(4).TransientSolution(end).Base),... + (md.results(4).TransientSolution(end).MaskGroundediceLevelset),... + (md.results(4).TransientSolution(end).FloatingiceMeltingRate),... + (melting_rate_4(:)),... + }; Index: /issm/trunk-jpl/test/NightlyRun/test4002.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test4002.m (revision 24455) +++ /issm/trunk-jpl/test/NightlyRun/test4002.m (revision 24456) @@ -3,72 +3,73 @@ % %Script control parameters -steps=[1 2 3 4 5 6 7 8 9 10 11 12]; +steps=1:12; final_time=1/365; %To download and recompile MITgcm from scratch: -!rm -rf $ISSM_DIR/test/MITgcm/install -!rm -rf $ISSM_DIR/test/MITgcm/build/* -!rm $ISSM_DIR/test/MITgcm/code/SIZE.h +!rm -rf ${ISSM_DIR}/test/MITgcm/install +!rm -rf ${ISSM_DIR}/test/MITgcm/build/* +!rm -f ${ISSM_DIR}/test/MITgcm/code/SIZE.h +!rm -rf Models %Organizer -mkdir Models -org=organizer('repository','Models/','prefix','IceOcean.','steps',steps); +!mkdir Models +org=organizer('repository','Models','prefix','IceOcean.','steps',steps); presentdirectory=pwd; -% {{{ Parameters: +% {{{ Parameters: if perform(org,'Parameters'), - Nx=20; %number of longitude cells - Ny=40; %number of latitude cells - Nz=30; %number of MITgcm vertical cells - nPx=2; %number of MITgcm processes to use in x direction - nPy=4; %number of MITgcm processes to use in y direction - xgOrigin=0; %origin of longitude - ygOrigin=-80; %origin of latitude - dLong=.25; %longitude grid spacing - dLat=.05; %latitude grid spacing - delZ=30; %thickness of vertical levels - icefront_position_ratio=.75; - ice_thickness=100; - 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: + Nx=20; %number of longitude cells + Ny=40; %number of latitude cells + Nz=30; %number of MITgcm vertical cells + nPx=2; %number of MITgcm processes to use in x direction + nPy=4; %number of MITgcm processes to use in y direction + xgOrigin=0; %origin of longitude + ygOrigin=-80; %origin of latitude + dLong=.25; %longitude grid spacing + dLat=.05; %latitude grid spacing + delZ=30; %thickness of vertical levels + icefront_position_ratio=.75; + ice_thickness=100; + 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'), @@ -95,47 +96,47 @@ end % }}} -% {{{ IceSheetGeometry: +% {{{ 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); + + 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: +% {{{ GetMITgcm: if perform(org,'GetMITgcm'), system([pwd '/../MITgcm/get_mitgcm.sh']); end % }}} -% {{{ BuildMITgcm: +% {{{ BuildMITgcm: if perform(org,'BuildMITgcm'), - %load data: + %load data: loaddata(org,'Parameters'); @@ -149,5 +150,5 @@ tline = fgetl(fidi); if ~ischar(tline), break, end - %do the change here: + %do the change here: if strcmpi(tline,' & sNx = 20,'), fprintf(fido,'%s%i%s\n',' & sNx = ',round(Nx/nPx),','); @@ -177,15 +178,15 @@ % }}} addpath(recursivepath([pwd '/../MITgcm'])); -% {{{ RunUncoupledMITgcm: +% {{{ RunUncoupledMITgcm: if perform(org,'RunUncoupledMITgcm'), - %load data: + %load data: loaddata(org,'Parameters'); loaddata(org,'Bathymetry'); loaddata(org,'IceSheetGeometry'); - endtime = round(MITgcmDeltaT * ... - floor(time_step*y2s*async_step_MITgcm_multiplier/MITgcmDeltaT)); - - % {{{ prepare MITgcm + 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') @@ -199,6 +200,6 @@ !\cp ../MITgcm/input/* run !\cp ../MITgcm/input/eedata_uncoupled run/eedata - - %load data: + + %load data: loaddata(org,'Parameters'); @@ -206,14 +207,14 @@ 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); @@ -236,5 +237,5 @@ U=obcUvel*ones(Ny,Nz); writebin('run/OBu.bin',U); - + % zero boundary conditions Z=zeros(Ny,Nz); @@ -249,5 +250,5 @@ tline = fgetl(fidi); if ~ischar(tline), break, end - %do the change here: + %do the change here: if strcmpi(tline,' OB_Iwest = 40*1,'), fprintf(fido,'%s%i%s\n',' OB_Iwest = ',Ny,'*1,'); @@ -269,8 +270,8 @@ % }}} - %start looping: + %start looping: for t=start_time:time_step:final_time, disp(['Year: ' num2str(t)]) - % {{{ generate MITgcm parameter file data + % {{{ generate MITgcm parameter file data fidi=fopen('../MITgcm/input/data','r'); fido=fopen('run/data','w'); @@ -280,5 +281,5 @@ tline = fgetl(fidi); if ~ischar(tline), break, end - %do the change here: + %do the change here: if strcmpi(tline,' xgOrigin = 0.0,'), fprintf(fido,'%s%i%s\n',' xgOrigin = ',xgOrigin,','); @@ -370,207 +371,207 @@ %Configure ISSM -% {{{ CreateMesh: +% {{{ CreateMesh: 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(:); - - %project lat,long: - [x,y]=ll2xy(lat,long,-1); - - 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: + + loaddata(org,'Parameters'); + loaddata(org,'Bathymetry'); + loaddata(org,'IceSheetGeometry'); + + %create model: + md=model(); + + %Grab lat,long from MITgcm: + lat=lat(:); + long=long(:); + + %project lat,long: + [x,y]=ll2xy(lat,long,-1); + + 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,'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(:); - - %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.surface=md.geometry.base+md.geometry.thickness; - - %nothing passes icefront: - pos=find(~ice_mask); - 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.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 - - savemodel(org,md); -end -% }}} -% {{{ ParameterizeIce: + + 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(:); + + %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.surface=md.geometry.base+md.geometry.thickness; + + %nothing passes icefront: + pos=find(~ice_mask); + 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.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 + + savemodel(org,md); +end +% }}} +% {{{ ParameterizeIce: if perform(org,'ParameterizeIce'), - - loaddata(org,'Parameters'); - loaddata(org,'CreateMesh'); - loaddata(org,'MeshGeometry'); - - %miscellaneous - md.miscellaneous.name='test4002'; - - %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=30*ones(md.mesh.numberofvertices,1); - pos=find(md.mask.groundedice_levelset<=0); - md.friction.coefficient(pos)=0; - 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-20)*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]; - - %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); - - %deal with water: - pos=find(md.mask.ice_levelset>0); - md.stressbalance.spcvx(pos)=0; - 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.lat==min(md.mesh.lat)); - md.stressbalance.spcvy(pos)=200; - - %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)=md.initialization.vx(pos); - md.stressbalance.spcvy(pos)=md.initialization.vy(pos); - md.stressbalance.spcvz(pos)=md.initialization.vz(pos); - md.masstransport.spcthickness(pos)=md.geometry.thickness(pos); - - 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: + + loaddata(org,'Parameters'); + loaddata(org,'CreateMesh'); + loaddata(org,'MeshGeometry'); + + %miscellaneous + md.miscellaneous.name='test4002'; + + %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=30*ones(md.mesh.numberofvertices,1); + pos=find(md.mask.groundedice_levelset<=0); + md.friction.coefficient(pos)=0; + 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-20)*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]; + + %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); + + %deal with water: + pos=find(md.mask.ice_levelset>0); + md.stressbalance.spcvx(pos)=0; + 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.lat==min(md.mesh.lat)); + md.stressbalance.spcvy(pos)=200; + + %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)=md.initialization.vx(pos); + md.stressbalance.spcvy(pos)=md.initialization.vy(pos); + md.stressbalance.spcvz(pos)=md.initialization.vz(pos); + md.masstransport.spcthickness(pos)=md.geometry.thickness(pos); + + 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=final_time; - md.timestepping.time_step=time_step; - md.transient.isgroundingline=1; - md.transient.isthermal=0; - md.groundingline.migration='SubelementMigration'; - md.groundingline.melt_interpolation='SubelementMelt2'; - md.groundingline.friction_interpolation='SubelementFriction2'; - - md.cluster=generic('name',oshostname(),'np',2); - md=solve(md,'Transient'); - - savemodel(org,md); + + loaddata(org,'Parameters'); + loaddata(org,'ParameterizeIce'); + + %timestepping: + md.timestepping.final_time=final_time; + md.timestepping.time_step=time_step; + md.transient.isgroundingline=1; + md.transient.isthermal=0; + md.groundingline.migration='SubelementMigration'; + md.groundingline.melt_interpolation='SubelementMelt2'; + md.groundingline.friction_interpolation='SubelementFriction2'; + + md.cluster=generic('name',oshostname(),'np',2); + md=solve(md,'Transient'); + + savemodel(org,md); end % }}} %Run MITgcm/ISSM -% {{{ RunCoupledMITgcmISSM: +% {{{ RunCoupledMITgcmISSM: if perform(org,'RunCoupledMITgcmISSM'), - %load data: - loaddata(org,'Parameters'); - loaddata(org,'ParameterizeIce'); - loaddata(org,'Bathymetry'); - loaddata(org,'IceSheetGeometry'); + %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 + % {{{ prepare MITgcm % rename previous run directory and create new one if exist ('run.old') @@ -631,5 +632,5 @@ tline = fgetl(fidi); if ~ischar(tline), break, end - %do the change here: + %do the change here: if strcmpi(tline,' OB_Iwest = 40*1,'), fprintf(fido,'%s%i%s\n',' OB_Iwest = ',Ny,'*1,'); @@ -649,31 +650,31 @@ 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='SubelementMelt2'; - md.groundingline.friction_interpolation='SubelementFriction2'; - - % }}} - - %start looping: - results=md.results; - - for t=start_time:time_step:final_time + + % {{{ 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='SubelementMelt2'; + 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; + %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]); @@ -681,49 +682,49 @@ 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); + % {{{ 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 + %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,' 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,','); @@ -734,12 +735,12 @@ continue; end - fprintf(fido,'%s\n',tline); - end - %close files - fclose(fidi); - fclose(fido); - % }}} - - % {{{ generate initial MITgcm conditions + fprintf(fido,'%s\n',tline); + end + %close files + fclose(fidi); + fclose(fido); + % }}} + + % {{{ generate initial MITgcm conditions ds=round(endtime/MITgcmDeltaT); if t>start_time @@ -760,5 +761,5 @@ 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) @@ -794,5 +795,5 @@ % }}} - % {{{ system call to run MITgcm + % {{{ system call to run MITgcm cd run eval(['!mpiexec -np ' int2str(nPx*nPy) ' ./mitgcmuv']); @@ -810,55 +811,55 @@ % }}} - %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: + %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'); - endtime = round(MITgcmDeltaT * floor(final_time*y2s/MITgcmDeltaT)); - outputtime = round(MITgcmDeltaT * floor(time_step*y2s/MITgcmDeltaT)); - - % {{{ prepare MITgcm + + loaddata(org,'Parameters'); + loaddata(org,'ParameterizeIce'); + loaddata(org,'Bathymetry'); + loaddata(org,'IceSheetGeometry'); + endtime = round(MITgcmDeltaT * floor(final_time*y2s/MITgcmDeltaT)); + outputtime = round(MITgcmDeltaT * floor(time_step*y2s/MITgcmDeltaT)); + + % {{{ prepare MITgcm % rename previous run directory and create new one if exist ('run.old') @@ -918,5 +919,5 @@ tline = fgetl(fidi); if ~ischar(tline), break, end - %do the change here: + %do the change here: if strcmpi(tline,' OB_Iwest = 40*1,'), fprintf(fido,'%s%i%s\n',' OB_Iwest = ',Ny,'*1,'); @@ -937,5 +938,5 @@ writebin('run/icetopo.bin',draft); % }}} - % {{{ generate MITgcm parameter file data + % {{{ generate MITgcm parameter file data fidi=fopen('../MITgcm/input/data','r'); fido=fopen('run/data','w'); @@ -945,5 +946,5 @@ tline = fgetl(fidi); if ~ischar(tline), break, end - %do the change here: + %do the change here: if strcmpi(tline,' xgOrigin = 0.0,'), fprintf(fido,'%s%i%s\n',' xgOrigin = ',xgOrigin,','); @@ -996,22 +997,22 @@ fclose(fido); % }}} - - 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.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); - -% %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']); + 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); + +% %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 % }}} @@ -1027,32 +1028,32 @@ melting_rate_4=readbin(fnm,[Nx Ny]); field_names ={'Base1','Melting1','Vx2','Vy2','Thickness2','Base2','MaskGroundediceLevelset2','FloatingiceMeltingRate2',... - 'Melting2','Vx3','Vy3','Thickness3','Base3','MaskGroundediceLevelset3','FloatingiceMeltingRate3',... - 'Melting3','Vx4','Vy4','Thickness4','Base4','MaskGroundediceLevelset4','FloatingiceMeltingRate4','Melting4'}; + 'Melting2','Vx3','Vy3','Thickness3','Base3','MaskGroundediceLevelset3','FloatingiceMeltingRate3',... + 'Melting3','Vx4','Vy4','Thickness4','Base4','MaskGroundediceLevelset4','FloatingiceMeltingRate4','Melting4'}; field_tolerances={2e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,... - 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13,... - 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13 }; + 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13,... + 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 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).MaskGroundediceLevelset),... - (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).MaskGroundediceLevelset),... - (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).MaskGroundediceLevelset),... - (md.results.TransientSolution(4).BasalforcingsFloatingiceMeltingRate),... - (melting_rate_4(:)),... - }; + (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).MaskGroundediceLevelset),... + (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).MaskGroundediceLevelset),... + (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).MaskGroundediceLevelset),... + (md.results.TransientSolution(4).BasalforcingsFloatingiceMeltingRate),... + (melting_rate_4(:)),... + }; Index: /issm/trunk-jpl/test/NightlyRun/test4003.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test4003.m (revision 24455) +++ /issm/trunk-jpl/test/NightlyRun/test4003.m (revision 24456) @@ -3,14 +3,15 @@ % %Script control parameters -steps=[1 2 3 4 5 6 7 8 9 10 11]; +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 ${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); +!mkdir Models +org=organizer('repository','Models','prefix','IceOcean.','steps',steps); presentdirectory=pwd;