SUBROUTINE GEOM INCLUDE 'PARAMETERS.INC' INCLUDE 'PLANET.INC' COMMON & /BLOCK10/ Z(NALT) & /BLOCK27/ PI & /ONED15/ THOUR & /ONED32/ SUNRISE(NALT) , SUNSET(NALT) & /ONED33/ DS , DSRAY(NALT) , KTAN(NALT) , & RAYTAN(NALT,NALT) , ZTAN(NALT) , & CHIRAY(NALT) , CHITAN(NALT,NALT) , & SINSZA & /ONED34/ ANG1 , ANG2 & /ONED41/ SZA & /ONED50/ COSSZA C THIS SUBROUTINE COMPUTES SOLAR ZENITH ANGLES AND SLANT DISTANCES. C DIST(X,Y,SINSZA,COSSZA) IS THE DISTANCE COVERED BY A RAY AS IT TRAVERSES C THE LAYER BETWEEN ALTITUDES X AND Y AND INCLINED AT A SOLAR ZENITH ANGLE C WHOSE COSINE AND SINE ARE GIVEN BY COSSZA AND SINSZA. THIS FORMULA CAN BE C DERIVED BY USING THE LAW OF COSINES. DIST(X,Y,SINSZA,COSSZA) = -(RADIUS + X)*COSSZA + & DSQRT( (RADIUS + Y)*(RADIUS + Y) - & (RADIUS + X)*(RADIUS + X)*SINSZA*SINSZA ) COSSZA = ANG1 + ANG2*COS( (2.0D0*PI/PDAY)*(THOUR - 0.5D0*PDAY) ) SINSZA = DSQRT(1.0D0 - COSSZA*COSSZA) SZA = ACOS(COSSZA) C COMPUTE SLANT DISTANCES C THE SLANT DISTANCE IS ONLY WEAKLY DEPENDENT UPON THE ALTITUDE FOR SOLAR ZENITH C ANGLES LESS THAN 75 DEGREES. THIS DEPENDENCY WILL BE IGNORED. ALSO, THE SLANT C DISTANCE WILL BE CONVERTED TO CENTIMETERS. KMID = NALT/2 IF( SZA .LT. 75.0d0*PI/180.0d0 ) THEN DS = DIST(Z(KMID),Z(KMID+1),SINSZA,COSSZA)*1d5 ENDIF C COMPUTE ACCURATE DISTANCE FUNCTIONS USING RAY-TRACING. THESE VALUES C ARE ONLY NEEDED FOR A SOLAR ZENITH ANGLE > 75 DEGREES. WE MAKE THE C FOLLOWING APPROXIMATION: THE SLANT PATH FOR ALTITUDE KLAYER > K IS GIVEN C BY DSRAY(KLAYER-K+1). THUS, THE SLANT PATHS FOR A RAY REACHING C ANY ALTITUDE CAN BE OBTAINED BY A SIMPLE TRANSLATION OF DSRAY. IF( ( SZA .GE. 75.0*PI/180.0 ) .AND. ( SZA .LE. PI/2.0) ) THEN CHI = SZA DO K = 1 , NALT-1 CHIRAY(K) = CHI X = Z(K) Y = Z(K+1) DSRAY(K) = DIST(X,Y,SIN(CHI),COS(CHI)) DELTA = ASIN( DSRAY(K)*SIN(CHI)/(RADIUS + Y) ) DSRAY(K) = DSRAY(K)*1d5 ! DISTANCE IN CM CHI = CHI - DELTA ENDDO CHIRAY(NALT) = CHI ENDIF IF( (SZA .GE. PI/2.0d0) .AND. (ISPHERE .EQ. 1) ) THEN C DETERMINE TANGENT POINT ALTITUDE LAYER FOR SOLAR ZENITH ANGLES > PI/2 DO K = 1 , NALT C TEST TO SEE IF THE ALTITUDE POINT IS IN DAYLIGHT IF( (THOUR .GE. SUNRISE(K)) .AND. & (THOUR .LT. SUNSET(K)) )THEN C ZTAN IS THE TANGENT POINT ALTITUDE C THE TANGENT POINT IS IN THE (KTAN-1)TH ALTITUDE LAYER ZTAN(K) = (RADIUS + Z(K))*SINSZA - RADIUS KTANG = K DO WHILE( ZTAN(K) .LT. Z(KTANG) & .AND. ( KTANG .GT. 1 ) ) KTANG = KTANG - 1 ENDDO KTAN(K) = KTANG + 1 KTAN(K) = MIN0(NALT,KTAN(K)) C CALCULATE SLANT PATH DISTANCES FOR TANGENT POINT CALCULATIONS. C CALCULATE DISTANCE TRAVELED TO THE TANGENT POINT IN THE LAST LAYER C (THE ONE CONTAINING THE TANGENT POINT). K1 = KTAN(K) - 1 K2 = KTAN(K) CHI = PI/2.0d0 RAYTAN(K1,K) = SQRT( (RADIUS + Z(K2))*(RADIUS + Z(K2)) - & (RADIUS + ZTAN(K))*(RADIUS + ZTAN(K)) ) DELTA = ASIN( RAYTAN(K1,K)*SIN(CHI)/ & (RADIUS + Z(K2)) ) RAYTAN(K1,K) = 1d5*RAYTAN(K1,K) ! DIST IN CM CHI = CHI - DELTA C ADD ON REMAINING LAYERS DO KLAYER = K2 , NALT-1 CHITAN(KLAYER,K) = CHI X = Z(KLAYER) Y = Z(KLAYER+1) RAYTAN(KLAYER,K) = DIST(X,Y,SIN(CHI),COS(CHI)) DELTA = ASIN( RAYTAN(KLAYER,K)*SIN(CHI)/ & (RADIUS + Y) ) RAYTAN(KLAYER,K) = 1d5*RAYTAN(KLAYER,K) CHI = CHI - DELTA ENDDO CHITAN(NALT,K) = CHI ENDIF ! END TEST ON DAYLIGHT ENDDO ! END LOOP OVER ALTITUDE C CALCULATE COMPLEMENTARY SLANT PATH DISTANCES. THE COMPLEMENTARY C SOLAR ZENITH ANGLE IS PI - CHI. CHI = PI - ACOS(COSSZA) ! COMPLEMENTARY SZA DO K = 1 , NALT-1 IF( (THOUR .GE. SUNRISE(K)) .AND. (THOUR .LT. SUNSET(K)) )THEN CHIRAY(K) = CHI X = Z(K) Y = Z(K+1) DSRAY(K) = DIST(X,Y,SIN(CHI),COS(CHI)) DELTA = ASIN( DSRAY(K)*SIN(CHI)/ & (RADIUS + Y) ) DSRAY(K) = 1d5*DSRAY(K) ! DISTANCE IN CM CHI = CHI - DELTA ENDIF ENDDO CHIRAY(NALT) = CHI ENDIF CALL SLANTDIU ! CALCULATE SLANT PATH COLUMN DENSITIES RETURN END C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ SUBROUTINE GEOMAVG INCLUDE 'PARAMETERS.INC' INCLUDE 'PLANET.INC' COMMON & /BLOCK10/ Z(NALT) & /BLOCK18/ CONC(MB,NALT) & /BLOCK16/ COSSZA(NTRAP) & /BLOCK17/ NANGLE(NALT) & /BLOCK20/ DS(NTRAP) , DSRAY(NALT,NTRAP) , & KTAN(NTRAP,NALT) , & RAYTAN(NALT,NTRAP,NALT) , & ZTAN(NTRAP,NALT) , & SINSZA(NTRAP) , CHIRAY(NALT,NTRAP) , & CHITAN(NALT,NTRAP,NALT) & /BLOCK27/ PI & /BLOCK30/ NPI2 & /ONED34/ ANG1 , ANG2 & /ONED49/ N75 C THIS SUBROUTINE COMPUTES SOLAR ZENITH ANGLES AND SLANT DISTANCES. C THE DIURNAL AVERAGING PORTION OF THE MAIN PROGRAM REQUIRES THESE VALUES. C BY PRECOMPUTING THESE VALUES AT THE INITIALIZATION PHASE OF THE PROGRAM C WE AVOID HAVING TO MAKE REDUNTANT (AND COSTLY) COMPUTATIONS. C COSSZA(KOUNT) AND SINSZA(KOUNT) ARE THE COSINES AND SINES OF THE SOLAR C ZENITH ANGLE. C DIST(X,Y,SINSZA,COSSZA) IS THE DISTANCE COVERED BY A RAY AS IT TRAVERSES C THE LAYER BETWEEN ALTITUDES X AND Y AND INCLINED AT A SOLAR ZENITH ANGLE C WHOSE COSINE AND SINE ARE GIVEN BY COSSZA AND SINSZA. THIS FORMULA CAN BE C DERIVED BY USING THE LAW OF COSINES. DIST(Xa,Ya,SINSZAa,COSSZAa) = -(RADIUS + Xa)*COSSZAa + & DSQRT( (RADIUS + Ya)*(RADIUS + Ya) - & (RADIUS + Xa)*(RADIUS + Xa)*SINSZAa*SINSZAa ) C COMPUTE SOLAR ZENITH ANGLES DO NANG = 1 , NTRAP OMEGA = (NANG-1.0d0)*(PI/(NTRAP-1.0d0)) COSSZA(NANG) = ANG1 + ANG2*COS(OMEGA) C ADJUST FOR POSSIBLE ROUNDOFF ERRORS. IT IS WISE TO CHOOSE NTRAP TO BE C AN EVEN NUMBER IN ORDER TO MINIMIZE THIS TYPE OF PROBLEM. IF( ABS(COSSZA(NANG)) .LT. 1d-7 )THEN COSSZA(NANG) = 0.0d0 ENDIF SINSZA(NANG) = SQRT(1.0d0 - COSSZA(NANG)*COSSZA(NANG)) ENDDO C FIND NUMBER OF SOLAR ZENITH ANGLES THAT ARE LESS THAN PI/2 NANG = 1 DO WHILE( (NANG .LE. NTRAP) .AND. & (COSSZA(NANG) .GE. 0.0d0) ) NANG = NANG + 1 ENDDO NPI2 = NANG -1 C DETERMINE THE NUMBER OF SOLAR ZENITH ANGLES LESS THAN 75 DEGREES. C ANGLES LESS THAN 75 DEGREES WILL USE THE SECANT APPROXIMATION. NANG = 1 DO WHILE( (NANG .LE. NPI2) .AND. & (SINSZA(NANG) .LT. 0.966d0) ) NANG = NANG + 1 ENDDO N75 = NANG - 1 C COMPUTE NUMBER OF ANGLES BEYOND PI/2 WHICH ARE NOT IN DARKNESS DO K = 1 , NALT NANG = NPI2 + 1 DO WHILE( (SINSZA(NANG) .GE. (RADIUS/(RADIUS + Z(K))) ) & .AND. (NANG .LE. NTRAP) ) NANG = NANG + 1 ENDDO NANGLE(K) = NANG - 1 ! TOTAL NUMBER OF DAYLIGHT ANGLES ENDDO C COMPUTE SLANT DISTANCES AT ANGLE NANG C THE SLANT DISTANCE IS ONLY WEAKLY DEPENDENT UPON THE ALTITUDE. THIS C DEPENDENCY WILL BE IGNORED. ALSO, THE SLANT DISTANCE WILL BE CONVERTED C TO CENTIMETERS. KMID = NALT/2 IF( ISPHERE .EQ. 0 ) THEN NEND = NPI2 ELSE NEND = N75 ENDIF DO NANG = 1 , NEND DS(NANG) = (1d5)*DIST(Z(KMID),Z(KMID+1), & SINSZA(NANG),COSSZA(NANG)) ENDDO C COMPUTE ACCURATE DISTANCE FUNCTIONS USING RAY-TRACING. THESE VALUES C ARE ONLY NEEDED FOR SOLAR ZENITH ANGLE > 75 DEGREES. WE MAKE THE C FOLLOWING APPROXIMATION: THE SLANT PATH FOR ALTITUDE KLAYER > K IS GIVEN C BY DSRAY(KLAYER-K+1,NANG). THUS, THE SLANT PATHS FOR A RAY REACHING C ANY ALTITUDE CAN BE OBTAINED BY A SIMPLE TRANSLATION OF DSRAY. IF( ISPHERE .EQ. 1) THEN ! RAY-TRACE CALCULATION N1 = N75 + 1 N2 = NPI2 DO NANG = N1 , N2 CHI = ASIN(SINSZA(NANG)) ! SZA IN RADIANS DO K = 1 , NALT-1 CHIRAY(K,NANG) = CHI X = Z(K) Y = Z(K+1) DSRAY(K,NANG) = DIST(X,Y,SIN(CHI),COS(CHI)) DELTA = ASIN( DSRAY(K,NANG)*SIN(CHI)/ & (RADIUS + Y) ) DSRAY(K,NANG) = (1d5)*DSRAY(K,NANG) ! DISTANCE IN CM CHI = CHI - DELTA ENDDO CHIRAY(NALT,NANG) = CHI ENDDO C DETERMINE TANGENT POINT ALTITUDE LAYER FOR SOLAR ZENITH ANGLES > PI/2 DO K = 1 , NALT N1 = NPI2 + 1 N2 = NANGLE(K) DO NANG = N1 , N2 C ZTAN IS THE TANGENT POINT ALTITUDE C THE TANGENT POINT IS IN THE (KTAN-1)TH ALTITUDE LAYER ZTAN(NANG,K) = (RADIUS + Z(K))*SINSZA(NANG) - RADIUS KTANG = K DO WHILE( ZTAN(NANG,K) .LT. Z(KTANG) & .AND. ( KTANG .GT. 1 ) ) KTANG = KTANG - 1 ENDDO KTAN(NANG,K) = KTANG + 1 KTAN(NANG,K) = MIN0(NALT,KTAN(NANG,K)) C CALCULATE SLANT PATH DISTANCES FOR TANGENT POINT CALCULATIONS. C CALCULATE DISTANCE TRAVELED TO THE TANGENT POINT IN THE LAST LAYER C (THE ONE CONTAINING THE TANGENT POINT) K1 = KTAN(NANG,K) - 1 K2 = KTAN(NANG,K) CHI = PI/2.0d0 RAYTAN(K1,NANG,K) = SQRT( & (RADIUS + Z(K2))*(RADIUS + Z(K2)) - & (RADIUS + ZTAN(NANG,K))*(RADIUS + ZTAN(NANG,K))) DELTA = ASIN( RAYTAN(K1,NANG,K)*SIN(CHI)/ & (RADIUS + Z(K2)) ) RAYTAN(K1,NANG,K) = (1d5)*RAYTAN(K1,NANG,K) ! DIST IN CM CHI = CHI - DELTA C ADD ON REMAINING LAYERS DO KLAYER = K2 , NALT-1 CHITAN(KLAYER,NANG,K) = CHI X = Z(KLAYER) Y = Z(KLAYER+1) RAYTAN(KLAYER,NANG,K) = DIST(X,Y,SIN(CHI),COS(CHI)) DELTA = ASIN( RAYTAN(KLAYER,NANG,K)*SIN(CHI)/ & (RADIUS + Y) ) RAYTAN(KLAYER,NANG,K) = (1d5)*RAYTAN(KLAYER,NANG,K) CHI = CHI - DELTA ENDDO CHITAN(NALT,NANG,K) = CHI ENDDO ENDDO C CALCULATE COMPLEMENTARY SLANT PATH DISTANCES. THE COMPLEMENTARY C SOLAR ZENITH ANGLE IS PI - CHI. N1 = NPI2 + 1 N2 = NANGLE(NALT) ! MAXIMUM NUMBER OF DAYLIGHT ANGLES DO NANG = N1 , N2 CHI = PI - ACOS(COSSZA(NANG)) ! COMPLEMENTARY SZA DO K = 1 , NALT-1 CHIRAY(K,NANG) = CHI X = Z(K) Y = Z(K+1) DSRAY(K,NANG) = DIST(X,Y,SIN(CHI),COS(CHI)) DELTA = ASIN( DSRAY(K,NANG)*SIN(CHI)/ & (RADIUS + Y) ) DSRAY(K,NANG) = (1d5)*DSRAY(K,NANG) ! DISTANCE IN CM CHI = CHI - DELTA ENDDO CHIRAY(NALT,NANG) = CHI ENDDO ENDIF ! END RAY-TRACE COMPUTATIONS RETURN END C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ SUBROUTINE SLANTDIU INCLUDE 'PARAMETERS.INC' COMMON & /BLOCK10/ Z(NALT) & /BLOCK11/ T(NALT) & /BLOCK14/ MTAU(MB) & /BLOCK15/ NTAU & /BLOCK18/ CONC(MB,NALT) & /BLOCK27/ PI & /BLOCK28/ SCALE(MB) & /BLOCK31/ SCALEH(MB) & /ONED15/ THOUR & /ONED32/ SUNRISE(NALT) , SUNSET(NALT) & /ONED33/ DS , DSRAY(NALT) , KTAN(NALT) , & RAYTAN(NALT,NALT) , ZTAN(NALT) , & CHIRAY(NALT) , CHITAN(NALT,NALT) , & SINSZA & /ONED41/ SZA & /ONED47/ COLUMN(NALT,MB) & /ONED50/ COSSZA C THIS SUBROUTINE COMPUTES THE COLUMN DEPTHS FOR THE OPTICALLY THICK SPECIES. C CALCULATE OPTICAL DEPTH AT TOP OF MODEL. THE NUMBER DENSITY C ABOVE THE TOP LEVEL IS APPROXIMATED BY N(Z) = N(NALT)*EXP(-Z/H) ; H = KT/MG C SINCE THE ATMOSPHERE ABOVE THE UPPER LIMIT IS ASSUMED TO BE EXPONENTIALY C DISTRIBUTED, THE CHAPMAN FUNCTION CAN BE USED TO COMPUTE THE OPTICAL DEPTH. C COMPUTE SCALE HEIGHTS AT TOP OF MODEL DO N = 1 , NTAU M = MTAU(N) SCALEH(M) = SCALE(M)*T(NALT) ! SCALE HEIGHT C CONVERT TO CENTIMETERS (CONC IS IN UNITS OF CM^-3). 1 AMU = 1.66E-27 KG SCALEH(M) = (100.0d0/1.66d-27)*SCALEH(M) ENDDO C PLANAR GEOMETRY CALCULATION IF( SZA .LT. 75.0d0*PI/180.0d0 ) THEN COSTOP = COSSZA ! SZA AT TOP SINTOP = SINSZA ! SZA AT TOP COLUMN(NALT,MTAU(1)) = CONC(MTAU(1),NALT)*SCALEH(MTAU(1))* & CHAP(COSTOP,SINTOP,SCALEH(MTAU(1))) COLUMN(NALT,MTAU(2)) = CONC(MTAU(2),NALT)*SCALEH(MTAU(2))* & CHAP(COSTOP,SINTOP,SCALEH(MTAU(2))) COLUMN(NALT,MTAU(3)) = CONC(MTAU(3),NALT)*SCALEH(MTAU(3))* & CHAP(COSTOP,SINTOP,SCALEH(MTAU(3))) DO K = NALT-1 , 1 , -1 COLUMN(K,MTAU(1)) = COLUMN(K+1,MTAU(1)) + (0.5D0*DS)* & (CONC(MTAU(1),K) + CONC(MTAU(1),K+1)) COLUMN(K,MTAU(2)) = COLUMN(K+1,MTAU(2)) + (0.5D0*DS)* & (CONC(MTAU(2),K) + CONC(MTAU(2),K+1)) COLUMN(K,MTAU(3)) = COLUMN(K+1,MTAU(3)) + (0.5D0*DS)* & (CONC(MTAU(3),K) + CONC(MTAU(3),K+1)) ENDDO ENDIF C COMPUTE COLUMN DEPTHS FOR 75 < SZA < 90 USING AN ACCURATE FORMULA IF( (SZA .GE. 75.0*PI/180.0 ) .AND. ( SZA .LE. PI/2.0 ) ) THEN C COMPUTE COLUMN DEPTHS AT TOP OF MODEL. COSTOP = COS(CHIRAY(NALT)) ! SZA AT TOP SINTOP = SIN(CHIRAY(NALT)) ! SZA AT TOP COLUMN(NALT,MTAU(1)) = CONC(MTAU(1),NALT)*SCALEH(MTAU(1))* & CHAP(COSTOP,SINTOP,SCALEH(MTAU(1))) COLUMN(NALT,MTAU(2)) = CONC(MTAU(2),NALT)*SCALEH(MTAU(2))* & CHAP(COSTOP,SINTOP,SCALEH(MTAU(2))) COLUMN(NALT,MTAU(3)) = CONC(MTAU(3),NALT)*SCALEH(MTAU(3))* & CHAP(COSTOP,SINTOP,SCALEH(MTAU(3))) DO K = 1 , NALT-1 C SUM UP CONTRIBUTIONS DUE TO DIFFERENTIAL COLUMN DEPTHS SUM1 = 0.0D0 SUM2 = 0.0D0 SUM3 = 0.0D0 DO KLAYER = K , NALT-1 DIST = 0.5D0*DSRAY(KLAYER-K+1) SUM1 = SUM1 + DIST*(CONC(MTAU(1),KLAYER) + & CONC(MTAU(1),KLAYER+1)) SUM2 = SUM2 + DIST*(CONC(MTAU(2),KLAYER) + & CONC(MTAU(2),KLAYER+1)) SUM3 = SUM3 + DIST*(CONC(MTAU(3),KLAYER) + & CONC(MTAU(3),KLAYER+1)) ENDDO C ADD ON CONTRIBUTION FROM TOP OF ATMOSPHERE COLUMN(K,MTAU(1)) = SUM1 + COLUMN(NALT,MTAU(1)) COLUMN(K,MTAU(2)) = SUM2 + COLUMN(NALT,MTAU(2)) COLUMN(K,MTAU(3)) = SUM3 + COLUMN(NALT,MTAU(3)) ENDDO ENDIF C CALCULATE COLUMN DEPTH FOR SZA > PI/2 IF( (SZA .GT. PI/2.0d0) .AND. ( ISPHERE .EQ. 1 ) ) THEN C COMPUTE COLUMN DEPTHS AT TOP OF MODEL FOR THE COMPLEMENTARY RAY AT PI - CHI COSTOP = COS(CHIRAY(NALT)) ! SZA AT TOP SINTOP = SIN(CHIRAY(NALT)) ! SZA AT TOP COLUMN(NALT,MTAU(1)) = CONC(MTAU(1),NALT)*SCALEH(MTAU(1))* & CHAP(COSTOP,SINTOP,SCALEH(MTAU(1))) COLUMN(NALT,MTAU(2)) = CONC(MTAU(2),NALT)*SCALEH(MTAU(2))* & CHAP(COSTOP,SINTOP,SCALEH(MTAU(2))) COLUMN(NALT,MTAU(3)) = CONC(MTAU(3),NALT)*SCALEH(MTAU(3))* & CHAP(COSTOP,SINTOP,SCALEH(MTAU(3))) DO K = 1 , NALT IF( (THOUR .GE. SUNRISE(K)) .AND. & (THOUR .LE. SUNSET(K)) ) THEN C COMPUTE COMPLEMENTARY COLUMN DEPTH AT PI - CHI SUM1 = 0.0D0 SUM2 = 0.0D0 SUM3 = 0.0D0 DO KLAYER = K , NALT-1 DIST = 0.5D0*DSRAY(KLAYER-K+1) SUM1 = SUM1 + DIST*(CONC(MTAU(1),KLAYER) + & CONC(MTAU(1),KLAYER+1)) SUM2 = SUM2 + DIST*(CONC(MTAU(2),KLAYER) + & CONC(MTAU(2),KLAYER+1)) SUM3 = SUM3 + DIST*(CONC(MTAU(3),KLAYER) + & CONC(MTAU(3),KLAYER+1)) ENDDO C ADD ON CONTRIBUTION FROM TOP OF ATMOSPHERE COLUMN(K,MTAU(1)) = SUM1 + COLUMN(NALT,MTAU(1)) COLUMN(K,MTAU(2)) = SUM2 + COLUMN(NALT,MTAU(2)) COLUMN(K,MTAU(3)) = SUM3 + COLUMN(NALT,MTAU(3)) K2 = KTAN(K) ! TANGENT POINT ALTITUDE K1 = K2 - 1 C SUM UP CONTRIBUTIONS DUE TO DIFFERENTIAL COLUMN DEPTHS FOR THE TANGENT RAY SUM1 = 0.0D0 SUM2 = 0.0D0 SUM3 = 0.0D0 DO KLAYER = K2 , NALT-1 DIST = 0.5D0*RAYTAN(KLAYER,K) SUM1 = SUM1 + DIST*(CONC(MTAU(1),KLAYER) + & CONC(MTAU(1),KLAYER+1)) SUM2 = SUM2 + DIST*(CONC(MTAU(2),KLAYER) + & CONC(MTAU(2),KLAYER+1)) SUM3 = SUM3 + DIST*(CONC(MTAU(3),KLAYER) + & CONC(MTAU(3),KLAYER+1)) ENDDO C CALCULATE COLUMN DEPTHS IN TANGENT LAYER AND INCLUDE IN TOTAL COLUMN DEPTH CONC1 = ( CONC(MTAU(1),K2) - CONC(MTAU(1),K1) )* & ( ZTAN(K) - Z(K1) )/(Z(K2) - Z(K1)) + & CONC(MTAU(1),K1) CONC2 = ( CONC(MTAU(2),K2) - CONC(MTAU(2),K1) )* & ( ZTAN(K) - Z(K1) )/(Z(K2) - Z(K1)) + & CONC(MTAU(2),K1) CONC3 = ( CONC(MTAU(3),K2) - CONC(MTAU(3),K1) )* & ( ZTAN(K) - Z(K1) )/(Z(K2) - Z(K1)) + & CONC(MTAU(3),K1) SUM1 = SUM1 + 0.5D0*RAYTAN(K1,K)*(CONC1 + CONC(MTAU(1),K2)) SUM2 = SUM2 + 0.5D0*RAYTAN(K1,K)*(CONC2 + CONC(MTAU(2),K2)) SUM3 = SUM3 + 0.5D0*RAYTAN(K1,K)*(CONC3 + CONC(MTAU(3),K2)) C ADD ON CONTRIBUTION FROM TOP OF ATMOSPHERE COSTOP = COS(CHITAN(NALT,K)) ! SZA AT TOP SINTOP = SIN(CHITAN(NALT,K)) ! SZA AT TOP SUM1 = SUM1 + CONC(MTAU(1),NALT)*SCALEH(MTAU(1))* & CHAP(COSTOP,SINTOP,SCALEH(MTAU(1))) SUM2 = SUM2 + CONC(MTAU(2),NALT)*SCALEH(MTAU(2))* & CHAP(COSTOP,SINTOP,SCALEH(MTAU(2))) SUM3 = SUM3 + CONC(MTAU(3),NALT)*SCALEH(MTAU(3))* & CHAP(COSTOP,SINTOP,SCALEH(MTAU(3))) C CALCULATE COLUMN DEPTHS FOR SOLAR ZENITH ANGLES > PI/2. THE COLUMN DEPTH IS C COMPUTED BY SUBTRACTING THE COLUMN DEPTH AT ANGLE PI - CHI FROM TWICE THE C COLUMN DEPTH TO THE TANGENT POINT. COLUMN(K,MTAU(1)) = 2.0D0*SUM1 - COLUMN(K,MTAU(1)) COLUMN(K,MTAU(2)) = 2.0D0*SUM2 - COLUMN(K,MTAU(2)) COLUMN(K,MTAU(3)) = 2.0D0*SUM3 - COLUMN(K,MTAU(3)) ENDIF ENDDO ENDIF RETURN END C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ C------------------------------------------------------------------------------ SUBROUTINE SLANTAVG INCLUDE 'PARAMETERS.INC' COMMON & /BLOCK10/ Z(NALT) & /BLOCK11/ T(NALT) & /BLOCK14/ MTAU(MB) & /BLOCK15/ NTAU & /BLOCK16/ COSSZA(NTRAP) & /BLOCK17/ NANGLE(NALT) & /BLOCK18/ CONC(MB,NALT) & /BLOCK20/ DS(NTRAP) , DSRAY(NALT,NTRAP) , & KTAN(NTRAP,NALT) , & RAYTAN(NALT,NTRAP,NALT) , & ZTAN(NTRAP,NALT) , & SINSZA(NTRAP) , CHIRAY(NALT,NTRAP) , & CHITAN(NALT,NTRAP,NALT) & /BLOCK28/ SCALE(MB) & /BLOCK30/ NPI2 & /BLOCK31/ SCALEH(MB) & /ONED46/ COLUMN(NTRAP,NALT,MB) & /ONED49/ N75 C THIS SUBROUTINE COMPUTES THE COLUMN DEPTHS FOR THE OPTICALLY THICK SPECIES. C CALCULATE OPTICAL DEPTH AT TOP OF MODEL. THE NUMBER DENSITY C ABOVE THE TOP LEVEL IS APPROXIMATED BY N(Z) = N(NALT)*EXP(-Z/H) ; H = KT/MG C SINCE THE ATMOSPHERE ABOVE THE UPPER LIMIT IS ASSUMED TO BE EXPONENTIALY C DISTRIBUTED, THE CHAPMAN FUNCTION CAN BE USED TO COMPUTE THE OPTICAL DEPTH. IF( ISPHERE .EQ. 1 ) THEN NEND = N75 ELSE NEND = NPI2 ENDIF DO N = 1 , NTAU M = MTAU(N) SCALEH(M) = SCALE(M)*T(NALT) ! SCALE HEIGHT C CONVERT TO CENTIMETERS (CONC IS IN UNITS OF CM^-3). 1 AMU = 1.66E-27 KG SCALEH(M) = (100.0d0/1.66d-27)*SCALEH(M) DO NANG = 1 , NEND COSTOP = COSSZA(NANG) ! SZA AT TOP SINTOP = SINSZA(NANG) ! SZA AT TOP COLUMN(NANG,NALT,M) = CONC(M,NALT)*SCALEH(M)* & CHAP(COSTOP,SINTOP,SCALEH(M)) ENDDO ENDDO C PLANAR GEOMETRY CALCULATION DO K = NALT-1 , 1 , -1 DO NANG = 1 , NEND DIST = 0.5D0*DS(NANG) COLUMN(NANG,K,MTAU(1)) = COLUMN(NANG,K+1,MTAU(1)) + & DIST* & (CONC(MTAU(1),K) + CONC(MTAU(1),K+1)) COLUMN(NANG,K,MTAU(2)) = COLUMN(NANG,K+1,MTAU(2)) + & DIST* & (CONC(MTAU(2),K) + CONC(MTAU(2),K+1)) COLUMN(NANG,K,MTAU(3)) = COLUMN(NANG,K+1,MTAU(3)) + & DIST* & (CONC(MTAU(3),K) + CONC(MTAU(3),K+1)) ENDDO ENDDO IF( ISPHERE .EQ. 1 ) THEN C COMPUTE COLUMN DEPTHS AT TOP OF MODEL. THIS BLOCK ALSO CALCULATES C THE INITIAL COLUMN DEPTHS USED FOR THE CASE WHERE SZA > PI/2 N1 = N75 + 1 N2 = NANGLE(NALT) DO NANG = N1 , N2 COSTOP = COS(CHIRAY(NALT,NANG)) ! SZA AT TOP SINTOP = SIN(CHIRAY(NALT,NANG)) ! SZA AT TOP COLUMN(NANG,NALT,MTAU(1)) = CONC(MTAU(1),NALT)* & SCALEH(MTAU(1))* & CHAP(COSTOP,SINTOP,SCALEH(MTAU(1))) COLUMN(NANG,NALT,MTAU(2)) = CONC(MTAU(2),NALT)* & SCALEH(MTAU(2))* & CHAP(COSTOP,SINTOP,SCALEH(MTAU(2))) COLUMN(NANG,NALT,MTAU(3)) = CONC(MTAU(3),NALT)* & SCALEH(MTAU(3))* & CHAP(COSTOP,SINTOP,SCALEH(MTAU(3))) ENDDO C COMPUTE COLUMN DEPTHS FOR 75 < SZA < 90 AND FOR THE COMPLEMENTARY ANGLES IF C SZA > 90 DO K = 1 , NALT-1 DO NANG = N75 + 1 , NANGLE(K) C SUM UP CONTRIBUTIONS DUE TO DIFFERENTIAL COLUMN DEPTHS SUM1 = 0.0D0 SUM2 = 0.0D0 SUM3 = 0.0D0 DO KLAYER = K , NALT-1 DIST = 0.5D0*DSRAY(KLAYER-K+1,NANG) SUM1 = SUM1 + DIST*(CONC(MTAU(1),KLAYER) + & CONC(MTAU(1),KLAYER+1)) SUM2 = SUM2 + DIST*(CONC(MTAU(2),KLAYER) + & CONC(MTAU(2),KLAYER+1)) SUM3 = SUM3 + DIST*(CONC(MTAU(3),KLAYER) + & CONC(MTAU(3),KLAYER+1)) ENDDO C ADD ON CONTRIBUTION FROM TOP OF ATMOSPHERE COLUMN(NANG,K,MTAU(1)) = SUM1 + COLUMN(NANG,NALT,MTAU(1)) COLUMN(NANG,K,MTAU(2)) = SUM2 + COLUMN(NANG,NALT,MTAU(2)) COLUMN(NANG,K,MTAU(3)) = SUM3 + COLUMN(NANG,NALT,MTAU(3)) ENDDO ENDDO C COMPUTES COLUMN DEPTHS FOR SZA > 90 DO K = 1 , NALT-1 DO NANG = NPI2 + 1 , NANGLE(K) K2 = KTAN(NANG,K) K1 = K2 - 1 C CALCULATE COLUMN DEPTH FOR TANGENT RAY C SUM UP CONTRIBUTIONS DUE TO DIFFERENTIAL COLUMN DEPTHS SUM1 = 0.0D0 SUM2 = 0.0D0 SUM3 = 0.0D0 DO KLAYER = K2 , NALT-1 DIST = 0.5D0*RAYTAN(KLAYER,NANG,K) SUM1 = SUM1 + DIST*(CONC(MTAU(1),KLAYER) + & CONC(MTAU(1),KLAYER+1)) SUM2 = SUM2 + DIST*(CONC(MTAU(2),KLAYER) + & CONC(MTAU(2),KLAYER+1)) SUM3 = SUM3 + DIST*(CONC(MTAU(3),KLAYER) + & CONC(MTAU(3),KLAYER+1)) ENDDO C CALCULATE COLUMN DEPTHS IN TANGENT LAYER AND INCLUDE IN TOTAL COLUMN DEPTH CONC1 = ( CONC(MTAU(1),K2) - CONC(MTAU(1),K1) )* & ( ZTAN(NANG,K) - Z(K1) )/(Z(K2) - Z(K1)) + & CONC(MTAU(1),K1) CONC2 = ( CONC(MTAU(2),K2) - CONC(MTAU(2),K1) )* & ( ZTAN(NANG,K) - Z(K1) )/(Z(K2) - Z(K1)) + & CONC(MTAU(2),K1) CONC3 = ( CONC(MTAU(3),K2) - CONC(MTAU(3),K1) )* & ( ZTAN(NANG,K) - Z(K1) )/(Z(K2) - Z(K1)) + & CONC(MTAU(3),K1) SUM1 = SUM1 + 0.5D0*RAYTAN(K1,NANG,K)* & (CONC1 + CONC(MTAU(1),K2)) SUM2 = SUM2 + 0.5D0*RAYTAN(K1,NANG,K)* & (CONC2 + CONC(MTAU(2),K2)) SUM3 = SUM3 + 0.5D0*RAYTAN(K1,NANG,K)* & (CONC3 + CONC(MTAU(3),K2)) C ADD ON CONTRIBUTION FROM TOP OF ATMOSPHERE COSTAN = COS(CHITAN(NALT,NANG,K)) ! SZA AT TOP SINTAN = SIN(CHITAN(NALT,NANG,K)) ! SZA AT TOP SUM1 = SUM1 + CONC(MTAU(1),NALT)*SCALEH(MTAU(1))* & CHAP(COSTAN,SINTAN,SCALEH(MTAU(1))) SUM2 = SUM2 + CONC(MTAU(2),NALT)*SCALEH(MTAU(2))* & CHAP(COSTAN,SINTAN,SCALEH(MTAU(2))) SUM3 = SUM3 + CONC(MTAU(3),NALT)*SCALEH(MTAU(3))* & CHAP(COSTAN,SINTAN,SCALEH(MTAU(3))) C CALCULATE COLUMN DEPTHS FOR SOLAR ZENITH ANGLES > PI/2. THE COLUMN DEPTH IS C COMPUTED BY SUBTRACTING THE COLUMN DEPTH AT ANGLE PI - CHI FROM TWICE THE C COLUMN DEPTH TO THE TANGENT POINT. COLUMN(NANG,K,MTAU(1)) = 2.0D0*SUM1 - COLUMN(NANG,K,MTAU(1)) COLUMN(NANG,K,MTAU(2)) = 2.0D0*SUM2 - COLUMN(NANG,K,MTAU(2)) COLUMN(NANG,K,MTAU(3)) = 2.0D0*SUM3 - COLUMN(NANG,K,MTAU(3)) ENDDO ENDDO ENDIF RETURN END