Warning, /verification/global_ocean.cs32x15/input.thsice/check_mom_budget.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 3158eb90 on 2019-08-09 23:24:37 UTC
3158eb90e6 Jean*0001 %- Uses few MITgcm matlab scripts that can be found in MITgcm/utils/matlab/
0002 % (rdmds.m & load_grid.m) and in MITgcm/utils/matlab/cs_grid/ (split_C_cub.m)
0003 % Suggestion: uncomment the 2 folowing lines to access them:
0004 % addpath ../../../utils/matlab/
0005 % addpath ../../../utils/matlab/cs_grid
a093b9b1a0 Jean*0006
0007 if size(who('kpr'),1) > 0,
0008 fprintf('kpr is defined and = %i \n',kpr);
0009 else
0010 fprintf('kpr undefined ; set to 1 \n'); kpr=1 ;
0011 end
0012
0013 rDir='../tr_run.thsice/'; nit=36010; deltaT=86400;
0014 %rDir='../tr_run.viscA4/'; nit=86405; deltaT=3600;
0015
0016 %-----
0017 %gDir='grid64/';
0018 gDir=rDir; %- for grid files
0019
0020 ii=strfind(rDir,'/'); i1=1; i2=length(rDir)-1;
0021 if length(ii) > 1, i1=1+ii(1); i2=ii(2)-1; end
0022 titexp=[' Glob.Oce ',rDir(i1:i2),' ; ']; titexp=strrep(titexp,'_','\_');
0023 MnIter=30*86400/deltaT;
0024 tit0=sprintf(' mn=%2i',nit/MnIter);
0025
0026 namF='momDiag'; namFs='srfDiag';
0027 %namFs='surfDiag';
0028 tit0=namF;
0029
0030 %--------------------------------------------------
0031 if kpr > 0,
0032 %-- Read in grif-files. Needed with older output to:
0033 % a) compute surface pressure gradient which was not accounted for in 'Um_dPHdx'
0034 % but is now included in 'Um_dPhiX'
0035 % b) compute Non-Hydrostatic pressure gradient (if Non-Hyd) which was not accounted
0036 % for in 'Um_dPHdx' but is now included in 'Um_dPhiX'
0037 % c) form vertical viscous flux 'VISrI_Um', to compute tendency from implicit vertical
0038 % viscosity (with explicit bottom drag) since diagnostic 'Um_ImplD' was not there
0039 G=load_grid(gDir,0,180,360);
0040 % fieldnames(G)
0041 nc=G.dims(2); nr=G.dims(3); nPxy=G.dims(1)*G.dims(2); nPp2=nPxy+2;
0042 ncx=6*nc; np1=nc+1; globArea=sum(G.rAc(:));
0043 %--
0044 yg2=zeros(nPp2,1); yg2([1:nPxy],1)=reshape(G.yG,[nPxy 1]);
0045 xg2=zeros(nPp2,1); xg2([1:nPxy],1)=reshape(G.xG,[nPxy 1]);
0046 rAz2=zeros(nPp2,1); rAz2([1:nPxy],1)=reshape(G.rAz,[nPxy 1]);
0047 %-- add missing corner:
0048 xg2(nPxy+1)=xg2(1); yg2(nPxy+1)=yg2(1+2*nc); rAz2(nPxy+1)=rAz2(1);
0049 xg2(nPxy+2)=xg2(1+3*nc); yg2(nPxy+2)=yg2(1); rAz2(nPxy+2)=rAz2(1);
0050
0051 mskW=min(ceil(G.hFacW),1); mskS=min(ceil(G.hFacS),1);
0052 %-------
0053 end
0054
0055 %-----------------------
0056 %- set constant: gravity, rhoConst
0057 rhoConst=1035.; gravity=9.81;
0058 %--------------------------------------------------
0059
0060 n3d=nr;
0061 if kpr > 0,
0062
0063 %-- Read in 2-D diagnostics file "namFs" and 3-D diagnostics file "namF":
0064 if isempty(namFs), namS=''; nV2d=0;
0065 else
0066 [v3d,iter,M]=rdmds([rDir,namFs],nit);
0067 eval(M) ; f2dList=fldList; namS=char(f2dList) ; nV2d=size(namS,1);
0068 end
0069 [v4d,iter,M]=rdmds([rDir,namF],nit);
0070 eval(M) ; namV=char(fldList) ; nV=size(namV,1);
0071 if nV2d == 0, f2dList=fldList; end
0072
0073 jdps=find(strcmp(f2dList,'PHI_SURF'));
0074 if isempty(jdps), jdps=0;
0075 jdps=find(strcmp(f2dList,'ETAN '));
0076 if isempty(jdps), jdps=0; else
0077 eta=v3d(:,:,jdps);
0078 %b0s=rdmds([gDir,'Bo_surf']); var=b0s.*eta;
0079 var=gravity*eta;
0080 end
0081 else var=v3d(:,:,jdps); end
0082 if jdps ~= 0,
0083 v6t=split_C_cub(var);
0084 var=v6t(2:np1,:,:)-v6t(1:nc,:,:); dpx=var(:,2:np1,:);
0085 dpx=reshape(permute(dpx,[1 3 2]),[ncx nc]);
0086 var=v6t(:,2:np1,:)-v6t(:,1:nc,:); dpy=var(2:np1,:,:);
0087 dpy=reshape(permute(dpy,[1 3 2]),[ncx nc]);
0088 dpx=-dpx./G.dXc; dpy=-dpy./G.dYc;
0089 end
0090
0091 fileName=sprintf('%s.%10.10i.%s',[rDir,'pnhDiag'],nit+1,'data');
0092 if isfile(fileName),
0093 fprintf(' -- loading file: %s ...',fileName);
0094 var=rdmds([rDir,'pnhDiag'],nit+1);
0095 fprintf(' done\n');
0096 else
0097 jnh=find(strcmp(fldList,'PHI_NH '));
0098 if isempty(jnh), jnh=0; else var=v4d(:,:,:,jnh); end
0099 end
0100 if jnh ~= 0,
0101 v6t=split_C_cub(var);
0102 var=v6t(2:np1,:,:,:)-v6t(1:nc,:,:,:); dpNHx=var(:,2:np1,:,:);
0103 dpNHx=reshape(permute(dpNHx,[1 4 2 3]),[ncx nc nr]);
0104 var=v6t(:,2:np1,:,:)-v6t(:,1:nc,:,:); dpNHy=var(2:np1,:,:,:);
0105 dpNHy=reshape(permute(dpNHy,[1 4 2 3]),[ncx nc nr]);
0106 var=repmat(G.dXc,[1 1 nr]);
0107 dpNHx=-dpNHx./var; dpNHx=dpNHx.*mskW;
0108 var=repmat(G.dYc,[1 1 nr]);
0109 dpNHy=-dpNHy./var; dpNHy=dpNHy.*mskS;
0110 end
0111
0112 jeta=find(strcmp(f2dList,'ETAN '));
0113 if isempty(jeta), jeta=0; else
0114 var=v3d(:,:,jeta).*G.rAc;
0115 v6t=split_C_cub(var);
0116 var=v6t(2:np1,:,:)+v6t(1:nc,:,:); vbx=0.5*var(:,2:np1,:);
0117 vbx=reshape(permute(vbx,[1 3 2]),[ncx nc]);
0118 var=v6t(:,2:np1,:)+v6t(:,1:nc,:); vby=0.5*var(2:np1,:,:);
0119 vby=reshape(permute(vby,[1 3 2]),[ncx nc]);
0120 vbx=vbx./G.rAw; vby=vby./G.rAs;
0121 var=G.depth;
0122 v6t=split_C_cub(var);
0123 var=min(v6t(2:np1,:,:),v6t(1:nc,:,:)); hhx=var(:,2:np1,:);
0124 hhx=reshape(permute(hhx,[1 3 2]),[ncx nc]);
0125 var=min(v6t(:,2:np1,:),v6t(:,1:nc,:)); hhy=var(2:np1,:,:);
0126 hhy=reshape(permute(hhy,[1 3 2]),[ncx nc]);
0127 %-- when using z* with older output, need to account for column vertical streaching
0128 % in computation of vertical viscosity tendency form vertical viscous flux 'VISrI_Um'
0129 rFacW=hhx; rFacW(find(hhx==0.))=-1; rFacW=vbx./rFacW;
0130 rFacW(find(hhx==0.))=0; rFacW=rFacW+ones(ncx,nc);
0131 rFacS=hhy; rFacS(find(hhy==0.))=-1; rFacS=vby./rFacS;
0132 rFacS(find(hhy==0.))=0; rFacS=rFacS+ones(ncx,nc);
0133 end
0134
0135 end
0136
0137 %-----------------------------------------------------------------------------------------
0138 fmtStats='Var= %s : Min,Max,Avr,StD= %12.5e %12.5e %12.5e %12.5e\n';
0139 %fmtSum0=' Sum Tend : Min,Max,Avr,StD= %12.5e %12.5e %12.5e %12.5e\n';
0140 fmtSum1=' Sum Tend: Avr,StD= %12.5e %12.5e ';
0141 fmtSum2='; Residual= %12.5e %12.5e +\n';
0142
0143 gUdp=zeros(ncx,nc,nr); gVdp=gUdp; titUdp=' ? '; titVdp=' ? ';
0144 j1=find(strcmp(fldList,'Um_dPhiX'));
0145 j2=find(strcmp(fldList,'Vm_dPhiY'));
0146 if isempty(j1) & isempty(j2), jdph=0;
0147 j1=find(strcmp(fldList,'Um_dPHdx'));
0148 j2=find(strcmp(fldList,'Vm_dPHdy'));
0149 if ~isempty(j1), jdph=j1; else
0150 if ~isempty(j2), jdph=j2; end; end
0151 if jdph ~= 0 & jdps ~= 0,
0152 gUdp=repmat(dpx,[1 1 nr]).*mskW;
0153 gVdp=repmat(dpy,[1 1 nr]).*mskS;
0154 end
0155 if jnh ~= 0, gUdp=gUdp+dpNHx; gVdp=gVdp+dpNHy; end
0156 else
0157 if isempty(j1), jdph=j2; else jdph=j1; end
0158 end
0159 if jdph ~= 0,
0160 if ~isempty(j1), gUdp=gUdp+squeeze(v4d(:,:,:,j1)); titUdp=char(fldList(j1)); end
0161 if ~isempty(j2), gVdp=gVdp+squeeze(v4d(:,:,:,j2)); titVdp=char(fldList(j2)); end
0162 if jdps ~= 0,
0163 titUdp=[titUdp(1:end-1),upper(titUdp(end))];
0164 titVdp=[titVdp(1:end-1),upper(titVdp(end))];
0165 %fprintf(' titUdp: >%s< ; titVdp: >%s<\n',titUdp,titVdp);
0166 end
0167 end
0168
0169 %-- Tendencies from implicit vertical viscous fluxes
0170 % Note: will be used to close momentum budget
0171 % a) if using older output (since 'Um_ImplD' was not there);
0172 % b) and using implicit viscosity but without implicit bottom friction
0173 % In the latest case (selectImplicitDrag=2,) cannot close the budget using older output
0174 j=find(strcmp(fldList,'VISrI_Um'));
0175 if isempty(j), juNz=0; else
0176 fprintf(' -- Tendencies from vertically visc. fluxes --\n');
0177 juNz=j; var=squeeze(v4d(:,:,:,j));
0178 %- compute tendency from minus div of vert. fluxes:
0179 div=var; div(:,:,1:nr-1)=div(:,:,1:nr-1)-var(:,:,2:nr);
0180 ddz=G.hFacW.*repmat(reshape(G.dRf,[1 1 nr]),[ncx nc 1]);
0181 rdz=ddz; rdz(find(ddz==0))=-1;
0182 rdz=1./rdz; rdz(find(ddz==0))=0;
0183 div=div.*rdz;
0184 div=div./repmat(G.rAw,[1 1 nr]);
0185 div=div./repmat(rFacW,[1 1 nr]);
0186 gUnuZ=-div;
0187 titv=char(fldList(j)); var=gUnuZ;
0188 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0189 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0190 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0191 %--
0192 j=find(strcmp(fldList,'Um_ImplD'));
0193 if ~isempty(j),
0194 titv=char(fldList(j)); var=squeeze(v4d(:,:,:,j));
0195 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0196 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0197 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0198 var=var-gUnuZ; titv='Diff:2-1';
0199 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0200 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0201 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0202 end
0203 fprintf('\n');
0204 end
0205
0206 j=find(strcmp(fldList,'VISrI_Vm'));
0207 if isempty(j), jvNz=0; else
0208 jvNz=j; var=squeeze(v4d(:,:,:,j));
0209 %- compute tendency from minus div of vert. fluxes:
0210 div=var; div(:,:,1:nr-1)=div(:,:,1:nr-1)-var(:,:,2:nr);
0211 ddz=G.hFacS.*repmat(reshape(G.dRf,[1 1 nr]),[ncx nc 1]);
0212 rdz=ddz; rdz(find(ddz==0))=-1;
0213 rdz=1./rdz; rdz(find(ddz==0))=0;
0214 div=div.*rdz;
0215 div=div./repmat(G.rAs,[1 1 nr]);
0216 div=div./repmat(rFacS,[1 1 nr]);
0217 gVnuZ=-div;
0218 titv=char(fldList(j)); var=gVnuZ;
0219 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0220 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0221 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0222 %--
0223 j=find(strcmp(fldList,'Vm_ImplD'));
0224 if ~isempty(j),
0225 titv=char(fldList(j)); var=squeeze(v4d(:,:,:,j));
0226 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0227 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0228 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0229 var=var-gVnuZ; titv='Diff:2-1';
0230 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0231 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0232 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0233 end
0234 fprintf('\n');
0235 end
0236 %juNz=0; jvNz=0;
0237
0238 %- Here we check that vertical integral of implicit vertical viscous tendency
0239 % match either bottom drag (if using implicit bottom drag) or simply zero.
0240 j1=find(strcmp(fldList,'Um_ImplD'));
0241 j2=find(strcmp(f2dList,'botTauX '));
0242 if ~isempty(j1) & ~isempty(j2),
0243 fprintf(' -- Vertically integrated tendencies --\n');
0244 j=j2;
0245 titv=char(f2dList(j)); var=v3d(:,:,j); bTauX=var;
0246 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0247 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0248 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0249 j=j1;
0250 titv=char(fldList(j1)); var=squeeze(v4d(:,:,:,j));
0251 %- vertically integrated:
0252 ddz=G.hFacW.*repmat(reshape(G.dRf,[1 1 nr]),[ncx nc 1]);
0253 var=var.*ddz; var=sum(var,3); var=var*rhoConst;
0254 var=var.*rFacW;
0255 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0256 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0257 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0258 var=var-bTauX; titv='Diff:2-1';
0259 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0260 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0261 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0262 fprintf('\n');
0263 end
0264
0265 j1=find(strcmp(fldList,'Vm_ImplD'));
0266 j2=find(strcmp(f2dList,'botTauY '));
0267 if ~isempty(j1) & ~isempty(j2),
0268 j=j2;
0269 titv=char(f2dList(j)); var=v3d(:,:,j); bTauY=var;
0270 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0271 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0272 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0273 j=j1;
0274 titv=char(fldList(j1)); var=squeeze(v4d(:,:,:,j));
0275 %- vertically integrated:
0276 ddz=G.hFacS.*repmat(reshape(G.dRf,[1 1 nr]),[ncx nc 1]);
0277 var=var.*ddz; var=sum(var,3); var=var*rhoConst;
0278 var=var.*rFacS;
0279 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0280 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0281 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0282 var=var-bTauY; titv='Diff:2-1';
0283 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0284 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0285 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0286 fprintf('\n');
0287 end
0288
0289 %fprintf('\n');
0290 %------
0291 fprintf(' -- Check Mom budget, exp: %s, files: %s & %s, it= %i\n',rDir,namF,namFs,nit);
0292
0293 j=find(strcmp(fldList,'TOTUTEND'));
0294 if ~isempty(j),
0295 titv=char(fldList(j)); var=squeeze(v4d(:,:,:,j))/86400;
0296 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0297 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0298 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0299 dUtot=var;
0300 end
0301
0302 %- For each term: a) print some stats of this term ;
0303 % b) add to other tendency and print stats of the sum
0304 % c) substract the sum from total tendency (-> residual) and print stats
0305
0306 j=find(strcmp(fldList,'Um_dPhiX'));
0307 if isempty(j), j=find(strcmp(fldList,'Um_dPHdx')); end
0308 if ~isempty(j),
0309 titv=titUdp; var=gUdp;
0310 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0311 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0312 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0313 gUtot=var;
0314 end
0315
0316 j=find(strcmp(fldList,'Um_Advec'));
0317 if ~isempty(j),
0318 titv=char(fldList(j)); var=squeeze(v4d(:,:,:,j));
0319 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0320 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0321 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0322 gUtot=gUtot+var; var=gUtot;
0323 Avr=mean(var(:)); vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0324 fprintf(fmtSum1,Avr,StD); var=dUtot-var;
0325 Avr=mean(var(:)); vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0326 fprintf(fmtSum2,Avr,StD);
0327 end
0328
0329 j=find(strcmp(fldList,'Um_Ext '));
0330 if ~isempty(j),
0331 titv=char(fldList(j)); var=squeeze(v4d(:,:,:,j));
0332 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0333 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0334 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0335 gUtot=gUtot+var; var=gUtot;
0336 Avr=mean(var(:)); vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0337 fprintf(fmtSum1,Avr,StD); var=dUtot-var;
0338 Avr=mean(var(:)); vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0339 fprintf(fmtSum2,Avr,StD);
0340 end
0341
0342 j=find(strcmp(fldList,'Um_Diss '));
0343 if ~isempty(j), jDexp=j;
0344 titv=char(fldList(j)); var=squeeze(v4d(:,:,:,j));
0345 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0346 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0347 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0348 gUtot=gUtot+var; var=gUtot;
0349 Avr=mean(var(:)); vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0350 fprintf(fmtSum1,Avr,StD); var=dUtot-var;
0351 Avr=mean(var(:)); vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0352 fprintf(fmtSum2,Avr,StD);
0353 end
0354 j=find(strcmp(fldList,'Um_ImplD'));
0355 if ~isempty(j),
0356 var=squeeze(v4d(:,:,:,j));
0357 elseif juNz ~=0,
0358 j=juNz; var=gUnuZ;
0359 else j=0; end
0360 if j ~= 0,
0361 titv=char(fldList(j));
0362 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0363 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0364 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0365 gUtot=gUtot+var; var=gUtot;
0366 Avr=mean(var(:)); vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0367 fprintf(fmtSum1,Avr,StD); var=dUtot-var;
0368 Avr=mean(var(:)); vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0369 fprintf(fmtSum2,Avr,StD);
0370 end
0371
0372 j=find(strcmp(fldList,'AB_gU '));
0373 if ~isempty(j),
0374 titv=char(fldList(j)); var=squeeze(v4d(:,:,:,j));
0375 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0376 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0377 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0378 gUtot=gUtot+var; var=gUtot;
0379 Avr=mean(var(:)); vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0380 fprintf(fmtSum1,Avr,StD); var=dUtot-var;
0381 Avr=mean(var(:)); vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0382 fprintf(fmtSum2,Avr,StD);
0383 end
0384
0385 fprintf('\n');
0386
0387 j=find(strcmp(fldList,'TOTVTEND'));
0388 if ~isempty(j),
0389 titv=char(fldList(j)); var=squeeze(v4d(:,:,:,j))/86400;
0390 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0391 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0392 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0393 dVtot=var;
0394 end
0395
0396 j=find(strcmp(fldList,'Vm_dPhiY'));
0397 if isempty(j), j=find(strcmp(fldList,'Vm_dPHdy')); end
0398 if ~isempty(j),
0399 %titv=char(fldList(j)); var=squeeze(v4d(:,:,:,j));
0400 titv=titVdp; var=gVdp;
0401 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0402 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0403 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0404 gVtot=var;
0405 end
0406
0407 j=find(strcmp(fldList,'Vm_Advec'));
0408 if ~isempty(j),
0409 titv=char(fldList(j)); var=squeeze(v4d(:,:,:,j));
0410 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0411 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0412 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0413 gVtot=gVtot+var; var=gVtot;
0414 Avr=mean(var(:)); vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0415 fprintf(fmtSum1,Avr,StD); var=dVtot-var;
0416 Avr=mean(var(:)); vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0417 fprintf(fmtSum2,Avr,StD);
0418 end
0419
0420 j=find(strcmp(fldList,'Vm_Ext '));
0421 if ~isempty(j),
0422 titv=char(fldList(j)); var=squeeze(v4d(:,:,:,j));
0423 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0424 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0425 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0426 gVtot=gVtot+var; var=gVtot;
0427 Avr=mean(var(:)); vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0428 fprintf(fmtSum1,Avr,StD); var=dVtot-var;
0429 Avr=mean(var(:)); vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0430 fprintf(fmtSum2,Avr,StD);
0431 end
0432
0433 j=find(strcmp(fldList,'Vm_Diss '));
0434 if ~isempty(j),
0435 titv=char(fldList(j)); var=squeeze(v4d(:,:,:,j));
0436 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0437 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0438 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0439 gVtot=gVtot+var; var=gVtot;
0440 Avr=mean(var(:)); vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0441 fprintf(fmtSum1,Avr,StD); var=dVtot-var;
0442 Avr=mean(var(:)); vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0443 fprintf(fmtSum2,Avr,StD);
0444 end
0445 j=find(strcmp(fldList,'Vm_ImplD'));
0446 if ~isempty(j),
0447 var=squeeze(v4d(:,:,:,j));
0448 elseif jvNz ~=0,
0449 j=jvNz; var=gVnuZ;
0450 else j=0; end
0451 if j ~= 0,
0452 titv=char(fldList(j));
0453 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0454 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0455 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0456 gVtot=gVtot+var; var=gVtot;
0457 Avr=mean(var(:)); vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0458 fprintf(fmtSum1,Avr,StD); var=dVtot-var;
0459 Avr=mean(var(:)); vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0460 fprintf(fmtSum2,Avr,StD);
0461 end
0462
0463 j=find(strcmp(fldList,'AB_gV '));
0464 if ~isempty(j),
0465 titv=char(fldList(j)); var=squeeze(v4d(:,:,:,j));
0466 mnV=min(var(:)); MxV=max(var(:)); Avr=mean(var(:));
0467 vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0468 fprintf(fmtStats,titv,mnV,MxV,Avr,StD);
0469 gVtot=gVtot+var; var=gVtot;
0470 Avr=mean(var(:)); vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0471 fprintf(fmtSum1,Avr,StD); var=dVtot-var;
0472 Avr=mean(var(:)); vv2=var.*var; StD=mean(vv2(:))-Avr*Avr; StD=sqrt(StD);
0473 fprintf(fmtSum2,Avr,StD);
0474 end
0475
0476 return
