Commit 922c3a10 authored by Pierre Cazenave's avatar Pierre Cazenave

Add the TMD v2.04 code.

parent 496d3f4d
% Bilinear interpolation (neigbouring NaNs avoided)
% in point xt,yt
% x(n),y(m) - coordinates for h(n,m)
% xt(N),yt(N) - coordinates for interpolated values
% Global case is considered ONLY for EXACTLY global solutions,
% i.e. given in lon limits, satisfying:
% ph_lim(2)-ph_lim(1)==360
% (thus for C-grid: x(end)-x(1)==360-dx )
%
% usage:
% [hi]=BLinterp(x,y,h,xt,yt,km);
% km=1 if model is on Cartesian grid (km) otherwise 0
%
function [hi]=BLinterp(x,y,h,xt,yt,km);
%
if nargin<5,km=0;end
%
dx=x(2)-x(1);dy=y(2)-y(1);
glob=0;
if km==0 & x(end)-x(1)==360-dx,glob=1;end
inan=find(isnan(h)>0);
h(inan)=0;
mz=(h~=0);
n=length(x);m=length(y);
[n1,m1]=size(h);
if n~=n1 | m~=m1,
fprintf('Check dimensions\n');
hi=NaN;
return
end
% extend solution both sides for global
if glob==1,
h0=h;mz0=mz;
[k1,k2]=size(x);if k1==1,x=x';end
[k1,k2]=size(y);if k1==1,y=y';end % just for consistency
x=[x(1)-2*dx;x(1)-dx;x;x(end)+dx;x(end)+2*dx];
h=[h(end-1,:);h(end,:);h;h(1,:);h(2,:)];
mz=[mz(end-1,:);mz(end,:);mz;mz(1,:);mz(2,:)];
end
% Adjust lon convention
% THIS should be done only if km==0 !!!!!!!
xti=xt;
if km==0,
ik=find(xti<x(1) | xti>x(end));
if x(end)>180,xti(ik)=xti(ik)+360;end
if x(end)<0, xti(ik)=xti(ik)-360;end
ik=find(xti>360);
xti(ik)=xti(ik)-360;
ik=find(xti<-180);
xti(ik)=xti(ik)+360;
end
%
[X,Y]=meshgrid(x,y);
q=1/(4+2*sqrt(2));q1=q/sqrt(2);
h1=q1*h(1:end-2,1:end-2)+q*h(1:end-2,2:end-1)+q1*h(1:end-2,3:end)+...
q1*h(3:end,1:end-2)+q*h(3:end,2:end-1)+q1*h(3:end,3:end)+...
q*h(2:end-1,1:end-2)+q*h(2:end-1,3:end);
mz1=q1*mz(1:end-2,1:end-2)+q*mz(1:end-2,2:end-1)+q1*mz(1:end-2,3:end)+...
q1*mz(3:end,1:end-2)+q*mz(3:end,2:end-1)+q1*mz(3:end,3:end)+...
q*mz(2:end-1,1:end-2)+q*mz(2:end-1,3:end);
mz1(find(mz1==0))=1;
h2=h;
h2(2:end-1,2:end-1)=h1./mz1;
ik=find(mz==1);
h2(ik)=h(ik);
h2(find(h2==0))=NaN;
hi=interp2(X,Y,h2',xti,yt);hi=conj(hi);
if glob==1, h=h0;mz=mz0;end
return
% Bilinear interpolation (neigbouring NaNs avoided)
% in point xt,yt
% x(n),y(m) - coordinates for h(n,m)
% xt(N),yt(N) - coordinates for interpolated values
% Global case is considered ONLY for EXACTLY global solutions,
% i.e. given in lon limits, satisfying:
% ph_lim(2)-ph_lim(1)==360
% (thus for C-grid: x(end)-x(1)==360-dx )
%
function [hi]=BLinterp(x,y,h,xt,yt);
%
dx=x(2)-x(1);dy=y(2)-y(1);
glob=0;
if x(end)-x(1)==360-dx,glob=1;end
inan=find(isnan(h)>0);
h(inan)=0;
mz=(h~=0);
n=length(x);m=length(y);
[n1,m1]=size(h);
if n~=n1 | m~=m1,
fprintf('Check dimensions\n');
hi=NaN;
return
end
% extend solution both sides for global
if glob==1,
h0=h;mz0=mz;
[k1,k2]=size(x);if k1==1,x=x';end
[k1,k2]=size(y);if k1==1,y=y';end % just for consistency
x=[x(1)-2*dx;x(1)-dx;x;x(end)+dx;x(end)+2*dx];
h=[h(end-1,:);h(end,:);h;h(1,:);h(2,:)];
mz=[mz(end-1,:);mz(end,:);mz;mz(1,:);mz(2,:)];
end
% Adjust lon convention
xti=xt;
ik=find(xti<x(1) | xti>x(end));
if x(end)>180,xti(ik)=xti(ik)+360;end
if x(end)<0, xti(ik)=xti(ik)-360;end
ik=find(xti>360);
xti(ik)=xti(ik)-360;
ik=find(xti<-180);
xti(ik)=xti(ik)+360;
%
[X,Y]=meshgrid(x,y);
q=1/(4+2*sqrt(2));q1=q/sqrt(2);
h1=q1*h(1:end-2,1:end-2)+q*h(1:end-2,2:end-1)+q1*h(1:end-2,3:end)+...
q1*h(3:end,1:end-2)+q*h(3:end,2:end-1)+q1*h(3:end,3:end)+...
q*h(2:end-1,1:end-2)+q*h(2:end-1,3:end);
mz1=q1*mz(1:end-2,1:end-2)+q*mz(1:end-2,2:end-1)+q1*mz(1:end-2,3:end)+...
q1*mz(3:end,1:end-2)+q*mz(3:end,2:end-1)+q1*mz(3:end,3:end)+...
q*mz(2:end-1,1:end-2)+q*mz(2:end-1,3:end);
mz1(find(mz1==0))=1;
h2=h;
h2(2:end-1,2:end-1)=h1./mz1;
ik=find(mz==1);
h2(ik)=h(ik);
h2(find(h2==0))=NaN;
hi=interp2(X,Y,h2',xti,yt);hi=conj(hi);
if glob==1, h=h0;mz=mz0;end
return
function [hu,hv] = Huv(hz);
% Usage [hu,hv] = Huv(hz);
[n,m] = size(hz);
[mu,mv,mz] = Muv(hz);
indxm = [n,1:n-1];
indym = [m,1:m-1];
hu = mu.*(hz + hz(indxm,:))/2;
hv = mv.*(hz + hz(:,indym))/2;
% Lana Erofeeva, re-make for matlab OCT 2004
% usage:
% [dh]=InferMinor(zmaj,cid,SDtime);
%
% Based on Richard Ray's code perth2
% Return correction for the 16 minor constituents
% zeros, if not enough input constituents for inference
% Input:
% cid(ncon,4) - GIVEN constituents
% zmaj(ncon,... - Complex HC for GIVEN constituents/points
% SDtime(... - time expressed in Serial Days (see help datenum)
% Modes:
% Time series: zmaj(ncon,1), SDtime(nt,1) ->Output: dh(nt,1)
% Drift Track: zmaj(ncon,nt), SDtime(nt,1) ->Output: dh(nt,1)
% Map: zmaj(ncon,N,M),SDtime(1,1) ->Output: dh(N,M)
% -------------------------------------------------------------------
function [dh]=InferMinor(zmaj,cid,SDtime);
rad=pi/180;PP=282.8;dh=NaN;
% check size
[ncon,dum]=size(cid);
if dum~=4, fprintf('InferMinor call: Wrong constituents\n');return;end
[n1,n2,n3]=size(zmaj);
if n1~=ncon,zmaj=conj(zmaj');[n1,n2,n3]=size(zmaj);end
if n1~=ncon,fprintf('InferMinor call: Wrong zmaj size\n');return;end
TimeSeries=0;DriftTrack=0;TMap=0;
nt=length(SDtime);
%[n1 n2 n3 nt]
if n2==1 & n3==1,
TimeSeries=1;N=nt;M=1;
elseif n3==1 & n2>1 & nt==n2,
DriftTrack=1;N=nt;M=1;
elseif n2>1 & nt==1,
TMap=1; N=n2;M=n3;
else
fprintf('InferMinor call: Wrong zmaj/SDtime size\n');
return;
end
%
dh=zeros(N,M);
d0=datenum(1992,1,1); % corresponds to 48622mjd
d1=SDtime;
time=d1-d0; % time (days) relatively Jan 1 1992 GMT
time_mjd=48622+time;
[ncon,dum]=size(cid);
cid8=['q1 ';'o1 ';'p1 ';'k1 ';'n2 ';'m2 ';'s2 ';'k2 '];
% re-order zmaj to correspond to cid8
z8=zeros(8,N,M);
ni=0;
for i1=1:8
for j1=1:ncon
if cid(j1,:)==cid8(i1,:),
z8(i1,:,:)=zmaj(j1,:,:);
if(i1~=3 & i1~=8), ni=ni+1; end
end
end
end
if ni<6,return;end % not enough constituents for inference
%
zmin=zeros(18,N,M);
zmin(1,:,:) = 0.263 *z8(1,:,:) - 0.0252*z8(2,:,:); %2Q1
zmin(2,:,:) = 0.297 *z8(1,:,:) - 0.0264*z8(2,:,:); %sigma1
zmin(3,:,:) = 0.164 *z8(1,:,:) + 0.0048*z8(2,:,:); %rho1 +
zmin(4,:,:) = 0.0140*z8(2,:,:) + 0.0101*z8(4,:,:); %M1
zmin(5,:,:) = 0.0389*z8(2,:,:) + 0.0282*z8(4,:,:); %M1
zmin(6,:,:) = 0.0064*z8(2,:,:) + 0.0060*z8(4,:,:); %chi1
zmin(7,:,:) = 0.0030*z8(2,:,:) + 0.0171*z8(4,:,:); %pi1
zmin(8,:,:) =-0.0015*z8(2,:,:) + 0.0152*z8(4,:,:); %phi1
zmin(9,:,:) =-0.0065*z8(2,:,:) + 0.0155*z8(4,:,:); %theta1
zmin(10,:,:) =-0.0389*z8(2,:,:) + 0.0836*z8(4,:,:); %J1 +
zmin(11,:,:) =-0.0431*z8(2,:,:) + 0.0613*z8(4,:,:); %OO1 +
zmin(12,:,:) = 0.264 *z8(5,:,:) - 0.0253*z8(6,:,:); %2N2 +
zmin(13,:,:) = 0.298 *z8(5,:,:) - 0.0264*z8(6,:,:); %mu2 +
zmin(14,:,:) = 0.165 *z8(5,:,:) + 0.00487*z8(6,:,:); %nu2 +
zmin(15,:,:) = 0.0040*z8(6,:,:) + 0.0074*z8(7,:,:); %lambda2
zmin(16,:,:) = 0.0131*z8(6,:,:) + 0.0326*z8(7,:,:); %L2 +
zmin(17,:,:) = 0.0033*z8(6,:,:) + 0.0082*z8(7,:,:); %L2 +
zmin(18,:,:) = 0.0585*z8(7,:,:); %t2 +
%
hour = (time - floor(time))*24.D0;
t1 = 15.*hour;
t2 = 30.*hour;
[S,H,P,omega]=astrol(time_mjd); % nsites x nt
%
arg=zeros(18,N,M);
arg(1,:,:) = t1 - 4.*S + H + 2.*P - 90.; % 2Q1
arg(2,:,:) = t1 - 4.*S + 3.*H - 90.; % sigma1
arg(3,:,:) = t1 - 3.*S + 3.*H - P - 90.; % rho1
arg(4,:,:) = t1 - S + H - P + 90.; % M1
arg(5,:,:) = t1 - S + H + P + 90.; % M1
arg(6,:,:) = t1 - S + 3.*H - P + 90.; % chi1
arg(7,:,:) = t1 - 2.*H + PP - 90.; % pi1
arg(8,:,:) = t1 + 3.*H + 90.; % phi1
arg(9,:,:) = t1 + S - H + P + 90.; % theta1
arg(10,:,:) = t1 + S + H - P + 90.; % J1
arg(11,:,:) = t1 + 2.*S + H + 90.; % OO1
arg(12,:,:) = t2 - 4.*S + 2.*H + 2.*P; % 2N2
arg(13,:,:) = t2 - 4.*S + 4.*H; % mu2
arg(14,:,:) = t2 - 3.*S + 4.*H - P; % nu2
arg(15,:,:) = t2 - S + P + 180.D0; % lambda2
arg(16,:,:) = t2 - S + 2.*H - P + 180.D0; % L2
arg(17,:,:) = t2 - S + 2.*H + P; % L2
arg(18,:,:) = t2 - H + PP; % t2
%
% determine nodal corrections f and u
sinn = sin(omega*rad);
cosn = cos(omega*rad);
sin2n = sin(2.*omega*rad);
cos2n = cos(2.*omega*rad);
f = ones(18,N,M);
f(1,:,:) = sqrt((1.0 + 0.189*cosn - 0.0058*cos2n).^2 +...
(0.189*sinn - 0.0058*sin2n).^2);
f(2,:,:) = f(1,:,:);
f(3,:,:) = f(1);
f(4,:,:) = sqrt((1.0 + 0.185*cosn).^2 + (0.185*sinn).^2);
f(5,:,:) = sqrt((1.0 + 0.201*cosn).^2 + (0.201*sinn).^2);
f(6,:,:) = sqrt((1.0 + 0.221*cosn).^2 + (0.221*sinn).^2);
f(10,:,:) = sqrt((1.0 + 0.198*cosn).^2 + (0.198*sinn).^2);
f(11,:,:) = sqrt((1.0 + 0.640*cosn + 0.134*cos2n).^2 +...
(0.640*sinn + 0.134*sin2n).^2 );
f(12,:,:) = sqrt((1.0 - 0.0373*cosn).^2 + (0.0373*sinn).^2);
f(13,:,:) = f(12,:,:);
f(14,:,:) = f(12,:,:);
f(16,:,:) = f(12,:,:);
f(17,:,:) = sqrt((1.0 + 0.441*cosn).^2 + (0.441*sinn).^2);
%
u = zeros(18,N,M);
u(1,:,:) = atan2(0.189*sinn - 0.0058*sin2n,...
1.0 + 0.189*cosn - 0.0058*sin2n)/rad;
u(2,:,:) = u(1,:,:);
u(3,:,:) = u(1,:,:);
u(4,:,:) = atan2( 0.185*sinn, 1.0 + 0.185*cosn)/rad;
u(5,:,:) = atan2(-0.201*sinn, 1.0 + 0.201*cosn)/rad;
u(6,:,:) = atan2(-0.221*sinn, 1.0 + 0.221*cosn)/rad;
u(10,:,:) = atan2(-0.198*sinn, 1.0 + 0.198*cosn)/rad;
u(11,:,:) = atan2(-0.640*sinn - 0.134*sin2n,...
1.0 + 0.640*cosn + 0.134*cos2n)/rad;
u(12,:,:) = atan2(-0.0373*sinn, 1.0 - 0.0373*cosn)/rad;
u(13,:,:) = u(12,:,:);
u(14,:,:) = u(12,:,:);
u(16,:,:) = u(12,:,:);
u(17,:,:) = atan2(-0.441*sinn, 1.0 + 0.441*cosn)/rad;
% sum over all tides
% ------------------
for k=1:18
tmp=squeeze(real(zmin(k,:,:)).*f(k,:,:).*...
cos((arg(k,:,:) + u(k,:,:))*rad)-...
imag(zmin(k,:,:)).*f(k,:,:).*...
sin((arg(k,:,:)+u(k,:,:))*rad));
[n1,n2]=size(tmp);if n1==1,tmp=tmp';end
dh(:,:) = dh(:,:) + tmp;
end
return
function [mz,mu,mv] = Muv(hz)
% Given a rectangular bathymetry grid, construct masks for zeta,
% u and v nodes on a C-grid
% USAGE: [mz,mu,mv] = Muv(hz);
[n,m] = size(hz);
mz = hz > 0;
indx = [2:n 1];
indy = [2:m 1];
mu(indx,:) = mz.*mz(indx,:);
mv(:,indy) = mz.*mz(:,indy);
% calculates tidal ellipse parameters for the arrays of
% u and v - COMPLEX amplitudes of EW and NS currents of
% a given tidal constituent
% land should be set to 0 or NaN in u,v prior to calling tideEl
% usage: [umajor,uminor,uincl,uphase]=tideEl(u,v);
function [umajor,uminor,uincl,uphase]=tideEl(u,v);
% change to polar coordinates
% in Robin's was - + + -, this is as in Foreman's
t1p = (real (u) - imag(v));
t2p = (real (v) + imag(u));
t1m = (real (u) + imag(v));
t2m = (real (v) - imag(u));
% ap, am - amplitudes of positevely and negatively
% rotated vectors
ap = sqrt( (t1p.^2 + t2p.^2)) / 2.;
am = sqrt( (t1m.^2 + t2m.^2)) / 2.;
% ep, em - phases of positively and negatively rotating vectors
ep = atan2( t2p, t1p);
ep = ep + 2 * pi * (ep < 0.0);
ep = 180. * ep / pi;
em = atan2( t2m, t1m);
em = em + 2 * pi * (em < 0.0);
em = 180. * em / pi;
% determine the major and minor axes, phase and inclination using Foreman's formula
umajor = (ap + am);
uminor = (ap - am);
uincl = 0.5 * (em + ep);
uincl = uincl - 180. * (uincl > 180);
uphase = - 0.5*(ep-em) ;
uphase = uphase + 360. * (uphase < 0);
uphase = uphase - 360. * (uphase >= 360);
return
function [x,y] = XY(ll_lims,n,m);
% Usage: [x,y] = XY(ll_lims,n,m);
dx = (ll_lims(2)-ll_lims(1))/n;
dy = (ll_lims(4)-ll_lims(3))/m;
x = ll_lims(1)+dx/2:dx:ll_lims(2)-dx/2;
y = ll_lims(3)+dy/2:dy:ll_lims(4)-dy/2;
% Computes the basic astronomical mean longitudes s, h, p, N.
% Note N is not N', i.e. N is decreasing with time.
% These formulae are for the period 1990 - 2010, and were derived
% by David Cartwright (personal comm., Nov. 1990).
% time is UTC in decimal MJD.
% All longitudes returned in degrees.
% R. D. Ray Dec. 1990
% Non-vectorized version. Re-make for matlab by Lana Erofeeva, 2003
% usage: [s,h,p,N]=astrol(time)
% time, MJD
function [s,h,p,N]=astrol(time);
circle=360;
T = time - 51544.4993;
% mean longitude of moon
% ----------------------
s = 218.3164 + 13.17639648 * T;
% mean longitude of sun
% ---------------------
h = 280.4661 + 0.98564736 * T;
% mean longitude of lunar perigee
% -------------------------------
p = 83.3535 + 0.11140353 * T;
% mean longitude of ascending lunar node
% --------------------------------------
N = 125.0445D0 - 0.05295377D0 * T;
%
s = mod(s,circle);
h = mod(h,circle);
p = mod(p,circle);
N = mod(N,circle);
%
if s<0, s = s + circle; end
if h<0, h = h + circle; end
if p<0, p = p + circle; end
if N<0, N = N + circle; end
return
% returns Flag=0, if ModName, GridName & type OK
% Flag>0, if they are not or files do not exist
% Flag<0, if not possible to define
% usage: [Flag]=checkTypeName(ModName,GridName,type);
function [Flag]=checkTypeName(ModName,GridName,type);
Flag=0;
% check if files exist
if exist(ModName,'file')==0,
fprintf('File %s NOT found\n',ModName); Flag=1;
return
end
if exist(GridName,'file')==0,
fprintf('File %s NOT found\n',GridName); Flag=1;
return
end
type=deblank(type);btype=deblank(type(end:-1:1));
type=btype(end:-1:1);type=type(1:1);
% Check type
if type~='z' & type~='U'& type~='u'& type~='V' & type~='v',
fprintf('WRONG TYPE %s: should be one of: ''z'',''U'',''u'',''V'',''v''\n',type);
Flag=2;
return
end
% Check type/name correspondence
i1=findstr(ModName,'/');
if isempty(i1)>0,i1=1;else i1=i1(end)+1;end
if ModName(i1:i1)=='h',
if type=='U'| type=='u'|type=='V'| type=='v',
fprintf('For file %s only legal type is ''z'' (%s given)\n',...
ModName(i1:end),type);
Flag=3;
return
end
elseif ModName(i1:i1)=='u' | ModName(i1:i1)=='U',
if type=='z',
fprintf('For file %s legal types are: ''U'',''u'',''V'',''v'' (%s given)\n',...
ModName(i1:end),type);
Flag=4;
return
end
else
fprintf('WARNING: Model name %s does not correspond TMD convention:\n',ModName(i1:end));
fprintf('Can not check type %s:...\n',type);
Flag=-1;
end
return
function [l]=check_dim(ModName,n,m,n1,m1);
l=1;
if n~=n1 | m~=m1,
fprintf('ERROR:\n');
fprintf('Control file %s contains inconsistent files:\n',ModName);
fprintf('Grid size: %d x %d \n',n,m);
fprintf('Elevation/transport size: %d x %d \n',n1,m1);
l=0;
end
return;
% constit returns amplitude, phase, frequency,alpha, species for
% a tidal constituent identified by a 4 character string
% This version returns zeros for all phases
% Usage: [ispec,amp,ph,omega,alpha,constitNum] = constit(c);
function [ispec,amp,ph,omega,alpha,constitNum] = constit(c);
ispec=-1;
amp=0;ph=0;omega=0;alpha=0;constitNum=0;
if( nargout < 6)
fprintf('Number of output arguments to constit.m does not agree with current usage\n');
end
% c : character string ID of constituents that the program knows about
% all other variables are in the same order
c_data = ['m2 ';'s2 ';'k1 ';'o1 '; ...
'n2 ';'p1 ';'k2 ';'q1 '; ...
'2n2 ';'mu2 ';'nu2 ';'l2 '; ...
't2 ';'j1 ';'no1 ';'oo1 '; ...
'rho1';'mf ';'mm ';'ssa ';'m4 ';'ms4 ';'mn4 ' ];
constitNum_data = [1,2,5,6,...
3,7,4,8,...
0,0,0,0,...
0,0,0,0,...
0,0,0,0,0,0,0];
% ispec : species type (spherical harmonic dependence of quadropole potential)
ispec_data = [2;2;1;1; ...
2;1;2;1; ...
2;2;2;2; ...
2;1;1;1; ...
1;0;0;0;0;0;0];
% alpha : loading love number ... frequncy dependance here is suspect
alpha_data = [ 0.693;0.693;0.736;0.695; ...
0.693;0.706;0.693;0.695; ...
0.693;0.693;0.693;0.693; ...
0.693;0.693;0.693;0.693; ...
0.693;0.693;0.693;0.693;0.693;0.693;0.693];
% omega : frequencies
omega_data = [1.405189e-04;1.454441e-04;7.292117e-05;6.759774e-05; ...
1.378797e-04;7.252295e-05;1.458423e-04;6.495854e-05; ...
1.352405e-04;1.355937e-04;1.382329e-04;1.431581e-04; ...
1.452450e-04;7.556036e-05;7.028195e-05;7.824458e-05; ...
6.531174e-05;0.053234e-04;0.026392e-04;0.003982e-04; ...
2.810377e-04;2.859630e-04;2.783984e-04];
% phase : Astronomical arguments (relative to t0 = 1 Jan 0:00 1992]
% Richrad Ray subs: "arguments" and "astrol"
phase_data = [ 1.731557546;0.000000000;0.173003674;1.558553872;...
6.050721243;6.110181633;3.487600001;5.877717569;
4.086699633;3.463115091;5.427136701;0.553986502;
0.052841931;2.137025284;2.436575100;1.929046130;
5.254133027;1.756042456;1.964021610;3.487600001;
3.463115091;1.731557546;1.499093481];
% amp : amplitudes
%amp_data = [0.242334;0.112743;0.141565;0.100661; ...
amp_data = [0.2441 ;0.112743;0.141565;0.100661; ...
0.046397;0.046848;0.030684;0.019273; ...
0.006141;0.007408;0.008811;0.006931; ...
0.006608;0.007915;0.007915;0.004338; ...
0.003661;0.042041;0.022191;0.019567;0;0;0];
nspecies = length(ispec_data);
nl = length(c);
kk = 0;
for k = 1:nspecies
if(lower(c) == c_data(k,1:nl)) ,
kk = k;
break
end
end
if(kk == 0 ),
ispec = -1;
return
else
constitNum = constitNum_data(kk);
ispec = ispec_data(kk);
amp = amp_data(kk);
alpha = alpha_data(kk);
omega = omega_data(kk);
% ph = 0;
ph = phase_data(kk);
end
return
% reads a grid file in matlab
% USAGE: [ll_lims,hz,mz,iob,dt] = grd_in(cfile);
% Location: C:\Toolboxes\TMD\FUNCTIONS
function [ll_lims,hz,mz,iob,dt] = grd_in(cfile);
fid = fopen(cfile,'r','b');
fseek(fid,4,'bof');
n = fread(fid,1,'long');
m = fread(fid,1,'long');
lats = fread(fid,2,'float');
lons = fread(fid,2,'float');
dt = fread(fid,1,'float');
if(lons(1) < 0) & (lons(2) < 0 ) & dt>0,
lons = lons + 360;
end
ll_lims = [lons ; lats ];
%fprintf('Time step (sec): %10.1f\n',dt);
nob = fread(fid,1,'long');
if nob == 0,
fseek(fid,20,'cof');
iob = [];
else
fseek(fid,8,'cof');
iob = fread(fid,[2,nob],'long');
fseek(fid,8,'cof');
end
hz = fread(fid,[n,m],'float');
fseek(fid,8,'cof');
mz = fread(fid,[n,m],'long');
fclose(fid);
% outputs a grid file in matlab
% USAGE: grd_out(cfile,ll_lims,hz,mz,iob,dt);
function grd_out(cfile,ll_lims,hz,mz,iob,dt);
% open this way for files to be read on Unix machine
fid = fopen(cfile,'w','b');
[dum,nob] = size(iob);
[n,m] = size(hz);
reclen = 32;
fwrite(fid,reclen,'long');
fwrite(fid,n,'long');
fwrite(fid,m,'long');
fwrite(fid,ll_lims(3:4),'float');
fwrite(fid,ll_lims(1:2),'float');
fwrite(fid,dt,'float');
fwrite(fid,nob,'long');
fwrite(fid,reclen,'long');
%
if nob == 0,
fwrite(fid,4,'long');
fwrite(fid,0,'long');
fwrite(fid,4,'long');
else
reclen=8*nob;
fwrite(fid,reclen,'long');
fwrite(fid,iob,'long');
fwrite(fid,reclen,'long');
end
%
reclen = 4*n*m;
fwrite(fid,reclen,'long');
fwrite(fid,hz,'float');
fwrite(fid,reclen,'long');
fwrite(fid,reclen,'long');
fwrite(fid,mz,'long');
fwrite(fid,reclen,'long');
fclose(fid);
return
function [h,th_lim,ph_lim] = h_in(cfile,ic)
%USAGE: [h,th_lim,ph_lim] = h_in(cfile,ic);
% reads in elevation for constituent # ic in file cfile
fid = fopen(cfile,'r','b');
ll = fread(fid,1,'long');
nm = fread(fid,3,'long');
n=nm(1);
m = nm(2);
nc = nm(3);
th_lim = fread(fid,2,'float');
ph_lim = fread(fid,2,'float');
%nskip = (ic-1)*(nm(1)*nm(2)*8+8) + 8 + 4*nc;
nskip = (ic-1)*(nm(1)*nm(2)*8+8) + 8 + ll - 28;
fseek(fid,nskip,'cof');
htemp = fread(fid,[2*n,m],'float');
h = htemp(1:2:2*n-1,:)+i*htemp(2:2:2*n,:);
fclose(fid);
function [] = h_out(cfile,h,th_lim,ph_lim,c_ids)
% Usage: [] = h_out(cfile,h,th_lim,ph_lim,c_ids)
%WARNING!!!! C_IDS HAVE TO BE 4 CHARACTER STRINGS !!!!!!!!!
% writes elevation file in standard format
fid = fopen(cfile,'w','b');
[n,m,nc] = size(h);
% length of header: allow for 4 character long c_id strings
llHead = 4*(7+nc);
fwrite(fid,llHead,'long');
fwrite(fid,n,'long');
fwrite(fid,m,'long');
fwrite(fid,nc,'long');
fwrite(fid,th_lim,'float');
fwrite(fid,ph_lim,'float');
fwrite(fid,c_ids,'char');
fwrite(fid,llHead,'long');
htemp = zeros(2*n,m);
llConstit = 8*n*m;
for ic = 1:nc
fwrite(fid,llConstit,'long');
htemp(1:2:2*n-1,:) = real(h(:,:,ic));
htemp(2:2:2*n,:) = imag(h(:,:,ic));
fwrite(fid,htemp,'float');
fwrite(fid,llConstit,'long');
end
fclose(fid);
% function to predict tidal elevation fields at "time"
% using harmonic constants
% INPUT: time (days) relatively Jan 1, 1992 (48622mjd)
% con(nc,4) - char*4 tidal constituent IDs
% hc(n,m,nc) - harmonic constant vector (complex)
% OUTPUT:hhat - tidal field at time
%
% Nodal corrections included
%
% usage: [hhat]=harp(time,hc,con);