cc--- boundary layer calculation

      parameter (NP1=124,NP2=34)
      
      DIMENSION DENS(NP1,NP2), XMOM(NP1,NP2), YMOM(NP1,NP2),
     *     ENER(NP1,NP2),
     *     ADENS(NP1),AXMOM(NP1),AYMOM(NP1),AENER(NP1),
     *     ASDENS(NP1),ASXMOM(NP1),ASYMOM(NP1),ASENER(NP1),
     *     FM(NP1),FP(NP1),FPX(NP1),FE(NP1),AAD(4),AAX(4),AAY(4),
     *     AAE(4),AAAD(NP1),AAAX(NP1),AAAY(NP1),AAAE(NP1),
     *     XX(NP1),YY(NP2),DYY(NP2), DXX(NP1)


      PI=4.0*atan(1.0)
      
      IXEND=NP1-2
      IYEND=NP2-2
      
      CK=3.0
      CR2 = 1.414213562
      PR  = 1.772453851
      GAM = 1.4
      RGAMH = SQRT(GAM/2.)
      PRM=(CK+1)/2.0
      PROT1  = 1.0E-10
      EPS=0.3
      AA=1.0
      EP=1.0E-7


      OPEN(UNIT=1,FILE='data1.out',STATUS='unknown')
      OPEN(UNIT=2,FILE='data2.out',STATUS='unknown')
      OPEN(UNIT=3,FILE='data3.out',STATUS='unknown')
      OPEN(UNIT=4,FILE='data4.out',STATUS='unknown')
      OPEN(UNIT=5,FILE='da.out',STATUS='old')
      OPEN(UNIT=7,FILE='data.out',STATUS='unknown')
      OPEN(UNIT=11,FILE='xy.out',STATUS='unknown')

      TIME   = 0.0
      TAUMAX = 10.0
      UBD    = 0.0
      IT     = 1



      NIN=70
      IXS=40
c--- viscosity coefficient to get Re= 1.0E5      
      VIS=0.0026395*1.18980/3.0

      TSTOP=400.0

      XMACH=0.30
      UU=1.0
      VV=0.0
      PPP=1.0/(GAM*XMACH*XMACH)
      EE=0.5*UU*UU+2.5*PPP
      CCC=SQRT(GAM*PPP)
      
      XS=0.0

      ETAXL=1.10
      ETAXR=1.05

      DXX(IXS+1)=0.100
      XX(IXS+1)=0.5*DXX(IXS+1)

      DO IX = IXS+2, IXEND+1
         DXX(IX)=ETAXR*DXX(IX-1)
         XX(IX)=XX(IX-1)+0.5*(DXX(IX)+DXX(IX-1))
      ENDDO

      DO IX=IXS, 2, -1
         DXX(IX)=ETAXL*DXX(IX+1)
         XX(IX) = XX(IX+1)-0.5*(DXX(IX)+DXX(IX+1)) 
      ENDDO

      DO IX=2, IXEND
         WRITE(*,*) 'DXX(IX), XX(IX)', IX, DXX(IX), XX(IX)
         
      ENDDO

      
      ETAY=1.2
      DYY(1)=0.0175
      YY(1)=-(ETAY*DYY(1)+0.5*DYY(1))

      DO IY=2,IYEND+1
         DYY(IY)=ETAY*DYY(IY-1)
         YY(IY)=YY(IY-1)+0.5*(DYY(IY)+DYY(IY-1))
         write(*,*) 'DYY(IY), YY(IY)', IY, DYY(IY), YY(IY)
      ENDDO


C     INITIAL CONDITIONS:
      DO 102 IY=3,IYEND
         DO 101 IX = 3, IXEND
c            DENS(IX,IY) = 1.0
c            XMOM(IX,IY) = UU
c            YMOM(IX,IY) = 0.0
c            ENER(IX,IY) = EE

      read(5,81) DENS(IX,IY),
     *           XMOM(IX,IY),YMOM(IX,IY),ENER(IX,IY)
 81   FORMAT(E15.7,1X,E15.7,1X,E15.7,1X,E15.7)

 101     CONTINUE
 102  CONTINUE
