close all
clear all
% temperatur og trykk for en standard atmosf?re:
load atmos.dat
temp = atmos(:,1);
trykk = atmos(:,2);
% blandingsforholdet av gassene i atmosf?ren, brukes kun for SW
load r_gass.dat
r_h2o=r_gass(:,1);
r_co2=r_gass(:,2);
r_o3=r_gass(:,3);
%absorpsjonskonstanter [m^2/kg]:
load abs_kons.dat
k_h2o_1=abs_kons(:,1);
k_h2o_2=abs_kons(:,2);
k_h2o_3=abs_kons(:,3);
k_hartley=abs_kons(:,4);
k_huggins=abs_kons(:,5);
k_chappuis=abs_kons(:,6);

%emissivitetene for hver gass [m^-1]:
load emiss.dat
emiss_h2o = emiss(:,1)*10;
emiss_co2 = emiss(:,2)*2; 
emiss_o3  = emiss(:,3);


R=287.05;                           % gasskonstant t?rr luft (J/(K*kg))
S0 = 1367;                          % solkonstanten, Wm-2
breddegrader = 45;                  % 
deklinasjonsvinkel = 0;             % 0=21 mars
sigma = 5.67e-8;                    % Stefan-Boltzmann konstant [W/(m2*K4)] 
Cp = 1004;                          % Spesifikk varme for t?rr luft ved konstant trykk [J/K*kg]
Cs = 1e+6;                          % Spesifikk varme (Cs=3e+5 tilsvarer 10 cm jord) [J/K*m2]
Alb = 0.3;                          % Bakkealbedo
In=10;                              % langb?lget str?ling ved 51km [W/m2]
T_s(1)=288.2;                       % Bakketemp, Kelvin

t_max=2400;                           % hvor mange timer vi ser p? frem i tid
lag=1:50;                           % lagene i modellen
temp=temp(:,1);                     % Kelvin [K]
trykk=trykk(:,1);                   % [hPa]
%bruker gassloven for ? finne tettheten:
tetthet=trykk(1:50)./(R*temp(1:50))*100; %gange med 100 for ? f? pascal! [kg/m3]
dz=1000;                                 % tykkelsen p? hvert lag, [m]
H_SH=4000;                                 % skalah?yde SH, [m]

%H2O: tilgjengelig energi for 3 b?nd for vanndamp:
delband1 = 0.73;
delband2 = 0.20;
delband3 = 0.07;
del_h2o = 0.346*S0;

%O3: tilgjengelig energi for 3 b?nd for ozon:
del_hartley = 0.002;
del_huggins = 0.01;
del_chappuis =.988;
del_o3=0.654*S0;

%NB! n?r du skal finne optisk tykkelse for hvert lag trenger du ikke en
%forl?kke! Det holder med matrisemultiplisering. Om du liker l?kker, s? er
%det ogs? helt ok ? lage:)
%optisk tykkelse H2O:
tau1 = 0;             %oppgave1
tau2 = 0;             %oppgave1
tau3 = 0;             %oppgave1

% optisk tykkelse O3:
tauH=0;               %oppgave1
tauHU=0;              %oppgave1
tauC=0;               %oppgave1


%deklarasjon av variabler, disse skal IKKE endres!
% dette er kun for raskere kj?ring av programmet
mu=zeros(t_max,1);
S_ned_1=zeros(51,t_max);
S_ned_2=zeros(51,t_max);
S_ned_3=zeros(51,t_max);
DeltaF1=zeros(50,t_max);
DeltaF2=zeros(50,t_max);
DeltaF3=zeros(50,t_max);
Temp1=zeros(50,t_max);
Temp2=zeros(50,t_max);
Temp3=zeros(50,t_max);
Temp_h2o=zeros(50,t_max);
S_ned_H=zeros(51,t_max);
S_ned_HU=zeros(51,t_max);
S_ned_C=zeros(51,t_max);
DeltaF_H=zeros(50,t_max);
DeltaF_HU=zeros(50,t_max); 
DeltaF_C=zeros(50,t_max);
Temp_H=zeros(50,t_max);
Temp_HU=zeros(50,t_max);
Temp_C=zeros(50,t_max);
Temp_o3=zeros(50,t_max);
Temp_kort=zeros(50,t_max);
totTempkort=zeros(50,1);
totTemp=zeros(50,1);
tempO3=zeros(50,1);
tempH2O=zeros(50,1);
T=zeros(50,t_max);
Temp_lang=zeros(50,t_max);
F_oppH2o=zeros(51,t_max);
F_oppCo2=zeros(51,t_max);
F_oppO3=zeros(51,t_max);
DeltaF_oppH2o=zeros(50,t_max);
DeltaF_oppCo2=zeros(50,t_max);
DeltaF_oppO3=zeros(50,t_max);
F_nedH2o=zeros(51,t_max);
F_nedCo2=zeros(51,t_max);
F_nedO3=zeros(51,t_max);
DeltaF_nedH2o=zeros(50,t_max);
DeltaF_nedCo2=zeros(50,t_max);
DeltaF_nedO3=zeros(5,t_max);
DeltaF_H2o=zeros(50,t_max);
DeltaF_Co2=zeros(50,t_max);
DeltaF_O3=zeros(50,t_max);
DeltaF_SH=zeros(50,t_max);
Temp_lang_H2o = zeros(50,t_max);
Temp_lang_Co2 = zeros(50,t_max);
Temp_lang_O3  = zeros(50,t_max);
Temp_SH  = zeros(50,t_max);
totTemp_H2o=zeros(50,1);
totTemp_Co2=zeros(50,1);
totTemp_O3=zeros(50,1);
totTemplang=zeros(50,1);
Total_temp=zeros(50,1);
S_ned_s = zeros(t_max);
F_ned_s = zeros(t_max);

