/[MITgcm]/MITgcm_contrib/osse/utils/cyl2cartuv.m
ViewVC logotype

Annotation of /MITgcm_contrib/osse/utils/cyl2cartuv.m

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


Revision 1.3 - (hide annotations) (download)
Thu Jan 13 15:50:36 2005 UTC (20 years, 6 months ago) by afe
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +1 -10 lines
removed commented cruft from m-files, added handy-dandy ready-made test data

1 afe 1.1 function [u,v] = cyl2cartuv(thetav,rhov,xi,yi,varargin)
2     % [u,v]=cyl2cartuv(thetav,rhov,xi,yi);
3     %
4     % Re-grids model output in cylindrical coords to cartesian.
5     % c is a 2-D or 3-D scalar or z-vector field
6     % xi,yi are vectors of the new regular lat-lon grid to interpolate to.
7     % z is the interpolated data with dimensions of size(xi) by size(yi).
8 afe 1.2 % theta=0 is at 12 o'clock.
9 afe 1.1 %
10     % e.g.
11     % >> t=rdmds('Ttave.0000513360');
12     % >> xi=-179:2:180;yi=-89:2:90;
13     % >> ti=cyl2cart(x,y,t,xi,yi);
14     %
15 afe 1.3 % $Header: /u/gcmpack/MITgcm_contrib/osse/utils/cyl2cartuv.m,v 1.2 2004/07/12 21:43:16 afe Exp $
16 afe 1.1
17     if ~isequal(size(thetav),size(rhov))
18     error('Theta and rho vector arrays must be same size');
19     end
20    
21     %work out mappings of polar to cartesian
22     NN=size(thetav);
23     [ntheta nrho nz]=size(thetav);
24     [RHO,THETA,NZ] = meshgrid(1:nrho,-pi+2*pi/ntheta:2*pi/ntheta:pi,1:nz);
25     [x,y] = pol2cart(THETA(:,:,1),RHO(:,:,1));
26     nx=ntheta;ny=nrho;
27    
28     % break out components
29     uv=thetav.*cos(THETA)+rhov.*sin(THETA);
30     vv=thetav.*-sin(THETA)+rhov.*cos(THETA);
31    
32     X=reshape(x,[1 nx*ny]);
33     Y=reshape(y,[1 nx*ny]);
34     del=griddata_preprocess(Y,X,yi,xi',varargin{:});
35    
36     for k=1:nz;
37     UV=reshape(uv(:,:,k),[1 nx*ny]);
38     VV=reshape(vv(:,:,k),[1 nx*ny]);
39     u(:,:,k)=griddata(Y,X,UV,yi,xi',varargin{:});
40     v(:,:,k)=griddata(Y,X,VV,yi,xi',varargin{:});
41     end % k
42    
43     % Split vertical and time dimensions
44     if size(NN,2)>2
45     u=reshape(u,[size(u,1) size(u,2) NN(3:end)]);
46     v=reshape(v,[size(v,1) size(v,2) NN(3:end)]);
47     end

  ViewVC Help
Powered by ViewVC 1.1.22