'Energy Balance Model
'copied from pages 77-81 of "A Climate Modelling Primer" second edition
' by K. McGuffie & A. Henderson-Sellers
'converted form BASICA to QBasic by Daniel Wedul


'declare subs
DECLARE SUB bad ()
DECLARE SUB intro ()
DECLARE SUB waitforspace ()
DECLARE SUB slow ()
DECLARE SUB mainmenu ()
DECLARE SUB albedomenu ()
DECLARE SUB lattrans ()
DECLARE SUB longwave ()
DECLARE SUB runmodel ()
DECLARE SUB alicetemp ()
DECLARE SUB alalbedo ()
DECLARE SUB alground ()
DECLARE SUB bad ()
DECLARE SUB showresults ()
DECLARE SUB initialize ()
DECLARE SUB save ()

'make all variables common
COMMON SHARED LATICE, A, B, C, AICE, TCRIT, IN, PIBY2, sx

DIM SHARED LATZ$(9), S(9), ALAND(9), TSTARTorig(9), TEMP(9), ALBEDO(9)
DIM SHARED ALANDorig(9), TSTART(9), TM(9)

'initialize variables
DEF FNR (x) = INT(100 * x) / 100

'constants
IN = 3.14159 / 18
PIBY2 = 3.14159 / 2

'arrays
FOR i = 1 TO 9
   READ LATZ$(i)
NEXT i
FOR lat = 1 TO 9
   READ S(lat)
NEXT lat
FOR i = 1 TO 9
   READ ALANDorig(i)
NEXT i
FOR lat = 1 TO 9
   READ TSTARTorig(lat)
NEXT lat

'parameters
initialize

'start of program

intro
mainmenu

'end of program



'begin data
DATA "80-90", "70-80", "60-70", "50-60", "40-50", "30-40", "20-30", "10-20", " 0-10"
DATA 0.5, 0.531, 0.624, 0.77, 0.892, 1.021, 1.12, 1.189, 1.219

DATA 0.5, 0.5, 0.452, 0.407, 0.357, 0.309, 0.272, 0.248, 0.254

DATA -16.9, -12.3, -5.1, 2.2, 8.8, 16.2, 22.9, 26.1, 26.4
'end data

SUB alalbedo
'shared AICE
CLS
PRINT
PRINT "     The default albedo of the ice is 0.62."
PRINT
DO
   INPUT "     What is the new value you want for this albedo ? ", AICE
   IF AICE < 0 OR AICE > .99 THEN bad ELSE EXIT DO
LOOP

END SUB

SUB albedomenu
DO
   CLS
   PRINT "* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *"
   PRINT "               P A R A M E T E R I Z A T I O N   O F"
   PRINT "                       A L B E D O"
   PRINT "* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *"
   PRINT
   PRINT
   PRINT "          There are three things which you may alter"
   PRINT
   PRINT "            1.   The temperature at which the surface"
   PRINT "                 becomes ice covered."
   PRINT
   PRINT "            2.   The albedo of this ice covered surface."
   PRINT
   PRINT "            3.   The albedo of the underlying ground."
   PRINT
   PRINT "            4.   Return to the main menu."
   PRINT
   PRINT " Choose which one you want"
   DO
      al$ = INPUT$(1)
   LOOP UNTIL al$ = "1" OR al$ = "2" OR al$ = "3" OR al$ = "4"
   SELECT CASE al$
      CASE "1"
         alicetemp
      CASE "2"
         alalbedo
      CASE "3"
         alground
   END SELECT
LOOP UNTIL al$ = "4"
END SUB

SUB alground
'shared ALAND(i)
DO
   CLS
   PRINT
   PRINT "     The albedos look like this from north to equator"
   PRINT
   FOR i = 1 TO 9
      PRINT "     ("; LTRIM$(RTRIM$(STR$(i))); ")        "; LATZ$(i); "   "; ALAND(i)
   NEXT i
   PRINT
   PRINT "     Which one do you want to alter ( zero for none of them )"
   PRINT
   INPUT "      Enter the number ", i
   IF i = 0 THEN EXIT DO
   PRINT
   PRINT "      The old value in band"; i; "is "; ALAND(i)
   DO
      INPUT "      What is your new value ? ", ALAND(i)
      IF ALAND(i) < 0 OR ALAND(i) > 1 THEN bad ELSE EXIT DO
   LOOP
LOOP
END SUB

SUB alicetemp
'shared TCRIT
CLS
PRINT
PRINT
PRINT "     The default value is -10 deg.C"
PRINT
INPUT "     What is the new value you want for TCRIT ? ", TCRIT
END SUB

SUB bad
PRINT "     Illegal response try again"
END SUB

SUB initialize
'shared A, B, C, AICE, TCRITE, ALAND(i), ALANDorig(i),
'shared TSTART(i), TSTARTorig(i)
A = 204
B = 2.17
C = 3.8
AICE = .62
TCRIT = -10
FOR i = 1 TO 9
   ALAND(i) = ALANDorig(i)
   TSTART(i) = TSTARTorig(i)
