      SUBROUTINE GFLUXV(NTODO,WAVEN,DTDEL,TDEL,WDEL,CDEL
     &          , F0PI,RSF,BTOP,BSURF,FP,FM,FMIDP,FMIDM,IPRINT)
c      PARAMETER (NL=101)
        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,UBARI,UBARV,UBAR0,CMMID,TAUMID,RSF,BTOP
      REAL EM,EP,PI,WAVEN,CPMID,BSURF,G4
      COMMON /UBARED/ UBARI,UBARV,UBAR0
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
