C CADRCL    SOURCE    PV090527  25/04/10    21:15:01     12230          
C  TRACE DE DEFORMES
C  CALCUL DU CADRE GLOBAL
C  PROJECTION SUCCESSIVE DE CHAQUE DEFORME (SELON ICLE)
C  1995 option DIRE P.PEGON JRC-ISPRA
C
C PP option DIRE
      SUBROUTINE CADRCL(KABCOR,LABCO2,IOEIL,XPROJ,
     #  ICLE,XMINT,YMINT,XMAXT,YMAXT,zmint,zmaxt,cgrav,diloc,LDIRE,
     >    axez)

      IMPLICIT INTEGER(I-N)
***   IMPLICIT REAL*8(A-H,O-Z)

-INC PPARAM
-INC CCOPTIO
-INC CCREEL
-INC SMCOORD

      dimension cgrav(*),axez(*)
C+PP option DIRE
      dimension diloc(*)
      LOGICAL LDIRE
C+PP
      DIMENSION XO(3),XP(3),XN(3),XG(3),XV(3),UI(3),UJ(3)
      COMMON /CCADRC/XN,XG,XO,UI,UJ
C  ATTENTION
      SEGMENT KABCOR(0)
      SEGMENT KABCO2(2,0)
      SEGMENT LABCO2(3,0)
      SEGMENT ICOR2(0)
      SEGMENT KXPRO2(NVEC)
      SEGMENT SXCORD
        REAL   XCORD(IDIM,ITE)
      ENDSEGMENT
      SEGMENT SXCOR2
        REAL   XCOR2(IDIM,ITE)
      ENDSEGMENT
      SEGMENT XPROJ(3,ITE)
      SEGMENT XPRO2(3,ITE)
C+PP option DIRE
      REAL*8 SSN

*dbg      write(6,*) 'coucou cadrcl' 
C+PP
C  ICLE = 0 RECHERCHE DES MIN MAX SUR TOUTES LES DEFORMES
C  CALCUL DE LA PROJECTION
      XMINT=1E30
      YMINT=XMINT
      ZMINT=XMINT
      XMAXT=-1E30
      YMAXT=XMAXT
      ZMAXT=XMAXT
      NCORD=KABCOR(/1)

      IF (IDIM.EQ.3) GOTO 50
      
      ZMINT=0.
      ZMAXT=0.
      DO ICORD=1,NCORD
         SXCORD=KABCOR(ICORD)
*        write(6,*) '1 cadrcl icord sxcord',icord,sxcord
         ITE=XCORD(/2)
         DO J=1,ITE
            XMINT=MIN(XMINT,XCORD(1,J))
            XMAXT=MAX(XMAXT,XCORD(1,J))
            YMINT=MIN(YMINT,XCORD(2,J))
            YMAXT=MAX(YMAXT,XCORD(2,J))
         ENDDO
      ENDDO
      
      IF (ICLE.NE.0) GOTO 1000
      IF (LABCO2.EQ.0) GOTO 60
      DO ICORD=1,NCORD
        IF (LABCO2(3,ICORD).EQ.0) GOTO 61
        KABCO2=LABCO2(1,ICORD)
        SXCORD=KABCOR(ICORD)
*        write(6,*) '2 cadrcl icord sxcord',icord,sxcord
        ITE=XCORD(/2)
        NVEC=KABCO2(/2)
        IF (NVEC.EQ.0) GOTO 61
        DO IVEC=1,NVEC
          SXCOR2=KABCO2(1,IVEC)
          ICOR2=KABCO2(2,IVEC)
          ITE=XCOR2(/2)
          DO J=1,ITE
            IF (ICOR2(J).EQ.0) GOTO 18
            UX=XCOR2(1,J)-XCORD(1,J)
            UY=XCOR2(2,J)-XCORD(2,J)
            U1=XCOR2(1,J)-UX/3-UY/5
            V1=XCOR2(2,J)-UY/3+UX/5
            XMINT=MIN(XMINT,U1)
            XMAXT=MAX(XMAXT,U1)
            YMINT=MIN(YMINT,V1)
            YMAXT=MAX(YMAXT,V1)
            U1=XCOR2(1,J)-UX/3+UY/5
            V1=XCOR2(2,J)-UY/3-UX/5
            XMINT=MIN(XMINT,U1)
            XMAXT=MAX(XMAXT,U1)
            YMINT=MIN(YMINT,V1)
            YMAXT=MAX(YMAXT,V1)
            XMINT=MIN(XMINT,XCOR2(1,J))
            XMAXT=MAX(XMAXT,XCOR2(1,J))
            YMINT=MIN(YMINT,XCOR2(2,J))
            YMAXT=MAX(YMAXT,XCOR2(2,J))
          ENDDO
  18      CONTINUE
        ENDDO
  61    CONTINUE
      ENDDO
  60  CONTINUE
