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

Annotation of /MITgcm_contrib/gmaze_pv/compute_alpha.m

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


Revision 1.1 - (hide annotations) (download)
Mon Jul 10 15:09:01 2006 UTC (17 years, 10 months ago) by gmaze
Branch: MAIN
add PVflux diag and update

1 gmaze 1.1 %
2     % [] = compute_alpha(SNAPSHOT)
3     %
4     % This function computes the thermal expansion coefficient from
5     % files of potential temperature THETA and salinity anomaly
6     % SALTanom.
7     %
8     % Files name are:
9     % INPUT:
10     % ./netcdf-files/<SNAPSHOT>/<netcdf_THETA>.<netcdf_domain>.<netcdf_suff>
11     % ./netcdf-files/<SNAPSHOT>/<netcdf_SALTanom>.<netcdf_domain>.<netcdf_suff>
12     % OUTPUT:
13     % ./netcdf-files/<SNAPSHOT>/ALPHA.<netcdf_domain>.<netcdf_suff>
14     %
15     % with: netcdf_* as global variables
16     %
17     % Alpha is computed with the subroutine sw_alpha from package SEAWATER
18     %
19     % 06/27/06
20     % gmaze@mit.edu
21    
22     function varargout = compute_alpha(snapshot)
23    
24     global sla toshow
25     global netcdf_suff netcdf_domain
26     global netcdf_SALTanom netcdf_THETA
27     pv_checkpath
28    
29    
30     % Path and extension to find netcdf-files:
31     pathname = strcat('netcdf-files',sla);
32     ext = netcdf_suff;
33    
34     % Load files:
35     ferfile = strcat(pathname,sla,snapshot,sla,netcdf_THETA,'.',netcdf_domain,'.',ext);
36     ncT = netcdf(ferfile,'nowrite');
37     [Tlon Tlat Tdpt] = coordfromnc(ncT);
38    
39     ferfile = strcat(pathname,sla,snapshot,sla,netcdf_SALTanom,'.',netcdf_domain,'.',ext);
40     ncS = netcdf(ferfile,'nowrite');
41     [Slon Slat Sdpt] = coordfromnc(ncS); % but normaly is the same grid as T
42    
43    
44    
45     %%%%%%%%%%%%%%%%%%%%%%%%%%%%
46     % surface PV flux
47     %%%%%%%%%%%%%%%%%%%%%%%%%%%%
48    
49     % Define axis:
50     nx = length(Tlon) ;
51     ny = length(Tlat) ;
52     nz = length(Tdpt) ;
53    
54    
55     % Pre-allocation:
56     if toshow,disp('Pre-allocate');end
57     ALPHA = zeros(nz,ny,nx).*NaN;
58    
59     % Compute alpha:
60     for iz = 1 : nz
61     if toshow,disp(strcat('Compute alpha for level:',num2str(iz),'/',num2str(nz)));end
62     TEMP = ncT{4}(iz,:,:);
63     SALT = ncS{4}(iz,:,:) + 35;
64     PRES = (0.09998*9.81*Tdpt(iz))*ones(ny,nx);
65     ALPHA(iz,:,:) = sw_alpha(SALT,TEMP,PRES,'ptmp');
66     end %for iz
67    
68    
69     %%%%%%%%%%%%%%%%%%%%%%%%%%%%
70     % Record
71     %%%%%%%%%%%%%%%%%%%%%%%%%%%%
72     if toshow, disp('record'), end
73    
74     % General informations:
75     netfil = 'ALPHA';
76     units = '1/K';
77     ncid = 'ALPHA';
78     longname = 'Thermal expansion coefficient';
79     uniquename = 'ALPHA';
80    
81     % Open output file:
82     nc = netcdf(strcat(pathname,sla,snapshot,sla,netfil,'.',netcdf_domain,'.',ext),'clobber');
83    
84     % Define axis:
85     nx = length(Tlon) ;
86     ny = length(Tlat) ;
87     nz = length(Tdpt) ;
88    
89     nc('X') = nx;
90     nc('Y') = ny;
91     nc('Z') = nz;
92    
93     nc{'X'} = ncfloat('X');
94     nc{'X'}.uniquename = ncchar('X');
95     nc{'X'}.long_name = ncchar('longitude');
96     nc{'X'}.gridtype = nclong(0);
97     nc{'X'}.units = ncchar('degrees_east');
98     nc{'X'}(:) = Tlon;
99    
100     nc{'Y'} = ncfloat('Y');
101     nc{'Y'}.uniquename = ncchar('Y');
102     nc{'Y'}.long_name = ncchar('latitude');
103     nc{'Y'}.gridtype = nclong(0);
104     nc{'Y'}.units = ncchar('degrees_north');
105     nc{'Y'}(:) = Tlat;
106    
107     nc{'Z'} = ncfloat('Z');
108     nc{'Z'}.uniquename = ncchar('Z');
109     nc{'Z'}.long_name = ncchar('depth');
110     nc{'Z'}.gridtype = nclong(0);
111     nc{'Z'}.units = ncchar('m');
112     nc{'Z'}(:) = Tdpt;
113    
114     % And main field:
115     nc{ncid} = ncfloat('Z', 'Y', 'X');
116     nc{ncid}.units = ncchar(units);
117     nc{ncid}.missing_value = ncfloat(NaN);
118     nc{ncid}.FillValue_ = ncfloat(NaN);
119     nc{ncid}.longname = ncchar(longname);
120     nc{ncid}.uniquename = ncchar(uniquename);
121     nc{ncid}(:,:,:) = ALPHA;
122    
123     nc=close(nc);
124    
125    

  ViewVC Help
Powered by ViewVC 1.1.22