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

Contents 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 - (show 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 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 loc.lon(loc.lon<0)=loc.lon(loc.lon<0)+360;
20
21 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 if any(data.Grid_0001(:)>1e10), continue, end
37 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