*  pour eviter des calculs sur des quantites non significatives,
*  on agrandit un peu le cadre
      XDIFF=XMAXT-XMINT
      YDIFF=YMAXT-YMINT
c*      ZDIFF=ZMAXT-ZMINT
      ZDIFF=0.
      XSOMM=XMAXT+XMINT
      YSOMM=YMAXT+YMINT
c*      ZSOMM=ZMAXT+ZMINT
      ZSOMM=0.
      XDIFF=MAX(XDIFF,YDIFF/10,ZDIFF/10)
      YDIFF=MAX(XDIFF/10,YDIFF,ZDIFF/10)
      ZDIFF=MAX(XDIFF/10,YDIFF/10,ZDIFF)
      XMINT=(XSOMM-XDIFF)/2
      XMAXT=(XSOMM+XDIFF)/2
      YMINT=(YSOMM-YDIFF)/2
      YMAXT=(YSOMM+YDIFF)/2
      ZMINT=(ZSOMM-ZDIFF)/2
      ZMAXT=(ZSOMM+ZDIFF)/2
*     write(6,*) '1 cadrcl ',xmint,xmaxt,ymint,ymaxt,zmint,zmaxt
      IF (ICLE.NE.0) GOTO 1000
      RETURN

C* Cas DIME = 3 :
  50  CONTINUE
      SEGACT,MCOORD
      IREF=(IOEIL-1)*4
      XO(1)=XCOOR(IREF+1)
      XO(2)=XCOOR(IREF+2)
      XO(3)=XCOOR(IREF+3)
C+PP option DIRE
      IF (LDIRE)THEN
         DO J=1,3
            XG(J)=cgrav(J)
            XN(J)=XO(J)-XG(J)
         ENDDO
C+PP option DIRE
C+PP
      ELSE
C POINT MOYEN
        XG(1)=0.D0
        XG(2)=0.D0
        XG(3)=0.D0
        NPTOT=0
        DO ICORD=1,NCORD
          SXCORD=KABCOR(ICORD)
*      write(6,*) '3 cadrcl icord sxcord',icord,sxcord
          ITE=XCORD(/2)
          NPTOT=NPTOT+ITE
          DO I=1,ITE
            XG(1)=XG(1)+XCORD(1,I)
            XG(2)=XG(2)+XCORD(2,I)
            XG(3)=XG(3)+XCORD(3,I)
          ENDDO
        ENDDO
        DO J=1,3
          XG(J)=XG(J)/NPTOT
          XN(J)=XO(J)-XG(J)
          cgrav(j)=XG(j)
        ENDDO
      ENDIF
C+PP
C  NORMALE AU PLAN
      SNI=MAX(ABS(XN(1)),ABS(XN(2)),ABS(XN(3)))+XSPETI
      SN=SQRT((XN(1)/SNI)**2+(XN(2)/SNI)**2+(XN(3)/SNI)**2)
      SN=SN*SNI
      IF (SN.EQ.0.) THEN
        CALL ERREUR(21)
        RETURN
      ENDIF
C+PP option DIRE
      SSN=SN
C+PP
      XN(1)=XN(1)/SN
      XN(2)=XN(2)/SN
      XN(3)=XN(3)/SN

C  REPERE LOCAL SUR LE PLAN
C+PP option DIRE
      IF (LDIRE)THEN
         UJ(1)=diloc(1)
         UJ(2)=diloc(2)
         UJ(3)=diloc(3)
      ELSE
C+PP
         UJ(1)=0.D0
         UJ(2)=0.D0
         UJ(3)=1.D0
C+PP option DIRE
      ENDIF
C+PP
  21  CONTINUE
      SV=UJ(1)*XN(1)+UJ(2)*XN(2)+UJ(3)*XN(3)
      DO 20 J=1,3
         UJ(J)=UJ(J)-SV*XN(J)
  20  CONTINUE
      SV=UJ(1)**2+UJ(2)**2+UJ(3)**2
      IF (ABS(SV).LT.0.1) THEN
         UJ(1)=0.D0
         UJ(2)=1.D0
         UJ(3)=1.D0
         GOTO 21
      ENDIF
      SV=SQRT(SV)
      UJ(1)=UJ(1)/SV
      UJ(2)=UJ(2)/SV
      UJ(3)=UJ(3)/SV
      UI(1)=UJ(2)*XN(3)-UJ(3)*XN(2)
      UI(2)=UJ(3)*XN(1)-UJ(1)*XN(3)
      UI(3)=UJ(1)*XN(2)-UJ(2)*XN(1)
