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

  • TabularUnified ../trunk-jpl/examples/EsaGRACE/grace.m

     
    1 % map grace loads in meters [m] of water equivalent height
    2 function water_load = grace(index,lat,lon,tmin,tmax);
    3 
    4         %compute centroids using (lat,lon) data {{{
    5         ne = length(index); % number of elements
    6         % lat -> [0,180]; long -> [0,360] to compute centroids
    7         lat=90-lat;
    8         lon(lon<0)=180+(180+lon(lon<0));
    9        
    10         ax_0=lat(index(:,1)); ay_0=lon(index(:,1));
    11         bx_0=lat(index(:,2)); by_0=lon(index(:,2));
    12         cx_0=lat(index(:,3)); cy_0=lon(index(:,3));
    13         % find whether long is 0 or 360! This is important to compute centroids as well as elemental area
    14         for ii=1:ne
    15                 if (min([ay_0(ii),by_0(ii),cy_0(ii)])==0 && max([ay_0(ii),by_0(ii),cy_0(ii)])>180)
    16                         if ay_0(ii)==0
    17                                 ay_0(ii)=360;
    18                         end
    19                         if by_0(ii)==0
    20                                 by_0(ii)=360;
    21                         end
    22                         if cy_0(ii)==0
    23                                 cy_0(ii)=360;
    24                         end
    25                 end
    26         end
    27         % correction at the north pole
    28         ay_0(ax_0==0)=(by_0(ax_0==0)+cy_0(ax_0==0))./2;
    29         by_0(bx_0==0)=(cy_0(bx_0==0)+ay_0(bx_0==0))./2;
    30         cy_0(cx_0==0)=(ay_0(cx_0==0)+by_0(cx_0==0))./2;
    31         % correction at the south pole
    32         ay_0(ax_0==180)=(by_0(ax_0==180)+cy_0(ax_0==180))./2;
    33         by_0(bx_0==180)=(cy_0(bx_0==180)+ay_0(bx_0==180))./2;
    34         cy_0(cx_0==180)=(ay_0(cx_0==180)+by_0(cx_0==180))./2;
    35         %
    36         lat_element=(ax_0+bx_0+cx_0)/3;
    37         lon_element=(ay_0+by_0+cy_0)/3;
    38 
    39         % [lat,lon] \in [-90:90,-180,180];
    40         lat_element_0 = 90-lat_element;         lon_element_0 = lon_element;
    41         lon_element_0(lon_element>180) = (lon_element(lon_element>180)-180) - 180;
    42         % }}}
    43 
    44         % Monthly GRACE data
    45         filename=['GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc'];
    46         time_0=ncread(filename,'time'); % days since 2002-01-01 00:00:00 UTC
    47         long_0=ncread(filename,'lon'); % longitudes: 0.27-359.75
    48         lati_0=ncread(filename,'lat'); % latitudes: -89.75:89.75
    49         rec=ncread(filename,'lwe_thickness');% rec_ensemble_mean [cm]
    50 
    51         time_yr = 2002+time_0/365; % [yr]
    52 
    53         [nn_0,mm_0] = size(squeeze(rec(:,:,1)));
    54         for jj=1:mm_0     % chose a latitude
    55                 for kk=1:nn_0
    56                         ii=nn_0*(jj-1)+kk;
    57                         lat_0(ii)=lati_0(jj);
    58                         lon_0(ii)=long_0(kk);
    59                         tws_monthly(:,ii) = rec(kk,jj,:);
    60                 end
    61         end
    62         %
    63         %% translate variables as grace-related variables -- so I do not need to do too much editing
    64         grace_monthly=tws_monthly;
    65         grace_monthly(grace_monthly<-1344.6016 | grace_monthly>858.5046)=0;
    66 
    67         % detrend over the entire time period
    68         grace_monthly = detrend(grace_monthly); % 159X64800 => remove trends from each column!
    69 
    70         % fill out the blanks {{{
    71        
    72         lat_grace=lat_0;
    73         lon_grace=lon_0;
    74         num_org=length(lon_grace);
    75 
    76         qq=1;                   mm=1;
    77         for jj=2:num_org-1
    78                 if (lat_grace(jj)~=lat_grace(jj+1))
    79                         lat_new(qq)=lat_grace(jj);
    80                         lon_new(qq)=lon_grace(jj)+(lon_grace(jj)-lon_grace(jj-1));
    81                         load_new(:,qq)=grace_monthly(:,mm);
    82                         lat_new(qq+1)=lat_grace(jj);
    83                         lon_new(qq+1)=lon_grace(mm)-(lon_grace(jj)-lon_grace(jj-1));
    84                         load_new(:,qq+1)=grace_monthly(:,jj);
    85                         qq=qq+2;
    86                         mm=jj+1; % to find out the value for monthly data
    87                 end
    88         end
    89        
    90         num_add=length(lat_new);
    91         num_plus=num_org+num_add;
    92 
    93         lat_grace_plus=zeros(num_plus,1);
    94         lon_grace_plus=zeros(num_plus,1);
    95         load_grace_plus=zeros(length(time_0),num_plus);
    96        
    97         lat_grace_plus(1:num_org)=lat_grace;
    98         lat_grace_plus(1+num_org:num_plus)=lat_new;
    99         lon_grace_plus(1:num_org)=lon_grace;
    100         lon_grace_plus(1+num_org:num_plus)=lon_new;
    101         load_grace_plus(:,1:num_org)=grace_monthly;
    102         load_grace_plus(:,1+num_org:num_plus)=load_new; %(:,:);
    103         % }}}
    104 
    105         % choose selected months ONLY
    106         [diff1,pos1] = min(abs(tmin-time_yr));
    107         [diff2,pos2] = min(abs(tmax-time_yr));
    108 
    109         time_yr=time_yr(pos1:pos2);
    110         load_grace_plus=load_grace_plus(pos1:pos2,:);
    111         num_yr=length(time_yr);
    112         water_load_0=zeros(ne,num_yr);
    113 
    114         for jj=1:num_yr
    115                 water_load_0(:,jj) = griddata(lat_grace_plus,lon_grace_plus,load_grace_plus(jj,:),lat_element_0,lon_element);
    116                 disp([num2str(jj),' of ',num2str(num_yr),' months done!']);
    117         end
    118 
    119         water_load=water_load_0/100;            % cm -> meters of water
    120         water_load(isnan(water_load))=0;
    121 

  • TabularUnified ../trunk-jpl/examples/EsaWahr/wahr.m 

    Property changes on: ../trunk-jpl/examples/EsaGRACE/grace.m
