cat EBM.f
      PROGRAM EBM
      IMPLICIT NONE

c	Daniel Wedul
c	CS 471
c	November 6, 2001
      
ccc	This program runs a 1 dimentional energy balance model.
c	Since I couldn't figure out how to correctly get fortran77 to
c	share variables between subprograms, I will use files to store
c	the values
c	the common variables are a, b, c, aice, tcrit
c	the common arrays are aland, tstart
c	the files are in folder SHARED and are named variablename.val

ccc	declare variables
      REAL	a,		!base longwave loss
     &		b,		!coefficient of lonwave loss
     &		c,		!coefficient of heat flux
     &		aice,		!albedo of ice
     &		tcrit		!temperature where ocean turns to ice
c	initialize variable files
      CALL initialize

c	display intro
      CALL intro

c	start the model
      CALL mainmenu

      END	!end program test



ccc	A sub to initialize the variables
      SUBROUTINE initialize
      REAL i			!temp variable
      REAL val(9)		!temp variable
      REAL val2(9)		!temp variable
      !a = 204
      i = 204
      OPEN(10, FILE="./SHARED/a.val", FORM='FORMATTED')
      WRITE(10,*) i
      CLOSE(10) 
      !b = 2.17
      i = 2.17
      OPEN(10, FILE="./SHARED/b.val", FORM='FORMATTED')
      WRITE(10,*) i
      CLOSE(10) 
      !c = 3.8
      i = 3.8
      OPEN(10, FILE="./SHARED/c.val", FORM='FORMATTED')
      WRITE(10,*) i
      CLOSE(10)
      !aice = .62
      i = .62
      OPEN(10, FILE="./SHARED/aice.val", FORM='FORMATTED')
      WRITE(10,*) i
      CLOSE(10) 
      !tcrit = -10
      i = -10
      OPEN(10, FILE="./SHARED/tcrit.val", FORM='FORMATTED')
      WRITE(10,*) i
      CLOSE(10) 
      !aland(i) = (.5, .5, .452, .407, .357, .309, .272, .248, .254)
      DATA val /.5, .5, .452, .407, .357, .309, .272, .248, .254/
      OPEN(10, FILE="./SHARED/aland.val", FORM='FORMATTED')
      DO i=1,9
      	WRITE(10,*) val(i)
      ENDDO
      CLOSE(10) 
      END	!end initialize

ccc     sub to display the intro
      SUBROUTINE intro
      INTEGER i
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*) ("*", i=1,30)
      WRITE(6,*) "	ENERGY"
      WRITE(6,*) "	BALANCE"
      WRITE(6,*) "	MODEL"
      WRITE(6,*) ("*", i=1,30)
      WRITE(6,*)
      WRITE(6,*) "   Copyright A. Henderson-Sellers"
      WRITE(6,*) "	     & K. McKuffie"
      WRITE(6,*) "			(1986)"
      WRITE(6,*)
      WRITE(6,*) " This model is similar to those of Budyko ",
     &		  "and Sellers."
      WRITE(6,*) " You will be offered the opportunity to alter the"
      WRITE(6,*) " values of the parameters which control the ",
     &		  "model climate."
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*) "   Press return to continue"
      READ(5,*)
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*) " There are arious possibilities for changing"
      WRITE(6,*)
      WRITE(6,*) " the model climate.  You can then test the"
      WRITE(6,*)
      WRITE(6,*) " sensitivity of this climate to changes in the"
      WRITE(6,*)
      WRITE(6,*) " solar constant.  That you should observe the"
      WRITE(6,*)
      WRITE(6,*) " changes due to your changing of the model"
      WRITE(6,*)
      WRITE(6,*) " parameters is also of importance in understanding"
      WRITE(6,*)
      WRITE(6,*) " the nature of this model."
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*) " Press return to continue"
      READ(5,*)
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*)
      END 	!end intro