%NB! temperaturen i atmosf?ren kalles T
% Initialiserer med T=250K i alle lag, alternativt initialiserer
% med standard atm.
%T(:,1)=temp(1:50);  
T(:,1)=250;

%her begynner tidsl?kka!
for t=1:t_max;

%KORTB?LGET STR?LING:

%cosinus til sentitvinkelen
mu(t)= sind(breddegrader)*sind(deklinasjonsvinkel)...
    + cosd(breddegrader)*cosd(deklinasjonsvinkel)*cosd((t+18)*15);  

%H2O:
if (mu(t) > 0)                     % Senitvinkelen st?rre enn null = sola er oppe
    S0_1=del_h2o*delband1;
    S0_2=del_h2o*delband2;
    S0_3=del_h2o*delband3;
    S_ned_1(51,t) = S0_1*mu(t);   % Instr?ling avhengig av senitvinkelen
    S_ned_2(51,t) = S0_2*mu(t);
    S_ned_3(51,t) = S0_3*mu(t);
else 
    S_ned_1(51,t) = 0;   % Sola g?tt ned
    S_ned_2(51,t) = 0;
    S_ned_3(51,t) = 0;
end

% Str?ling nedover i hvert lag avhengig av str?ling i lag over og
% eksponensielt av optisk tykkelse og senitvinkelen; ulik veilengde for str?lingen ? passere
for n=1:50
    S_ned_1(51-n,t) =   0; %oppgave 2
    S_ned_2(51-n,t) =   0; %oppgave 2
    S_ned_3(51-n,t) =   0; %oppgave 2
end

% Forskjell: hva kommer inn i laget minus det som g?r ut av laget. Str?ling
% blir absorbert og kan ha en effekt p? oppvarmingsrate
for n=1:50
    DeltaF1(n,t) = S_ned_1(n+1,t)-S_ned_1(n,t); 
    DeltaF2(n,t) = S_ned_2(n+1,t)-S_ned_2(n,t);
    DeltaF3(n,t) = S_ned_3(n+1,t)-S_ned_3(n,t);
end

%------ Oppvarmingsrate ----K/time
for n=1:50
    Temp1(n,t) = 0; %oppgave 3
    Temp2(n,t) = 0; %oppgave 3
    Temp3(n,t) = 0; %oppgave 3
end

%total oppvarmingsrate grunnet vanndamp:
for n=1:50
     Temp_h2o(n,t) = Temp1(n,t) + Temp2(n,t) + Temp3(n,t);
end

%O3:
if (mu(t) > 0)                         % Senitvinkelen st?rre enn null = sola er oppe
    S0_H=del_o3*del_hartley;           % hvor stor del av energien fra sola som er 
    S0_HU=del_o3*del_huggins;          % tilgjengelig for hvert b?nd
    S0_C=del_o3*del_chappuis;
    S_ned_H(51,t) = S0_H*mu(t);         % Innstr?ling avhengig av senitvinkelen
    S_ned_HU(51,t) = S0_HU*mu(t);
    S_ned_C(51,t) = S0_C*mu(t);
else 
    S_ned_H(51,t) = 0;   % Sola g?tt ned
    S_ned_HU(51,t) = 0;
    S_ned_C(51,t) = 0;
                    
end

