Skip to content

Commit ea71aa0

Browse files
gaelforgetgaelforget
authored
Merge pull request #11 from gaelforget/develstuff01
lagrangian, vorticity, and divergence tools
2 parents d8feb4e + f606a96 commit ea71aa0

File tree

9 files changed

+5287
-9
lines changed

9 files changed

+5287
-9
lines changed

‎gcmfaces_devel/budget_mom_vort.STDOUT‎

Lines changed: 4826 additions & 0 deletions
Large diffs are not rendered by default.

‎gcmfaces_devel/budget_mom_vort.m‎

Lines changed: 190 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,190 @@
1+
function [out]=budget_mom_vort(nmRun,tim);
2+
% BUDGET_MOM_VORT performs a simple check of momentum budget closure for output
3+
% in nmRun ('diags.20160531' by default) and time tim (62772 by default)
4+
%
5+
% example 1: budget_mom_vort('diags.20160531',62772);
6+
%
7+
% example 2: execute commands below to compute time mean BPV budget
8+
%
9+
%list0=dir('diags.20160531/budg_tendU.*.meta'); nt=240;
10+
%list1=zeros(nt,1); for tt=1:nt; list1(tt)=str2num(list0(tt).name(12:end-5)); end;
11+
%
12+
%tic;
13+
%for tt=1:nt;
14+
% disp(tt);
15+
% [out]=budget_mom_vort('diags.20160531',list1(tt));
16+
% if tt==1;
17+
% list2=fieldnames(out);
18+
% for nn=1:length(list2);
19+
% ave.(list2{nn})=out.(list2{nn})/nt;
20+
% end;
21+
% else;
22+
% for nn=1:length(list2);
23+
% ave.(list2{nn})=ave.(list2{nn})+out.(list2{nn})/nt;
24+
% end;
25+
% end;
26+
%end;
27+
%toc;
28+
%
29+
%%save budget_mom_vort_ave.mat ave;
30+
%
31+
%rot.TOTUTEND=calc_UV_curl(ave.TOTUTEND,ave.TOTVTEND,1);
32+
%rot.AB_gU=calc_UV_curl(ave.AB_gU,ave.AB_gV,1);
33+
%rot.Um_Diss=calc_UV_curl(ave.Um_Diss,ave.Vm_Diss,1);
34+
%rot.Um_Advec=calc_UV_curl(ave.Um_Advec,ave.Vm_Advec,1);
35+
%rot.Um_dPHdx=calc_UV_curl(ave.Um_dPHdx,ave.Vm_dPHdx,1);
36+
%rot.Um_Ext=calc_UV_curl(ave.Um_Ext,ave.Vm_Ext,1);
37+
%rot.VISrI_Um=calc_UV_curl(ave.VISrI_Um,ave.VISrI_Vm,1);
38+
%rot.Um_dETANdx=calc_UV_curl(ave.Um_dETANdx,ave.Vm_dETANdy,1);
39+
%
40+
%%msk=mygrid.mskC(:,:,1);
41+
%msk=exch_T_N(mygrid.mskC(:,:,1));
42+
%for ff=1:msk.nFaces;
43+
% tmp1=msk{ff}(2:end-1,2:end-1);
44+
% for ii=-1:1; for jj=-1:1; tmp1=tmp1.*msk{ff}(2+ii:end-1+ii,2+jj:end-1+jj); end; end;
45+
% msk{ff}=tmp1;
46+
%end;
47+
%
48+
%list2=fieldnames(rot);
49+
%distXC=3*mygrid.DXC; distYC=3*mygrid.DYC;
50+
%for nn=1:length(list2);
51+
% tmp1=msk.*rot.(list2{nn});
52+
% tmp1=diffsmooth2D(tmp1,distXC,distYC);
53+
% eval([list2{nn} '=tmp1;']);
54+
%end;
55+
%
56+
%figureL; fac2=2;
57+
%subplot(3,2,1); qwckplot(TOTUTEND); caxis(fac2*[-1 1]*2e-10); colorbar; title('TEND');
58+
%subplot(3,2,2); qwckplot(Um_dETANdx+Um_dPHdx); caxis(fac2*[-1 1]*2e-10); colorbar; title('PRESS');
59+
%subplot(3,2,3); qwckplot(AB_gU+Um_Diss+VISrI_Um); caxis(fac2*[-1 1]*2e-10); colorbar; title('VARIOUS');
60+
%subplot(3,2,4); qwckplot(Um_Advec); caxis(fac2*[-1 1]*2e-10); colorbar; title('ADVEC');
61+
%subplot(3,2,5); qwckplot(Um_Ext); caxis(fac2*[-1 1]*2e-10); colorbar; title('EXT');
62+
%subplot(3,2,6); qwckplot(TOTUTEND-Um_dETANdx-Um_dPHdx-AB_gU...
63+
% -Um_Diss-VISrI_Um-Um_Advec-Um_Ext); caxis(fac2*[-1 1]*2e-10); colorbar; title('RESIDUAL');
64+
65+
gcmfaces_global;
66+
67+
doPlot=0;
68+
if nargin<1; nmRun='diags.20160531'; end;
69+
if nargin<2; tim=62772; end;
70+
71+
dir00='eccov4_release2_mom/';
72+
dir0=[dir00 nmRun '/'];
73+
dirGrid=[dir00 'GRID/'];
74+
75+
if isempty(mygrid); grid_load(dirGrid,5,'compact'); end;
76+
77+
RAW=rdmds2gcmfaces([dirGrid 'RAW']);
78+
RAS=rdmds2gcmfaces([dirGrid 'RAS']);
79+
nr=length(mygrid.RC); kk=1; fac0=1e3; fac1=1e3;
80+
doIce=0;
81+
82+
listVars={'TOTUTEND','AB_gU ','Um_Diss ','Um_Advec','Um_dPHdx','Um_Ext ','VISrI_Um'};
83+
listVars=deblank(listVars);
84+
for vv=1:7;
85+
fld=rdmds2gcmfaces([dir0 'budg_tendU'],tim,'rec',vv);
86+
%fld=fld(:,:,:,end);
87+
if vv==1; fld=fld/86400;
88+
elseif vv==7;
89+
fld=fld./repmat(RAW,[1 1 nr]);
90+
tmp1=mygrid.hFacW.*mk3D(mygrid.DRF,fld);
91+
tmp2=fld; tmp2(:,:,nr)=fld(:,:,nr);
92+
tmp2(:,:,1:nr-1)=fld(:,:,1:nr-1)-fld(:,:,2:nr);
93+
fld=tmp2./tmp1;
94+
end;
95+
tmp1=mygrid.hFacW.*mk3D(mygrid.DRF,fld);
96+
fld=nansum(tmp1.*fld,3);
97+
eval([listVars{vv} '=fld;']);
98+
eval(['out.' listVars{vv} '=fld;']);
99+
end;
100+
101+
listVars_V={'TOTVTEND','AB_gV ','Vm_Diss ','Vm_Advec','Vm_dPHdx','Vm_Ext ','VISrI_Vm'};
102+
listVars_V=deblank(listVars_V);
103+
for vv=1:7;
104+
fld=rdmds2gcmfaces([dir0 'budg_tendV'],tim,'rec',vv);
105+
%fld=fld(:,:,:,end);
106+
if vv==1; fld=fld/86400;
107+
elseif vv==7;
108+
fld=fld./repmat(RAS,[1 1 nr]);
109+
tmp1=mygrid.hFacS.*mk3D(mygrid.DRF,fld);
110+
tmp2=fld; tmp2(:,:,nr)=fld(:,:,nr);
111+
tmp2(:,:,1:nr-1)=fld(:,:,1:nr-1)-fld(:,:,2:nr);
112+
fld=tmp2./tmp1;
113+
end;
114+
tmp1=mygrid.hFacS.*mk3D(mygrid.DRF,fld);
115+
fld=nansum(tmp1.*fld,3);
116+
eval([listVars_V{vv} '=fld;']);
117+
eval(['out.' listVars_V{vv} '=fld;']);
118+
end;
119+
120+
ETAN=1/9.81*rdmds2gcmfaces([dir0 'budg_aveSURF'],tim,'rec',6);
121+
[Um_dETANdx,Vm_dETANdy]=calc_T_grad(-9.81*mygrid.mskC(:,:,1).*ETAN,0);
122+
123+
tmp1=nansum(mygrid.hFacW.*mk3D(mygrid.DRF,mygrid.hFacW),3);
124+
Um_dETANdx=Um_dETANdx.*tmp1;
125+
tmp1=nansum(mygrid.hFacS.*mk3D(mygrid.DRF,mygrid.hFacS),3);
126+
Vm_dETANdy=Vm_dETANdy.*tmp1;
127+
128+
out.Um_dETANdx=Um_dETANdx;
129+
out.Vm_dETANdy=Vm_dETANdy;
130+
131+
%%
132+
133+
if nargout==0;
134+
135+
m=mygrid.mskC(:,:,kk);
136+
137+
figureL;
138+
subplot(3,2,1); qwckplot(m.*VISrI_Um(:,:,kk));
139+
title('VISrI_Um (d/dr/rac of VISC)'); caxis(fac0*[-1 1]*1e-8); colorbar;
140+
subplot(3,2,2); qwckplot(m.*Um_dETANdx);
141+
title('Um_dETANdx (d/dx of PHI_SURF/g)'); caxis(fac0*[-1 1]*1e-5); colorbar;
142+
for jj=1:6;
143+
vv=listVars{jj+1};
144+
subplot(3,3,jj+3);
145+
eval(['fld=' vv ';']);
146+
qwckplot(m.*fld(:,:,kk));
147+
if strcmp(vv,'Um_Advec')|strcmp(vv,'Um_dPHdx')|strcmp(vv,'Um_Ext');
148+
caxis(fac0*[-1 1]*1e-5); colorbar;
149+
else;
150+
caxis(fac0*[-1 1]*1e-8); colorbar;
151+
end;
152+
title(listVars{jj+1});
153+
end;
154+
155+
%compute R.H.S.
156+
fld=Um_dPHdx(:,:,kk)+Um_dETANdx+Um_Advec(:,:,kk);
157+
fld=fld+Um_Diss(:,:,kk)+Um_Ext(:,:,kk)+AB_gU(:,:,kk);
158+
fld=fld-VISrI_Um(:,:,kk);
159+
160+
figureL;
161+
subplot(3,1,1); qwckplot(m.*TOTUTEND(:,:,kk));
162+
title('TOTUTEND'); caxis(fac1*[-1 1]*1e-8); colorbar;
163+
subplot(3,1,2); qwckplot(m.*fld);
164+
title('R.H.S.'); caxis(fac1*[-1 1]*1e-8); colorbar;
165+
subplot(3,1,3); qwckplot(m.*TOTUTEND(:,:,kk)-fld);
166+
title('residual'); caxis(fac1*[-1 1]*1e-11); colorbar;
167+
168+
%compute R.H.S.
169+
fld=Vm_dPHdx(:,:,kk)+Vm_dETANdy+Vm_Advec(:,:,kk);
170+
fld=fld+Vm_Diss(:,:,kk)+Vm_Ext(:,:,kk)+AB_gV(:,:,kk);
171+
fld=fld-VISrI_Vm(:,:,kk);
172+
173+
figureL;
174+
subplot(3,1,1); qwckplot(m.*TOTVTEND(:,:,kk));
175+
title('TOTVTEND'); caxis(fac1*[-1 1]*1e-8); colorbar;
176+
subplot(3,1,2); qwckplot(m.*fld);
177+
title('R.H.S.'); caxis(fac1*[-1 1]*1e-8); colorbar;
178+
subplot(3,1,3); qwckplot(m.*TOTVTEND(:,:,kk)-fld);
179+
title('residual'); caxis(fac1*[-1 1]*1e-11); colorbar;
180+
181+
if 0;
182+
figureL;
183+
subplot(2,1,1); qwckplot(m.*TOTUTEND(:,:,kk)-fld);
184+
title('residual'); caxis(fac1*[-1 1]*1e-8); colorbar;
185+
subplot(2,1,2); qwckplot(m.*VISrI_Um(:,:,kk));
186+
title('residual'); caxis(fac1*[-1 1]*1e-8); colorbar;
187+
end;
188+
189+
end;%if nargout==0;
190+

