program RMCM;

{Version 15: October 2012}
{4 Short Wave intervals (channels)}
{Ozone concentration prescribed}
{4 infrared intervals (bands), using information from Kondratyev (1969) and Bliss (1961)}
{Water vapour absorption cross sections are shown in figure 2.51 of the lecture notes}
{CO2 absorption band 13-17 micron using method described in Box 2.6 of lecture notes}
{Atmospheric window band 8-13 micron}

{KMAX-Layer Radiative-Moist-Convective-Model}

{Program written by Aarnout van Delden, Utrecht University, Netherlands}
{a.j.vandelden@uu.nl}

const
{+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++}
{localDir1 = '~/Desktop/WorkFolderRMCM15/Input/';}
localDir1 = '/Users/Shared/WorkFolderRMCM15/Input/';
{localDir1 = '/Users/Shared/WorkFolderRMCM15/Input/';}
{Path to folder where input is located}
{ The following files should be in Input: InputRMCM15 (f5), O3mixingratio_25x31 (f4) and CIRAandSA.txt (optional) }
{localDir2 = '~/Desktop/WorkFolderRMCM15/Results/';}
localDir2 = '/Users/Shared/WorkFolderRMCM15/Results/';
{localDir2 = '/Users/Delden/Desktop/WorkFolderRMCM15/Results/';}
{localDir2 = '/Users/Shared/WorkFolderRMCM15/Results/';}
{Path to folder where output is saved}
KMAX = 25;{number of layers}
dimi=31;{number of latitudes}{needed to read Ozone-file}
binmin=150;{minimum temperature in Planck energy partition}
binmax=350;{maximum temperature in Planck energy partition}

channelmax = 4;{SW}
bandmax=4;{LW}
pi = 3.1415927;
sigma = 0.000000056704;

JulianDay0=1.0;{1 corresponds to 1 January}
dm = 149598000000.0; {averaged distance sun_earth in meters; it is 1 AU}

KH = 400.0;{100.0 if KMAX = 100}{reference turbulent/convective exchange coefficient [Wm^-2K^-1]}
KE = 0.25;

C=1000000.0;{10000000.0}{heat capacity of the ground}

ps = 100000.0;{surface pressure}
cp = 1004.0;{specific heat cap. at const. pressure}
Rd = 287.0;{gas constant for dry air [J/(kg K)]}
RCO2 = 188.9;{gas constant for CO2 [J/(kg K)]}
Rv = 461.5;{gas constant for water vapour [J/(kg K)]}
g = 9.8;{acceleration due to gravity}
kBoltzman=1.38*exp(ln(10.0)*(-23));{Boltzmans constant JK^-1}

RHc0=0.98;{0.999}{0.9}{0.99}{parameter related to condensation scheme (see SPEEDY)}
DRHc=0.10;{0.20}{0.1}{0.65}{parameter related to condensation scheme (see SPEEDY)}
RH0=0.34;{0.3777: De Bilt. Badajoz:    0.34}{parameters for parametrization of cloud cover}
Ccc=1.28;{1.503: De Bilt.  Badajoz: 1.28}{parameters for parametrization of cloud cover}

{Other parameters are specified in the Inputfile: InputRMCM}
{+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++}

