C     Last change:  LE   15 Oct 2001   12:45 pm
      SUBROUTINE CalcShad(task,
     &     sol_dir_deg,sol_elev_deg,
     &     Coord_Surface,
     &     N_Shades, N_Points, Coord_Shades, Transp_Shades,
     &     Transp_Result,Three,Four)

* Updates
* 980315 AB Sun's position [deg & rad], in separate variables
* 980306 ML Major revision

C
      DOUBLE PRECISION pi
      PARAMETER(pi = 3.141592653589793)
C
cabg  DOUBLE PRECISION sol_dir, sol_elev
      DOUBLE PRECISION sol_dir_deg, sol_elev_deg
      DOUBLE PRECISION Coord_Surface(3,4)
      INTEGER          task, N_Shades, N_points, Three, Four
      DOUBLE PRECISION Coord_Shades(3,N_Shades*4)
      DOUBLE PRECISION Transp_Shades(N_Shades)
      DOUBLE PRECISION Transp_Result
C
C * Inputs
C     task           - 1 for direct,
C                      2 for diffuse, which allows sun below horizon
C     sol_dir_deg    - Direction to the sun [deg]
C     sol_elev_deg   - Elevation of the sun [deg]
C     Coord_Surface  - Coordinates of the surface to project upon
C     N_Shades       - Number of shades (rectangles)
C     Coord_Shades   - Coordinates for the shades
C     Transp_Shades  - Transparency of each shade [0.0, 1.0]
C                        1.0 = transparent
C * Outputs:
C     Transp_Result  - Result as equivalent transparency [0.0,1.0]
C                        1.0 = transparent
C
      INTEGER I, Ii, Iii, nn_shades
      INTEGER, DIMENSION(N_Shades) :: shade_index
      DOUBLE PRECISION sol_dir, sol_elev
      DOUBLE PRECISION, DIMENSION(4,N_Shades*4) ::  Coord_Proj
      DOUBLE PRECISION, DIMENSION(4,4) ::  Coord_Surf_Proj
      DOUBLE PRECISION Trans(4,4)
      DOUBLE PRECISION Transx(4,4), Transy(4,4), Transz(4,4)
      DOUBLE PRECISION Translat(4,4)
      DOUBLE PRECISION X_Rot, Y_rot, Z_rot
      DOUBLE PRECISION Sol_Vec(4), R, Normal_to_surf(4)
      DOUBLE PRECISION Surf_X(4), Surf_Y(4), Surf_tmp(4)
      DOUBLE PRECISION Scalar_Product
      DOUBLE PRECISION C1(3),C2(3),C3(3),C4(3)
      LOGICAL Debug,debug1