C BOUNDARY CONDITIONS:
      
      
 601  CONTINUE
      TAU=100.0

      DO 201 IY = 3,IYEND
      DO 202 IX=3,IXEND
         U=SQRT(XMOM(IX,IY)*XMOM(IX,IY)+
     *        YMOM(IX,IY)*YMOM(IX,IY))/DENS(IX,IY)
         ELAM=0.25*(CK + 2.0)*DENS(IX,IY)/(ENER(IX,IY) 
     *        -0.5*XMOM(IX,IY)*XMOM(IX,IY)/DENS(IX,IY)
     *        -0.5 *YMOM(IX,IY)*YMOM(IX,IY)/DENS(IX,IY))
         IF(ELAM.LE.0.0) THEN 
            print *, 'ELAM NEGATIVE',DENS(IX,IY),XMOM(IX,IY),
     *           YMOM(IX,IY),ENER(IX,IY), IX,IY
         ENDIF
         S    = ABS(U) + RGAMH/SQRT(ELAM)
         UBD  = MAX(UBD,S)
         
         TAU1=DYY(IY)/UBD

         TAU=MIN(TAU1,TAU)

 202  CONTINUE
 201  CONTINUE

      TAU=(1.0-EPS)*TAU
      UBD = 0.0
      call flush(6)
      IF(MOD(IT,50).EQ.0) THEN
         print *, ' IT, TIME, TAU: ',IT,TIME,TAU
      ENDIF


      IF(TIME.GT.TSTOP) then


      DO IY=3,IYEND
      DO IX=3,IXEND 

      WRITE(7,10) DENS(IX,IY),XMOM(IX,IY),YMOM(IX,IY),ENER(IX,IY)
 10   FORMAT(E15.7,1X,E15.7,1X,E15.7,1X,E15.7)


      ENDDO
      ENDDO

      DO IY=3,IYEND
         DO IX=3,IXEND 
            WRITE(11,33) XX(IX), YY(IY), XMOM(IX,IY)/DENS(IX,IY)
 33         FORMAT(E15.7,1X,E15.7,1X,E15.7)
         ENDDO
      ENDDO


      IX=NIN

      DO IY=3,IYEND
      XL=XX(IX)-0.5*(XX(IXS+1)+XX(IXS))
      YV=YY(IY)-0.5*(YY(2)+YY(3))
      ETA=YV*SQRT(UU/(VIS*XL))
      UX=XMOM(IX,IY)/(UU*DENS(IX,IY))
      VY=(YMOM(IX,IY)/(DENS(IX,IY)*UU))*SQRT(UU*XL/VIS)
      WRITE(1,11) ETA, UX, VY
 11   FORMAT(E10.4,2X,E10.4,2X,E10.4)

      ENDDO

      IX=NIN+10
      DO IY=3,IYEND
      XL=XX(IX)-0.5*(XX(IXS+1)+XX(IXS))
      YV=YY(IY)-0.5*(YY(2)+YY(3))
      ETA=YV*SQRT(UU/(VIS*XL))
      UX=XMOM(IX,IY)/(UU*DENS(IX,IY))
      VY=(YMOM(IX,IY)/(DENS(IX,IY)*UU))*SQRT(UU*XL/VIS)


      WRITE(2,12) ETA, UX, VY
 12   FORMAT(E10.4,2X,E10.4,2X,E10.4)


      ENDDO

      IX=NIN+20
      DO IY=3,IYEND
      XL=XX(IX)-0.5*(XX(IXS+1)+XX(IXS))
      YV=YY(IY)-0.5*(YY(2)+YY(3))
      ETA=YV*SQRT(UU/(VIS*XL))
      UX=XMOM(IX,IY)/(UU*DENS(IX,IY))
      VY=(YMOM(IX,IY)/(DENS(IX,IY)*UU))*SQRT(UU*XL/VIS)
      WRITE(3,13) ETA, UX, VY
 13   FORMAT(E10.4,2X,E10.4,2X,E10.4)


      ENDDO

      IX=NIN+30
      DO IY=3,IYEND
      XL=XX(IX)-0.5*(XX(IXS+1)+XX(IXS))
      YV=YY(IY)-0.5*(YY(2)+YY(3))
      ETA=YV*SQRT(UU/(VIS*XL))
      UX=XMOM(IX,IY)/(UU*DENS(IX,IY))
      VY=(YMOM(IX,IY)/(DENS(IX,IY)*UU))*SQRT(UU*XL/VIS)
      WRITE(4,14) ETA, UX, VY
 14   FORMAT(E10.4,2X,E10.4,2X,E10.4)

      ENDDO

      RE=UU*FLOAT(XX(IXEND)-0.5*(XX(IXS+1)+XX(IXS)))/VIS
      write(*,*) 'Reynold= , Mach= , U= '
      write(*,*) RE, XMACH, UU

      STOP

      ENDIF


      TIME=TIME+TAU
      IT=IT+1


      DO 161 IY=3,IYEND
         DO 162 IX=3,IXEND
            ADENS(IX)=DENS(IX,IY)
            AXMOM(IX)=XMOM(IX,IY)
            AYMOM(IX)=YMOM(IX,IY)
            AENER(IX)=ENER(IX,IY)
 162     CONTINUE
         

         ADENS(1)= 1.0
         AXMOM(1)= UU
         AYMOM(1)=0.0
         AENER(1)=EE

         ADENS(2)=1.0
         AXMOM(2)=UU
         AYMOM(2)=0.0
         AENER(2)=EE

         D1=ADENS(3)
         X1=AXMOM(3)
         Y1=AYMOM(3)
         E1=AENER(3)

         PP1=0.4*(E1-0.5*(X1*X1+Y1*Y1)/D1)
         CC1=SQRT(GAM*PP1/D1)
         R2=X1/D1-2.0*CC1/(GAM-1)
         R1=UU+2.0*CCC/(GAM-1)
         UU1=0.5*(R1+R2)
         CS=0.25*(GAM-1.0)*(R1-R2)
         PP=CS*CS*D1/GAM

         ADENS(2)= 1.0
         AXMOM(2)= UU1
         AYMOM(2)= Y1
         AENER(2)= 0.5*(UU1*UU1+Y1*Y1)+ 2.5*PP

         ADENS(1)= 1.0
         AXMOM(1)= UU1
         AYMOM(1)= Y1
         AENER(1)= 0.5*(UU1*UU1+Y1*Y1)+ 2.5*PP


         D1=2.0*ADENS(IXEND)-ADENS(IXEND-1)
         X1=2.0*AXMOM(IXEND)-AXMOM(IXEND-1)
         Y1=2.0*AYMOM(IXEND)-AYMOM(IXEND-1)
         E1=2.0*AENER(IXEND)-AENER(IXEND-1)

         ADENS(IXEND+1)= D1
         AXMOM(IXEND+1)= X1
         AYMOM(IXEND+1)= Y1
         AENER(IXEND+1)= E1


         ADENS(IXEND+2)=2.0*ADENS(IXEND+1)-ADENS(IXEND)
         AXMOM(IXEND+2)=2.0*AXMOM(IXEND+1)-AXMOM(IXEND)
         AYMOM(IXEND+2)=2.0*AYMOM(IXEND+1)-AYMOM(IXEND)
         AENER(IXEND+2)=2.0*AENER(IXEND+1)-AENER(IXEND)


         DO II=2,IXEND+1

            DFD1=ADENS(II)-ADENS(II-1)
            DFD2=ADENS(II+1)-ADENS(II)

            DFX1=AXMOM(II)-AXMOM(II-1)
            DFX2=AXMOM(II+1)-AXMOM(II)

            DFY1=AYMOM(II)-AYMOM(II-1)
            DFY2=AYMOM(II+1)-AYMOM(II)

            DFE1=AENER(II)-AENER(II-1)
            DFE2=AENER(II+1)-AENER(II)
            
            AAAD(II)=(SIGN(AA,DFD1)+SIGN(AA,DFD2))
     *           *(ABS(DFD1)*ABS(DFD2))/(ABS(DFD1)+ABS(DFD2)+EP)
            AAAX(II)=(SIGN(AA,DFX1)+SIGN(AA,DFX2))
     *           *(ABS(DFX1)*ABS(DFX2))/(ABS(DFX1)+ABS(DFX2)+EP)
            AAAY(II)=(SIGN(AA,DFY1)+SIGN(AA,DFY2))
     *           *(ABS(DFY1)*ABS(DFY2))/(ABS(DFY1)+ABS(DFY2)+EP)
            AAAE(II)=(SIGN(AA,DFE1)+SIGN(AA,DFE2))
     *           *(ABS(DFE1)*ABS(DFE2))/(ABS(DFE1)+ABS(DFE2)+EP)

         ENDDO

         DO 999 I= 2,IXEND
            I1=I+1

            DX=XX(I1)-XX(I)

            AAD(1)=AAAD(I)/DX
            AAD(2)=ADENS(I)
            AAD(3)=ADENS(I1)
            AAD(4)=AAAD(I1)/DX

            AAX(1)=AAAX(I)/DX
            AAX(2)=AXMOM(I)
            AAX(3)=AXMOM(I1)
            AAX(4)=AAAX(I1)/DX

            AAY(1)=AAAY(I)/DX
            AAY(2)=AYMOM(I)
            AAY(3)=AYMOM(I1)
            AAY(4)=AAAY(I1)/DX
               
            AAE(1)=AAAE(I)/DX
            AAE(2)=AENER(I)
            AAE(3)=AENER(I1)
            AAE(4)=AAAE(I1)/DX

            TIE=TAU
            CALL DXEE(AAD,AAX,AAY,AAE,DX,TIE,VIS,AFM,AFP,AFX,AFE)
            FM(I)=AFM
            FP(I)=AFP
            FPX(I)=AFX
            FE(I)=AFE
               
 999     CONTINUE


         DO 502 IX=3,IXEND
            DX=0.5*(XX(IX+1)+XX(IX))
     *           -0.5*(XX(IX)+XX(IX-1))
            DENS(IX,IY) = ADENS(IX)+(FM(IX-1)-FM(IX))/DX
            XMOM(IX,IY) = AXMOM(IX)+(FP(IX-1)-FP(IX))/DX
            YMOM(IX,IY) = AYMOM(IX)+(FPX(IX-1)-FPX(IX))/DX
            ENER(IX,IY) = AENER(IX)+(FE(IX-1)-FE(IX))/DX
 502     CONTINUE
 161  CONTINUE

         
      DO 171 IX=3,IXEND
         DO 172 IY=3,IYEND
            ADENS(IY)=DENS(IX,IY)
            AYMOM(IY)=-XMOM(IX,IY)
            AXMOM(IY)=YMOM(IX,IY)
            AENER(IY)=ENER(IX,IY)
 172     CONTINUE

         D1=2.0*ADENS(IYEND)-ADENS(IYEND-1)
         X1=2.0*AXMOM(IYEND)-AXMOM(IYEND-1)
         Y1=2.0*AYMOM(IYEND)-AYMOM(IYEND-1)
         E1=2.0*AENER(IYEND)-AENER(IYEND-1)

         PP1=0.4*(E1-0.5*(X1*X1+Y1*Y1)/D1)
         CC1=SQRT(GAM*PP1/D1)
         R1=X1/D1+2.0*CC1/(GAM-1)
         R2=-2.0*CCC/(GAM-1)
         UU1=0.5*(R1+R2)
         CS=0.25*(GAM-1.0)*(R1-R2)
         PP=CS*CS*D1/GAM
 
         ADENS(IYEND+1)= D1
         AXMOM(IYEND+1)= D1*UU1
         AYMOM(IYEND+1)= Y1
         AENER(IYEND+1)= 0.5*D1*(UU1*UU1+Y1*Y1/(D1*D1))
     *        + 2.5*PP

         ADENS(IYEND+2)= 2.0*ADENS(IYEND+1)-ADENS(IYEND)
         AXMOM(IYEND+2)= 2.0*AXMOM(IYEND+1)-AXMOM(IYEND)
         AYMOM(IYEND+2)= 2.0*AYMOM(IYEND+1)-AYMOM(IYEND)
         AENER(IYEND+2)= 2.0*AENER(IYEND+1)-AENER(IYEND)
         
         IF(IX.GT.(IXS+1)) THEN
            ADENS(2)= ADENS(3)
            AXMOM(2)=-AXMOM(3)
            AYMOM(2)=-AYMOM(3)
            AENER(2)= AENER(3)

            ADENS(1)=ADENS(4)
            AXMOM(1)=-AXMOM(4)
            AYMOM(1)=-AYMOM(4)
            AENER(1)=AENER(4)
         ELSE
            ADENS(2)= ADENS(3)
            AXMOM(2)= -AXMOM(3)
            AYMOM(2)= AYMOM(3)
            AENER(2)= AENER(3)

            ADENS(1)=ADENS(4) 
            AXMOM(1)= -AXMOM(4)
            AYMOM(1)= AYMOM(4)
            AENER(1)= AENER(4)
         ENDIF



         DO II=2,IYEND+1
            DFD1=ADENS(II)-ADENS(II-1)
            DFD2=ADENS(II+1)-ADENS(II)

            DFX1=AXMOM(II)-AXMOM(II-1)
            DFX2=AXMOM(II+1)-AXMOM(II)

            DFY1=AYMOM(II)-AYMOM(II-1)
            DFY2=AYMOM(II+1)-AYMOM(II)

            DFE1=AENER(II)-AENER(II-1)
            DFE2=AENER(II+1)-AENER(II)
            
            AAAD(II)=(SIGN(AA,DFD1)+SIGN(AA,DFD2))
     *           *(ABS(DFD1)*ABS(DFD2))/(ABS(DFD1)+ABS(DFD2)+EP)
            AAAX(II)=(SIGN(AA,DFX1)+SIGN(AA,DFX2))
     *           *(ABS(DFX1)*ABS(DFX2))/(ABS(DFX1)+ABS(DFX2)+EP)
            AAAY(II)=(SIGN(AA,DFY1)+SIGN(AA,DFY2))
     *           *(ABS(DFY1)*ABS(DFY2))/(ABS(DFY1)+ABS(DFY2)+EP)
            AAAE(II)=(SIGN(AA,DFE1)+SIGN(AA,DFE2))
     *           *(ABS(DFE1)*ABS(DFE2))/(ABS(DFE1)+ABS(DFE2)+EP)
               
         ENDDO

 979     DO 998  I =  2 , IYEND

            I1=I+1

            DY=YY(I1)-YY(I)

            AAD(1)=AAAD(I)/DY
            AAD(2)=ADENS(I)
            AAD(3)=ADENS(I1)
            AAD(4)=AAAD(I1)/DY

            AAX(1)=AAAX(I)/DY
            AAX(2)=AXMOM(I)
            AAX(3)=AXMOM(I1)
            AAX(4)=AAAX(I1)/DY

            AAY(1)=AAAY(I)/DY
            AAY(2)=AYMOM(I)
            AAY(3)=AYMOM(I1)
            AAY(4)=AAAY(I1)/DY

            AAE(1)=AAAE(I)/DY
            AAE(2)=AENER(I)
            AAE(3)=AENER(I1)
            AAE(4)=AAAE(I1)/DY

            TIE=TAU
            CALL DXEE(AAD,AAX,AAY,AAE,DY,TIE,VIS,AFM,AFP,AFX,AFE)

            FM(I)=AFM
            FP(I)=AFP
            FPX(I)=AFX
            FE(I)=AFE

 998     CONTINUE

         DO II=3,IYEND
            DY=0.5*(YY(II+1)+YY(II))
     *           -0.5*(YY(II)+YY(II-1))
            ASDENS(II) = (FM(II-1)-FM(II))/DY+ADENS(II)
            ASXMOM(II) = (FP(II-1)-FP(II))/DY+AXMOM(II)
            ASYMOM(II) = (FPX(II-1)-FPX(II))/DY+AYMOM(II)
            ASENER(II) = (FE(II-1)-FE(II))/DY+AENER(II)
         ENDDO
         
         DO 176 IY=3,IYEND
            DENS(IX,IY)=ASDENS(IY)
            XMOM(IX,IY)=-ASYMOM(IY)
            YMOM(IX,IY)=ASXMOM(IY)
            ENER(IX,IY)=ASENER(IY)
 176     CONTINUE
 171  CONTINUE

      GOTO 601

      END
         



      SUBROUTINE DXE(AXU,AYU,AL,A,B1,B2,C,XQ,Y1Q,Y2Q,ZQ)
      
      CC=2.0*C-(AXU*AXU+AYU*AYU+2.5/AL)*A
      BB=B2-AYU*A
      AA=B1-AXU*A
      ZQ=0.400*AL*AL*(CC-2.0*AXU*AA-2.0*AYU*BB)
      Y1Q=2.0*AL*(AA-AXU*ZQ/AL)
      Y2Q=2.0*AL*(BB-AYU*ZQ/AL)
      XQ=A-Y1Q*AXU-Y2Q*AYU-ZQ*(AXU*AXU+AYU*AYU+2.5/AL)
      RETURN
      END
      
      SUBROUTINE DXEE(AD,AX,AY,AE,DD,TAU,VIS,AFM,AFP,AFPX,AFE)
      DIMENSION AD(4),AX(4),AY(4),AE(4)
      CK = 3.0
      PI=4.0*atan(1.0)
      
      ADE1=AD(2)+0.5*DD*AD(1)
      AXM1=AX(2)+0.5*DD*AX(1)
      AYM1=AY(2)+0.5*DD*AY(1)
      AEN1=AE(2)+0.5*DD*AE(1)
      
      ADE2=AD(3)-0.5*DD*AD(4)
      AXM2=AX(3)-0.5*DD*AX(4)
      AYM2=AY(3)-0.5*DD*AY(4)
      AEN2=AE(3)-0.5*DD*AE(4)

      RADE1=1.0/ADE1
      RADE2=1.0/ADE2

      AE1=0.25*(CK+2)*ADE1/(AEN1
     *     -0.50*AXM1*AXM1/ADE1-0.50*AYM1*AYM1/ADE1)
      AXU1=AXM1*RADE1
      AYU1=AYM1*RADE1

      AE2=0.25*(CK+2)*ADE2/(AEN2
     *     -0.50*AXM2*AXM2/ADE2-0.50*AYM2*AYM2/ADE2)
      AXU2=AXM2*RADE2
      AYU2=AYM2*RADE2

      AW11=AD(1)
      BXW11=AX(1)
      BYW11=AY(1)
      CW11=AE(1)

      AW22=AD(4)
      BXW22=AX(4)
      BYW22=AY(4)
      CW22=AE(4)



      IF(AE1.LE.0.0) THEN
         ADE1=AD(2)
         AXM1=AX(2)
         AYM1=AY(2)
         AEN1=AE(2)
         AW11=0.0
         BXW11=0.0
         BYW11=0.0
         CW11=0.0
         RADE1=1.0/ADE1
         AE1=0.25*(CK+2)*ADE1/(AEN1
     *        -0.50*AXM1*AXM1/ADE1
     *        -0.50*AYM1*AYM1/ADE1)
         AXU1=AXM1/ADE1
         AYU1=AYM1/ADE1
      endif

      IF(AE2.LE.0.0) THEN
         ADE2=AD(3)
         AXM2=AX(3)
         AYM2=AY(3)
         AEN2=AE(3)
         AW22=0.0
         BXW22=0.0
         BYW22=0.0
         CW22=0.0
         RADE2=1.0/ADE2

         AE2=0.25*(CK+2)*ADE2/(AEN2
     *        -0.50*AXM2*AXM2/ADE2
     *        -0.50*AYM2*AYM2/ADE2)
         AXU2=AXM2/ADE2
         AYU2=AYM2/ADE2
      ENDIF

      CALL DXE(AXU1,AYU1,AE1,AW11,BXW11,BYW11,CW11,X1,Y1,YX1,Z1)
      CALL DXE(AXU2,AYU2,AE2,AW22,BXW22,BYW22,CW22,X2,Y2,YX2,Z2)
      
            
      RAE1=1.0/AE1
      RAE2=1.0/AE2
      
      TEU0=0.5*DERFC(-AXU1*SQRT(AE1))
      TEU1=AXU1*TEU0+0.50*EXP(-AE1*AXU1*AXU1)/SQRT(AE1*PI)
      TEU2=AXU1*TEU1+0.5*TEU0*RAE1
      TEU3=AXU1*TEU2+1.0*TEU1*RAE1
      TEU4=AXU1*TEU3+1.5*TEU2*RAE1
      TEU5=AXU1*TEU4+2.0*TEU3*RAE1
      TEU6=AXU1*TEU5+2.5*TEU4*RAE1
          
      TGU0=0.5*(DERFC(AXU2*SQRT(AE2)))
      TGU1=AXU2*TGU0-0.50*EXP(-AE2*AXU2*AXU2)/SQRT(AE2*PI)
      TGU2=AXU2*TGU1+0.5*TGU0*RAE2
      TGU3=AXU2*TGU2+1.0*TGU1*RAE2
      TGU4=AXU2*TGU3+1.5*TGU2*RAE2
      TGU5=AXU2*TGU4+2.0*TGU3*RAE2
      TGU6=AXU2*TGU5+2.5*TGU4*RAE2


      EY0=1.0          
      EY1=AYU1
      EY2=AYU1*EY1+0.5*EY0*RAE1
      EY3=AYU1*EY2+1.0*EY1*RAE1
      EY4=AYU1*EY3+1.5*EY2*RAE1
      EY5=AYU1*EY4+2.0*EY3*RAE1
      EI2=0.5*CK*RAE1
      EI4=(0.75*CK+0.25*CK*(CK-1.0))*RAE1*RAE1


      GY0=1.0
      GY1=AYU2
      GY2=AYU2*GY1+0.5*GY0*RAE2
      GY3=AYU2*GY2+1.0*GY1*RAE2
      GY4=AYU2*GY3+1.5*GY2*RAE2
      GY5=AYU2*GY4+2.0*GY3*RAE2
      GI2=0.5*CK*RAE2
      GI4=(0.75*CK+0.25*CK*(CK-1.0))*RAE2*RAE2
            