% Str?ling nedover i hvert lag avhengig av str?ling i lag over og
% eskponensilt av optisk tykkelse og senitvinkelen; ulik veilengde for str?lingen ? passere
for n=1:50
    S_ned_H(51-n,t) = 0;      %oppgave 2
    S_ned_HU(51-n,t) = 0;      %oppgave 2
    S_ned_C(51-n,t) = 0;      %oppgave 2
end

% For bakketemperaturen
    S_ned_s(t) = S_ned_1(1,t) + S_ned_2(1,t) + S_ned_3(1,t) +S_ned_H(1,t) + S_ned_HU(1,t) + S_ned_C(1,t) ; 

% Forskjell: hva kommer inn i laget minus det som g?r ut av laget. Str?ling
% blir absorbert og kan ha en effekt p? oppvarmingsrate
for n=1:50
    DeltaF_H(n,t) = S_ned_H(n+1,t)-S_ned_H(n,t); 
    DeltaF_HU(n,t) = S_ned_HU(n+1,t)-S_ned_HU(n,t); 
    DeltaF_C(n,t) = S_ned_C(n+1,t)-S_ned_C(n,t); 
end

%------ Oppvarmingsrate ----K/time
for n=1:50
    Temp_H(n,t) = 0;   %oppgave 3
    Temp_HU(n,t) = 0;  %oppgave 3
    Temp_C(n,t) = 0;   %oppgave 3
end

%total oppvarmingsrate grunnet ozon:
for n=1:50
     Temp_o3(n,t) = Temp_H(n,t) + Temp_HU(n,t) + Temp_C(n,t);
 end
%total oppvarmingsrate kortb?lget (K/time):
 for n=1:50
     Temp_kort(n,t) =Temp_o3(n,t) + Temp_h2o(n,t);% 
 end

%LANGB?LGET STR?LING:
% %-------- str?ling oppover ---------------------------------------
%bruker Stefan-Boltzmann lov, antar at jorda er et svartlegeme
 F_oppH2o(1,t) = sigma*(T_s(t))^4;  
 F_oppCo2(1,t) = sigma*(T_s(t))^4;
 F_oppO3(1,t)  = sigma*(T_s(t))^4;
 
 %langb?lget str?ling oppover for hvert lag
 for n=1:50
     F_oppH2o(n+1,t) = 0; %oppgave 4
     F_oppCo2(n+1,t) = 0; %oppgave 4
     F_oppO3(n+1,t)  = 0; %oppgave 4
 end

% Forskjell: hva kommer inn i laget minus det som g?r ut av laget. Str?ling
% blir absorbert og kan ha en effekt p? oppvarmingsrate
for n=1:50
% TB2011, feil. Rettet 16.02 2012
    DeltaF_oppH2o(n,t) =  F_oppH2o(n,t)-F_oppH2o(n+1,t);                  
    DeltaF_oppCo2(n,t) =  F_oppCo2(n,t)-F_oppCo2(n+1,t);
    DeltaF_oppO3(n,t)  =   F_oppO3(n,t)-F_oppO3(n+1,t);
end

% %-------- str?ling nedover ---------------------------------------
% langb?lget str?ling ned fra 51 km
F_nedH2o(51,t) = In;
F_nedCo2(51,t) = In;
F_nedO3(51,t)  = In;
%langb?lget str?ling nedover for hvert lag
for n=1:50
    F_nedH2o(51-n,t) = 0; %oppgave 5
    F_nedCo2(51-n,t) = 0; %oppgave 5
    F_nedO3(51-n,t)  = 0; %oppgave 5
end

% For bakketemperaturen
    F_ned_s(t) = F_nedH2o(1,t) + F_nedCo2(1,t) + F_nedO3(1,t) ; 

    % Forskjell: hva kommer inn i laget minus det som g?r ut av laget. Str?ling
% blir absorbert og kan ha en effekt p? oppvarmingsrate
for n=1:50
    DeltaF_nedH2o(n,t) = F_nedH2o(n+1,t) - F_nedH2o(n,t);
    DeltaF_nedCo2(n,t) = F_nedCo2(n+1,t) - F_nedCo2(n,t);
    DeltaF_nedO3(n,t)  = F_nedO3(n+1,t) - F_nedO3(n,t);
end 

%---- Total langb?lget fluks ----------------------------------  
for n=1:50
    DeltaF_H2o(n,t) =  DeltaF_nedH2o(n,t) + DeltaF_oppH2o(n,t) ;
    DeltaF_Co2(n,t) =  DeltaF_nedCo2(n,t) + DeltaF_oppCo2(n,t) ;
    DeltaF_O3(n,t)  =  DeltaF_nedO3(n,t) + DeltaF_oppO3(n,t) ;
end

