| 1 |
gforget |
1.1 |
|
| 2 |
|
|
layersOffline=0; |
| 3 |
gforget |
1.4 |
layersEulerian=0; |
| 4 |
gforget |
1.1 |
integrateFromBottom=1; |
| 5 |
|
|
if isempty(whos('setDiagsParams')); setDiagsParams={}; end; |
| 6 |
|
|
|
| 7 |
gforget |
1.3 |
if length(setDiagsParams)==2; |
| 8 |
gforget |
1.1 |
layersName=setDiagsParams{1}; |
| 9 |
gforget |
1.3 |
layersLims=setDiagsParams{2}; |
| 10 |
gforget |
1.1 |
layersOffline=1; |
| 11 |
gforget |
1.4 |
layersEulerian=1;%hacked : should be 0 |
| 12 |
gforget |
1.1 |
elseif length(setDiagsParams)==1; |
| 13 |
gforget |
1.3 |
layersName=setDiagsParams{1}; |
| 14 |
gforget |
1.1 |
else; |
| 15 |
gforget |
1.3 |
layersName='1SLT'; |
| 16 |
gforget |
1.1 |
end; |
| 17 |
|
|
|
| 18 |
|
|
layersFileSuff=['LAYERS_' layersName]; |
| 19 |
|
|
if layersOffline; layersFileSuff=['LAYERS_offline_' layersName]; end; |
| 20 |
gforget |
1.3 |
|
| 21 |
|
|
%override default file name: |
| 22 |
|
|
%--------------------------- |
| 23 |
|
|
fileMat=['diags_set_' layersFileSuff]; |
| 24 |
gforget |
1.1 |
|
| 25 |
|
|
if userStep==1;%diags to be computed |
| 26 |
gforget |
1.3 |
listDiags=['layersParams gloOV gloThick gloDelThick gloBAR']; |
| 27 |
|
|
if sum([90 1170]~=mygrid.ioSize)==0; |
| 28 |
|
|
listDiags=[listDiags ' atlOV atlThick atlDelThick atlBAR']; |
| 29 |
|
|
listDiags=[listDiags ' pacindOV pacindThick pacindDelThick pacindBAR']; |
| 30 |
|
|
end; |
| 31 |
gforget |
1.1 |
elseif userStep==2&~layersOffline;%input files and variables |
| 32 |
|
|
listFlds={ ['LaUH' layersName],['LaVH' layersName],... |
| 33 |
|
|
['LaHw' layersName],['LaHs' layersName]}; |
| 34 |
|
|
listFldsNames=deblank(listFlds); |
| 35 |
|
|
listFiles={'layersDiags'}; |
| 36 |
|
|
listSubdirs={[dirModel 'diags/LAYERS/' ]}; |
| 37 |
|
|
%load layersLims consistent with online MITgcmpkg/layers |
| 38 |
|
|
layersLims=squeeze(rdmds([dirModel 'diags/LAYERS/layers' layersName])); |
| 39 |
|
|
elseif userStep==2&layersOffline;%input files and variables |
| 40 |
|
|
listFlds={'THETA','SALT','UVELMASS','VVELMASS','GM_PsiX','GM_PsiY'}; |
| 41 |
|
|
listFldsNames=deblank(listFlds); |
| 42 |
|
|
listFiles={'state_3d_set1','trsp_3d_set1','trsp_3d_set2'}; |
| 43 |
gforget |
1.3 |
elseif userStep==3; |
| 44 |
gforget |
1.1 |
|
| 45 |
|
|
%override default file name: |
| 46 |
|
|
%--------------------------- |
| 47 |
gforget |
1.3 |
fileMat=['diags_set_' layersFileSuff '_' num2str(tt) '.mat']; |
| 48 |
gforget |
1.1 |
|
| 49 |
|
|
if ~layersOffline; |
| 50 |
|
|
|
| 51 |
|
|
eval(['fldU=LaUH' layersName '; fldV=LaVH' layersName ';']); |
| 52 |
|
|
if ~isempty(whos(['LaHw' layersName])); |
| 53 |
|
|
eval(['fldHw=LaHw' layersName '; fldHs=LaHs' layersName ';']); |
| 54 |
|
|
else; |
| 55 |
|
|
fldHw=[]; fldHs=[]; |
| 56 |
|
|
end; |
| 57 |
|
|
layersGrid=(layersLims(1:end-1)+layersLims(2:end))'/2; |
| 58 |
|
|
|
| 59 |
|
|
else; |
| 60 |
|
|
|
| 61 |
|
|
if ~layersEulerian; |
| 62 |
|
|
[fldUbolus,fldVbolus,fldWbolus]=calc_bolus(GM_PsiX,GM_PsiY); |
| 63 |
|
|
fldUbolus=fldUbolus.*mygrid.mskW; fldVbolus=fldVbolus.*mygrid.mskS; |
| 64 |
|
|
UVELMASS=UVELMASS+fldUbolus; VVELMASS=VVELMASS+fldVbolus; |
| 65 |
|
|
end; |
| 66 |
|
|
|
| 67 |
|
|
U=UVELMASS.*mk3D(mygrid.DRF,UVELMASS); |
| 68 |
|
|
V=VVELMASS.*mk3D(mygrid.DRF,VVELMASS); |
| 69 |
|
|
|
| 70 |
|
|
if strcmp(layersName,'2TH'); |
| 71 |
|
|
tracer=THETA; |
| 72 |
gforget |
1.4 |
if isempty(whos('layersLims')); layersLims=[-2:35]; end; |
| 73 |
gforget |
1.1 |
elseif strcmp(layersName,'1SLT'); |
| 74 |
|
|
tracer=SALT; |
| 75 |
gforget |
1.4 |
if isempty(whos('layersLims')); layersLims=[33:0.1:38]; end; |
| 76 |
gforget |
1.1 |
elseif strcmp(layersName,'3RHO'); |
| 77 |
|
|
P=mk3D(-mygrid.RC,THETA); |
| 78 |
|
|
t=convert2vector(THETA); |
| 79 |
|
|
s=convert2vector(SALT); |
| 80 |
|
|
p=convert2vector(P); |
| 81 |
|
|
pref=2000+0*p;%could be made an extra param |
| 82 |
|
|
[rhop,rhpis,rhor] = density(t,s,p,pref); |
| 83 |
|
|
tracer=convert2vector(rhor); |
| 84 |
gforget |
1.4 |
if isempty(whos('layersLims')); layersLims=1000+[30:0.1:38]'; end; |
| 85 |
gforget |
1.1 |
end |
| 86 |
|
|
|
| 87 |
|
|
layersGrid=(layersLims(1:end-1)+layersLims(2:end))'/2; |
| 88 |
|
|
tic; [fldU,fldV]=layers_remap({U,V},'extensive',tracer,layersGrid,1); toc;%hacked: should be 2 |
| 89 |
|
|
fldHw=[]; fldHs=[]; |
| 90 |
|
|
|
| 91 |
|
|
end;%if ~layersOffline; |
| 92 |
gforget |
1.3 |
|
| 93 |
|
|
if sum([90 1170]~=mygrid.ioSize)==0; |
| 94 |
|
|
listBasins=[1:3]; |
| 95 |
|
|
else; |
| 96 |
|
|
listBasins=1; |
| 97 |
|
|
end; |
| 98 |
|
|
|
| 99 |
|
|
for bb=listBasins; |
| 100 |
|
|
|
| 101 |
|
|
%mask : global, atlantic or Pac+Ind |
| 102 |
|
|
if bb==1; mskC=mygrid.mskC(:,:,1); mskW=mygrid.mskW(:,:,1); mskS=mygrid.mskS(:,:,1); |
| 103 |
|
|
elseif bb==2; [mskC,mskW,mskS]=v4_basin({'atlExt'}); |
| 104 |
|
|
elseif bb==3; [mskC,mskW,mskS]=v4_basin({'pacExt','indExt'}); |
| 105 |
gforget |
1.1 |
end; |
| 106 |
|
|
mskC3d=1*(mskC>0); mskW3d=1*(mskW>0); mskS3d=1*(mskS>0); |
| 107 |
|
|
mskC3d=mk3D(mskC3d,fldU); mskW3d=mk3D(mskW3d,fldU); mskS3d=mk3D(mskS3d,fldU); |
| 108 |
gforget |
1.2 |
|
| 109 |
gforget |
1.3 |
%the overturning streamfunction computation itself |
| 110 |
|
|
layersOV=calc_overturn(fldU.*mskW3d,fldV.*mskS3d,1,{'dh'}); |
| 111 |
|
|
|
| 112 |
gforget |
1.1 |
%the associated barotropic streamfunction (for checking) |
| 113 |
gforget |
1.3 |
layersBAR=calc_barostream(fldU.*mskW3d,fldV.*mskS3d,1,{'dh'}); |
| 114 |
gforget |
1.1 |
|
| 115 |
|
|
%the associated thickness |
| 116 |
|
|
if isempty(fldHw); |
| 117 |
|
|
layersThick=[]; |
| 118 |
|
|
layersDelThick=[]; |
| 119 |
|
|
else; |
| 120 |
|
|
%interpolate to tracer points |
| 121 |
|
|
fldH=0*fldHw; |
| 122 |
gforget |
1.3 |
[fldHwp,fldHsp]=exch_UV(LaHw1SLT.*mskW3d,LaHs1SLT.*mskS3d); |
| 123 |
gforget |
1.1 |
for iF=1:fldH.nFaces; |
| 124 |
|
|
tmpA=fldHwp{iF}(2:end,:,:); |
| 125 |
|
|
tmpB=fldHwp{iF}(1:end-1,:,:); |
| 126 |
|
|
tmpC=fldHsp{iF}(:,2:end,:); |
| 127 |
|
|
tmpD=fldHsp{iF}(:,1:end-1,:); |
| 128 |
|
|
tmpTot=tmpA+tmpB+tmpC+tmpD; |
| 129 |
|
|
tmpNb=1*(tmpA~=0)+1*(tmpB~=0)+1*(tmpC~=0)+1*(tmpD~=0); |
| 130 |
|
|
jj=find(tmpNb>0); tmpTot(jj)=tmpTot(jj)./tmpNb(jj); |
| 131 |
|
|
jj=find(isnan(tmpTot)); tmpTot(jj)=0; |
| 132 |
|
|
fldH{iF}=tmpTot; |
| 133 |
|
|
end; |
| 134 |
|
|
%compute zonal mean |
| 135 |
|
|
fldH=calc_zonmean_T(fldH.*mskC3d); |
| 136 |
|
|
%integrate to overturning points |
| 137 |
|
|
if integrateFromBottom; |
| 138 |
|
|
layersThick=[flipdim(cumsum(flipdim(fldH,2),2),2) zeros(size(fldH,1),1)]; |
| 139 |
|
|
else; |
| 140 |
|
|
layersThick=[zeros(size(fldH,1),1) cumsum(fldH,2)]; |
| 141 |
|
|
end; |
| 142 |
|
|
%compute dH/dLayer |
| 143 |
|
|
layersDelThick=[zeros(size(fldH,1),1) fldH./( ones(size(fldH,1),1)*diff(layersLims') )]; |
| 144 |
|
|
end; |
| 145 |
gforget |
1.3 |
|
| 146 |
|
|
%store to global, atlantic or Pac+Ind arrays: |
| 147 |
|
|
if bb==1; |
| 148 |
|
|
gloOV=layersOV; gloThick=layersThick; gloDelThick=layersDelThick; gloBAR=layersBAR; |
| 149 |
|
|
elseif bb==2; |
| 150 |
|
|
kk=find(mygrid.LATS<-35|mygrid.LATS>70); |
| 151 |
|
|
layersOV(kk,:)=NaN; if ~isempty(layersThick); layersThick(kk,:)=NaN; layersDelThick(kk,:)=NaN; end; |
| 152 |
|
|
atlOV=layersOV; atlThick=layersThick; atlDelThick=layersDelThick; atlBAR=layersBAR; |
| 153 |
|
|
elseif bb==3; |
| 154 |
|
|
kk=find(mygrid.LATS<-35|mygrid.LATS>65); |
| 155 |
|
|
layersOV(kk,:)=NaN; if ~isempty(layersThick); layersThick(kk,:)=NaN; layersDelThick(kk,:)=NaN; end; |
| 156 |
|
|
pacindOV=layersOV; pacindThick=layersThick; pacindDelThick=layersDelThick; pacindBAR=layersBAR; |
| 157 |
|
|
end; |
| 158 |
|
|
end; |
| 159 |
gforget |
1.1 |
|
| 160 |
|
|
layersParams.name=layersName; |
| 161 |
|
|
layersParams.LC=layersGrid; |
| 162 |
|
|
layersParams.LF=layersLims; |
| 163 |
|
|
layersParams.isOffline=layersOffline; |
| 164 |
|
|
layersParams.suffix=layersFileSuff; |
| 165 |
|
|
layersParams.isEulerian=layersEulerian; |
| 166 |
gforget |
1.3 |
|
| 167 |
|
|
|
| 168 |
|
|
%===================== COMPUTATIONAL SEQUENCE ENDS =========================% |
| 169 |
|
|
%===================== PLOTTING SEQUENCE BEGINS =========================% |
| 170 |
|
|
|
| 171 |
|
|
elseif userStep==-1;%plotting |
| 172 |
|
|
|
| 173 |
|
|
|
| 174 |
|
|
%meridional streamfunction (Eulerian) : |
| 175 |
|
|
fld=mean(alldiag.gloOV(:,:,tt),3); fld(fld==0)=NaN; |
| 176 |
|
|
X=mygrid.LATS*ones(1,length(alldiag.layersParams.LF)); Y=ones(length(mygrid.LATS),1)*(alldiag.layersParams.LF'); |
| 177 |
|
|
cc=[[-50:10:-30] [-24:3:24] [30:10:50]]; title0='Meridional Stream Function'; |
| 178 |
|
|
if doAnomalies; cc=scaleAnom*[-1:0.1:1]; end; |
| 179 |
|
|
figureL; set(gcf,'Renderer','zbuffer'); %set(gcf,'Units','Normalized','Position',[0.05 0.1 0.4 0.8]); |
| 180 |
|
|
pcolor(X,Y,fld); shading interp; cbar=gcmfaces_cmap_cbar(cc); title(title0); |
| 181 |
|
|
myCaption={myYmeanTxt,'mean -- overturning streamfunction (Sv)'}; |
| 182 |
|
|
if addToTex; write2tex(fileTex,2,myCaption,gcf); end; |
| 183 |
|
|
|
| 184 |
|
|
if multiBasins; |
| 185 |
|
|
%meridional streamfunction (Eulerian): |
| 186 |
|
|
fld=mean(alldiag.atlOV(:,:,tt),3); fld(fld==0)=NaN; |
| 187 |
|
|
X=mygrid.LATS*ones(1,length(alldiag.layersParams.LF)); Y=ones(length(mygrid.LATS),1)*(alldiag.layersParams.LF'); |
| 188 |
|
|
cc=[[-50:10:-30] [-24:3:24] [30:10:50]]; title0='Atlantic Meridional Stream Function'; |
| 189 |
|
|
if doAnomalies; cc=scaleAnom*[-1:0.1:1]; end; |
| 190 |
|
|
figureL; set(gcf,'Renderer','zbuffer'); %set(gcf,'Units','Normalized','Position',[0.05 0.1 0.4 0.8]); |
| 191 |
|
|
pcolor(X,Y,fld); shading interp; cbar=gcmfaces_cmap_cbar(cc); title(title0); |
| 192 |
|
|
myCaption={myYmeanTxt,'mean -- Atlantic overturning streamfunction (Sv)'}; |
| 193 |
|
|
if addToTex; write2tex(fileTex,2,myCaption,gcf); end; |
| 194 |
|
|
%meridional streamfunction (residual): |
| 195 |
|
|
fld=mean(alldiag.pacindOV(:,:,tt),3); fld(fld==0)=NaN; |
| 196 |
|
|
X=mygrid.LATS*ones(1,length(alldiag.layersParams.LF)); Y=ones(length(mygrid.LATS),1)*(alldiag.layersParams.LF'); |
| 197 |
|
|
cc=[[-50:10:-30] [-24:3:24] [30:10:50]]; title0='Pac+Ind Meridional Stream Function'; |
| 198 |
|
|
if doAnomalies; cc=scaleAnom*[-1:0.1:1]; end; |
| 199 |
|
|
figureL; set(gcf,'Renderer','zbuffer'); %set(gcf,'Units','Normalized','Position',[0.05 0.1 0.4 0.8]); |
| 200 |
|
|
pcolor(X,Y,fld); shading interp; cbar=gcmfaces_cmap_cbar(cc); title(title0); |
| 201 |
|
|
myCaption={myYmeanTxt,'mean -- Pac+Ind overturning streamfunction (Sv)'}; |
| 202 |
|
|
if addToTex; write2tex(fileTex,2,myCaption,gcf); end; |
| 203 |
|
|
end; |
| 204 |
|
|
|
| 205 |
gforget |
1.1 |
|
| 206 |
|
|
end; |
Parent Directory
| 
Revision Log
| 
Revision Graph