c--- average state
      ADE=0.5*(AD(2)+AD(3))
      AXM=0.5*(AX(2)+AX(3))
      AYM=0.5*(AY(2)+AY(3))
      AEN=0.5*(AE(2)+AE(3))

      
      RADE=1.0/ADE
      AE0=0.25*(CK+2)*ADE/(AEN
     *     -0.50*AXM*AXM/ADE-0.50*AYM*AYM/ADE)
      AXU0=AXM*RADE
      AYU0=AYM*RADE
      RAE0=1.0/AE0

      AA1=ADE1/AE1
      AA2=ADE2/AE2

      TE=2.0*AE0*VIS+TAU*ABS(AA1-AA2)/(AA1+AA2)
      
      SP=-TAU/TE
      SW=EXP(SP)
      SE1=TE*(1.0-SW)
      SET0=TE*(-TAU+SE1)
      RSET0=1.0/SET0
      SET1=-TE*TAU*SW+TE*TE*(1.0-SW)
      SET2=TE*(-TAU+SE1)+SET1

      SETV=SET1

      FMI1=SE1*ADE1*TEU0
     *     -SETV*
     *     (X1*TEU1+Y1*TEU2+YX1*TEU1*EY1+
     *     Z1*(TEU3+TEU1*EY2+TEU1*EI2))

      FPI1=SE1*ADE1*TEU1
     *     -SETV*
     *     (X1*TEU2+Y1*TEU3+YX1*TEU2*EY1+
     *     Z1*(TEU4+TEU2*EY2+TEU2*EI2))

      FPIX1=SE1*ADE1*TEU0*EY1
     *     -SETV*
     *     (X1*TEU1*EY1+Y1*TEU2*EY1
     *     +YX1*TEU1*EY2+
     *     Z1*(TEU3*EY1+TEU1*EY3+TEU1*EY1*EI2))

      FEI1=0.5*SE1*ADE1*(TEU2+TEU0*EY2+TEU0*EI2)
     *     -0.5*SETV*
     *     (X1*(TEU3+TEU1*EY2+TEU1*EI2)+
     *     Y1*(TEU4+TEU2*EY2+TEU2*EI2)+
     *     YX1*(TEU3*EY1+TEU1*EY3+TEU1*EY1*EI2)+
     *     Z1*(TEU5+TEU1*EY4+TEU1*EI4
     *     +2.0*TEU3*EY2+2.0*TEU3*EI2+2.0*TEU1*EY2*EI2))
      
      FMI2=SE1*ADE2*TGU0
     *     -SETV*
     *     (X2*TGU1+Y2*TGU2+YX2*TGU1*GY1+
     *     Z2*(TGU3+TGU1*GY2+TGU1*GI2))

      FPI2=SE1*ADE2*TGU1
     *     -SETV*
     *     (X2*TGU2+Y2*TGU3+YX2*TGU2*GY1+
     *     Z2*(TGU4+TGU2*GY2+TGU2*GI2))

      FPIX2=SE1*ADE2*TGU0*GY1
     *     -SETV*
     *     (X2*TGU1*GY1+Y2*TGU2*GY1
     *     +YX2*TGU1*GY2+
     *     Z2*(TGU3*GY1+TGU1*GY3+TGU1*GY1*GI2))

      FEI2=0.5*SE1*ADE2*(TGU2+TGU0*GY2+TGU0*GI2)
     *     -0.5*SETV*
     *     (X2*(TGU3+TGU1*GY2+TGU1*GI2)+
     *     Y2*(TGU4+TGU2*GY2+TGU2*GI2)+
     *     YX2*(TGU3*GY1+TGU1*GY3+TGU1*GY1*GI2)+
     *     Z2*(TGU5+TGU1*GY4+TGU1*GI4
     *     +2.0*TGU3*GY2+2.0*TGU3*GI2+2.0*TGU1*GY2*GI2))
      
      AMI=FMI1+FMI2
      API=FPI1+FPI2
      APIX=FPIX1+FPIX2
      AEI=FEI1+FEI2


      FM1=ADE1*(SE1*TEU1)-(SETV)*
     *     (X1*TEU2+Y1*TEU3+YX1*TEU2*EY1
     *     +Z1*(TEU4+TEU2*EY2+TEU2*EI2))
      FP1=ADE1*(SE1*TEU2)-(SETV)*
     *     (X1*TEU3+Y1*TEU4+YX1*TEU3*EY1
     *     +Z1*(TEU5+TEU3*EY2+TEU3*EI2))
      FPX1=ADE1*(SE1*TEU1*EY1)-(SETV)*
     *     (X1*TEU2*EY1+Y1*TEU3*EY1
     *     +YX1*TEU2*EY2
     *     +Z1*(TEU4*EY1+TEU2*EY3+TEU2*EY1*EI2))
      FE1=ADE1*(0.5*SE1*(TEU3+TEU1*EI2+TEU1*EY2))
     *     -0.5*(SETV)*
     *     (X1*(TEU4+TEU2*EY2+TEU2*EI2)+
     *     Y1*(TEU5+TEU3*EY2+TEU3*EI2)+
     *     YX1*(TEU4*EY1+TEU2*EY3+TEU2*EY1*EI2)+
     *     Z1*(TEU6+TEU2*EY4+TEU2*EI4+2.0*TEU4*EY2
     *     +2.0*TEU4*EI2+2.0*TEU2*EY2*EI2))
      
      FM2=ADE2*(SE1*TGU1)-(SETV)*
     *     (X2*TGU2+Y2*TGU3+YX2*TGU2*GY1
     *     +Z2*(TGU4+TGU2*GY2+TGU2*GI2))
      FP2=ADE2*(SE1*TGU2)-(SETV)*
     *     (X2*TGU3+Y2*TGU4+YX2*TGU3*GY1
     *     +Z2*(TGU5+TGU3*GY2+TGU3*GI2))
      FPX2=ADE2*(SE1*TGU1*GY1)-(SETV)*
     *     (X2*TGU2*GY1+Y2*TGU3*GY1
     *     +YX2*TGU2*GY2
     *     +Z2*(TGU4*GY1+TGU2*GY3+TGU2*GY1*GI2))
      FE2=ADE2*(0.5*SE1*(TGU3+TGU1*GI2+TGU1*GY2))
     *     -0.5*(SETV)*
     *     (X2*(TGU4+TGU2*GY2+TGU2*GI2)+
     *     Y2*(TGU5+TGU3*GY2+TGU3*GI2)+
     *     YX2*(TGU4*GY1+TGU2*GY3+TGU2*GY1*GI2)+
     *     Z2*(TGU6+TGU2*GY4+TGU2*GI4+2.0*TGU4*GY2+2.0*TGU4*GI2+
     *     2.0*TGU2*GY2*GI2))
            


