Warning, /verification/fizhi-gridalt-hs/matlab/fromPhys2Dyn.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
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