ccc	The sub to govern main flow of the program
      SUBROUTINE mainmenu

      CHARACTER*10 choice
      INTEGER i
      choice = "0"
      DOWHILE(choice .ne. "6")
      	WRITE(6,*)
      	WRITE(6,*)
      	WRITE(6,*)
      	WRITE(6,*) ("* ", i=1,12)
      	WRITE(6,*) "   M A I N   M E N U"
      	WRITE(6,*) ("* ", i=1,12)
      	WRITE(6,*)
      	WRITE(6,*) " There are 3 main parameterization schemes within the"
      	WRITE(6,*) " model. You may make alterations to any or all of them"
      	WRITE(6,*) " at any one time. Any which you chose not to alter"
      	WRITE(6,*) " will be filled by default values."
      	WRITE(6,*)
      	WRITE(6,*)
      	WRITE(6,*) " Choice	Parameterization"
      	WRITE(6,*) ("-", i=1,40)
      	WRITE(6,*) " (1)		ALBEDO"
      	WRITE(6,*) " (2)		LATITUDINAL TRANSPORT"
      	WRITE(6,*) " (3)		LONGWAVE RADIATION TO SPACE"
      	WRITE(6,*) " (4)		RUN"
      	WRITE(6,*) " (5)		RESET VARIABLES TO DEFAULT"
      	WRITE(6,*) " (6)		QUIT"
      	WRITE(6,*)
      	WRITE(6,*) " Enter the number of your choice"
	READ(5,*) choice
      	WRITE(6,*)
      	WRITE(6,*)
      	WRITE(6,*)
	IF (choice .eq. "1") THEN 
	  CALL albedomenu
	ENDIF
	IF (choice .eq. "2") THEN 
	  CALL lattrans
	ENDIF
	IF (choice .eq. "3") THEN 
	  CALL longwave
	ENDIF
	IF (choice .eq. "4") THEN 
	  CALL runmodel
	ENDIF        
	IF (choice .eq. "5") THEN 
	  CALL initialize
	ENDIF        
	IF (choice .eq. "6") THEN 
	  WRITE(6,*) "GOOD BYE"
	ENDIF
      ENDDO
      END	!end mainmenu

ccc	sub to change albedo options
      SUBROUTINE albedomenu
      INTEGER i
      CHARACTER*10 choice
      choice = "0"
      DOWHILE(choice .ne. "4")
      	WRITE(6,*)
      	WRITE(6,*)
      	WRITE(6,*)
      	WRITE(6,*)
      	WRITE(6,*) ("* ", i=1,22)
      	WRITE(6,*) "   P A R A M E T E R I Z A T I O N   O F "
      	WRITE(6,*) "             A L B E D O"
      	WRITE(6,*) ("* ", i=1,22)
      	WRITE(6,*)
      	WRITE(6,*)
      	WRITE(6,*) "  There are three things which you may alter"
      	WRITE(6,*)
      	WRITE(6,*) "	1. The temperature at which the surface"
      	WRITE(6,*) "	   becomes ice covered."
      	WRITE(6,*)
      	WRITE(6,*) "	2. The albedo of this ice covered surface"
      	WRITE(6,*)
      	WRITE(6,*) "	3. the albedo of the underlying ground"
      	WRITE(6,*)
      	WRITE(6,*) "	4. Return to the main menu"
      	WRITE(6,*)
      	WRITE(6,*) " Choose which one you want"
	READ(5,*) choice
	IF (choice .eq. "1") THEN
	  CALL albedoone
	ENDIF
	IF (choice .eq. "2") THEN
	  CALL albedotwo
	ENDIF
	IF (choice .eq. "3") THEN
	  CALL albedothree
	ENDIF
      ENDDO
      END	!end albedomenu

ccc	sub to change temp where surface becomes ice
      SUBROUTINE albedoone
      REAL tcrit
      ! get current value of tcrit
      OPEN(10, FILE="./SHARED/tcrit.val", FORM='FORMATTED')
      READ(10,*) tcrit
      CLOSE(10)

      WRITE(6,*) " The current value of tcrit is ", tcrit
      WRITE(6,*)
      WRITE(6,*) " What is the new value you want for tcrit?"
      READ(5,*) tcrit

      ! save new value of tcrit
      OPEN(10, FILE="./SHARED/tcrit.val", FORM='FORMATTED')
      WRITE(10,*) tcrit
      CLOSE(10)
      END	!end albedoone

ccc	sub to change albedo of ice
      SUBROUTINE albedotwo
      REAL aice
      ! get current value of aice
      OPEN(10, FILE="./SHARED/aice.val", FORM='FORMATTED')
      READ(10,*) aice
      CLOSE(10)

      WRITE(6,*) " The current albedo of the ice is ", aice
      WRITE(6,*)
      WRITE(6,*) " What is the new value you want for this albedo?"
      READ(5,*) aice

      ! save new value of aice
      OPEN(10, FILE="./SHARED/aice.val", FORM='FORMATTED')
      WRITE(10,*) aice
      CLOSE(10)
      END	!end albedotwo

