--- MITgcm_contrib/dfer/matlab_stuff/calcHeatTransDirect.m	2013/03/05 18:08:06	1.1
+++ MITgcm_contrib/dfer/matlab_stuff/calcHeatTransDirect.m	2018/03/07 22:01:18	1.2
@@ -1,25 +1,23 @@
 function HT = calcHeatTransDirect(varargin)
 
-% HT = calcHeatTransDirect(d,g,time,flu,blkFile,[optional]);
-% HT = calcHeatTransDirect(d,g,time,flu,blkFile,'grav',9.81);
+% HT = calcHeatTransDirect(d,g,dE,GM,time,flu,blkFile,[mask,phi_ref]);
 %
 % Input arguements:  
 %   The incoming field data (d) and grid data (g) must be in a structured
 %   array format (which is the format that comes from rdmnc):
-%       d  [Field data]  hUtave,hVtave,uVeltave,vVeltave,Ttave,UTtave,
-%                        VTtave,(Stave,UStave,VStave for atm)
-%       g  [Grid data ]  drF,dxG,dyG,dxC,dyC,HFacW,HFacS,rA
+%       d  [Structure]    Various inputs depending on oceanic or 
+%                         atmospheric computation, rstar or not,
+%                         and GM form/ouputs
+%       g  [Structure]    drF,dxG,dyG,dxC,dyC,HFacW,HFacS,rA
+%       dE [Structure]    More field data
+%       GM [Structure]    Gent anf McWilliams data
 %   Other input parameters:
-%       time     (vec)  Time levels to analyze ([] for all)
-%       flu      (str)  'O' or 'A' for ocean or atmosphere
-%       blkFile  (str)  Broken line file (eg 'isoLat_cs32_59.mat')
+%       time [vector]     Time levels to analyze ([] for all)
+%       flu  [string]     'O' or 'A' for ocean or atmosphere
+%       blkFile [string]  Broken line file (eg 'isoLat_cs32_59.mat')
 %   Optional parameters:
-%       'grav'   (num, default 9.81)  Acceleration due to gravity
-%       'LhVap'  (num, default 2501)  Latent heat of vaporization
-%       'CpO'    (num, default 3994)  Specific heat capacity of water
-%       'RhoO'   (num, default 1030)  Density of sea water
-%       'CpA'    (num, default 1004)  Specific heat capacity of water
-%       %'DiffKh' (num, default 800)   Horizontal diffusivity
+%       mask [Structure]  W and S mask           
+%       phi_ref [vector]  Reference geopotential height 
 %
 % Output:
 %   HT_Out is a structured array with the following fields:
@@ -79,9 +77,6 @@
 nt= varargin{5};
 flu = varargin{6};
 blkFile = varargin{7};
-%for ii = 8:2:length(varargin)
-%    temp = varargin{ii+1}; eval([varargin{ii},' = temp;']);
-%end
 if length(varargin) >= 8 & ~isempty(varargin{8})
    mask = varargin{8};
    masking = 1;
@@ -530,64 +525,3 @@
     HT.PotHT = reshape(PotHT,ydim+2,nt,nout);
 end
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%       Stuff that might need to be added back in later (for basins)      %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-% Block for constructing mskG for different basins:
-% if nBas > 0,
-%     mskBc=rdda([Rac,'maskC_bas.bin'],[6*nc*nc 3],1,'real*4','b');
-%     mskBw=rdda([Rac,'maskW_bas.bin'],[6*nc*nc 3],1,'real*4','b');
-%     mskBs=rdda([Rac,'maskS_bas.bin'],[6*nc*nc 3],1,'real*4','b');
-%     if nBas==2,
-%         mskBc(:,2)=mskBc(:,2)+mskBc(:,3);
-%         mskBw(:,2)=mskBw(:,2)+mskBw(:,3);
-%         mskBs(:,2)=mskBs(:,2)+mskBs(:,3);
-%         mskBc=min(1,mskBc); mskBw=min(1,mskBw); mskBs=min(1,mskBs);
-%     end
-%     %- load: np_Sep, ij_Sep, tp_Sep:
-%     sep_lineF=[Rac,'sepBas_cs32_60'];
-%     load(sep_lineF);
-%     fprintf([' + bassin mask & Sep.line:',sep_lineF,' \n']);
-%     %- compute mask for each bassin: -----------
-%     kMsep=1;
-%     if kMsep,
-%         mskW=1+min(1,ceil(hw(:,1)));
-%         mskS=1+min(1,ceil(hs(:,1)));
-%         for b=1:nBas,
-%             bs=b; b1=1+bs; b2=2+rem(bs,nBas);
-%             if nBas == 2, bs=b+b-1; b1=2; b2=3 ; end
-%             for j=1:ydim+1,
-%                 for i=1:np_Sep(bs,j),
-%                     ij=ij_Sep(bs,j,i); typ=abs(tp_Sep(bs,j,i));
-%                     if typ == 1,
-%                         mskG(j,b1)=mskG(j,b1)*mskW(ij);
-%                         mskG(j,b2)=mskG(j,b2)*mskW(ij);
-%                     elseif typ == 2,
-%                         mskG(j,b1)=mskG(j,b1)*mskS(ij);
-%                         mskG(j,b2)=mskG(j,b2)*mskS(ij);
-%                     end
-%                 end
-%             end
-%         end
-%         mskG=2-min(2,mskG); 
-%     end
-% end
-
-% % Masking for different basin.
-% mskZ = ones(ydim+2,1+nBas);
-% mskZ([2:ydim+1],:) = mskG([1:ydim],:) + mskG([2:ydim+1],:);
-% mskZ = reshape(min(mskZ,1),(ydim+2)*(1+nBas),1);
-% I = find(mskZ == 0);
-% mskZ = reshape(mskZ,(ydim+2),1+nBas);
-% for n=1:min(6,nout)
-%     var=SenHT(:,:,n); 
-%     var=reshape(var,(ydim+2)*(1+nBas),1);  var(I) = NaN;
-%     var=reshape(var,(ydim+2),1+nBas);  SenHT(:,:,n) = var;
-% end
-
-% Makes sure that there are no zeros in area AreaZon:
-% - set AreaZon to 1 if = zero
-% AreaZon=reshape(AreaZon,(ydim+1)*(1+nBas),1);
-% [I]=find(AreaZon==0); AreaZon(I)=1; 
-% AreaZon=reshape(AreaZon,ydim+1,1+nBas);