C PROJECTION CONIQUE SUR LE PLAN
      DO 170 ICORD=1,NCORD
         SXCORD=KABCOR(ICORD)
*        write(6,*) '4 cadrcl icord sxcord',icord,sxcord
         ITE=XCORD(/2)
         DO 12 I=1,ITE
            DO J=1,3
               XP(J)=XCORD(J,I)
               XV(J)=XP(J)-XO(J)
            ENDDO
            SV=XV(1)*XN(1)+XV(2)*XN(2)+XV(3)*XN(3)
C+PP option DIRE
            IF ((LDIRE.AND.-SV.GE.SSN).OR.(.NOT.LDIRE))THEN
C+PP
               SN=(XP(1)-XG(1))*XN(1)+(XP(2)-XG(2))*XN(2)+(XP(3)-XG(3))
     $              *XN(3)
               DO 14 J=1,3
                  XP(J)=XP(J)-(SN/SV)*XV(J)-XG(J)
 14            CONTINUE
               XPRO=XP(1)*UI(1)+XP(2)*UI(2)+XP(3)*UI(3)
               YPRO=XP(1)*UJ(1)+XP(2)*UJ(2)+XP(3)*UJ(3)
               ZPRO=-SN
               XMINT=MIN(XMINT,XPRO)
               XMAXT=MAX(XMAXT,XPRO)
               YMINT=MIN(YMINT,YPRO)
               YMAXT=MAX(YMAXT,YPRO)
               ZMINT=MIN(ZMINT,ZPRO)
               ZMAXT=MAX(ZMAXT,ZPRO)
C+PP option DIRE
            ENDIF
C+PP
 12      CONTINUE
 170  CONTINUE
*     write(6,*) 'cadrcl en 239 ',xmint,xmaxt,ymint,ymaxt,zmint,zmaxt
      IF (LABCO2.EQ.0) GOTO 11
      DO 171 ICORD=1,NCORD
         IF (LABCO2(3,ICORD).EQ.0) GOTO 171
         KABCO2=LABCO2(1,ICORD)
         SXCORD=KABCOR(ICORD)
*        write(6,*) '5 cadrcl icord sxcord',icord,sxcord
         ITE=XCORD(/2)
         NVEC=KABCO2(/2)
         IF (NVEC.EQ.0) GOTO 171
         DO 6 IVEC=1,NVEC
            SXCOR2=KABCO2(1,IVEC)
            ICOR2=KABCO2(2,IVEC)
            if (icor2.le.0) goto 6
            DO 17 I=1,ITE
               IF(I.GT.ICOR2(/1)) goto 17
               IF (ICOR2(I).EQ.0) GOTO 17
C+PP a faire meme ss DIRE car sinon que signifie XPRO,YPRO?
               DO J=1,3
                  XP(J)=XCORD(J,I)
                  XV(J)=XP(J)-XO(J)
               ENDDO
               SV=XV(1)*XN(1)+XV(2)*XN(2)+XV(3)*XN(3)
C+PP option DIRE
               IF (LDIRE.AND.-SV.LT.SSN) THEN
                  ICOR2(I)=0
                  GOTO 17
               ENDIF
C+PP a faire meme ss DIRE car sinon que signifie XPRO,YPRO?
               SN=(XP(1)-XG(1))*XN(1)+(XP(2)-XG(2))*XN(2)+(XP(3)-XG(3))
     $              *XN(3)
               DO J=1,3
                  XP(J)=XP(J)-(SN/SV)*XV(J)-XG(J)
               ENDDO
               XPRO=XP(1)*UI(1)+XP(2)*UI(2)+XP(3)*UI(3)
               YPRO=XP(1)*UJ(1)+XP(2)*UJ(2)+XP(3)*UJ(3)
C+PP
               DO 15 J=1,3
                  XP(J)=XCOR2(J,I)
                  XV(J)=XP(J)-XO(J)
 15            CONTINUE
               SV=XV(1)*XN(1)+XV(2)*XN(2)+XV(3)*XN(3)
C+PP option DIRE
               IF (LDIRE.AND.-SV.LT.SSN) THEN
                  ICOR2(I)=0
                  GOTO 17
               ENDIF