ccc	sub to change albedo of ground
      SUBROUTINE albedothree
      REAL aland(9)
      INTEGER i
      !get aland
      OPEN(10, FILE="./SHARED/aland.val", FORM='FORMATTED')
      DO i=1,9
	READ(10,*) aland(i)
      ENDDO

      DOWHILE(i .ne. 0)
	WRITE(6,*)
	WRITE(6,*)
	WRITE(6,*)
	WRITE(6,*) " The albedos look like this from north to equator"
	WRITE(6,*)
        DO i=1,9
	  WRITE(6,*) "   ",i,"   ",(9-i)*10,"-",(10-i)*10,"    ",aland(i)
	ENDDO
	WRITE(6,*) " Which one do you want to alter ",
     &  	 	     "(zero for none of them)"        
	WRITE(6,*)
	READ(5,*), i
	IF (i .GT. 0 .AND. i .LE. 9) THEN
	  WRITE(6,*) "The old value in band ",i," is ",aland(i)
	  WRITE(6,*) "The new value needs to be between 0 and 1."
	  WRITE(6,*) "What is the new value?"
	  !get input
	  aland(i)=-1
	  DOWHILE(aland(i) <0 .OR. aland(i) >1)
	    READ(5,*) aland(i)
	  ENDDO
	  OPEN(10, FILE="./SHARED/aland.val", FORM='FORMATTED')
 	  DO i=1,9
	    WRITE(10,*), aland(i)
	  ENDDO
	ENDIF	  
      ENDDO
      END	!end albedothree

ccc	sub to change latitudinal transport
      SUBROUTINE lattrans
      INTEGER i
      ! read the parameter
      REAL c
      OPEN(10, FILE="./SHARED/c.val", FORM='FORMATTED')
      READ(10,*) c
      CLOSE(10)

      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*) ("* ", i=1,20)
      WRITE(6,*) "    T R A N S P O R T "
      WRITE(6,*) ("* ", i=1,20)
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*) "  In this case you can alter the rate at which ",
     &		 "heat is"
      WRITE(6,*) "  transported around the model by varying the ",
     &		 "value of C"
      WRITE(6,*) "  in the following eqation."
      WRITE(6,*)
      WRITE(6,*) "	Heat Flux = C x (T(mean)-T(zone))"
      WRITE(6,*)
      WRITE(6,*) "	The current C is ",c
      WRITE(6,*)
      WRITE(6,*) "Enter a new value of c from 0 to 50"
      c = -10
      DOWHILE(c .LE. 0 .OR. c .GE. 50)
	READ(5,*) c
      ENDDO
      !save new c value
      OPEN(10, FILE="./SHARED/c.val", FORM='FORMATTED')
      WRITE(10,*) c
      CLOSE(10)
      END	!end lattrans

ccc	sub to chang longwave out
      SUBROUTINE longwave
      REAL a, b
      INTEGER i
      !get a and b
      OPEN(10, FILE="./SHARED/a.val", FORM='FORMATTED')
      READ(10,*) a
      CLOSE(10)
      OPEN(10, FILE="./SHARED/b.val", FORM='FORMATTED')
      READ(10,*) b
      CLOSE(10)	
      !change the values
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*) ("* ", i=1,25)
      WRITE(6,*) "    L O N G W A V E   L O S S   T O   S P A C E"
      WRITE(6,*) ("* ", i=1,25)
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*) "  The longwave loss to space is determined by"
      WRITE(6,*) "  the following equation."
      WRITE(6,*)
      WRITE(6,*) "	R = A + B x T(zone)"
      WRITE(6,*)
      WRITE(6,*) "  currently A = ", a
      WRITE(6,*) "	      B = ", b    
      WRITE(6,*)
      WRITE(6,*) "enter a 1 to change these and 0 to keep them"
      READ(5,*) i
      IF (i .eq. 1) THEN
	WRITE(6,*) "Enter new value of A ";
	READ(5,*) a
	WRITE(6,*) "Enter new value of b ";
	READ(5,*) b
      ENDIF
      !save values
      OPEN(10, FILE="./SHARED/a.val", FORM='FORMATTED')
      WRITE(10,*) a
      CLOSE(10)
      OPEN(10, FILE="./SHARED/b.val", FORM='FORMATTED')
      WRITE(10,*) b
      CLOSE(10)	
      END	!end longwavye

