| 1 |
% |
| 2 |
% [JFz] = compute_JFz(SNAPSHOT) |
| 3 |
% |
| 4 |
% Here we compute the PV flux due to frictionnal forces as |
| 5 |
% JFz = ( TAUx * dSIGMATHETA/dy - TAUy * dSIGMATHETA/dx ) / RHO / EKL |
| 6 |
% |
| 7 |
% where: |
| 8 |
% TAU is the surface wind-stress (N/m2) |
| 9 |
% SIGMATHETA is the potential density (kg/m3) |
| 10 |
% RHO is the density (kg/m3) |
| 11 |
% EKL is the Ekman layer depth (m, positive) |
| 12 |
% |
| 13 |
% Files names are: |
| 14 |
% INPUT: |
| 15 |
% ./netcdf-files/<SNAPSHOT>/<netcdf_SIGMATHETA>.<netcdf_domain>.<netcdf_suff> |
| 16 |
% ./netcdf-files/<SNAPSHOT>/<netcdf_TAUX>.<netcdf_domain>.<netcdf_suff> |
| 17 |
% ./netcdf-files/<SNAPSHOT>/<netcdf_TAUY>.<netcdf_domain>.<netcdf_suff> |
| 18 |
% ./netcdf-files/<SNAPSHOT>/<netcdf_RHO>.<netcdf_domain>.<netcdf_suff> |
| 19 |
% ./netcdf-files/<SNAPSHOT>/<netcdf_EKL>.<netcdf_domain>.<netcdf_suff> |
| 20 |
% OUTPUT: |
| 21 |
% ./netcdf-files/<SNAPSHOT>/JFz.<netcdf_domain>.<netcdf_suff> |
| 22 |
% |
| 23 |
% with netcdf_* as global variables |
| 24 |
% |
| 25 |
% 06/27/06 |
| 26 |
% gmaze@mit.edu |
| 27 |
|
| 28 |
function varargout = compute_JFz(snapshot) |
| 29 |
|
| 30 |
global sla toshow |
| 31 |
global netcdf_suff netcdf_domain |
| 32 |
global netcdf_TAUX netcdf_TAUY netcdf_SIGMATHETA netcdf_EKL netcdf_RHO |
| 33 |
pv_checkpath |
| 34 |
|
| 35 |
|
| 36 |
% NETCDF file name: |
| 37 |
filST = netcdf_SIGMATHETA; |
| 38 |
filTx = netcdf_TAUX; |
| 39 |
filTy = netcdf_TAUY; |
| 40 |
filRHO = netcdf_RHO; |
| 41 |
filH = netcdf_EKL; |
| 42 |
|
| 43 |
% Path and extension to find them: |
| 44 |
pathname = strcat('netcdf-files',sla); |
| 45 |
ext = netcdf_suff; |
| 46 |
|
| 47 |
% Load files: |
| 48 |
ferfile = strcat(pathname,sla,snapshot,sla,filST,'.',netcdf_domain,'.',ext); |
| 49 |
ncST = netcdf(ferfile,'nowrite'); |
| 50 |
[STlon STlat STdpt] = coordfromnc(ncST); |
| 51 |
|
| 52 |
ferfile = strcat(pathname,sla,snapshot,sla,filTx,'.',netcdf_domain,'.',ext); |
| 53 |
ncTx = netcdf(ferfile,'nowrite'); |
| 54 |
ferfile = strcat(pathname,sla,snapshot,sla,filTy,'.',netcdf_domain,'.',ext); |
| 55 |
ncTy = netcdf(ferfile,'nowrite'); |
| 56 |
|
| 57 |
ferfile = strcat(pathname,sla,snapshot,sla,filRHO,'.',netcdf_domain,'.',ext); |
| 58 |
ncRHO = netcdf(ferfile,'nowrite'); |
| 59 |
RHO = ncRHO{4}(1,:,:); |
| 60 |
|
| 61 |
ferfile = strcat(pathname,sla,snapshot,sla,filH,'.',netcdf_domain,'.',ext); |
| 62 |
ncH = netcdf(ferfile,'nowrite'); |
| 63 |
EKL = ncH{4}(1,:,:); |
| 64 |
|
| 65 |
|
| 66 |
|
| 67 |
|
| 68 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 69 |
% First term |
| 70 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 71 |
|
| 72 |
% Dim: |
| 73 |
if toshow, disp('dim'), end |
| 74 |
nx = length(STlon) ; |
| 75 |
ny = length(STlat) - 1 ; |
| 76 |
nz = length(STdpt); |
| 77 |
|
| 78 |
% Pre-allocate: |
| 79 |
if toshow, disp('pre-allocate'), end |
| 80 |
dSIGMATHETAdy = zeros(nz,ny-1,nx).*NaN; |
| 81 |
dy = zeros(1,ny).*NaN; |
| 82 |
STup = zeros(nz,ny); |
| 83 |
STdw = zeros(nz,ny); |
| 84 |
|
| 85 |
% Meridional gradient of SIGMATHETA: |
| 86 |
if toshow, disp('grad'), end |
| 87 |
% Assuming the grid is regular, dy is independent of x: |
| 88 |
[dy b] = meshgrid( m_lldist([1 1]*STlon(1),STlat(1:ny+1) ), STdpt ) ; clear b |
| 89 |
for ix = 1 : nx |
| 90 |
if toshow, disp(strcat(num2str(ix),'/',num2str(nx))), end |
| 91 |
STup = squeeze(ncST{4}(:,2:ny+1,ix)); |
| 92 |
STdw = squeeze(ncST{4}(:,1:ny,ix)); |
| 93 |
dSTdy = ( STup - STdw ) ./ dy; |
| 94 |
% Change horizontal grid point definition to fit with SIGMATHETA: |
| 95 |
dSTdy = ( dSTdy(:,1:ny-1) + dSTdy(:,2:ny) )./2; |
| 96 |
dSIGMATHETAdy(:,:,ix) = dSTdy; |
| 97 |
end %for iy |
| 98 |
|
| 99 |
% Make TAUx having same limits: |
| 100 |
TAUx = ncTx{4}(1,2:ny,:); |
| 101 |
|
| 102 |
% Compute first term: TAUx * dSIGMATHETA/dy |
| 103 |
iz = 1; |
| 104 |
JFz_a = TAUx .* squeeze(dSIGMATHETAdy(iz,:,:)) ; |
| 105 |
|
| 106 |
|
| 107 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 108 |
% Second term |
| 109 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 110 |
|
| 111 |
% Dim: |
| 112 |
if toshow, disp('dim'), end |
| 113 |
nx = length(STlon) - 1; |
| 114 |
ny = length(STlat) ; |
| 115 |
nz = length(STdpt) ; |
| 116 |
|
| 117 |
% Pre-allocate: |
| 118 |
if toshow, disp('pre-allocate'), end |
| 119 |
dSIGMATHETAdx = zeros(nz,ny,nx-1).*NaN; |
| 120 |
dx = zeros(1,nx).*NaN; |
| 121 |
STup = zeros(nz,nx); |
| 122 |
STdw = zeros(nz,nx); |
| 123 |
|
| 124 |
% Zonal gradient of SIGMATHETA |
| 125 |
if toshow, disp('grad'), end |
| 126 |
for iy = 1 : ny |
| 127 |
if toshow, disp(strcat(num2str(iy),'/',num2str(ny))), end |
| 128 |
[dx b] = meshgrid( m_lldist(STlon(1:nx+1),[1 1]*STlat(iy)), STdpt ) ; clear b |
| 129 |
STup = squeeze(ncST{4}(:,iy,2:nx+1)); |
| 130 |
STdw = squeeze(ncST{4}(:,iy,1:nx)); |
| 131 |
dSTdx = ( STup - STdw ) ./ dx; |
| 132 |
% Change horizontal grid point definition to fit with SIGMATHETA: |
| 133 |
dSTdx = ( dSTdx(:,1:nx-1) + dSTdx(:,2:nx) )./2; |
| 134 |
dSIGMATHETAdx(:,iy,:) = dSTdx; |
| 135 |
end %for iy |
| 136 |
|
| 137 |
% Make TAUy having same limits: |
| 138 |
TAUy = ncTy{4}(1,:,2:nx); |
| 139 |
|
| 140 |
% Compute second term: TAUy * dSIGMATHETA/dx |
| 141 |
iz = 1; |
| 142 |
JFz_b = TAUy .* squeeze(dSIGMATHETAdx(iz,:,:)) ; |
| 143 |
|
| 144 |
|
| 145 |
|
| 146 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 147 |
% Finish ... |
| 148 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 149 |
% Then make all terms having same limits: |
| 150 |
nx = length(STlon) ; |
| 151 |
ny = length(STlat) ; |
| 152 |
nz = length(STdpt) ; |
| 153 |
JFz_a = squeeze(JFz_a(:,2:nx-1)); |
| 154 |
JFz_b = squeeze(JFz_b(2:ny-1,:)); |
| 155 |
delta_e = squeeze(EKL(2:ny-1,2:nx-1)); |
| 156 |
rho = squeeze(RHO(2:ny-1,2:nx-1)); |
| 157 |
|
| 158 |
% and finish: |
| 159 |
JFz = (JFz_a - JFz_b)./delta_e./rho; |
| 160 |
|
| 161 |
|
| 162 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 163 |
% Record |
| 164 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 165 |
if toshow, disp('record'), end |
| 166 |
|
| 167 |
% General informations: |
| 168 |
netfil = 'JFz'; |
| 169 |
units = 'kg/m3/s2'; |
| 170 |
ncid = 'JFz'; |
| 171 |
longname = 'Vertical PV flux due to frictional forces'; |
| 172 |
uniquename = 'JFz'; |
| 173 |
|
| 174 |
% Open output file: |
| 175 |
nc = netcdf(strcat(pathname,sla,snapshot,sla,netfil,'.',netcdf_domain,'.',ext),'clobber'); |
| 176 |
|
| 177 |
% Define axis: |
| 178 |
nx = length(STlon) ; |
| 179 |
ny = length(STlat) ; |
| 180 |
nz = 1 ; |
| 181 |
|
| 182 |
nc('X') = nx-2; |
| 183 |
nc('Y') = ny-2; |
| 184 |
nc('Z') = nz; |
| 185 |
|
| 186 |
nc{'X'} = ncfloat('X'); |
| 187 |
nc{'X'}.uniquename = ncchar('X'); |
| 188 |
nc{'X'}.long_name = ncchar('longitude'); |
| 189 |
nc{'X'}.gridtype = nclong(0); |
| 190 |
nc{'X'}.units = ncchar('degrees_east'); |
| 191 |
nc{'X'}(:) = STlon(2:nx-1); |
| 192 |
|
| 193 |
nc{'Y'} = ncfloat('Y'); |
| 194 |
nc{'Y'}.uniquename = ncchar('Y'); |
| 195 |
nc{'Y'}.long_name = ncchar('latitude'); |
| 196 |
nc{'Y'}.gridtype = nclong(0); |
| 197 |
nc{'Y'}.units = ncchar('degrees_north'); |
| 198 |
nc{'Y'}(:) = STlat(2:ny-1); |
| 199 |
|
| 200 |
nc{'Z'} = ncfloat('Z'); |
| 201 |
nc{'Z'}.uniquename = ncchar('Z'); |
| 202 |
nc{'Z'}.long_name = ncchar('depth'); |
| 203 |
nc{'Z'}.gridtype = nclong(0); |
| 204 |
nc{'Z'}.units = ncchar('m'); |
| 205 |
nc{'Z'}(:) = STdpt(1); |
| 206 |
|
| 207 |
% And main field: |
| 208 |
nc{ncid} = ncfloat('Z', 'Y', 'X'); |
| 209 |
nc{ncid}.units = ncchar(units); |
| 210 |
nc{ncid}.missing_value = ncfloat(NaN); |
| 211 |
nc{ncid}.FillValue_ = ncfloat(NaN); |
| 212 |
nc{ncid}.longname = ncchar(longname); |
| 213 |
nc{ncid}.uniquename = ncchar(uniquename); |
| 214 |
nc{ncid}(:,:,:) = JFz; |
| 215 |
|
| 216 |
nc=close(nc); |
| 217 |
|
| 218 |
|
| 219 |
% Output: |
| 220 |
output = struct('JFz',JFz,'lat',STlat(2:ny-1),'lon',STlon(2:nx-1)); |
| 221 |
switch nargout |
| 222 |
case 1 |
| 223 |
varargout(1) = {output}; |
| 224 |
end |
Parent Directory
| 
Revision Log
| 
Revision Graph