C+PP
               SN=(XP(1)-XG(1))*XN(1)+(XP(2)-XG(2))*XN(2)+(XP(3)-XG(3))
     $              *XN(3)
               DO 16 J=1,3
                  XP(J)=XP(J)-(SN/SV)*XV(J)-XG(J)
 16            CONTINUE
               XPROO=XP(1)*UI(1)+XP(2)*UI(2)+XP(3)*UI(3)
               YPROO=XP(1)*UJ(1)+XP(2)*UJ(2)+XP(3)*UJ(3)
               XMINT=MIN(XMINT,XPROO)
               XMAXT=MAX(XMAXT,XPROO)
               YMINT=MIN(YMINT,YPROO)
               YMAXT=MAX(YMAXT,YPROO)
               UX=XPROO-XPRO
               UY=YPROO-YPRO
               U1=XPROO-UX/3-UY/5
               V1=YPROO-UY/3+UX/5
               XMINT=MIN(XMINT,U1)
               XMAXT=MAX(XMAXT,V1)
               U1=XPROO-UX/3+UY/5
               V1=YPROO-UY/3-UX/5
               XMINT=MIN(XMINT,U1)
               XMAXT=MAX(XMAXT,V1)
 17         CONTINUE
 6       CONTINUE
 171  CONTINUE
 11   CONTINUE
*  pour eviter des calculs sur des quantites non significatives,
*  on agrandi un peu le cadre
*     write(6,*) '2 cadrcl ',xmint,xmaxt,ymint,ymaxt,zmint,zmaxt
      XDIFF=XMAXT-XMINT
      YDIFF=YMAXT-YMINT
      ZDIFF=ZMAXT-ZMINT
      XSOMM=XMAXT+XMINT
      YSOMM=YMAXT+YMINT
      ZSOMM=ZMAXT+ZMINT
*     write(6,*) '4 cadrcl ',xdiff,ydiff,zdiff,xsomm,ysomm,zsomm
      XDIFF=MAX(XDIFF,YDIFF/10,ZDIFF/10)
      YDIFF=MAX(XDIFF/10,YDIFF,ZDIFF/10)
      ZDIFF=MAX(XDIFF/10,YDIFF/10,ZDIFF)
*     write(6,*) '5 cadrcl ',xdiff,ydiff,zdiff
      XMINT=(XSOMM-XDIFF)/2
      XMAXT=(XSOMM+XDIFF)/2
      YMINT=(YSOMM-YDIFF)/2
      YMAXT=(YSOMM+YDIFF)/2
      ZMINT=(ZSOMM-ZDIFF)/2
      ZMAXT=(ZSOMM+ZDIFF)/2
*     write(6,*) '3 cadrcl ',xmint,xmaxt,ymint,ymaxt,zmint,zmaxt
*  axez pour tourner correctement avec opengl
      axez(1)=0
      axez(2)=uj(3)
      axez(3)=sqrt(1-uj(3)**2)
      if (xn(3).lt.0) axez(3)=-axez(3)
*     write (6,*) ' axez ',axez(1),axez(2),axez(3)
      IF (ICLE.NE.0) GOTO 1000
      RETURN
 1000 CONTINUE
C  ON REMPLIT XPROJ POUR LA DEFORME CONCERNEE
      SXCORD=KABCOR(ICLE)
*        write(6,*) '6 cadrcl icle  sxcord',icord,sxcord
      ITE=XCORD(/2)
      IF (IDIM.EQ.2) GOTO 1100
C+PP option DIRE
      SNI=MAX(ABS(XN(1)),ABS(XN(2)),ABS(XN(3)))+XSPETI
      SN=SQRT((XN(1)/SNI)**2+(XN(2)/SNI)**2+(XN(3)/SNI)**2)
      SN=SN*SNI
      SSN=0.95*SSN
C+PP
C PROJECTION CONIQUE SUR LE PLAN
      DO 612 I=1,ITE
         DO 613 J=1,3
            XP(J)=XCORD(J,I)
            XV(J)=XP(J)-XO(J)
 613     CONTINUE
         SV=XV(1)*XN(1)+XV(2)*XN(2)+XV(3)*XN(3)
         SN=(XP(1)-XG(1))*XN(1)+(XP(2)-XG(2))*XN(2)+(XP(3)-XG(3))*XN(3)
         XPROJ(3,I)=-SN
         DO 614 J=1,3
            XP(J)=XP(J)-(SN/SV)*XV(J)-XG(J)
 614     CONTINUE
