/[MITgcm]/MITgcm_contrib/gmaze_pv/C_compute_potential_vorticity.m
ViewVC logotype

Diff of /MITgcm_contrib/gmaze_pv/C_compute_potential_vorticity.m

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

revision 1.5 by gmaze, Thu Feb 1 17:02:02 2007 UTC revision 1.6 by gmaze, Wed Sep 19 14:45:59 2007 UTC
# Line 1  Line 1 
1  %  %
2  % [Q] = C_compute_potential_vorticity(SNAPSHOT,[WANTSPLPV])  % [Q] = C_compute_potential_vorticity(SNAPSHOT,[WANTSPLPV])
3    % [Q1,Q2,Q3] = C_compute_potential_vorticity(SNAPSHOT,[WANTSPLPV])
4  %  %
5  % This file computes the potential vorticity Q from  % This file computes the potential vorticity Q from
6  % netcdf files of relative vorticity (OMEGAX, OMEGAY, ZETA)  % netcdf files of relative vorticity (OMEGAX, OMEGAY, ZETA)
# Line 49  if nargin==2  % case of optional flag pr Line 50  if nargin==2  % case of optional flag pr
50      FLpv3 = 2;      FLpv3 = 2;
51    end    end
52  end %if  end %if
53    %[FLpv1 FLpv2 FLpv3]
54    
55    
56  %% Optionnal flags:  %% Optionnal flags:
57  global toshow % Turn to 1 to follow the computing process  global toshow % Turn to 1 to follow the computing process
# Line 58  global toshow % Turn to 1 to follow the Line 61  global toshow % Turn to 1 to follow the
61    
62  % Path and extension to find them:  % Path and extension to find them:
63  pathname = strcat('netcdf-files',sla,snapshot,sla);  pathname = strcat('netcdf-files',sla,snapshot,sla);
64    %pathname = '.';
65  ext      = strcat('.',netcdf_suff);  ext      = strcat('.',netcdf_suff);
66    
67  % Names:  % Names:
# Line 245  dSIGMATHETAdz = zeros(nz-1,ny,nx).*NaN; Line 249  dSIGMATHETAdz = zeros(nz-1,ny,nx).*NaN;
249             dz = zeros(1,nz).*NaN;             dz = zeros(1,nz).*NaN;
250    
251  % Vertical grid differences:  % Vertical grid differences:
252           dz = diff(STdpt);  % STdpt contains negative values with STdpt(1) at the surface
253    % and STdpt(end) at the bottom of the ocean.
254    % So dz is positive with respect to z axis upward:
255             dz = -diff(STdpt);
256  [a dz_3D c] = meshgrid(STlat,dz,STlon); clear a c  [a dz_3D c] = meshgrid(STlat,dz,STlon); clear a c
257    
258  % Vertical gradient:  % Vertical gradient:
259  if toshow,disp('Vertical gradient of SIGMATHETA'), end  if toshow,disp('Vertical gradient of SIGMATHETA'), end
260             ST = ncST{4}(:,:,:);             ST = ncST{4}(:,:,:);
261  dSIGMATHETAdz = ( ST(2:nz+1,:,:) - ST(1:nz,:,:) ) ./ dz_3D;             % Z axis upward, so vertical derivative is upper-part
262               % minus lower-part:
263    dSIGMATHETAdz = ( ST(1:nz,:,:) - ST(2:nz+1,:,:) ) ./ dz_3D;
264  clear dz_3D ST  clear dz_3D ST
265    
266  % Change vertical gridding:  % Change vertical gridding:
# Line 403  nc{ncid}.uniquename    = ncchar(uniquena Line 412  nc{ncid}.uniquename    = ncchar(uniquena
412  nc{ncid}(:,:,:)        = PV;  nc{ncid}(:,:,:)        = PV;
413    
414  nc=close(nc);  nc=close(nc);
415  close(ncOx);  if FLpv3 ~= 2
416  close(ncOy);     close(ncOx);
417  close(ncOz);     close(ncOy);
418       close(ncOz);
419    end
420  close(ncST);  close(ncST);
421    
422  % Outputs:  % Outputs:
# Line 413  OUT = struct('PV',PV,'dpt',PV_dpt,'lat', Line 424  OUT = struct('PV',PV,'dpt',PV_dpt,'lat',
424  switch nargout  switch nargout
425   case 1   case 1
426    varargout(1) = {OUT};    varargout(1) = {OUT};
427     case 2
428      varargout(1) = {struct('PV1',PV1,'dpt',PV1_dpt,'lat',PV1_lat,'lon',PV1_lon)};
429      varargout(2) = {struct('PV2',PV2,'dpt',PV2_dpt,'lat',PV2_lat,'lon',PV2_lon)};
430     case 3
431      varargout(1) = {struct('PV1',PV1,'dpt',PV1_dpt,'lat',PV1_lat,'lon',PV1_lon)};
432      varargout(2) = {struct('PV2',PV2,'dpt',PV2_dpt,'lat',PV2_lat,'lon',PV2_lon)};
433      varargout(3) = {struct('PV3',PV3,'dpt',PV3_dpt,'lat',PV3_lat,'lon',PV3_lon)};
434  end  end
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.22