___________________________________________________________________
Deleted: svn:executable
## -1 +0,0 ##
-*
\ No newline at end of property
     
    1 % compute semi-analytic solutions for a disc load
    2 function [vert, horz] = wahr(disc_rad,xi,love_h,love_l);
    3 
    4         disc_rad = disc_rad/1000; % km
    5         % compute P(x), dP(x)/dx, d2P(x)/dx2
    6         %---------------------------------------------------------------------
    7         % compute p_value
    8         theta=km2deg(xi/1000)';
    9         ang = theta/180*pi;
    10         alpha=cos(ang);
    11         m=length(alpha);
    12         n=length(love_h)-1;
    13         p_value = p_polynomial_value(m,n,alpha);
    14         p_prime = p_polynomial_prime(m,n,alpha);
    15         %---------------------------------------------------------------------
    16         nn=[0:n];
    17         nn_plus_1=nn+1;
    18 
    19         % disc radius in degree
    20         disc_rad = km2deg(disc_rad)/180*pi;
    21         tau=zeros(size(love_h));
    22         tau(1) = 0.5*(1-cos(disc_rad)); % tau_0
    23         p_value_disc = p_polynomial_value(1,n+1,cos(disc_rad));
    24         p_prime_disc = p_polynomial_prime(1,n,cos(disc_rad));
    25         for jj=2:n+1
    26                 nnn = jj-1;
    27                 tau(jj) = 0.5 * (p_value_disc(jj-1) - p_value_disc(jj+1));
    28         end
    29 
    30         const=zeros(size(love_h));
    31         for jj=1:n+1
    32                 const(jj) = 1/(2*(jj-1)+1);
    33         end
    34 
    35         disc=sum(bsxfun(@times,p_value,tau'),2);
    36 
    37         g1 = -sum(bsxfun(@times,p_value,(tau.*love_h.*const)'),2);
    38         g5 = -sum(bsxfun(@times,-sin(ang),bsxfun(@times,p_prime,(tau.*love_l.*const)')),2);
    39 
    40         % coeff
    41         coeff = 1000*4*pi*(6.67408*10^-11)*(6.3781*10^6)/9.81;
    42 
    43         % vertical and horizontal solutions in mm
    44         vert = g1*coeff*1000; % mm
    45         horz = g5*coeff*1000; % mm
    46 
Note: See TracBrowser for help on using the repository browser.