fizhi-gridalt-hs/matlab/fromPhys2Dyn.m

cc5fb01035 Jean*0001 function [vDy,klev,dpD]=fromPhys2Dyn(vPh,dpP,inpPres,delR);
                0002 % [vDy,klev]=fromPhys2Dyn(vPh,dpP,Dp,delR); (1rst call)
                0003 % [vDy,klev]=fromPhys2Dyn(vPh,dpP,klev);    (next calls)
                0004 %
                0005 %  vPh :: 3.D input field on Physics grid, to put on Dynamics grid.
                0006 %  dpP :: delta.P field (3D) for Physics grid.
                0007 %  Dp  :: Total air-column thickness (2.D, pressure units)
                0008 % delR :: Dynamics delta.P (vector) (pressure units)
                0009 % klev :: index of the Dynamics level for each physics grid level (integer array, 3.D)
                0010 %  vDy :: 3.D output field on Dynamics grid.
                0011 %  dpD :: delta.P field (3D) for dynamics grid (~hFacC*delR) [usefull for checking]
                0012 %
                0013 %- 1rst call:
                0014 % Dp=rdmds('Depth'); delR (from file "data"); dpP from diagnostics output ('DPPHYS  ')
                0015 % [vDy,klev]=fromPhys2Dyn(vPh,dpP,Dp,delR);
                0016 %- next calls:
                0017 % [vDy]=fromPhys2Dyn(vPh,dpP,klev);
                0018 if nargin < 4,
                0019  first=0;
                0020  klev=inpPres;
                0021  if size(klev) ~= size(vPh),
                0022    fprintf(' Mismatch in size vPh:');fprintf(' %i',size(vPh));
                0023    fprintf(' & klev:');fprintf(' %i',size(klev));
                0024    fprintf('\n => stop'); vDy=0; klev=0; return
                0025  end
                0026  nr=max(max(max(klev)));
                0027 else
                0028  first=1;
                0029  Dp=inpPres;
                0030  nr=length(delR);
                0031 end
                0032
                0033  ncx=size(dpP,1); nc=size(dpP,2); nP=size(dpP,3); nPg=ncx*nc;
                0034  epsil=1.e-6;
                0035
                0036 %- find corresponding levels:
                0037 if first == 1,
                0038  fprintf(' Compute klev: ');
                0039  TimeT0=clock;
                0040  var=Dp./sum(dpP,3);
                0041  for k=1:nP, dpP(:,:,k)=var.*dpP(:,:,k); end
                0042  pp1=cumsum(dpP,3); pp1=reshape(pp1,nPg,nP);
                0043  pDw=zeros(1,nr+1); pDw(nr:-1:1)=cumsum(delR(nr:-1:1));
                0044  klev=ones(nPg,nP);
                0045  for ij=1:nPg,
                0046   for k=1:nP,
                0047    [K]=find(pDw>pp1(ij,k)-epsil);
                0048    klev(ij,k)=max(K);
                0049   end
                0050  end
                0051  klev=reshape(klev,ncx,nc,nP);
                0052  TimeT1=clock;
                0053  fprintf(' <- done (Time= %8.3f s)\n',etime(TimeT1,TimeT0));
                0054 else
                0055  if size(klev) ~= size(vPh),
                0056    fprintf(' Mismatch in size vPh:');fprintf(' %i',size(vPh));
                0057    fprintf(' & klev:');fprintf(' %i',size(klev));
                0058    fprintf('\n => stop'); vDy=0; klev=0; return
                0059  end
                0060 end
                0061
                0062  fprintf(' Average on Dyn.Grid: ');
                0063  TimeT2=clock;
                0064  vDy=zeros(ncx,nc,nr);
                0065  dpD=zeros(ncx,nc,nr);
                0066  for k=1:nr,
                0067    var=dpP; var(find(klev~=k))=0;
                0068    dpDyn=sum(var,3);
                0069    dpD(:,:,k)=dpDyn;
                0070    dpDyn(find(dpDyn==0))=-1;
                0071    var=var.*vPh;
                0072    vDyn=sum(var,3)./dpDyn;
                0073    vDyn(find(dpDyn==-1))=0;
                0074    vDy(:,:,k)=vDyn;
                0075  end
                0076  TimeT3=clock;
                0077  fprintf(' <- done (Time= %8.3f s)\n',etime(TimeT3,TimeT2));
                0078
                0079 return