      SUBROUTINE DELTAE(ROS,G,WC1,TO1,RE1,TR1,RMUE,RE2,TR2,IL1,IL2,LON) 
  
C........  DELTA EDDINGTON APPROXIMATION  (VECTORIZED VERSION)
C 
C  INPUT
C     ROS: REFLECTIVITY OF UNDERLYING LAYER 
C     G  : ASSYMETRY FACTOR 
C     TO1 : OPTICAL THICKNESS 
C     GP,WCP,TOP: EQUIVALENT OPTICAL PARAMETERS 
C  OUTPUT 
C     RE1 : REFLECT. OF LAYER WITHOUT REF. FROM UNDERLYING LAYER
C     TR1 : TRANS.       "              "                 " 
C     RE2 : REFLECT.     "      WITH     "               "
C     TR2 : TRANS.       "       "        "             " 
  
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)

      COMMON /RADSIZE/ KMX, LMX, KMXP, NGLP1, NGL 
  
      DIMENSION ROS(LMX),G(LMX),WC1(LMX),TO1(LMX),RE1(LMX),TR1(LMX),
     1          RMUE(LMX),RE2(LMX),TR2(LMX) 
  
  
C---------------------------------------------------------------------- 
  
      DO 1 I=IL1,IL2
        IF (TO1(I) .EQ. 0.) THEN
          RE1(I) = 0. 
          TR1(I) = 1. 
          RE2(I) = 0. 
          TR2(I) = 1. 
        ELSE
          FF=G(I)*G(I)
          GP=G(I)/(1.+G(I)) 
          TOP=(1.-WC1(I)*FF)*TO1(I) 
          WCP=(1-FF)*WC1(I)/(1.-WC1(I)*FF)
          DT=2./3.
          X1=1.-WCP*GP
          WM=1.-WCP 
          RM2=RMUE(I)*RMUE(I) 
          RK=SQRT(3.*WM*X1) 
          X2=4.*(1.-RK*RK*RM2)
          RP=SQRT(3.*WM/X1) 
          ALPHA=3.*WCP*RM2*(1.+GP*WM)/X2
          BETA=3.*WCP*RMUE(I)*(1.+3.*GP*RM2*WM)/X2
          EXMUO=EXP(-TOP/RMUE(I)) 
          EXKP=EXP(RK*TOP)
          EXKM=1./EXKP
          XP2P=1.+DT*RP 
          XM2P=1.-DT*RP 
          AP2B=ALPHA+DT*BETA
          AM2B=ALPHA-DT*BETA
  
C........  WITHOUT REFLEXION FROM THE UNDERLYING LAYER
  
          A11=XP2P
          A12=XM2P
          A13=AP2B
          A22=XP2P*EXKP 
          A21=XM2P*EXKM 
          A23=AM2B*EXMUO
          D=A11*A22-A21*A12 
          C1=(A22*A13-A12*A23)/D
          C2=(A11*A23-A21*A13)/D
          RI0=C1+C2-ALPHA 
          RI1=RP*(C1-C2)-BETA 
          RE1(I)=(RI0-DT*RI1)/RMUE(I) 
          RI0D=C1*EXKM+C2*EXKP-ALPHA*EXMUO
          RI1D=RP*(C1*EXKM-C2*EXKP)-BETA*EXMUO
          TR1(I)=EXMUO+(RI0D+DT*RI1D)/RMUE(I) 
  
C........  WITH REFLEXION FROM THE UNDERLYING LAYER 
  
          A21=A21-ROS(I)*XP2P*EXKM
          A22=A22-ROS(I)*XM2P*EXKP
          A23=A23-ROS(I)*EXMUO*(AP2B-RMUE(I)) 
          D=A11*A22-A21*A12 
          C1=(A22*A13-A12*A23)/D
          C2=(A11*A23-A21*A13)/D
          RI0=C1+C2-ALPHA 
          RI1=RP*(C1-C2)-BETA 
          RE2(I)=(RI0-DT*RI1)/RMUE(I) 
          RI0D=C1*EXKM+C2*EXKP-ALPHA*EXMUO
          RI1D=RP*(C1*EXKM-C2*EXKP)-BETA*EXMUO
          TR2(I)=EXMUO+(RI0D+DT*RI1D)/RMUE(I) 
        ENDIF 
    1 CONTINUE
  
C---------------------------------------------------------------------- 
  
      RETURN
      END 
