matlab_class/gcmfaces_calc/gcmfaces_phitheta.m

function [PHI,phi,varargout]=gcmfaces_phitheta(T,dT,Theta,varargin);
% GCMFACES_PHITHETA parameterized probability distribution diagnostic
%
%     [PHI,phi]=gcmfaces_phitheta(T,dT,Theta) computes the 
%     probabilities that theta<=Theta (PHI) and theta=Theta (phi) 
%     assuming that the theta probably density is homogeneous around T 
%     (i.e., equal to 1/2/dT in the [T-dT T+dT] range).
%
%     [PHI,phi,dPHIdt]=gcmfaces_phitheta(T,dT,Theta,H) additionally
%     computes the -phi*H rate of change in PHI associated with the 
%     action of H.
%
%     The T input can be from the double or gcmfaces class whereas dT and
%     Theta must be of the double class. While dT must be a scalar value, 
%     T and/or Theta can have one or several non-singleton dimensions. 
%     In these cases T and Theta dimensions get coumpounded except 
%     for singleton dimensions when T or Theta is a column vector.
%
% examples:
% 
%     T=[3:0.01:10]; dT=0.5; Theta=6; 
%     [PHI,phi]=gcmfaces_phitheta(T,dT,Theta);
%
%     T=[3:0.01:10]'; dT=0.5; Theta=[-3:0.1:30]'; 
%     [PHI,phi]=gcmfaces_phitheta(T,dT,Theta);
%
%     dir0='release2_climatology/nctiles_climatology/'; rhocp=1029.*3994;
%     SST=read_nctiles([dir0 'THETA/THETA'],'THETA',[],1);
%     oceQnet=read_nctiles([dir0 'oceQnet/oceQnet'],'oceQnet');
%     [PHI,phi,dPHIdt]=gcmfaces_phitheta(SST,0.5,6,1/rhocp*oceQnet);

if ~isempty(which('gcmfaces_global')); gcmfaces_global; end;
if isempty(which('gcmfaces_global'))&isa(T,'gcmfaces');
    error('gcmfaces is missing from Matlab path');
end;

if nargin<3; error('incorrect input parameter specification'); end;

%identify input dimensions
if isa(T,'gcmfaces'); ndimT=size(T{1}); else; ndimT=size(T); end;
if ndimT(end)==1; ndimT=length(ndimT)-1; else; ndimT=length(ndimT); end;
ndimTheta=length(size(Theta));

%replicate T if needed to coumpound dimensions
tmp1=[ones(1,ndimT) size(Theta)];
T=repmat(T,tmp1);
for ii=1:nargin-3;
    varargout{ii}=repmat(varargin{ii},tmp1);
end;

%replicate Theta if needed to coumpound dimensions
tmp1=[[1:ndimT]+ndimTheta [1:ndimTheta]];
Theta=permute(Theta,tmp1);
if isa(T,'gcmfaces');
    tmp1=size(T{1});
    tmp1=[mygrid.ioSize tmp1(3:ndimT) ones(1,ndimTheta)];
    Theta=convert2gcmfaces(repmat(Theta,tmp1));
else;
    tmp1=size(T);
    tmp1=[tmp1([1:ndimT]) ones(1,ndimTheta)];
    Theta=repmat(Theta,tmp1);
end;

%compte PHI and phi
PHI=(Theta-T+dT)/(2*dT);
phi=0*T+1/(2*dT);

PHI(PHI<0)=0; PHI(PHI>1)=1;
phi(PHI==0)=0; phi(PHI==1)=0;

%compute transformation rates
for ii=1:nargin-3;
    varargout{ii}=-varargout{ii}.*phi;
end;


1	gforget	1.1	function [PHI,phi,varargout]=gcmfaces_phitheta(T,dT,Theta,varargin);
2			% GCMFACES_PHITHETA parameterized probability distribution diagnostic
3			%
4			% [PHI,phi]=gcmfaces_phitheta(T,dT,Theta) computes the
5			% probabilities that theta<=Theta (PHI) and theta=Theta (phi)
6			% assuming that the theta probably density is homogeneous around T
7			% (i.e., equal to 1/2/dT in the [T-dT T+dT] range).
8			%
9			% [PHI,phi,dPHIdt]=gcmfaces_phitheta(T,dT,Theta,H) additionally
10			% computes the -phi*H rate of change in PHI associated with the
11			% action of H.
12			%
13			% The T input can be from the double or gcmfaces class whereas dT and
14			% Theta must be of the double class. While dT must be a scalar value,
15			% T and/or Theta can have one or several non-singleton dimensions.
16			% In these cases T and Theta dimensions get coumpounded except
17			% for singleton dimensions when T or Theta is a column vector.
18			%
19			% examples:
20			%
21			% T=[3:0.01:10]; dT=0.5; Theta=6;
22			% [PHI,phi]=gcmfaces_phitheta(T,dT,Theta);
23			%
24			% T=[3:0.01:10]'; dT=0.5; Theta=[-3:0.1:30]';
25			% [PHI,phi]=gcmfaces_phitheta(T,dT,Theta);
26			%
27			% dir0='release2_climatology/nctiles_climatology/'; rhocp=1029.*3994;
28			% SST=read_nctiles([dir0 'THETA/THETA'],'THETA',[],1);
29			% oceQnet=read_nctiles([dir0 'oceQnet/oceQnet'],'oceQnet');
30			% [PHI,phi,dPHIdt]=gcmfaces_phitheta(SST,0.5,6,1/rhocp*oceQnet);
31
32			if ~isempty(which('gcmfaces_global')); gcmfaces_global; end;
33			if isempty(which('gcmfaces_global'))&isa(T,'gcmfaces');
34			error('gcmfaces is missing from Matlab path');
35			end;
36
37			if nargin<3; error('incorrect input parameter specification'); end;
38
39			%identify input dimensions
40			if isa(T,'gcmfaces'); ndimT=size(T{1}); else; ndimT=size(T); end;
41			if ndimT(end)==1; ndimT=length(ndimT)-1; else; ndimT=length(ndimT); end;
42			ndimTheta=length(size(Theta));
43
44			%replicate T if needed to coumpound dimensions
45			tmp1=[ones(1,ndimT) size(Theta)];
46			T=repmat(T,tmp1);
47			for ii=1:nargin-3;
48			varargout{ii}=repmat(varargin{ii},tmp1);
49			end;
50
51			%replicate Theta if needed to coumpound dimensions
52			tmp1=[[1:ndimT]+ndimTheta [1:ndimTheta]];
53			Theta=permute(Theta,tmp1);
54			if isa(T,'gcmfaces');
55			tmp1=size(T{1});
56			tmp1=[mygrid.ioSize tmp1(3:ndimT) ones(1,ndimTheta)];
57			Theta=convert2gcmfaces(repmat(Theta,tmp1));
58			else;
59			tmp1=size(T);
60			tmp1=[tmp1([1:ndimT]) ones(1,ndimTheta)];
61			Theta=repmat(Theta,tmp1);
62			end;
63
64			%compte PHI and phi
65			PHI=(Theta-T+dT)/(2*dT);
66			phi=0T+1/(2dT);
67
68			PHI(PHI<0)=0; PHI(PHI>1)=1;
69			phi(PHI==0)=0; phi(PHI==1)=0;
70
71			%compute transformation rates
72			for ii=1:nargin-3;
73			varargout{ii}=-varargout{ii}.*phi;
74			end;
75
76