Index: /issm/trunk-jpl/test/NightlyRun/test2021.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2021.m	(revision 25428)
+++ /issm/trunk-jpl/test/NightlyRun/test2021.m	(revision 25428)
@@ -0,0 +1,87 @@
+%Test Name: SESAWslr 
+% SESAW method of solving GRD slr 
+% reference: Adhikari et al., 2016, GMD, https://doi.org/10.5194/gmd-9-1087-2016 
+
+md=model;
+md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',500); %500 km resolution mesh
+
+% read in love numbers. 
+love_numbers = lovenumbers('maxdeg',10000);
+
+% compute Green's functions 
+disp(['Computing Greens functions...']); 
+[Grigid,Gelastic,Uelastic]=greensfunctions(md.mesh.elements,md.mesh.lat,md.mesh.long,love_numbers); 
+greens.Grigid = Grigid;
+greens.Gelast = Gelastic;
+greens.Uelast = Uelastic;
+clearvars Grigid Gelastic Uelastic
+
+% compute lat,long at elemental centroids. 
+[late,longe] = latelonge(md.mesh.elements,md.mesh.lat,md.mesh.long); 
+
+% load GRACE data. Same as used in test2020.  
+load('../Data/GRACE_JPL_April2002_WEH.mat'); 
+lat = repmat(lat',720,1);
+lon = repmat(lon,1,360);
+F = scatteredInterpolant(lat(:),lon(:),weh(:)); 
+
+% map GRACE data onto elemental centorids. 
+loads_element = F(late,longe);
+loads_element(isnan(loads_element))=0;
+
+% ocean mask mapped onto the elemental centroids. 
+ocean_element = gmtmask(late,longe); 
+
+% Area of individual elements
+area_element=GetAreasSphericalTria(md.mesh.elements,md.mesh.lat,md.mesh.long,md.solidearth.planetradius); 
+
+% Parameters input for SESAWslr solver. 
+para.ocean_element = ocean_element; 
+para.loads_element = loads_element; 
+para.area_element = area_element; 
+para.earth_density = md.materials.earth_density; 
+para.ocean_density = md.materials.rho_water;
+para.loads_density = md.materials.rho_freshwater; % if land loads are ice, use ice density. 
+
+para.rel_tol = 1e-5; 
+
+% solid earth rheology. 
+para.solidearth = 'rigid'; % 'rigid' or 'elastic'; 
+
+% rotational feedbacks. 
+para.rotational.flag = 0; % Rotational flag on (1) or off (0) 
+para.rotational.earth_radius = md.solidearth.planetradius; 
+para.rotational.load_love_k2 = love_numbers.k(3); 
+para.rotational.tide_love_k2 = love_numbers.tk(3); 
+para.rotational.tide_love_h2 = love_numbers.th(3);  
+para.rotational.tide_love_k2secular = love_numbers.tk2secular;  
+para.rotational.moi_p = md.solidearth.rotational.polarmoi; 
+para.rotational.moi_e = md.solidearth.rotational.equatorialmoi; 
+para.rotational.omega = md.solidearth.rotational.angularvelocity;
+
+% solve: Rigid without rotational feedbacks. 
+disp(['Solving sesaw-slr for Rigid Earth WITHOUT rotational feedback...']); 
+[eus_rigid,rsl_rigid] = SESAWslr(md.mesh.elements,md.mesh.lat,md.mesh.long,greens,para); 
+
+% solve: Rigid with rotational feedbacks. 
+para.rotational.flag = 1; 
+disp(['Solving sesaw-slr for Rigid Earth WITH rotational feedback...']); 
+[eus_rigid_rot,rsl_rigid_rot] = SESAWslr(md.mesh.elements,md.mesh.lat,md.mesh.long,greens,para); 
+
+% solve: Elastic with rotational feedbacks. 
+para.solidearth = 'elastic'; 
+disp(['Solving sesaw-slr for Elastic Earth WITH rotational feedback...']); 
+[eus_elast_rot,rsl_elast_rot] = SESAWslr(md.mesh.elements,md.mesh.lat,md.mesh.long,greens,para); 
+
+% solve: Elastic with rotational feedbacks. 
+para.rotational.flag = 0; 
+disp(['Solving sesaw-slr for Elastic Earth WITHOUT rotational feedback...']); 
+[eus_elast,rsl_elast] = SESAWslr(md.mesh.elements,md.mesh.lat,md.mesh.long,greens,para); 
+
+%Fields and tolerances to track changes
+field_names={'eus_rigid','eus_rigid_rot','eus_elast','eus_elast_rot',...
+	          'rsl_rigid','rsl_rigid_rot','rsl_elast','rsl_elast_rot'};
+field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13};
+field_values={eus_rigid,eus_rigid_rot,eus_elast,eus_elast_rot,...
+	          rsl_rigid,rsl_rigid_rot,rsl_elast,rsl_elast_rot};
+