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	roquet	1.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	roquet	1.2
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	roquet	1.1
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