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

Annotation of /MITgcm_contrib/gmaze_pv/compute_JBz.m

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


Revision 1.5 - (hide annotations) (download)
Wed Sep 19 15:37:38 2007 UTC (16 years, 7 months ago) by gmaze
Branch: MAIN
CVS Tags: HEAD
Changes since 1.4: +0 -0 lines
General Update

1 gmaze 1.1 %
2 gmaze 1.3 % [JBz] = compute_JBz(SNAPSHOT)
3 gmaze 1.1 %
4     % Here we compute the PV flux due to diabatic processes as
5     % JFz = - alpha * f * Qnet / MLD / Cw
6     % where:
7     % alpha = 2.5*E-4 1/K is the thermal expansion coefficient
8     % f = 2*OMEGA*sin(LAT) is the Coriolis parameter
9     % Qnet is the net surface heat flux (W/m^2), positive downward
10 gmaze 1.2 % MLD is the mixed layer depth (m, positive)
11 gmaze 1.1 % Cw = 4187 J/kg/K is the specific heat of seawater
12     %
13     % Files names are:
14     % INPUT:
15     % ./netcdf-files/<SNAPSHOT>/<netcdf_Qnet>.<netcdf_domain>.<netcdf_suff>
16     % ./netcdf-files/<SNAPSHOT>/<netcdf_MLD>.<netcdf_domain>.<netcdf_suff>
17     % OUTPUT:
18     % ./netcdf-files/<SNAPSHOT>/JBz.<netcdf_domain>.<netcdf_suff>
19     %
20     % with: netcdf_* as global variables
21     %
22     % 06/27/06
23     % gmaze@mit.edu
24    
25 gmaze 1.3 function varargout = compute_JBz(snapshot)
26 gmaze 1.1
27     global sla toshow
28     global netcdf_suff netcdf_domain
29     global netcdf_Qnet netcdf_MLD
30     pv_checkpath
31    
32    
33     % Path and extension to find netcdf-files:
34     pathname = strcat('netcdf-files',sla);
35     ext = netcdf_suff;
36    
37     % Load files:
38     ferfile = strcat(pathname,sla,snapshot,sla,netcdf_Qnet,'.',netcdf_domain,'.',ext);
39     ncQ = netcdf(ferfile,'nowrite');
40     [Qlon Qlat Qdpt] = coordfromnc(ncQ);
41    
42     ferfile = strcat(pathname,sla,snapshot,sla,netcdf_MLD,'.',netcdf_domain,'.',ext);
43     ncMLD = netcdf(ferfile,'nowrite');
44     [MLDlon MLDlat MLDdpt] = coordfromnc(ncMLD);
45    
46    
47    
48     %%%%%%%%%%%%%%%%%%%%%%%%%%%%
49     % surface PV flux
50     %%%%%%%%%%%%%%%%%%%%%%%%%%%%
51    
52     % Define axis:
53     nx = length(Qlon) ;
54     ny = length(Qlat) ;
55     nz = length(Qdpt) ;
56    
57    
58     % Planetary vorticity:
59     f = 2*(2*pi/86400)*sin(Qlat*pi/180);
60     [a f] = meshgrid(Qlon,f); clear a c
61    
62    
63     % Net surface heat flux:
64     Qnet = ncQ{4}(:,:,:);
65    
66    
67     % Mixed layer Depth:
68     MLD = ncMLD{4}(:,:,:);
69    
70    
71     % Coefficient:
72     alpha = 2.5*10^(-4); % Surface average value
73     Cw = 4187;
74     coef = - alpha / Cw;
75    
76    
77     % JBz:
78     JBz = zeros(nz,ny,nx).*NaN;
79     JBz(1,:,:) = coef*f.*Qnet./MLD;
80    
81    
82    
83     %%%%%%%%%%%%%%%%%%%%%%%%%%%%
84     % Record
85     %%%%%%%%%%%%%%%%%%%%%%%%%%%%
86     if toshow, disp('record'), end
87    
88     % General informations:
89     netfil = 'JBz';
90     units = 'kg/m3/s2';
91     ncid = 'JBz';
92     longname = 'Vertical PV flux due to diabatic processes';
93     uniquename = 'JBz';
94    
95     % Open output file:
96     nc = netcdf(strcat(pathname,sla,snapshot,sla,netfil,'.',netcdf_domain,'.',ext),'clobber');
97    
98     % Define axis:
99     nx = length(Qlon) ;
100     ny = length(Qlat) ;
101     nz = 1 ;
102    
103     nc('X') = nx;
104     nc('Y') = ny;
105     nc('Z') = nz;
106    
107     nc{'X'} = ncfloat('X');
108     nc{'X'}.uniquename = ncchar('X');
109     nc{'X'}.long_name = ncchar('longitude');
110     nc{'X'}.gridtype = nclong(0);
111     nc{'X'}.units = ncchar('degrees_east');
112     nc{'X'}(:) = Qlon;
113    
114     nc{'Y'} = ncfloat('Y');
115     nc{'Y'}.uniquename = ncchar('Y');
116     nc{'Y'}.long_name = ncchar('latitude');
117     nc{'Y'}.gridtype = nclong(0);
118     nc{'Y'}.units = ncchar('degrees_north');
119     nc{'Y'}(:) = Qlat;
120    
121     nc{'Z'} = ncfloat('Z');
122     nc{'Z'}.uniquename = ncchar('Z');
123     nc{'Z'}.long_name = ncchar('depth');
124     nc{'Z'}.gridtype = nclong(0);
125     nc{'Z'}.units = ncchar('m');
126 gmaze 1.3 nc{'Z'}(:) = Qdpt(1);
127 gmaze 1.1
128     % And main field:
129     nc{ncid} = ncfloat('Z', 'Y', 'X');
130     nc{ncid}.units = ncchar(units);
131     nc{ncid}.missing_value = ncfloat(NaN);
132     nc{ncid}.FillValue_ = ncfloat(NaN);
133     nc{ncid}.longname = ncchar(longname);
134     nc{ncid}.uniquename = ncchar(uniquename);
135     nc{ncid}(:,:,:) = JBz;
136    
137     nc=close(nc);
138 gmaze 1.3 close(ncQ);
139     close(ncMLD);
140 gmaze 1.1
141    
142 gmaze 1.3
143     % Output:
144     output = struct('JBz',JBz,'lat',Qlat,'lon',Qlon);
145     switch nargout
146     case 1
147     varargout(1) = {output};
148     end

  ViewVC Help
Powered by ViewVC 1.1.22