| 1 |
gforget |
1.1 |
function [Ug,Vg]=calc_UV_geos(P); |
| 2 |
|
|
%[Ug, Vg]=calc_UV_geos(P); |
| 3 |
|
|
% computes geostrophic velocities from horizontal pressure anomaly field P |
| 4 |
|
|
% |
| 5 |
|
|
% example: |
| 6 |
|
|
% |
| 7 |
|
|
% %read actual model velocity field |
| 8 |
|
|
% U=mean(read_nctiles('release2_climatology/nctiles_climatology/UVELMASS/UVELMASS'),4); |
| 9 |
|
|
% V=mean(read_nctiles('release2_climatology/nctiles_climatology/VVELMASS/VVELMASS'),4); |
| 10 |
|
|
% %read hydrostatic pressure gradient computed along R* surface: |
| 11 |
|
|
% P=mean(read_nctiles('release2_climatology/nctiles_climatology/PHIHYD/PHIHYD'),4); |
| 12 |
|
|
% %correct, approximately, for slope of R* surface: |
| 13 |
|
|
% ETAN=mean(read_nctiles('release2_climatology/nctiles_climatology/ETAN/ETAN'),3); |
| 14 |
|
|
% tmp1=mygrid.mskC.*mk3D(ETAN,P);%free surface height |
| 15 |
|
|
% tmp2=mk3D(mygrid.DRF,P).*(1-mygrid.hFacC);%grounded thickness |
| 16 |
|
|
% tmp3=mygrid.mskC.*mk3D(-mygrid.RC,P)-1/2*tmp2;%depth of R* center points |
| 17 |
|
|
% tmp4=mygrid.mskC.*mk3D(mygrid.Depth,P);%bottom depth |
| 18 |
|
|
% tmp5=1-(tmp1+tmp3)./(tmp1+tmp4); |
| 19 |
|
|
% P=P+tmp5.*mk3D(9.81*ETAN,P); |
| 20 |
|
|
% |
| 21 |
|
|
% %compute geostrophy |
| 22 |
|
|
% [Ug,Vg]=calc_UV_geos(P); |
| 23 |
|
|
% |
| 24 |
|
|
% %display results |
| 25 |
|
|
% kk=10; cc=[-1 1]*0.1; |
| 26 |
|
|
% kk=30; cc=[-1 1]*0.05; |
| 27 |
|
|
% kk=40; cc=[-1 1]*0.02; |
| 28 |
|
|
% figure; orient tall; m=mygrid.mskC(:,:,kk); |
| 29 |
|
|
% [tmpu,tmpv]=calc_UEVNfromUXVY(U(:,:,kk),V(:,:,kk)); |
| 30 |
|
|
% [tmpug,tmpvg]=calc_UEVNfromUXVY(Ug(:,:,kk),Vg(:,:,kk)); |
| 31 |
|
|
% subplot(3,1,1); qwckplot(m.*tmpu); caxis(cc); colorbar; title('zonal flow'); |
| 32 |
|
|
% subplot(3,1,2); qwckplot(m.*tmpug); caxis(cc); colorbar; title('geostrophic flow'); |
| 33 |
|
|
% subplot(3,1,3); qwckplot(m.*tmpu-tmpug); caxis(cc); colorbar; title('difference'); |
| 34 |
|
|
% |
| 35 |
|
|
|
| 36 |
|
|
% development notes: |
| 37 |
|
|
% is P assumed NaN-masked? |
| 38 |
|
|
% add to example_transports? Along with Ekman transport? |
| 39 |
|
|
% hard-coded rhoconst, g, gamma should be stated in help section |
| 40 |
|
|
% k loop should be possible to remove; once done also for calc_T_grad |
| 41 |
|
|
% should allow for application to ETAN x appropriate constant as 2D input |
| 42 |
|
|
|
| 43 |
|
|
gcmfaces_global; |
| 44 |
|
|
|
| 45 |
|
|
% mask=squeeze(mygrid.hFacC(:,:,1)); |
| 46 |
|
|
% dzf=mygrid.DRF; |
| 47 |
|
|
nz=length(mygrid.RC); |
| 48 |
|
|
|
| 49 |
|
|
%constants |
| 50 |
|
|
rhoconst=1029; |
| 51 |
|
|
g=9.81; |
| 52 |
|
|
omega=7.27e-5; |
| 53 |
|
|
|
| 54 |
|
|
[dxC, dyC]=exch_UV_N(mygrid.DXC, mygrid.DYC); |
| 55 |
|
|
|
| 56 |
|
|
%f=2*omega*sind(mygrid.YC); |
| 57 |
|
|
%fmat=exch_T_N(f); |
| 58 |
|
|
f=2*omega*sind(mygrid.YG); |
| 59 |
|
|
fmat=exch_Z(f); |
| 60 |
|
|
|
| 61 |
|
|
%replace NaN/1 mask with 0/1 mask: |
| 62 |
|
|
P(isnan(P))=0; |
| 63 |
|
|
mskW=mygrid.mskW; mskW(isnan(mskW))=0; |
| 64 |
|
|
mskS=mygrid.mskS; mskS(isnan(mskS))=0; |
| 65 |
|
|
%%weight average by portion that is filled |
| 66 |
|
|
%mskW=mygrid.hFacW; |
| 67 |
|
|
%mskS=mygrid.hFacS; |
| 68 |
|
|
|
| 69 |
|
|
%mask out near-equator values: |
| 70 |
|
|
tmp1=1+0*mygrid.YC; |
| 71 |
|
|
tmp1(abs(mygrid.YC)<10)=NaN; |
| 72 |
|
|
P=P.*repmat(tmp1,[1 1 nz]); |
| 73 |
|
|
|
| 74 |
|
|
%main computational loop: |
| 75 |
|
|
Ug=zeros(P,nz); |
| 76 |
|
|
Vg=zeros(P,nz); |
| 77 |
|
|
|
| 78 |
|
|
for iz=1:nz |
| 79 |
|
|
[dpdx,dpdy]=calc_T_grad(squeeze(P(:,:,iz)), 0); |
| 80 |
|
|
dpdx=dpdx.*mskW(:,:,iz); |
| 81 |
|
|
dpdy=dpdy.*mskS(:,:,iz); |
| 82 |
|
|
[dpdx, dpdy]=exch_UV_N(dpdx, dpdy); %add extra points |
| 83 |
|
|
|
| 84 |
|
|
[mask_u, mask_v]=exch_UV_N(squeeze(mskW(:,:,iz)), squeeze(mskS(:,:,iz))); |
| 85 |
|
|
mask_u=abs(mask_u); mask_v=abs(mask_v); %mask is always positive |
| 86 |
|
|
|
| 87 |
|
|
%average up to 4 points to get value at correct location |
| 88 |
|
|
ucur=Ug(:,:,iz); |
| 89 |
|
|
vcur=Vg(:,:,iz); |
| 90 |
|
|
|
| 91 |
|
|
for iF=1:mygrid.nFaces |
| 92 |
|
|
[nx,ny]=size(mygrid.XC{iF}); |
| 93 |
|
|
[nx2,ny2]=size(mask_u{iF}); |
| 94 |
|
|
|
| 95 |
|
|
cur_dpdx=dpdx{iF}(2:nx2,1:ny2-1); |
| 96 |
|
|
cur_dpdy=dpdy{iF}(1:nx2-1,2:ny2); |
| 97 |
|
|
cur_masku=mask_u{iF}(2:nx2,1:ny2-1); |
| 98 |
|
|
cur_maskv=mask_v{iF}(1:nx2-1,2:ny2); |
| 99 |
|
|
|
| 100 |
|
|
f=(fmat{iF}(1:nx,1:ny)+fmat{iF}(2:nx+1,1:ny))/2; |
| 101 |
|
|
npts=cur_masku(1:nx,1:ny)+cur_masku(2:nx+1,1:ny)+cur_masku(1:nx,2:ny+1)+cur_masku(2:nx+1,2:ny+1); |
| 102 |
|
|
ii=find(npts==0); |
| 103 |
|
|
npts(ii)=NaN; |
| 104 |
|
|
vcur{iF}=(cur_dpdx(1:nx,1:ny)+cur_dpdx(2:nx+1,1:ny)+cur_dpdx(1:nx,2:ny+1)+cur_dpdx(2:nx+1,2:ny+1))./(f.*npts); |
| 105 |
|
|
vcur{iF}(ii)=NaN; |
| 106 |
|
|
|
| 107 |
|
|
f=(fmat{iF}(1:nx,1:ny)+fmat{iF}(1:nx,2:ny+1))/2; |
| 108 |
|
|
npts=cur_maskv(1:nx,1:ny)+cur_maskv(2:nx+1,1:ny)+cur_maskv(1:nx,2:ny+1)+cur_maskv(2:nx+1,2:ny+1); |
| 109 |
|
|
ii=find(npts==0); |
| 110 |
|
|
npts(ii)=NaN; |
| 111 |
|
|
ucur{iF}=-1*(cur_dpdy(1:nx,1:ny)+cur_dpdy(2:nx+1,1:ny)+cur_dpdy(1:nx,2:ny+1)+cur_dpdy(2:nx+1,2:ny+1))./(f.*npts); |
| 112 |
|
|
ucur{iF}(ii)=NaN; |
| 113 |
|
|
end; |
| 114 |
|
|
|
| 115 |
|
|
Ug(:,:,iz)=ucur; |
| 116 |
|
|
Vg(:,:,iz)=vcur; |
| 117 |
|
|
end |
| 118 |
|
|
|
| 119 |
|
|
mskW=mygrid.mskW; ii=find(isnan(mskW)); mskW(ii)=0; |
| 120 |
|
|
mskS=mygrid.mskS; ii=find(isnan(mskS)); mskS(ii)=0; |
| 121 |
|
|
Vg=Vg.*mskS; |
| 122 |
|
|
Ug=Ug.*mskW; |
| 123 |
|
|
|
| 124 |
|
|
|
| 125 |
|
|
|
Parent Directory
| 
Revision Log
| 
Revision Graph