#!/usr/bin/perl5
#ebm.pl
#Energy Balance Model
# Robert R. Gotwals, September 11, 1997
#
#Initialize Variables
$A=204;
$B=2.17;
$C=3.87;
$Aice=0.62;
$smalla=0.3;
$Tcrit=-10;
$SCfrac=1;
$SC=1370;
$iterations=24;
@Zones=("80-90","70-80","60-70","50-60","40-50","30-40","0-30","0-20","0-10");
@ZoneLat=(85,75,65,55,45,35,25,15,6);
@SunWt=(0.5,0.531,0.624,0.77,0.892,1.021,1.12,1.189,1.219);
@Init_T=(-16.9,-12.3,-5.1,2.2,8.8,16.2,22.9,26.1,26.4);

#Calculate Init_a
for ($lat=0;$lat<=8;$lat+=1){
if ($Init_T[$lat] < $Tcrit) {
		$Init_a[$lat]=$Aice
} else {
			$Init_a[$lat]=$smalla
}
}
#
#  Calculate sumcos from Cell L17
$sumcos=0;
for ($lat=0;$lat<=8;$lat+=1){
$coz[$lat]=cos($ZoneLat[$lat]*3.1415/180);
$sumcos=$sumcos+$coz[$lat];
}
#Calculate R_in[$lat]
for ($lat=0;$lat<=8;$lat+=1){
				$R_in[$lat]=$SC/4*$SCfrac*$SunWt[$lat];
}
#
#STEP ONE
#for ($step=1;$step<=1;$step +=1){
#Calculate T-cos
		for ($lat=0;$lat<=8;$lat +=1){
						$Tcos[$lat]=$Init_T[$lat]*$coz[$lat];
}
# Calculate Mean_T
					$Mean_T[1]=0;
				for ($lat=0;$lat<=8;$lat +=1){
								$Mean_T[1]=$Mean_T[1]+$Tcos[$lat];						
				}
$Mean_T[1]=$Mean_T[1]/$sumcos;
#
#Calculate Temp and Albedo for Step one
		for ($lat=0;$lat<=8;$lat +=1){
						$Temp[$lat]=($R_in[$lat]*(1-$Init_a[$lat])+$C*$Mean_T[1]-$A)/($B+$C);
      if ($Temp[$lat] < $Tcrit) {
											$Albedo[$lat]=$Aice
							} else {
										$Albedo[$lat]=$smalla
							}
		}
#
#
#Calculate intermediate steps
for ($step=1;$step<=$iterations;$step +=1){
#Calculate T-cos
		for ($lat=0;$lat<=8;$lat +=1){
						$Tcos[$lat]=$Temp[$lat]*$coz[$lat];
}
# Calculate Mean_T
					$Mean_T[$step]=0;
				for ($lat=0;$lat<=8;$lat +=1){
								$Mean_T[$step]=$Mean_T[$step]+$Tcos[$lat];						
				}
$Mean_T[$step]=$Mean_T[$step]/$sumcos;
#
#Calculate Temp and Albedo for Step one
		for ($lat=0;$lat<=8;$lat +=1){
						$Temp[$lat]=($R_in[$lat]*(1-$Albedo[$lat])+$C*$Mean_T[$step]-$A)/($B+$C);
      if ($Temp[$lat] < $Tcrit) {
											$Albedo[$lat]=$Aice
							} else {
										$Albedo[$lat]=$smalla
							}
		}
}

#FINAL STEP
$GMT=$Mean_T[$iterations];
for ($lat=0;$lat<=8;$lat+=1) {
						$Final_T[$lat]=($R_in[$lat]*(1-$Init_a[$lat])+$C*$GMT-$A)/($B+$C);
			      if ($Final_T[$lat] < $Tcrit) {
											$Final_a[$lat]=$Aice
							} else {
										$Final_a[$lat]=$smalla
							}
$R_out[$lat]=$A+$B*$Final_T[$lat];
print "$ZoneLat[$lat], $Final_a[$lat], $R_in[$lat], $R_out[$lat], $Final_T[$lat]\n";					
}