Index: /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/SESAWslr.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/SESAWslr.m (revision 25426) +++ /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/SESAWslr.m (revision 25427) @@ -1,3 +1,3 @@ -function [eust,rsl] = SESAWslr(index,greens,para) +function [eust,rsl] = SESAWslr(index,lat,long,greens,para) %SESAWslr :: computes GRD slr due to applied surface loads based on SESAW method. % @@ -5,5 +5,5 @@ % % Usage: -% [eus,rsl]=greensfunctions(md.mesh.elements,greens,para) +% [eus,rsl]=greensfunctions(md.mesh.elements,md.mesh.lat,md.mesh.long,greens,para) % % greens.Grigid = Gravitational Green's function for the rigid planet. @@ -18,6 +18,16 @@ % para.loads_density = land loads (ice or freshwater) density [kg/m^3] % para.rel_tol = relative tolerance for iterative SLR solver +% % para.solidearth = rheological model for solid Earth: 'rigid' or 'elastic' -% para.isgrd_R = Do we want to account for Rotational feedback? 0 or 1. +% +% para.rotational.flag = Do we want to accoutn for rotational feedback? 1 (yes) or 0 (no) +% para.rotational.earth_radius = mean radius of the Earth [m] +% para.rotational.load_love_k2 = degree 2 load love number k +% para.rotational.tide_love_k2 = degree 2 tide love number k +% para.rotational.tide_love_h2 = degree 2 tide love number h +% para.rotational.tide_love_k2secular = degree 2 secular tide love number k +% para.rotational.moi_p = polar moment of inertia [kg m^2] +% para.rotational.moi_e = mean equatorial moment of inertia [kg m^2] +% para.rotational.omega = mean rotational velocity of earth [rad per second] area_element = para.area_element; @@ -29,37 +39,59 @@ rel_tol = para.rel_tol; +rot_flag = para.rotational.flag; + % total Green's function for elastic earth, i.e. (1+k_l-h_l) if strcmpi(para.solidearth,'rigid') Galpha = greens.Grigid; elseif strcmpi(para.solidearth,'elastic') - Galpha = greens.Grigid + greens.Gelast- greens.Uelast; + Galpha = greens.Grigid + greens.Gelast - greens.Uelast; else error(['Unknown solidearth model:' para.solidearth, '; Should be rigid or elastic']); end +if rot_flag + % set up parameters for rotationalfeedback function. + para_rot = para.rotational; + para_rot.area_element = para.area_element; + para_rot.loads_element = para.loads_element + para_rot.loads_density = para.loads_density; +end + +% densitity ratios density_o_e = rho_o/rho_e; density_l_e = rho_l/rho_e; density_l_o = rho_l/rho_o; +% ocean and earth's surface areas: +ocean_area = sum(ocean_element.*area_element); +earth_area = sum(area_element); + % eustatic term -eust=-density_l_o*sum(loads_element.*area_element)/sum(ocean_element.*area_element); +eust = -density_l_o*sum(loads_element.*area_element)/ocean_area; -term1=3*density_l_e.*sum(bsxfun(@times,Galpha,(loads_element.*area_element)'),2)./sum(area_element); -func3=mean(term1(index),2).*ocean_element; -term3=sum(func3.*area_element)./sum(ocean_element.*area_element); +term1 = 3*density_l_e.*sum(bsxfun(@times,Galpha,(loads_element.*area_element)'),2)./earth_area; +if rot_flag + term1 = term1 + rotationalfeedback(index,lat,long,para_rot); +end +func3 = mean(term1(index),2).*ocean_element; +term3 = sum(func3.*area_element)./ocean_area; % computation of sea level change -rsl=eust+term1-term3; +rsl = eust+term1-term3; -norm_diff = 10; p = 0; - +norm_diff = 10; p = 0; while norm_diff > rel_tol norm_old = sqrt(sum(rsl.^2)); + % + term2 = 3*density_o_e.*sum(bsxfun(@times,Galpha,(mean(rsl(index),2).*ocean_element.*area_element)'),2)./earth_area; + if rot_flag + para_rot.loads_element = mean(rsl(index),2).*ocean_element; + para_rot.loads_density = para.ocean_density; + term2 = term2 + rotationalfeedback(index,lat,long,para_rot); + end + func4 = mean(term2(index),2).*ocean_element; + term4 = sum(func4.*area_element)./ocean_area; % - term2=3*density_o_e.*sum(bsxfun(@times,Galpha,(mean(rsl(index),2).*ocean_element.*area_element)'),2)./sum(area_element); - func4=mean(term2(index),2).*ocean_element; - term4=sum(func4.*area_element)./sum(ocean_element.*area_element); - % - rsl=eust+term1+term2-term3-term4; + rsl = eust+term1+term2-term3-term4; norm_new = sqrt(sum(rsl.^2)); norm_diff = abs(norm_new-norm_old)./norm_old; Index: /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/greensfunctions.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/greensfunctions.m (revision 25426) +++ /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/greensfunctions.m (revision 25427) @@ -1,13 +1,13 @@ -function [Grigid,Gelastic,Uelastic] = greensfunctions(md_mesh,love_numbers) +function [Grigid,Gelastic,Uelastic] = greensfunctions(mesh_elements,mesh_lat,mesh_long,love_numbers) %greensfunctions :: computes Greens funtions for rigid and elastic Earth. % % Usage: -% [Grigid,Gelastic,Uelastic]=greensfunctions(md_mesh,love_numbers) +% [Grigid,Gelastic,Uelastic]=greensfunctions(md.mesh.elements,md.mesh.lat,md.mesh.long,md.solidearth.lovenumbers); -nv=md_mesh.numberofvertices; -ne=md_mesh.numberofelements; +ne = length(mesh_elements); +nv = length(mesh_lat); % compute lat,long at elemental centroids. -[late,longe] = latelonge(md_mesh); +[late,longe] = latelonge(mesh_elements,mesh_lat,mesh_long); % love numbers. @@ -36,5 +36,5 @@ qq=1; for j=1:nv - phi1(:,1)=md_mesh.lat(j)./180.*pi; lambda1(:,1)=md_mesh.long(j)./180.*pi; % size #vertices + phi1(:,1)=mesh_lat(j)./180.*pi; lambda1(:,1)=mesh_long(j)./180.*pi; % size #vertices delPhi=abs(phi2-phi1); delLambda=abs(lambda2-lambda1); alpha=2.*asin(sqrt(sin(delPhi./2).^2+cos(phi1).*cos(phi2).*sin(delLambda./2).^2)); @@ -43,5 +43,5 @@ Uelastic(j,:)=interp1(xx,elast_loveH,cos(alpha)); if (j==500*qq) - display([num2str(j),' of ', num2str(md_mesh.numberofvertices),' vertices done!']); + display([num2str(j),' of ', num2str(nv),' vertices done!']); qq=qq+1; end Index: /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/latelonge.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/latelonge.m (revision 25426) +++ /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/latelonge.m (revision 25427) @@ -1,22 +1,20 @@ -function [late,longe] = latelonge(md_mesh) +function [late,longe] = latelonge(mesh_elements,mesh_lat,mesh_long) %latelonge :: computes lat,long at elemental centroids on spherical surface % % Usage: -% [late,longe]=latelonge(md.mesh) +% [late,longe]=latelonge(md.mesh.elements,md.mesh.lat,md.mesh.long) -nv = md_mesh.numberofvertices; -if length(md_mesh.lat)~=nv | length(md_mesh.long)~=nv - error('lat,long not defined properly at vertices.'); -end +ne = length(mesh_elements); +nv = length(mesh_lat); % lat -> [0,180]; long -> [0,360] to compute centroids -lat=90-md_mesh.lat; lon=md_mesh.long; +lat=90-mesh_lat; lon=mesh_long; lon(lon<0)=180+(180+lon(lon<0)); -ax_0=lat(md_mesh.elements(:,1)); ay_0=lon(md_mesh.elements(:,1)); -bx_0=lat(md_mesh.elements(:,2)); by_0=lon(md_mesh.elements(:,2)); -cx_0=lat(md_mesh.elements(:,3)); cy_0=lon(md_mesh.elements(:,3)); +ax_0=lat(mesh_elements(:,1)); ay_0=lon(mesh_elements(:,1)); +bx_0=lat(mesh_elements(:,2)); by_0=lon(mesh_elements(:,2)); +cx_0=lat(mesh_elements(:,3)); cy_0=lon(mesh_elements(:,3)); % find whether long is 0 or 360! This is important to compute centroids as well as elemental area -for ii=1:md_mesh.numberofelements +for ii=1:ne if (min([ay_0(ii),by_0(ii),cy_0(ii)])==0 && max([ay_0(ii),by_0(ii),cy_0(ii)])>180) if ay_0(ii)==0 Index: /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/rotationalfeedback.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/rotationalfeedback.m (revision 25427) +++ /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/rotationalfeedback.m (revision 25427) @@ -0,0 +1,53 @@ +function [rsl_rot] = rotationalfeedback(index,lat,long,para_rot) +%rotationalfeedback :: computes rotational feedback for slr based on SESAW method. +% +% Usage: +% rsl = otationalfeedback(md.mesh.elements,md.mesh.lat,md.mesh.long,para_rot) +% +% para_rot.area_element = para.area_element [m^2]; +% para_rot.loads_element = surface loads (ice, terrestrial water storae, sea level) [m] +% para_rot.loads_density = surface loads (ice, freshwater, ocean water) density [kg.m^3]; +% +% para_rot.earth_radius = mean radius of the Earth [m] +% para_rot.load_love_k2 = degree 2 load love number k +% para_rot.tide_love_k2 = degree 2 tide love number k +% para_rot.tide_love_h2 = degree 2 tide love number h +% para_rot.tide_love_k2secular = degree 2 secular tide love number k +% para_rot.moi_p = polar moment of inertia [kg m^2] +% para_rot.moi_e = mean equatorial moment of inertia [kg m^2] +% para_rot.omega = mean rotational velocity of earth [rad per second] + +area_element = para_rot.area_element; +loads = para_rot.loads_element; +rho_loads = para_rot.loads_density; +re = para_rot.earth_radius; +load_love_k2 = para_rot.load_love_k2; +tide_love_k2 = para_rot.tide_love_k2; +tide_love_h2 = para_rot.tide_love_h2; +tide_love_k2secular =para_rot.tide_love_k2secular; +moi_p = para_rot.moi_p; +moi_e = para_rot.moi_e; +omega = para_rot.omega; + +% compute lat,long at elemental centroids. +[late,longe] = latelonge(index,lat,long); +lat = lat*pi/180; long = long*pi/180; +late= late*pi/180; longe = longe*pi/180; + +% Perturbation terms for moment of inertia (moi_list): +% computed analytically (see Wu & Peltier, eqs 10 & 32); also consistent with my GMD formulation! +moi_list_1 = -re^2 *rho_loads*sum((loads.*area_element).*(sin(late).*cos(late).*cos(longe))); +moi_list_2 = -re^2 *rho_loads*sum((loads.*area_element).*(sin(late).*cos(late).*sin(longe))); +moi_list_3 = re^2 *rho_loads*sum((loads.*area_element).*(1-sin(late).^2)); + +% compute perturbation terms for angular velocity vector: +m1 = 1/(1-tide_love_k2/tide_love_k2secular) * (1+load_love_k2)/(moi_p-moi_e) * moi_list_1; +m2 = 1/(1-tide_love_k2/tide_love_k2secular) * (1+load_love_k2)/(moi_p-moi_e) * moi_list_2; +m3 = -(1+load_love_k2)/moi_p * moi_list_3; % term associated with fluid number (3-order-of-magnitude smaller) is negelected + +% Green's function (1+k_2-h_2/g): checked against Glenn Milne's thesis Chapter 3 (eqs: 3.3-4, 3.10-11) +% Perturbation terms for angular velocity vector (m1, m2, m3): checked against Mitrovica (2005 Appendix) & Jensen et al (2013 Appendix A3) +% Sea level rotational feedback: checked against GMD eqs 8-9 (only first order terms, i.e., degree 2 order 0 & 1 considered) +% only first order terms are considered now: +rsl_rot = ((1.0+tide_love_k2-tide_love_h2)/9.81)*(omega*re)^2*(-m3/6.0 + 0.5*m3*cos(2*lat) - 0.5*sin(2*lat).*(m1*cos(long)+m2*sin(long))); +