| 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