matlab_class/gcmfaces_diags/diags_set_B.m


if userStep==1;%diags to be computed
    %here _s stands for cumulative sum and _2s for cumulative sum of squares
    %the _s should be stated first
    listDiags=['THETA_s SALT_s SIatmQnt_s SIatmFW_s oceQnet_s oceFWflx_s ' ...
               'fldTZ_s fldTM_s curlTau_s fldETAN_s fldETANLEADS_s fldMLD_s ' ...
               'THETA_2s SALT_2s SIatmQnt_2s SIatmFW_2s oceQnet_2s oceFWflx_2s ' ...
               'fldTZ_2s fldTM_2s curlTau_2s fldETAN_2s fldETANLEADS_2s fldMLD_2s'];
elseif userStep==2;%input files and variables
    listFlds={    'THETA','SALT','SIatmQnt','SIatmFW ','oceQnet ','oceFWflx','oceTAUX','oceTAUY','ETAN','sIceLoad','MXLDEPTH'};
    listFldsNames=deblank(listFlds);
    listFiles={'state_3d_set1','state_2d_set1','other_2d_set1'};
    listSubdirs={[dirModel 'diags/OTHER/' ],[dirModel 'diags/STATE/' ]};
elseif userStep==3;%computational part;
    %mask fields:
    SIatmQnt=SIatmQnt.*mygrid.mskC(:,:,1);
    SIatmFW=SIatmFW.*mygrid.mskC(:,:,1);
    oceQnet=oceQnet.*mygrid.mskC(:,:,1);
    oceFWflx=oceFWflx.*mygrid.mskC(:,:,1);
    fldTX=oceTAUX.*mygrid.mskW(:,:,1);
    fldTY=oceTAUY.*mygrid.mskS(:,:,1);
    %compute Eastward/Northward wind stresses:
    [fldTZ,fldTM]=calc_UEVNfromUXVY(fldTX,fldTY);
    %compute wind stress curl:
    curlTau=calc_UV_curl(fldTX, fldTY,1 );%the doMask argument should not matter as msk was already applied
    %mask and re-arrange fields:
    fldETAN=ETAN.*mygrid.mskC(:,:,1);
    fldETANLEADS=(ETAN+sIceLoad/myparms.rhoconst).*mygrid.mskC(:,:,1);
    fldMLD=MXLDEPTH.*mygrid.mskC(:,:,1);
    %
    THETA=THETA.*mygrid.mskC;
    SALT=SALT.*mygrid.mskC;
    %
    if ii>1;
      fileMatPrev=['diags_set_' tmp1 '_' num2str(listTimes(ii-1)) '.mat'];
      listTimesBak=listTimes; 
      load([dirMat fileMatPrev]);
      listTimes=listTimesBak; 
    end;
    for jj=1:length(listDiags)/2;
        myDiag=listDiags{jj}(1:end-2);
        if ii==1;
            eval([myDiag '_s=0*' myDiag ';']);
            eval([myDiag '_2s=0*' myDiag '.^2;']);
        end;
        eval([myDiag '_s=' myDiag '_s+' myDiag ';']);
        eval([myDiag '_2s=' myDiag '_2s+' myDiag '.^2;']);
    end;

end;
1	gforget	1.1
2			if userStep==1;%diags to be computed
3	gforget	1.2	%here _s stands for cumulative sum and _2s for cumulative sum of squares
4			%the _s should be stated first
5			listDiags=['THETA_s SALT_s SIatmQnt_s SIatmFW_s oceQnet_s oceFWflx_s ' ...
6			'fldTZ_s fldTM_s curlTau_s fldETAN_s fldETANLEADS_s fldMLD_s ' ...
7			'THETA_2s SALT_2s SIatmQnt_2s SIatmFW_2s oceQnet_2s oceFWflx_2s ' ...
8			'fldTZ_2s fldTM_2s curlTau_2s fldETAN_2s fldETANLEADS_2s fldMLD_2s'];
9	gforget	1.1	elseif userStep==2;%input files and variables
10			listFlds={ 'THETA','SALT','SIatmQnt','SIatmFW ','oceQnet ','oceFWflx','oceTAUX','oceTAUY','ETAN','sIceLoad','MXLDEPTH'};
11			listFldsNames=deblank(listFlds);
12			listFiles={'state_3d_set1','state_2d_set1','other_2d_set1'};
13			listSubdirs={[dirModel 'diags/OTHER/' ],[dirModel 'diags/STATE/' ]};
14			elseif userStep==3;%computational part;
15			%mask fields:
16			SIatmQnt=SIatmQnt.*mygrid.mskC(:,:,1);
17			SIatmFW=SIatmFW.*mygrid.mskC(:,:,1);
18			oceQnet=oceQnet.*mygrid.mskC(:,:,1);
19			oceFWflx=oceFWflx.*mygrid.mskC(:,:,1);
20			fldTX=oceTAUX.*mygrid.mskW(:,:,1);
21			fldTY=oceTAUY.*mygrid.mskS(:,:,1);
22			%compute Eastward/Northward wind stresses:
23			[fldTZ,fldTM]=calc_UEVNfromUXVY(fldTX,fldTY);
24			%compute wind stress curl:
25			curlTau=calc_UV_curl(fldTX, fldTY,1 );%the doMask argument should not matter as msk was already applied
26			%mask and re-arrange fields:
27			fldETAN=ETAN.*mygrid.mskC(:,:,1);
28			fldETANLEADS=(ETAN+sIceLoad/myparms.rhoconst).*mygrid.mskC(:,:,1);
29			fldMLD=MXLDEPTH.*mygrid.mskC(:,:,1);
30			%
31			THETA=THETA.*mygrid.mskC;
32			SALT=SALT.*mygrid.mskC;
33	gforget	1.2	%
34			if ii>1;
35			fileMatPrev=['diags_set_' tmp1 '_' num2str(listTimes(ii-1)) '.mat'];
36	gforget	1.3	listTimesBak=listTimes;
37	gforget	1.2	load([dirMat fileMatPrev]);
38	gforget	1.3	listTimes=listTimesBak;
39	gforget	1.2	end;
40			for jj=1:length(listDiags)/2;
41			myDiag=listDiags{jj}(1:end-2);
42			if ii==1;
43			eval([myDiag '_s=0*' myDiag ';']);
44			eval([myDiag '_2s=0*' myDiag '.^2;']);
45			end;
46			eval([myDiag '_s=' myDiag '_s+' myDiag ';']);
47			eval([myDiag '_2s=' myDiag '_2s+' myDiag '.^2;']);
48			end;
49
50	gforget	1.1	end;