SUBROUTINE DIFFLUX INCLUDE 'PARAMETERS.INC' INCLUDE 'PLANET.INC' PARAMETER (SMALL = 1E-20) COMMON & /BLOCK6/ MCROSS(MB) & /BLOCK10/ Z(NALT) & /BLOCK12/ CROSS(NWAVE,MB) & /BLOCK14/ MTAU(MB) & /BLOCK18/ CONC(MB,NALT) & /BLOCK31/ SCALEH(MB) ! COMPUTED IN SLANTAVG/SLANTDIU & /BLOCK35/ INDX(NALT,NDIFWAV) , & A(NALT,NALT,NDIFWAV) & /ONED51/ ALBEDO(NWAVE) & /ONED52/ TAUVERT(NDIFWAV,NALT) , TOTALV(NDIFWAV) DIMENSION TAURAY(NALT,NDIFWAV) , TAUABS(NALT,NDIFWAV) , & COLVERT(NALT,MB) , DTAU(NALT,NDIFWAV) , & DTAURAY(NALT,NDIFWAV) , SCALE(NALT) , & RAYV(NDIFWAV) , ABSV(NDIFWAV) C THIS SUBROUTINE COMPUTES THE DIFFUSE FLUX DUE TO GROUND ALBEDO AND RAYLEIGH C SCATTERING. THIS VALUE IS ADDED ONTO THE SOLAR FLUX IN THE MAIN PROGRAM. C THE INDICES OF THE ABSORBING SPECIES HAVE BEEN EXPLICITLY WRITTEN OUT IN C ORDER TO IMPROVE THE EXECUTION SPEED. C THE FORMULAS WHICH ARE USED BELOW WERE DERIVED IN "Radiation Field in the C Troposphere and Stratosphere..." , C Planetary Space Science , Vol 30 , No. 9 , pp 923-933 C THESE FORMULAS ASSUME THAT THE SOURCE FUNCTION IS CONSTANT BETWEEN LAYERS C AND THUS EVALUATE THE SOURCE FUNCTION AT THE *MIDPOINT* OF A LAYER. C THUS, THE VERTICAL OPTICAL DEPTHS AND SCATTERING MATRIX WHICH ARE DERIVED C BELOW ALL APPLY TO THE *MIDPOINTS*. C COMPUTE COLUMN DEPTHS AT TOP OF MODEL FOR RAYLEIGH SCATTERING SPECIES AND C ABSORBING SPECIES COLVERT(NALT,MTAU(1)) = CONC(MTAU(1),NALT)*SCALEH(MTAU(1)) COLVERT(NALT,MTAU(2)) = CONC(MTAU(2),NALT)*SCALEH(MTAU(2)) COLVERT(NALT,MTAU(3)) = CONC(MTAU(3),NALT)*SCALEH(MTAU(3)) DO K = NALT-1 , 1 , -1 DIST = (0.5D5)*( Z(K+1) - Z(K) ) COLVERT(K,MTAU(1)) = DIST*(CONC(MTAU(1),K+1) + & CONC(MTAU(1),K)) + & COLVERT(K+1,MTAU(1)) COLVERT(K,MTAU(2)) = DIST*(CONC(MTAU(2),K+1) + & CONC(MTAU(2),K)) + & COLVERT(K+1,MTAU(2)) COLVERT(K,MTAU(3)) = DIST*(CONC(MTAU(3),K+1) + & CONC(MTAU(3),K)) + & COLVERT(K+1,MTAU(3)) ENDDO C COMPUTE DIFFERENTIAL OPTICAL DEPTHS AND VERTICAL OPTICAL DEPTHS. C NOTICE THAT THE VERTICAL OPTICAL DEPTHS ARE EVALUATED AT THE MIDPOINT OF THE C LAYERS. C TAURAY IS THE OPTICAL DEPTH DUE TO RAYLEIGH SCATTERING ONLY DO LAMBDIF = 1 , NDIFWAV LAMBDA = LOWDIF + LAMBDIF - 1 DO K = 1 , NALT-1 DTAURAY(K,LAMBDIF) = COLVERT(K,MTAU(3))* & CROSS(LAMBDA,MCROSS(3)) - & COLVERT(K+1,MTAU(3))* & CROSS(LAMBDA,MCROSS(3)) DTAUABS = ( COLVERT(K,MTAU(1))*CROSS(LAMBDA,MCROSS(1)) + & COLVERT(K,MTAU(2))*CROSS(LAMBDA,MCROSS(2)) ) - & ( COLVERT(K+1,MTAU(1))*CROSS(LAMBDA,MCROSS(1)) + & COLVERT(K+1,MTAU(2))*CROSS(LAMBDA,MCROSS(2)) ) TAURAY(K,LAMBDIF) = COLVERT(K,MTAU(3))* & CROSS(LAMBDA,MCROSS(3)) - & 0.5D0*DTAURAY(K,LAMBDIF) TAUABS(K,LAMBDIF) = COLVERT(K,MTAU(1))* & CROSS(LAMBDA,MCROSS(1)) + & COLVERT(K,MTAU(2))* & CROSS(LAMBDA,MCROSS(2)) - & 0.5D0*DTAUABS TAUVERT(LAMBDIF,K) = TAURAY(K,LAMBDIF) + TAUABS(K,LAMBDIF) ENDDO ENDDO C EXTRAPOLATE TAURAY(1) AND TAUABS(1) TO THE GROUND DO LAMBDIF = 1 , NDIFWAV Z1 = 0.5d0*Z(2) Z2 = 0.5d0*( Z(3) + Z(2) ) RAYV(LAMBDIF) = (Z2*TAURAY(1,LAMBDIF) - & Z1*TAURAY(2,LAMBDIF))/(Z2 - Z1) ABSV(LAMBDIF) = (Z2*TAUABS(1,LAMBDIF) - & Z1*TAUABS(2,LAMBDIF))/(Z2 - Z1) TOTALV(LAMBDIF) = RAYV(LAMBDIF) + ABSV(LAMBDIF) ! TOTAL VERTICAL OPTICAL DEPTH ENDDO ! AT GROUND C SET UP THE MULTIPLE SCATTERING MATRIX DO LAMBDIF = 1 , NDIFWAV LAMBDA = LOWDIF + LAMBDIF - 1 TAUV = RAYV(LAMBDIF) TV = ABSV(LAMBDIF) DO K2 = 1 , NALT-1 TAUN = TAURAY(K2,LAMBDIF) TN = TAUABS(K2,LAMBDIF) DTAUN = 0.5D0*DTAURAY(K2,LAMBDIF) DO K1 = 1 , NALT-1 TAUK = TAURAY(K1,LAMBDIF) TK = TAUABS(K1,LAMBDIF) ARG1 = ABS( TAUK - TAUN + DTAUN ) + & ABS( TK - TN ) ARG2 = ABS( TAUK - TAUN - DTAUN ) + & ABS( TK - TN ) ARG3 = TAUV - TAUK + TV - TK ARG3 = DMAX1(ARG3,0.0D0) ! AVOID ROUND-OFF PROBLEMS ( ARG3 > 0 !! ) ARG4 = TAUV - TAUN + DTAUN + TV - TN ARG5 = TAUV - TAUN - DTAUN + TV - TN A(K1,K2,LAMBDIF) = -(0.5D0)*ABS( E(2,ARG1) - E(2,ARG2) ) - & ALBEDO(LAMBDA)*E(2,ARG3)* & ABS(E(3,ARG4) - E(3,ABS(ARG5))) IF( K2 .EQ. K1 ) THEN A(K1,K2,LAMBDIF) = E(2,ABS(DTAUN)) ! ACCOUNT FOR SELF SCATTERING ENDIF ENDDO ENDDO ENDDO ! END LOOP OVER LAMBDIF C NOW PERFORM LU DECOMPOSITON OF A. THIS ALGORITHM IS ADAPTED FROM "Numerical C Recipes in Fortran" WHERE IT IS LISTED AS 'LUDCMP'. THIS DECOMPOSITION C IS NEEDED IN ORDER TO SOLVE A SET OF MATRIX EQUATIONS IN SUBROUTINE DIURNL. C NOTICE THAT LAMBDIF IS A "DUMMY" INDEX IN THE FOLLOWING ALGORITHM. DO LAMBDIF = 1 , NDIFWAV DO I = 1 , NALT-1 AAMAX = 0.0D0 DO J = 1 , NALT-1 IF( ABS(A(J,I,LAMBDIF)) .GT. AAMAX ) THEN AAMAX = ABS(A(J,I,LAMBDIF)) ENDIF ENDDO SCALE(I) = 1.0D0/AAMAX ENDDO DO J = 1 , NALT-1 DO I = 1 , J-1 SUM = A(J,I,LAMBDIF) DO K = 1 , I-1 SUM = SUM - A(K,I,LAMBDIF)*A(J,K,LAMBDIF) ENDDO A(J,I,LAMBDIF) = SUM ENDDO AAMAX = 0.0D0 DO I = J , NALT-1 SUM = A(J,I,LAMBDIF) DO K = 1 , J-1 SUM = SUM - A(K,I,LAMBDIF)*A(J,K,LAMBDIF) ENDDO A(J,I,LAMBDIF) = SUM DUM = SCALE(I)*ABS(SUM) IF( DUM .GE. AAMAX) THEN IMAX = I AAMAX = DUM ENDIF ENDDO IF( J .NE. IMAX ) THEN DO K = 1 , NALT-1 DUM = A(K,IMAX,LAMBDIF) A(K,IMAX,LAMBDIF) = A(K,J,LAMBDIF) A(K,J,LAMBDIF) = DUM ENDDO SCALE(IMAX) = SCALE(J) ENDIF INDX(J,LAMBDIF) = IMAX IF( A(J,J,LAMBDIF) .EQ. 0.0d0 ) A(J,J,LAMBDIF) = SMALL IF( J .NE. NALT-1 ) THEN DUM = 1.0D0/A(J,J,LAMBDIF) DO I = J+1 , NALT-1 A(J,I,LAMBDIF) = DUM*A(J,I,LAMBDIF) ENDDO ENDIF ENDDO ENDDO ! END LOOP OVER DUMMY INDICES RETURN END C---------------------------------------------------------------------------- C---------------------------------------------------------------------------- C---------------------------------------------------------------------------- SUBROUTINE ALBAVG INCLUDE 'PARAMETERS.INC' COMMON & /BLOCK16/ COSSZA(NTRAP) & /BLOCK30/ NPI2 & /ONED51/ ALBEDO(NWAVE) & /ONED52/ TAUVERT(NDIFWAV,NALT) , TOTALV(NDIFWAV) & /BLOCK41/ S(NTRAP,NWAVE,NALT) & /ONED58/ S0(NTRAP,NDIFWAV,NALT) DIMENSION TGROUND(NTRAP,NDIFWAV) C THIS SUBROUTINE INCLUDES GROUND REFLECTION IN THE SOURCE TERM C COMPUTE FLUX DUE TO DIRECT ATTENUATION, GROUND REFLECTANCE DO LAMBDIF = 1 , NDIFWAV LAMBDA = LOWDIF + LAMBDIF - 1 A = 2.0D0*ALBEDO(LAMBDA) DO NANG = 1 , NPI2 B = DMAX1( COSSZA(NANG) , 0.0D0 ) TGROUND(NANG,LAMBDIF) = A*B*S(NANG,LAMBDA,1) ENDDO ENDDO C THE MULTIPLE SCATTERING MATRIX REQUIRES THAT THE SOURCE FUNCTION - S0 BE C EVALUATED AT THE MIDPOINT OF THE LAYERS. THIS ACCOUNT FOR THE C S0 = 0.5*(TRANS(K) + TRANS(K+1)) TERM BELOW DO K = 1 , NALT-1 DO LAMBDIF = 1 , NDIFWAV LAMBDA = LOWDIF + LAMBDIF - 1 TVERT = TOTALV(LAMBDIF) - TAUVERT(LAMBDIF,K) IF( ALBEDO(LAMBDA) .GT. 0.001D0 )THEN SGROUND = E(2,TVERT) ELSE SGROUND = 0.0D0 ENDIF C ADD REFLECTION FROM GROUND ONTO SOURCE FUNCTION DO NANG = 1 , NPI2 S0(NANG,LAMBDIF,K) = 0.5D0*(S(NANG,LAMBDA,K) + & S(NANG,LAMBDA,K+1)) S0(NANG,LAMBDIF,K) = S0(NANG,LAMBDIF,K) + & SGROUND*TGROUND(NANG,LAMBDIF) ENDDO ENDDO ENDDO RETURN END C---------------------------------------------------------------------------- C---------------------------------------------------------------------------- C---------------------------------------------------------------------------- SUBROUTINE ALBDIU INCLUDE 'PARAMETERS.INC' COMMON & /BLOCK19/ TRANS(NWAVE,NALT) & /ONED50/ COSSZA & /ONED51/ ALBEDO(NWAVE) & /ONED52/ TAUVERT(NDIFWAV,NALT) , TOTALV(NDIFWAV) & /ONED57/ S0(NDIFWAV,NALT) DIMENSION TGROUND(NDIFWAV) C THIS SUBROUTINE INCLUDES GROUND REFLECTION IN THE SOURCE TERM C COMPUTE FLUX DUE TO DIRECT ATTENUATION, GROUND REFLECTANCE DO LAMBDIF = 1 , NDIFWAV LAMBDA = LOWDIF + LAMBDIF - 1 A = 2.0D0*ALBEDO(LAMBDA) B = DMAX1( COSSZA , 0.0D0 ) TGROUND(LAMBDIF) = A*B*TRANS(LAMBDA,1) ENDDO C THE MULTIPLE SCATTERING MATRIX REQUIRES THAT THE SOURCE FUNCTION - S0 BE C EVALUATED AT THE MIDPOINT OF THE LAYERS. THIS ACCOUNT FOR THE C S0 = 0.5*(TRANS(K) + TRANS(K+1)) TERM BELOW. TOTALV IS THE VERTICAL C OPTICAL DEPTH AT THE GROUND. DO K = 1 , NALT-1 DO LAMBDIF = 1 , NDIFWAV LAMBDA = LOWDIF + LAMBDIF - 1 S0(LAMBDIF,K) = 0.5D0*(TRANS(LAMBDA,K+1) + TRANS(LAMBDA,K)) TVERT = TOTALV(LAMBDIF) - TAUVERT(LAMBDIF,K) IF( ALBEDO(LAMBDA) .GT. 0.0D0 )THEN SGROUND = E(2,TVERT) ELSE SGROUND = 0.0D0 ENDIF C ADD REFLECTION FROM GROUND ONTO SOURCE FUNCTION S0(LAMBDIF,K) = S0(LAMBDIF,K) + SGROUND*TGROUND(LAMBDIF) ENDDO ENDDO RETURN END C---------------------------------------------------------------------------- C---------------------------------------------------------------------------- C---------------------------------------------------------------------------- SUBROUTINE MSAVG INCLUDE 'PARAMETERS.INC' COMMON & /BLOCK30/ NPI2 & /BLOCK35/ INDX(NALT,NDIFWAV) , & A(NALT,NALT,NDIFWAV) & /ONED58/ S0(NTRAP,NDIFWAV,NALT) C INCLUDE MULTIPLE SCATTERING EFFECTS: C THIS BLOCK SOLVES A MATRIX EQUATION A*S = S0 , WHERE A IS THE MULTIPLE C SCATTERING MATRIX FROM SUBROUTINE DIFFUSE AND S0 IS THE SOURCE FUNCTION C S0 IS OVERWRITTEN AND THE SOLUTION IS RETURNED AS S0. C RECALL THAT THE SOURCE FUNCTION "S0" AND THE MULTIPLE SCATTERING C MATRIX "A" ARE EVALUATED AT THE MIDPOINT OF THE LAYERS. THUS, AFTER C SOLVING FOR "S", WE WILL HAVE TO INTERPOLATE THE VALUES ONTO THE GRID C DEFINED BY Z(K). NEND = NPI2 DO LAMBDIF = 1 , NDIFWAV II = 0 DO K1 = 1 , NALT-1 KL = INDX(K1,LAMBDIF) DO NANG = 1 , NEND SUM = S0(NANG,LAMBDIF,KL) S0(NANG,LAMBDIF,KL) = S0(NANG,LAMBDIF,K1) IF( II .NE. 0 ) THEN DO K2 = II , K1-1 SUM = SUM - A(K1,K2,LAMBDIF)*S0(NANG,LAMBDIF,K2) ENDDO ELSE IF (SUM .NE. 0 ) THEN II = K1 ENDIF S0(NANG,LAMBDIF,K1) = SUM ENDDO ENDDO DO K1 = NALT-1 , 1 , -1 DO NANG = 1 , NEND SUM = S0(NANG,LAMBDIF,K1) DO K2 = K1+1 , NALT-1 SUM = SUM - A(K1,K2,LAMBDIF)*S0(NANG,LAMBDIF,K2) ENDDO S0(NANG,LAMBDIF,K1) = SUM/A(K1,K1,LAMBDIF) ENDDO ENDDO ENDDO ! END LOOP OVER LAMBDIF CALL INTERAVG ! INTERPOLATE ONTO CORRECT GRID RETURN END C---------------------------------------------------------------------------- C---------------------------------------------------------------------------- C---------------------------------------------------------------------------- SUBROUTINE MSDIU INCLUDE 'PARAMETERS.INC' COMMON & /BLOCK35/ INDX(NALT,NDIFWAV) , & A(NALT,NALT,NDIFWAV) & /ONED57/ S0(NDIFWAV,NALT) C INCLUDE MULTIPLE SCATTERING EFFECTS: C THIS BLOCK SOLVES A MATRIX EQUATION A*S = S0 , WHERE A IS THE MULTIPLE C SCATTERING MATRIX FROM SUBROUTINE DIFFUSE AND S0 IS THE SOURCE FUNCTION. C S0 IS OVERWRITTEN AND THE SOLUTION IS RETURNED AS S0. C RECALL THAT THE SOURCE FUNCTION "S0" AND THE MULTIPLE SCATTERING C MATRIX "A" ARE EVALUATED AT THE MIDPOINT OF THE LAYERS. THUS, AFTER C SOLVING FOR "S", WE WILL HAVE TO INTERPOLATE THE VALUES ONTO THE GRID C DEFINED BY Z(K). DO LAMBDIF = 1 , NDIFWAV II = 0 DO K1 = 1 , NALT-1 KL = INDX(K1,LAMBDIF) SUM = S0(LAMBDIF,KL) S0(LAMBDIF,KL) = S0(LAMBDIF,K1) IF( II .NE. 0 ) THEN DO K2 = II , K1-1 SUM = SUM - A(K1,K2,LAMBDIF)*S0(LAMBDIF,K2) ENDDO ELSE IF (SUM .NE. 0 ) THEN II = K1 ENDIF S0(LAMBDIF,K1) = SUM ENDDO DO K1 = NALT-1 , 1 , -1 SUM = S0(LAMBDIF,K1) DO K2 = K1+1 , NALT-1 SUM = SUM - A(K1,K2,LAMBDIF)*S0(LAMBDIF,K2) ENDDO S0(LAMBDIF,K1) = SUM/A(K1,K1,LAMBDIF) ENDDO ENDDO ! END LOOP OVER LAMBDIF CALL INTERDIU ! INTERPOLATE S0 ONTO CORRECT GRID RETURN END C---------------------------------------------------------------------------- C---------------------------------------------------------------------------- C---------------------------------------------------------------------------- SUBROUTINE INTERAVG INCLUDE 'PARAMETERS.INC' COMMON & /BLOCK30/ NPI2 & /BLOCK41/ S(NTRAP,NWAVE,NALT) & /ONED58/ S0(NTRAP,NDIFWAV,NALT) C THIS SUBROUTINE INTERPOLATES THE MULTIPLY SCATTERED SOURCE FUNCTION C BACK ONTO THE GRID USED IN CHEM1D. DO K = 2 , NALT-1 DO LAMBDIF = 1 , NDIFWAV LAMBDA = LOWDIF + LAMBDIF - 1 DO NANG = 1 , NPI2 S(NANG,LAMBDA,K) = 0.5D0*( S0(NANG,LAMBDIF,K) + & S0(NANG,LAMBDIF,K-1) ) ENDDO ENDDO ENDDO C SET UPPER AND LOWER BOUNDRIES TO NEAREST VALUE DO LAMBDIF = 1 , NDIFWAV LAMBDA = LOWDIF + LAMBDIF - 1 DO NANG = 1 , NPI2 S(NANG,LAMBDA,1) = S0(NANG,LAMBDIF,1) S(NANG,LAMBDA,NALT) = S0(NANG,LAMBDIF,NALT-1) ENDDO ENDDO RETURN END C---------------------------------------------------------------------------- C---------------------------------------------------------------------------- C---------------------------------------------------------------------------- SUBROUTINE INTERDIU INCLUDE 'PARAMETERS.INC' COMMON & /BLOCK19/ TRANS(NWAVE,NALT) & /ONED57/ S0(NDIFWAV,NALT) C THIS SUBROUTINE INTERPOLATES THE MULTIPLY SCATTERED SOURCE FUNCTION C BACK ONTO THE GRID USED IN CHEM1D. DO K = 2 , NALT-1 DO LAMBDIF = 1 , NDIFWAV LAMBDA = LOWDIF + LAMBDIF - 1 TRANS(LAMBDA,K) = 0.5D0*( S0(LAMBDIF,K) + & S0(LAMBDIF,K-1) ) ENDDO ENDDO C SET UPPER AND LOWER BOUNDRIES TO NEAREST VALUE DO LAMBDIF = 1 , NDIFWAV LAMBDA = LOWDIF + LAMBDIF - 1 TRANS(LAMBDA,1) = S0(LAMBDIF,1) TRANS(LAMBDA,NALT) = S0(LAMBDIF,NALT-1) ENDDO RETURN END