matlab_class/gcmfaces_diags/diags_set_MLD.m



if userStep==1;%diags to be computed
    listDiags='fldMldBoyer fldMldSuga fldMldKara';
elseif userStep==2;%input files and variables
    listFlds={    'THETA','SALT'};
    listFldsNames=deblank(listFlds);
    listFiles={'monthly_2d_set1','monthly_3d_set1','state_2d_set1','other_2d_set1','state_3d_set1'};
elseif userStep==3;%computational part;
        fldT=THETA.*mygrid.mskC; fldS=SALT.*mygrid.mskC;
        %
        %prepare to compute potential density:
        fldP=0*mygrid.mskC; for kk=1:length(mygrid.RC); fldP(:,:,kk)=-mygrid.RC(kk); end;
        T=convert2vector(fldT);
        S=convert2vector(fldS);
        msk=convert2vector(mygrid.mskC);
        P=convert2vector(fldP);
        %compute potential density:
        RHO=0*msk; alpha=0*msk;
        tmp1=find(~isnan(msk));
        RHO(tmp1) = density(T(tmp1),S(tmp1),P(tmp1));
        fldRhoPot=convert2vector(RHO);
        alpha(tmp1) = density(T(tmp1)+1e-4,S(tmp1),P(tmp1));
        fldAlpha=(convert2vector(alpha)-fldRhoPot)/1e-4;

        clear T S P msk RHO RHOis tmp1;

        %compute mld:
        tmp1=NaN*mygrid.mskC(:,:,1);
        for kk=1:50;
          tmp2=fldRhoPot(:,:,kk)-fldRhoPot(:,:,1);
          %if we pass RHO(1)+0.03 for the first time (or we reach the bottom)
          %then mld is the velocity point above RC(kk), which is RF(kk)
          jj=find((tmp2>0.03|isnan(tmp2))&isnan(tmp1));
         tmp1(jj)=-mygrid.RF(kk);
        end;
        fldMldBoyer=tmp1;

        %compute mld:
        tmp1=NaN*mygrid.mskC(:,:,1);
        for kk=1:50;
          tmp2=fldRhoPot(:,:,kk)-fldRhoPot(:,:,1);
          %if we pass RHO(1)+0.125 for the first time (or we reach the bottom)
          %then mld is the velocity point above RC(kk), which is RF(kk)
          jj=find((tmp2>0.125|isnan(tmp2))&isnan(tmp1));
         tmp1(jj)=-mygrid.RF(kk);
        end;
        fldMldSuga=tmp1;

        %compute mld:
        tmp1=NaN*mygrid.mskC(:,:,1);
        fldRhoPotMax=fldRhoPot(:,:,1)-0.8*fldAlpha(:,:,1);
        for kk=1:50;
          tmp2=fldRhoPot(:,:,kk)-fldRhoPotMax;
          %if we pass RHO(1)+0.8*alpha(1) for the first time (or we reach the bottom)
          %then mld is the velocity point above RC(kk), which is RF(kk)
          jj=find((tmp2>0|isnan(tmp2))&isnan(tmp1));
         tmp1(jj)=-mygrid.RF(kk);
        end;
        fldMldKara=tmp1;

elseif userStep==-1;%computational part;
                
        fld=mean(alldiag.fldMldKara(:,:,tt),3);
        cc=[[0:20:100] [150:50:300] 400 [500:200:1100] [1500:500:2000]]; title0='Mixed Layer Depth (Kara)';
        if doAnomalies; cc=scaleAnom*[-5:0.5:5]; end;
        figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
        myCaption={myYmeanTxt,'mean -- Kara mld (m)'};
        if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

        fld=mean(alldiag.fldMldSuga(:,:,tt),3);
        cc=[[0:20:100] [150:50:300] 400 [500:200:1100] [1500:500:2000]]; title0='Mixed Layer Depth (Suga)';
        if doAnomalies; cc=scaleAnom*[-5:0.5:5]; end;
        figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
        myCaption={myYmeanTxt,'mean -- Suga mld (m)'};
        if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

        fld=mean(alldiag.fldMldBoyer(:,:,tt),3);
        cc=[[0:20:100] [150:50:300] 400 [500:200:1100] [1500:500:2000]]; title0='Mixed Layer Depth (Boyer)';
        if doAnomalies; cc=scaleAnom*[-5:0.5:5]; end;
        figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
        myCaption={myYmeanTxt,'mean -- Boyer mld (m)'};
        if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

end;