%----  Fluks (W/m2) av SH fra bakken, ser bort fra LH ----------------------------------  
if (mu(t) > 0)                     % Senitvinkelen st?rre enn null = sola er oppe
  SH0 = Cp*tetthet(1)*1.5e-1*(T_s(t)-T(1,t));   % Hartmann likning 4.26
                                              % (CD=3E-3, U=5m/s)
else
  SH0 = 0;
end
    
%Energifluksen ved bakken gitt ved SH0 fordeles vertikalt
for n=1:50
    DeltaF_SH(n,t) = SH0*(exp(-(n-1)*dz/H_SH)-exp(-n*dz/H_SH));
end


    
%------ Oppvarmingsrate ----K/time
for n=1:50
    Temp_lang_H2o(n,t) = 0; %oppgave 6 
    Temp_lang_Co2(n,t) = 0; %oppgave 6
    Temp_lang_O3(n,t) = 0; %oppgave 6
    Temp_SH(n,t) = 0; %oppgave 6   % Hint: DeltaF_SH og DeltaF_H2o
                      %er begge energiflukser
    % total oppvarmingsrate for langb?lget:
    Temp_lang(n,t)=Temp_lang_H2o(n,t)+Temp_lang_Co2(n,t)+Temp_lang_O3(n,t);
end
%------ Total Oppvarmingsrate, b?de kort og lang ----K/time
for n=1:50
    Total_temp(n,t) = Temp_kort(n,t) + Temp_lang(n,t);
end   
% Oppdaterer bakketemperaturen
    T_s(t+1) = T_s(t)+((1-Alb)*S_ned_s(t)+F_ned_s(t)-sigma*(T_s(t))^4-SH0)*3600/Cs;
    
% Oppdaterer temperaturen i atmosf?ren
for n=1:50
T(n,t+1) = T(n,t)+ Total_temp(n,t)+Temp_SH(n,t);
end

end %her ender tidsl?kka:

%finner ut total oppvarmingsrate for det siste d?gnet vi har sett p?
%brukes i figurene nedenfor
if t_max>=24
    t_min=t_max-23;
else
    t_min=1;
end
for n=1:50
    totTemp_H2o(n)=sum(Temp_lang_H2o(n,t_min:t_max));
    totTemp_Co2(n)=sum(Temp_lang_Co2(n,t_min:t_max));
    totTemp_O3(n)=sum(Temp_lang_O3(n,t_min:t_max));
    totTemplang(n)=sum(Temp_lang(n,t_min:t_max));
    totTempkort(n)=sum(Temp_kort(n,t_min:t_max));
    tempO3(n)=sum(Temp_o3(n,t_min:t_max));
    tempH2O(n)=sum(Temp_h2o(n,t_min:t_max));
    totTemp(n)=sum(Total_temp(n,t_min:t_max));
end

T_s(t_max)
T(:,t_max)

%figur over oppvarmingsraten pga. absorbsjon av kortb?lget str?ling, [Kelvin/d?gn]
figure(1)
plot(totTempkort,lag,'Linewidth',3)
hold on
plot(tempO3,lag,'r')
plot(tempH2O,lag,'g')
legend('total', 'O3', 'H2O','Location','SE')
title('Oppvarming kortb?lget pr d?gn')
ylabel('km')
xlabel('Kelvin pr. d?gn')
ylim([1 50])

%figur over oppvarmingsraten pga. absorbsjon av langb?lget str?ling, [Kelvin/d?gn]
figure(2)
plot(totTemplang,lag,'LineWidth',4)
hold on
plot(totTemp_H2o,lag,'g-','LineWidth',3)
plot(totTemp_Co2,lag,'k-','LineWidth',3)
plot(totTemp_O3,lag,'r-','LineWidth',3)
ylim([1 50])
title('Oppvarming langb?lget pr d?gn')
ylabel('km')
xlabel('Kelvin pr. d?gn')
legend('Total', 'H2O', 'CO2','O3','Location','SW')

%figur over total oppvarmingsrate pga. str?ling, [kelvin/d?gn]
figure(3)
plot(totTemp,lag,'b','LineWidth',3)
hold on
plot(totTemplang,lag,'k','Linewidth',3)
plot(totTempkort,lag,'r','LineWidth',3)
legend('Total', 'Lang', 'Kort','Location','SW')
ylim([1 50])
title('Oppvarming pr d?gn')
ylabel('km')
xlabel('Kelvin pr. d?gn')

%figur over temperaturprofil, [kelvin]
figure(4)
plot(T(:,t_max),lag,'b','LineWidth',3)
hold on
legend('Temp','Location','SW')
ylim([1 50])
title('Temperaturprofil')
ylabel('km')
xlabel('Kelvin')