{USEFUL UNIX COMMANDS:
 gpc RCM.p
 gpc -O3 -funroll-loops RMCM15.p                  (MAC-G5-processor)
 gpc -O3 -funroll-loops -msse3 RMCM15.p		     (Intel-processor)
 time ./a.out
}
{#L#}{<------long wave radiation budget is computed}
{#S#}{<------Solar radiation budget is computed}


type
		vector1 = array[1..KMAX] of double;
		vector2 = array[0..4] of double;
		vector3 = array[0..channelmax] of double;{for channelmax channels of Solar radiation}
		vector4 = array[1..bandmax] of double;{for bandmax bands of long-wave radiation}
		vector5 = array[1..dimi] of double;
		matrix1 = array[1..KMAX,0..4] of double;
		matrix2 = array[1..KMAX,1..channelmax] of double;{for channelmax channels of Solar radiation}
		matrix3 = array[1..dimi,1..KMAX,1..12] of double;
		matrix4 = array[1..KMAX,1..bandmax,1..2] of double;{for bandmax channels of LW radiation}
		matrix5 = array[binmin..binmax,1..bandmax] of double;{for bandmax bands of long-wave radiation}


var
{Input-Variables}
SUN:integer;{0: then Solar irradiance=constant in time and equal to global average; 1: Solar irradiance follow yearly cycle depending on JulianDay and latitude}
INITIALPROFILE:string[13];
Runnumber,s2:string[2];
tmax,ozone,soilwetness,grey:integer;
S0,PWInit,CCinit,Bowen,Hv0,Lv,pCO2,qstrat,epsg,alfag,alfac,LWP,Tauc:real;
{End Input-Variables}

sigmaO3,sigmaCO2,sigmaL,sigmav,TotalOpticalPath:vector4;
eps,OptPathLW:matrix4;
LWEnergy:matrix5;
bin:integer;
LWbeam,SWbeam:integer;

F,OptPathSW,abscoeff,p,H,IRa,Ua,SR,k,TauOpt,ro,rov,qO3,qCO2,qv,qL,dz,z,NBV,RH,rs,qvs,dqv,RHc,FLatentHeat,Theta,mrO3Mean,opticalpathCO2:vector1;
Tg,PW:vector2;
T:matrix1;
mrO3:matrix3;
sum:vector5;

OpticalPathSW:matrix2;{Optical path=0 if beam is not absorbed}
gamma,SgD,SgU,sigmaSWO3,sigmaSWwv:vector3;

albedoCT,Emission:real;
CT,nCT:integer;
EF:real;
modeEF,modeCC:integer;

SumGamma,OpticalPathSW0,Zenit,Fg,FPW,albedo,ER,e,Ug,Sgtot,Sabsorbed,SRD,SRU,SRC,B,CounterRadiation,Hs,Ts,dt,time,OLRTOA,kappa,omega,QH,QE,dp,dummy,dummy1,dummy2,mu,TauS,Hopt,UgOut,IRnet,eg,esg,es,Evaporationenergyflux,rg,RHg,Hv,LR,DH:double;{Hv is the water vapour scale height}
thetaupper,thetalower,dEXNER,phi,daymonth,Absphi:real;

q0:real;{q0 is the basic state minimum specific humidity. Important to have some water vapour in the stratosphere}
PW0,PWprime:real;{PW=PW0+PWprime; PW0 is the precipitable water if the specific humidity over the full depth of the atmosphere were the same as the specific humidity in the stratosphere}
Rovg,Rovg0,Rovgprime,Rov0,Rovprime:real;

CC:real;{cloud cover, defined using the max RH}
i,j,run,channel,band:integer;{i is used to denote layer-number; j is used to denote substep in Runge-Kutta scheme}
n,m,mm,mon,level,CloudTopLevel:integer;
nt,dimt,tf:longint;

delta:real;{solar declination angle}
Q:real;{Solar irradiance: function of time and latitude}
JulianDay:double;

nrun,runmin,runmax:integer;

sumSR,SumLatentHeat:real;{#####}
f1,f2,f3,f4,f5,f6:text;
label
1,2;

{#############################################}

{BEGIN: SOLAR INPUT AS A FUNCTION OF LATITUDE AND JULIAN DAY}
procedure SOLARINPUT(phi,Day:real);
{This procedure has been written by Niels Zweers (Feb. 2007)}

{phi is latitude in degrees}
{JulianDay is Day of the year JulianDay=1 on January 1; JulianDay=80 on March 21}
var
time2, d, Hday:real;

begin
time2:=trunc(Day);{no daily cycle!!!}

d:=dm*(1+0.0178*sin((2*pi/365.0)*(285.0+time2))); {distance between the earth and the sun parameterized as a function of time}
{d[t]:=dm;} {distance between the earth and the sun if orbit is circular, with radius 1 AU}

delta:=(23.45*pi/180)*sin((2*pi/365.0)*(285.0+time2)); {declination angle of the sun as function of time}

	if phi=-90 then begin
	if (delta>-0.00001) and (delta<0.00001) then Hday:=pi/2;
	if delta<-0.00001 then Hday:=pi;
	if delta>0.00001 then Hday:=0;
	end;{if}
	
	if phi=90 then begin
	if (delta>-0.00001) and (delta<0.00001) then Hday:=pi/2;
	if delta<-0.00001 then Hday:=0;
	if delta>0.00001 then Hday:=pi;
	end;{if}
	
	if (phi<>-90) and (phi<>90) then begin
	if (sin((pi/180)*phi)/cos((pi/180)*phi))*(sin(delta)/cos(delta)) > 1 then Hday:=pi; {24 hours of sunlight on a day}
	if (sin((pi/180)*phi)/cos((pi/180)*phi))*(sin(delta)/cos(delta)) < -1 then Hday:=0; {24 hours of darkness on a day}
	if ((sin((pi/180)*phi)/cos((pi/180)*phi))*(sin(delta)/cos(delta))<=1) and ((sin((pi/180)*phi)/cos((pi/180)*phi))*(sin(delta)/cos(delta))>=-1) then Hday:=arccos(-1*(sin((pi/180)*phi)/cos((pi/180)*phi))*(sin(delta)/cos(delta)));
	end;{if}

	Q:=(S0/pi)*sqr(dm/d)*((Hday*sin((pi/180)*phi)*sin(delta))+(cos((pi/180)*phi)*cos(delta)*sin(Hday))); {insolation Q at the TOA as a function of time and phi}
{Q=average insolation during one day as a function of Julian day and latitude}

Hday:=24*Hday/pi;{length of the day in hours}
if Hday>0 then Zenit:=arccos((Q*24/(S0*Hday))*sqr(d/dm));{in radians}
if Hday=0 then Zenit:=pi/2+0.0000000000000000000000000001;{in radians}
{Zenit=daily average Zenit angle: needed to calculate the aborption of Solar radiation by ozone}

if SUN=0 then Q:=S0/4;{global average constant insolation}
if SUN=0 then Zenit:=pi/3;{pi/3}
{See Ramanathan&Coakley (1978), p. 474: cosine of the solar zeith angle is 1/2, i.e. Zenit angle is 60� or pi/3}

{global average and time average cos(Zenit)=0.387??? }

end;
{END: SOLAR INPUT AS A FUNCTION OF LATITUDE AND JULIAN DAY}
{#############################################}
procedure SATVAPOURPRESSURE(Temp:real);
var
aa,bb:real;
begin
aa:=17.269388;
bb:=35.86;
es:=610.78*exp((aa*(Temp-273.16))/(Temp-bb));
end;
{end: procedure SATVAPOURPRESSURE}
{#############################################}
procedure PLANCK;

const
maxintervals=22;
h=6.63*exp(ln(10.0)*(-34));
c=299792458.0;
pi=3.14159;

type
vector9 = array[0..maxintervals] of double;

var
Temp:real;
B:vector9;
Bintegrated:vector4;
labda,dlabda,Btot,Tblackbody,sumBintegrated:real;
interval:integer;
band:integer;

begin

for bin:=binmin to binmax do begin

dlabda:=0.000001;
labda:=0.0;
Temp:=bin+0.5;
Btot:=sigma*sqr(Temp)*sqr(Temp);{Stefan-Boltzman}

for band:=1 to bandmax do begin
Bintegrated[band]:=0.0;
end;

B[0]:=0;
labda:=0.0000065;
for interval:= 7 to 22 do begin
labda:=labda+dlabda;
B[interval]:=(2*pi*h*sqr(c)/(labda*labda*labda*labda*labda))/(exp(h*c/(labda*kBoltzman*Temp))-1);
end;

labda:=0.0000065;
for interval:= 7 to 20 do begin
labda:=labda+dlabda;

{LW}
if (labda>=0.000008) and (labda<0.000013) then Bintegrated[1]:=Bintegrated[1]+(0.5*(B[interval]+B[interval+1])*dlabda);{window}
if (labda>=0.000013) and (labda<0.000017) then Bintegrated[2]:=Bintegrated[2]+(0.5*(B[interval]+B[interval+1])*dlabda);{CO2-absorbption band}
if (labda>=0.000017) and (labda<0.000021) then Bintegrated[3]:=Bintegrated[3]+(0.5*(B[interval]+B[interval+1])*dlabda);{CO2-absorbption band}
end;{interval}

sumBintegrated:=0;
for band:=1 to bandmax-1 do begin
sumBintegrated:=sumBintegrated+Bintegrated[band];
end;
Bintegrated[bandmax]:=Btot-sumBintegrated;

for band:=1 to bandmax do begin
LWEnergy[bin,band]:=Bintegrated[band]/Btot;{fraction of total energy that is emitted in a certain long wave band}
end;

end;{bin}

end;
{end: procedure PLANCK}

{#############################################}

{#############################################}
procedure Stringof1(NUMBER:integer);
var
			f6: text;
begin
rewrite(f6,'HELPFILE');
if NUMBER<100 then begin write(f6,NUMBER:2);close(f6);end;
reset(f6,'HELPFILE');
if NUMBER<100 then read(f6, Runnumber);
close(f6);
end;
{#############################################}
{%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%}
Function TransmissivityCO2(W:double):double;
{This procedure calculates the optical path of in the CO2 absorption band due to CO2 as a function the column abundance, W in kg per sq. m}
{First it calculates the transmissivity  using eq. 4 in box 2.5}

const
nu0=66750.0;{m^^-1}
Avogadro=6.02214129*exp(ln(10.0)*(23));
max=10000;
{sigmazeroCO2=3.71*exp(ln(10.0)*(-23));}
sigmazeroCO2=0.025*exp(ln(10.0)*(-23));{0.02}
rplus=0.00038;
rmin=0.00044;
Dlabda=0.001*0.000001;

type
		vector = array[0..max] of double;

var
ACS:vector;
labda,labda1,labda2,Tau:double;
n,labdaintervals:integer;

begin
labda1:=13.0*0.000001;labda2:=17.0*0.000001;Tau:=0.0;

labdaintervals:=trunc((labda2-labda1)/Dlabda);
labda:=labda1-(0.5*Dlabda);
{labdaintervals intervals in wavelength}
for n:=1 to labdaintervals do begin
labda:=labda+Dlabda;
if labda>=1/nu0 then ACS[n]:=(Avogadro/0.044)*sigmazeroCO2*exp(-rmin*abs((1/labda)-nu0));
if labda<1/nu0 then ACS[n]:=(Avogadro/0.044)*sigmazeroCO2*exp(-rplus*abs((1/labda)-nu0));
end;{n}

{compute transmission coefficient with equation 4 of box 2.5}
Tau:=0.0;
for n:=1 to labdaintervals do begin
Tau:=Tau+(Dlabda*(exp(-ACS[n]*W)));
end;
TransmissivityCO2:=Tau/(labda2-labda1);
end;
{%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%}

begin
rewrite(f6,'RMCM_INIT');
close(f6);

writeln('');
writeln('=======================================');
writeln('RADIATIVE-MOIST-CONVECTIVE MODEL (RMCM)');
writeln;
writeln('Author: Aarnout van Delden (a.j.vandelden@uu.nl)');
writeln('Version 15: October 2012');
writeln('Path of workfolder, which contains input and output:');
writeln('/Users/Shared/WorkFolderRMCM15/');
writeln('');
writeln('Number of layers is',KMAX:4);
writeln('Global and yearly average Insolation is ',SUN:2,'(1: constant; 0: seasonal cycle)');
writeln('Greenhouse "gases" are water vapour, CO2, and clouds ');
writeln('CO2- and Ozone-distribution are prescribed');
writeln;
writeln('Solar radiation: 4 channels');
writeln('channel 1 is Hartley and Huggins band (<0.34 micrometre)');
writeln('channel 2 radiation in the wavelength range: 0.34-0.5 micrometre');
writeln('channel 3 is Chappuis band (0.5-0.7 micrometre)');
writeln('channel 4 is Near infrared (>0.7 micrometre)');
writeln;
writeln('Long wave (terrestrial) radiation: 4 bands');
writeln('Band 1 is window band: 8-13 micron');
writeln('Band 2 is "CO2" absorption band: 13-17 micron');
writeln('Band 3 is "H2O" absorption band: 17-21 micron');
writeln('Band 4 is aggregate of remaining radiation absorbed by water-vapour & clouds');
writeln;
writeln('Initial temperature: Isothermal');
writeln('');
writeln('Output can be found in the following three files:');
writeln(concat(localDir1,'EndstateRun..'));
writeln(concat(localDir2,'TimeevolutionRun..'));
writeln(concat(localDir2,'EnergyfluxesRun..'));
writeln('');
writeln('Run:',run:5);
writeln('');
writeln('press RETURN to continue');
readln;

PLANCK;{calculate energy partition of long-wave spectrum in dependence of temperature}
{begin read input}

reset(f5,concat(localDir1,'InputRMCM15'));

readln(f5);
readln(f5,tmax,SUN,phi,run);
readln(f5);
readln(f5,S0,PWInit,CCinit,modeCC,EF,modeEF);
readln(f5);
for channel:=1 to channelmax do begin
read(f5,gamma[channel]);
end;
readln(f5);readln(f5);
readln(f5,Ozone,Soilwetness,grey);
readln(f5);
readln(f5,Bowen,Hv0,Lv);
readln(f5);
readln(f5,pCO2,qstrat,epsg,Tauc);
Tauc:=Tauc*3600.0;{hrs to seconds}
readln(f5);
for band:=1 to bandmax do begin
readln(f5,sigmaO3[band],sigmaCO2[band],sigmaL[band],sigmav[band]);
end;
readln(f5);
for channel:=1 to channelmax do begin
read(f5,sigmaSWwv[channel]);{m^2 kg^-1}{absorption crosssection water vapour}
end;
readln(f5);readln(f5);
readln(f5,alfag,alfac);
close(f5);
{end read input}

ER:=1/(Bowen+1);{Evaporative ratio}
dp:=ps/KMAX;
Bowen:=(1/ER)-1;

SumGamma:=0;
for channel:=1 to channelmax do begin
SumGamma:=SumGamma+gamma[channel];
end;

{****************************************************************************************}
{Specify absorption cross-section O3 Solar radiation in each channel}
sigmaSWO3[1] := 200.0*exp(ln(10.0)*(-24));{m^2}{weighted absorption crosssection ozone in the Hartley band 200-300 nm and Huggins band 300-340 nm }
sigmaSWO3[2] := 0.0;{no absorption due to ozone in channel 2}
sigmaSWO3[3] := 0.285*exp(ln(10.0)*(-24));{m^2}{weighted absorption crosssection ozone 500-700 nm (Chappius band)}{Lindzen&Will, JAS, 1973; D.F.Strobel, JGR, 1978}
sigmaSWO3[4] := 0;{no absorption due to ozone in channel 4}
{****************************************************************************************}

if (KMAX<5) then begin
writeln('Too few layers to capture inhomogeneous atmosphere');
writeln('');
goto 2;
end;{if}





{Specify tracer distribution as a molecule mixing ratio by volume}
if KMAX=400 then reset(f4,concat(localDir1,'O3mixingratio_400x31'));
if KMAX=100 then reset(f4,concat(localDir1,'O3mixingratio_100x31'));
if KMAX=50 then reset(f4,concat(localDir1,'O3mixingratio_50x31'));
if KMAX=25 then reset(f4,concat(localDir1,'O3mixingratio_25x31'));

for mon:=1 to 12 do begin
readln(f4);
for n:=1 to KMAX do begin
for i:=1 to dimi do begin
read(f4,mrO3[i,n,mon]);{i is index of latitude (runs from north to south); n is index of level}
{i=1 is North Pole}
mrO3[i,n,mon]:=Ozone*mrO3[i,n,mon]/1000000.0;
if ozone=0 then mrO3[i,n,mon]:=0.0;
end;end;end;
if ozone=0 then gamma[1]:=1-SumGamma;{if ozone=0 then radiation in channel 1 is not absorbed above model-top}
close(f4);

{Tg, Ta, Te from analytical expressions for one layer model}

JulianDay:=JulianDay0;
if phi<900 then SOLARINPUT(phi,JulianDay);
if phi>900 then Q:=S0/4;


kappa:=Rd/cp;

dt:=300;{timestep in seconds}
{dt:=0.0001;}{&}
dimt:=trunc(tmax*24*3600/dt);{number of  time steps}
tf:=288;{192}{8}{96}{24}{frequency of output in units of dt}
{ settings}{********************************************************}

PW[0]:=PWInit;{initialize precipitable water}
{pCO2:=pCO2/2;}
{Bowen:=0.0;}CCinit:=CCinit-0.1;
runmin:=run;
runmax:=run;
for nrun:=runmin to runmax do begin{<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<nrun}
{pCO2:=pCO2*2;}
{Bowen:=Bowen+0.05;}CCinit:=CCinit+0.1;CC:=CCinit;{initialize cloud cover}

stringof1(nrun);
rewrite(f1,concat(localDir1,'EndstateRun',Runnumber));
rewrite(f2,concat(localDir2,'TimeevolutionRun',Runnumber));
rewrite(f3,concat(localDir2,'EnergyfluxesRun',Runnumber));

{initial condition temperature, pressure and threshold mixing ratio}
p[1]:=0.5*dp;
T[1,0]:=255.0;
RHc[1]:=RHc0+DRHc*(sqr(p[1]/ps)-1);{predefined maximum mixing ratio}
for n:=2 to KMAX do begin
T[n,0]:=255.0;
p[n]:=p[n-1]+(ps/KMAX);
RHc[n]:=RHc0+DRHc*(sqr(p[n]/ps)-1);{predefined maximum mixing ratio}
RH[n]:=0.0;{initialize relative humidity}
dqv[n]:=0.0;
Theta[n]:=(T[n,0]*exp(ln(ps/p[n])*kappa));
end;
Tg[0]:=255.0;

{Brunt Vaisala frequency}
{writeln('    p[hPa]     Theta[K]     NBV[s^-1] ');}
for n:=3 to KMAX-1 do begin
dz[n]:=(Rd*T[n,0]/g)*(ln(p[n]+0.5*dp)-ln(p[n]-0.5*dp));
{NBV[n]:=sqrt((g/Theta[n])*(Theta[n-1]-Theta[n+1])/(2*dz[n]));}
{writeln(p[n]/100:10:1,Theta[n]:10:2,NBV[n]:10:5);}
end;

CloudTopLevel:=KMAX;

{Heat transfer coefficient}
for n:=KMAX downto 1 do begin
k[n]:=KH*KMAX/100;{The value of KH is the value of the transfer coefficient if KMAX=100. Transfer coeffiecient is proportional to KMAX}
end;

{thickness in metres of each layer}
for n:=2 to KMAX do begin
dz[n]:=(Rd*T[n,0]/g)*(ln(p[n]+0.5*dp)-ln(p[n]-0.5*dp));
end;

{height each layer}
z[KMAX]:=0.5*dz[KMAX];{z[KMAX]:=dz[KMAX];CORRECTION:08122010}
for n:=KMAX-1 downto 1 do begin
z[n]:=z[n+1]+dz[n];
end;

{set specific concentrations of each radiatively active constituent}
Hv:=Hv0;
Rovg:=PW[0]/Hv;{Density of water vapour at the ground}
e:=Rovg*Rv*T[KMAX,0];{vapour pressure at level n}
SATVAPOURPRESSURE(T[KMAX,0]);
RHg:=e/es;{relative humidity at level n}

for n:=1 to KMAX do begin
Rov[n]:=Rovg*exp(-z[n]/Hv);
Ro[n]:=p[n]/(Rd*T[n,0]);
qv[n]:=Rov[n]/Ro[n];
e:=Rov[n]*Rv*T[n,0];{vapour pressure at level n}
SATVAPOURPRESSURE(T[n,0]);
RH[n]:=e/es;{relative humidity at level n}
qCO2[n]:=(Rd/RCO2)*pCO2;
qL[n]:=0.0;
{ if (p[n]<85000) and (p[n]>20000) then qL[n]:=0.00001
 else qL[n]:=0.0;}{###}
end;

for n:=1 to KMAX do begin
opticalpathCO2[n]:=-ln(TransmissivityCO2(qCO2[n]*dp/g));
end;

{writeln('TransmissivityCO2=',TransmissivityCO2(qCO2[KMAX]*dp/g):10:5);
writeln('optical path=',-ln(TransmissivityCO2(qCO2[KMAX]*dp/g)):10:5);}


{set LW-optical path and LW-emission coefficient each layer}{**%**}
for band:=1 to bandmax do begin
for n:=1 to KMAX do begin

if grey=1 then begin
OptPathLW[n,band,1]:=((1-CC)*((qO3[n]*sigmaO3[band])+(qCO2[n]*sigmaCO2[band])+(qv[n]*sigmav[band]))*dp/g);
OptPathLW[n,band,2]:=(CC*((qO3[n]*sigmaO3[band])+(qCO2[n]*sigmaCO2[band])+(qv[n]*sigmav[band])+((qL[n]/CC)*sigmaL[band]))*dp/g);
end;

if grey=0 then begin
OptPathLW[n,band,1]:=(1-CC)*((((qO3[n]*sigmaO3[band])+(qv[n]*sigmav[band]))*dp/g)+(sigmaCO2[band]*OpticalPathCO2[n]));{sigmaCO2 is either 1 or 0}
OptPathLW[n,band,2]:=(CC)*((((qO3[n]*sigmaO3[band])+(qv[n]*sigmav[band])+((qL[n]/CC)*sigmaL[band]))*dp/g)+(sigmaCO2[band]*OpticalPathCO2[n]));{sigmaCO2 is either 1 or 0}
end;
{qL is weighted with cloud cover!! because qL is average over cloudy and clear area}

eps[n,band,1]:=1-exp(-OptPathLW[n,band,1]);
eps[n,band,2]:=1-exp(-OptPathLW[n,band,2]);
end;
end;{band}


if SUN=0 then begin
{average area weighted ozone mixing ratio}

{average cosine of latitude}
sum[i]:=0;
for i:=1 to dimi do begin
phi:=90.0-(i-1)*180/(dimi-1);
sum[i]:=sum[i]+cos(phi*pi/180);
end;

for n:=1 to KMAX do begin
mrO3Mean[n]:=0.0;
for i:=1 to dimi do begin
phi:=90.0-(i-1)*180/(dimi-1);
for mon:=1 to 12 do begin
mrO3Mean[n]:=mrO3Mean[n]+(cos(phi*pi/180)*mrO3[i,n,mon]/(12*dimi*sum[i]));
end;{mon}
end;{i}
end;{n}

{specific concentration O3 (qO3) is approximately the same as the mixing ratio}
for n:=1 to KMAX do begin
qO3[n]:=mrO3Mean[n];
end;

end;{SUN=0}
{writeln(' p[hPa]  O3_mixingratio  OptPath[ch1]  OptPath[ch2]  OptPath[ch3]  es[hPa]  RelHum[%] ');}

{absorption depends on optical path which depends on absorber mixing ratio}
for n:=1 to KMAX do begin
if Zenit<pi/2 then begin

for channel:=1 to channelmax do begin
OpticalPathSW[n,channel]:=(sigmaSWO3[channel]*mrO3[i,n,mon]*(Rd/kBoltzman)+(qv[n]*sigmaSWwv[channel]))*(1/cos(Zenit))*dp/g;
if SUN=0 then OpticalPathSW[n,channel]:=(sigmaSWO3[channel]*mrO3Mean[n]*(Rd/kBoltzman)+(qv[n]*sigmaSWwv[channel]))*(1/cos(Zenit))*dp/g;
if Q=0 then OpticalPathSW[n,channel]:=0.0;
end;

end;

end;
{##}

{^settings}{********************************************************}
write(f2,' time[days] ');
for n:=1 to KMAX do begin
write(f2,' ',n:5,' ');
end;writeln(f2,'     day  Tg[K]  albedo OLRTOA[Wm^-2] ASR[Wm^-2] HeatFlux_g[Wm^-2]  PW[kgm^-2]  RHg   Ac   Prec[kgm^-2d^-1]   E[kgm^-2d^-1]');
writeln('     day   Tg[K]  albedo OLRTOA[Wm^-2] ASR[Wm^-2] HeatFlux_g[Wm^-2]  PW[kgm^-2]  RHg   Ac   Prec[kgm^-2d^-1]   E[kgm^-2d^-1]');

albedo:=(alfag*(1-CC))+(alfac*CC);{this is the weighted albedo of the Earth's surface}
albedoCT:=(alfac*CC);{albedo at cloud top, radiation in clear atmosphere is not reflected}

{*********************************************************************}



{#L#}
{energy balance at the earth's surface}
Ug:=0;

{Long wave radiation (emitted by the atmosphere) received at Earth's surface}
for LWbeam:=1 to 2 do begin
for n:=1 to KMAX do begin

{PLANCK1(T[n,0]);}
bin:=trunc(T[n,0]);
for band:=1 to bandmax do begin
if LWbeam=1 then Emission:=(1-CC)*(LWEnergy[bin,band]*eps[n,band,LWbeam]*sigma*sqr(T[n,0])*sqr(T[n,0]));{Energy emitted by layer n in particular band}
if LWbeam=2 then Emission:=CC*(LWEnergy[bin,band]*eps[n,band,LWbeam]*sigma*sqr(T[n,0])*sqr(T[n,0]));{Energy emitted by layer n in particular band}


for m:=n+1 to KMAX do begin
Emission:=Emission*(1-eps[m,band,LWbeam]);
end;{m}

Ug:=Ug+(epsg*Emission);{absorption coefficient of Earth's surface=epsg}
end;{band}end;{n}end;{LWbeam}
{#L#}


{Solar radiation received at Earth's surface} {#S#}
for channel:=1 to channelmax do begin

SgD[channel]:=gamma[channel]*Q;

{SgD: downwelling radiation; SgU: upwelling radiation}

CT:=0;
for n:=1 to KMAX do begin
if (qL[n]>0.0) and (CT=0) then begin
{then cloud top}
CT:=1;
SgD[channel]:=(1-albedoCT)*SgD[channel];
end;{if}
SgD[channel]:=SgD[channel]*exp(-OpticalPathSW[n,channel]);
end;{n}
end;{channel}

{absorption of Solar radiation by the earth's surface}
for channel:=1 to channelmax do begin
Ug:=Ug+((1-alfag)*SgD[channel]);
end;


{Tg[0]:=sqrt(sqrt(Ug/sigma));}{&}
{*********************************************************************}


{Convective heat flux}
for n:=1 to KMAX do begin
H[n]:=0;
end;

time:=0;

write(time/(24*3600):8:1);
write(f2,time/(24*3600):11:3);
for n:=1 to KMAX do begin
{write(T[n,0]:7:1);}
write(f2,T[n,0]:7:1);
end;write(Tg[0]:7:1);write(albedo:9:2);write(OLRTOA:10:1);write((1-albedo)*Q:12:1,'    0');
write(f2,Tg[0]:7:2);write(f2,albedo:7:2);write(f2,OLRTOA:9:2);write(f2,(1-albedo)*Q:9:2,'    0.0');
writeln;writeln(f2);

Hv:=Hv0;FPW:=0.0;

{dimt:=2;}{&}
{>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>TIME LOOP:}
for nt:=1 to dimt do begin
time:=time+dt;

Hv:=Hv0-(3456000*4*FPW);Hv:=Hv0;
if Hv>3000 then Hv:=3000;
if Hv<1700 then Hv:=1700;

JulianDay:=JulianDay+dt/(3600*24);
SOLARINPUT(phi,JulianDay);

{absorption depends on optical path which depends on absorber mixing ratio}
for n:=1 to KMAX do begin
if Zenit<pi/2 then begin
for channel:=1 to channelmax do begin
OpticalPathSW[n,channel]:=(sigmaSWO3[channel]*mrO3[i,n,mon]*(Rd/kBoltzman)+(qv[n]*sigmaSWwv[channel]))*(1/cos(Zenit))*dp/g;
if SUN=0 then OpticalPathSW[n,channel]:=(sigmaSWO3[channel]*mrO3Mean[n]*(Rd/kBoltzman)+(qv[n]*sigmaSWwv[channel]))*(1/cos(Zenit))*dp/g;
end;
end;{if}

for channel:=1 to channelmax do begin
if Q=0 then OpticalPathSW[n,channel]:=0.0;
end;

end;{n}
{##}

{prescribe the minimum specific humidity (important for the stratosphere)}
q0:=qStrat;
PW0:=q0*ps/g;
if PW0>PW[0] then begin
PW0:=0.0;
q0:=0.0;
end;

{Determine cloud cover; depends on RH at the ground }
Rovg0:=ps*q0/(Rd*Tg[0]);
PWprime:=PW[0]-PW0;
Rovgprime:=PWprime/Hv;
Rovg:=Rovg0+Rovgprime;{Density of water vapour at the ground}
e:=Rovg*Rv*Tg[0];{vapour pressure at the ground}
SATVAPOURPRESSURE(Tg[0]);
RHg:=e/es;{relative humidity at level KMAX}
{if modeCC=0 then CC:=CCinit;}{&}
if modeCC=1 then CC:=Ccc*(RHg-RH0);{&}
If CC>=1 then CC:=0.99;
if CC<0 then CC:=0.0;

{Liquid water content}
for n:=1 to KMAX do begin
Rov[n]:=(p[n]*q0/(Rd*T[n,0]))+(Rovgprime*exp(-z[n]/Hv));
Ro[n]:=p[n]/(Rd*T[n,0]);
qv[n]:=Rov[n]/Ro[n];
if (dqv[n]>0) and (CC>0) then qL[n]:=dqv[n]*Tauc/dt;
if dqv[n]<=0 then qL[n]:=0.0;{###}
if CC=0 then qL[n]:=0.0;
SATVAPOURPRESSURE(T[n,0]);
e:=Ro[n]*Rv*T[n,0];{vapour pressure}
RH[n]:=e/es;
end;


{albedo depends on cloud cover:}
albedo:=(alfag*(1-CC))+(alfac*CC);

{set LW-optical path and LW-emission coefficient each layer}{**%**}
for band:=1 to bandmax do begin
for n:=1 to KMAX do begin

if grey=1 then begin
OptPathLW[n,band,1]:=((1-CC)*((qO3[n]*sigmaO3[band])+(qCO2[n]*sigmaCO2[band])+(qv[n]*sigmav[band]))*dp/g);
if (dqv[n]>0) then OptPathLW[n,band,2]:=(CC*((qO3[n]*sigmaO3[band])+(qCO2[n]*sigmaCO2[band])+((RHc[n]*qvs[n])*sigmav[band])+((qL[n]/CC)*sigmaL[band]))*dp/g);
{if condensation then part of PW as cloud(liquid) water}
{qL is weighted with cloud cover!! because qL itself is average over cloudy and clear area}
if (dqv[n]<=0) then OptPathLW[n,band,2]:=(CC*((qO3[n]*sigmaO3[band])+(qCO2[n]*sigmaCO2[band])+(qv[n]*sigmav[band]))*dp/g);
{if no condensation (above & below cloud base: then no cloud water}
end;

if grey=0 then begin
OptPathLW[n,band,1]:=(1-CC)*((((qO3[n]*sigmaO3[band])+(qv[n]*sigmav[band]))*dp/g)+(sigmaCO2[band]*OpticalPathCO2[n]));{sigmaCO2 is either 1 or 0}
if (dqv[n]>0) then OptPathLW[n,band,2]:=(CC*((qO3[n]*sigmaO3[band])+((RHc[n]*qvs[n])*sigmav[band])+((qL[n]/CC)*sigmaL[band]))*dp/g)+(CC*sigmaCO2[band]*OpticalPathCO2[n]);
{if condensation then part of PW as cloud(liquid) water}
{qL is weighted with cloud cover!! because qL itself is average over cloudy and clear area}
if (dqv[n]<=0) then OptPathLW[n,band,2]:=(CC*((qO3[n]*sigmaO3[band])+(qv[n]*sigmav[band]))*dp/g)+(CC*sigmaCO2[band]*OpticalPathCO2[n]);
{if no condensation (above & below cloud base: then no cloud water}
end;{if grey}


eps[n,band,1]:=1-exp(-OptPathLW[n,band,1]);
eps[n,band,2]:=1-exp(-OptPathLW[n,band,2]);
end;
end;{band}


for j:=1 to 4 do begin {substeps in R-K scheme}

albedo:=(alfag*(1-CC))+(alfac*CC);{this is the weighted albedo}
albedoCT:=(alfac*CC);{albedo at cloud top, radiation in clear atmosphere is not reflected}

{*********************************************************************}
for n:=1 to KMAX do begin

end;

{#L#}
{tendency at the earth's surface}
CounterRadiation:=0.0;

{Long wave radiation (emitted by the atmosphere) received at Earth's surface}
for LWbeam:=1 to 2 do begin
for n:=1 to KMAX do begin

{PLANCK1(T[n,j-1]);}
bin:=trunc(T[n,j-1]);
for band:=1 to bandmax do begin
if LWbeam=1 then Emission:=(1-CC)*LWEnergy[bin,band]*eps[n,band,LWbeam]*sigma*sqr(sqr(T[n,j-1]));{Energy emitted by layer n in particular band}
if LWbeam=2 then Emission:=(CC)*LWEnergy[bin,band]*eps[n,band,LWbeam]*sigma*sqr(sqr(T[n,j-1]));{Energy emitted by layer n in particular band}

for m:=n+1 to KMAX do begin
Emission:=Emission*(1-eps[m,band,LWbeam]);{transmission through each layer to Earth's surface}
end;
CounterRadiation:=CounterRadiation+Emission;{radiation emitted by atmosphere reaching Earth}

end;{band}end;{n}end;{LWbeam}

B:=(epsg*CounterRadiation);{B=Counter radiation that is absorbed by Earth}
{a fraction (1-epsg) of B is reflected and must later be added to flux coming from the Earth's surface}{^^}
{absorption coefficient Earth' surface=epsg}

{#L#}

{solar radiation received at Earth's surface}{#S#}
for channel:=1 to channelmax do begin

SgD[channel]:=gamma[channel]*Q;

{SgD: downwelling radiation; SgU: upwelling radiation}

CT:=0;
for n:=1 to KMAX do begin
if (qL[n]>0.0) and (CT=0) then begin
{then cloud top}
CT:=1;
SgD[channel]:=(1-albedoCT)*SgD[channel];
end;{if}
SgD[channel]:=SgD[channel]*exp(-OpticalPathSW[n,channel]);
end;{n}
end;{channel}

{add contribution of each channel}
Sgtot:=0;
for channel:=1 to channelmax do begin
Sgtot:=Sgtot+((1-alfag)*SgD[channel]);
end;

n:=KMAX;
H[KMAX]:=k[n]*(((Tg[j-1]*exp(ln(ps/ps)*kappa))-(T[n,0]*exp(ln(ps/p[KMAX])*kappa))));
if H[KMAX]<0 then H[KMAX]:=0;

for n:=1 to KMAX-1 do begin
H[n]:=k[n]*((T[n+1,j-1]*exp(ln(ps/p[n+1])*kappa))-(T[n,j-1]*exp(ln(ps/p[n])*kappa)));
if H[n]<0 then H[n]:=0;
end;

if modeEF=2 then begin
Rovg:=PW[j-1]/Hv;
e:=Rovg*Rv*Tg[j-1];
SATVAPOURPRESSURE(Tg[j-1]);
RHg:=e/es;
Evaporationenergyflux:=Soilwetness*KE*es*(1-RHg);
end;

If modeEF=1 then Evaporationenergyflux:=Soilwetness*H[KMAX]/Bowen;
If modeEF=0 then Evaporationenergyflux:=EF;                    
if Lv=0 then Evaporationenergyflux:=0;

Fg:=(1/C)*(B+Sgtot-H[KMAX]-Evaporationenergyflux-(epsg*sigma*sqr(Tg[j-1])*sqr(Tg[j-1])));{Tendency of surface temperature}
if Lv>0 then FPW:=Evaporationenergyflux/Lv;{tendency of PW due to evaporation}{kg/(s m^2)}
if Lv=0 then FPW:=0.0;{tendency of PW due to evaporation}{kg/(s m^2)}


{end tendency at the earth's surface}
{*********************************************************************}

{BEGIN: Calculate vertical distribution of water vapour and calculate condensation and latent heat release}

{First: calculate the thickness of the subcloud+cloud layer, DH}
DH:=0;
for n:=KMAX downto cloudTopLevel do begin
dz[n]:=(Rd*T[n,0]/g)*(ln(p[n]+0.5*dp)-ln(p[n]-0.5*dp));
if (n=KMAX) or (n=cloudTopLevel) then DH:=DH+(0.5*dz[n])
else
DH:=DH+(dz[n])
end;

{Calculate the lapse rate within this layer}
if cloudTopLevel=KMAX then LR:=g/cp;
if cloudTopLevel<KMAX then LR:=-(T[n-1,j-1]-T[KMAX,j-1])/DH;
{Use this to calculate Hv}
{Hv:=Rv*sqr(Tg[j-1])/(Lv*LR);if Hv>3000.0 then Hv:=3000.0;}{DYNAMIC WATER VAPOUR SCALE HEIGHT}


{Rovg:=PW[j-1]/Hv0;}Rovg:=PW[j-1]/Hv;{Density of water vapour at the ground}
if Rovg>0.0 then eg:=Rovg*Rv*Tg[j-1];{vapour pressure at the ground}
if Rovg=0.0 then eg:=0;{vapour pressure at the ground}

SATVAPOURPRESSURE(Tg[j-1]);
RHg:=eg/es;

CloudTopLevel:=KMAX;
{>>}for n:=KMAX downto 1 do begin
if n=KMAX then Rov[n]:=Rovg*exp(((Rd*T[KMAX,j-1]/(g*Hv))*(ln(ps-(0.5*dp))-ln(ps))));
if n<KMAX then Rov[n]:=Rov[n+1]*exp((Rd*0.5*(T[n,j-1]+T[n+1,j-1])/(g*Hv))*(ln(p[n])-ln(p[n+1])));
e:=Rov[n]*Rv*T[n,j-1];{vapour pressure at level n}
SATVAPOURPRESSURE(T[n,j-1]);
RH[n]:=e/es;{relative humidity at level n}
qv[n]:=e*Rd/(Rv*(p[n]));
qvs[n]:=es*Rd/(Rv*(p[n]));

dqv[n]:=(qv[n]-(RHc[n]*qvs[n]))*dt/Tauc;{tendency of mixing ratio towards a predefined maximum value, as in SPEEDY}
If (CC=0.0) or (CC=1.0) then dqv[n]:=0.0;{if no clouds then no precipitation; if there is nevertheless supersaturation, i.e. dqv>0, we put CC=1 (see above: !!!!)}

if (dqv[n]>0) and (CC>0) then begin
CloudTopLevel:=n;
FLatentHeat[n]:=Lv*dqv[n]*dp/(g*dt);
FPW:=FPW-(dqv[n]*dp/(dt*g));{tendency is PW due to condensation: SINK of WATER VAPOUR}
end;{if}
if dqv[n]<=0 then begin
dqv[n]:=0.0;
FLatentHeat[n]:=0.0;
FPW:=FPW;
end;{if}
if CC=0 then begin
dqv[n]:=0.0;
FLatentHeat[n]:=0.0;
FPW:=FPW;
end;{if}

{>>}end;{n}
Rov[1]:=0.0;

if soilwetness=0.0 then FPW:=0.0;
{END: Calculate vertical distribution of water vapour and calculate condensation and latent heat release}


{compute tendencies in each layer (n) due to infrared radiation received from each other layer (m) and from the Earth's surface}
sumSR:=0;{#####}

for n:=1 to KMAX do begin{>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>begin radiation budget of each layer}


{#L#}

{radiation from the atmosphere}
IRa[n]:=0;
if n=1 then OLRTOA:=0;
for LWbeam:=1 to 2 do begin
for m:=1 to KMAX do begin{summation over all layers, m}

if m=n then goto 1;

{PLANCK1(T[m,j-1]);}
bin:=trunc(T[m,j-1]);
for band:=1 to bandmax do begin
if LWbeam=1 then Emission:=(1-CC)*LWEnergy[bin,band]*eps[m,band,LWbeam]*sigma*sqr(sqr(T[m,j-1]));{Planck}
if LWbeam=2 then Emission:=(CC)*LWEnergy[bin,band]*eps[m,band,LWbeam]*sigma*sqr(sqr(T[m,j-1]));{Planck}

{Upwelling: IR radiation from below:}
if m>n then begin
for mm:=m-1 downto n+1 do begin
Emission:=Emission*(1-eps[mm,band,LWbeam]);
end;
end;{if}
{Downwelling: IR radiation from above:}
if m<n then begin
for mm:=m+1 to n-1 do begin
Emission:=Emission*(1-eps[mm,band,LWbeam]);
end;
end;{if}
if n=1 then OLRTOA:=OLRTOA+Emission*(1-eps[n,band,LWbeam]);
IRa[n]:=IRa[n]+(Emission*eps[n,band,LWbeam]);{sum of contributions of all layers multiplied by absorption coeficient of layer n}
end;{band}


1:
end;{m}end;{LWbeam}

{net absorbed radiation from the Earth's surface (by layer n) }
{emission coefficient Earth's surface does not depend on spectral band}
Ua[n]:=0;
if n=1 then UgOut:=0;
for LWbeam:=1 to 2 do begin
{PLANCK1(Tg[j-1]);}
bin:=trunc(Tg[j-1]);
for band:=1 to bandmax do begin
if LWbeam=1 then Emission:=(1-CC)*LWEnergy[bin,band]*(((1-epsg)*CounterRadiation)+(epsg*sigma*sqr(sqr(Tg[j-1]))));{^^}{term with B is atmsopheric radiation that is reflected by Earth's surface}
if LWbeam=2 then Emission:=CC*LWEnergy[bin,band]*(((1-epsg)*CounterRadiation)+(epsg*sigma*sqr(sqr(Tg[j-1]))));{^^}{term with B is atmsopheric radiation that is reflected by Earth's surface}
{Emission coefficient of Earth is epsg is all bands; 
there is also reflection of radiation that is received from the atmosphere: reflection coefficient is (1-epsg)}
for m:=(KMAX) downto (n+1) do begin
Emission:=Emission*(1-eps[m,band,LWbeam]);
end;{m}
if n=1 then OLRTOA:=OLRTOA+Emission*(1-eps[n,band,LWbeam]);
if n=1 then UgOut:=UgOut+Emission*(1-eps[n,band,LWbeam]);
Ua[n]:=Ua[n]+Emission*eps[n,band,LWbeam];
end;{band}end;{LWbeam}

{#L#}


{#S#}
for channel:=1 to channelmax do begin

SRD:=gamma[channel]*Q;

CT:=0;
for m:=1 to n-1 do begin
if (qL[m]>0.0) and (CT=0) then begin
{then cloud top}
CT:=1;SRD:=SRD*(1-albedoCT);{Solar radiation partly reflected at cloud top}
end;{at cloud top}
SRD:=SRD*exp(-OpticalPathSW[m,channel]);{part of downward beam that passes through layer m}
end;{m}
SRD:=SRD*(1-exp(-OpticalPathSW[n,channel]));{part of downward beam that is absorbed by layer n}


{radiation from the sun which reaches Earth's surface and is reflected at Earth's surface and absorbed in layer n }

SRU:=gamma[channel]*Q;

{Downward beam}
CT:=0;
nCT:=0;
for m:=1 to KMAX do begin
if (qL[m]>0.0) and (CT=0) then begin
{then cloud top}
CT:=1;nCT:=m;SRU:=SRU*(1-albedoCT);
end;{if}{at cloud top}
SRU:=SRU*exp(-OpticalPathSW[m,channel]);
end;{m}
SRU:=(alfag)*SRU;{radiation that is reflected at Earth's surface}
{Upward beam}
for m:=KMAX downto n+1 do begin
SRU:=SRU*exp(-OpticalPathSW[m,channel]);
end;{m}
SRU:=SRU*(1-exp(-OpticalPathSW[n,channel]));{This part of upward beam of reflected Solar radiation is absorbed by layer n}

{radiation from the sun which reaches cloud top and is reflected at cloud top and absorbed in layer above cloud top }
if (nCT>0) and (n<nCT) then begin

SRC:=gamma[channel]*Q;

{Downward beam}
CT:=0;
for m:=1 to nCT-1 do begin
SRC:=SRC*exp(-OpticalPathSW[m,channel]);{above cloud top}
end;{now the beam is at the top of the first (highest) cloud layer, n=nCT}
SRC:=SRC*(albedoCT);{beam is reflected}
for m:=nCT-1 downto n+1 do begin
SRC:=SRC*exp(-OpticalPathSW[m,channel]);{above cloud top}
end;
SRC:=SRC*(1-exp(-OpticalPathSW[n,channel]));{This part of upward beam of reflected Solar radiation is absorbed by layer n}
end;{if nCT and n}
if (nCT>0) and (n>=nCT) then SRC:=0;{level n is below cloud top, then this level does not receive reflected radiation}
if (nCT=0) then SRC:=0;{no clouds}

if channel=1 then SR[n]:=SRD+SRU+SRC;
if (channel>1) then SR[n]:=SR[n]+SRD+SRU+SRC;

sumSR:=sumSR+SRD+SRU+SRC;{#####}

end;{channel}{#S#}

if n<KMAX then begin
H[n]:=k[n]*((T[n+1,j-1]*exp(ln(ps/p[n+1])*kappa))-(T[n,j-1]*exp(ln(ps/p[n])*kappa)));
if H[n]<0 then H[n]:=0;
end;

{#L#}
{Planck emission}
{PLANCK1(T[n,j-1]);}
bin:=trunc(T[n,j-1]);
Emission:=0;
for LWbeam:=1 to 2 do begin
for band:=1 to bandmax do begin
if LWbeam=1 then Emission:=Emission+((1-CC)*LWEnergy[bin,band]*eps[n,band,LWbeam]*sigma*sqr(sqr(T[n,j-1])));{Planck emission}
if LWbeam=2 then Emission:=Emission+(CC*LWEnergy[bin,band]*eps[n,band,LWbeam]*sigma*sqr(sqr(T[n,j-1])));{Planck emission}
end;end;
{#L#}
if n=1 then OLRTOA:=OLRTOA+Emission;

if n>1 then F[n]:=(KMAX*g/(cp*ps))*(FLatentHeat[n]+Ua[n]+IRa[n]+SR[n]+H[n]-H[n-1]-(2*Emission));
if n=1 then F[n]:=(KMAX*g/(cp*ps))*(Ua[n]+IRa[n]+SR[n]+H[n]-(2*Emission));

end;{n}{<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<end radiation budget of each layer}

{Compute total absorbed radiation by atmosphere}
Sabsorbed:=0;
for n:=1 to KMAX do begin
Sabsorbed:=Sabsorbed+SR[n];
end;

if j=1 then PW[j]:=PW[0]+(FPW*dt/2);
if j=2 then PW[j]:=PW[0]+(FPW*dt/2);
if j=3 then PW[j]:=PW[0]+(FPW*dt);
if j=4 then PW[j]:=PW[0]-(FPW*dt/2);

if j=1 then Tg[j]:=Tg[0]+(Fg*dt/2);
if j=2 then Tg[j]:=Tg[0]+(Fg*dt/2);
if j=3 then Tg[j]:=Tg[0]+(Fg*dt);
if j=4 then Tg[j]:=Tg[0]-(Fg*dt/2);

for n:=1 to KMAX do begin
if j=1 then T[n,j]:=T[n,0]+(F[n]*dt/2);
if j=2 then T[n,j]:=T[n,0]+(F[n]*dt/2);
if j=3 then T[n,j]:=T[n,0]+(F[n]*dt);
if j=4 then T[n,j]:=T[n,0]-(F[n]*dt/2);
end;

end;{j}{end substeps in R-K scheme}

{final step}
PW[0]:=(PW[1]+(2*PW[2])+PW[3]-PW[4])/3;
Tg[0]:=(Tg[1]+(2*Tg[2])+Tg[3]-Tg[4])/3;

for n:=1 to KMAX do begin
T[n,0]:=(T[n,1]+(2*T[n,2])+T[n,3]-T[n,4])/3;
end;{n}

{prescribe the minimum specific humidity (important for the stratosphere)}
q0:=qStrat;
PW0:=q0*ps/g;
if PW0>PW[0] then begin
PW0:=0.0;
q0:=0.0;
end;

{BEGIN: Calculate vertical distribution of water vapour}
Rovg0:=ps*q0/(Rd*Tg[0]);
PWprime:=PW[0]-PW0;
Rovgprime:=PWprime/Hv;
Rovg:=Rovg0+Rovgprime;{Density of water vapour at the ground}
eg:=Rovg*Rv*Tg[0];{vapour pressure at the ground}
SATVAPOURPRESSURE(Tg[0]);
RHg:=eg/es;

{#}
for n:=KMAX downto 1 do begin
Rov0:=p[n]*q0/(Rd*T[n,0]);

if n=KMAX then Rovprime:=Rovgprime*exp(((Rd*T[KMAX,0]/(g*Hv))*(ln(p[KMAX])-ln(ps))));
if n<KMAX then Rovprime:=Rovprime*exp((Rd*0.5*(T[n,0]+T[n+1,0])/(g*Hv))*(ln(p[n])-ln(p[n+1])));

Rov[n]:=Rov0+Rovprime;
if Rovg=0.0 then Rov[n]:=0.0;

e:=Rov[n]*Rv*T[n,0];{vapour pressure at level n}
SATVAPOURPRESSURE(T[n,0]);
RH[n]:=e/es;{relative humidity at level n}
qv[n]:=e*Rd/(Rv*(p[n]));
qvs[n]:=es*Rd/(Rv*(p[n]));
end;

{END: Calculate vertical distribution of water vapour}

{Compute total absorbed radiation by atmosphere}
Sabsorbed:=0;
for n:=1 to KMAX do begin
Sabsorbed:=Sabsorbed+SR[n];
end;

{output}
if nt=trunc(nt/tf)*tf then begin
write(time/(24*3600):8:1);
write(f2,time/(24*3600):11:3);
for n:=1 to KMAX do begin
{write(T[n,0]:7:1);}
write(f2,T[n,0]:7:1);
end;write(Tg[0]:7:1);write(albedo:9:2);write(OLRTOA:10:2);write(Sabsorbed+Sgtot:13:2,H[KMAX]:7:1,PW[0]:21:1,RHg:7:2,CC:7:2,24*3600*(FPW-(Evaporationenergyflux/Lv)):12:2,24*3600*Evaporationenergyflux/Lv:14:2);
write(f2,Tg[0]:7:1);write(f2,albedo:7:2);write(f2,OLRTOA:9:2);write(f2,Sabsorbed+Sgtot:13:2,H[KMAX]:7:1,PW[0]:21:1,RHg:7:2,CC:7:2,24*3600*(FPW-(Evaporationenergyflux/Lv)):12:2,24*3600*Evaporationenergyflux/Lv:14:2);
writeln;writeln(f2);
end;

end;{time loop}

{Integrated liquid water}
LWP:=0;
for n:=1 to KMAX do begin
LWP:=LWP+(dp*qL[n]/g);
end;

for n:=KMAX downto 1 do begin
Rov0:=p[n]*q0/(Rd*T[n,0]);

if n=KMAX then Rovprime:=Rovgprime*exp(((Rd*T[KMAX,0]/(g*Hv))*(ln(p[KMAX])-ln(ps))));
if n<KMAX then Rovprime:=Rovprime*exp((Rd*0.5*(T[n,0]+T[n+1,0])/(g*Hv))*(ln(p[n])-ln(p[n+1])));

Rov[n]:=Rov0+Rovprime;
if Rovg=0.0 then Rov[n]:=0.0;

e:=Rov[n]*Rv*T[n,0];{vapour pressure at level n}
SATVAPOURPRESSURE(T[n,0]);
RH[n]:=e/es;{relative humidity at level n}
qv[n]:=e*Rd/(Rv*(p[n]));
qvs[n]:=es*Rd/(Rv*(p[n]));
rs[n]:=es*Rd/(Rv*(p[n]-es));
Theta[n]:=(T[n,0]*exp(ln(ps/p[n])*kappa));
end;

writeln(f1,'pressure[hPa] Temperature[K] potTemp[K]    HeatFlux    qv      RH[%]  RHc    FLatentHeat[K/day]      qL');
writeln('  n  pressure[hPa] Temperature[K] potTemp[K]    HeatFlux    qv      RH[%]  RHc    FLatentHeat[K/day]       qL');

for n:=1 to KMAX do begin
writeln(f1,p[n]/100:10:1,T[n,0]:15:2,Theta[n]:15:2,H[n]:12:5,qv[n]:9:5,RH[n]*100:7:1,RHc[n]*100:7:1,24*3600*(KMAX*g/(cp*ps))*FLatentHeat[n]:10:2,qL[n]:12:7);
writeln(n:4,p[n]/100:10:1,T[n,0]:15:2,Theta[n]:15:2,H[n]:12:5,qv[n]:9:5,RH[n]*100:7:1,RHc[n]*100:7:1,24*3600*(KMAX*g/(cp*ps))*FLatentHeat[n]:10:2,qL[n]:12:7);
end;
e:=Rovg*Rv*Tg[0];{vapour pressure at level n}
SATVAPOURPRESSURE(Tg[0]);
rg:=e*Rd/(Rv*(ps-e));

writeln('  n  pressure[hPa] Temperature[K] potTemp[K]    HeatFlux    r      RH[%]  PW[kg/m^2]  Evaporation_rate[kg/(m^2.day)] Precipitation_rate[kg/(m^2.day)]    ');
if Lv>0 then writeln(f1,ps/100:10:1,Tg[0]:15:2,Tg[0]:15:2,-H[KMAX]:12:5,rg:9:5,RH[n]*100:7:1,PW[0]:8:1,24*3600*Evaporationenergyflux/Lv:20:3,24*3600*(FPW-(Evaporationenergyflux/Lv)):20:3);
if Lv>0 then writeln((KMAX+1):4,ps/100:10:1,Tg[0]:15:2,Tg[0]:15:2,-H[KMAX]:12:5,rg:9:5,RH[n]*100:7:1,PW[0]:8:1,24*3600*Evaporationenergyflux/Lv:20:3,24*3600*(FPW-(Evaporationenergyflux/Lv)):20:3);

if Lv=0 then writeln(f1,ps/100:10:1,Tg[0]:15:2,Tg[0]:15:2,-H[KMAX]:12:5,rg:9:5,RH[n]*100:7:1,PW[0]:8:1,FPW:20:3,FPW:20:3);
if Lv=0 then writeln((KMAX+1):4,ps/100:10:1,Tg[0]:15:2,Tg[0]:15:2,-H[KMAX]:12:5,rg:9:5,RH[n]*100:7:1,PW[0]:8:1,FPW:20:3,FPW:20:3);



writeln;
writeln(' n=layer index;');
writeln(' T=temperature of layer;');
writeln(' Th=potential temperature of layer;');
writeln(' The=equivalent potential temperature of layer;');
writeln(' H=convective heat flux from layer below(n+1)');
writeln(' U=absorbed terrestrial radiation;');
writeln(' IR=absorbed radiation coming from other layers in the atmosphere');
writeln(' B=emitted radiation by layer n');
writeln(' L=latent heat release in layer n/Energy flux due to evaporation at the ground');
writeln(' F=Total diabatic heating');
writeln(' NBV=Brunt-Vaisala frequency (if NBV=0: Neutral or unstable stratification)');
writeln;
writeln('   n  p[hPa]  T[K]  Th[K] SW[W^m^-2] U[W^m^-2] IR[W^m^-2] B[W^m^-2] H[W^m^-2] L[W^m^-2] F[KDay^-1] NBV[s^-1]   RH[%]  qO3[ppm]   qL[ppm]     qv[per mil]');
writeln(f3,'   n  p[hPa]  T[K]  Th[K] SW[W^m^-2] U[W^m^-2] IR[W^m^-2] B[W^m^-2] H[W^m^-2] L[W^m^-2] F[KDay^-1] NBV[s^-1]   RH[%]  qO3[ppm]   qL[ppm]     qv[per mil]');


{compute Brunt-V�is�l� frequency, NBV}
for n:=1 to KMAX do begin


if (n>1) and (n<KMAX) then begin
thetaupper:=(0.5*(T[n-1,0]+T[n,0])*exp(ln(ps/(p[n]-(0.5*dp)))*kappa));
thetalower:=(0.5*(T[n,0]+T[n+1,0])*exp(ln(ps/(p[n]+(0.5*dp)))*kappa));
end;

if (n=KMAX) then begin
thetaupper:=(0.5*(T[n-1,0]+T[n,0])*exp(ln(ps/(p[n]-(0.5*dp)))*kappa));
thetalower:=Tg[0];
end;

if (n>1) and (n<KMAX) then dEXNER:=cp*((0.5*(T[n-1,0]+T[n,0])/thetaupper)-(0.5*(T[n,0]+T[n+1,0])/thetalower));
if (n=KMAX) then dEXNER:=cp*((0.5*(T[n-1,0]+T[n,0])/thetaupper)-1);
if thetaupper>thetalower then NBV[n]:=(g/(0.5*(thetaupper+thetalower)))*sqrt(-(thetaupper-thetalower)/dEXNER);
if thetaupper<=thetalower then NBV[n]:=0.0;
end;
NBV[1]:=NBV[2];

SumLatentHeat:=0.0;
for n:=1 to KMAX do begin
SumLatentHeat:=SumLatentHeat+FLatentHeat[n];end;

for n:=1 to KMAX do begin

{PLANCK1(T[n,0]);}
bin:=trunc(T[n,0]);
Emission:=0;
for LWbeam:=1 to 2 do begin
for band:=1 to bandmax do begin
if LWbeam=1 then Emission:=Emission+((1-CC)*LWEnergy[bin,band]*2*eps[n,band,LWbeam]*sigma*sqr(T[n,0])*sqr(T[n,0]));
if LWbeam=2 then Emission:=Emission+(CC*LWEnergy[bin,band]*2*eps[n,band,LWbeam]*sigma*sqr(T[n,0])*sqr(T[n,0]));
end;{band}end;{LWbeam}





writeln(n:4,p[n]/100:7:1,T[n,0]:7:1,Theta[n]:7:1,SR[n]:11:2,Ua[n]:10:2,IRa[n]:11:2,-Emission:10:2,H[n]:10:2,FLatentHeat[n]:10:2,3600*24*F[n]:10:3,NBV[n]:10:4,RH[n]*100:7:1,1000000*qO3[n]:8:2,1000000*qL[n]:11:2,1000*qv[n]:11:2);{Planck}
writeln(f3,n:4,p[n]/100:7:1,T[n,0]:7:1,Theta[n]:7:1,SR[n]:11:2,Ua[n]:10:2,IRa[n]:11:2,-Emission:10:2,H[n]:10:2,FLatentHeat[n]:10:2,3600*24*F[n]:10:3,NBV[n]:10:4,RH[n]*100:7:1,1000000*qO3[n]:8:2,1000000*qL[n]:11:2,1000*qv[n]:11:2);{Planck}
end;
writeln('   g',ps/100:7:1,Tg[0]:7:1,Tg[0]:7:1,Sgtot:11:2,-(epsg*sigma*sqr(Tg[0])*sqr(Tg[0])):10:2,B:11:2,'     99999',-H[KMAX]:10:2,-Evaporationenergyflux:10:2,3600*24*Fg:10:3,'     99999',RHg:7:2,'    9999','       0.00','       9999');{Planck}
writeln(f3,'   g',ps/100:7:1,Tg[0]:7:1,Tg[0]:7:1,Sgtot:11:2,-(epsg*sigma*sqr(Tg[0])*sqr(Tg[0])):10:2,B:11:2,'     99999',-H[KMAX]:10:2,-Evaporationenergyflux:10:2,3600*24*Fg:10:3,'     99999',RHg:7:2,'    9999','       0.00','       9999');{Planck}
writeln('');

Bowen:=Soilwetness*H[KMAX]/Evaporationenergyflux;

writeln('');
writeln('NOTE: Relative humidity is calculated from the total water content of the atmosphere');
writeln('      Part of this water is in the form of "cloud water"');

writeln('');
writeln('Run:',nrun:5);
writeln('tmax:',tmax:6,' days');
writeln('SUN:',SUN:5,' (0:constant insolation; 1:seasonal cycle)');
writeln('phi:',phi:8:1,' degrees');
writeln('Sabsorbed-A=',Sabsorbed:9:2,' W m^-2',' (net absorbed Solar radiation by the atmosphere)',' (fraction of incoming irradiance below highest level=',Sabsorbed/(SumGamma*S0/4):7:3,')');
writeln('Sabsorbed-A+g=',Sabsorbed+Sgtot:9:2,' W m^-2',' (net absorbed Solar radiation by the planet)',' (fraction of incoming irradiance below highest level=',(Sabsorbed+Sgtot)/(SumGamma*S0/4):7:3,')');
writeln('Long-wave back-radiation',CounterRadiation:10:2);
writeln('Net radiation at the ground=',(Sgtot+B-(epsg*sigma*sqr(Tg[0])*sqr(Tg[0]))):9:2,' W m^-2');
if ozone=1 then writeln('Fraction of incoming irradiance in ultraviolet that is absorbed above highest level=',(1-SumGamma):7:3);
writeln('Total latent heat released=',SumLatentHeat:9:2,' W m^-2');
writeln('U-TOA=',UgOut:10:2,' W m^-2',' (net transmitted radiation from the surface of planet)');
writeln('U-TOA=',UgOut/(((1-epsg)*CounterRadiation)+(epsg*sigma*sqr(Tg[0])*sqr(Tg[0]))):9:3,' (net fraction transmitted radiation from the surface of planet)');{Planck}
writeln('OLR-TOA=',OLRTOA:9:2,' W m^-2');
writeln('Precipitable water=',PW[0]:11:2,'[ kg m^-2]');
if Lv>0 then writeln('Evaporation rate=',24*3600*Evaporationenergyflux/Lv:14:2,'[ kg m^-2 d^-1]');
if Lv>0 then writeln('Precipitation rate=',24*3600*(FPW-(Evaporationenergyflux/Lv)):12:2,'[ kg m^-2 d^-1]');
if Lv=0 then writeln('Evaporation rate=0');
if Lv=0 then writeln('Precipitation rate=0');
writeln('Cloud cover=',CC:10:2);
if CC>0 then writeln('Cloud top level=',p[CloudTopLevel]/100:10:2,' [hPa]');
if CC=0 then writeln('Cloud top level:',' no clouds');
{writeln('Weighted albedo=',albedo:10:2);}
if ozone=0 then writeln('Planetary or Bond albedo=',1-((Sabsorbed+Sgtot)/(S0/4)):7:2);
if ozone=1 then writeln('Planetary or Bond albedo=',1-((Sabsorbed+Sgtot+((1-SumGamma)*S0/4))/(S0/4)):7:2);
writeln('Bowen=',Bowen:12:4,' (Evaporative Ratio=',ER:12:4,')');
writeln('Ozone=',Ozone:5,' (0: no ozone; 1: ozone prescribed according to observations)');
writeln('Hv=',Hv:10:1,' [m]');
writeln('epsg=',epsg:10:2);
writeln('pCO2=',pCO2:10:6,' ppmv');
writeln('alfag=',alfag:10:3);
writeln('alfac=',alfac:10:2);
writeln('Tauc=',Tauc/3600:10:3,' hrs');
writeln('LWP=',1000*LWP:10:3,' g m^-2');

for channel:=1 to channelmax do begin
writeln('SW-absorption crosssection water vapour, channel',channel:5,sigmaSWwv[channel]:10:4,' m^2 kg^-1');
end;

writeln(f3,'');
writeln(f3,'Run:',nrun:5);
writeln(f3,'tmax:',tmax:6,' days');
writeln(f3,'SUN:',SUN:5,' (0:constant insolation; 1:seasonal cycle)');
writeln(f3,'phi:',phi:8:1,' degrees');
writeln(f3,'Sabsorbed-A=',Sabsorbed:9:2,' W m^-2',' (net absorbed Solar radiation by the atmosphere)',' (fraction of incoming irradiance below highest level=',Sabsorbed/(SumGamma*S0/4):7:3,')');
writeln(f3,'Sabsorbed-A+g=',Sabsorbed+Sgtot:9:2,' W m^-2',' (net absorbed Solar radiation by the planet)',' (fraction of incoming irradiance below highest level=',(Sabsorbed+Sgtot)/(SumGamma*S0/4):7:3,')');
writeln(f3,'Long-wave back-radiation',CounterRadiation:10:2);
writeln(f3,'Net radiation at the ground=',(Sgtot+B-(epsg*sigma*sqr(Tg[0])*sqr(Tg[0]))):9:2,' W m^-2');
if ozone=1 then writeln(f3,'Fraction of incoming irradiance in ultraviolet that is absorbed above highest level=',gamma[1]:7:3);
writeln(f3,'Total latent heat released=',SumLatentHeat:9:2,' W m^-2');
writeln(f3,'U-TOA=',UgOut:10:2,' W m^-2',' (net transmitted radiation from the surface of planet)');
writeln(f3,'U-TOA=',UgOut/(((1-epsg)*CounterRadiation)+(epsg*sigma*sqr(Tg[0])*sqr(Tg[0]))):9:3,' (net fraction transmitted radiation from the surface of planet)');{Planck}
writeln(f3,'OLR-TOA=',OLRTOA:9:2,' W m^-2');
writeln(f3,'Precipitable water=',PW[0]:11:2,'[ kg m^-2]');
if Lv>0 then writeln(f3,'Evaporation rate=',24*3600*Evaporationenergyflux/Lv:14:2,'[ kg m^-2 d^-1]');
if Lv>0 then writeln(f3,'Precipitation rate=',24*3600*(FPW-(Evaporationenergyflux/Lv)):12:2,'[ kg m^-2 d^-1]');
if Lv=0 then writeln(f3,'Evaporation rate=0');
if Lv=0 then writeln(f3,'Precipitation rate=0');
writeln(f3,'Cloud cover=',CC:10:2);
if CC>0 then writeln(f3,'Cloud top level=',p[CloudTopLevel]/100:10:2,' [hPa]');
if CC=0 then writeln(f3,'Cloud top level:',' no clouds');
writeln(f3,'Weighted albedo=',albedo:10:2);
if ozone=0 then writeln(f3,'Planetary or Bond albedo=',1-((Sabsorbed+Sgtot)/(S0/4)):7:2);
if ozone=1 then writeln(f3,'Planetary or Bond albedo=',1-((Sabsorbed+Sgtot+((1-SumGamma)*S0/4))/(S0/4)):7:2);
writeln(f3,'Bowen=',Bowen:12:4,' (Evaporative Ratio=',ER:12:4,')');
writeln(f3,'Ozone=',Ozone:5,' (0: no ozone; 1: ozone prescribed according to observations)');
writeln(f3,'Hv=',Hv:10:1,' [m]');
writeln(f3,'epsg=',epsg:10:2);
writeln(f3,'pCO2=',pCO2:10:6,' ppmv');
writeln(f3,'alfag=',alfag:10:3);
writeln(f3,'alfac=',alfac:10:3);
writeln(f3,'Tauc=',Tauc/3600:10:3,' hrs');
writeln(f3,'LWP=',1000*LWP:10:3,' g m^-2');

for channel:=1 to channelmax do begin
writeln(f3,'SW-absorption crosssection water vapour, channel',channel:5,sigmaSWwv[channel]:10:4,' m^2 kg^-1');
end;

writeln(f3,'LW absorption crossections,      O3       CO2  H2O-vapour  cloud');
writeln('LW absorption crossections,      O3       CO2  H2O-vapour  cloud');
for band:=1 to bandmax do begin
write('band=',band:4,sigmaO3[band]:26:3,sigmaCO2[band]:10:3,sigmav[band]:10:3,sigmaL[band]:10:3);writeln;
write(f3,'band=',band:4,sigmaO3[band]:26:3,sigmaCO2[band]:10:3,sigmav[band]:10:3,sigmaL[band]:10:3);writeln(f3);
end;

close(f1);close(f2);close(f3);

end;{runnumber}
2:
writeln('END RADIATIVE-CONVECTIVE MODEL (RCM)');
writeln('===================================');
writeln('press RETURN to finish');
readln;
end.