‎gcmfaces_devel/calc_UV_geos.m‎

Lines changed: 6 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -46,7 +46,7 @@
4646

4747
% mask=squeeze(mygrid.hFacC(:,:,1));
4848
% dzf=mygrid.DRF;
49-
nz=length(mygrid.RC);
49+
nz=size(P.f1,3);
5050

5151
%constants
5252
rhoconst=1029;
@@ -62,8 +62,8 @@
6262

6363
%replace NaN/1 mask with 0/1 mask:
6464
P(isnan(P))=0;
65-
mskW=mygrid.mskW; mskW(isnan(mskW))=0;
66-
mskS=mygrid.mskS; mskS(isnan(mskS))=0;
65+
mskW=mygrid.mskW(:,:,1:nz); mskW(isnan(mskW))=0;
66+
mskS=mygrid.mskS(:,:,1:nz); mskS(isnan(mskS))=0;
6767
%%weight average by portion that is filled
6868
%mskW=mygrid.hFacW;
6969
%mskS=mygrid.hFacS;
@@ -118,10 +118,9 @@
118118
Vg(:,:,iz)=vcur;
119119
end
120120

121-
mskW=mygrid.mskW; ii=find(isnan(mskW)); mskW(ii)=0;
122-
mskS=mygrid.mskS; ii=find(isnan(mskS)); mskS(ii)=0;
123-
Vg=Vg.*mskS;
124-
Ug=Ug.*mskW;
121+
mskW=mygrid.mskW(:,:,1:nz); ii=find(isnan(mskW)); mskW(ii)=0;
122+
mskS=mygrid.mskS(:,:,1:nz); ii=find(isnan(mskS)); mskS(ii)=0;
123+
Vg=Vg.*mskS; Ug=Ug.*mskW;
125124

