% Script to compute column integrated hydrometeor mass from WRF mixing ratio data.

% 12 Feb 2014
% John D McMillen
% john.mcmillen@utah.edu

% This script will compute column integrated hydrometeor mass from WRF 
% mixing ratio data.  The integration I will use is adapted from Grasso 
% and Greenwald 2004, published in Monthly Weather Review.

% CIRx = summation(Qx*rho*dz) summed from 0 to z top

% x here represents the specific hydrometeor
% I will calculate the dry air rho from the ideal gas law.
% I will solve for dz using z computed from geopotential at w grid points.
% This will center the vertical difference around the Qx and rho values 
% which are solved by WRF at mass points.

close all
clear all

% Constants
r_d = 287;          % Dry gas constant
cp = 1004;          % Specific heat of dry air at const pressure
g   = 9.81;         % gravity
t0  = 300;          % Base state potential temperature, K
p00 = 100000;       % Surface pressure, Pa

% Import data

% Geopotential (perturbation and base state)
PH = double(ncread('/uufs/chpc.utah.edu/common/home/steenburgh-group3/johnmc/WRF_3.4/WRF_Real/27oct2010/controlmpvar/wrfout/wrfout_d03_2010-10-27_00:00:00',...
    'PH'));

PHB = double(ncread('/uufs/chpc.utah.edu/common/home/steenburgh-group3/johnmc/WRF_3.4/WRF_Real/27oct2010/controlmpvar/wrfout/wrfout_d03_2010-10-27_00:00:00',...
    'PHB'));

% Pressure (perturbation and base state)
P = double(ncread('/uufs/chpc.utah.edu/common/home/steenburgh-group3/johnmc/WRF_3.4/WRF_Real/27oct2010/controlmpvar/wrfout/wrfout_d03_2010-10-27_00:00:00',...
    'P'));

PB = double(ncread('/uufs/chpc.utah.edu/common/home/steenburgh-group3/johnmc/WRF_3.4/WRF_Real/27oct2010/controlmpvar/wrfout/wrfout_d03_2010-10-27_00:00:00',...
    'PB'));

% Perturbation Potential Temperature
T = double(ncread('/uufs/chpc.utah.edu/common/home/steenburgh-group3/johnmc/WRF_3.4/WRF_Real/27oct2010/controlmpvar/wrfout/wrfout_d03_2010-10-27_00:00:00',...
    'T'));

% Cloud Liquid Water
Qcw = double(ncread('/uufs/chpc.utah.edu/common/home/steenburgh-group3/johnmc/WRF_3.4/WRF_Real/27oct2010/controlmpvar/wrfout/wrfout_d03_2010-10-27_00:00:00',...
    'QCLOUD'));

% Determine height in meters at each model level. 

GEO = PH + PHB;

Z = GEO./g;

% Compute dz

for i=1:(size(Z,3) - 1)
    dz(:,:,i) = Z(:,:,i+1) - Z(:,:,i);
end

% Determine the pressure at each model level.

pres = P + PB;

% Determine the potential temerature at each model level.

theta = T + t0;

% Compute the temperature at each model layer.

temp = theta./((p00./pres).^(r_d/cp));

% Compute rho.

rho = pres./(r_d.*temp);

% Calculate column integral.


for i=1:(size(Z,3) - 1)
    CIRx(:,:,i) = Qcw(:,:,i).*rho(:,:,i).*dz(:,:,i);
end

CIRx = sum(CIRx,3);

% Visualize the result
lake = double(ncread('/uufs/chpc.utah.edu/common/home/steenburgh-group3/johnmc/WRF_3.4/WRF_Real/27oct2010/controlmpvar/wrfout/wrfout_d03_2010-10-27_02:30:00'...
,'LANDMASK'));

CIRx = CIRx';
lake = lake';

figure(1)
contour(CIRx,30,'k')
hold on
contour(lake,1,'k','linewidth',2)
hold off

print -dpng -r200 intcloud.png;

