Dynamical properties analysis

gistfile1.matlab

classdef DynamicsAnalysis < handle
    properties (SetAccess = public, GetAccess = public)
        filename
        t_start = 300;
        t_end
        numLat
        latDeg % latitudes in degree
        pfull % pressure in mb
        refP % (1:2) don't use surface or top wind
        refLat % reference latitude for f and beta
%        refLatRange %(1:2), index of latDeg latitutde range for averaging, for tropopause height and criticality
        midlatRange %(1:2) EKE within 30% of its maximum
        
        meanZonalWind % (lat,plevel)
        meanZonalTemp
        meanEKE
        potentialTemp
        meanPs % converted to mb
        
        RossbyRadius % unit: km
        RhineScale
        
        EARTHRADIUS = 6400; % in km
        ROTATION_PERIOID = 24*3600; % in seconds
        F %unit: 1/s
        BETA %unit: 1/s/km
    end
    
    properties 
        ncid
    end
    
    methods
        function obj = DynamicsAnalysis(filename,t_start,EARTHRADIUS,ROTATION_PERIOID)
            obj.ncid = netcdf.open(filename,'NC_NOWRITE');
            if nargin > 1
                obj.t_start = t_start;
                obj.EARTHRADIUS = EARTHRADIUS;
                obj.ROTATION_PERIOID = ROTATION_PERIOID;
            end
            var_id_lat   = netcdf.inqVarID(obj.ncid,'lat');
            obj.latDeg   = netcdf.getVar(obj.ncid,var_id_lat);
            var_id_pfull = netcdf.inqVarID(obj.ncid,'pfull');
            obj.pfull    = netcdf.getVar(obj.ncid,var_id_pfull);
        end
        
        function [U,EKE] = calZonalMeanKinetics(obj)
            var_id_u = netcdf.inqVarID(obj.ncid,'ucomp');
            ucomp = netcdf.getVar(obj.ncid,var_id_u);
            numPlevels = size(ucomp,3);
            obj.numLat = size(ucomp,2);
            obj.t_end = size(ucomp,4);
            U = zeros(size(ucomp,2),size(ucomp,3));
            U3d = zeros(size(ucomp,1),size(ucomp,2),size(ucomp,3));
            
            U3d = sum(ucomp(:,:,:,obj.t_start:obj.t_end),4)/(obj.t_end-obj.t_start+1);
            U = reshape(sum(U3d,1)/size(U3d,1),size(U3d,2),size(U3d,3));
            obj.meanZonalWind = U;
            
            EKE = zeros(size(ucomp,2),size(ucomp,3));
            for i = 1:size(ucomp,1)
                for t = obj.t_start:obj.t_end
                    EKE = EKE + (reshape(ucomp(i,:,:,t),size(ucomp,2),size(ucomp,3))-U).^2;
                end
            end
            
            clear ucomp
            var_id_v = netcdf.inqVarID(obj.ncid,'vcomp');
            vcomp = netcdf.getVar(obj.ncid,var_id_v);
            for i = 1:size(vcomp,1)
                for t = obj.t_start:obj.t_end
                    EKE = EKE + reshape(vcomp(i,:,:,t),size(vcomp,2),size(vcomp,3)).^2;
                end
            end
            
            COS_WEIGHT = cos(pi*obj.latDeg(:)/180);
            COS_WEIGHT = repmat(COS_WEIGHT,1,size(EKE,2));
            EKE = EKE.*COS_WEIGHT;
            
            obj.meanEKE = EKE;
            clear vcomp
            
            obj.refP = [round(0.2*numPlevels), round(0.8*numPlevels)];
            EKE1d = reshape(sum(EKE(:,obj.refP(1):obj.refP(2)),2)./(obj.refP(2)-obj.refP(1)+1),...
                1,[]);
            obj.refLat = find(EKE1d == max(EKE1d(1:round(obj.numLat/2))));
            obj.midlatRange(1) = min(find(EKE1d(1:round(obj.numLat/2))>=0.3*EKE1d(obj.refLat)));
            obj.midlatRange(2) = max(find(EKE1d(1:round(obj.numLat/2))>=0.3*EKE1d(obj.refLat)));
            
            obj.F = abs(sin(pi/180*obj.latDeg(obj.refLat))*4*pi/obj.ROTATION_PERIOID);
            obj.BETA = cos(pi/180*obj.latDeg(obj.refLat))*4*pi/obj.ROTATION_PERIOID/obj.EARTHRADIUS;
            
            var_id_ps = netcdf.inqVarID(obj.ncid,'ps');
            ps = netcdf.getVar(obj.ncid,var_id_ps);
            timeMeanPs = reshape(sum(ps(:,:,obj.t_start:obj.t_end),3)/(obj.t_end-obj.t_start+1),...
                size(ps,1),size(ps,2));
            obj.meanPs = reshape(sum(timeMeanPs,1)/size(timeMeanPs,1),1,[])/100; % converts to mb
        end
        
        function T = calZonalMeanTemp(obj,timeMeanFile)
            if nargin == 1
            var_id_Temp = netcdf.inqVarID(obj.ncid,'temp');
            temp = netcdf.getVar(obj.ncid,var_id_Temp);
            
            temp3d = sum(temp(:,:,:,obj.t_start:obj.t_end),4)/(obj.t_end-obj.t_start+1);
            T = reshape(sum(temp3d,1)/size(temp3d,1),size(temp3d,2),size(temp3d,3));
            obj.meanZonalTemp = T;
            clear temp
            end
            
            if nargin == 2
            ncid_aveTemp = netcdf.open(timeMeanFile,'NC_NOWRITE');
            var_id_Temp = netcdf.inqVarID(ncid_aveTemp,'temp');
            temp = netcdf.getVar(ncid_aveTemp,var_id_Temp);
            netcdf.close(ncid_aveTemp);
            T = reshape(sum(temp,1)/size(temp,1),size(temp,2),size(temp,3));
            obj.meanZonalTemp = T;
            clear temp
            end
            
            pFactor = zeros(size(T));
            for i=1:size(pFactor,1)
                pFactor(i,:) = (1000./obj.pfull(:)).^(2/7);
            end
            obj.potentialTemp = T.*pFactor;
        end
        
        function P_trop = tropopauseP(obj,whichLat)
            % use WMO criterion for tropopause
            pfull_range=size(obj.pfull); % levels of pressure in modelling
            lat_point=whichLat;
            pfull = obj.pfull;
            T_p=obj.meanZonalTemp(lat_point,:)';
            d_T=T_p(1:pfull_range-2)-T_p(3:pfull_range);
            d_pfull=pfull(1:pfull_range-2)-pfull(3:pfull_range);
            dT_dz=-9.8/287*pfull(2:pfull_range-1).*d_T./d_pfull./T_p(2:pfull_range-1);
         
            n_lower=2;
            n_reverse=find(dT_dz==min(dT_dz));
            P_trop=interp1(1000*dT_dz(n_lower:n_reverse),pfull(1+n_lower:(1+n_reverse)),-2,'pchip');
        end
        
        function [dtheta_dp  Np_2] = partial_potentialT_p(obj,whichLat)
            % unit K/mb
            pfull_range = size(obj.pfull,1);
            y_p_temp=obj.potentialTemp(whichLat,round(pfull_range/2):pfull_range-2);
            x_pfull=obj.pfull(round(pfull_range/2):pfull_range-2);
            
            % linear fit to get d_theta/d_p           
            [fitresult, gof] = polyfit( x_pfull(:), y_p_temp(:), 1);
            dtheta_dp = fitresult(1);
            
            pRef = 700;
            Np_2=-287/pRef*(pRef/1000)^(2/7)*dtheta_dp;
        end
        
        function RossbyR = rossbyRadius(obj,whichLat)
            P_trop = obj.tropopauseP(whichLat);
            Ps = obj.meanPs(whichLat);
            [dtheta_dp Np_2] = obj.partial_potentialT_p(whichLat);
            RossbyR = sqrt(Np_2)*(Ps-P_trop)/obj.F/1000;
        end
        
        function dtheta_dy = partialThetaY(obj,whichLat,pLevel)
            if nargin == 2
                pLevel = round(0.75*size(obj.pfull));
            end
            theta1d = obj.potentialTemp(:,pLevel);
            dtheta_di = (1/12)*(-theta1d(whichLat+2)+8*theta1d(whichLat+1)...
                -8*theta1d(whichLat-1)+theta1d(whichLat-2));
            dtheta_dy = dtheta_di*size(obj.latDeg,1)/pi/obj.EARTHRADIUS;
        end
        
        function Sc = criticality(obj,whichLat)
            P_trop = tropopauseP(obj,whichLat);
            Ps = obj.meanPs(whichLat);
            [dtheta_dp  Np_2] = partial_potentialT_p(obj,whichLat);
            dtheta_dy = partialThetaY(obj,whichLat);
            Sc = 0.5*obj.F/obj.BETA*dtheta_dy/((Ps-P_trop)*(-dtheta_dp));
        end
        
        function delete(obj)
            netcdf.close(obj.ncid);
            clear obj
        end
    end
end