Changeset 25937

Timestamp:

01/14/21 16:11:11 (4 years ago)

Author:

adhikari

Message:

CHG: upgraded slr solver parameters etc to comply with the new formulation

File:

• issm/trunk-jpl/examples/SlrGRACE/runme.m (modified) (10 diffs)

issm/trunk-jpl/examples/SlrGRACE/runme.m

@@ -2 +2 @@
 addpath('../Data','../Functions');
 
-steps=[1:7];
-
-try
- % [1:7]
+steps=[6];
 
 if any(steps==1) %{{{
@@ -11 +8 @@
 
         numrefine=1;
-        resolution=150*1e3;                     % inital resolution [m]
+        resolution=300*1e3;                     % inital resolution [m]
         radius = 6.371012*10^6;         % mean radius of Earth, m
         mindistance_coast=150*1e3;      % coastal resolution [m]
@@ -22 +19 @@
         for i=1:numrefine,
 
-                md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long);
+                ocean_mask=gmtmask(md.mesh.lat,md.mesh.long);
 
                 distance=zeros(md.mesh.numberofvertices,1);
 
-                pos=find(~md.mask.ocean_levelset);      coaste.lat=md.mesh.lat(pos);    coaste.long=md.mesh.long(pos); 
-                pos=find(md.mask.ocean_levelset);       coasto.lat=md.mesh.lat(pos);    coasto.long=md.mesh.long(pos); 
+                pos=find(~ocean_mask);  coaste.lat=md.mesh.lat(pos);    coaste.long=md.mesh.long(pos); 
+                pos=find(ocean_mask);   coasto.lat=md.mesh.lat(pos);    coasto.long=md.mesh.long(pos); 
 
                 for j=1:md.mesh.numberofvertices
                         phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi;
-                        if md.mask.ocean_levelset(j),
+                        if ocean_mask(j),
                                 phi2=coaste.lat/180*pi; lambda2=coaste.long/180*pi;
                                 deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1);
@@ -44 +41 @@
                 pos=find(distance<mindistance_coast); distance(pos)=mindistance_coast;
                 
-                pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land);
+                pos2=find(ocean_mask~=1 & distance>mindistance_land);
                 distance(pos2)=mindistance_land;
 
@@ -51 +48 @@
         end
 
-        md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long);
+        ocean_mask=gmtmask(md.mesh.lat,md.mesh.long);
+        pos = find(ocean_mask==0); 
+        md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1); 
+        md.mask.ocean_levelset(pos)=1; 
 
         save ./Models/SlrGRACE_Mesh md;
@@ -67 +67 @@
         water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2));
 
-        md.slr.deltathickness = water_load*md.materials.rho_freshwater/md.materials.rho_ice;
+        md.solidearth.surfaceload.icethicknesschange = water_load*md.materials.rho_freshwater/md.materials.rho_ice;
 
         save ./Models/SlrGRACE_Loads md;
 
-        plotmodel (md,'data',md.slr.deltathickness,'view',[90 -90],'caxis',[-.1 .1],'title','Ice height equivalent [m]');
+        plotmodel (md,'data',md.solidearth.surfaceload.icethicknesschange,'view',[90 -90],'caxis',[-.1 .1],'title','Ice height equivalent [m]');
 end %}}}
 
@@ -80 +80 @@
         md.solidearth.lovenumbers = lovenumbers('maxdeg',10000);
 
-        md.mask.ocean_levelset = ones(md.mesh.numberofvertices,1);
-        md.mask.ice_levelset = ones(md.mesh.numberofvertices,1);
-        pos=find(md.slr.deltathickness~=0);
-        md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
-        md.mask.land_levelset = 1-md.mask.ocean_levelset;
-
-        md.slr.sealevel=zeros(md.mesh.numberofvertices,1);
-        md.slr.steric_rate=zeros(md.mesh.numberofvertices,1);
-        md.slr.ocean_area_scaling = 1;
-
-        di=md.materials.rho_ice/md.materials.rho_water;
-        md.geometry.thickness=ones(md.mesh.numberofvertices,1);
-        md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1);
-        md.geometry.base=md.geometry.surface-md.geometry.thickness;
-        md.geometry.bed=md.geometry.base;
+        md.mask.ice_levelset=-md.mask.ocean_levelset; 
+
+        md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
+        md.dsl.global_average_thermosteric_sea_level_change=[0;0];
+        md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);
+        md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
+
+        md.solidearth.settings.ocean_area_scaling=1;
+
+        % arbitary to pass consistency check. 
+        md.geometry.bed=-ones(md.mesh.numberofvertices,1);
+        md.geometry.surface=ones(md.mesh.numberofvertices,1);
+        md.geometry.base=md.geometry.bed; 
+        md.geometry.thickness=md.geometry.surface-md.geometry.base; 
 
         md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
@@ -109 +108 @@
         md = loadmodel('./Models/SlrGRACE_Parameterization');
 
-        md.slr.rigid=1;
-        md.slr.elastic=1;
-        md.slr.rotation=1;
+        md.solidearth.settings.rigid=1;
+        md.solidearth.settings.elastic=1;
+        md.solidearth.settings.rotation=1;
+
+        md.solidearth.requested_outputs = {'Sealevel','SealevelRSL'};
 
         md.cluster=generic('name',oshostname(),'np',3);
@@ -125 +126 @@
         md = loadmodel('./Models/SlrGRACE_Solution');
 
-        sol1 = md.slr.deltathickness*100;                                               % [cm]
-        sol2 = md.results.SealevelriseSolution.Sealevel*1000;   % [mm]
+        sol1 = md.solidearth.surfaceload.icethicknesschange*100;                                                % [cm]
+        sol2 = md.results.SealevelriseSolution.SealevelRSL*1000;        % [mm]
 
         sol_name={'Change in water equivalent height [cm]', 'Relative sea-level [mm]'};
@@ -190 +191 @@
 
         [ne,nt]=size(water_load);
-        md.slr.deltathickness = zeros(ne+1,nt);
-        md.slr.deltathickness(1:ne,:) = water_load*md.materials.rho_freshwater/md.materials.rho_ice;
-        md.slr.deltathickness(ne+1,:)=[1:nt]; % times
+        md.solidearth.surfaceload.icethicknesschange = zeros(ne+1,nt);
+        md.solidearth.surfaceload.icethicknesschange(1:ne,:) = water_load*md.materials.rho_freshwater/md.materials.rho_ice;
+        md.solidearth.surfaceload.icethicknesschange(ne+1,:)=[1:nt]; % times
 
         md.transient.issmb=0;