SUBROUTINE GFLUXV(NTODO,WAVEN,DTDEL,TDEL,WDEL,CDEL & , F0PI,RSF,BTOP,BSURF,FP,FM,FMIDP,FMIDM,IPRINT) c PARAMETER (NL=101) USE TGMDAT_MOD, ONLY: UBARI,UBARV,UBAR0 IMPLICIT NONE INTEGER IPRINT,I,J,K,IDELTA,NAYER,NL,NTODO PARAMETER (NL=401) C ??FLAG? THIS VALUE (101) MUST BE .GE. NTODO C* THIS SUBROUTINE TAKES THE OPTICAL CONSTANTS AND BOUNDARY CONDITONS C FOR THE VISIBLE FLUX AT ONE WAVELENGTH AND SOLVES FOR THE FLUXES AT C THE LEVELS. THIS VERSION IS SET UP TO WORK WITH LAYER OPTICAL DEPTHS C MEASURED FROM THE TOP OF EACH LAEYER. (DTAU) TOP OF EACH LAYER HAS C OPTICAL DEPTH TAU(N).IN THIS SUB LEVEL N IS ABOVE LAYER N. THAT IS LAYER N C HAS LEVEL N ON TOP AND LEVEL N+1 ON BOTTOM. OPTICAL DEPTH INCREASES C FROM TOP TO BOTTOM. SEE C.P. MCKAY, TGM NOTES. C THIS SUBROUTINE DIFFERS FROM ITS IR CONTERPART IN THAT HERE WE SOLVE FOR C THE FLUXES DIRECTLY USING THE GENERALIZED NOTATION OF MEADOR AND WEAVOR C J.A.S., 37, 630-642, 1980. REAL B81,B82,R81 REAL AP,AM,F0PI,DENOM,CMMID,TAUMID,RSF,BTOP REAL EM,EP,PI,WAVEN,CPMID,BSURF,G4 C THIS NEXT ROW OF VARIABLES ARE THOSE ACTUALLY USED IN THE C ROUTINE REAL W0(NL), COSBAR(NL),DTAU(NL),TAU(NL) REAL wdel(ntodo-1),cdel(ntodo-1),dtdel(ntodo-1) REAL tdel(ntodo),fm(ntodo),fp(ntodo) REAL fmidm(ntodo),fmidp(ntodo) REAL ALPHA(NL),LAMDA(NL),XK1(NL),XK2(NL) REAL G1(NL),G2(NL),G3(NL) REAL GAMA(NL),CP(NL),CM(NL),CPM1(NL),CMM1(NL) &,E1(NL),E2(NL),E3(NL),E4(NL),EXPTRM(NL) DATA PI/3.14159265358979323846/ DATA IDELTA/1/ NAYER=NTODO-1 505 CONTINUE C TURN ON THE DELTA-FUNCTION IF REQUIRED HERE c PRINT*,' UBAR0 ',UBARI,UBARV,UBAR0 IF (IDELTA .EQ. 0) THEN DO 101 J=1,NAYER W0(J)=WDEL(J) COSBAR(J)=CDEL(J) DTAU(J)=DTDEL(J) TAU(J)=TDEL(J) 101 CONTINUE TAU(NTODO)=TDEL(NTODO) ELSE C FOR THE DELTA FUNCTION HERE... TAU(1)=TDEL(1)*(1.-WDEL(1)*CDEL(1)**2) DO 102 J=1,NAYER W0(J)=WDEL(J)*(1.-CDEL(J)**2)/(1.-WDEL(J)*CDEL(J)**2) COSBAR(J)=CDEL(J)/(1.+CDEL(J)) DTAU(J)=DTDEL(J)*(1.-WDEL(J)*CDEL(J)**2) TAU(J+1)=TAU(J)+DTAU(J) c print*,'W0=',W0(J),' COSBAR=',COSBAR(J) c print*,'DTAU=',DTAU(J),' TAU=',TAU(J) c print*,'WDEL=',WDEL(J),' CDEL=',CDEL(J),' DTDEL=',DTDEL(J) 102 CONTINUE ENDIF C WE GO WITH THE HEMISPHERIC CONSTANT APPROACH C AS DEFINED BY M&W - THIS IS THE WAY THE IR IS DONE DO 20 J=1,NAYER ALPHA(J)=SQRT( (1.-W0(J))/(1.-W0(J)*COSBAR(J)) ) C THIS SET OF G'S IS FOR THE TWO STREAM HEMI-CONSTANT C G1(J)=(1.-W0(J)*0.5*(1.+COSBAR(J)))/UBARV C G2(J)=W0(J)*0.5*(1.-COSBAR(J))/UBARV C G3(J)=0.5*(1.-COSBAR(J)) C SET OF CONSTANTS DETERMINED BY DOM G1(J)= (SQRT(3.)*0.5)*(2. - W0(J)*(1.+COSBAR(J))) G2(J)= (SQRT(3.)*W0(J)*0.5)*(1.-COSBAR(J)) G3(J)=0.5*(1.-SQRT(3.)*COSBAR(J)*UBAR0) LAMDA(J)=SQRT(G1(J)**2 - G2(J)**2) GAMA(J)=(G1(J)-LAMDA(J))/G2(J) 20 CONTINUE DO 7 J=1,NAYER G4=1.-G3(J) DENOM=LAMDA(J)**2 - 1./UBAR0**2 c print*,'G1(J)**2- G2(J)**2',G1(J)**2 - G2(J)**2 c print*,'SQRT(G1(J)**2 - G2(J)**2)',SQRT(G1(J)**2 c & - G2(J)**2) c print*,'DENOM=',LAMDA(J)**2,1./UBAR0**2 C THERE IS A POTENTIAL PROBLEM HERE IF W0=0 AND UBARV=UBAR0 C THEN DENOM WILL VANISH. THIS ONLY HAPPENS PHYSICALLY WHEN C THE SCATTERING GOES TO ZERO C PREVENT THIS WITH A IF STATEMENT IF ( DENOM .EQ. 0.) THEN PRINT*,UBAR0 DENOM=1.E-4 PRINT*,'G1(J)**2- G2(J)**2',G1(J)**2 - G2(J)**2 PRINT*,'SQRT(G1(J)**2 - G2(J)**2)',SQRT(G1(J)**2 & - G2(J)**2) PRINT*,'DENOM=',LAMDA(J)**2,1./UBAR0**2 WRITE (6,99) 99 FORMAT (' DENOM ZERO; RESET # 1IN GFLUXV, W0=0?') END IF AM=F0PI*W0(J)*(G4 *(G1(J)+1./UBAR0) +G2(J)*G3(J) )/DENOM AP=F0PI*W0(J)*(G3(J)*(G1(J)-1./UBAR0) +G2(J)*G4 )/DENOM C* CPM1 AND CMM1 ARE THE CPLUS AND CMINUS TERMS EVALUATED C AT THE TOP OF THE LAYER, THAT IS LOWER OPTICAL DEPTH TAU(J) CPM1(J)=AP*EXP(-TAU(J)/UBAR0) CMM1(J)=AM*EXP(-TAU(J)/UBAR0) C* CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE C BOTTOM OF THE LAYER. THAT IS AT HIGHER OPTICAL DEPTH TAU(J+1) CP(J)=AP*EXP(-TAU(J+1)/UBAR0) CM(J)=AM*EXP(-TAU(J+1)/UBAR0) c print*,'AM=',AM,' AP=',AP c print*,'CPM1(J)=',CPM1(J),' CMM1(J)=',CMM1(J) c print*,'CP(J)=',CP(J),' CM(J)=',CM(J) 7 CONTINUE C* C* NOW CALCULATE THE EXPONENTIAL TERMS NEEDED C* FOR THE TRIDIAGONAL ROTATED LAYERED METHOD C* WARNING IF DTAU(J) IS GREATER THAN ABOUT 35 C* WE CLIPP IT TO AVOID OVERFLOW. C* EXP (TAU) - EXP(-TAU) WILL BE NONSENSE THIS IS C* CORRECTED IN THE DSOLVER ROUTINE. ??FLAG? DO 103 J=1,NAYER EXPTRM(J)=35. IF ( LAMDA(J)*DTAU(J) .LT. 35. ) EXPTRM(J)=LAMDA(J)*DTAU(J) 103 CONTINUE C* NEED TO CLIPP THE EXPONENTIAL HERE. DO 8 J=1,NAYER EP=EXP(EXPTRM(J)) EM=1.0/EP E1(J)=EP+GAMA(J)*EM E2(J)=EP-GAMA(J)*EM E3(J)=GAMA(J)*EP+EM E4(J)=GAMA(J)*EP-EM 8 CONTINUE C B81=BTOP B82=BSURF R81=RSF CALL DSOLVER(NAYER,GAMA,CP,CM,CPM1,CMM1 & ,E1,E2,E3,E4,B81,B82,R81,XK1,XK2) C C EVALUATE THE NTODO FLUXES THROUGH THE NAYER LAYERS C USE THE TOP (TAU=0) OPTICAL DEPTH EXPRESSIONS TO EVALUATE FP AND FM C AT THE THE TOP OF EACH LAYER,J = LEVEL J DO 46 J=1,NAYER FP(J)= XK1(J) + GAMA(J)*XK2(J) + CPM1(J) FM(J)=GAMA(J)*XK1(J) + XK2(J) + CMM1(J) 46 CONTINUE C USE EXPRESSION FOR BOTTOM FLUX TO GET THE FP AND FM AT NTODO J=NAYER EP=EXP(EXPTRM(J)) EM=1.0/EP FP(J+1)=XK1(J)*EP + GAMA(J)*XK2(J)*EM + CP(J) FM(J+1)=XK1(J)*EP*GAMA(J) + XK2(J)*EM + CM(J) C NOTE THAT WE HAVE SOLVED FOR THE FLUXES DIRECTLY AND NO C FURTHER INTEGRATION IS NEEDED, THE UBARV TERM IS ABSORBED C INTO THE DEFINITION OF G1 THU G4 C UBARV = .5 IS HEMISPHERIC CONSTANT C UBARV = SQRT(3) IS GAUSS QUADRATURE C AND OTHER CASES AS PER MEADOR AND WEAVOR, JAS, 37, 630-643,1980. C C ADD THE DIRECT FLUX TERM TO THE DOWNWELLING RADIATION, LIOU 182 DO 80 J=1,NTODO FM(J)=FM(J)+UBAR0*F0PI*EXP(-TAU(J)/UBAR0) 80 CONTINUE C C NOW WE CALCULATE THE FLUXES AT THE MIDPOINTS OF THE LAYERS. C EXACLTY ANALOGOUS TO THE ABOVE COMPUTATION C DO 1982 J=1,NAYER EP=EXP(0.5*EXPTRM(J)) EM=1.0/EP G4=1.-G3(J) DENOM=LAMDA(J)**2 - 1./UBAR0**2 C THERE IS A POTENTIAL PROBLEM HERE IF W0=0 AND UBARV=UBAR0 C THEN DENOM WILL VANISH. THIS ONLY HAPPENS PHYSICALLY WHEN C THE SCATTERING GOES TO ZERO C PREVENT THIS WITH A IF STATEMENT IF ( DENOM .EQ. 0.) THEN DENOM=1.E-4 PRINT*,'G1(J)**2- G2(J)**2',G1(J)**2 - G2(J)**2 PRINT*,'SQRT(G1(J)**2 - G2(J)**2)',SQRT(G1(J)**2 & - G2(J)**2) PRINT*,'DENOM=',LAMDA(J)**2,1./UBAR0**2 WRITE (6,78) 78 FORMAT (' DENOM ZERO; RESET # 2 IN GFLUXV, W0=0?') END IF AM=F0PI*W0(J)*(G4 *(G1(J)+1./UBAR0) +G2(J)*G3(J) )/DENOM AP=F0PI*W0(J)*(G3(J)*(G1(J)-1./UBAR0) +G2(J)*G4 )/DENOM C* CPMID AND CMMID ARE THE CPLUS AND CMINUS TERMS EVALUATED C AT THE MIDDLE OF THE LAYER, THAT IS LOWER OPTICAL DEPTH TAU(J)+ TAUMID= (TAU(J)+0.5*DTAU(J) ) CPMID=AP*EXP(-TAUMID/UBAR0) CMMID=AM*EXP(-TAUMID/UBAR0) FMIDP(J)=XK1(J)*EP + GAMA(J)*XK2(J)*EM + CPMID FMIDM(J)=XK1(J)*EP*GAMA(J) + XK2(J)*EM + CMMID C ADD THE DIRECT FLUX TO THE DOWNWELLING TERM FMIDM(J)= FMIDM(J) +UBAR0*F0PI*EXP(-TAUMID/UBAR0) 1982 CONTINUE C ** NOW PRINTOUT IF NECESSARY IF (IPRINT .LT. 9) RETURN WRITE(6,601) F0PI,RSF,BTOP,BSURF DO 120 J=1,NAYER WRITE (6,301) TAU(J),FP(J),FM(J),DTAU(J),W0(J),COSBAR(J) & ,ALPHA(J), LAMDA(J),G1(J),G2(J),G3(J) 120 CONTINUE J=NTODO WRITE (6,301) TAU(J),FP(J),FM(J) WRITE (6,602) DO 130 J=1,NAYER WRITE (6,301) CP(J),CM(J),E1(J),E2(J),E3(J),E4(J),CPM1(J) & ,CMM1(J),GAMA(J) 130 CONTINUE 301 FORMAT(1X,1P13E10.3) 602 FORMAT(' CP(J),CM(J),E1(J),E2(J),EM1(J),EM2(J),CPM1(J),CMM1(J)' & ,',GAMA(J)') 601 FORMAT(1X,'F0PI,RSF,BTOP,BSURF= ',1P4E10.3,/ & ' TAU,FUP,FDOWN,DTAU(J),W0(J),COSBAR(J)' & ,',ALPHA(J),LAMDA(J),G1,G2,G3') C *******************************************************************012 RETURN END