profilesMatlabProcessing/profiles_misc/model_interp.m

function [MITprofCur]=model_interp(MITprofCur,model,varargin)
% function [MITprofCur]=model_interp(MITprofCur,model,varargin)
%   interpolate model hydrographic fields at the location of given profiles
%
%   MITprofCur: structure containing profile information
%   model is a string used to select the model to be loaded
%       'OCCA' : ECCOv4 grid + OCCA atlas
%       'SOSE59' : SOSE59 grid + atlas

% load model
[grid,atlas]=model_load(model);

% protect global grid variables
global mygrid mytri MYBASININDEX
mygrid1=mygrid; mytri1=mytri; MYBASININDEX1=MYBASININDEX;
mygrid=grid.mygrid; mytri=grid.mytri; MYBASININDEX=grid.MYBASININDEX;

% compute index of closest profile
MITprofCur.prof_point=gcmfaces_bindata(MITprofCur.prof_lon,MITprofCur.prof_lat);
[MITprofCur.ii,MITprofCur.jj]=gcmfaces_bindata(MITprofCur.prof_lon,MITprofCur.prof_lat);

% temporal index
switch model
    case 'SOSE59'
        MITprofCur.imonth=MITprofCur.ii*0+1;
    case 'OCCA'
        MITprofCur.imonth=floor(MITprofCur.prof_YYYYMMDD/1e2)-1e2*floor(MITprofCur.prof_YYYYMMDD/1e4);
    otherwise
        error('not a valid model string');
end

% extract profiles
nk=length(mygrid.RC);kk=ones(1,nk);
np=length(MITprofCur.ii);pp=ones(np,1);
ind2prof=sub2ind(size(atlas.T{1}),MITprofCur.ii*kk,MITprofCur.jj*kk,pp*[1:nk],MITprofCur.imonth*kk);

I=find(~isnan(ind2prof));
Tatlas=ind2prof*NaN;Tatlas(I)=atlas.T{1}(ind2prof(I)); Tatlas(Tatlas==0)=NaN;
Satlas=ind2prof*NaN;Satlas(I)=atlas.S{1}(ind2prof(I)); Satlas(Satlas==0)=NaN;

warning off
Tatlas=interp1(-mygrid.RC',Tatlas',MITprofCur.prof_depth)';
Satlas=interp1(-mygrid.RC',Satlas',MITprofCur.prof_depth)';
warning on

% record values
MITprofCur=setfield(MITprofCur,['prof_T_' model],Tatlas');
MITprofCur=setfield(MITprofCur,['prof_S_' model],Satlas');

% put back global grid values
mygrid=mygrid1; mytri=mytri1; MYBASININDEX=MYBASININDEX1;

1	function [MITprofCur]=model_interp(MITprofCur,model,varargin)
2	% function [MITprofCur]=model_interp(MITprofCur,model,varargin)
3	% interpolate model hydrographic fields at the location of given profiles
4	%
5	% MITprofCur: structure containing profile information
6	% model is a string used to select the model to be loaded
7	% 'OCCA' : ECCOv4 grid + OCCA atlas
8	% 'SOSE59' : SOSE59 grid + atlas
9
10	% load model
11	[grid,atlas]=model_load(model);
12
13	% protect global grid variables
14	global mygrid mytri MYBASININDEX
15	mygrid1=mygrid; mytri1=mytri; MYBASININDEX1=MYBASININDEX;
16	mygrid=grid.mygrid; mytri=grid.mytri; MYBASININDEX=grid.MYBASININDEX;
17
18	% compute index of closest profile
19	MITprofCur.prof_point=gcmfaces_bindata(MITprofCur.prof_lon,MITprofCur.prof_lat);
20	[MITprofCur.ii,MITprofCur.jj]=gcmfaces_bindata(MITprofCur.prof_lon,MITprofCur.prof_lat);
21
22	% temporal index
23	switch model
24	case 'SOSE59'
25	MITprofCur.imonth=MITprofCur.ii*0+1;
26	case 'OCCA'
27	MITprofCur.imonth=floor(MITprofCur.prof_YYYYMMDD/1e2)-1e2*floor(MITprofCur.prof_YYYYMMDD/1e4);
28	otherwise
29	error('not a valid model string');
30	end
31
32	% extract profiles
33	nk=length(mygrid.RC);kk=ones(1,nk);
34	np=length(MITprofCur.ii);pp=ones(np,1);
35	ind2prof=sub2ind(size(atlas.T{1}),MITprofCur.iikk,MITprofCur.jjkk,pp[1:nk],MITprofCur.imonthkk);
36
37	I=find(~isnan(ind2prof));
38	Tatlas=ind2prof*NaN;Tatlas(I)=atlas.T{1}(ind2prof(I)); Tatlas(Tatlas==0)=NaN;
39	Satlas=ind2prof*NaN;Satlas(I)=atlas.S{1}(ind2prof(I)); Satlas(Satlas==0)=NaN;
40
41	warning off
42	Tatlas=interp1(-mygrid.RC',Tatlas',MITprofCur.prof_depth)';
43	Satlas=interp1(-mygrid.RC',Satlas',MITprofCur.prof_depth)';
44	warning on
45
46	% record values
47	MITprofCur=setfield(MITprofCur,['prof_T_' model],Tatlas');
48	MITprofCur=setfield(MITprofCur,['prof_S_' model],Satlas');
49
50	% put back global grid values
51	mygrid=mygrid1; mytri=mytri1; MYBASININDEX=MYBASININDEX1;
52