/[MITgcm]/MITgcm_contrib/gael/matlab_class/gcmfaces_diags/diags_set_LAYERS.m
ViewVC logotype

Annotation of /MITgcm_contrib/gael/matlab_class/gcmfaces_diags/diags_set_LAYERS.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (hide annotations) (download)
Thu Nov 1 19:07:17 2012 UTC (12 years, 8 months ago) by gforget
Branch: MAIN
- more modular, and suposedly more user-friendly, version
  of ecco_v4/comp_driver.m and basic_diags_ecco.m

1 gforget 1.1
2     layersOffline=0;
3     layersEulerian=NaN;
4     integrateFromBottom=1;
5     if isempty(whos('setDiagsParams')); setDiagsParams={}; end;
6    
7     if length(setDiagsParams)==3;
8     layersName=setDiagsParams{1};
9     nameBasin=setDiagsParams{2};
10     layersLims=setDiagsParams{3};
11     layersOffline=1;
12     layersEulerian=1;%hacked : should be 0
13     elseif length(setDiagsParams)==2;
14     layersName=setDiagsParams{1}; nameBasin=setDiagsParams{2};
15     elseif length(setDiagsParams)==1;
16     layersName=setDiagsParams{1}; nameBasin='';
17     else;
18     layersName='1SLT'; nameBasin='';
19     end;
20    
21     layersFileSuff=['LAYERS_' layersName];
22     if layersOffline; layersFileSuff=['LAYERS_offline_' layersName]; end;
23     if ~isempty(nameBasin); layersFileSuff=[layersFileSuff nameBasin]; end;
24    
25     if userStep==1;%diags to be computed
26     listDiags=['layersParams layersOV layersThick layersDelThick layersBAR'];
27     elseif userStep==2&~layersOffline;%input files and variables
28     listFlds={ ['LaUH' layersName],['LaVH' layersName],...
29     ['LaHw' layersName],['LaHs' layersName]};
30     listFldsNames=deblank(listFlds);
31     listFiles={'layersDiags'};
32     listSubdirs={[dirModel 'diags/LAYERS/' ]};
33     %load layersLims consistent with online MITgcmpkg/layers
34     layersLims=squeeze(rdmds([dirModel 'diags/LAYERS/layers' layersName]));
35     elseif userStep==2&layersOffline;%input files and variables
36     listFlds={'THETA','SALT','UVELMASS','VVELMASS','GM_PsiX','GM_PsiY'};
37     listFldsNames=deblank(listFlds);
38     listFiles={'state_3d_set1','trsp_3d_set1','trsp_3d_set2'};
39     else;
40    
41     %override default file name:
42     %---------------------------
43     tmp1=layersFileSuff;
44     fileMat=['diags_set_' tmp1 '_' num2str(tt) '.mat'];
45    
46     if ~layersOffline;
47    
48     eval(['fldU=LaUH' layersName '; fldV=LaVH' layersName ';']);
49     if ~isempty(whos(['LaHw' layersName]));
50     eval(['fldHw=LaHw' layersName '; fldHs=LaHs' layersName ';']);
51     else;
52     fldHw=[]; fldHs=[];
53     end;
54     layersGrid=(layersLims(1:end-1)+layersLims(2:end))'/2;
55    
56     else;
57    
58     if ~layersEulerian;
59     [fldUbolus,fldVbolus,fldWbolus]=calc_bolus(GM_PsiX,GM_PsiY);
60     fldUbolus=fldUbolus.*mygrid.mskW; fldVbolus=fldVbolus.*mygrid.mskS;
61     UVELMASS=UVELMASS+fldUbolus; VVELMASS=VVELMASS+fldVbolus;
62     end;
63    
64     U=UVELMASS.*mk3D(mygrid.DRF,UVELMASS);
65     V=VVELMASS.*mk3D(mygrid.DRF,VVELMASS);
66    
67     if strcmp(layersName,'2TH');
68     tracer=THETA;
69     if isempty(layersLims); layersLims=[-2:35]; end;
70     elseif strcmp(layersName,'1SLT');
71     tracer=SALT;
72     if isempty(layersLims); layersLims=[33:0.1:38]; end;
73     elseif strcmp(layersName,'3RHO');
74     P=mk3D(-mygrid.RC,THETA);
75     t=convert2vector(THETA);
76     s=convert2vector(SALT);
77     p=convert2vector(P);
78     pref=2000+0*p;%could be made an extra param
79     [rhop,rhpis,rhor] = density(t,s,p,pref);
80     tracer=convert2vector(rhor);
81     if isempty(layersLims); layersLims=1000+[30:0.1:38]'; end;
82     end
83    
84     layersGrid=(layersLims(1:end-1)+layersLims(2:end))'/2;
85     tic; [fldU,fldV]=layers_remap({U,V},'extensive',tracer,layersGrid,1); toc;%hacked: should be 2
86     fldHw=[]; fldHs=[];
87    
88     end;%if ~layersOffline;
89    
90     %the very overturning streamfunction computation
91     if ~isempty(nameBasin); [mskC,mskW,mskS]=v4_basin(nameBasin);
92     else; mskC=mygrid.mskC(:,:,1); mskW=mygrid.mskW(:,:,1); mskS=mygrid.mskS(:,:,1);
93     end;
94     mskC3d=1*(mskC>0); mskW3d=1*(mskW>0); mskS3d=1*(mskS>0);
95     mskC3d=mk3D(mskC3d,fldU); mskW3d=mk3D(mskW3d,fldU); mskS3d=mk3D(mskS3d,fldU);
96     fldU=fldU.*mskW3d; fldV=fldV.*mskS3d;
97     layersOV=calc_overturn(fldU,fldV,1,{'dh'});
98    
99     %the associated barotropic streamfunction (for checking)
100     layersBAR=calc_barostream(fldU,fldV,1,{'dh'});
101    
102     %the associated thickness
103     if isempty(fldHw);
104     layersThick=[];
105     layersDelThick=[];
106     else;
107     %interpolate to tracer points
108     fldH=0*fldHw;
109     [fldHwp,fldHsp]=exch_UV(LaHw1SLT,LaHs1SLT);
110     for iF=1:fldH.nFaces;
111     tmpA=fldHwp{iF}(2:end,:,:);
112     tmpB=fldHwp{iF}(1:end-1,:,:);
113     tmpC=fldHsp{iF}(:,2:end,:);
114     tmpD=fldHsp{iF}(:,1:end-1,:);
115     tmpTot=tmpA+tmpB+tmpC+tmpD;
116     tmpNb=1*(tmpA~=0)+1*(tmpB~=0)+1*(tmpC~=0)+1*(tmpD~=0);
117     jj=find(tmpNb>0); tmpTot(jj)=tmpTot(jj)./tmpNb(jj);
118     jj=find(isnan(tmpTot)); tmpTot(jj)=0;
119     fldH{iF}=tmpTot;
120     end;
121     %compute zonal mean
122     fldH=calc_zonmean_T(fldH.*mskC3d);
123     %integrate to overturning points
124     if integrateFromBottom;
125     layersThick=[flipdim(cumsum(flipdim(fldH,2),2),2) zeros(size(fldH,1),1)];
126     else;
127     layersThick=[zeros(size(fldH,1),1) cumsum(fldH,2)];
128     end;
129     %compute dH/dLayer
130     layersDelThick=[zeros(size(fldH,1),1) fldH./( ones(size(fldH,1),1)*diff(layersLims') )];
131     end;
132    
133     layersParams.name=layersName;
134     layersParams.LC=layersGrid;
135     layersParams.LF=layersLims;
136     layersParams.isOffline=layersOffline;
137     layersParams.suffix=layersFileSuff;
138     layersParams.isEulerian=layersEulerian;
139    
140     end;

  ViewVC Help
Powered by ViewVC 1.1.22