1	gforget	1.1
2
3			if userStep==1;%diags to be computed
4			listDiags='fldMldBoyer fldMldSuga fldMldKara';
5			elseif userStep==2;%input files and variables
6			listFlds={ 'THETA','SALT'};
7			listFldsNames=deblank(listFlds);
8			listFiles={'monthly_2d_set1','monthly_3d_set1','state_2d_set1','other_2d_set1','state_3d_set1'};
9			elseif userStep==3;%computational part;
10			fldT=THETA.mygrid.mskC; fldS=SALT.mygrid.mskC;
11			%
12			%prepare to compute potential density:
13			fldP=0*mygrid.mskC; for kk=1:length(mygrid.RC); fldP(:,:,kk)=-mygrid.RC(kk); end;
14			T=convert2vector(fldT);
15			S=convert2vector(fldS);
16			msk=convert2vector(mygrid.mskC);
17			P=convert2vector(fldP);
18			%compute potential density:
19			RHO=0msk; alpha=0msk;
20			tmp1=find(~isnan(msk));
21			RHO(tmp1) = density(T(tmp1),S(tmp1),P(tmp1));
22			fldRhoPot=convert2vector(RHO);
23			alpha(tmp1) = density(T(tmp1)+1e-4,S(tmp1),P(tmp1));
24			fldAlpha=(convert2vector(alpha)-fldRhoPot)/1e-4;
25
26			clear T S P msk RHO RHOis tmp1;
27
28			%compute mld:
29			tmp1=NaN*mygrid.mskC(:,:,1);
30			for kk=1:50;
31			tmp2=fldRhoPot(:,:,kk)-fldRhoPot(:,:,1);
32			%if we pass RHO(1)+0.03 for the first time (or we reach the bottom)
33			%then mld is the velocity point above RC(kk), which is RF(kk)
34			jj=find((tmp2>0.03\|isnan(tmp2))&isnan(tmp1));
35			tmp1(jj)=-mygrid.RF(kk);
36			end;
37			fldMldBoyer=tmp1;
38
39			%compute mld:
40			tmp1=NaN*mygrid.mskC(:,:,1);
41			for kk=1:50;
42			tmp2=fldRhoPot(:,:,kk)-fldRhoPot(:,:,1);
43			%if we pass RHO(1)+0.125 for the first time (or we reach the bottom)
44			%then mld is the velocity point above RC(kk), which is RF(kk)
45			jj=find((tmp2>0.125\|isnan(tmp2))&isnan(tmp1));
46			tmp1(jj)=-mygrid.RF(kk);
47			end;
48			fldMldSuga=tmp1;
49
50			%compute mld:
51			tmp1=NaN*mygrid.mskC(:,:,1);
52			fldRhoPotMax=fldRhoPot(:,:,1)-0.8*fldAlpha(:,:,1);
53			for kk=1:50;
54			tmp2=fldRhoPot(:,:,kk)-fldRhoPotMax;
55			%if we pass RHO(1)+0.8*alpha(1) for the first time (or we reach the bottom)
56			%then mld is the velocity point above RC(kk), which is RF(kk)
57			jj=find((tmp2>0\|isnan(tmp2))&isnan(tmp1));
58			tmp1(jj)=-mygrid.RF(kk);
59			end;
60			fldMldKara=tmp1;
61
62			elseif userStep==-1;%computational part;
63
64			fld=mean(alldiag.fldMldKara(:,:,tt),3);
65			cc=[[0:20:100] [150:50:300] 400 [500:200:1100] [1500:500:2000]]; title0='Mixed Layer Depth (Kara)';
66			if doAnomalies; cc=scaleAnom*[-5:0.5:5]; end;
67			figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
68			myCaption={myYmeanTxt,'mean -- Kara mld (m)'};
69			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
70
71			fld=mean(alldiag.fldMldSuga(:,:,tt),3);
72			cc=[[0:20:100] [150:50:300] 400 [500:200:1100] [1500:500:2000]]; title0='Mixed Layer Depth (Suga)';
73			if doAnomalies; cc=scaleAnom*[-5:0.5:5]; end;
74			figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
75			myCaption={myYmeanTxt,'mean -- Suga mld (m)'};
76			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
77
78			fld=mean(alldiag.fldMldBoyer(:,:,tt),3);
79			cc=[[0:20:100] [150:50:300] 400 [500:200:1100] [1500:500:2000]]; title0='Mixed Layer Depth (Boyer)';
80			if doAnomalies; cc=scaleAnom*[-5:0.5:5]; end;
81			figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
82			myCaption={myYmeanTxt,'mean -- Boyer mld (m)'};
83			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
84
85			end;
86
87