matlab_class/gcmfaces_diags/diags_set_MLD.m

3.3


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;

  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)'};

  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:240),3); mon='March';
    else; fld_seas=nanmedian(fld(:,:,9:12:240),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;


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