C
C * Initializations
C
      Debug = .FALSE.
      Debug1 = .FALSE.
      IF(DEBUG) THEN
        WRITE(*,*) '==============='
        Write(*,*) 'S ',sol_dir_deg, sol_elev_deg
      END IF

      IF(DEBUG1)
     & Write(127,'(a,2f10.2,a)') '''(',sol_dir_deg, sol_elev_deg,')'
C
C * Transform Sun position to a vector
C
      Sol_dir = sol_dir_deg * PI / 180.0
      sol_elev = sol_elev_deg * PI / 180.0
      Sol_Vec =
     &  (/ Sin(Sol_Dir) * Cos(Sol_Elev),
CCCC switch!
     &     Cos(Sol_Dir) * Cos(Sol_Elev),
     &     Sin(Sol_Elev),
     &     1.0D0 /)
c>      Debug = .TRUE.
      Transp_Result = 0.0D0
cle>  IF(sol_elev_deg .LE. 0.0D0 .AND. task .EQ. 0) RETURN
      IF(sol_elev_deg .LE. 0.0D0 .AND. task .EQ. 1) RETURN
C
C 1 3, 4 3
C It is only interesting to align window with XY, hence the Abs
C ?? Very strange order of points in Coord_Surf:
C ??  4 2
C ??  3 1
C
      Surf_X(1:3) = Abs(Coord_Surface(:,1) - Coord_Surface(:,3))
      Surf_X(4) = 1.0
      Surf_Y(1:3) = Abs(Coord_Surface(:,4) - Coord_Surface(:,3))
      Surf_Y(4) = 1.0
C
      CALL Cross_product(Normal_To_Surf, Surf_Y, Surf_X)
C fix
C ?? this is wrong if the vector has coordinates with different signs ??
      Normal_to_surf = ABS(Normal_to_surf)
      IF(DEBUG) THEN
        Write(*,'(1x,A,1x,3F6.2)') 'Sol_Vec',Sol_Vec(1:3)
        Write(*,'(1x,A,1x,3F6.2)') 'Coord_surface',Coord_surface(1:3,1)
        Write(*,'(1x,A,1x,3F6.2)') 'Coord_surface',Coord_surface(1:3,2)
        Write(*,'(1x,A,1x,3F6.2)') 'Coord_surface',Coord_surface(1:3,3)
        Write(*,'(1x,A,1x,3F6.2)') 'Coord_surface',Coord_surface(1:3,4)
        Write(*,'(1x,A,1x,3F6.2)') 'Normal_To_Surf',Normal_To_Surf(1:3)
        Write(*,'(1x,A,1x,3F6.2)') 'Scalar_product',
     & Scalar_product(Sol_Vec, Normal_To_Surf)
      ENDIF
      IF (Scalar_product(Sol_Vec, Normal_To_Surf) .LE. 0.0D0) RETURN
C
C * Make transformation matrix for the surface on the XY-plane
C
C     * Normalize the surface x and y
C
      R = Sqrt(
     &     Surf_X(1)*Surf_X(1) +
     &     Surf_X(2)*Surf_X(2) +
     &     Surf_X(3)*Surf_X(3))
      Surf_X(1:3) = Surf_X(1:3)/R
      R = Sqrt(
     &     Surf_Y(1)*Surf_Y(1) +
     &     Surf_Y(2)*Surf_Y(2) +
     &     Surf_Y(3)*Surf_Y(3))
      Surf_Y(1:3) = Surf_Y(1:3)/R
C
C     * Initialize the transformation matrix
C
      Translat = 0.0
      DO i = 1,4
         Translat(i,i) = 1.0
      END DO
      Transx = Translat
      Transz = Translat
      Transy = Translat
C
C     * Translate
C
      Translat(1:3,4) = - Coord_Surface(:,3)
C
C     * Rotate so that the Surf_x is pointing in the X -direction (1,0,0)
C
      IF(Debug)THEN
         Write(*,*) 'Surf_X'
         Write(*,'(4f6.2)') Surf_X
      END IF
      IF(Surf_X(2) .EQ. 0.0) THEN
         X_rot = 0.0
      ELSE
CC         X_Rot = Atan2(Surf_X(2),Surf_X(3))
         X_Rot = Atan2(Surf_X(3),Surf_X(2))
         IF(DEBUG)THEN
            write(*,'(3f6.2)') X_Rot, Surf_X(3), Surf_X(2)
         END IF
      END IF
      Transx(2,2) = Cos(X_Rot)
      Transx(3,3) = Cos(X_Rot)
      Transx(2,3) = -Sin(X_Rot)
      Transx(3,2) = Sin(X_Rot)
      Surf_X=MATMUL(Transx, Surf_X)
      Surf_Y=MATMUL(Transx, Surf_Y)
      IF(Debug)THEN
         Write(*,'(A,4f6.2)') ' transx on X ',Surf_X
      END IF
      IF(Surf_X(3) .EQ. 0.0) THEN
         Y_rot = 0.0
      ELSE
         Y_rot = Atan2(Surf_X(3),Surf_X(1))
      END IF
      Transz(1,1) = Cos(Y_rot)
      Transz(3,3) = Cos(Y_rot)
      Transz(1,3) = Sin(Y_rot)
      Transz(3,1) = -Sin(Y_rot)
      Surf_tmp=MATMUL(Transz, Surf_X)
      IF(Surf_tmp(1) .LT. 0) THEN
         Y_rot = Y_rot+pi
         Transz(1,1) = Cos(Y_rot)
         Transz(3,3) = Cos(Y_rot)
         Transz(1,3) = Sin(Y_rot)
         Transz(3,1) = -Sin(Y_rot)
      END IF
      Surf_X=MATMUL(Transz, Surf_X)
      Surf_Y=MATMUL(Transz, Surf_Y)
      IF(Debug)THEN
         Write(*,'(A,4f6.2)') ' transz onX ',Surf_X
      END IF
C
C     * Rotate so that the Surf_y is pointing in the y -direction (0,1,0)
C
      IF(Debug)THEN
         Write(*,*) 'Surf_y'
         Write(*,'(4f6.2)') Surf_y
      END IF
      IF(Surf_y(3) .EQ. 0.0) THEN
         Z_rot = 0.0
      ELSE
         Z_Rot = Atan2(Surf_Y(3),Surf_Y(2))
      END IF
      Transy(2,2) = Cos(Z_Rot)
      Transy(3,3) = Cos(Z_Rot)
      Transy(2,3) = Sin(Z_Rot)
CCC * Fix!
      Transy(3,2) = Abs(-Sin(Z_Rot))
C>      Transy(2,3) = -Sin(Z_Rot)
C>      Transy(3,2) = Sin(Z_Rot)
C>
      Surf_x=MATMUL(Transy, Surf_x)
      Surf_Y=MATMUL(Transy, Surf_Y)
      IF(Debug)THEN
         Write(*,'(A,4f6.2)') ' transy on X ',Surf_x
         Write(*,'(A,4f6.2/)') ' transy on Y ',Surf_Y
      END IF
C
      Trans = MATMUL(Transx,Translat)
      Trans = MATMUL(Transz, Trans)
      Trans = MATMUL(Transy, Trans)
C
      IF(Debug)THEN
         Write(*,'(A,/4(4f6.2/)/)') ' Translat ',(Translat(:,i),i=1,4)
         Write(*,'(A,/4(4f6.2/)/)') ' Transx   ',(Transx(:,i),i=1,4)
         Write(*,'(A,/4(4f6.2/)/)') ' Transy   ',(Transy(:,i),i=1,4)
         Write(*,'(A,/4(4f6.2/)/)') ' Trans    ',(Trans(:,i),i=1,4)
      END IF
C
C     * Transform the surface
C
      Coord_Surf_Proj(1:3,:) = Coord_Surface
      Coord_Surf_Proj(4,:) = 1.0
      Coord_Surf_Proj = MATMUL(Trans, Coord_Surf_Proj)
      IF(Debug)THEN
         Write(*,*) 'Projected surface'
         Write(*,*) Nint(Coord_Surf_Proj(1:2,1)*100)/100.0
         Write(*,*) Nint(Coord_Surf_Proj(1:2,2)*100)/100.0
         Write(*,*) Nint(Coord_Surf_Proj(1:2,3)*100)/100.0
         Write(*,*) Nint(Coord_Surf_Proj(1:2,4)*100)/100.0
      END IF
C
C * Transform every Shade by transforming the coordinate vector
C
      Coord_Proj(1,:) = Coord_Shades(1,:)
      Coord_Proj(2,:) = Coord_Shades(2,:)
      Coord_Proj(3,:) = Coord_Shades(3,:)
      Coord_Proj(4,:) = 1.0
      IF(Debug)THEN
        write(*,*) 'Coord_proj before trans'
        do i=1,n_shades*4
          Write(*,*) Coord_proj(1:3,I)
        end do
      END IF
C
      Coord_Proj = MATMUL(Trans, Coord_Proj)
C
      IF(Debug)THEN
        write(*,*) 'Coord_proj after trans'
        do i=1,n_shades*4
          Write(*,*) Coord_proj(1:3,I)
        end do
      END IF
C
C * Cut with XY-plane
C
C ?? The intersection may be pentagonal ??

      nn_shades = 0
      DO I=1,n_shades
        Ii = (I-1)*4
        C1=Coord_proj(1:3,Ii+1)
        C2=Coord_proj(1:3,Ii+2)
        C3=Coord_proj(1:3,Ii+3)
        C4=Coord_proj(1:3,Ii+4)
Cabg?? Ok to update on the fly?
        call CutZZero(C1,C2)
        call CutZZero(C2,C3)
        call CutZZero(C3,C4)
        call CutZZero(C4,C1)
        Coord_proj(1:3,Ii+1)=C1
        Coord_proj(1:3,Ii+2)=C2
        Coord_proj(1:3,Ii+3)=C3
        Coord_proj(1:3,Ii+4)=C4
        IF(C1(3) .LE. 0 .AND. C2(3) .LE. 0 .AND.
     &     C3(3) .LE. 0 .AND. C4(3) .LE. 0) THEN
C * Discard
        ELSE
C * Cut
           nn_shades=nn_shades+1 
           shade_index(nn_shades)=I
        END IF
      end do

      do i=1,nn_shades
        Ii = (i-1)*4
        Iii = (shade_index(i)-1)*4
        Coord_proj(1:3,Ii+1)=Coord_proj(1:3,Iii+1)
        Coord_proj(1:3,Ii+2)=Coord_proj(1:3,Iii+2)
        Coord_proj(1:3,Ii+3)=Coord_proj(1:3,Iii+3)
        Coord_proj(1:3,Ii+4)=Coord_proj(1:3,Iii+4)
      end do
      n_shades = nn_shades
C
      IF(Debug)THEN
        write(*,*) 'Coord_proj efter CutZZero'
        do i=1,n_shades*4
          Write(*,*) Coord_proj(1:3,I)
        end do
      END IF
C
C * Only rotation of the solar direction
C
      Sol_vec(1:3) = MATMUL(Trans(1:3,1:3),Sol_vec(1:3))
      IF(Debug)THEN
         Write(*,*) 'X_Rot,Y_rot, Sol_Dir, Sol_Elev'
         Write(*,*) X_Rot,Y_rot, Sol_Dir, Sol_Elev
         Write(*,*) 'Sol_Vec'
         Write(*,'(4f6.2)') Sol_Vec
      END IF
C
C * Trans is reused for transformation without translation
CC
CC Sol_dir is in the right direction
CC
CC      Trans = MATMUL(Transz, Transx)
CC      Trans = MATMUL(Transy, Trans)
CC      Sol_Vec = MATMUL(Trans, Sol_Vec)
CC      IF(Debug)THEN
CC         Write(*,*) Sol_Vec
CC      END IF
C
C * Project every point in the shades to the transformed surface (Z=0)
C
      IF(Debug)THEN
         Write(*,*) 'Coord_Proj'
         IF(debug1) write(127,*) '(mapcar #''draw-1 ''('
      END IF
      DO I = 1, 4*N_Shades
         R=Coord_Proj(3,I)/Sol_Vec(3)
         Coord_Proj(1:3,I) = Coord_Proj(1:3,I) - R * Sol_Vec(1:3)
         Coord_Proj(1:3,I) = Nint(Coord_Proj(1:3,I)*100) / 100.0
         IF(Debug)THEN
            Write(*,'(1x,2F7.2)') Coord_Proj(1:2,I)
         END IF
         IF(Debug1)THEN
            if(mod(I,4).eq.1) write(127,*) ' ('
            write(127,'(3x,a1,3F6.2,a1)') '(',Coord_Proj(1:3,I),')'
            if(mod(I,4).eq.0) write(127,*) ' )'
         END IF
      END DO
      IF(Debug1)THEN
         write(127,*) '))'
      END IF
C
C * Scanning
C
      CALL Scan_and_cut(
     &     Coord_surf_proj, N_Shades, Coord_proj,
     &     Transp_Shades, Transp_Result)
C>      Debug = .TRUE.
      IF(DEBUG) WRITE(*,*) Transp_Result
      RETURN
      END

      SUBROUTINE Scan_and_cut(
     &     Coord_surface, N_Shades, Coord_shades,
     &     Transp_Shades, Transp_Result)
C
      DOUBLE PRECISION Coord_Surface(4,4)
      INTEGER N_Shades
      DOUBLE PRECISION Coord_Shades(4,N_Shades*4)
      DOUBLE PRECISION Transp_Shades(N_Shades)
      DOUBLE PRECISION Transp_Result
C
C * Inputs
C     Coord_Surface    - Coordinates of the surface to project upon
C     N_Shades         - Number of shades (rectangles)
C     Coord_Shades     - Coordinates for the shades
C     Transp_Shades    - Transparency of each shade [0.0, 1.0]
C                        1.0 = transparent
C * Outputs:
C     Transp_Result    - Result as equivalent transparency [0.0,1.0]
C
      LOGICAL Debug
      INTEGER Ip1, Ip2, Ip3, Ip4
      DOUBLE PRECISION Xmin, Xmax, Ymin, Ymax, Total_area
      DOUBLE PRECISION Xp1, Xp2, Xp3, Xp4, Yp1, Yp2, Yp3, Yp4
      DOUBLE PRECISION Eps
      PARAMETER (Eps=1.0D-8)
C
C * Pointers etc.
C
      INTEGER Nullp, SUCC, PRE, DAT, TRANS, AREA
      PARAMETER (Nullp=-1,
     &           SUCC=1, PRE=2, DAT=3, TRANS=1, AREA=2)
      INTEGER MAX_Surf, MAX_Pnt, Max_Coord
      PARAMETER (MAX_Surf=1000,
     &           MAX_Pnt=4000,
     &           Max_Coord=4000)
      INTEGER M_Surf, M_Pnt, M_Coord
      INTEGER Slistp(Max_Surf), Slist(DAT,Max_Surf)
      INTEGER Clistp(Max_Coord)
      DOUBLE PRECISION Clist(2,Max_Coord)
      DOUBLE PRECISION RSlist(AREA,Max_Surf)
      INTEGER Plistp(Max_Pnt), Plist(3,Max_Pnt)
C Functions
      INTEGER New_Coord, New_Surf
      DOUBLE PRECISION Trans_surfaces
C
      INTEGER Root_out, Rootp, Root_scan
      INTEGER I_SurfP, I_Coord, I_Surfout, I_Scan
      INTEGER I_P1, I_P2, I_PS, I, I_PX
      DOUBLE PRECISION Coord(2), Coord1(2), Coord2(2), Coordx(2), R
      DOUBLE PRECISION CoordQ1(2), CoordQ2(2)
      INTEGER I_SurfQ, I_Q1, I_Q2, I_PScan
      LOGICAL intersect, Lmin
      DOUBLE PRECISION Trans_w_size
C
C * SList (1,:) = SUCC,
C *       (2,:) = PRE,
C *       (3,:) = DAT -> PList,  NullP = -1
C
C * PList (1,:) = SUCC,
C *       (2,:) = PRE,
C *       (3,:) = DAT -> CList,  NullP = -1
C
C * CList (1,:) = X
C * CList (2,:) = Y
C
C *RSList (1,:) = TRANSP,
C *       (2,:) = AREA,
C
C * Initializations
C
      Debug = .FALSE.
      IF (DEBUG) THEN
          WRITE(*,*) 'Scan_and_cut'
      ENDIF
C
      Rootp = Nullp
      Root_out = Nullp
      Root_scan = Nullp
C      
      DO I = 1, Max_Surf
          Slistp(I)=I
      END DO
      M_Surf=0

      DO I = 1, Max_Coord
          Clistp(I)=I
      END DO
      M_Coord=0

      DO I = 1, Max_Pnt
          Plistp(I)=I
      END DO
      M_Pnt=0
C
      Trans_w_size = 0.0D0
C ======================================================================
C
C * Move every shade that is potentially inside the surface to
C * a list of shadeareas SList
C
      Xmin = Coord_Surface(1,1)
      Xmax = Coord_Surface(1,3)
      Ymin = Coord_Surface(2,1)
      Ymax = Coord_Surface(2,3)
      DO I=1,4
         Xmin = Min(Xmin,Coord_Surface(1,I))
         Xmax = Max(Xmax,Coord_Surface(1,I))
         Ymin = Min(Ymin,Coord_Surface(2,I))
         Ymax = Max(Ymax,Coord_Surface(2,I))
      END DO
      Total_area = (Xmax - Xmin) * (Ymax - Ymin)
      DO I=1,N_Shades
        Ip1 = (I-1)*4+1
        Ip2 = (I-1)*4+2
        Ip3 = (I-1)*4+3
        Ip4 = (I-1)*4+4
        Xp1 = Coord_Shades(1,Ip1)
        Yp1 = Coord_Shades(2,Ip1)
        Xp2 = Coord_Shades(1,Ip2)
        Yp2 = Coord_Shades(2,Ip2)
        Xp3 = Coord_Shades(1,Ip3)
        Yp3 = Coord_Shades(2,Ip3)
        Xp4 = Coord_Shades(1,Ip4)
        Yp4 = Coord_Shades(2,Ip4)
        IF (
     &     (Xmin .GT. XP1 .AND. Xmin .GT. XP2 .AND.
     &      Xmin .GT. XP3 .AND. Xmin .GT. XP4)) THEN
C Discard left
CC           Write(*,*) 'l',Xp1,yp1
CC           Write(*,*) Xp2,yp2
CC           Write(*,*) Xp3,yp3
CC           Write(*,*) Xp4,yp4
        ELSE IF (
     &     (Xmax .LT. XP1 .AND. Xmax .LT. XP2 .AND.
     &      Xmax .LT. XP3 .AND. Xmax .LT. XP4)) THEN
C Discard right
CC           Write(*,*) 'r',Xp1,yp1
CC           Write(*,*) Xp2,yp2
CC           Write(*,*) Xp3,yp3
CC           Write(*,*) Xp4,yp4
        ELSE IF (
     &     (Ymin .GT. YP1 .AND. Ymin .GT. YP2 .AND.
     &      Ymin .GT. YP3 .AND. Ymin .GT. YP4)) THEN
C Discard below
CC           Write(*,*) 'b',Xp1,yp1
CC           Write(*,*) Xp2,yp2
CC           Write(*,*) Xp3,yp3
CC           Write(*,*) Xp4,yp4
        ELSE IF (
     &     (Ymax .LT. YP1 .AND. Ymax .LT. YP2 .AND.
     &      Ymax .LT. YP3 .AND. Ymax .LT. YP4)) THEN
C Discard above
CC           Write(*,*) 'a',Xp1,yp1
CC           Write(*,*) Xp2,yp2
CC           Write(*,*) Xp3,yp3
CC           Write(*,*) Xp4,yp4
        ELSE
C Do not discard!
          I_SurfP = New_surf(M_Surf,Rootp,Slistp,Slist,RSlist,
     &          Transp_Shades(I))
          I_Coord = New_coord(M_Coord,Clistp,Clist,Coord_Shades(1,Ip1))
          Call Add_coord(I_SurfP,I_Coord,Slist,M_Pnt,Plistp,Plist)
          I_Coord = New_coord(M_Coord,Clistp,Clist,Coord_Shades(1,Ip2))
          Call Add_coord(I_SurfP,I_Coord,Slist,M_Pnt,Plistp,Plist)
          I_Coord = New_coord(M_Coord,Clistp,Clist,Coord_Shades(1,Ip3))
          Call Add_coord(I_SurfP,I_Coord,Slist,M_Pnt,Plistp,Plist)
          I_Coord = New_coord(M_Coord,Clistp,Clist,Coord_Shades(1,Ip4))
          Call Add_coord(I_SurfP,I_Coord,Slist,M_Pnt,Plistp,Plist)
          IF (debug) THEN
            WRITE(*,*) 'Added shade'
            CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
          END IF
        END IF
      END DO
C
C * Now we have all surfaces stored in Slist beginning from Rootp
C
C ======================================================================
C
C * Start by cutting against the window surface
C * 1st left & right against Xmin & Xmax
C
      I_SurfP = Rootp
 2000 CONTINUE
      IF(I_SurfP .EQ. Nullp) GOTO 2050
      I_P1=Slist(DAT,I_SurfP)
C>ML:980904 moved down>      Coord1(:) = Clist(:,Plist(DAT,I_P1))
      IF(I_P1 .EQ. Nullp) THEN
         I_SurfP = Slist(SUCC,I_SurfP)
         GOTO 2000
      END IF
C>>
      Coord1(:) = Clist(:,Plist(DAT,I_P1))
 2100 CONTINUE
      I_P2=Plist(SUCC,I_P1)
      IF (I_P2 .EQ. Nullp) THEN
         Coord2(:) = Clist(:,Plist(DAT,Slist(DAT,I_SurfP)))
      ELSE
         Coord2(:) = Clist(:,Plist(DAT,I_P2))
      ENDIF
C cut Left and Right
      IF( Abs(Coord1(1)-Coord2(1)) .GT. Eps ) THEN
        Lmin = Coord1(1) .LE. Coord2(1)
        I_PX = I_P1
CC      write(*,*) 'xm',Xmin,Xmax
CC      write(*,*) 'coord1',coord1
CC      write(*,*) 'coord2',coord2
        R=(Xmin-Coord1(1))/(coord2(1)-Coord1(1))
        IF (R .GE. 0.0 .AND. R+eps .LT. 1.0) THEN
          Coordx(1)=Xmin+eps
          Coordx(2)=R*(Coord2(2)-Coord1(2))+Coord1(2)
          IF(Lmin)THEN
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coordx)
             Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
             Coord(:)=Coordx(:)-(/2*eps,0.0D0/)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
             Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
             I_PX=Plist(SUCC,Plist(SUCC,I_P1))
          ELSE
             Coord(:)=Coordx(:)-(/2*eps,0.0D0/)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
             Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coordx)
             Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
          END IF
          IF (debug) THEN
            WRITE(*,*) '2100 cut left'
            CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
          END IF
        END IF       
        R=(Xmax-Coord1(1))/(coord2(1)-Coord1(1))
        IF (R .GE. 0.0 .AND. R+eps .LT. 1.0) THEN
          Coordx(1)=Xmax-eps
          Coordx(2)=R*(Coord2(2)-Coord1(2))+Coord1(2)
          IF(.NOT. Lmin)THEN
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coordx)
             Call Add_coord_after(I_PX,I_Coord,Slist,M_Pnt,Plistp,Plist)
             Coord(:)=Coordx(:)+(/2*eps,0.0D0/)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
             Call Add_coord_after(I_PX,I_Coord,Slist,M_Pnt,Plistp,Plist)
          ELSE
             Coord(:)=Coordx(:)+(/2*eps,0.0D0/)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
             Call Add_coord_after(I_PX,I_Coord,Slist,M_Pnt,Plistp,Plist)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coordx)
             Call Add_coord_after(I_PX,I_Coord,Slist,M_Pnt,Plistp,Plist)
          END IF       
          IF (debug) THEN
            WRITE(*,*) '2100 cut right'
            CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
          END IF
        END IF
CC        write(*,*) 'coords lr',coordx
      END IF
      Coord1(:)=Coord2(:)
      I_P1 = I_P2
      IF (I_P1 .NE. Nullp) GOTO 2100
      I_SurfP = Slist(SUCC,I_SurfP)
      GOTO 2000
 2050 CONTINUE
C>
CC      CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
C ----------------------------------------------------------
C Remove right and left points
      I_SurfP = Rootp
 3000 CONTINUE
      IF(I_SurfP .EQ. Nullp) GOTO 3050
      I_P1=Slist(DAT,I_SurfP)
 3100 CONTINUE
      IF(I_P1 .EQ. Nullp) THEN
         I_SurfP = Slist(SUCC,I_SurfP)
         GOTO 3000
      END IF
      I_Coord = Plist(DAT,I_P1)
      Coord1(:) = Clist(:,I_Coord)
C cut Left and Right
      IF( Coord1(1) .LT. Xmin ) THEN
C        CALL Remove_coord(I_Coord,M_Coord,Clistp)
        I_PS = Plist(PRE,I_P1)
        CALL Remove_point(I_P1,I_SurfP,M_Pnt,Plistp,Plist,Slist)
        IF( I_PS .NE. Nullp)THEN
          I_P1 = I_PS
        ELSE
          I_P1=Slist(DAT,I_SurfP)
          GOTO 3100
        END IF
      ELSE IF( Coord1(1) .GT. Xmax ) THEN
C        CALL Remove_coord(I_Coord,M_Coord,Clistp)
        I_PS = Plist(PRE,I_P1)
        CALL Remove_point(I_P1,I_SurfP,M_Pnt,Plistp,Plist,Slist)
        IF( I_PS .NE. Nullp)THEN
          I_P1 = I_PS
        ELSE
          I_P1=Slist(DAT,I_SurfP)
          GOTO 3100
        END IF
      ELSE
      END IF
      I_P1 = Plist(SUCC,I_P1)
      IF (I_P1 .NE. Nullp) GOTO 3100
      I_SurfP = Slist(SUCC,I_SurfP)
      GOTO 3000
 3050 CONTINUE
      IF (debug) THEN
        WRITE(*,*) '3050 removed left & right'
        CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
      END IF
C ======================================================================
C
C * Now bottom & top against Ymin & Ymax
C
      I_SurfP = Rootp
 2200 CONTINUE
      IF(I_SurfP .EQ. Nullp) GOTO 2250
      I_P1=Slist(DAT,I_SurfP)
      IF(I_P1 .EQ. Nullp) THEN
         I_SurfP = Slist(SUCC,I_SurfP)
         GOTO 2200
      END IF
      Coord1(:) = Clist(:,Plist(DAT,I_P1))
 2300 CONTINUE
      I_P2=Plist(SUCC,I_P1)
      IF (I_P2 .EQ. Nullp) THEN
         Coord2(:) = Clist(:,Plist(DAT,Slist(DAT,I_SurfP)))
      ELSE
         Coord2(:) = Clist(:,Plist(DAT,I_P2))
      ENDIF
C cut Up and down
      IF( Abs(Coord1(2)-Coord2(2)) .GT. Eps ) THEN
        Lmin = Coord1(2) .LE. Coord2(2)
        I_PX = I_P1
C      write(*,*) 'ym',Ymin,Ymax
C      write(*,*) 'coord1',coord1
C      write(*,*) 'coord2',coord2
        R=(Ymin-Coord1(2))/(coord2(2)-Coord1(2))
        IF (R .GE. 0.0 .AND. R+eps .LT. 1.0) THEN
          Coordx(2)=Ymin+eps
          Coordx(1)=R*(Coord2(1)-Coord1(1))+Coord1(1)
          IF(Lmin)THEN
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coordx)
             Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
             Coord(:)=Coordx(:)-(/0.0D0,2*eps/)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
             Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
             I_PX=Plist(SUCC,Plist(SUCC,I_P1))
          ELSE
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coordx)
             Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
             Coord(:)=Coordx(:)-(/0.0D0,2*eps/)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
             Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
          END IF
          IF (debug) THEN
            WRITE(*,*) '2300 cut bottom'
            CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
          END IF
        END IF       
        R=(Ymax-Coord1(2))/(coord2(2)-Coord1(2))
        IF (R .GE. 0.0 .AND. R+eps .LT. 1.0) THEN
          Coordx(2)=Ymax-eps
          Coordx(1)=R*(Coord2(1)-Coord1(1))+Coord1(1)
          IF(.NOT. Lmin)THEN
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coordx)
             Call Add_coord_after(I_PX,I_Coord,Slist,M_Pnt,Plistp,Plist)
             Coord(:)=Coordx(:)+(/0.0D0,2*eps/)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
             Call Add_coord_after(I_PX,I_Coord,Slist,M_Pnt,Plistp,Plist)
          ELSE
             Coord(:)=Coordx(:)+(/0.0D0,2*eps/)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
             Call Add_coord_after(I_PX,I_Coord,Slist,M_Pnt,Plistp,Plist)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coordx)
             Call Add_coord_after(I_PX,I_Coord,Slist,M_Pnt,Plistp,Plist)
          END IF
          IF (debug) THEN
            WRITE(*,*) '2300 cut top'
            CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
          END IF
        END IF       
CC        write(*,*) 'coords bt',coordx
      END IF
C      
      Coord1(:)=Coord2(:)
      I_P1 = I_P2
      IF (I_P1 .NE. Nullp) GOTO 2300
      I_SurfP = Slist(SUCC,I_SurfP)
      GOTO 2200
 2250 CONTINUE
C>
C      CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
C ----------------------------------------------------------
C
C Remove top and bottom points
C
      I_SurfP = Rootp
 3200 CONTINUE
      IF(I_SurfP .EQ. Nullp) GOTO 3250
      I_P1=Slist(DAT,I_SurfP)
 3300 CONTINUE
      IF(I_P1 .EQ. Nullp) THEN
         I_SurfP = Slist(SUCC,I_SurfP)
         GOTO 3200
      END IF
      I_Coord = Plist(DAT,I_P1)
      Coord1(:) = Clist(:,I_Coord)
C cut Top and Down
      IF( Coord1(2) .LT. Ymin ) THEN
C        CALL Remove_coord(I_Coord,M_Coord,Clistp)
        I_PS = Plist(PRE,I_P1)
        CALL Remove_point(I_P1,I_SurfP,M_Pnt,Plistp,Plist,Slist)
        IF( I_PS .NE. Nullp)THEN
          I_P1 = I_PS
        ELSE
          I_P1=Slist(DAT,I_SurfP)
          GOTO 3300
        END IF
      ELSE IF( Coord1(2) .GT. Ymax ) THEN
C        CALL Remove_coord(I_Coord,M_Coord,Clistp)
        I_PS = Plist(PRE,I_P1)
        CALL Remove_point(I_P1,I_SurfP,M_Pnt,Plistp,Plist,Slist)
        IF( I_PS .NE. Nullp)THEN
          I_P1 = I_PS
        ELSE
          I_P1=Slist(DAT,I_SurfP)
          GOTO 3300
        END IF
      ELSE
      END IF
      I_P1 = Plist(SUCC,I_P1)
      IF (I_P1 .NE. Nullp) GOTO 3300
      I_SurfP = Slist(SUCC,I_SurfP)
      GOTO 3200
 3250 CONTINUE
C
C>
      IF(Debug)THEN
        write(*,*) 'Rootp scan'
        CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
        CALL Remove_duplicates (Rootp,
     &  Slistp,Slist,Plistp,Plist,Clist,Clistp,M_Coord,M_Pnt,eps)
      END IF
C
C ======================================================================
C
C * Remove void surfaces
C * Add intersects between surfaces
C
      I_SurfP = Rootp
 7400 CONTINUE
      IF(I_SurfP .EQ. Nullp) GOTO 7450
      I_P1=Slist(DAT,I_SurfP)
      IF(I_P1 .EQ. Nullp) THEN
         I_SurfP = Slist(SUCC,I_SurfP)
         GOTO 7400
      END IF
      Coord1(:) = Clist(:,Plist(DAT,I_P1))
 7500 CONTINUE
      I_P2=Plist(SUCC,I_P1)
      IF (I_P2 .EQ. Nullp) THEN
         Coord2(:) = Clist(:,Plist(DAT,Slist(DAT,I_SurfP)))
      ELSE
         Coord2(:) = Clist(:,Plist(DAT,I_P2))
      ENDIF

cle>  Coord2(:) = Clist(:,Plist(DAT,I_P2))
      I_SurfQ = Slist(SUCC,I_SurfP)
 7600 CONTINUE
      IF(I_SurfQ .EQ. Nullp) GOTO 7650
      I_Q1=Slist(DAT,I_SurfQ)
      IF (I_Q1 .EQ. Nullp) THEN
        CALL Remove_surf(I_SurfQ,Rootp,M_Surf,Slistp,Slist)
        I_SurfQ = Slist(SUCC,I_SurfQ)
        GOTO 7600
      END IF
      CoordQ1(:) = Clist(:,Plist(DAT,I_Q1))
 7700 CONTINUE
      I_Q2=Plist(SUCC,I_Q1)
      IF (I_Q2 .EQ. Nullp) THEN
        CoordQ2(:) = Clist(:,Plist(DAT,Slist(DAT,I_SurfQ)))
      ELSE
        CoordQ2(:) = Clist(:,Plist(DAT,I_Q2))
      ENDIF
      IF (intersect(Coord1,Coord2,CoordQ1,CoordQ2,Coord)) THEN
        I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
        Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
        I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
        Call Add_coord_after(I_Q1,I_Coord,Slist,M_Pnt,Plistp,Plist)
        IF (debug) THEN
          WRITE(*,*) '7700 intersect surfaces'
          CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
        END IF
      END IF
      IF (I_Q2 .NE. Nullp) THEN
        I_Q1=I_Q2
        CoordQ1 = CoordQ2
        GOTO 7700
      END IF
      I_SurfQ = Slist(SUCC,I_SurfQ)
      GOTO 7600
 7650 CONTINUE
C
      IF (I_P2 .NE. Nullp) THEN
        I_P1=I_P2
        Coord1 = Coord2
        GOTO 7500
      END IF
      I_SurfP = Slist(SUCC,I_SurfP)
      GOTO 7400
 7450 CONTINUE
C
C ----------------------------------------------------------
C
C * Combine all surfaces to one scan list
C
      I_SurfP = Rootp
      I_Scan = New_surf(M_Surf,Root_scan,Slistp,Slist,RSlist,0.0D0)
 7000 CONTINUE
      IF(I_SurfP .EQ. Nullp) GOTO 7050
      I_P1=Slist(DAT,I_SurfP)
      IF(I_P1 .EQ. Nullp) THEN
         I_SurfP = Slist(SUCC,I_SurfP)
         IF(I_SurfP .EQ. Nullp) GOTO 7050
C> ????
C>         GOTO 7400
         GOTO 7000
      END IF
      Coord1(:) = Clist(:,Plist(DAT,I_P1))
      I_Coord = New_coord(M_Coord,Clistp,Clist,Coord1)
      Call Add_coord(I_Scan,I_Coord,Slist,M_Pnt,Plistp,Plist)
      IF (debug) THEN
        WRITE(*,*) '7000 combine'
        CALL Write_surfaces (Root_scan,Slist,Plist,Clist,RSlist)
      END IF
 7100 CONTINUE
      I_P1=Plist(SUCC,I_P1)
      IF (I_P1 .NE. Nullp) THEN
        Coord1(:) = Clist(:,Plist(DAT,I_P1))
        I_Coord = New_coord(M_Coord,Clistp,Clist,Coord1)
        Call Add_coord(I_Scan,I_Coord,Slist,M_Pnt,Plistp,Plist)
        IF (debug) THEN
          WRITE(*,*) '7100 combine'
          CALL Write_surfaces (Root_scan,Slist,Plist,Clist,RSlist)
        END IF
        GOTO 7100
      END IF
      I_SurfP = Slist(SUCC,I_SurfP)
      GOTO 7000
 7050 CONTINUE

      IF (debug) THEN
        WRITE(*,*) 'before sort'
        CALL Write_surfaces (Root_scan,Slist,Plist,Clist,RSlist)
      END IF
C
C ======================================================================
C
C Sort points and remove duplicate x-values
C
 7900 CONTINUE
      I_SurfP = Root_scan
 8000 CONTINUE
      IF(I_SurfP .EQ. Nullp) GOTO 8050
      I_P1=Slist(DAT,I_SurfP)
      IF(I_P1 .EQ. Nullp) THEN
         I_SurfP = Slist(SUCC,I_SurfP)
         GOTO 8000
      END IF
      Coord1(:) = Clist(:,Plist(DAT,I_P1))
 8100 CONTINUE
      I_P2=Plist(SUCC,I_P1)
      IF (I_P2 .NE. Nullp) THEN
        Coord2(:) = Clist(:,Plist(DAT,I_P2))
        IF (Coord2(1) .LT. Coord1(1))THEN
          I_Coord = Plist(DAT,I_P1)
          CALL Remove_point(I_P1,I_SurfP,M_Pnt,Plistp,Plist,Slist)
          Call Add_coord_after(I_P2,I_Coord,Slist,M_Pnt,Plistp,Plist)
          GOTO 7900
        ELSE IF(Coord2(1) - Coord1(1) .LT. eps)THEN
          CALL Remove_point(I_P1,I_SurfP,M_Pnt,Plistp,Plist,Slist)
          GOTO 7900
        END IF
          IF (debug) THEN
            WRITE(*,*) '8100 sort'
            CALL Write_surfaces (Root_scan,Slist,Plist,Clist,RSlist)
          END IF
        Coord1 = Coord2
        I_P1 = I_P2
        GOTO 8100
      END IF
      I_SurfP = Slist(SUCC,I_SurfP)
      GOTO 8000
 8050 CONTINUE
C>
      IF(Debug)THEN
        Write(*,*) 'Scan'
        CALL Write_surfaces (Root_scan,Slist,Plist,Clist,RSlist)
        Write(*,*) 'Rootp before'
        CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
      END IF
C
C ======================================================================
C
C * Start pruning by scanlines
C
      I_Scan = Root_scan
 9400 CONTINUE
      IF(I_Scan .EQ. Nullp) GOTO 9450
      I_PScan=Slist(DAT,I_Scan)
      IF (I_PScan .EQ. Nullp) THEN
        GOTO 9450
      END IF
 9500 CONTINUE
      Xmin = Clist(1,Plist(DAT,I_PScan))
C
      CALL Remove_duplicates (Rootp,
     &  Slistp,Slist,Plistp,Plist,Clist,Clistp,M_Coord,M_Pnt,eps)
C
      I_SurfP = Rootp
      CALL Erase_list(Root_out,
     &  Slistp,Slist,Plistp,Plist,Clistp,M_Coord,M_Pnt,M_Surf)
 9000 CONTINUE
      IF(I_SurfP .EQ. Nullp) GOTO 9050
      I_P1=Slist(DAT,I_SurfP)
      IF (I_P1 .EQ. Nullp) THEN
        I_SurfP = Slist(SUCC,I_SurfP)
        GOTO 9000
      END IF
      Coord1(:) = Clist(:,Plist(DAT,I_P1))
 9100 CONTINUE
      I_P2=Plist(SUCC,I_P1)
      IF (I_P2 .EQ. Nullp) THEN
         Coord2(:) = Clist(:,Plist(DAT,Slist(DAT,I_SurfP)))
      ELSE
         Coord2(:) = Clist(:,Plist(DAT,I_P2))
      ENDIF
C cut Left
      IF( Abs(Coord1(1)-Coord2(1)) .GT. Eps ) THEN
        Lmin = Coord1(1) .LE. Coord2(1)
        R=(Xmin-Coord1(1))/(coord2(1)-Coord1(1))
CC      write(*,*) I_P1, Coord1
CC      write(*,*) I_P2, Coord2
CC      Write(*,*) R, Xmin
        IF (R .GE. 0.0 .AND. R .LE. 1.0) THEN
          Coord1(1)=Xmin+eps
          Coord1(2)=R*(Coord2(2)-Coord1(2))+Coord1(2)
          IF(Lmin)THEN
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coord1)
             Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
             Coord(:)=Coord1(:)-(/2*eps,0.0D0/)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
             Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
          ELSE
             Coord(:)=Coord1(:)-(/2*eps,0.0D0/)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
             Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coord1)
             Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
          END IF
          IF (debug) THEN
            WRITE(*,*) '9100 cut left'
            CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
          END IF
        END IF       
      END IF
      Coord1(:)=Coord2(:)
      I_P1 = I_P2
      IF (I_P1 .NE. Nullp) GOTO 9100
      I_SurfP = Slist(SUCC,I_SurfP)
      GOTO 9000
 9050 CONTINUE
C
C ----------------------------------------------------------
C
C * Cut points to the left of current Xmin
C
      I_SurfP = Rootp
 9200 CONTINUE
      I_Surfout = Nullp
      IF(I_SurfP .EQ. Nullp) GOTO 9250
      I_P1=Slist(DAT,I_SurfP)
 9300 CONTINUE
      IF (I_P1 .EQ. Nullp) THEN
        I_SurfP = Slist(SUCC,I_SurfP)
        GOTO 9200
      END IF
      I_Coord = Plist(DAT,I_P1)
      Coord1(:) = Clist(:,I_Coord)

      IF( Coord1(1) .LT. Xmin ) THEN
C        CALL Remove_coord(I_Coord,M_Coord,Clistp)
        I_PS = Plist(PRE,I_P1)
        CALL Remove_point(I_P1,I_SurfP,M_Pnt,Plistp,Plist,Slist)
        IF( I_Surfout .EQ. Nullp) THEN
          I_Surfout = New_surf(M_Surf,Root_out,Slistp,Slist,RSlist,
     &        RSlist(TRANS,I_SurfP))
        END IF
        Call Add_coord(I_Surfout,I_Coord,Slist,M_Pnt,Plistp,Plist)
        IF (debug) THEN
          WRITE(*,*) '9300 removed left'
          CALL Write_surfaces (Root_out,Slist,Plist,Clist,RSlist)
        END IF
        IF( I_PS .NE. Nullp)THEN
          I_P1 = I_PS
        ELSE
          I_P1=Slist(DAT,I_SurfP)
          GOTO 9300
        END IF
      ELSE
C Keep the point
      END IF
      I_P1 = Plist(SUCC,I_P1)
      IF (I_P1 .NE. Nullp) GOTO 9300
      IF(I_Surfout .NE. Nullp) THEN
        CALL Remove_duplicates (I_Surfout,
     &    Slistp,Slist,Plistp,Plist,Clist,Clistp,M_Coord,M_Pnt,eps)
          IF(Debug)THEN
            Write(*,*) 'removed dupl root_out'
            CALL Write_surfaces (Root_out,Slist,Plist,Clist,RSlist)
          END IF
        CALL Compute_area (I_Surfout,Slist,Plist,Clist,RSlist,eps)
      END IF
      I_SurfP = Slist(SUCC,I_SurfP)
      GOTO 9200
 9250 CONTINUE
C>
      IF(Debug)THEN
        Write(*,*) 'root_out'
        CALL Write_surfaces (Root_out,Slist,Plist,Clist,RSlist)
        Write(*,*) 'rootp'
        CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
      END IF
C
C ======================================================================
C
C * Sort by area
C
      CALL Sort_by_area (Root_out,  Slist, RSlist)
      Trans_w_size = Trans_w_size +
     &    Trans_surfaces(Root_out,Slist,Plist,Clist,RSlist,eps)
C
      IF(Debug)THEN
        Write(*,*) "Trans_w_size:",Trans_w_size
      END IF
      I_PScan=Plist(SUCC,I_PScan)
      IF (I_PScan .NE. Nullp) THEN
        GOTO 9500
      END IF
 9450 CONTINUE
      IF(Debug)THEN
        Write(*,*) Trans_w_size, 1.0 - Trans_w_size / Total_area
      END IF
      Transp_result = 1.0 - Trans_w_size / Total_area
C
      RETURN
      END
C
      INTEGER FUNCTION New_surf(M_Surf,Rootp,Slistp,Slist,RSlist,Transp)
      INTEGER Nullp, SUCC, PRE, DAT, TRANS
      PARAMETER (Nullp=-1, SUCC=1, PRE=2, DAT=3, TRANS=1)
      INTEGER Slistp(*),Slist(3,*),M_Surf,Rootp,Ip
      DOUBLE PRECISION Transp, RSlist(2,*)
      M_Surf=M_Surf+1
      Ip=Slistp(M_Surf)
      IF( Rootp .NE. Nullp) THEN
        Slist(PRE,Rootp) = Ip
      ENDIF
      Slist(DAT,Ip) = NullP
      RSlist(TRANS,Ip) = Transp
      Slist(SUCC,Ip) = Rootp
      Slist(PRE,Ip) = Nullp
      Rootp = Ip
      New_surf = Ip
      RETURN
      END
C
      SUBROUTINE Remove_surf(ISurf,Rootp,M_Surf,Slistp,Slist)
      INTEGER Nullp, SUCC, PRE
      PARAMETER (Nullp=-1, SUCC=1, PRE=2)
      INTEGER Slistp(*),Slist(3,*),M_Surf,ISurf,Rootp
      INTEGER ISs
C
      IF( Rootp .EQ. ISurf)THEN
         Rootp = Slist(SUCC,ISurf)
      END IF
      Slistp(M_Surf)=ISurf
      M_Surf=M_Surf-1
      ISs = Slist(SUCC,ISurf)
      IF (ISs .NE. Nullp) THEN
        Slist(PRE,ISs)=Slist(PRE,ISurf)
      END IF
      ISs = Slist(PRE,ISurf)
      IF (ISs .NE. Nullp) THEN
        Slist(SUCC,ISs)=Slist(SUCC,ISurf)
      END IF
      RETURN
      END
C
      SUBROUTINE Add_coord(I_Surf,I_Coord,Slist,M_Pnt,Plistp,Plist)
      INTEGER Nullp, SUCC, PRE, DAT
      PARAMETER (Nullp=-1, SUCC=1, PRE=2, DAT=3)
      INTEGER I_Surf, Slist(3,*),Plistp(*),Plist(3,*),I_Coord,M_Pnt
      INTEGER Ip, Ip1
      M_Pnt=M_Pnt+1
      Ip=Plistp(M_Pnt)
      Ip1=Slist(DAT,I_Surf)
      IF( Ip1 .NE. Nullp) THEN
        Plist(PRE,Ip1) = Ip
      ENDIF
      IF(Ip .EQ. Nullp) Write(*,*) 'Error in add_coord',M_Pnt
      Slist(DAT,I_Surf)=Ip
      Plist(PRE,Ip)=Nullp
      Plist(SUCC,Ip)=Ip1
      Plist(DAT,IP)=I_Coord
      RETURN
      END
C
      SUBROUTINE Add_coord_after(IC,I_Coord,Slist,M_Pnt,Plistp,Plist)
      INTEGER Nullp, SUCC, PRE, DAT
      PARAMETER (Nullp=-1, SUCC=1, PRE=2, DAT=3)
      INTEGER Slist(3,*),Plistp(*),Plist(3,*), IC, I_Coord, M_Pnt
      INTEGER Ip, Ips
      IF(IC .EQ. Nullp) THEN
        write(*,*) 'Error Null in add_coord_after'
      END IF
      M_Pnt=M_Pnt+1
      Ip=Plistp(M_Pnt)
      Ips=Plist(SUCC,IC)
      Plist(SUCC,Ip)=Ips
      Plist(PRE,Ip)=IC
      IF( Ips .NE. Nullp) THEN
         Plist(PRE,Ips)=Ip
      END IF
      Plist(SUCC,IC)=Ip
      Plist(DAT,IP)=I_Coord
      RETURN
      END
C
      SUBROUTINE Remove_point(IP,IS,M_Pnt,Plistp,Plist,Slist)
      INTEGER Nullp, SUCC, PRE, DAT
      PARAMETER (Nullp=-1, SUCC=1, PRE=2, DAT=3)
      INTEGER Plistp(*),Plist(3,*),M_Pnt,Ip,IS,Slist(3,*)
      INTEGER Ips
      IF(IP .EQ. Nullp) THEN
        write(*,*) 'ERROR Remove_point'
        return
      END IF
      IF( Slist(DAT,IS) .EQ. IP)THEN
         Slist(DAT,IS) = Plist(SUCC,Ip)
      END IF
      Plistp(M_Pnt)=IP
      M_Pnt=M_Pnt-1
      Ips = Plist(SUCC,Ip)
      IF (Ips .NE. Nullp) THEN
        Plist(PRE,Ips)=Plist(PRE,Ip)
      END IF
      Ips = Plist(PRE,Ip)
      IF (Ips .NE. Nullp) THEN
        Plist(SUCC,Ips)=Plist(SUCC,Ip)
      END IF
      RETURN
      END
C
      SUBROUTINE Remove_coord(IC,M_Coord,Clistp)
C      INTEGER Nullp, SUCC, PRE, DAT
C      PARAMETER (Nullp=-1, SUCC=1, PRE=2, DAT=3)
      INTEGER Clistp(*),M_Coord,IC
      Clistp(M_Coord)=IC
      M_Coord=M_Coord-1
      RETURN
      END
C
      INTEGER FUNCTION New_coord(M_Coord,Clistp,Clist,Coords)
C      INTEGER Nullp, SUCC, PRE, DAT
C      PARAMETER (Nullp=-1, SUCC=1, PRE=2, DAT=3)
      DOUBLE PRECISION Coords(2),Clist(2,*)
      INTEGER Clistp(*),M_Coord,Ip
      M_Coord=M_Coord+1
      Ip=Clistp(M_Coord)
      Clist(:,Ip)=Coords(:)
      New_coord = Ip
      RETURN
      END

      DOUBLE PRECISION FUNCTION dist(Coord1, Coord2)
      DOUBLE PRECISION Coord1(2), Coord2(2)
      DOUBLE PRECISION x, y
C
      X=Coord1(1)-Coord2(1)
      Y=Coord1(2)-Coord2(2)
C
      dist = sqrt( x*x + y*y )
      RETURN
      END
      

      LOGICAL FUNCTION  solve2 (A,X,B)
C
C Solve a 2x2 equation system A*X=B
C
C * Inputs
C    A                - Left hand side
C    B                - Right hand side
C
C * Outputs:
C     solve2        - .TRUE. if successfull
C    X                - solution
C
      DOUBLE PRECISION A(2,2), X(2), B(2)
      DOUBLE PRECISION epsilon, SLASK, R
      PARAMETER (epsilon = 1D-10)
      INTEGER I
      LOGICAL LPIV
C
      LPIV = .FALSE.
C
      IF( ABS(A(1,1)) .LE. epsilon)THEN
        IF( ABS(A(2,1)) .LE. epsilon)THEN
          solve2 = .FALSE.
          RETURN
        ELSE
          LPIV = .TRUE.
        END IF
      END IF
C
      IF( ABS(A(1,1)) .LE. ABS(A(2,1)) )THEN
        LPIV = .TRUE.
      END IF
      IF (LPIV) THEN
C Row Pivot
        DO I = 1, 2
             SLASK = A(1,I)
          A(1,I) = A(2,I)
          A(2,I) = SLASK
        END DO
        SLASK = B(1)
        B(1) = B(2)
        B(2) = SLASK
      END IF
C
      R = A(2,1) / A(1,1)
C Not used      A(2,1) = A(2,1) - A(1,1) * R (= 0.0)
      A(2,2) = A(2,2) - A(1,2) * R
      B(2) = B(2) - B(1) * R
      IF( ABS(A(2,2)) .LE. epsilon)THEN
        solve2 = .FALSE.
        RETURN
      END IF
      R = A(1,2) / A(2,2)
C Not used      A(1,2) = A(1,2) - A(2,2) * R (= 0.0)
      B(1) = B(1) - B(2) * R
      X(1) = B(1) / A(1,1)
      X(2) = B(2) / A(2,2)
      solve2 = .TRUE.
      RETURN      
      END
C
      LOGICAL FUNCTION intersect (P1, P2, Q1, Q2, S)
C
C Compute the intersection of the line between P1 and P2 and the
C line Q1 and Q2.
C
C * Inputs
C    P1,P2            - Points of the first line
C    Q1,Q2            - Points of the second line
C
C * Outputs:
C     intersect    - .TRUE. if successfull
C    S                - The point of intersection
C
      DOUBLE PRECISION P1(2), P2(2), Q1(2), Q2(2), S(2)
      DOUBLE PRECISION A(2,2), X(2), B(2), eps
      PARAMETER (eps = 1.0D-8)
      LOGICAL solve2, Debug
C
      Debug = .FALSE.
C
      A(1,1) = P2(1) - P1(1)
      A(1,2) = Q1(1) - Q2(1)
      A(2,1) = P2(2) - P1(2)
      A(2,2) = Q1(2) - Q2(2)
      B(1) = Q1(1) - P1(1)
      B(2) = Q1(2) - P1(2)
      IF (solve2(A,X,B)) THEN
         IF (X(1) .GE. +eps .AND. X(1) .LE. 1.0D0-eps .AND.
     &       X(2) .GE. +eps .AND. X(2) .LE. 1.0D0-eps) THEN
           intersect = .TRUE.
           S(1) = X(1) * (P2(1) - P1(1)) + P1(1)
           S(2) = X(1) * (P2(2) - P1(2)) + P1(2)
           IF(Debug)THEN
             Write(*,*) P1(1), S(1), P2(1)
             Write(*,*) P1(2), S(2), P2(2)
             Write(*,*) Q1(1), X(2) * (Q2(1) - Q1(1)) + Q1(1), Q2(1)
             Write(*,*) Q1(2), X(2) * (Q2(2) - Q1(2)) + Q1(2), Q2(2)
           END IF
         ELSE
           intersect = .FALSE.
         END IF
      ELSE
         intersect = .FALSE.
      ENDIF
      RETURN
      END
C
      SUBROUTINE Remove_duplicates (Rootp,
     &  Slistp,Slist,Plistp,Plist,Clist,Clistp,M_Coord,M_Pnt,eps)
      INTEGER Nullp, SUCC, DAT
      PARAMETER (Nullp=-1, SUCC=1, DAT=3)
      INTEGER Slistp(*),Slist(3,*),Plistp(*),Plist(3,*),Clistp(3,*)
      INTEGER Rootp,M_Coord,M_Pnt
      DOUBLE PRECISION Clist(2,*),eps, dist
C
C Remove duplicate points in all surfaces in the list Rootp
C
      INTEGER I_SurfP, I_P1, I_P2
      DOUBLE PRECISION Coord1(2), Coord2(2)
C
      I_SurfP = Rootp
 4000 CONTINUE
      IF(I_SurfP .EQ. Nullp) GOTO 4050
      I_P1=Slist(DAT,I_SurfP)
      IF (I_P1 .EQ. Nullp) THEN
        GOTO 4300
      END IF 
      Coord1(:) = Clist(:,Plist(DAT,I_P1))
 4100 CONTINUE
      I_P2=Plist(SUCC,I_P1)
      IF (I_P2 .EQ. Nullp) THEN
         Coord2(:) = Clist(:,Plist(DAT,Slist(DAT,I_SurfP)))
      ELSE
         Coord2(:) = Clist(:,Plist(DAT,I_P2))
      ENDIF
C cut Top and Down
      IF( 2*eps .GT. dist(Coord1, Coord2) ) THEN
        CALL Remove_coord(Plist(DAT,I_P1),M_Coord,Clistp)
        CALL Remove_point(I_P1,I_SurfP,M_Pnt,Plistp,Plist,Slist)
      END IF
      Coord1(:)=Coord2(:)
      I_P1 = I_P2
      IF (I_P1 .NE. Nullp) GOTO 4100
 4300 CONTINUE
      I_SurfP = Slist(SUCC,I_SurfP)
      GOTO 4000
 4050 CONTINUE
C
      RETURN      
      END
C
      SUBROUTINE Compute_area (Rootp,
     &  Slist,Plist,Clist,RSlist,eps)
      INTEGER Nullp, SUCC, DAT, AREA
      PARAMETER (Nullp=-1, SUCC=1, DAT=3, AREA=2)
      INTEGER Slist(3,*),Plist(3,*)
      INTEGER Rootp
      DOUBLE PRECISION Clist(2,*),RSlist(2,*),eps
C
C Remove duplicate points in all surfaces integer the list Rootp
C
      INTEGER I_SurfP, I_P1, I_P2, Ic, I
cle>  DOUBLE PRECISION Coord(2,20), Xmin, Xmax, Ymin, Ymax
      DOUBLE PRECISION Coord(2,100), Xmin, Xmax, Ymin, Ymax
C
cle>  Alternative way of computing area. Commented and kept
C     for later use?????
C      INTEGER SurfP
C      DOUBLE PRECISION AB(2), AC(2), ASurf, ATriangle,
C     &  Coord1(2), Coord2(2)
C
C      SurfP = Rootp
CC     Fetch 1st point A
C      I_P1=Slist(DAT,SurfP)
C      IF(I_P1 == NullP) THEN
C        RSlist(AREA,SurfP) = 0.0
C        RETURN
C      END IF
C      Coord1(:) = Clist(:,Plist(DAT,I_P1))
CC     Fetch 2nd point B
C      I_P1=Plist(SUCC,I_P1)
C      IF(I_P1 == NullP) THEN
C        RSlist(AREA,SurfP) = 0.0
C        RETURN
C      END IF
C      Coord2(:) = Clist(:,Plist(DAT,I_P1))
C      AB = Coord2 - Coord1
CC Repeat from here for each C point (nr 3,..)
C      ASurf = 0.0
C 7200 CONTINUE
C      I_P1=Plist(SUCC,I_P1)
C      IF(I_P1 /= NullP) THEN
C        Coord2(:) = Clist(:,Plist(DAT,I_P1))
C        AC = Coord2 - Coord1
C        ATriangle = 0.5 * (AB(1)*AC(2) - AC(1)*AB(2))
C        ASurf = ASurf + Atriangle
C
CC Next point
C        AB = AC
C        GOTO 7200
C      END IF
C      RSlist(AREA,SurfP) = ABS(ASurf)
C      RETURN
cle>
      Xmin=1.0D10
      Xmax=-1.0D10
      Ymin=1.0D10
      Ymax=-1.0D10
C
      I_SurfP = Rootp
 4000 CONTINUE
      Ic=0
      IF(I_SurfP .EQ. Nullp) GOTO 4050
      I_P1=Slist(DAT,I_SurfP)
      IF (I_P1 .EQ. Nullp) THEN
        GOTO 4300
      END IF 
      Ic=Ic+1
      Coord(:,Ic) = Clist(:,Plist(DAT,I_P1))
 4100 CONTINUE
      I_P2=Plist(SUCC,I_P1)
      IF (I_P2 .NE. Nullp) THEN
        Ic=Ic+1
        Coord(:,Ic) = Clist(:,Plist(DAT,I_P2))
      ENDIF
      I_P1 = I_P2
      IF (I_P1 .NE. Nullp) GOTO 4100
 4300 CONTINUE
CC        write(*,*) 'Ic = ',Ic      
      IF (Ic .LE. 2)THEN
        RSlist(AREA,I_SurfP) = 0.0D0
      ELSE IF (Ic .EQ. 3)THEN
        DO I = 1, 3
          Xmin=Min(Xmin,Coord(1,I))
          Xmax=Max(Xmax,Coord(1,I))
          Ymin=Min(Ymin,Coord(2,I))
          Ymax=Max(Ymax,Coord(2,I))
        END DO
CC        write(*,*) xmin,xmax,ymin,ymax
        RSlist(AREA,I_SurfP) = (Xmax-Xmin)*(Ymax-Ymin)/2.0
      ELSE IF (Ic .EQ. 4)THEN
        DO I = 1, 3
          Xmin=Min(Xmin,Coord(1,I))
          Xmax=Max(Xmax,Coord(1,I))
          Ymin=Min(Ymin,Coord(2,I))
          Ymax=Max(Ymax,Coord(2,I))
        END DO
CC        write(*,*) xmin,xmax,ymin,ymax
        RSlist(AREA,I_SurfP) = (Xmax-Xmin)*(Ymax-Ymin)/2.0
        Xmin=1.0D10
        Xmax=-1.0D10
        Ymin=1.0D10
        Ymax=-1.0D10
        DO I = 3, 4
          Xmin=Min(Xmin,Coord(1,I))
          Xmax=Max(Xmax,Coord(1,I))
          Ymin=Min(Ymin,Coord(2,I))
          Ymax=Max(Ymax,Coord(2,I))
        END DO
        I=1
        Xmin=Min(Xmin,Coord(1,I))
        Xmax=Max(Xmax,Coord(1,I))
        Ymin=Min(Ymin,Coord(2,I))
        Ymax=Max(Ymax,Coord(2,I))
CC        write(*,*) xmin,xmax,ymin,ymax
        RSlist(AREA,I_SurfP) = RSlist(AREA,I_SurfP) +
     &     (Xmax-Xmin)*(Ymax-Ymin)/2.0
      ELSE
CC        Write(*,*) 'Error Ic',Ic
        DO I = 1, Ic
CC           Write(*,*) Coord(:,I)
        END DO
      END IF
C      I_SurfP = Slist(SUCC,I_SurfP)
C      GOTO 4000
 4050 CONTINUE
C
      RETURN
      END
C
      SUBROUTINE Write_surfaces (Rootp,
     &  Slist,Plist,Clist,RSlist)
      INTEGER Nullp, SUCC, DAT, AREA
      PARAMETER (Nullp=-1, SUCC=1, DAT=3, AREA=2)
      INTEGER Slist(3,*),Plist(3,*)
      INTEGER Rootp
      DOUBLE PRECISION Clist(2,*),RSlist(2,*)
C
      INTEGER I_SurfP, I_P1
      DOUBLE PRECISION Coord1(2)
C
      I_SurfP = Rootp
 7200 CONTINUE
      IF(I_SurfP .EQ. Nullp) GOTO 7250
      Write(*,*) 'Area:',RSlist(AREA,I_SurfP)
      I_P1=Slist(DAT,I_SurfP)
      if(i_p1 .ne. nullp) then
        Coord1(:) = Clist(:,Plist(DAT,I_P1))
        Write(*,*) I_SurfP, I_P1, Coord1
      end if
 7300 CONTINUE
      IF (I_P1 .NE. Nullp) THEN
         I_P1=Plist(SUCC,I_P1)
         IF (I_P1 .NE. Nullp) THEN
            Coord1(:) = Clist(:,Plist(DAT,I_P1))
            Write(*,*) I_SurfP, I_P1, Coord1
            GOTO 7300
         END IF
      END IF
      I_SurfP = Slist(SUCC,I_SurfP)
      GOTO 7200
 7250 CONTINUE
      RETURN
      END
C
      SUBROUTINE Erase_list (Rootp,
     &  Slistp,Slist,Plistp,Plist,Clistp,M_Coord,M_Pnt,M_Surf)
      INTEGER Nullp, SUCC, DAT
      PARAMETER (Nullp=-1, SUCC=1, DAT=3)
      INTEGER Slistp(*),Slist(3,*),Plistp(*),Plist(3,*),Clistp(3,*)
      INTEGER Rootp,M_Coord,M_Pnt,M_Surf
C
      INTEGER I_SurfP, I_P1, I_P2
C
      I_SurfP = Rootp
 7200 CONTINUE
      IF(I_SurfP .EQ. Nullp) GOTO 7250
      I_P1=Slist(DAT,I_SurfP)
 7300 CONTINUE
      IF (I_P1 .NE. Nullp) THEN
        I_P2=Plist(SUCC,I_P1)
        CALL Remove_coord(Plist(DAT,I_P1),M_Coord,Clistp)
        CALL Remove_point(I_P1,I_SurfP,M_Pnt,Plistp,Plist,Slist)
        I_P1=I_P2
        GOTO 7300
      END IF
      CALL Remove_surf(I_SurfP,Rootp,M_Surf,Slistp,Slist)
      I_SurfP = Slist(SUCC,I_SurfP)
      GOTO 7200
 7250 CONTINUE
      RETURN
      END
C
      SUBROUTINE Sort_by_area (Rootp,  Slist, RSlist)
      INTEGER Nullp, SUCC, PRE, AREA
      PARAMETER (Nullp=-1, SUCC=1, PRE=2, AREA=2)
      INTEGER Rootp, Slist(3,*)
      DOUBLE PRECISION RSlist(2,*)
C
      INTEGER I_S1, I_S2
 7900 CONTINUE
      I_S1 = Rootp
 8000 CONTINUE
      IF(I_S1 .EQ. Nullp) GOTO 8050
      I_S2 = Slist(SUCC,I_S1)
      IF(I_S2 .EQ. Nullp) GOTO 8050
      IF( RSlist(AREA,I_S1) .GT. RSlist(AREA,I_S2))THEN
        IF( I_S1 .EQ. Rootp)THEN
          Rootp = I_S2
        END IF
        IF( Slist(SUCC,I_S2) .NE. Nullp) THEN
          Slist(PRE,Slist(SUCC,I_S2)) = I_S1
          Slist(SUCC,I_S1) = Slist(SUCC,I_S2)
        ELSE
          Slist(SUCC,I_S1) = Nullp
        END IF
        IF( Slist(PRE,I_S1) .NE. Nullp) THEN
          Slist(SUCC,Slist(PRE,I_S1)) = I_S2
          Slist(PRE,I_S2) = Slist(PRE,I_S1)
        ELSE
          Slist(PRE,I_S2) = Nullp
        END IF
        Slist(SUCC,I_S2) = I_S1
        Slist(PRE,I_S1) = I_S2
        GOTO 7900
      END IF
      I_S1 = I_S2
      GOTO 8000
 8050 CONTINUE
      RETURN
      END
C
      DOUBLE PRECISION FUNCTION Trans_surfaces (Rootp,
     &  Slist,Plist,Clist,RSlist,eps)
      INTEGER Nullp, SUCC, DAT
      PARAMETER (Nullp=-1, SUCC=1, DAT=3)
      INTEGER Slist(3,*),Plist(3,*)
      INTEGER Rootp
      DOUBLE PRECISION Clist(2,*),RSlist(2,*),eps
C
      INTEGER I_SurfP
      DOUBLE PRECISION Tloc, Xmin, Xmax, Y1min, Y1max, Y2min, Y2max
C>ML:990108
      DOUBLE PRECISION Xmin1, Xmax1
C>
cle>  DOUBLE PRECISION Coord(2,20)
      DOUBLE PRECISION Coord(2,100)
      INTEGER IC, I_P1, I_P2, I, Icmax
      INTEGER N_points
      LOGICAL UNDEFX,UNDEF1,UNDEF2
C
      UNDEFX=.TRUE.
      UNDEF1=.TRUE.
      UNDEF2=.TRUE.
      Xmin=1.0D10
      Xmax=-1.0D10
      Y1min=1.0D10
      Y2min=1.0D10
      Y1max=-1.0D10
      Y2max=-1.0D10
      I_SurfP = Rootp
      Tloc = 0.0D0
      Icmax=0
C      
 7200 CONTINUE
      Ic=0
      IF(I_SurfP .EQ. Nullp) GOTO 7250
      IF(N_points(I_Surfp,Slist,Plist) .LE. 2)THEN
        I_SurfP = Slist(SUCC,I_SurfP)
        GOTO 7200
      END IF
      I_P1=Slist(DAT,I_SurfP)
      IF (I_P1 .EQ. Nullp) THEN
        GOTO 4300
      END IF 
      Ic=Ic+1
      Coord(:,Ic) = Clist(:,Plist(DAT,I_P1))
 4100 CONTINUE
CC      write(*,*) 'I_P1',I_P1
      I_P2=Plist(SUCC,I_P1)
      IF (I_P2 .NE. Nullp) THEN
        Ic=Ic+1
        Coord(:,Ic) = Clist(:,Plist(DAT,I_P2))
      ENDIF
      I_P1 = I_P2
      IF (I_P1 .NE. Nullp) GOTO 4100
 4300 CONTINUE
      Icmax= Max(Ic,Icmax)
C>ML:990108
      Xmin1=1.0D10
      Xmax1=-1.0D10
      DO I=1,IC
        Xmin1=Min(Xmin1,Coord(1,I))
        Xmax1=Max(Xmax1,Coord(1,I))
      END DO
      IF(IC .GT. 2 .AND. Xmax1-Xmin1 > 2*eps) THEN
C>
        DO I=1,IC
          UNDEFX=.FALSE.
          Xmin=Min(Xmin,Coord(1,I))
          Xmax=Max(Xmax,Coord(1,I))
        END DO
        DO I=1,IC
          IF(Abs(Coord(1,I)-Xmin) .LE. 2*eps)THEN
            UNDEF1=.FALSE.
            Y1min=Min(Y1min,Coord(2,I))
            Y1max=Max(Y1max,Coord(2,I))
          ELSE
            UNDEF2=.FALSE.
            Y2min=Min(Y2min,Coord(2,I))
            Y2max=Max(Y2max,Coord(2,I))
          END IF
        END DO
C>ML:990108
      END IF
C>
C>      Tloc = Tloc +
CA>    & RSlist(AREA,I_SurfP) * (1.0 - RSlist(TRANS,I_SurfP))
      I_SurfP = Slist(SUCC,I_SurfP)
      GOTO 7200
 7250 CONTINUE
      IF(Icmax .LE. 2) THEN
        UNDEFX = .TRUE.
        UNDEF1 = .TRUE.
        UNDEF2 = .TRUE.
      END IF
      IF(.NOT. (UNDEFX .OR. UNDEF1 .OR. UNDEF2))THEN
CC        write(*,*) 'icmax',icmax
CC        write(*,*) xmax,xmin
CC        write(*,*) y1max,y1min
CC        write(*,*) y2max,y2min
        Tloc = (Y1Max - Y1Min) * (Xmax - Xmin) / 2.0
        Tloc = Tloc + (Y2Max - Y2Min) * (Xmax - Xmin) / 2.0
      END IF
      Trans_surfaces = Tloc
      RETURN
      END

      INTEGER FUNCTION  N_Points (I_Surfp,Slist,Plist)
      INTEGER Nullp, SUCC, DAT
      PARAMETER (Nullp=-1, SUCC=1, DAT=3)
      INTEGER Slist(3,*),Plist(3,*)
C
      INTEGER I_SurfP
      INTEGER IC, I_P1
C
 7200 CONTINUE
      Ic=0
      IF(I_SurfP .EQ. Nullp) GOTO 7250
      I_P1=Slist(DAT,I_SurfP)
      IF (I_P1 .EQ. Nullp) THEN
        GOTO 7250
      END IF 
 4100 CONTINUE
      Ic=Ic+1
      I_P1=Plist(SUCC,I_P1)
      IF (I_P1 .NE. Nullp) THEN
        GOTO 4100
      ENDIF
 7250 CONTINUE
      N_Points = Ic
      RETURN
      END

      DOUBLE PRECISION FUNCTION Scalar_product(X,Y)
      DOUBLE PRECISION X(*),Y(*)
C
C
      Scalar_product=X(1)*Y(1)+X(2)*Y(2)+X(3)*Y(3)
      RETURN
      END

      SUBROUTINE Cross_product(res,vec1,vec2)
      DOUBLE PRECISION res(*), vec1(*), vec2(*)
      res(1)=vec1(2)*vec2(3) - vec2(2)*vec1(3)
      res(2)=vec1(3)*vec2(1) - vec2(3)*vec1(1)
      res(3)=vec1(1)*vec2(2) - vec2(1)*vec1(2)
      RETURN
      END

      SUBROUTINE CutZZero(C1,C2)
      DOUBLE PRECISION C1(3), C2(3)
      DOUBLE PRECISION R
C * Cut the line C1C2 with Z=0 Return the part that is positive Z
      IF(C1(3) .LE. 0)THEN
         IF(C2(3) .LE. 0)THEN
C     * Do nothing
         ELSE 
            R=C1(3)/(C2(3)-C1(3))
            C1=C1-R*(C2-C1)
         END IF
      ELSE IF(C2(3) .LE. 0)THEN
         R=C2(3)/(C1(3)-C2(3))
         C2=C2-R*(C1-C2)
      END IF
      RETURN
      END