Warning, /utils/matlab/cs_grid/calcHeatTransFluxes.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 4cf7d6e4 on 2006-06-29 18:52:07 UTC
d43e8328b9 Dani*0001 function HT = calcHeatTransFluxes(varargin);
32ec7e18d0 Dani*0002
d43e8328b9 Dani*0003 % HT = calcHeatTransFluxes(od,ad,sd,og,ag,time,blkFile,[optional]);
0004 % HT = calcHeatTransFluxes(od,ad,sd,og,ag,time,blkFile,'grav',9.81);
32ec7e18d0 Dani*0005 %
0006 % Input:
0007 % The incoming data must be in a structured array format (which is the
d43e8328b9 Dani*0008 % format that comes from rdmnc) having the following fields:
32ec7e18d0 Dani*0009 % od [Ocean tave data] T,Ttave
0010 % ad [Atmos tave data] T,Ttave
0011 % sd [Aim tave data ] T,TSRtave,OLRtave,SSRtave,SLRtave,
d43e8328b9 Dani*0012 % EVAPtave,SHFtave,EnFxPrtave
32ec7e18d0 Dani*0013 % og [Ocean grid data] drF,HFacC,rA
0014 % ag [Atmos grid data] drF,HFacC
d43e8328b9 Dani*0015 % Other input parameters:
0016 % time (vec) Time levels to analyze ([] for all, miminum 3)
0017 % blkFile (str) Broken line file (eg 'isoLat_cs32_59.mat')
0018 % Optional parameters:
0019 % 'grav' (num, default 9.81) Acceleration due to gravity
0020 % 'LhVap' (num, default 2501) Latent heat of vaporization
0021 % 'CpO' (num, default 3994) Specific heat capacity of water
0022 % 'RhoO' (num, default 1035) Density of sea water
0023 % 'CpA' (num, default 1004) Specific heat capacity of water
32ec7e18d0 Dani*0024 %
0025 % Output:
0026 % HT is a structured array with the oceanic, atmpspheric, and total heat
0027 % transport at time(2:end-1) in PW, the calculated adjustments, heat
0028 % content (and its change) in the bands, and appropriate axis
0029 % information.
0030 %
d43e8328b9 Dani*0031 % Description:
0032 % Calculate heat transport from fluxes. Using a heat content budget
0033 % equation that equates the change in heat content of a zonal band with
0034 % the fluxes in and out of the top and bottom, find the poleward heat
0035 % transport by integrating from the north pole (HT = 0) southward. The
0036 % adjustment to the flux fields required for the heat transport at the
0037 % south pole to be 0 is also tabulated. In the atmosphere this
0038 % adjustment is associated with the frictional damping not coded to
0039 % increase the temperature.
0040 %
0041 % This function is built to be used with cubed-sphere output from the
0042 % Aim-ocean coupled model of the MITgcm. The zonal bands are defined by
0043 % a set of broken lines, the information for which is contain in the
0044 % broken line file "blkFile".
0045 %
0046 % Original Author: Daniel Enderton
0047 %
32ec7e18d0 Dani*0048 % Usage:
0049 %
0050 % >> og = rdmnc('/GridOcn/grid.*');
0051 % >> ag = rdmnc('/GridAtm/grid.*');
0052 % >> od = rdmnc('/DataOcn/tave.*');
0053 % >> ad = rdmnc('/DataAtm/tave.*');
0054 % >> sd = rdmnc('/DataAtm/aim_tave.*');
0055 %
0056 % >> LhVap = 2500;
0057 % >> aHeat = 1004./9.81; % aCp/g
0058 % >> oHeat = 3994.*1035; % oCp*oRho
0059 %
0060 % >> HT = CalcHT_Fluxes(od,ad,sd,og,ag,[],LhVap,aHeat,oHeat,blkFile)
0061 %
0062 % >> plot(HT.ylatHT,HT.HTocn(:,1),'g',...
0063 % HT.ylatHT,HT.HTatm(:,1),'b',...
0064 % HT.ylatHT,HT.HTtot(:,1),'r'); grid on;
0065
4cf7d6e47a Dani*0066 StorageAtm = 1;
0067 StorageOcn = 1;
0068
d43e8328b9 Dani*0069 % Default constants (can be overriden).
0070 LhVap = 2501;
0071 grav = 9.81;
0072 CpO = 3994;
0073 RhoO = 1035;
0074 CpA = 1004;
0075
0076 % Read input parameters.
0077 od = varargin{1};
0078 ad = varargin{2};
0079 sd = varargin{3};
0080 og = varargin{4};
0081 ag = varargin{5};
0082 time = varargin{6};
0083 blkFile = varargin{7};
0084 for ii = 8:2:length(varargin)
0085 temp = varargin{ii+1}; eval([varargin{ii},' = temp;']);
0086 end
0087
0088 aHeat = CpA./grav;
0089 oHeat = RhoO.*CpO;
0090
32ec7e18d0 Dani*0091
0092 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0093 % Prepare / reform incoming data %
0094 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0095
0096 % Determine time indecies.
0097 if isequal(od.T,ad.T) & isequal(od.T,sd.T)
0098 if isempty(time), time = od.T; i_time = 1:length(time);
0099 else [dump,i_time] = ismember(time,od.T); end
0100 else
0101 error('Time axes must match for ocean, atmospheric, and aim data!');
0102 end
0103
0104 % Grid information.
0105 nt = length(time);
0106 nc = size(og.HFacC,2);
0107 onr = size(og.HFacC,3);
0108 anr = size(ag.HFacC,3);
0109
0110 odrf = og.drF;
0111 adrf = ag.drF;
0112 ac = reshape(og.rA ,[6*nc*nc,1]);
0113 ohc = reshape(og.HFacC(1:6*nc,1:nc,1:onr),[6*nc*nc,onr]);
0114 ahc = reshape(ag.HFacC(1:6*nc,1:nc,1:anr),[6*nc*nc,anr]);
0115
0116 % Temperature fields.
4cf7d6e47a Dani*0117 if StorageOcn
0118 ot = reshape(od.Ttave(1:6*nc,1:nc,1:onr,i_time),[6*nc*nc,onr,nt]);
0119 end
0120 if StorageAtm
0121 at = reshape(ad.Ttave(1:6*nc,1:nc,1:anr,i_time),[6*nc*nc,anr,nt]);
0122 end
32ec7e18d0 Dani*0123
0124 % Aim fields.
0125 toaSW = reshape( sd.TSRtave(1:6*nc,1:nc,i_time),[6*nc*nc,nt]);
0126 toaLW = reshape( sd.OLRtave(1:6*nc,1:nc,i_time),[6*nc*nc,nt]);
0127 srfSW = reshape( sd.SSRtave(1:6*nc,1:nc,i_time),[6*nc*nc,nt]);
0128 srfLW = reshape( sd.SLRtave(1:6*nc,1:nc,i_time),[6*nc*nc,nt]);
0129 srfEvp = reshape( sd.EVAPtave(1:6*nc,1:nc,i_time),[6*nc*nc,nt]);
0130 srfSH = reshape( sd.SHFtave(1:6*nc,1:nc,i_time),[6*nc*nc,nt]);
0131 srfEFP = reshape(sd.EnFxPrtave(1:6*nc,1:nc,i_time),[6*nc*nc,nt]);
0132
0133 % Load broken line information. Determine latitude axis of a broken line
0134 % by the mean of the latitude of the edges along the broken line. These
0135 % are the latitude coordinates for the poleward heat transport calculations
0136 % and the surface flux / heat contents values are the mean between two
0137 % broken lines.
0138 load(blkFile);
0139 ydim=length(bkl_Ylat);
0140 YlatAv=sum(bkl_Ysg,1)./(1+bkl_Npts');
0141 ylatHT = [-90,YlatAv,90];
0142 ylatSF = ( ylatHT(2:ydim+2) + ylatHT(1:ydim+1) )./2;
0143
0144
0145 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0146 % Make heat transport calculations %
0147 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0148
0149 % Top of atmosphere and surface fluxes.
0150 Ft = toaSW-toaLW;
fbb955f500 Dani*0151 Fs = srfSW-srfLW-LhVap.*srfEvp-srfSH+srfEFP;
32ec7e18d0 Dani*0152
0153 % Loop over bands and calculate (per band) (1) surface area, (2) top of
0154 % atmpsphere and surface fluxes, and (3) oceanic and atmospheric heat
0155 % content. There are nbkl broken lines and nbkl+1 bands between broken
0156 % lines. The variable "bkl_Zon" gives the band number (nbkl+1 total) for a
0157 % given index between 0 and nbkl, that is, nbkl+1 total zones.
0158 AreaZon = NaN.*zeros([ydim+1,1]);
0159 Ft_Zon = NaN.*zeros([ydim+1,nt]);
0160 Fs_Zon = NaN.*zeros([ydim+1,nt]);
0161 HA = NaN.*zeros([ydim+1,nt]);
0162 HO = NaN.*zeros([ydim+1,nt]);
0163 dHAdt = NaN.*zeros([ydim+1,nt-2]);
0164 dHOdt = NaN.*zeros([ydim+1,nt-2]);
0165 for izon = 1:ydim+1
0166 ii = find(bkl_Zon == izon-1);
0167 AreaZon(izon,1) = sum(ac(ii));
0168 for it = 1:length(time)
0169 % ??? HFacC for surface fluxes? How to deal with runoff / land?
0170 Ft_Zon(izon,it) = sum(Ft(ii,it).*ac(ii));
0171 Fs_Zon(izon,it) = sum(Fs(ii,it).*ac(ii));
4cf7d6e47a Dani*0172 if StorageAtm
0173 HA(izon,it) = sum(aHeat.*((at(ii,:,it).*ahc(ii,:))*adrf).*ac(ii));
0174 end
32ec7e18d0 Dani*0175 HO(izon,it) = sum(oHeat.*((ot(ii,:,it).*ohc(ii,:))*odrf).*ac(ii));
0176 end
0177 end
0178
0179 % Calculate change in heat content (this means that the heat transport is
0180 % calculated for time(2:end-1) only, hence subsetting fluxes as well).
0181 for it = 2:length(time)-1
4cf7d6e47a Dani*0182 if StorageAtm
0183 dHAdt(:,it-1) = (HA(:,it+1)-HA(:,it-1))./(time(it+1)-time(it-1));
0184 end
32ec7e18d0 Dani*0185 dHOdt(:,it-1) = (HO(:,it+1)-HO(:,it-1))./(time(it+1)-time(it-1));
0186 end
0187 Ft_Zon = Ft_Zon(:,2:end-1);
0188 Fs_Zon = Fs_Zon(:,2:end-1);
0189
0190 % Compute offsets.
0191 % ??? However small, why is there an oceanic offset?
0192 ocnadj_Zon = Fs_Zon - dHOdt;
0193 atmadj_Zon = Ft_Zon - Fs_Zon - dHAdt;
0194 for it = 1:size(ocnadj_Zon,2)
0195 ocnadj(it) = sum(ocnadj_Zon(:,it))./sum(AreaZon);
0196 atmadj(it) = sum(atmadj_Zon(:,it))./sum(AreaZon);
0197 end
0198
0199 % Using a heat content budget equation that equates the change in heat
0200 % content of a band with the fluxes in and out of the top and bottom, find
0201 % the poleward heat transport by integrating from the north pole (HT = 0)
0202 % southward. The above computed adjustment is required for the heat
0203 % transport at the south pole to be 0.
0204 HTocn = zeros(size(Ft_Zon)+[1,0]);
0205 HTatm = zeros(size(Ft_Zon)+[1,0]);
0206 for ilat = ydim+2:-1:2
0207 HTocn(ilat-1,:) = HTocn(ilat,:) ...
0208 - Fs_Zon(ilat-1,:) ...
0209 + dHOdt(ilat-1,:) ...
0210 + ocnadj.*AreaZon(ilat-1);
0211 HTatm(ilat-1,:) = HTatm(ilat,:) ...
0212 - Ft_Zon(ilat-1,:) ...
0213 + Fs_Zon(ilat-1,:) ...
0214 + dHAdt(ilat-1,:) ...
0215 + atmadj.*AreaZon(ilat-1);
0216 end
0217 HTtot = HTocn + HTatm;
0218
0219
0220 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0221 % Assign outputs %
0222 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0223
0224 HT.time = time(2:end-1);
0225 HT.ylatSF = ylatSF;
0226 HT.ylatHT = ylatHT;
0227 HT.AreaZon = AreaZon;
0228 HT.Fs = Fs_Zon;
0229 HT.Ft = Ft_Zon;
0230 HT.HTtot = HTtot*1e-15;
0231 HT.HTocn = HTocn*1e-15;
0232 HT.HTatm = HTatm*1e-15;
0233 HT.HO = HO(:,2:end-1);
0234 HT.HA = HA(:,2:end-1);
0235 HT.dHOdt = dHOdt;
0236 HT.dHAdt = dHAdt;
0237 HT.ocnadj = ocnadj;
0238 HT.atmadj = atmadj;
