Homework+1+Matlab+code

clear all; close all; % If you want to use Brian's codes, you need to add path where you put his codes, using addpath(''), for example, addpath('/users/ejung/MPO663/hw1/MATLAB_SOUNDINGS_FUNCS/'); %

% MMMMM Code to read the netcdf files with Matlab MMMMM % Read Netcdf File FileName = '/Users/mulatemedrano/Documents/MATLAB/Courses/MPO663/HM1/COARE_soundings_goodset.nc';

% Display the information of the file ncdump (FileName)

% Extract the variables of interest (these are the variables in the COARE data set) p = nc_varget (FileName, 'p',:); T = nc_varget (FileName, 'T',:); RH = nc_varget (FileName, 'RH',:); IRTEMP = nc_varget (FileName, 'IRTEMP',:); % MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM

( For HW: Replace this simple crude sounding with a data grab of the COARE soundings .nc file) % Play with a sounding using mse as thermo variable CRUDE SOUNDING % yi = interp1(x,Y,xi)

p1 = [1000.00, 850.000, 700.000, 500.000, 400.000, 300.000, 200.000, 100.000]; T1 = [27.6167, 18.4951, 10.2421, -4.81654, -14.9950, -29.9441, -52.9454, -82.0287]; q1 = [0.0178057,0.0114316,0.00704032,0.00311171,0.00150254,0.000455441,5.48021e-05,2.48807e-06];

% interpolate linearly in pressure p = 1000:-5:100 ; % target pressures every 5mb T = interp1(p1,T1,p); q = interp1(p1,q1,p);

This code should get you pretty far. But revisit it for your purposes.

% Get a hypsometric z profile for ent z = z_hypsometric (p, T, q);

% ent of the actual sounding ent = entropy(T,p,q);

% SATURATION ent of the sounding entsat = entropy(T,p,mixingratio(esat(T),p));

% dry version entdry = entropy(T,p,1e-7*q);

% PLOT figure(8) plot(entsat,z/1000.); title('Moist entropy version') hold on plot(ent,z/1000.) plot(entdry,z/1000.)

%%%%%%%%%%%%%%%%%%% UNDILUTE PARCEL ASCENT %%%%%%%%%%%%%%%%%%% entp = ent*0 + ent(1); qp = q*0 + q(1);

plot (entp, z/1000.,'r') hold off

%%%%%%%%%%%%%%%%%%%% Invert parcel ent %%%% [Tppp,qvppp] = invert_entropy(entp, qp, p);

%%% NOW BUOYANCY Tref = 273.15; eps = 0.622; % Wallace and Hobbs 3.60 -- use Kelvin temperature Tv = (T+Tref) .* (q + eps) ./ eps ./ (1+q); Tvppp=(Tppp+Tref) .* (qvppp+ eps) ./ eps ./ (1+qvppp);

%% Virtual buoyancy of parcel figure(9) plot(Tvppp-Tv,z/1000); title('ent based parcel T_rho excess (in K)') hold on %% WITH condensate loading %plot(Tvppp-Tv-Tvp.*(qp-qvppp),z/1000,'r') % original. just typo plot(Tvppp-Tv-Tvppp.*(qp-qvppp),z/1000,'r') hold off