source: issm/oecreview/Archive/22755-22818/ISSM-22809-22810.diff@ 22819

Last change on this file since 22819 was 22819, checked in by Mathieu Morlighem, 7 years ago

CHG: added Archive/22755-22818
File size: 6.2 KB
RevLine 
[22819]1Index: ../trunk-jpl/examples/EsaGRACE/grace.m
2===================================================================
3--- ../trunk-jpl/examples/EsaGRACE/grace.m (revision 22809)
4+++ ../trunk-jpl/examples/EsaGRACE/grace.m (nonexistent)
5@@ -1,121 +0,0 @@
6-% map grace loads in meters [m] of water equivalent height
7-function water_load = grace(index,lat,lon,tmin,tmax);
8-
9- %compute centroids using (lat,lon) data {{{
10- ne = length(index); % number of elements
11- % lat -> [0,180]; long -> [0,360] to compute centroids
12- lat=90-lat;
13- lon(lon<0)=180+(180+lon(lon<0));
14-
15- ax_0=lat(index(:,1)); ay_0=lon(index(:,1));
16- bx_0=lat(index(:,2)); by_0=lon(index(:,2));
17- cx_0=lat(index(:,3)); cy_0=lon(index(:,3));
18- % find whether long is 0 or 360! This is important to compute centroids as well as elemental area
19- for ii=1:ne
20- if (min([ay_0(ii),by_0(ii),cy_0(ii)])==0 && max([ay_0(ii),by_0(ii),cy_0(ii)])>180)
21- if ay_0(ii)==0
22- ay_0(ii)=360;
23- end
24- if by_0(ii)==0
25- by_0(ii)=360;
26- end
27- if cy_0(ii)==0
28- cy_0(ii)=360;
29- end
30- end
31- end
32- % correction at the north pole
33- ay_0(ax_0==0)=(by_0(ax_0==0)+cy_0(ax_0==0))./2;
34- by_0(bx_0==0)=(cy_0(bx_0==0)+ay_0(bx_0==0))./2;
35- cy_0(cx_0==0)=(ay_0(cx_0==0)+by_0(cx_0==0))./2;
36- % correction at the south pole
37- ay_0(ax_0==180)=(by_0(ax_0==180)+cy_0(ax_0==180))./2;
38- by_0(bx_0==180)=(cy_0(bx_0==180)+ay_0(bx_0==180))./2;
39- cy_0(cx_0==180)=(ay_0(cx_0==180)+by_0(cx_0==180))./2;
40- %
41- lat_element=(ax_0+bx_0+cx_0)/3;
42- lon_element=(ay_0+by_0+cy_0)/3;
43-
44- % [lat,lon] \in [-90:90,-180,180];
45- lat_element_0 = 90-lat_element; lon_element_0 = lon_element;
46- lon_element_0(lon_element>180) = (lon_element(lon_element>180)-180) - 180;
47- % }}}
48-
49- % Monthly GRACE data
50- filename=['GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc'];
51- time_0=ncread(filename,'time'); % days since 2002-01-01 00:00:00 UTC
52- long_0=ncread(filename,'lon'); % longitudes: 0.27-359.75
53- lati_0=ncread(filename,'lat'); % latitudes: -89.75:89.75
54- rec=ncread(filename,'lwe_thickness');% rec_ensemble_mean [cm]
55-
56- time_yr = 2002+time_0/365; % [yr]
57-
58- [nn_0,mm_0] = size(squeeze(rec(:,:,1)));
59- for jj=1:mm_0 % chose a latitude
60- for kk=1:nn_0
61- ii=nn_0*(jj-1)+kk;
62- lat_0(ii)=lati_0(jj);
63- lon_0(ii)=long_0(kk);
64- tws_monthly(:,ii) = rec(kk,jj,:);
65- end
66- end
67- %
68- %% translate variables as grace-related variables -- so I do not need to do too much editing
69- grace_monthly=tws_monthly;
70- grace_monthly(grace_monthly<-1344.6016 | grace_monthly>858.5046)=0;
71-
72- % detrend over the entire time period
73- grace_monthly = detrend(grace_monthly); % 159X64800 => remove trends from each column!
74-
75- % fill out the blanks {{{
76-
77- lat_grace=lat_0;
78- lon_grace=lon_0;
79- num_org=length(lon_grace);
80-
81- qq=1; mm=1;
82- for jj=2:num_org-1
83- if (lat_grace(jj)~=lat_grace(jj+1))
84- lat_new(qq)=lat_grace(jj);
85- lon_new(qq)=lon_grace(jj)+(lon_grace(jj)-lon_grace(jj-1));
86- load_new(:,qq)=grace_monthly(:,mm);
87- lat_new(qq+1)=lat_grace(jj);
88- lon_new(qq+1)=lon_grace(mm)-(lon_grace(jj)-lon_grace(jj-1));
89- load_new(:,qq+1)=grace_monthly(:,jj);
90- qq=qq+2;
91- mm=jj+1; % to find out the value for monthly data
92- end
93- end
94-
95- num_add=length(lat_new);
96- num_plus=num_org+num_add;
97-
98- lat_grace_plus=zeros(num_plus,1);
99- lon_grace_plus=zeros(num_plus,1);
100- load_grace_plus=zeros(length(time_0),num_plus);
101-
102- lat_grace_plus(1:num_org)=lat_grace;
103- lat_grace_plus(1+num_org:num_plus)=lat_new;
104- lon_grace_plus(1:num_org)=lon_grace;
105- lon_grace_plus(1+num_org:num_plus)=lon_new;
106- load_grace_plus(:,1:num_org)=grace_monthly;
107- load_grace_plus(:,1+num_org:num_plus)=load_new; %(:,:);
108- % }}}
109-
110- % choose selected months ONLY
111- [diff1,pos1] = min(abs(tmin-time_yr));
112- [diff2,pos2] = min(abs(tmax-time_yr));
113-
114- time_yr=time_yr(pos1:pos2);
115- load_grace_plus=load_grace_plus(pos1:pos2,:);
116- num_yr=length(time_yr);
117- water_load_0=zeros(ne,num_yr);
118-
119- for jj=1:num_yr
120- water_load_0(:,jj) = griddata(lat_grace_plus,lon_grace_plus,load_grace_plus(jj,:),lat_element_0,lon_element);
121- disp([num2str(jj),' of ',num2str(num_yr),' months done!']);
122- end
123-
124- water_load=water_load_0/100; % cm -> meters of water
125- water_load(isnan(water_load))=0;
126-
127
128Property changes on: ../trunk-jpl/examples/EsaGRACE/grace.m
129___________________________________________________________________
130Deleted: svn:executable
131## -1 +0,0 ##
132-*
133\ No newline at end of property
134Index: ../trunk-jpl/examples/EsaWahr/wahr.m
135===================================================================
136--- ../trunk-jpl/examples/EsaWahr/wahr.m (revision 22809)
137+++ ../trunk-jpl/examples/EsaWahr/wahr.m (nonexistent)
138@@ -1,46 +0,0 @@
139-% compute semi-analytic solutions for a disc load
140-function [vert, horz] = wahr(disc_rad,xi,love_h,love_l);
141-
142- disc_rad = disc_rad/1000; % km
143- % compute P(x), dP(x)/dx, d2P(x)/dx2
144- %---------------------------------------------------------------------
145- % compute p_value
146- theta=km2deg(xi/1000)';
147- ang = theta/180*pi;
148- alpha=cos(ang);
149- m=length(alpha);
150- n=length(love_h)-1;
151- p_value = p_polynomial_value(m,n,alpha);
152- p_prime = p_polynomial_prime(m,n,alpha);
153- %---------------------------------------------------------------------
154- nn=[0:n];
155- nn_plus_1=nn+1;
156-
157- % disc radius in degree
158- disc_rad = km2deg(disc_rad)/180*pi;
159- tau=zeros(size(love_h));
160- tau(1) = 0.5*(1-cos(disc_rad)); % tau_0
161- p_value_disc = p_polynomial_value(1,n+1,cos(disc_rad));
162- p_prime_disc = p_polynomial_prime(1,n,cos(disc_rad));
163- for jj=2:n+1
164- nnn = jj-1;
165- tau(jj) = 0.5 * (p_value_disc(jj-1) - p_value_disc(jj+1));
166- end
167-
168- const=zeros(size(love_h));
169- for jj=1:n+1
170- const(jj) = 1/(2*(jj-1)+1);
171- end
172-
173- disc=sum(bsxfun(@times,p_value,tau'),2);
174-
175- g1 = -sum(bsxfun(@times,p_value,(tau.*love_h.*const)'),2);
176- g5 = -sum(bsxfun(@times,-sin(ang),bsxfun(@times,p_prime,(tau.*love_l.*const)')),2);
177-
178- % coeff
179- coeff = 1000*4*pi*(6.67408*10^-11)*(6.3781*10^6)/9.81;
180-
181- % vertical and horizontal solutions in mm
182- vert = g1*coeff*1000; % mm
183- horz = g5*coeff*1000; % mm
184-
185
186Property changes on: ../trunk-jpl/examples/EsaWahr/wahr.m
187___________________________________________________________________
188Deleted: svn:executable
189## -1 +0,0 ##
190-*
191\ No newline at end of property
Note: See TracBrowser for help on using the repository browser.