Changeset 26259

Timestamp:

05/12/21 11:41:59 (4 years ago)

Author:

adhikari

Message:

CHG: steps 6 and 7 are updated. some aspects are still unclear, e.g., whether GRACE monthly solutions should be cumsum'ed

File:

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

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

@@ -2 +2 @@
 addpath('../Data','../Functions');
 
-steps=[1:7];
+steps=[6:7];
 
 if any(steps==1) % {{{ 
@@ -142 +142 @@
         sol2 = (md.results.TransientSolution.Sealevel-md.results.TransientSolution.Bed)*1000;   % [mm]
 
-        sol_name={'Change in water equivalent height [cm]', 'Relative sea-level [mm]'};
+        sol_name={'Change in water equivalent height [cm]', 'Relative sea level [mm]'};
         fig_name={'Fig_dH.pdf','Fig_slr.pdf'};
 
@@ -200 +200 @@
         disp('Projecting  loads onto the mesh...');
         time_range = 2007 + [15 350]/365;
-        water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2));
-
-        [ne,nt]=size(water_load);
-        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
-
+        onvertex = 1; % map data on vertex. If 0, it maps on the elemental centroid. 
+        water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2),onvertex);
+        rho_water2ice = md.materials.rho_freshwater/md.materials.rho_ice; 
+        ice_load = water_load*rho_water2ice; % ice height equivalent. 
+
+        % masstransport evalulates diff between the successive times, so we should cumsum. 
+        num_time = size(ice_load,2); 
+        md.masstransport.spcthickness = [md.geometry.thickness + ice_load];
+        md.masstransport.spcthickness(md.mesh.numberofvertices+1,:)=[0:num_time-1]; 
+
+        %Physics 
         md.transient.issmb=0;
-        md.transient.ismasstransport=0;
         md.transient.isstressbalance=0;
         md.transient.isthermal=0;
-        md.transient.isgia=1;   md.transient.isslr=1;
-
-        md.timestepping.start_time=-1;
-        md.timestepping.final_time=nt;
+        md.transient.ismasstransport=1;
+        md.transient.isslc=1;
+
+        md.solidearth.settings.rigid=1;
+        md.solidearth.settings.elastic=1;
+        md.solidearth.settings.rotation=1;
+        
+        %time stepping:
+        md.timestepping.start_time=0;
         md.timestepping.time_step=1;
-        md.timestepping.interp_forcing=0;
-        md.settings.output_frequency=1;
+        md.timestepping.final_time=num_time; 
 
         md.cluster=generic('name',oshostname(),'np',3);
         md.verbose=verbose('111111111');
 
-        md.solidearth.requested_outputs = {'Sealevel','SealevelRSL'};
+        md.solidearth.requested_outputs = {'Sealevel','Bed'};
 
         md=solve(md,'tr');
@@ -232 +239 @@
         md = loadmodel('./Models/SlrGRACE_Transient');
 
-        time = md.solidearth.surfaceload.icethicknesschange(end,:);
+        time = md.masstransport.spcthickness(end,:);
 
         for tt=1:length(time)
-                gmsl(tt) = md.results.TransientSolution(tt).Bslr*1000; % GMSL rate mm/yr
-                sol1(:,tt) = md.solidearth.surfaceload.icethicknesschange(1:end-1,tt)*100;      % ice equivalent height [cm/yr]
-                sol2(:,tt) = md.results.TransientSolution(tt+1).SealevelRSL*1000;       % mm/yr
+                gmsl(tt) = md.results.TransientSolution(tt).Bslc*1000; % [mm]
+                sol1(:,tt) = (md.masstransport.spcthickness(1:end-1,tt)-md.geometry.thickness)*100; % [cm]
+                sol2(:,tt) = (md.results.TransientSolution(tt).Sealevel-md.results.TransientSolution(tt).Bed)*1000;     % [mm]
         end
-        sol_name = {'Change in water equivalent height [cm]', 'Relative sea-level [mm]'};
+        sol_name = {'Change in water equivalent height [cm]', 'Relative sea level [mm]'};
         movie_name = {'Movie_dH.mp4','Movie_slr.mp4'};
-
+        
         res = 1.0;
 
@@ -284 +291 @@
                                 geoshow(flipud(coastlat),flipud(coastlon),'DisplayType','polygon','FaceColor','white');
                         else
-                                geoshow(coastlat,coastlon,'DisplayType','polygon','FaceColor','white');
+                                geoshow(coastlat,coastlon,'DisplayType','polygon','FaceColor','none');
                         end
                         plot(coastlon,coastlat,'k'); hold off;