C+PP option DIRE
         IF ((LDIRE.AND.-SV.GE.SSN).OR.(.NOT.LDIRE))THEN
C+PP
            XPROJ(1,I)=XP(1)*UI(1)+XP(2)*UI(2)+XP(3)*UI(3)
            XPROJ(2,I)=XP(1)*UJ(1)+XP(2)*UJ(2)+XP(3)*UJ(3)
C+PP option DIRE
         ELSE
            XPROJ(1,I)=(XMINT+XMAXT)/2
            XPROJ(2,I)=(YMINT+YMAXT)/2
         ENDIF
C+PP
 612  CONTINUE
      IF (LABCO2.EQ.0) GOTO 618
      IF (LABCO2(3,ICLE).EQ.0) GOTO 618
      KABCO2=LABCO2(1,ICLE)
      NVEC=KABCO2(/2)
      SEGINI KXPRO2
      LABCO2(2,ICLE)=KXPRO2
      IF (NVEC.EQ.0) GOTO 618
      DO 8 IVEC=1,NVEC
         SXCOR2=KABCO2(1,IVEC)
         ICOR2=KABCO2(2,IVEC)
*     ajout SG 2016/07/16 apparemment le ITE de XCOR2 n'est pas
*     forcement celui de XCORD en cas de coupe
         ITE=XCOR2(/2)
         SEGINI XPRO2
         KXPRO2(IVEC)=XPRO2
         DO 617 I=1,ITE
            IF (ICOR2(I).EQ.0) GOTO 617
            DO 615 J=1,3
               XP(J)=XCOR2(J,I)
               XV(J)=XP(J)-XO(J)
 615        CONTINUE
            SV=XV(1)*XN(1)+XV(2)*XN(2)+XV(3)*XN(3)
            SN=(XP(1)-XG(1))*XN(1)+(XP(2)-XG(2))*XN(2)+(XP(3)-XG(3))
     $           *XN(3)
            XPRO2(3,I)=-SN
            DO 616 J=1,3
               XP(J)=XP(J)-(SN/SV)*XV(J)-XG(J)
 616        CONTINUE
*           write(6,*) ' cadrcl on passe la 353 ',i
            XPRO2(1,I)=XP(1)*UI(1)+XP(2)*UI(2)+XP(3)*UI(3)
            XPRO2(2,I)=XP(1)*UJ(1)+XP(2)*UJ(2)+XP(3)*UJ(3)
 617     CONTINUE
**       SEGSUP SXCOR2
 8    CONTINUE
 618  CONTINUE
**    SEGSUP SXCORD
*     write(6,*) '6 cadrcl ',xmint,xmaxt,ymint,ymaxt,zmint,zmaxt
      RETURN
 1100 CONTINUE
      DO J=1,ITE
        DO I=1,IDIM
          XPROJ(I,J)=XCORD(I,J)
        ENDDO
      ENDDO
      IF (LABCO2.EQ.0) GOTO 1111
      IF (LABCO2(3,ICLE).EQ.0) GOTO 1111
      KABCO2=LABCO2(1,ICLE)
      NVEC=KABCO2(/2)
      SEGINI KXPRO2
      LABCO2(2,ICLE)=KXPRO2
      IF (NVEC.EQ.0) GOTO 1111
      DO 9 IVEC=1,NVEC
      SXCOR2=KABCO2(1,IVEC)
      ICOR2=KABCO2(2,IVEC)
      segact sxcor2
      ITE=XCOR2(/2)
      SEGINI XPRO2
      KXPRO2(IVEC)=XPRO2
*     write(6,*) ' cadrcl on passe aussi la 381',j
      DO 1113 J=1,ITE
        IF (ICOR2(J).EQ.0) GOTO 1113
        DO I=1,IDIM
          XPRO2(I,J)=XCOR2(I,J)
        ENDDO
*     write(6,*) j,xcor2,(xcor2(i,j),i=1,idim)
 1113 CONTINUE
**    SEGSUP SXCOR2
    9 CONTINUE
 1111 CONTINUE
**    SEGSUP SXCORD
*  axez pour tourner correctement avec opengl
      axez(1)=0
      axez(2)=uj(3)
      axez(3)=sqrt(1-uj(3)**2)
      if (xn(3).lt.0) axez(3)=-axez(3)
*     write (6,*) ' axez ',axez(1),axez(2),axez(3)
*     write(6,*) '7 cadrcl ',xmint,xmaxt,ymint,ymaxt,zmint,zmaxt

c      RETURN
      END