C This program calculates the sensitivity of the climate system to competing feedback C processes by allowing for various distributions of land and sea C REAL TL(1:9), CIL, CIW(1:9), S(1:9), ALPHAL(1:9), ALPHAW(1:9) REAL K, KM, TW(1:9),LANDFRAC, WATERFRAC,ALPHAM(1:9) REAL TIMESTEP,DTL(1:9),DTW(1:9),WVFL,WVFW,PI,SF,TBAR ,T(1:9) REAL ALPHALI(1:9),ALPHALC(1:9), ALPHALB(1:9),ALPHAWI(1:9) REAL ALPHAWC(1:9), A, B, IF,TBARLANDSUM REAL TBARWATERSUM,TBARLAND,TBARWATER REAL INCREMENT,ECC,TILT,DN,DECANGLE,SRANGLE,SCON,X,ALBLAT REAL ALPHAMON(1:9),ECCRANGE,TAVEYR,TAVEYRSUM,ALBAVE,ALBAVESUM REAL DELTACO2,KORIG,KMORIG,KMULT,KMMULT LOGICAL VOLCANO,SEAICE(1:9) CHARACTER * 3 PERIHELION,peri INTEGER LAT PARAMETER (SCON = 1370.0) c PARAMETER (VOLCANO = .FALSE.) PARAMETER (VOLCANO = .TRUE.) PARAMETER (PI = 3.14159) PARAMETER (WVFL = 0.0067) PARAMETER (WVFW = 0.0134) PARAMETER (SF = 0.01) PARAMETER (IF = 0.02) PARAMETER (CIL = 4.186E7) PARAMETER (B = 2.17) PARAMETER (KORIG = 3.80) PARAMETER (KMORIG = 1.90) PARAMETER (TIMESTEP= 365.25/12.*86400.) write (*,*) write (*,*) 'This is an Energy Balance Climate Model designed to' write (*,*) ' test the role of changes in climatic variables.' write (*,*) ' You will be asked to choose values for various ' write (*,*) 'parameters, such as the tilt of the earth,' write (*,*) ' CO2 levels,etc.' write (*,*) write (*,*) 'Choose the date of perihelion (the date when the' write (*,*) 'earth and sun are closest together). If you choose' write (*,*) 'January, enter "JAN". If you choose July, ' write (*,*) 'enter "JUL". The modern condition is "JAN" ' read (*,*) peri if (peri .eq. 'JAN') perihelion = 'JAN' if (peri .eq. 'JUL') perihelion = 'JUL' write (*,*) write (*,*) 'Choose the tilt of the earth. Enter an angle' write (*,*) 'between 22.45 and 24.45. The modern value is 23.45' read (*,*) tilt write (*,*) write (*,*) 'Choose the eccentricity for the orbit of the earth.' write (*,*) 'Enter a number between 0 and 0.06. The modern value' write (*,*) 'is 0.016' read (*,*) eccrange write (*,*) write (*,*) 'Choose the carbon dioxide level of the atmosphere.' write (*,*) 'Enter a CO2 level as a multiple of the modern value' write (*,*) 'For example, enter "1" for the modern level,"2" for' write (*,*) 'an atmosphere with twice as much CO2 as today' read (*,*) deltaco2 A = 204.0 + -6.*alog(deltaco2) write (*,*) write (*,*) 'Choose the fraction of global land cover. Enter' write (*,*) 'a value between 0 (all-water planet) and 1 (all-' write (*,*) 'land planet). The modern value is about 0.33' read (*,*) landfrac waterfrac = 1 - landfrac write (*,*) write (*,*) 'Choose the magnitude of meridional (north-south)' write (*,*) 'atmospheric heat transport efficiency. Enter a' write (*,*) 'multiple of the modern efficiency, ranging from' write (*,*) '0.5 to 4 times the modern value.' read (*,*) kmult k = kmult * korig write (*,*) write (*,*) 'Choose the magnitude of zonal (east-west)' write (*,*) 'atmospheric heat transport efficiency. Enter a' write (*,*) 'multiple of the modern efficiency, ranging from' write (*,*) '0.5 to 4 times the modern value.' read (*,*) kmmult km = kmmult * kmorig C OPEN (UNIT = 1, FILE = 'mixed.out', FORM = 'FORMATTED', $ ACCESS = 'SEQUENTIAL', STATUS= 'old' ) C C INITIALIZE TEMPERATURES AND ALBEDOS AT EACH LATITUDE C M = 1 DO 5 I = 1, 9 TL(I) = 288.0 TW(I) = 288.0 T(I) = 288.0 TBAR = 288.0 TBARLAND = 288.0 TBARWATER = 288.0 TAVEYRSUM = 0.0 ALPHALI(I) = 0.0 ALPHALC(I) = 0.0 ALPHAWI(I) = 0.0 ALPHAWC(I) = 0.0 ALPHAL(I) = 0.18 + 0.04*I ALPHALB(I) = ALPHAL(I) CIW(I) = 4.186E8 SEAICE(I) = .FALSE. 5 CONTINUE C C ALLOW FOR CHANGES IN TEMPERATURE AND ALBEDO DUE TO DIFFERENTIAL SOLAR FORCING C DO 6 L = 1,20 DO 7 J = 1,12 ALBAVESUM = 0.0 TBARSUM = 0.0 TBARVSUM = 0.0 TBARLANDSUM = 0.0 TBARWATERSUM = 0.0 IF (L.EQ.20) then WRITE (1,*) 'YEAR = ', L, ' MONTH = ', J WRITE (1,*) 'LATITUDE',' LAND ',' WATER ', ' AVERAGE ', $ ' LAND ',' WATER ', ' AVERAGE ', 'SOLAR ',' DEC ',' ECC ' WRITE (1,*) ' ',' TEMP ','TEMP ',' TEMP ', $ ' ALBEDO ','ALBEDO ','ALBEDO ',' RAD ' ENDIF C******** Tom *** day #, eccentricity, perhihelion, tilt****** DN = 365.25/12.*(J-1) + 15./30.*365.25/12. DECANGLE = TILT*SIN((DN+284)*PI/180.*360./365.25) IF (PERIHELION .EQ. 'JAN') THEN ECC = 1. + ECCRANGE*COS(2.*PI*DN/365.25) ELSE ECC = 1. + ECCRANGE*COS(2.*PI*(DN+365.25/2.)/365.25) ENDIF C************ ******** DO 10 I = 1,9 LAT = I*10 - 5 ALBLAT = (LAT-DECANGLE)/1.2*PI/180. DENCOS = COS(ALBLAT) IF (DENCOS .LT. 0.15) DENCOS = 0.15 ALPHAMON(I) = 0.06/DENCOS ALPHAW(I) = I*0.04+ALPHAMON(I)+ALPHAWI(I)+ALPHAWC(I) ALPHAL(I) = ALPHALB(I) + ALPHALI(I) + ALPHALC(I) C IF ((L.EQ.14).AND.(VOLCANO .EQ. .TRUE.)) THEN IF (((L.EQ.14).AND.(J.GE.10)).OR.((L.EQ.15) $ .AND.(J.LE.9)).AND.(VOLCANO.EQ. .TRUE.)) THEN INCREMENT = 0.20 - (M -1)*.018 C INCREMENT = 0.20 - (J-1)*.018 ALPHAL(I) = ALPHAL(I) + INCREMENT ALPHAW(I) = ALPHAW(I) + INCREMENT IF (ALPHAL(I).GT.(0.70+INCREMENT)) ALPHAL(I) = 0.70 + $ INCREMENT IF (ALPHAW(I).GT.(0.60+INCREMENT)) ALPHAW(I) = 0.60 + $ INCREMENT ELSE IF (ALPHAL(I).GT. 0.70) ALPHAL(I) = 0.70 IF (ALPHAW(I).GT. 0.60) ALPHAW(I) = 0.60 ENDIF ALPHAM(I) = ALPHAL(I)*LANDFRAC + ALPHAW(I)*WATERFRAC X = -TAN(LAT*PI/180.)*TAN(DECANGLE*PI/180.) IF (ABS(X) .LE. 1.0) THEN SRANGLE = ACOS(X) ELSE IF (X .LT. -1.0) THEN SRANGLE = PI ELSE SRANGLE = 0.0 ENDIF S(I) = (1./PI)*SCON*ECC*SIN(LAT*PI/180.)* $ SIN(DECANGLE*PI/180.)*(SRANGLE-TAN(SRANGLE)) DTL(I) = TIMESTEP/CIL*(S(I)*(1-ALPHAL(I))- $ (A+B*(TL(I)-273.2))-K*(T(I)-TBAR)-KM*(TL(I)-T(I))) DTW(I) = TIMESTEP/CIW(I)*(S(I)*(1-ALPHAW(I))- $ (A+B*(TW(I)-273.2))-K*(T(I)-TBAR)-KM*(TW(I)-T(I))) TL(I) = TL(I) + DTL(I) TW(I) = TW(I) + DTW(I) IF ((TW(I).GT.273.2).AND.(SEAICE(I) .EQ. .TRUE.)) THEN TW(I) = 271.2 + (TW(I) - 273.2)*0.1 SEAICE(I) = .FALSE. CIW(I) = 4.186E8 ELSE IF ((TW(I).LT. 271.2).AND.(SEAICE(I).EQ. .FALSE.)) $ THEN TW(I) = 271.2 + (TW(I) - 271.2)*10. SEAICE(I) = .TRUE. CIW(I) = CIL ENDIF T(I) = TL(I)*LANDFRAC + TW(I)*WATERFRAC IF (L .EQ. 20) THEN WRITE (1,100) LAT, TL(I)-273.15,TW(I)-273.15,T(I)-273.15, $ ALPHAL(I), ALPHAW(I), ALPHAM(I),S(I),DECANGLE, ECC C WRITE(1, 101) DTL(I),DTW(I), TBARLAND, TBARWATER, TBAR ENDIF TBARSUM = TBARSUM + T(I)*(SIN(I*10*PI/180.) - $ SIN((I-1)*10*PI/180.)) TBARLANDSUM = TBARLANDSUM + TL(I)*(SIN(I*10*PI/180.)- $ SIN((I-1)*10*PI/180.)) TBARWATERSUM = TBARWATERSUM + TW(I)*(SIN(I*10*PI/180.)- $ SIN((I-1)*10*PI/180.)) ALBAVESUM = ALBAVESUM + ALPHAM(I)*(SIN(I*10*PI/180.)- $ SIN((I-1)*10*PI/180.)) IF (((L.GE.14).AND.(J.GE.10)).OR.(L.GE.15).AND. $ (VOLCANO.EQ. .TRUE.)) $ TBARVSUM = TBARVSUM + T(I)*(SIN(I*10*PI/180.) - $ SIN((I-1)*10*PI/180.)) IF (TL(I) .LT. 273.2) THEN ALPHALI(I) = (273.2-TL(I))*SF ELSE IF (TL(I) .GT. 293.0) THEN ALPHALC(I) = (TL(I)-293.)*WVFL ELSE ALPHALI(I) = 0.0 ALPHALC(I) = 0.0 ENDIF IF (SEAICE(I) .EQ. .TRUE.) THEN IF (TW(I) .LT. 271.2) THEN ALPHAWI(I) = (271.2-TW(I))*IF ELSE ALPHAWI(I) = (273.2-TW(I))*IF*0.5 ENDIF ELSE IF (TW(I).GT. 293.0) THEN ALPHAWC(I) = (TW(I) -293.)*WVFW ELSE ALPHAWI(I) = 0.0 ALPHAWC(I) = 0.0 ENDIF 10 CONTINUE TBAR = TBARSUM ALBAVE = ALBAVESUM TBARLAND = TBARLANDSUM TBARWATER = TBARWATERSUM IF (L .EQ. 20) THEN WRITE (1,*) 'GLOBAL AVERAGE TEMPERATURE = ', TBAR-273.15, $ ' AVERAGE ALBEDO = ', ALBAVE write (1,*) ENDIF TAVEYRSUM = TAVEYRSUM + TBAR ALBAVEYRSUM = ALBAVEYRSUM + ALBAVE IF (((L.GE.14).AND.(J.GE.10)).OR.(L.GE.15)) M = M + 1 TBARV = TBARVSUM TVAVEYRSUM = TVAVEYRSUM + TBARV if (volcano .eq. .true.) then IF (L.GE.14) WRITE (2,*) L, J, ' TBARV = ', TBARV IF (M .EQ. 13) THEN TVAVEYR = TVAVEYRSUM/12. TVAVEYRSUM = 0.0 TBARVSUM = 0.0 WRITE (2,*) 'YEAR = ', L, ' RUNNING AVERAGE TEMP = ', $ TVAVEYR M = 1 ENDIF endif 7 CONTINUE ALBAVEYR = ALBAVEYRSUM/12. TAVEYR = TAVEYRSUM/12. TAVEYRSUM = 0.0 ALBAVEYRSUM = 0.0 IF (L .EQ. 20) THEN write (*,*) write (*,*) 'GLOBAL ANNUAL AVERAGE TEMPERATURE=',TAVEYR- $ 273.15 write (*,*) 'GLOBAL ANNUAL AVERAGE ALBEDO=',ALBAVEYR write (*,*) WRITE (*,*) 'If you want to see the temperatures for' write (*,*) 'each latitude and for each month, then' write (*,*) 'type "vi mixed.out". For a printout of' write (*,*) 'this information, type "lpe mixed.out"' write (*,*) write (1,*) 'GLOBAL ANNUAL TEMP = ', TAVEYR-273.15, $ ' GLOBAL ANNUAL ALBEDO = ', ALBAVEYR ENDIF 6 CONTINUE 100 FORMAT (I5, 3X, 2(F5.1,1X),F7.1,3x,3(F5.2,2X), $ F6.1,2X,F6.2,2X,F6.4) 101 FORMAT (2(F6.2,2X),3(F6.2,2X)) 110 STOP END