| 1 |
% $Header: $ |
| 2 |
% $Name: $ |
| 3 |
|
| 4 |
%-- after running mitgcmuv, in ../tr_run.thsice: |
| 5 |
% 1) extract numerical values from ASCII stat-Diag output files: |
| 6 |
% ( script "extract_StD" & matlab script read_StD.m |
| 7 |
% are located in in MITgcm_contrib/jmc_script ) |
| 8 |
% > extract_StD instStDiag.0000036000.txt = txt |
| 9 |
% > extract_StD surfStDiag.0000036000.txt = txt |
| 10 |
% 2) start matlab and execute this script. |
| 11 |
|
| 12 |
rac='../tr_run.thsice/'; nam='check-Conserv'; dBug=0; |
| 13 |
nit0=36000; nit=nit0+[1:9:20]; |
| 14 |
i1=1; i2=3; % <- to check half of this period: (i1,i2)=(1,2) or (i1,i2)=(2,3) |
| 15 |
|
| 16 |
arc=rdmds([rac,'RAC']); |
| 17 |
drF=squeeze(rdmds([rac,'DRF'])); delR=drF'; |
| 18 |
[ncx,nc]=size(arc); nPg=ncx*nc; nr=size(drF,1); |
| 19 |
cpwater=3994; Lfresh=334000; |
| 20 |
rhosw=1035; rhofw=1000; |
| 21 |
rhoIc=900; rhos=330; |
| 22 |
salIce=4; salSW=35; |
| 23 |
epsAB=0.6; deltaT=86400; |
| 24 |
|
| 25 |
hFac=rdmds([rac,'hFacC']); |
| 26 |
msk1=squeeze(hFac(:,:,1)); msk1=ceil(msk1); |
| 27 |
|
| 28 |
ocArea=sum(sum(msk1.*arc)); |
| 29 |
Depth=rdmds([rac,'Depth']); |
| 30 |
var=msk1.*arc; var=var.*Depth; DpAv1=sum(var(:))/ocArea; |
| 31 |
rDepth=Depth+msk1-1; rDepth=msk1./rDepth; |
| 32 |
|
| 33 |
var=reshape(arc,nPg,1); |
| 34 |
vol=var*delR; var=reshape(vol,ncx,nc,nr); |
| 35 |
vol=hFac.*var; Volum=sum(sum(sum(vol))); DpAv2=Volum/ocArea; |
| 36 |
fprintf(['Exp: ',nam,' (dir=',rac(1:end-1),'),']); |
| 37 |
fprintf(' Check-Cons between iter: %i %i\n',nit(i1),nit(i2)); |
| 38 |
fprintf('Ocean Area= %g ; Volume= %g mean Depth= %11.6f \n',ocArea,Volum,DpAv1); |
| 39 |
fprintf(' mean Depth2= %11.6f (Diff: %10.3e)\n',DpAv2,DpAv2-DpAv1); |
| 40 |
|
| 41 |
%-------------------- |
| 42 |
[v2i,iti,M]=rdmds([rac,'surfInst'],nit-1); |
| 43 |
if dBug, |
| 44 |
eval(M) |
| 45 |
v4d=reshape(v2i,[dimList nFlds length(iti)]); |
| 46 |
for i=1:length(iti), |
| 47 |
for n=1:nFlds, |
| 48 |
if nDims == 2, var=v4d(:,:,n,i); else var=v4d(:,:,:,n,i); end |
| 49 |
fprintf(['it= %i %2i ',char(fldList(n)),': min,max = %e , %e\n'], ... |
| 50 |
iti(i),n,min(var(:)),max(var(:))); |
| 51 |
end |
| 52 |
end |
| 53 |
end |
| 54 |
[v2c,iti,M]=rdmds([rac,'etaInst'],nit); |
| 55 |
if dBug, |
| 56 |
eval(M) |
| 57 |
v4d=reshape(v2c,[dimList nFlds length(iti)]); |
| 58 |
for i=1:length(iti), |
| 59 |
for n=1:nFlds, |
| 60 |
if nDims == 2, var=v4d(:,:,n,i); else var=v4d(:,:,:,n,i); end |
| 61 |
fprintf(['it= %i %2i ',char(fldList(n)),': min,max = %e , %e\n'], ... |
| 62 |
iti(i),n,min(var(:)),max(var(:))); |
| 63 |
end |
| 64 |
end |
| 65 |
end |
| 66 |
|
| 67 |
[v3i,iti,M]=rdmds([rac,'dynInst'],nit-1); |
| 68 |
if dBug, |
| 69 |
eval(M) |
| 70 |
v4d=reshape(v3i,[dimList nFlds length(iti)]); |
| 71 |
for i=1:length(iti), |
| 72 |
for n=1:nFlds, |
| 73 |
if nDims == 2, var=v4d(:,:,n,i); else var=v4d(:,:,:,n,i); end |
| 74 |
fprintf(['it= %i %2i ',char(fldList(n)),': min,max = %e , %e\n'], ... |
| 75 |
iti(i),n,min(var(:)),max(var(:))); |
| 76 |
end |
| 77 |
end |
| 78 |
end |
| 79 |
|
| 80 |
%- vvA: kLev, time_rec, region_rec, [ave,std,min,max,vol], var_rec |
| 81 |
namF=[rac,'surfStDiag']; |
| 82 |
[nIt,rList,tim1,vvA,listV1]=read_StD(namF,'txt','all_flds'); |
| 83 |
flxAve=squeeze(vvA); |
| 84 |
namF=[rac,'instStDiag']; |
| 85 |
[nIt,rList,tim2,vvA,listV2]=read_StD(namF,'txt','all_flds'); |
| 86 |
dynAve=squeeze(vvA); |
| 87 |
|
| 88 |
delT=(nit(2)-nit(1))*deltaT; |
| 89 |
|
| 90 |
et1=v2i(:,:,1,i1); et2=v2i(:,:,1,i2); |
| 91 |
pr1=v2i(:,:,2,i1); pr2=v2i(:,:,2,i2); |
| 92 |
et1p=v2c(:,:,i1); et2p=v2c(:,:,i2); |
| 93 |
|
| 94 |
t1=v3i(:,:,:,1,i1); t2=v3i(:,:,:,1,i2); |
| 95 |
s1=v3i(:,:,:,2,i1); s2=v3i(:,:,:,2,i2); |
| 96 |
|
| 97 |
%--- to compute T* & S* (= conserved quantities with AB-2) at it1 & it2 : |
| 98 |
dh1=et1p-et1; |
| 99 |
dv1=reshape(dh1.*rDepth,[nPg 1])*ones(1,nr); dv1=reshape(dv1,[ncx nc nr]); |
| 100 |
dv1(:,:,1)=dv1(:,:,1)-pr1/rhosw*deltaT/delR(1); |
| 101 |
dh2=et2p-et2; |
| 102 |
dv2=reshape(dh2.*rDepth,[nPg 1])*ones(1,nr); dv2=reshape(dv2,[ncx nc nr]); |
| 103 |
dv2(:,:,1)=dv2(:,:,1)-pr2/rhosw*deltaT/delR(1); |
| 104 |
%--- |
| 105 |
|
| 106 |
H1=dynAve(1,i1,1,2)*dynAve(1,i1,5,2)/ocArea*rhosw*cpwater; |
| 107 |
H2=dynAve(1,i2,1,2)*dynAve(1,i2,5,2)/ocArea*rhosw*cpwater; |
| 108 |
var=epsAB*t1.*dv1; |
| 109 |
H3=sum(sum(sum(var.*vol)))/ocArea*rhosw*cpwater; |
| 110 |
var=epsAB*t2.*dv2; |
| 111 |
H4=sum(sum(sum(var.*vol)))/ocArea*rhosw*cpwater; |
| 112 |
%heating=flxAve(:,1,3)+flxAve(:,1,9); % = tRelax + Qnet (EnPrec missing) |
| 113 |
heating=flxAve(:,1,7); % = TFLUX |
| 114 |
heating=sum(heating(1+i1:i2))*delT; |
| 115 |
fprintf(' H2-H1 = %8.3e ; heating= %8.3e ; Diff: %11.4e (%11.4e)\n', ... |
| 116 |
H2-H1,heating,H2-H1-heating,(H2-H1-heating)/H2); |
| 117 |
fprintf(' h2-h1 = %10.3e (h1=%8.1e, h2=%8.1e): %11.4e (%11.4e)\n', ... |
| 118 |
H4-H3,H3,H4,H2+H4-H1-H3-heating,(H2+H4-H1-H3-heating)/H2); |
| 119 |
|
| 120 |
hS1=dynAve(1,i1,1,3)*dynAve(1,i1,5,3)/ocArea*rhosw; |
| 121 |
hS2=dynAve(1,i2,1,3)*dynAve(1,i2,5,3)/ocArea*rhosw; |
| 122 |
var=epsAB*s1.*dv1; |
| 123 |
hS3=sum(sum(sum(var.*vol)))/ocArea*rhosw; |
| 124 |
var=epsAB*s2.*dv2; |
| 125 |
hS4=sum(sum(sum(var.*vol)))/ocArea*rhosw; |
| 126 |
saltfx=flxAve(:,1,4)+flxAve(:,1,10); % = sRelax + oceSflux |
| 127 |
%saltfx=flxAve(:,1,8); % = SFLUX |
| 128 |
saltfx=sum(saltfx(1+i1:i2))*delT; |
| 129 |
fprintf(' S2-S1 = %8.6f ; saltfx = %8.6f ; Diff: %11.4e (%11.4e)\n', ... |
| 130 |
hS2-hS1,saltfx,hS2-hS1-saltfx,(hS2-hS1-saltfx)/hS2); |
| 131 |
fprintf(' s2-s1 = %10.3e (s1=%8.1e, s2=%8.1e): %11.4e (%11.4e)\n', ... |
| 132 |
hS4-hS3,hS3,hS4,hS2+hS4-hS1-hS3-saltfx,(hS2+hS4-hS1-hS3-saltfx)/hS2); |
| 133 |
|
| 134 |
etAv1=dynAve(1,i1,1,1)*rhosw; |
| 135 |
etAv2=dynAve(1,i2,1,1)*rhosw; |
| 136 |
%- more precised when recomputing average from from full precision 2D field: |
| 137 |
%var=et1p; var=msk1.*var; etAv1=sum(sum(var.*arc))/ocArea*rhosw; |
| 138 |
%var=et2p; var=msk1.*var; etAv2=sum(sum(var.*arc))/ocArea*rhosw; |
| 139 |
fresh=flxAve(:,1,11); fresh=sum(fresh(1+i1:i2))*delT; |
| 140 |
fprintf(' Et2-Et1 = %8.6f ; fresh= %8.6f ; Diff: %11.4e (%11.4e)\n', ... |
| 141 |
etAv2-etAv1,fresh,etAv2-etAv1-fresh,(etAv2-etAv1-fresh)/etAv2); |
| 142 |
|
| 143 |
return |
Parent Directory
| 
Revision Log
| 
Revision Graph