/[MITgcm]/MITgcm_contrib/gael/profilesMatlabProcessing/profiles_misc/aviso_get_madt.m
ViewVC logotype

Annotation of /MITgcm_contrib/gael/profilesMatlabProcessing/profiles_misc/aviso_get_madt.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.2 - (hide annotations) (download)
Wed Apr 20 20:37:22 2011 UTC (14 years, 3 months ago) by roquet
Branch: MAIN
CVS Tags: checkpoint65x, checkpoint65r, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, HEAD
Changes since 1.1: +3 -1 lines
- bug fix for longitudes (must be between 0 and 360 in opendap grid)
- create the variable prof_madt_aviso first, then fill it with extracted values.

1 roquet 1.1 function [madt]=aviso_get_madt(loc,grid_aviso)
2     % function [madt]=aviso_get_madt(loc,grid)
3     % extract madt value (mapped absolute dynamic topography) at 1 location
4     % linear interpolation in space and time at each location.
5     %
6     % loc is a structure with three array of size (N,1) :
7     % loc.lat : latitude [deg_north]
8     % loc.lon : longitude [deg_east]
9     % loc.time : date in Matlab julian date format
10     %
11     % grid is the structure returned by aviso_get_grid
12     %
13     % madt is the output array [unit: m], of size (N,1).
14     %
15    
16     N=length(loc.lat);
17     madt=NaN*zeros(N,1);
18    
19 roquet 1.2 loc.lon(loc.lon<0)=loc.lon(loc.lon<0)+360;
20    
21 roquet 1.1 Itime=interp1(grid_aviso.time,(1:grid_aviso.Ntime)'-1,loc.time);
22     Ilon=interp1(grid_aviso.NbLongitudes,(1:grid_aviso.Nlon)'-1,loc.lon);
23     Ilat=interp1(grid_aviso.NbLatitudes,(1:grid_aviso.Nlat)'-1,loc.lat);
24    
25     [X,Y,T]=meshgrid([0 1],[0 1],[0 1]);
26    
27     for kk=1:N,
28    
29     if isnan(Itime(kk)*Ilon(kk)*Ilat(kk)), continue, end
30     it=floor(Itime(kk));
31     il=floor(Ilon(kk));
32     iL=floor(Ilat(kk));
33     URL=sprintf('http://opendap.aviso.oceanobs.com/thredds/dodsC/dataset-duacs-dt-ref-global-merged-madt-h?Grid_0001[%d:%d][%d:%d][%d:%d]',...
34     it,it+1,il,il+1,iL,iL+1);
35     a=loaddap(URL);data=a.Grid_0001; % data.Grid_0001 : lon x lat x time
36 roquet 1.2 if any(data.Grid_0001(:)>1e10), continue, end
37 roquet 1.1 madt(kk)=interp3(X,Y,T,data.Grid_0001/100,Ilon(kk)-il,Ilat(kk)-iL,Itime(kk)-it);
38    
39     end

  ViewVC Help
Powered by ViewVC 1.1.22