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

Annotation of /MITgcm_contrib/gmaze_pv/D_compute_potential_vorticity.m

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


Revision 1.1 - (hide annotations) (download)
Thu Jun 22 18:15:24 2006 UTC (19 years, 1 month ago) by gmaze
Branch: MAIN
Update routines D

1 gmaze 1.1 %
2     % [] = D_COMPUTE_POTENTIAL_VORTICITY(SNAPSHOT,[WANTSPLPV])
3     %
4     % For a time snapshot, this program multiplies the potential
5     % vorticity computed with C_COMPUTE_POTENTIAL_VORTICITY by the
6     % coefficient: -1/RHO
7     % Optional flag WANTSPLPV is turn to 0 by default. Turn it to 1
8     % if the PV computed is the simple one (f.dSIGMATHETA/dz).
9     %
10     % 06/21/2006
11     % gmaze@mit.edu
12     %
13    
14    
15     function D_compute_potential_vorticity(snapshot,varargin)
16    
17    
18     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
19     %% Setup
20     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
21     global sla netcdf_RHO netcdf_PV netcdf_domain netcdf_suff
22     pv_checkpath
23    
24     %% Flags to choose which term to compute (by default, all):
25     FLpv3 = 1;
26     if nargin==2 % case of optional flag presents:
27     if varargin{1}(1) == 1 % Case of the simple PV:
28     FLpv3 = 0;
29     end
30     end %if
31    
32     %% PV and RHO netcdf-files:
33     filPV = strcat(netcdf_PV ,'.',netcdf_domain);
34     filRHO = strcat(netcdf_RHO,'.',netcdf_domain);
35    
36     %% Path and extension to find them:
37     pathname = strcat('netcdf-files',sla,snapshot);
38     ext = strcat('.',netcdf_suff);
39    
40     %% Load netcdf files:
41     ferfile = strcat(pathname,sla,filPV,ext);
42     ncPV = netcdf(ferfile,'nowrite');
43     [PV_lon PV_lat PV_dpt] = coordfromnc(ncPV);
44    
45     ferfile = strcat(pathname,sla,filRHO,ext);
46     ncRHO = netcdf(ferfile,'nowrite');
47     [RHO_lon RHO_lat RHO_dpt] = coordfromnc(ncRHO);
48    
49     %% Error check:
50     if PV_dpt~RHO_dpt | PV_lon~RHO_lon | PV_lat~RHO_lat
51     error('D_compute_potential_vorticity.m: RHO and PV fields must be defined on the same grid');
52     return;
53     end
54    
55     %% Flags:
56     global toshow % Turn to 1 to follow the computing process
57    
58    
59    
60     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
61     %% Apply the coefficient
62     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
63    
64     %% Pre-allocate:
65     if toshow,disp('Pre-allocate');end
66     PV = zeros(length(PV_dpt),length(PV_lat),length(PV_dpt)).*NaN;
67    
68     %% Apply:
69     if toshow,disp('Multiplying PV field by -1/RHO'),end
70     PV = - ncPV{4}(:,:,:) ./ ncRHO{4}(:,:,:) ;
71    
72    
73    
74     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
75     % Record:
76     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
77     if toshow,disp('Now reccording PV file ...'),end
78    
79     % General informations:
80     if FLpv3 == 1
81     netfil = strcat('PV','.',netcdf_domain,'.',netcdf_suff);
82     units = '1/s/m';
83     ncid = 'PV';
84     longname = 'Potential vorticity';
85     uniquename = 'potential_vorticity';
86     else
87     netfil = strcat('splPV','.',netcdf_domain,'.',netcdf_suff);
88     units = '1/s/m';
89     ncid = 'splPV';
90     longname = 'Simple Potential vorticity';
91     uniquename = 'simple_potential_vorticity';
92     end %if
93    
94     % Open output file:
95     nc = netcdf(strcat(pathname,sla,netfil),'clobber');
96    
97     % Define axis:
98     nc('X') = length(PV_lon);
99     nc('Y') = length(PV_lat);
100     nc('Z') = length(PV_dpt);
101    
102     nc{'X'} = 'X';
103     nc{'Y'} = 'Y';
104     nc{'Z'} = 'Z';
105    
106     nc{'X'} = ncfloat('X');
107     nc{'X'}.uniquename = ncchar('X');
108     nc{'X'}.long_name = ncchar('longitude');
109     nc{'X'}.gridtype = nclong(0);
110     nc{'X'}.units = ncchar('degrees_east');
111     nc{'X'}(:) = PV_lon;
112    
113     nc{'Y'} = ncfloat('Y');
114     nc{'Y'}.uniquename = ncchar('Y');
115     nc{'Y'}.long_name = ncchar('latitude');
116     nc{'Y'}.gridtype = nclong(0);
117     nc{'Y'}.units = ncchar('degrees_north');
118     nc{'Y'}(:) = PV_lat;
119    
120     nc{'Z'} = ncfloat('Z');
121     nc{'Z'}.uniquename = ncchar('Z');
122     nc{'Z'}.long_name = ncchar('depth');
123     nc{'Z'}.gridtype = nclong(0);
124     nc{'Z'}.units = ncchar('m');
125     nc{'Z'}(:) = PV_dpt;
126    
127     % And main field:
128     nc{ncid} = ncfloat('Z', 'Y', 'X');
129     nc{ncid}.units = ncchar(units);
130     nc{ncid}.missing_value = ncfloat(NaN);
131     nc{ncid}.FillValue_ = ncfloat(NaN);
132     nc{ncid}.longname = ncchar(longname);
133     nc{ncid}.uniquename = ncchar(uniquename);
134     nc{ncid}(:,:,:) = PV;
135    
136     nc=close(nc);
137    

  ViewVC Help
Powered by ViewVC 1.1.22