126125

127126

‎gcmfaces_devel/calc_uvw_nodiv.m‎

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
function [Unodiv,Vnodiv,Wnodiv]=calc_uvw_nodiv(UVELMASS,VVELMASS);
2+
%[Unodiv,Vnodiv,Wnodiv]=calc_uvw_nodiv(UVELMASS,VVELMASS);
3+
% calculates non-divergent U,V,W fields from UVELMASS,VVELMASS
4+
5+
gcmfaces_global;
6+
7+
drf=mk3D(mygrid.DRF,UVELMASS);
8+
dxg=mk3D(mygrid.DXG,drf); dyg=mk3D(mygrid.DYG,drf);
9+
facW=drf.*dyg; facS=drf.*dxg;
10+
11+
%integrate vertically and apply surface mask:
12+
tmpU=nansum(facW.*UVELMASS,3).*mygrid.mskW(:,:,1);
13+
tmpV=nansum(facS.*VVELMASS,3).*mygrid.mskS(:,:,1);
14+
15+
%compute divergent transport:
16+
[tmpUdiv,tmpVdiv,tmpDivPot]=diffsmooth2D_div_inv(tmpU,tmpV);
17+
18+
%compute divergent velocity:
19+
tmpU=sum(facW.*mygrid.hFacW,3).*mygrid.mskW(:,:,1);
20+
tmpV=sum(facS.*mygrid.hFacS,3).*mygrid.mskS(:,:,1);
21+
tmpUdiv=tmpUdiv./tmpU; tmpVdiv=tmpVdiv./tmpV;
22+
23+
Udiv=mk3D(tmpUdiv,drf).*mygrid.hFacW;
24+
Vdiv=mk3D(tmpVdiv,drf).*mygrid.hFacS;
25+
Unodiv=UVELMASS-Udiv;
26+
Vnodiv=VVELMASS-Vdiv;
27+
28+
%verify result:
29+
%tmpU=nansum(facW.*Unodiv,3).*mygrid.mskW(:,:,1);
30+
%tmpV=nansum(facS.*Vnodiv,3).*mygrid.mskS(:,:,1);
31+
%[tmpUdiv2,tmpVdiv2,tmpDivPot2]=diffsmooth2D_div_inv(tmpU,tmpV);
32+
33+
%compute vertical component:
34+
tmp=calc_UV_conv(facW.*Unodiv,facS.*Vnodiv);
35+
Wnodiv=cumsum(tmp,3,'reverse')./mk3D(mygrid.RAC,drf);
36+
37+