CCC---------one slope for g


      AW=(AD(3)-AD(2))/DD
      BXW=(AX(3)-AX(2))/DD
      BYW=(AY(3)-AY(2))/DD
      CW=(AE(3)-AE(2))/DD

      CALL DXE(AXU0,AYU0,AE0,AW,BXW,BYW,CW,XL,YL,YXL,ZL)

      T0=1.0
      T1=AXU0
      T2=AXU0*T1+0.5*T0*RAE0
      T3=AXU0*T2+1.0*T1*RAE0
      T4=AXU0*T3+1.5*T2*RAE0
      T5=AXU0*T4+2.0*T3*RAE0
      T6=AXU0*T5+2.5*T4*RAE0

      Y0=1.0          
      Y1=AYU0
      Y2=AYU0*Y1+0.5*Y0*RAE0
      Y3=AYU0*Y2+1.0*Y1*RAE0
      Y4=AYU0*Y3+1.5*Y2*RAE0
      Y5=AYU0*Y4+2.0*Y3*RAE0
      EI2=0.5*CK*RAE0
      EI4=(0.75*CK+0.25*CK*(CK-1.0))*RAE0*RAE0

            
      TRIU=(XL*T1+YL*T2+YXL*T1*Y1+
     *     ZL*(T3+T1*Y2+T1*EI2))
      TRIU2=(XL*T2+YL*T3+YXL*T2*Y1+
     *     ZL*(T4+T2*Y2+T2*EI2))
      TRIUY1=(XL*T1*Y1+YL*T2*Y1
     *     +YXL*T1*Y2+
     *     ZL*(T3*Y1+T1*Y3+T1*Y1*EI2))
      TRIEU=0.5*
     *     (XL*(T3+T1*Y2+T1*EI2)+
     *     YL*(T4+T2*Y2+T2*EI2)+
     *     YXL*(T3*Y1+T1*Y3+T1*Y1*EI2)+
     *     ZL*(T5+T1*Y4+T1*EI4
     *     +2.0*T3*Y2+2.0*T3*EI2+2.0*T1*Y2*EI2))
      

      TRIU3=(XL*T3+YL*T4+YXL*T3*Y1
     *     +ZL*(T5+T3*Y2+T3*EI2))
      TRIU2Y1=(XL*T2*Y1+YL*T3*Y1
     *     +YXL*T2*Y2
     *     +ZL*(T4*Y1+T2*Y3+T2*Y1*EI2))
      TRIEU2=0.5*
     *     (XL*(T4+T2*Y2+T2*EI2)+
     *     YL*(T5+T3*Y2+T3*EI2)+
     *     YXL*(T4*Y1+T2*Y3+T2*Y1*EI2)+
     *     ZL*(T6+T2*Y4+T2*EI4+2.0*T4*Y2+2.0*T4*EI2+
     *     2.0*T2*Y2*EI2))