ccc	sub to run the model
      SUBROUTINE runmodel
      REAL 	sx,
     &		a,
     &		b,
     &		c,
     &		aice,
     &		tcrit,
     &		IN,
     &		PIBY2,
     & 		solcon,
     &		latice,
     &		dp,
     &		nc,
     &		sm,
     &		tx,
     &		ma,
     &		am,
     &		ac,
     &		a3,a4,a5
     &		work1,
     &		work2,
     &		x

      INTEGER	lat,
     &		nl,
     &		h,
     &		ic

      REAL	aland(9),
     &		tstart(9),
     &		temp(9),
     &		albedo(9),
     &		tm(9),
     &		s(9)

      PARAMETER(IN = 3.14159/18, PIBY2=3.14159/2)

      !load shared variables
      OPEN(10, FILE="./SHARED/a.val", FORM='FORMATTED')
      READ(10,*) a
      CLOSE(10)
      OPEN(10, FILE="./SHARED/b.val", FORM='FORMATTED')
      READ(10,*) b
      CLOSE(10)
      OPEN(10, FILE="./SHARED/c.val", FORM='FORMATTED')
      READ(10,*) c
      CLOSE(10)
      OPEN(10, FILE="./SHARED/aice.val", FORM='FORMATTED')
      READ(10,*) aice
      CLOSE(10)
      OPEN(10, FILE="./SHARED/tcrit.val", FORM='FORMATTED')
      READ(10,*) tcrit
      CLOSE(10)
      OPEN(10, FILE="./SHARED/aland.val", FORM='FORMATTED')
      DO i=1,9
	READ(10,*) aland(i)
      ENDDO
      CLOSE(10)

      ! do any other initializations
      DATA tstart /-16.9, -12.3, -5.1, 2.2, 8.8, 
     &			16.2, 22.9, 26.1, 26.4 /
      DATA s / .5, .531, .624, .77, .892,
     &		1.021, 1.12, 1.189, 1.219/

      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*)
      WRITE(6,*) "What fraction of the solar constant would you like?"
      READ(5,*)	sx

      DO lat=1,9
	temp(lat)=tstart(lat)
      ENDDO

      ! start of routine to calculate temperatures
      h = 0	
      DOWHILE(h .LE. 50)
	h = h + 1
	solcon = sx*1370/4
	! calculate albedo of zones
	latice = 0
	DO lat=1,9
	  nl = 0
	  albedo(lat) = aland(lat)
	  IF( temp(lat) .LE. tcrit) THEN
	    albedo(lat) = aice
	    IF (lat .EQ. 9) THEN
	      nl = 0
	    ELSE
	      IF ( temp(lat+1) .GT. tcrit) THEN
		dp = (-tcrit + temp(lat+1))*.1745/(temp(lat+1)-temp(lat))
		work2 = (PIBY2-(lat+.5)*IN)
		latice = work2 + dp
		IF (dp .GT. .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
		ENDIF
	      ENDIF
	    ENDIF
	  ENDIF
	ENDDO	!end lat do loop
	IF (albedo(1) .EQ. aland(1)) THEN
	  NI=90/57.296
	ENDIF
	!calculate mean temperature
	sm = 0
	DO lat=1,9
	  work1 = PIBY2-(lat-1)*IN
	  work2 = work1 - IN
	  ac = albedo(lat)
	  IF (lat .EQ. NL) THEN
	    ac = nc
	  ENDIF
	  sm = sm+(SIN(work1)-SIN(work2))*ac*s(lat)
	ENDDO
	tx = (solcon * (1 - sm) - a) / b
	! end of mean temperature routine
	
	DO lat = 1,9
	  tm(lat) = temp(lat)
	  temp(lat) = (solcon*s(lat)*(1-albedo(lat))-a+c*tx)
	  temp(lat) = temp(lat)/(c+b)
	ENDDO

	am = 0
	ic = 0
	DO lat = 1,9
	  ma = ABS(temp(lat)-tm(lat))
	  IF (ma .GT. am) THEN
	    am = ma
	  ENDIF
	ENDDO

	IF (am .LT .01) THEN
	  ic = 1
	ENDIF
	IF (ic .EQ. 1) THEN
	  h = 51
	ENDIF
      ENDDO	!end of h loop 
      IF (ic .NE. 1) THEN
	WRITE(6,*) " Model has failed to converge ",
     &		      "-- Think about your input"
      ELSE
	! show the results
	WRITE(6,*)
	WRITE(6,*)
	WRITE(6,*)
	WRITE(6,*) ("-",h=1,60)
	WRITE(6,*) "		RESULTS"
	WRITE(6,*) ("-", h=1,60)
	WRITE(6,*)
	WRITE(6,*) "	Zone	Temperature	Albedo"
	WRITE(6,*) "    ====    ===========     ======"
	DO lat=1,9
	  WRITE(6,1000) (9-lat)*10, (10-lat)*10, temp(lat), albedo(lat)
	ENDDO
	latice = latice*57.296
	WRITE(6,1005) latice
	WRITE(6,*)
	WRITE(6,*) "	input parameters"
	WRITE(6,1010) sx
  	WRITE(6,1020) a, b
	WRITE(6,1030) c
	WRITE(6,1040) aice, tcrit
	WRITE(6,*)
	WRITE(6,*) "Press return to continue"
      ENDIF	
      READ(5,*)

1000  FORMAT("    ",I2.2,"-",I2.2,"      ",F5.1,"          ",F3.2)
1005  FORMAT(" The ice edge is at ", F5.2," Degrees North")
1010  FORMAT(" fraction of solar constant = ", F4.2)
1020  FORMAT("      A = ", F6.1, "   B = ",F5.2)
1030  FORMAT("      C = ", F5.2)
1040  FORMAT(" Ice albed = ", F4.2," changes at ", F5.1," Deg C")
      END	!end runmodel



[sck01@kluane MIDTERM]$ exit
logout