‎gcmfaces_devel/flt_traj_init.m‎

Lines changed: 57 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,57 @@
1+
function []=prep_flt_init(dirOut);
2+
% prep_flt_init(dirOut) creates initial conditions for MITgcm/pkg/flt
3+
% that will be stored into dirOut ('init_flt/' by default).
4+
5+
%in summary:
6+
%- get tile map
7+
%- skip blank tiles
8+
%- loop over tile
9+
%- select all ocean points
10+
%- add uniform noise to i,j (in the -0.5 to 0.5 I think)
11+
%- output to file (single prec for ini or double for pickup)
12+
13+
%to do list:
14+
%- add lon, lat (vector?) input arguments and get i,j, etc using `gcmfaces_interp_coeffs`
15+
16+
if isempty(whos('dirOut')); dirOut=[pwd filesep 'init_flt' filesep]; end;
17+
if ~isdir(dirOut); mkdir(dirOut); end;
18+
19+
%%
20+
21+
gcmfaces_global; if isempty(mygrid); grid_load; end;
22+
map_tile=gcmfaces_loc_tile(30,30);
23+
24+
tmp1=convert2vector(mygrid.mskC(:,:,1).*map_tile);
25+
tmp1=unique(tmp1);
26+
list_tile=tmp1(find(~isnan(tmp1)));
27+
28+
map_i=NaN*map_tile; map_j=NaN*map_tile;
29+
[tmp_j,tmp_i]=meshgrid(1:30,1:30);
30+
for ff=1:map_i.nFaces;
31+
[m,n]=size(map_i{ff});
32+
map_i{ff}=repmat(tmp_i,[m/30 n/30]);
33+
map_j{ff}=repmat(tmp_j,[m/30 n/30]);
34+
end;
35+
36+
%%
37+
38+
vec_i=convert2vector(map_i);
39+
vec_j=convert2vector(map_j);
40+
vec_tile=convert2vector(map_tile.*mygrid.mskC(:,:,1));
41+
42+
kk=0;
43+
for ii=1:length(list_tile);
44+
jj=find(vec_tile==list_tile(ii));
45+
tmp_i=vec_i(jj); tmp_j=vec_j(jj);
46+
nn=length(jj);
47+
%
48+
tmp1=[nn 1 3600 0 0 9000 9 0 0];
49+
tmp2=[kk+[1:nn]' -ones(nn,1) tmp_i tmp_j ones(nn,1) zeros(nn,3) -ones(nn,1)];
50+
arrOut=[tmp1;tmp2]';
51+
%
52+
filOut=sprintf('%s/init_flt.%03d.001.data',dirOut,ii);
53+
write2file(filOut,arrOut,32);
54+
kk=kk+nn;
55+
end;
56+
57+

‎gcmfaces_devel/flt_traj_plot.m‎

Lines changed: 49 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,49 @@
1+
function []=flt_traj_plot(flts,proj,bgrd);
2+
% flt_traj_plot(flts) plots trajectories generated using pkg/flt
3+
%
4+
%A possible calling sequence:
5+
%
6+
%grid_load;
7+
%dirRun='run/';
8+
%[flts,data,header]=flt_traj_read([dirRun 'float_trajectories'],4);
9+
%flt_traj_plot(flts);
10+
11+
gcmfaces_global;
12+
13+
if isempty(whos('proj')); proj=0; end;
14+
15+
if proj==0;
16+
figure; orient landscape;
17+
col='rgbrgb';
18+
for ii=1:5;
19+
plot(mygrid.XC{ii}(1:5:end,1:5:end),mygrid.YC{ii}(1:5:end,1:5:end),[col(ii) 'x']); hold on;
20+
end;
21+
for ii=1:length(flts);
22+
x=flts(ii).x; y=flts(ii).y;
23+
jj=find(abs(diff(x))>90); x(jj)=NaN; y(jj)=NaN;
24+
x(end)=NaN; y(end)=NaN;
25+
plot(x,y,'k-'); hold on;
26+
end;
27+
end;
28+
29+
if proj>0;
30+
31+
figure; orient landscape;
32+
if isempty(whos('bgrd')); bgrd=mygrid.Depth; end;
33+
m_map_gcmfaces(bgrd,proj,{'myCmap','gray'},{'doCbar',0});
34+
for jj=1:3;
35+
if jj==1; mrk='r.';
36+
elseif jj==2; mrk='g.';
37+
else; mrk='b.';
38+
end;
39+
mrk='k.'; mrkini='go'; mrkend='co';
40+
for ii=jj:3:length(flts);
41+
lon=flts(ii).x; lat=flts(ii).y;
42+
if proj==1.2; lon(lon<20)=lon(lon<20)+360; end;
43+
m_map_gcmfaces({'plot',lon,lat,mrk,'MarkerSize',0.5},proj,{'doHold',1});
44+
m_map_gcmfaces({'plot',lon(1),lat(1),mrkini,'MarkerSize',2},proj,{'doHold',1});
45+
m_map_gcmfaces({'plot',lon(end-1),lat(end-1),mrkend,'MarkerSize',2},proj,{'doHold',1});
46+
end;
47+
end;
48+
end;
49+

0 commit comments

Comments
 (0)