% %look at the who OA flux data from robert
% 
% clear all; close all;
% 
% LW=3;
% PW=10^(-15);
% 
% 
% 
% fil = 'C:\Users\Aaron Donohoe\Desktop\research\turbulent\whoi_OA\from_robert\dOA_ERAL_CERES.nc';   %%lh_oaflux_2000_2018.nc'
% 
% %ncdisp(fil);
% obs_lat=double(ncread(fil, 'lat')); NN=length(obs_lat);
% obs_lon=double(ncread(fil, 'lon')); N=length(obs_lon);
% OBS_LH=squeeze(mean(permute(double(ncread(fil, 'lh', [1 1 1], [Inf Inf Inf])), [3 2 1] ))); 
% OBS_SH=squeeze(mean(permute(double(ncread(fil, 'sh', [1 1 1], [Inf Inf Inf])), [3 2 1] ))); 
% 
%  load('C:\Users\Aaron Donohoe\Desktop\research\radiation\ceres_ebaf\4.0\CERES40_climatology_2020', 'lat', 'lon', 'CDT', 'CUT', 'CDS', 'CUS', 'COLR', 'CUS_clr', 'CDS_clr',  'CUT_clr', 'CLDS_clr','CLDS','CLUS', 'CLUS_clr', 'COLR_clr')
%  
% OBS_DS=squeeze(mean(CDS)); 
% OBS_US=squeeze(mean(CUS)); 
% OBS_LDS=squeeze(mean(CLDS)); 
% OBS_LUS=squeeze(mean(CLUS)); 
% OBS_DT=squeeze(mean(CDT)); 
% OBS_UT=squeeze(mean(CUT)); 
% OBS_OLR=squeeze(mean(COLR)); 
% 
% %%% make landmask
% 
% Lon=obs_lon;
% Lon(find(Lon > 180)) = Lon(find(Lon > 180))-360;  
% [LON LAT] = meshgrid(Lon,obs_lat);
% land = double(landmask(LAT,LON));
% 
% %sphcolor2jet(land,obs_lat,obs_lon,5);
% [LON LAT] = meshgrid(obs_lon, obs_lat);
% land(find(LAT<=-83)) = 1;
% 
% 
% % some misisng data at a single grid point
% OBS_SHF = OBS_DS-OBS_US-OBS_SH-OBS_LH+OBS_LDS-OBS_LUS;
% OBS_SRAD = OBS_DS-OBS_US+OBS_LDS-OBS_LUS;
% OBS_TURB = -OBS_SH-OBS_LH;
% OBS_ARAD = OBS_DT-OBS_UT-OBS_OLR-OBS_LDS+OBS_LUS-OBS_DS+OBS_US;
% 
% % % %% annual surface heat flux
% ocean=1-land;
% %  
% OBS_G_SHF = trapz(sind(obs_lat), nansum((OBS_SHF.*ocean)')/N) / trapz(sind(obs_lat), nansum((ocean)')/N); 
% OBS_G_SRAD = trapz(sind(obs_lat), nansum((OBS_SRAD.*ocean)')/N) / trapz(sind(obs_lat), nansum((ocean)')/N); 
% OBS_G_SRAD_L = trapz(sind(obs_lat), nansum((OBS_SRAD.*land)')/N) / trapz(sind(obs_lat), nansum((land)')/N);
% OBS_G_TURB = trapz(sind(obs_lat), nansum((OBS_TURB.*ocean)')/N) / trapz(sind(obs_lat), nansum((ocean)')/N); 
% OBS_G_E = trapz(sind(obs_lat), nansum((OBS_LH.*ocean)')/N) / trapz(sind(obs_lat), nansum((ocean)')/N); 
% OBS_G_SH = trapz(sind(obs_lat), nansum((OBS_SH.*ocean)')/N) / trapz(sind(obs_lat), nansum((ocean)')/N); 
% %%% ARAD over the global land and ocean
% OBS_G_ARAD = trapz(sind(obs_lat), mean(OBS_ARAD')) / trapz(sind(obs_lat), ones(size(obs_lat))); 
% 
% 
%  ocean2 = ocean;
%  ocean2(find(ocean == 0)) = NaN;
%  land2 = land;
%  land2(find(land == 0)) = NaN;
% 
% 
% 
% 
% % AND IMPLIES oht
% ae=6.36*10^6;
%  mlat=[-90; (lat(1:end-1)+lat(2:end))/2; 90];
%  dlat=diff(mlat);
%  dA = (ae*pi*dlat/180).*(2*pi*ae*cosd(lat)).*(mean(ocean')'); % area of the region centered around lat 
% dA_tot =  (ae*pi*dlat/180).*(2*pi*ae*cosd(lat));
% dA_L =  (ae*pi*dlat/180).*(2*pi*ae*cosd(lat)).*(mean(land')'); % area of the region centered around lat 
% 
% % 
% PW=10^-15;
% %  figure(3); clf;
% %  hold on; plot(mlat, PW*OBS_OHT_tot, 'b-', 'LineWidth', LW);
% % hold on; plot(mlat, PW*OBS_AHT_ARAD, 'r-', 'LineWidth', LW);
% 
% % 
% % %%%% repeat with models
% %annual and zonal mean ARAD on global domain
% load('C:\Users\Aaron Donohoe\Desktop\research\cmip6\climatologt_compare_scripts\cmip3to6andobs_commongrid.mat', 'DT', 'UT', 'OLR', 'DS', 'US', 'LDS', 'LUS');
% 
%  MOD_Z_ARAD = squeeze(mean((DT-UT-OLR-DS+US-LDS+LUS)));
%  MOD_Z_OLR = squeeze(mean((OLR)));
%  
% 
% 
% %%%%
% load('C:\Users\Aaron Donohoe\Desktop\research\cmip6\climatologt_compare_scripts\Surface_common_grid', 'clat', 'clon', 'MODS', 'CLH', 'CSH', 'CLUS','CLDS', 'CUS', 'CDS')
% load('C:\Users\Aaron Donohoe\Desktop\research\cmip6\climatologt_compare_scripts\CMIP3and5_Surface_common_grid', 'clat', 'clon', 'dirs', 'C_LH', 'C_SH', 'C_LUS','C_LDS', 'C_US', 'C_DS')
% %stack cmip3 and 5 on cmip 6
% [NM MN NN N]=size(CLH);
% [NM2 MN NN N]=size(C_LH);
% 
% CLH(NM+1:NM+NM2, 1:MN, 1:NN, 1:N)= C_LH;
% CSH(NM+1:NM+NM2, 1:MN, 1:NN, 1:N)= C_SH;
% CLUS(NM+1:NM+NM2, 1:MN, 1:NN, 1:N)= C_LUS;
% CLDS(NM+1:NM+NM2, 1:MN, 1:NN, 1:N)= C_LDS;
% CUS(NM+1:NM+NM2, 1:MN, 1:NN, 1:N)= C_US;
% CDS(NM+1:NM+NM2, 1:MN,1:NN, 1:N)= C_DS;
% 
% %stack cmip3 and 5 on cmip 6
% [NM MN NN N]=size(CLH);
% [NM2 MN NN N]=size(C_LH);
% 
% CLH(NM+1:NM+NM2, 1:MN, 1:NN, 1:N)= C_LH;
% CSH(NM+1:NM+NM2, 1:MN, 1:NN, 1:N)= C_SH;
% CLUS(NM+1:NM+NM2, 1:MN, 1:NN, 1:N)= C_LUS;
% CLDS(NM+1:NM+NM2, 1:MN, 1:NN, 1:N)= C_LDS;
% CUS(NM+1:NM+NM2, 1:MN, 1:NN, 1:N)= C_US;
% CDS(NM+1:NM+NM2, 1:MN, 1:NN, 1:N)= C_DS;
% 
% 
% 
% 
% 
% MOD_LH=squeeze(mean(squeeze(mean(CLH)))); 
% MOD_SH=squeeze(mean(squeeze(mean(CSH))));
% MOD_DS=squeeze(mean(squeeze(mean(CDS)))); 
% MOD_US=squeeze(mean(squeeze(mean(CUS)))); 
% MOD_LDS=squeeze(mean(squeeze(mean(CLDS)))); 
% MOD_LUS=squeeze(mean(squeeze(mean(CLUS))));  
% 
% 
% 
% % some misisng data at a single grid point
% MOD_SHF = MOD_DS-MOD_US-MOD_SH-MOD_LH+MOD_LDS-MOD_LUS;
% MOD_SRAD = MOD_DS-MOD_US+MOD_LDS-MOD_LUS;
% MOD_TURB = -MOD_SH-MOD_LH;
% %% annual surface heat flux
% ocean=1-land;
% 
% MOD_G_SHF = trapz(sind(obs_lat), nansum((MOD_SHF.*ocean)')/N) / trapz(sind(obs_lat), nansum((ocean)')/N) 
% MOD_G_SRAD = trapz(sind(obs_lat), nansum((MOD_SRAD.*ocean)')/N) / trapz(sind(obs_lat), nansum((ocean)')/N) ;
% MOD_G_SRAD_L = trapz(sind(obs_lat), nansum((MOD_SRAD.*land)')/N) / trapz(sind(obs_lat), nansum((land)')/N)
% MOD_G_TURB = trapz(sind(obs_lat), nansum((MOD_TURB.*ocean)')/N) / trapz(sind(obs_lat), nansum((ocean)')/N) ;
% MOD_G_E = trapz(sind(obs_lat), nansum((MOD_LH.*ocean)')/N) / trapz(sind(obs_lat), nansum((ocean)')/N) ;
% MOD_G_SH = trapz(sind(obs_lat), nansum((MOD_SH.*ocean)')/N) / trapz(sind(obs_lat), nansum((ocean)')/N); 
% MOD_G_ARAD = trapz(sind(obs_lat), MOD_Z_ARAD) / trapz(sind(obs_lat), ones(size(obs_lat)));
% 
% MOD_Z_SHF = nanmean((MOD_SHF.*ocean2)');
% MOD_Z_SRAD = nanmean((MOD_SRAD.*ocean2)');
% MOD_Z_TURB = nanmean((MOD_TURB.*ocean2)');
% MOD_Z_E = nanmean((MOD_LH.*ocean2)');
% MOD_Z_SH = nanmean((MOD_SH.*ocean2)');
% MOD_Z_SRAD_L = nanmean((MOD_SRAD.*land2)');
% 
% OBS_Z_SHF = nanmean((OBS_SHF.*ocean2)'); 
% OBS_Z_SRAD = nanmean((OBS_SRAD.*ocean2)');
% OBS_Z_TURB = nanmean((OBS_TURB.*ocean2)');
% OBS_Z_E = nanmean((OBS_LH.*ocean2)');
% OBS_Z_SH = nanmean((OBS_SH.*ocean2)');
% OBS_Z_ARAD = mean(OBS_ARAD');
% OBS_Z_SRAD_L = nanmean((OBS_SRAD.*land2)');
% 
% 
% 
% 
% 
% 
% green= [0 .5 0];
% orange=[1 .5 0];
% red=[1 0 0]; 
% black = [ 0 0 0];
% gray=.2*[1 1 1]; 
% purp = [75 0 130]/255;
% RM = 1; %remove the global mean
% FSz=36; FN='Arial';
% LWA=3;
% 
% LW=4.5;
% RM=1;
% figure(99);clf;
% hold on; plot(obs_lat, -(OBS_Z_E - RM*OBS_G_E), '-','Color', green, 'LineWidth', LW);
% hold on; plot(obs_lat, -(MOD_Z_E - RM*MOD_G_E), '--', 'Color', green,'LineWidth', LW);
% 
% %legend('Observed', 'Model Mean')
% xlabel('Latitude', 'FontSize', FSz, 'FontName', FN)
% ylabel('Evaporation (W m^-^2)', 'FontSize', FSz, 'FontName', FN)
% 
% 
%  hold on; plot(obs_lat, -(OBS_Z_SH - RM*OBS_G_SH), '-','Color', red, 'LineWidth', LW);
%  hold on; plot(obs_lat, -(MOD_Z_SH - RM*MOD_G_SH), '--', 'Color', red,'LineWidth', LW);
%  
%  hold on; plot(obs_lat, OBS_Z_SRAD - RM*OBS_G_SRAD, '-','Color', orange, 'LineWidth', LW);
%  hold on; plot(obs_lat, MOD_Z_SRAD - RM*MOD_G_SRAD, '--', 'Color', orange,'LineWidth', LW);
%  
%  hold on; plot(obs_lat, OBS_Z_SHF - RM*OBS_G_SHF, '-','Color', black, 'LineWidth', LW);
%  hold on; plot(obs_lat, MOD_Z_SHF - RM*MOD_G_SHF, '--', 'Color', black,'LineWidth', LW);
% 
% xmn=-90; xmx=90;
% ymn=-150; ymx=150;
% xtic=-80:20:80;
% ytic=ymn:50:ymx;
% hold on; plot([xmn xmx], [0 0], '-', 'Color', gray, 'LineWidth', LWA)
% title('Ocean Energy Input', 'FontSize', FSz)
% set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', LWA, 'XTick', xtic, 'YTick', ytic); grid on; box on;
% axis([xmn xmx ymn ymx]);
% 
% 
% SCL= mean(land'); %scaling to put the land zonal mean in comparable units to the ocean
% 
% figure(100);clf;
% hold on; plot(obs_lat, (OBS_Z_E) - RM*OBS_G_E, '-','Color', green, 'LineWidth', LW);
% hold on; plot(obs_lat, (MOD_Z_E) - RM*MOD_G_E, '--', 'Color', green,'LineWidth', LW);
% 
% hold on; plot(obs_lat, (OBS_Z_ARAD) - RM*OBS_G_ARAD, '-','Color', orange, 'LineWidth', LW);
% hold on; plot(obs_lat, (MOD_Z_ARAD) - RM*MOD_G_ARAD, '--', 'Color', orange,'LineWidth', LW);
% 
% hold on; plot(obs_lat, (OBS_Z_SH - RM*OBS_G_SH), '-','Color', red, 'LineWidth', LW);
% hold on; plot(obs_lat, (MOD_Z_SH - RM*MOD_G_SH), '--', 'Color', red,'LineWidth', LW);
% 
%  hold on; plot(obs_lat, SCL.*(OBS_Z_SRAD_L - RM*OBS_G_SRAD_L), '-','Color', purp, 'LineWidth', LW);
%  hold on; plot(obs_lat, SCL.*(MOD_Z_SRAD_L - RM*MOD_G_SRAD_L), '--', 'Color', purp,'LineWidth', LW);
% 
% title('Atmospheric Energy Input', 'FontSize', FSz)
% hold on; plot([xmn xmx], [0 0], '-', 'Color', gray, 'LineWidth', LWA);
% set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', LWA, 'XTick', xtic, 'YTick', ytic); grid on; box on;
% axis([xmn xmx ymn ymx]);
% 
% 
% figure(71); clf;
% hold on; plot(lat, mean(OBS_OLR') - RM*trapz(sind(lat), mean(OBS_OLR'))/2 , 'b-', 'LineWidth', LW );
% hold on; plot(lat, MOD_Z_OLR' - RM*trapz(sind(lat), MOD_Z_OLR)/2, 'b--', 'LineWidth', LW );
% 
% hold on; plot(lat, mean(OBS_OLR'-OBS_LDS'+OBS_LUS') - RM*trapz(sind(lat), mean(OBS_OLR'-OBS_LDS'+OBS_LUS') )/2 , 'r-', 'LineWidth', LW );
% hold on; plot(lat, (MOD_Z_OLR-mean(MOD_LDS') + mean(MOD_LUS')) - RM*trapz(sind(lat), (MOD_Z_OLR-mean(MOD_LDS') + mean(MOD_LUS')))/2 , 'r--', 'LineWidth', LW );
% 
% % and the inmplied biases
% PW=10^(-15);
% % figure(98);clf;
% % hold on; plot(lat, implied_MHT(OBS_Z_SHF',lat), 'k-', 'LineWidth', LW)
% % hold on; plot(lat, implied_MHT(MOD_Z_SHF',lat), 'k--', 'LineWidth', LW)
% % hold on; plot(lat, implied_MHT(-OBS_Z_E',lat), 'g-', 'Color', green, 'LineWidth', LW)
% % hold on; plot(lat, implied_MHT(-MOD_Z_E',lat), 'g--','Color', green, 'LineWidth', LW)


%replace ocean only nan with zeros
OBS_Z_SHF(find(isnan(OBS_Z_SHF) == 1)) = 0;
OBS_Z_E(find(isnan(OBS_Z_E) == 1)) = 0;
OBS_Z_SH(find(isnan(OBS_Z_SH) == 1)) = 0;

MOD_Z_SHF(find(isnan(MOD_Z_SHF) == 1)) = 0;
MOD_Z_E(find(isnan(MOD_Z_E) == 1)) = 0;
MOD_Z_SH(find(isnan(MOD_Z_SH) == 1)) = 0;

% AND IMPLIES oht
MOD_OHT_tot = [0; cumsum(dA.*( (MOD_Z_SHF'-MOD_G_SHF) ))];
MOD_OHT_E = [0; cumsum(dA.*( -(MOD_Z_E'-MOD_G_E) ))];
MOD_AHT_E = [0; cumsum(dA.*( (MOD_Z_E'-MOD_G_E) ))];
MOD_AHT_SH = [0; cumsum(dA.*( (MOD_Z_SH'-MOD_G_SH) ))];
MOD_OHT_SRAD = [0; cumsum(dA.*( (MOD_Z_SRAD'-MOD_G_SRAD) ))];
MOD_OHT_SH = [0; cumsum(dA.*( -(MOD_Z_SH'-MOD_G_SH) ))];
MOD_AHT_ARAD = [0; cumsum(dA_tot.*( (MOD_Z_ARAD'-MOD_G_ARAD) ))];
FL= (MOD_Z_SRAD_L'-MOD_G_SRAD_L);
FL(find(isnan(FL) ==1)) = 0'



MOD_AHT_SRAD = [0; cumsum(dA_L.*(FL  ))];

figure(32);clf; 
plot(mlat, MOD_AHT_SRAD, 'k-', 'LineWidth', LW)





%%% inmplied heat transports

OBS_OHT_tot = [0; cumsum(dA.*( (OBS_Z_SHF'-OBS_G_SHF) ))];
OBS_OHT_E = [0; cumsum(dA.*( -(OBS_Z_E'-OBS_G_E) ))];
OBS_AHT_E = [0; cumsum(dA.*( (OBS_Z_E'-OBS_G_E) ))];
OBS_AHT_SH = [0; cumsum(dA.*( (OBS_Z_SH'-OBS_G_SH) ))];
OBS_OHT_SRAD = [0; cumsum(dA.*( (OBS_Z_SRAD'-OBS_G_SRAD) ))];
OBS_OHT_SH = [0; cumsum(dA.*( -(OBS_Z_SH'-OBS_G_SH) ))];
OBS_AHT_ARAD = [0; cumsum(dA_tot.*( (OBS_Z_ARAD'-OBS_G_ARAD) ))];

FL= (OBS_Z_SRAD_L'-OBS_G_SRAD_L);
FL(find(isnan(FL) ==1)) = 0'

OBS_AHT_SRAD = [0; cumsum(dA_L.*(FL  ))];

figure(32); hold on; 
plot(mlat, OBS_AHT_SRAD, 'k--', 'LineWidth', LW)

%  figure(3); hold on;
%  %hold on; plot(mlat, PW*OBS_OHT, 'b-', 'LineWidth', LW);
% hold on; plot(mlat, PW*MOD_AHT_ARAD, 'r--', 'LineWidth', LW);
% set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', LWA); grid on; box on;
% axis([-90 90 -1 0.5]); 
% title('Implied AHT of radiation', 'FontSize', FSz, 'FontName', FN)

%start with the heat transport biases inferred from (CERES) TOA and
%reanlysis (ERA5) AHT compared to CMIP3-6 ensemble mean
load('C:\Users\Aaron Donohoe\Documents\my_papers\AHT_OHT\figs\different_obs\best_OBS_model.mat', 'mlat', 'OBS_MHT' , 'OBS_AHT', 'OBS_OHT','MOD_MHT' , 'MOD_AHT', 'MOD_OHT');

ymn=-0.92; ymx=0.8;
ytic=-0.8:.2:.8;
%atmospheric biases 
 figure(10); clf;

 %hold on; plot(mlat,PW*(OBS_AHT-MOD_AHT'), 'r-', 'LineWidth', LW); %calculated from residual of 
  hold on; plot(mlat, PW*(OBS_AHT_E-MOD_AHT_E), 'g-','Color', green, 'LineWidth', LW);
hold on; plot(mlat, PW*(OBS_AHT_ARAD-MOD_AHT_ARAD), 'r-', 'Color', orange, 'LineWidth', LW);
  hold on; plot(mlat, PW*(OBS_AHT_SH-MOD_AHT_SH), '-','Color', red, 'LineWidth', LW);
  hold on; plot(mlat, PW*(OBS_AHT_SRAD-MOD_AHT_SRAD), '-','Color', purp, 'LineWidth', LW);
  hold on; plot(mlat, PW*(OBS_AHT_ARAD-MOD_AHT_ARAD) +PW*(OBS_AHT_E-MOD_AHT_E)+PW*(OBS_AHT_SH-MOD_AHT_SH) + PW*(OBS_AHT_SRAD-MOD_AHT_SRAD), 'k-', 'Color', black, 'LineWidth', LW);
  figure(10); hold on; plot(mlat, PW*(OBS_AHT-MOD_AHT'), 'k--', 'LineWidth', LW)


set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', LWA, 'XTick', xtic, 'YTick', ytic); grid on; box on;
axis([-90 90 ymn ymx]); 
%title({'Bias (obs - models) in atmospheric heat transport'; 'and contribution of evaporation bias (dashed)'}, 'FontSize', FSz, 'FontName', FN)
hold on; plot([-90 90], [0 0], 'LineWidth', LWA, 'Color', .3*[1 1 1])
xlabel('Latitude', 'FontSize', FSz, 'FontName', FN)
ylabel('Northward heat transport bias (PW)', 'FontSize', FSz, 'FontName', FN)

 %Ocean biases
figure(11); clf;

hold on; plot(mlat, PW*(OBS_OHT_E-MOD_OHT_E), 'b-', 'Color', green,'LineWidth', LW);
hold on; plot(mlat, PW*(OBS_OHT_SRAD-MOD_OHT_SRAD), '-','Color', orange, 'LineWidth', LW);
hold on; plot(mlat, PW*(OBS_OHT_SH-MOD_OHT_SH), '-','Color', red, 'LineWidth', LW);
% %% the sume of the terms
 hold on; plot(mlat, PW*(OBS_OHT_E-MOD_OHT_E)+PW*(OBS_OHT_SRAD-MOD_OHT_SRAD)+PW*(OBS_OHT_SH-MOD_OHT_SH) , 'k-','LineWidth', LW);
hold on; plot(mlat, PW*(OBS_OHT-MOD_OHT), 'k--', 'LineWidth', LW);

set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', LWA, 'XTick', xtic, 'YTick', ytic); grid on; box on;
axis([-90 90 ymn ymx]); 
%title({'Bias (obs - models) in ocean heat transport'; 'and contribution of evaporation bias (dashed)'}, 'FontSize', FSz, 'FontName', FN)
hold on; plot([-90 90], [0 0], 'LineWidth', LWA, 'Color', .3*[1 1 1])
xlabel('Latitude', 'FontSize', FSz, 'FontName', FN)
ylabel('Northward heat transport bias (PW)', 'FontSize', FSz, 'FontName', FN)

%%% add in the full model biases



return
%%% biases in spatial anomaly of srad and turb

figure(2); clf;
hold on; plot(obs_lat, (OBS_Z_SHF - OBS_G_SHF)- (MOD_Z_SHF - MOD_G_SHF), 'b-', 'LineWidth', LW)
hold on; plot(obs_lat, (OBS_Z_SRAD - OBS_G_SRAD)- (MOD_Z_SRAD - MOD_G_SRAD), 'r-', 'LineWidth', LW)
hold on; plot(obs_lat, (OBS_Z_TURB - OBS_G_TURB)- (MOD_Z_TURB - MOD_G_TURB), 'k-', 'LineWidth', LW)

%% impact on OHT biases
OHT_BIAS = OBS_OHT - MOD_OHT;
RAD_OHT_BIAS = [0; cumsum(dA.*((OBS_Z_SRAD' - OBS_G_SRAD)- (MOD_Z_SRAD' - MOD_G_SRAD)))];
TURB_OHT_BIAS = [0; cumsum(dA.*((OBS_Z_TURB' - OBS_G_TURB)- (MOD_Z_TURB' - MOD_G_TURB)))]; 

figure(4); clf;;
hold on; plot(mlat, PW*OHT_BIAS, 'b-', 'LineWidth', LW);
hold on; plot(mlat, PW*RAD_OHT_BIAS, 'r-', 'LineWidth', LW);
hold on; plot(mlat, PW*TURB_OHT_BIAS, 'k-', 'LineWidth', LW);