NEXT i
END SUB

SUB intro
CLS
PRINT
PRINT
PRINT
PRINT "********************************************************************************"
PRINT "                                ENERGY"
PRINT "                                BALANCE"
PRINT "                                MODEL"
PRINT "********************************************************************************"
PRINT
PRINT "          Copyright   A. Henderson-Sellers"
PRINT "                       & K. McKuffie"
PRINT "                                            (1986)"
PRINT
PRINT "     This model is similar to those of Budyko and Sellers."
PRINT "     You will be offered the opportunity to alter the"
PRINT "     values of the parameters which control the model climate."
PRINT
PRINT
PRINT
PRINT "               Press the space bar to continue"
waitforspace
CLS
PRINT
PRINT
PRINT
PRINT "           There are various possibilities for changing"
slow
PRINT "           the model climate.     You can then test the"
slow
PRINT "           sensitivity of this climate to changes in the"
slow
PRINT "           solar constatn.    That you should observe the"
slow
PRINT "           changes due to your changing of the model"
slow
PRINT "           parameters is also of importance in understanding"
slow
PRINT "           the nature of this model."
PRINT
PRINT
slow
PRINT "           Press space to continue"
waitforspace
END SUB

SUB lattrans
'shared C
CLS
PRINT "* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *"
PRINT "                       T R A N S P O R T"
PRINT "* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *"
PRINT
PRINT "     In this case you can alter the rate at which heat is"
PRINT "     transorted around the model by varying the value of C"
PRINT "     in the following equation."
PRINT
PRINT "            Heat Flux = C x ( T(mean) - T(zone) )"
PRINT
PRINT "                 The current value is "; C
PRINT
DO
   INPUT "  What is the value you want to use ? "; C
   IF C <= 0 OR C >= 50 THEN bad ELSE EXIT DO
LOOP

END SUB

SUB longwave
'shared A, B
CLS
PRINT "* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *"
PRINT "                 L O N G W A V E  L O S S  T O  S P A C E"
PRINT "* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *"
PRINT
PRINT "        The longwave loss to space is determined by the"
PRINT "        following equation."
PRINT
PRINT "            R = A + B x T(zone) "
PRINT
PRINT "      currently       A="; A
PRINT "                      B="; B
PRINT
PRINT "  Enter 1 to change them 0 to keep them the same"
DO
   r$ = INPUT$(1)
LOOP UNTIL r$ = "1" OR r$ = "0"
IF r$ = "1" THEN
   PRINT ""
   INPUT " Enter new value of A"; A
   INPUT " Enter new value of B"; B
END IF
END SUB

SUB mainmenu
DO
   CLS
   PRINT
   PRINT "                    * * * * * * * * * * * * * * * * * * * * "
   PRINT "                             M A I N   M E N U"
   PRINT "                    * * * * * * * * * * * * * * * * * * * * "
   PRINT
   PRINT
   PRINT "        There are 3 main parameterization schemes within the "
   PRINT "        model.  You may make alterations to any or all of them"
   PRINT "        at any one time.  Any which you choose not to alter"
   PRINT "        will be filled by default values."
   PRINT
   PRINT
   PRINT "  Choice                     Parameterization"
   PRINT "--------------------------------------------------------------------------------"
   PRINT "   (1)                        ALBEDO"
   PRINT "   (2)                        LATITUDINAL TRANSPORT"
   PRINT "   (3)                        LONGWAVE RADIATION TO SPACE"
   PRINT "   (4)                        RUN"
   PRINT "   (5)                        RESET VARIABLES"
   PRINT "   (6)                        QUIT"
   PRINT
   PRINT "  Enter the number of your choice"
   DO
      n$ = INPUT$(1)
   LOOP UNTIL VAL(n$) >= 1 AND VAL(n$) <= 6
   SELECT CASE n$
      CASE "1"
         albedomenu
      CASE "2"
         lattrans
      CASE "3"
         longwave
      CASE "4"
         runmodel
      CASE "5"
         initialize
      CASE "6"
         PRINT "good bye"
   END SELECT
LOOP UNTIL n$ = "6"
END SUB

SUB runmodel
'shared TSTART(i), TEMP(i), TSTARTorig(i)
CLS
PRINT
PRINT
PRINT "     What fraction of the solar constant would you like ?"
INPUT "     Your choice >", sx
'  start of routine to calculate temperatures
FOR lat = 1 TO 9
   TSTART(lat) = TSTARTorig(lat)
   TEMP(lat) = TSTART(lat)