ccc------------------------time evolution part


      AS=(-(TRIU)*SET2+SE1*ADE-AMI)*RSET0
      BS=(-(TRIU2)*SET2+SE1*AXM-API)*RSET0
      BSX=(-(TRIUY1)*SET2+SE1*AYM-APIX)*RSET0
      CS=(-(TRIEU)*SET2+SE1*AEN-AEI)*RSET0

      CALL DXE(AXU0,AYU0,AE0,AS,BS,BSX,CS,A1,B1,C1,D1)
            
CCFC----------------------------------------------------

      ATRIU=A1*T1+B1*T2+C1*T1*Y1+D1*(T3+T1*EI2+T1*Y2)
      ATRIUY1=A1*T1*Y1+B1*T2*Y1+C1*T1*Y2
     *     +D1*(T3*Y1+T1*EI2*Y1+T1*Y3)
      ATRIU2=A1*T2+B1*T3+C1*T2*Y1+D1*(T4+T2*EI2+T2*Y2)
      ATRIEU=0.5*(
     *     A1*(T3+T1*EI2+T1*Y2)+
     *     B1*(T4+T2*EI2+T2*Y2)+
     *     C1*(T3*Y1+T1*EI2*Y1+T1*Y3)+
     *     D1*(T5+T1*Y4+T1*EI4+2.0*T3*Y2
     *     +2.0*T3*EI2+2.0*T1*EI2*Y2))
      
      A=ADE*(TAU-SE1)
      B=SET2
      D=(0.5*TAU*TAU
     *     -TE*TAU+TE*TE*(1.0-SW))
      
      AFM=A*T1+B*TRIU2+D*ATRIU+(FM1+FM2)
      AFP=A*T2+B*TRIU3+D*ATRIU2+(FP1+FP2)
      AFPX=A*T1*Y1+B*TRIU2Y1+D*ATRIUY1
     *     +(FPX1+FPX2)
      AFE=0.5*A*(T3+T1*EI2+T1*Y2)+B*TRIEU2+
     *     D*ATRIEU+(FE1+FE2)

      RETURN
      END



