Warning, /verification/global_ocean.cs32x15/input.thsice/check_OceConserve.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
7e42d13321 Jean*0001 %-- after running mitgcmuv, in ../tr_run.thsice:
0002 % 1) extract numerical values from ASCII stat-Diag output files:
0003 % ( script "extract_StD" & matlab script read_StD.m
0004 % are located in in MITgcm_contrib/jmc_script )
0005 % > extract_StD instStDiag.0000036000.txt = txt
0006 % > extract_StD surfStDiag.0000036000.txt = txt
0007 % 2) start matlab and execute this script.
0008
0009 rac='../tr_run.thsice/'; nam='check-Conserv'; dBug=0;
0010 nit0=36000; nit=nit0+[1:9:20];
0011 i1=1; i2=3; % <- to check half of this period: (i1,i2)=(1,2) or (i1,i2)=(2,3)
0012
0013 arc=rdmds([rac,'RAC']);
0014 drF=squeeze(rdmds([rac,'DRF'])); delR=drF';
0015 [ncx,nc]=size(arc); nPg=ncx*nc; nr=size(drF,1);
0016 cpwater=3994; Lfresh=334000;
0017 rhosw=1035; rhofw=1000;
0018 rhoIc=900; rhos=330;
0019 salIce=4; salSW=35;
0020 epsAB=0.6; deltaT=86400;
0021
0022 hFac=rdmds([rac,'hFacC']);
0023 msk1=squeeze(hFac(:,:,1)); msk1=ceil(msk1);
0024
0025 ocArea=sum(sum(msk1.*arc));
0026 Depth=rdmds([rac,'Depth']);
0027 var=msk1.*arc; var=var.*Depth; DpAv1=sum(var(:))/ocArea;
0028 rDepth=Depth+msk1-1; rDepth=msk1./rDepth;
0029
0030 var=reshape(arc,nPg,1);
0031 vol=var*delR; var=reshape(vol,ncx,nc,nr);
0032 vol=hFac.*var; Volum=sum(sum(sum(vol))); DpAv2=Volum/ocArea;
0033 fprintf(['Exp: ',nam,' (dir=',rac(1:end-1),'),']);
0034 fprintf(' Check-Cons between iter: %i %i\n',nit(i1),nit(i2));
0035 fprintf('Ocean Area= %g ; Volume= %g mean Depth= %11.6f \n',ocArea,Volum,DpAv1);
0036 fprintf(' mean Depth2= %11.6f (Diff: %10.3e)\n',DpAv2,DpAv2-DpAv1);
0037
0038 %--------------------
0039 [v2i,iti,M]=rdmds([rac,'surfInst'],nit-1);
0040 if dBug,
0041 eval(M)
0042 v4d=reshape(v2i,[dimList nFlds length(iti)]);
0043 for i=1:length(iti),
0044 for n=1:nFlds,
0045 if nDims == 2, var=v4d(:,:,n,i); else var=v4d(:,:,:,n,i); end
0046 fprintf(['it= %i %2i ',char(fldList(n)),': min,max = %e , %e\n'], ...
0047 iti(i),n,min(var(:)),max(var(:)));
0048 end
0049 end
0050 end
0051 [v2c,iti,M]=rdmds([rac,'etaInst'],nit);
0052 if dBug,
0053 eval(M)
0054 v4d=reshape(v2c,[dimList nFlds length(iti)]);
0055 for i=1:length(iti),
0056 for n=1:nFlds,
0057 if nDims == 2, var=v4d(:,:,n,i); else var=v4d(:,:,:,n,i); end
0058 fprintf(['it= %i %2i ',char(fldList(n)),': min,max = %e , %e\n'], ...
0059 iti(i),n,min(var(:)),max(var(:)));
0060 end
0061 end
0062 end
0063
0064 [v3i,iti,M]=rdmds([rac,'dynInst'],nit-1);
0065 if dBug,
0066 eval(M)
0067 v4d=reshape(v3i,[dimList nFlds length(iti)]);
0068 for i=1:length(iti),
0069 for n=1:nFlds,
0070 if nDims == 2, var=v4d(:,:,n,i); else var=v4d(:,:,:,n,i); end
0071 fprintf(['it= %i %2i ',char(fldList(n)),': min,max = %e , %e\n'], ...
0072 iti(i),n,min(var(:)),max(var(:)));
0073 end
0074 end
0075 end
0076
0077 %- vvA: kLev, time_rec, region_rec, [ave,std,min,max,vol], var_rec
0078 namF=[rac,'surfStDiag'];
0079 [nIt,rList,tim1,vvA,listV1]=read_StD(namF,'txt','all_flds');
0080 flxAve=squeeze(vvA);
0081 namF=[rac,'instStDiag'];
0082 [nIt,rList,tim2,vvA,listV2]=read_StD(namF,'txt','all_flds');
0083 dynAve=squeeze(vvA);
0084
0085 delT=(nit(2)-nit(1))*deltaT;
0086
0087 et1=v2i(:,:,1,i1); et2=v2i(:,:,1,i2);
0088 pr1=v2i(:,:,2,i1); pr2=v2i(:,:,2,i2);
0089 et1p=v2c(:,:,i1); et2p=v2c(:,:,i2);
0090
0091 t1=v3i(:,:,:,1,i1); t2=v3i(:,:,:,1,i2);
0092 s1=v3i(:,:,:,2,i1); s2=v3i(:,:,:,2,i2);
0093
0094 %--- to compute T* & S* (= conserved quantities with AB-2) at it1 & it2 :
0095 dh1=et1p-et1;
0096 dv1=reshape(dh1.*rDepth,[nPg 1])*ones(1,nr); dv1=reshape(dv1,[ncx nc nr]);
0097 dv1(:,:,1)=dv1(:,:,1)-pr1/rhosw*deltaT/delR(1);
0098 dh2=et2p-et2;
0099 dv2=reshape(dh2.*rDepth,[nPg 1])*ones(1,nr); dv2=reshape(dv2,[ncx nc nr]);
0100 dv2(:,:,1)=dv2(:,:,1)-pr2/rhosw*deltaT/delR(1);
0101 %---
0102
0103 H1=dynAve(1,i1,1,2)*dynAve(1,i1,5,2)/ocArea*rhosw*cpwater;
0104 H2=dynAve(1,i2,1,2)*dynAve(1,i2,5,2)/ocArea*rhosw*cpwater;
0105 var=epsAB*t1.*dv1;
0106 H3=sum(sum(sum(var.*vol)))/ocArea*rhosw*cpwater;
0107 var=epsAB*t2.*dv2;
0108 H4=sum(sum(sum(var.*vol)))/ocArea*rhosw*cpwater;
0109 %heating=flxAve(:,1,3)+flxAve(:,1,9); % = tRelax + Qnet (EnPrec missing)
0110 heating=flxAve(:,1,7); % = TFLUX
0111 heating=sum(heating(1+i1:i2))*delT;
0112 fprintf(' H2-H1 = %8.3e ; heating= %8.3e ; Diff: %11.4e (%11.4e)\n', ...
0113 H2-H1,heating,H2-H1-heating,(H2-H1-heating)/H2);
0114 fprintf(' h2-h1 = %10.3e (h1=%8.1e, h2=%8.1e): %11.4e (%11.4e)\n', ...
0115 H4-H3,H3,H4,H2+H4-H1-H3-heating,(H2+H4-H1-H3-heating)/H2);
0116
0117 hS1=dynAve(1,i1,1,3)*dynAve(1,i1,5,3)/ocArea*rhosw;
0118 hS2=dynAve(1,i2,1,3)*dynAve(1,i2,5,3)/ocArea*rhosw;
0119 var=epsAB*s1.*dv1;
0120 hS3=sum(sum(sum(var.*vol)))/ocArea*rhosw;
0121 var=epsAB*s2.*dv2;
0122 hS4=sum(sum(sum(var.*vol)))/ocArea*rhosw;
0123 saltfx=flxAve(:,1,4)+flxAve(:,1,10); % = sRelax + oceSflux
0124 %saltfx=flxAve(:,1,8); % = SFLUX
0125 saltfx=sum(saltfx(1+i1:i2))*delT;
0126 fprintf(' S2-S1 = %8.6f ; saltfx = %8.6f ; Diff: %11.4e (%11.4e)\n', ...
0127 hS2-hS1,saltfx,hS2-hS1-saltfx,(hS2-hS1-saltfx)/hS2);
0128 fprintf(' s2-s1 = %10.3e (s1=%8.1e, s2=%8.1e): %11.4e (%11.4e)\n', ...
0129 hS4-hS3,hS3,hS4,hS2+hS4-hS1-hS3-saltfx,(hS2+hS4-hS1-hS3-saltfx)/hS2);
0130
0131 etAv1=dynAve(1,i1,1,1)*rhosw;
0132 etAv2=dynAve(1,i2,1,1)*rhosw;
0133 %- more precised when recomputing average from from full precision 2D field:
0134 %var=et1p; var=msk1.*var; etAv1=sum(sum(var.*arc))/ocArea*rhosw;
0135 %var=et2p; var=msk1.*var; etAv2=sum(sum(var.*arc))/ocArea*rhosw;
0136 fresh=flxAve(:,1,11); fresh=sum(fresh(1+i1:i2))*delT;
0137 fprintf(' Et2-Et1 = %8.6f ; fresh= %8.6f ; Diff: %11.4e (%11.4e)\n', ...
0138 etAv2-etAv1,fresh,etAv2-etAv1-fresh,(etAv2-etAv1-fresh)/etAv2);
0139
0140 return
