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 |