fizhi-gridalt-hs/matlab/fromPhys2Dyn.m

function [vDy,klev,dpD]=fromPhys2Dyn(vPh,dpP,inpPres,delR);
% [vDy,klev]=fromPhys2Dyn(vPh,dpP,Dp,delR); (1rst call)
% [vDy,klev]=fromPhys2Dyn(vPh,dpP,klev);    (next calls)
%
%  vPh :: 3.D input field on Physics grid, to put on Dynamics grid.
%  dpP :: delta.P field (3D) for Physics grid.
%  Dp  :: Total air-column thickness (2.D, pressure units)
% delR :: Dynamics delta.P (vector) (pressure units)
% klev :: index of the Dynamics level for each physics grid level (integer array, 3.D)
%  vDy :: 3.D output field on Dynamics grid.
%  dpD :: delta.P field (3D) for dynamics grid (~hFacC*delR) [usefull for checking]
%
%- 1rst call:
% Dp=rdmds('Depth'); delR (from file "data"); dpP from diagnostics output ('DPPHYS  ')
% [vDy,klev]=fromPhys2Dyn(vPh,dpP,Dp,delR);
%- next calls:
% [vDy]=fromPhys2Dyn(vPh,dpP,klev);
% $Header:  $

if nargin < 4,
 first=0;
 klev=inpPres;
 if size(klev) ~= size(vPh),
   fprintf(' Mismatch in size vPh:');fprintf(' %i',size(vPh));
   fprintf(' & klev:');fprintf(' %i',size(klev)); 
   fprintf('\n => stop'); vDy=0; klev=0; return
 end
 nr=max(max(max(klev)));
else
 first=1;
 Dp=inpPres;
 nr=length(delR);
end

 ncx=size(dpP,1); nc=size(dpP,2); nP=size(dpP,3); nPg=ncx*nc;
 epsil=1.e-6;

%- find corresponding levels:
if first == 1,
 fprintf(' Compute klev: ');
 TimeT0=clock;
 var=Dp./sum(dpP,3);
 for k=1:nP, dpP(:,:,k)=var.*dpP(:,:,k); end
 pp1=cumsum(dpP,3); pp1=reshape(pp1,nPg,nP);
 pDw=zeros(1,nr+1); pDw(nr:-1:1)=cumsum(delR(nr:-1:1));
 klev=ones(nPg,nP);
 for ij=1:nPg,
  for k=1:nP,
   [K]=find(pDw>pp1(ij,k)-epsil);
   klev(ij,k)=max(K);
  end
 end
 klev=reshape(klev,ncx,nc,nP);
 TimeT1=clock;
 fprintf(' <- done (Time= %8.3f s)\n',etime(TimeT1,TimeT0));
else
 if size(klev) ~= size(vPh),
   fprintf(' Mismatch in size vPh:');fprintf(' %i',size(vPh));
   fprintf(' & klev:');fprintf(' %i',size(klev)); 
   fprintf('\n => stop'); vDy=0; klev=0; return
 end
end

 fprintf(' Average on Dyn.Grid: ');
 TimeT2=clock;
 vDy=zeros(ncx,nc,nr);
 dpD=zeros(ncx,nc,nr);
 for k=1:nr,
   var=dpP; var(find(klev~=k))=0;
   dpDyn=sum(var,3);
   dpD(:,:,k)=dpDyn;
   dpDyn(find(dpDyn==0))=-1;
   var=var.*vPh;
   vDyn=sum(var,3)./dpDyn;
   vDyn(find(dpDyn==-1))=0;
   vDy(:,:,k)=vDyn;
 end
 TimeT3=clock;
 fprintf(' <- done (Time= %8.3f s)\n',etime(TimeT3,TimeT2));

return
1	function [vDy,klev,dpD]=fromPhys2Dyn(vPh,dpP,inpPres,delR);
2	% [vDy,klev]=fromPhys2Dyn(vPh,dpP,Dp,delR); (1rst call)
3	% [vDy,klev]=fromPhys2Dyn(vPh,dpP,klev); (next calls)
4	%
5	% vPh :: 3.D input field on Physics grid, to put on Dynamics grid.
6	% dpP :: delta.P field (3D) for Physics grid.
7	% Dp :: Total air-column thickness (2.D, pressure units)
8	% delR :: Dynamics delta.P (vector) (pressure units)
9	% klev :: index of the Dynamics level for each physics grid level (integer array, 3.D)
10	% vDy :: 3.D output field on Dynamics grid.
11	% dpD :: delta.P field (3D) for dynamics grid (~hFacC*delR) [usefull for checking]
12	%
13	%- 1rst call:
14	% Dp=rdmds('Depth'); delR (from file "data"); dpP from diagnostics output ('DPPHYS ')
15	% [vDy,klev]=fromPhys2Dyn(vPh,dpP,Dp,delR);
16	%- next calls:
17	% [vDy]=fromPhys2Dyn(vPh,dpP,klev);
18	% $Header: $
19
20	if nargin < 4,
21	first=0;
22	klev=inpPres;
23	if size(klev) ~= size(vPh),
24	fprintf(' Mismatch in size vPh:');fprintf(' %i',size(vPh));
25	fprintf(' & klev:');fprintf(' %i',size(klev));
26	fprintf('\n => stop'); vDy=0; klev=0; return
27	end
28	nr=max(max(max(klev)));
29	else
30	first=1;
31	Dp=inpPres;
32	nr=length(delR);
33	end
34
35	ncx=size(dpP,1); nc=size(dpP,2); nP=size(dpP,3); nPg=ncx*nc;
36	epsil=1.e-6;
37
38	%- find corresponding levels:
39	if first == 1,
40	fprintf(' Compute klev: ');
41	TimeT0=clock;
42	var=Dp./sum(dpP,3);
43	for k=1:nP, dpP(:,:,k)=var.*dpP(:,:,k); end
44	pp1=cumsum(dpP,3); pp1=reshape(pp1,nPg,nP);
45	pDw=zeros(1,nr+1); pDw(nr:-1:1)=cumsum(delR(nr:-1:1));
46	klev=ones(nPg,nP);
47	for ij=1:nPg,
48	for k=1:nP,
49	[K]=find(pDw>pp1(ij,k)-epsil);
50	klev(ij,k)=max(K);
51	end
52	end
53	klev=reshape(klev,ncx,nc,nP);
54	TimeT1=clock;
55	fprintf(' <- done (Time= %8.3f s)\n',etime(TimeT1,TimeT0));
56	else
57	if size(klev) ~= size(vPh),
58	fprintf(' Mismatch in size vPh:');fprintf(' %i',size(vPh));
59	fprintf(' & klev:');fprintf(' %i',size(klev));
60	fprintf('\n => stop'); vDy=0; klev=0; return
61	end
62	end
63
64	fprintf(' Average on Dyn.Grid: ');
65	TimeT2=clock;
66	vDy=zeros(ncx,nc,nr);
67	dpD=zeros(ncx,nc,nr);
68	for k=1:nr,
69	var=dpP; var(find(klev~=k))=0;
70	dpDyn=sum(var,3);
71	dpD(:,:,k)=dpDyn;
72	dpDyn(find(dpDyn==0))=-1;
73	var=var.*vPh;
74	vDyn=sum(var,3)./dpDyn;
75	vDyn(find(dpDyn==-1))=0;
76	vDy(:,:,k)=vDyn;
77	end
78	TimeT3=clock;
79	fprintf(' <- done (Time= %8.3f s)\n',etime(TimeT3,TimeT2));
80
81	return