% read_radprof.m

clear
close all

% matlab script path
path_prog = '/Users/byronb/Documents/MATLAB/2020_ATOMIC';
cd(path_prog);                  % cd to matlab directory
addpath(genpath(path_prog));	% add scripts to path
rehash;

% read all radiative profiles
dfl = '/Users/byronb/Raw_Data_Files.noindex/ATOMIC_2020/rad_profiles/all_rad_profiles.nc';
I = ncinfo(dfl);

swdn = ncread(dfl,'sw_dn');
swup = ncread(dfl,'sw_up');
lwdn = ncread(dfl,'lw_dn');
lwup = ncread(dfl,'lw_up');
plev = ncread(dfl,'plev');
plev = double(plev);
play = ncread(dfl,'play');
play = double(play);
ta = ncread(dfl,'tlay');
Q = ncread(dfl,'h2o');

% tam = nanmean(ta,2);
% qam = nanmean(Q,2);
% rho = play/Rg./tam - 0.61*qam;
% figure; plot(nanmean(Q'),play); set(gca,'ydir','reverse');
% ylim([800 1015]*100);
% figure; plot(rho,play); set(gca,'ydir','reverse');
% ylim([800 1015]*100);

swnet = swdn-swup;
lwnet = lwdn-lwup;

netradg = diff(nanmean(swnet+lwnet,2))/100; % W/m2/Pa = J*s/kg/m
netradg(1:715) = despike2(netradg(1:715));  % remove low level spikes

% tdot = g/cpa * netradg * 86400 s/d  -  in K/d units
% K/d = (m/s2) / (J/kg/K) * (J*s/kg/m) * (s/d)
g = 9.80665; % m/s2
cpa = 1001;  % J/kg/K, ~mean over air column
tdot = g/cpa*netradg*86400;

figure; plot([0,0],[play(1),play(end)],'k-',tdot,play,'.'); set(gca,'ydir','reverse');
grid; xlabel('Q_r (deg/day)'); ylim([play(end),play(1)]);
ylabel('Pressure (pa)');title('Mean ATOMIC profiles');
ppath = '/Users/byronb/Raw_Data_Files.noindex/ATOMIC_2020/plots/heating_profile.png';
print(ppath,'-dpng');

sfl = '/Users/byronb/Raw_Data_Files.noindex/ATOMIC_2020/rad_profiles/tdot.mat';
save(sfl,'tdot','play');

