SUBROUTINE ORBIT INCLUDE 'PARAMETERS.INC' INCLUDE 'PLANET.INC' CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC O R B I T CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC CC CC PURPOSE: TO COMPUTE THE DISTANCE BETWEEN THE PLANET AND THE CC CC SUN, THE SOLAR DECLINATION ANGLE. CC CC CC CC INPUT: PERORB = PERIOD OF THE PLANET'S ORBIT (SEC) CC CC SMA = LENGTH OF SEMI-MAJOR AXIS OF THE ORBIT (M) CC CC ECC = ECCENTRICITY OF ORBIT CC CC DTPEQ = TIME BETWEEN THE PERIHELION CROSSING AND EQUI- CC CC NOX (SECONDS) CC CC = CC CC TE = TIME SINCE UERNAL EQUINOX CROSSING (SECONDS) CC CC PAXINC = TILT OF PLANETARY AXIS TO PLANE OF ORBIT (DEGREE) CC CC CC CC OUTPUT: DIST = DISTANCE TO SUN (METERS) CC CC DEC = SOLAR DECLINATION ANGLE (DEGREES) CC CC CC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC COMMON & /ONED3/ XAU & /ONED4/ DEC & /ONED36/ DAY DATA PI/3.1415926d0/ C Find time since perapsis crossing TIM = DAY + PERIH C Find eccentric anomoly using Newton's method FM = 2.0d0*PI*TIM/PYEAR E0 = ECC*SIN(FM) + FM DE = 1.0d0 DO WHILE( DE .GE. 1.0E-7) E = E0 + (FM - E0 + ECC*SIN(E0))/(1.0d0 - ECC*COS(E0)) DE = ABS(E - E0) E0 = E ENDDO C Find the radius and obrit angle between the planet and the sun DIST = AU*(1.0d0 - ECC*COS(E)) XAU = DIST/AU PSI = 2.0d0*ATAN(SQRT((1.0d0 + ECC)/(1.0d0 - ECC))*TAN(0.5d0*E)) C Find the angle between the vernal equinox and the perihelion crossing FM = 2.0d0*PI*PERIH/PYEAR E0 = ECC*SIN(FM) + FM DE = 1.0d0 DO WHILE( DE .GE. 1.0E-7) E = E0 + (FM - E0 + ECC*SIN(E0))/(1.0d0 - ECC*COS(E0)) DE = ABS(E - E0) E0 = E ENDDO ANG = 2.0d0*ATAN(SQRT((1.0d0 + ECC)/(1.0d0 - ECC))*TAN(0.5d0*E)) C Find the apparent declination angle in degrees DEC = PAXINC*SIN(PSI - ANG) RETURN END C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ SUBROUTINE DAYTIME INCLUDE 'PARAMETERS.INC' INCLUDE 'PLANET.INC' COMMON & /SWITCHES/ ZLAT , DELTAT , DELT1 , DAYSTART , & TIMESTART , TIMESTOP , OUTTIME , & MODE , NTIME , NITER , IDIFF , istepdump & /ONED4/ DEC & /ONED5/ SRTOP , SSTOP & /ONED21/ SRBOT , SSBOT & /BLOCK10/ Z(NALT) & /BLOCK27/ PI & /ONED32/ SUNRISE(NALT) , SUNSET(NALT) & /ONED34/ ANG1 , ANG2 C THIS SECTION COMPUTES THE LENGTH OF A DAY. IF ISPHERE = 1 , THEN THE DAY C IS A FUNCTION OF THE ALTITUDE DUE TO THE CURVATURE OF THE EARTH. PI = 4.0D0*DATAN(1.0D0) XDEC = DEC*(PI/180.0D0) XLAT = ZLAT*(PI/180.0D0) ANG1 = SIN(XDEC)*SIN(XLAT) ANG2 = COS(XDEC)*COS(XLAT) C THIS PORTION COMPUTES THE LENGTH OF DAY IF( ISPHERE .EQ. 0 ) THEN ! PLANE-PARALLEL ATMOSPHERE C GROUND LEVEL VALUES FOR LENGTH OF DAY (TDAY), ETC. IF (ABS(ANG1/ANG2) .LE. 1.0d0) THEN TDAY = ACOS(-ANG1/ANG2)/(PI/PDAY) ELSE IF (ANG1/ANG2 .GT. 1.0d0) TDAY = PDAY IF (ANG1/ANG2 .LT. -1.0d0) TDAY = 0.0d0 ENDIF SRBOT = 0.5d0*(PDAY - TDAY) SSBOT = SRBOT + TDAY SRBOT = JIDINT(SRBOT*60.0d0)/60.0d0 ! SUNRISE TO NEAREST MINUTE SSBOT = JIDINT(SSBOT*60.0d0)/60.0d0 SRBOT = DMAX1(SRBOT,0.0D0) IF ((SSBOT .LT. 0.0d0) .OR. (SSBOT .GT. PDAY)) SSBOT = PDAY SRTOP = SRBOT ! SUNRISE AT TOP OF MODEL SSTOP = SSBOT ! SUNSET AT TOP OF MODEL ELSE C CURVED ATMOSPHERE - LENGTH OF DAY VARIES WITH ALTITUDE. DO K = 1 , NALT CHI = ACOS(RADIUS/(RADIUS + Z(K))) + 0.5D0*PI COSSZA = (COS(CHI) - ANG1)/ANG2 IF (ABS(COSSZA) .LT. 1.0d0) THEN TDAY = ACOS((COS(CHI) - ANG1)/ANG2)/(PI/PDAY) ELSE IF (COSSZA .GE. 1.0d0) TDAY = 0.0d0 IF (COSSZA .LE. -1.0d0) TDAY = PDAY ENDIF SUNRISE(K) = 0.5d0*(PDAY - TDAY) SUNSET(K) = SUNRISE(K) + TDAY SUNRISE(K) = DMAX1(0.0D0,SUNRISE(K)) IF ((SUNSET(K) .LT. 0.0d0) .OR. (SUNSET(K) .GT. PDAY)) THEN SUNSET(K) = PDAY ENDIF ENDDO SRBOT = SUNRISE(1) SRTOP = SUNRISE(NALT) SSBOT = SUNSET(1) SSTOP = SUNSET(NALT) ENDIF RETURN END C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ SUBROUTINE STEP INCLUDE 'PARAMETERS.INC' INCLUDE 'PLANET.INC' COMMON & /SWITCHES/ ZLAT , DELTAT , DELT1 , DAYSTART , & TIMESTART , TIMESTOP , OUTTIME , & MODE , NTIME , NITER , IDIFF , istepdump & /ONED3/ XAU & /ONED5/ SRTOP , SSTOP & /ONED8/ TIME & /ONED14/ DELTNEW & /ONED15/ THOUR & /ONED21/ SRBOT , SSBOT & /ONED35/ FLUX(NWAVE) & /ONED36/ DAY & /ONED56/ DAYNEW & /BLOCK29/ FLUX1AU(NWAVE) ! FROM SUBROUTINE PHOTO C CALCULATES TEMPORAL PARAMETERS. THIS ROUTINE IS USED ONLY FOR DOING A C DIURNAL RUN. C THE FOLLOWING SECTION ADJUSTS THE SIZE OF DELTHR DELTHR = DELTNEW/3600.0d0 IF (DELTHR .GE. DELT1/3600.0d0) THEN ! STEP SIZE > TERMINATOR STEP SIZE TNOW = THOUR TNEXT = THOUR + DELTHR TNEXT = JIDINT(TNEXT*60.0d0)/60.0d0 ! TIME TO NEAREST MINUTE TNOW = JIDINT(TNOW*60.0d0)/60.0d0 C IF PRESENTLY WITHIN OR GOING TO BE WITHIN AN HOUR BEFORE OR AFTER C SUNSET OR AN HOUR AFTER SUNRISE, SET TIME STEP EQUAL TO DELT1. IF ( ((TNEXT .GT. SRTOP) .AND. (TNEXT .LE. (SRBOT+1.0d0))) & .OR. & ((TNOW .GT. SRTOP) .AND. (TNOW .LT. (SRBOT+1.0d0))) & .OR. & ((TNOW .GT. (SSBOT-1.0d0)) .AND. (TNOW .LT. (SSTOP+1.0d0))) & .OR. & ((TNEXT .GT. (SSBOT-1.0d0)) .AND. & (TNEXT .LE. (SSTOP+1.0d0))) ) THEN DELTHR = DELT1/3600.0d0 ENDIF ENDIF IF (DELTHR .GE. 1.0d0) DELTHR = JIDINT(DELTHR) ! TRUNCATE TO A ! MULTIPLE OF ONE HOUR C NOW UPDATE TEMPORAL VARIABLES USING THE NEWLY COMPUTED DELTHR DELTNEW = 3600.0d0*DELTHR ! DELTA TIME IN SECONDS TIME = TIME + DELTHR ! NEW TIME VALUE C EVALUATE SOLAR ZENITH ANGLE AT MIDPOINT BETWEEN OLD AND NEW TIME VALUES C TAKE THE PLANETARY DAY INTO ACCOUNT IF(MODE .NE. 2 )THEN THOUR = DMOD(TIME,DBLE(PDAY)) - 0.5d0*DELTHR ENDIF DELDAY = DELTHR/PDAY ! DELTA DAY IN HOURS DAY = DAY + DELDAY DAYNEW = DAYNEW + DELDAY C RESET TIME TO NEXT DAY(YEAR) IF ONE DAY(YEAR) HAS ELAPSED C UPDATE ORBIT, SOLAR RADIATION, SUNRISE, AND SUNSET AFTER ONE DAY IF (DAY .GT. PYEAR) DAY = DAY - PYEAR IF (DAYNEW .GT. PDAY) THEN DAYNEW = DAYNEW - PDAY IF( IORBIT .EQ. 1 ) THEN ! IORBIT = 0 IS A PERPETUAL SEASON CALL ORBIT CALL DAYTIME ! CALCULATE NEW LENGTH OF DAY. C SCALE FLUX FOR DISTANCE FROM SUN. FAC = 1.0D0/(XAU**2) DO LAMBDA = 1 , NWAVE FLUX(LAMBDA) = FLUX1AU(LAMBDA)*FAC ENDDO ENDIF IF( IDIFF .EQ. 1 ) CALL DIFFLUX ! COMPUTE NEW SCATTERING MATRIX ENDIF RETURN END C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ SUBROUTINE TIMEHALF INCLUDE 'PARAMETERS.INC' INCLUDE 'PLANET.INC' LOGICAL ITERATE COMMON & /SWITCHES/ ZLAT , DELTAT , DELT1 , DAYSTART , & TIMESTART , TIMESTOP , OUTTIME , & MODE , NTIME , NITER , IDIFF , istepdump & /ONED2/ ITERATE & /ONED8/ TIME & /ONED14/ DELTNEW & /ONED25/ JUPDATE & /ONED36/ DAY & /ONED38/ DELTOLD , NHALFS & /ONED56/ DAYNEW C THIS SUBROUTINE HALVES THE TIME STEP IN AN ATTEMPT TO GET THE SOLUTION C TO CONVERGE AND RESETS TIME BACK TO ITS PREVIOUS VALUE BEFORE IT RAN C INTO CONVERGENCE PROBLEMS. IF (NHALFS .NE. 0) THEN NHALFS = NHALFS - 1 JUPDATE = 1 ! CALC IS IN TROUBLE- UPDATE PHOTOLYSIS OFTEN ITERATE = .TRUE. ! CONTINUE ITERATING IF( MODE .EQ. 1 ) THEN DELTHR = DELTNEW/3600.0d0 ! DELTAT IN UNITS OF HOURS C GO BACK TO PREVIOUS VALUES OF TIME TIME = TIME - DELTHR IF (TIME .LE. 0.0d0) TIME = TIME + PDAY DELDAY = DELTHR/PDAY DAY = DAY - DELDAY DAYNEW = DAYNEW - DELDAY DAYNEW = DMAX1(DAYNEW,0.0D0) IF (DAY .LE. 0.0d0) DAY = DAY + PYEAR DELTNEW = 0.5d0*DELTNEW ! *REDUCE TIME STEP* CALL STEP ! TAKE NEW TIME STEP ELSE CCC DELTNEW = DELTOLD + 0.5d0*(DELTNEW - DELTOLD) ! DIURNAL AVERAGE deltnew = 0.5d0*deltnew ENDIF ELSE ITERATE = .FALSE. ENDIF RETURN END