NEXT lat
FOR H = 1 TO 50
   SOLCON = sx * 1370 / 4
   ' Calculate albedo of zones
   LATICE = 0
   FOR lat = 1 TO 9
      NL = 0
      ALBEDO(lat) = ALAND(lat)
      IF TEMP(lat) <= TCRIT THEN
         ALBEDO(lat) = AICE
         IF lat = 9 THEN
            NL = 0
         ELSE
            IF TEMP(lat + 1) > TCRIT THEN
               DP = (-TCRIT + TEMP(lat + 1)) * .1745 / (TEMP(lat + 1) - TEMP(lat))
               WORK2 = (PIBY2 - (lat + .5) * IN)
               LATICE = WORK2 + DP
               IF DP > .0872564 THEN
                  A3 = PIBY2 - lat * IN
                  A4 = PIBY2 - (lat - 1) * IN
                  A5 = (SIN(LATICE) - SIN(A3)) / (SIN(A4) - SIN(A3))
                  NC = AICE - (AICE - ALBEDO(lat)) * A5
                  NL = lat
               ELSE
                  A3 = PIBY2 - (lat + 1) * IN
                  A4 = A3 - IN
                  A5 = (SIN(A4) - SIN(LATICE)) / (SIN(A4) - SIN(A3))
                  NC = ALBEDO(lat + 1) * (1 - A5) + AICE * A5
                  NL = lat + 1
               END IF
            END IF
         END IF
      END IF
   NEXT lat
   IF ALBEDO(1) = ALAND(1) THEN NI = 90 / 57.296
   'calculate mean temperature
   SM = 0
   FOR lat = 1 TO 9
      WORK1 = PIBY2 - (lat - 1) * IN
      WORK2 = WORK1 - IN
      AC = ALBEDO(lat)
      IF lat = NL THEN AC = NC
      SM = SM + (SIN(WORK1) - SIN(WORK2)) * AC * S(lat)
   NEXT lat
   TX = (SOLCON * (1 - SM) - A) / B
   'end of mean temp routine
   FOR lat = 1 TO 9
      TM(lat) = TEMP(lat)
      TEMP(lat) = (SOLCON * S(lat) * (1 - ALBEDO(lat)) - A + C * TX)
      TEMP(lat) = FNR(TEMP(lat) / (C + B))
   NEXT lat
   ' now test for convergence
   AM = 0
   IC = 0
   FOR lat = 1 TO 9
      MA = ABS(TEMP(lat) - TM(lat))
      IF MA > AM THEN AM = MA
   NEXT lat
   IF AM < .01 THEN IC = 1
   IF IC = 1 THEN EXIT FOR
NEXT H
LATICE = FNR(LATICE * 57.296)
IF IC <> 1 THEN
   PRINT "  Model has failed to converge - - Think about your input "
ELSE
   'show the results
   CLS
   PRINT " ------------------------------------------------------------------------------ "
   PRINT "                                R E S U L T S"
   PRINT " ------------------------------------------------------------------------------ "
   PRINT "            Zone                 Temperature            Albedo"
   PRINT "         =============         ===========              ======"
   FOR lat = 1 TO 9
      PRINT "               "; LATZ$(lat); "            ";
      PRINT USING "###.#"; TEMP(lat);
      PRINT USING "                   ###.##"; ALBEDO(lat)
   NEXT lat
   PRINT "         The ice edge is at "; LATICE; " degrees north."
   PRINT
   PRINT "         Input parameters "
   PRINT "  Fraction of solar constant"; FNR(sx)
   PRINT "         A="; A; "      B="; B
   PRINT "         C="; C
   PRINT " Ice albedo= "; AICE; "Changes at "; TCRIT; " Deg C"
   PRINT
   PRINT "   Do you want to save these results(y/n)"
   DO
      ch$ = LCASE$(INPUT$(1))
   LOOP UNTIL ch$ = "y" OR ch$ = "n"
   IF ch$ = "y" THEN save
   PRINT ""
   PRINT "press any key to continue"
   A$ = INPUT$(1)
END IF

END SUB

SUB save
OPEN "results.txt" FOR OUTPUT AS #1
PRINT #1, " - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  "
PRINT #1, "                                R E S U L T S"
PRINT #1, " - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  "
PRINT #1, "            Zone                 Temperature            Albedo"
PRINT #1, "         =============         ===========              ======"
FOR lat = 1 TO 9
   PRINT #1, "               "; LATZ$(lat); "            ";
   PRINT #1, USING "###.#"; TEMP(lat);
   PRINT #1, USING "                   ###.##"; ALBEDO(lat)
NEXT lat
PRINT #1, "         The ice edge is at "; LATICE; " degrees north."
PRINT #1,
PRINT #1, "         Input parameters "
PRINT #1, "  Fraction of solar constant"; FNR(sx)
PRINT #1, "         A="; A; "      B="; B
PRINT #1, "         C="; C
PRINT #1, " Ice albedo= "; AICE; "Changes at "; TCRIT; " Deg C"
PRINT #1, ""
CLOSE #1
END SUB

SUB slow
FOR time = 1 TO 350000: NEXT time
PRINT
END SUB

SUB waitforspace
WHILE INPUT$(1) <> " ": WEND
PRINT
END SUB