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	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