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;%plotting;

  list_var={'fldMldKara','fldMldSuga','fldMldBoyer'};

  list_tit={' mixed layer depth per Kara formula (m)',...
            ' mixed layer depth per Suga formula (m)',...
            ' mixed layer depth per Boyer M. formula (m)'};

  if myparms.diagsAreMonthly==1;
  for seas=1:2;
  for vv=1:length(list_var);

    eval(['fld=alldiag.' list_var{vv} ';']);

    %compute mean march and september fields
    if seas==1; fld_seas=nanmedian(fld(:,:,3:12:nt),3); mon='March';
    else; fld_seas=nanmedian(fld(:,:,9:12:nt),3); mon='September';
    end;
    fld_seas(fld_seas==0)=NaN;
  
    %plot
    cc=[[0:20:100] [150:50:300] 400 [500:200:1100] [1500:500:2000]];
    if doAnomalies; cc=scaleAnom*[-5:0.5:5]; end;
    figureL; set(gcf,'Renderer','zbuffer');
    m_map_gcmfaces(fld_seas,0,{'myCaxis',cc});
    myCaption={myYmeanTxt,mon,' mean -- ',list_tit{vv}};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

  end;%for vv=1:length(list_var);
  end;%for seas=1:2;
  end;%if myparms.diagsAreMonthly==1

end;


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