%make vertically integrated atmospheric temperature
clear all; close all;
      load C:\Users\thedhoe\Desktop\CMIP_seasonal_lap\obs\data\T.mat

       tlat=lat;
       tlon=lon;
       [TLON, TLAT]=meshgrid(tlon, tlat);
       [NT, NL, NN, N]=size(T);      
       load('C:\Users\thedhoe\Desktop\aqua_depth\figures\phase_cross_sections\OBS_mean_PS', 'lat', 'lon','PS')
       pslon=double(lon);
       pslat=double(lat);
       [PSLON, PSLAT]=meshgrid(pslon, pslat);
      
      %interpolate the surface pressure to the temperature gbrid
     IPS=interp2(PSLON, PSLAT,PS, TLON, TLAT);
     MASK=ones(NL,NN,N);
       
       
       for nn=1:NN;
           for n=1:N;
               do=find(plev >= (squeeze(IPS(nn,n)))/100); %levels that are below the surface at the given grid point
               MASK(do,nn,n)=NaN; %set these levels equal to NaN in the Mask
           end
       end

   doz=find(plev>00); %troposphere

     NL=length(plev);
     Phalf=mean([plev(1:end-1) plev(2:end)]')';
     Pdiff=[Phalf;1000]-[0;Phalf]; %pressure differenc by layer
     W=zeros(NL,1);W(doz)=Pdiff(doz); W=W/sum(W); %normalized weighting vector for mass integrals over the troposphere

      for mn=1:NT;
          for nn=1:NN;
              for n=1:N;
                  w=W.*squeeze(MASK(:,nn,n));
                  w=w/nansum(w); %normalize the weighting vector (above the topography)
                  VT(mn,nn,n) = nansum(squeeze(T(mn,:,nn,n)).*w');
              end
          end
      end