C**********************************************************************C
      SUBROUTINE SOLVECZR

C     Provides roots by integration on a contour in the complex plane
C     Use of a circle new version for the number of modes
C     rectangles are used as contours
C     Can split intervals if there is too many modes
    
      INCLUDE 'comc1.f'

      PARAMETER ( MXCER = 20000, MAXSSP = 2001, MXINT = 100 )

      INTEGER*4 NINTER, MXINTMOD, NVITINT, NVBP, NT, NCER, NCERM
      REAL*8 CWAT, CMINI, HEAU, HENS, CIMX, TAILLCER, CSMIN
      REAL*8 KBP( MXINT ), VBP( MXINT ), CBP( MXINT ), OMEGA, CMAXI
      REAL*8 BPC( 5 ), U, V, DR( 2 ), ALPHAM, CL, CH
      REAL*8 VITINT( MAXM ), TMN, TMX, DTM, RAY, CREF
      REAL*8 ALPHAMAX, CENTR, ANG( MXCER ), KBPINC, KBP1, KBP2
      COMPLEX*16 CENT, SOL( MAXM )
      COMPLEX*16 ALP( MAXSSP ), BET( MAXSSP ), RH( MAXSSP )
      COMPLEX*16 POINT( MXINT ), PT, SOL2
      INTEGER*4 NBZ( MXINT ), NBZP, J, I, IPOW, IPOWM, IPOW1, NMODE
      INTEGER*4 REP( 5 ), NRECT, NMP, NCER1, MLT, KK, NCIR, ZINT
      INTEGER NCEL, NBPC
      LOGICAL FLAGRES
      CHARACTER*80 ERRMSG
      
      COMMON /DONNEES/ ALP, BET, RH
      
      NCEL = 2 * NMEDIA
      OMEGA = 2.0 * PI * FREQ
   
      CALL  GETMINFO( CWAT, CMINI, HEAU, HENS, CIMX, CSMIN, CMAXI)   
   
C     Number of points taken to discretize the contour
C     minimum for one side of the rectangle

      NCER = 500
      
C     Number of points taken to discretize a circle

      NCIR = 500

C     Maximum number of modes allowed

      MXINTMOD = 6

      DO K = 1, NCIR  
         ANG( K ) = 2.0D0 * PI * DBLE( K ) / DBLE( NCIR )
      END DO

C     Look for the total number of modes
C     and split all the intervals

C     First check for the branch points

      BPC( 1 ) = CICR( CPB )
      BPC( 2 ) = CICR( CSB )
      BPC( 3 ) = CWAT
      BPC( 4 ) = CICR( CPT )
      BPC( 5 ) = CICR( CST )
      NBPC = 5
      
      CL = MIN( CMINI * 0.9, CLOW )
      CH = MAX( CMAXI * 1.1, CHIGH )
      
c      write(*,*) 'Cmini, Cmaxi :', CMINI, CMAXI
                       
      CALL PKSORT( NBPC, BPC )

C     Check for the number of roots inside the specified 
C     interval CLOW CHIGH

       VBP( 1 ) = MIN( CMINI * 0.9, CLOW )
       VBP( 2 ) = CWAT
       VBP( 3 ) = MAX( CMAXI * 1.1, CHIGH )
       NVBP = 3
       
C     Loop to get an approximation of the number of modes
C     between intermediate points defined by vector VBP

      DO K = 1, NVBP - 1

         NBZ( K ) = 1

         IF ( ( VBP(K) .LT. CWAT ) .AND. ( NMEDIA .GE. 2 ) 
     &      .AND. ( CMINI .NE. CSB )  .AND. ( CSMIN .GT. 0.0 ))  THEN
           NBZ( K ) = INT( 2 * ( HENS - HEAU ) * FREQ / CMINI *
     &       SQRT( 1 - CMINI * CMINI / VBP( K + 1 ) / VBP( K + 1 ) ) )
         ENDIF

         IF ( VBP( K ) .GE. CWAT ) THEN
           NBZ( K ) = INT( 2 * HEAU * FREQ / CWAT *
     &       SQRT( 1 - CWAT * CWAT / VBP( K + 1 ) / VBP( K + 1 ) ) )
         ENDIF
c         write(*,*) NBZ( K )
      END DO 

C     If there is too many modes inside one contour split
C     intervals to get less modes inside

      NVITINT = 0

      DO I = 1, NVBP - 1
 
         IF ( NBZ( I ) .GT. MXINTMOD )  THEN
            IF ( VBP( I ) .GE. CWAT ) CREF = CWAT
            IF ( VBP( I ) .LT. CWAT ) CREF = CMINI

            NINTER = NINT( REAL( NBZ( I ) ) / MXINTMOD )
            IF ( NINTER .EQ. 1 ) NINTER = 2
c            write(*,*) 'Ninter:', NINTER
            TMN = ASIN( CREF / VBP( I + 1 ) )
c            write(*,*) 'Tmn :', TMN
            IF ( VBP( I ) .LE. CREF ) THEN
               TMX = PI / 2.0D0
            ELSE
               TMX = ASIN( CREF / VBP( I     ) )
            END IF
c     write(*,*) 'Tmx :', TMX   
            DTM = ( TMX - TMN ) / NINTER
            DO J = 1, NINTER - 1    
               NVITINT = NVITINT + 1
               VITINT( NVITINT ) = CREF / SIN(TMX - DTM*J)
c     write(*,*) 'vint:', VITINT(NVITINT)
            END DO
         ENDIF
      END DO
       
      VITINT( NVITINT +  1 ) = CICR( CPB )
      VITINT( NVITINT +  2 ) = CICR( CSB )
      VITINT( NVITINT +  3 ) = CICR( CPT ) 
      VITINT( NVITINT +  4 ) = CICR( CST ) 
      VITINT( NVITINT +  5 ) = CWAT 
      VITINT( NVITINT +  6 ) = CWAT * 0.99

      NT = NVITINT + 6
       
C     Now get all the intervalls where we are going to use a contour
C     between CLOW and CHIGH

      CBP( 1 ) = CLOW
      J = 1
      DO I = 1, NT
         IF ( ( VITINT( I ) .GT. CLOW ) .AND.
     &        ( VITINT( I ) .LT. CHIGH ) ) THEN
            J = J + 1
            CBP( J ) = VITINT( I )
         ENDIF
      END DO
      CBP( J + 1 ) = CHIGH
      
      NT = J + 1
       
      CALL PKSORT( NT, CBP )

      DO K = 1, NT
         KBP( K ) = OMEGA / CBP( K )
      END DO  


C      End for the process of initialization for the first mesh

       
C      Maximum value for the attenuation of a mode
        
       ALPHAM = 40.0 / 8.69 / RMIN / 1000.0
c        write(*,*) ' CIMX, ALPHAM :', CIMX, ALPHAM

C      Initialization of the total number of modes

       NMODE = 0

       DO K = 1, NT - 1

c         write(*,*) '************************************************'
         
         KBP1 = KBP( K )
         KBP2 = KBP( K + 1 )

         CALL ZONE( RAY, CENT, KBP1, KBP2, ALPHAM, MLT, OMEGA, U )

         CALL RESZONE( MLT, CENT, RAY, U, NCER, NCIR, ANG, FLAGRES,
     &        NBZP, SOL, NMODE, MXCER, MAXM )

         IF ( FLAGRES ) THEN
            
            ZINT = INT( NBZP / MXINTMOD )
c            write(*,*) 'Zint :', ZINT
            KBPINC = ( KBP( K + 1 ) - KBP( K ) ) / DBLE( ZINT )

            DO KK = 1, ZINT
c               write(*,*) '****************************************'
               KBP1 = KBP( K ) +  DBLE( KK - 1 ) * KBPINC
               KBP2 = KBP( K ) +  DBLE( KK ) * KBPINC
               CALL ZONE( RAY, CENT, KBP1, KBP2, ALPHAM, MLT, OMEGA, U )
               CALL RESZONE( MLT, CENT, RAY, U, NCER, NCIR, ANG,
     &              FLAGRES, NBZP, SOL, NMODE, MXCER, MAXM )
            END DO
         ENDIF         
      END DO
        
      M = NMODE

      CALL ZPKSORT( M, SOL )
         
      NMP = 0 
      DO I = 1, M
         IF (SOL( I ) .NE. (0.0,0.0)) THEN
            NMP = NMP + 1
            EVMAT( ISET, NMP ) = SOL( I )           
         ENDIF
c         write(*,*)'Sol:', I, OMEGA/REAL(SQRT(SOL(I)))  
      END DO
      M = NMP
c     write(*,*) 'Total number of modes : ', M                  
      
      RETURN
      END 
      
C**********************************************************************
      SUBROUTINE ZONE( RAY, CENT, KBP1, KBP2, ALPHAM, MLT, OMEGA, U )

      INTEGER*4 MLT
      REAL*8 RAY, RAY1, RAY2, C1, C2, OMEGA, TAILLCER
      REAL*8 KBP1, KBP2, ALPHAM, U
      COMPLEX*16 CENT

      TAILLCER = 0.9999999D0
 
      CENT = KBP1 - ( KBP1 - KBP2 ) / 2.0D0
      RAY  = ( KBP1 - KBP2 ) / 2.0D0 * TAILLCER

      RAY1 = CENT + RAY
      RAY2 = CENT - RAY

      C1 = OMEGA / RAY1
      C2 = OMEGA / RAY2

      U = 2.0D0 * RAY
         
      IF ( ALPHAM .GE. U ) THEN
         MLT = INT( ALPHAM / U )
c         write(*,*) 'MLt * : ', MLT
      ELSE
         MLT = 1
c         write(*,*) 'MLt \ : ', MLT
      ENDIF

c      write(*,*) 'Vit ex:', C1, C2
        
      RETURN
      END         
C************************************************************************
      SUBROUTINE RESZONE( MLT, CENT, RAY, U, NCER, NCIR, ANG, FLAGRES,
     &     NBZP, SOL, NMODE, MXCER, MAXM )
       
      PARAMETER( MXCER = 20000 )
      INTEGER*4 MLT, NCER, NCIR, IRECT, NBZP, NMODE, NRECT
      INTEGER*4 REP( 5 ), IPOW, NBZPC, NPM
      REAL*8 U, ANG( MXCER ), DR( 2 ), RAY
      COMPLEX*16 SOL( MAXM ), T( MXCER ), FY( MXCER ), FYC( MXCER )
      COMPLEX*16 DELTA, X, CENT, SOLP( 30 ), CI
      LOGICAL FLAGRES        
       
      FLAGRES = .FALSE.

      CI = (0.0D0, 1.0D0 )
       
      DO IRECT = 1, MLT
         
         IF ( IRECT .GT. 1 ) THEN         
            CENT = CENT - CI * U
         ENDIF

         CALL RECT( T, REP, DR, CENT, U, U, NCER )

C        Total number of points used to discretize the rectangle

         NRECT  = REP( 5 )
c         write(*,*) 'NRECT : ',NRECT
c         write(*,*)'Nmode2 :',NMODE

C        Determination of the values of the function along the contour

         DO I = 1, NRECT
            X = T( I ) * T( I )
            CALL FUNCTZ( X , DELTA, IPOW, CENT )
            FY( I ) = DELTA * 10.0D0**IPOW
         END DO

         CALL FNBZEROS( NBZP, FY, NRECT)

c         write(*,*) 'Nbzp : ', NBZP
         
         IF ( NBZP .EQ. 0 ) RETURN

         IF ( NBZP .GT. 15 ) THEN
            FLAGRES = .TRUE.
            RETURN
         ENDIF

         DO I = 1, NCIR
            X = EXP( CI * ANG( I ) ) * RAY + CENT
c            write(*,*)'X : ',X
            X = X * X
            CALL FUNCTZ( X, DELTA, IPOW, CENT )
            FYC( I ) = DELTA * 10.0D0**IPOW
c            write(*,*)'FYC :',FYC(I)
         END DO 
                       
         CALL FNBZEROSC( NBZPC, FYC, NCIR)
           
c             write(*,*) 'CENT : ', CENT
c         write(*,*) 'Nbzp, Nbzpc : ', NBZP,NBZPC
          
         IF ( NBZPC .EQ. NBZP ) THEN
c            write(*,*) 'Cercle'
            CALL FSOLANA( SOLP, ANG, FYC, RAY, CENT, NBZPC, NCIR )
         ELSE
c            write(*,*) 'Rectangle'
            CALL FSOLREC( SOLP, T, FY, REP, CENT, NBZP, NRECT, DR )
         ENDIF

         DO NPM = 1, NBZP
            NMODE = NMODE + 1
            SOL( NMODE ) = SOLP( NPM ) ** 2
c               PIM = ABS( DIMAG( CDSQRT( SOL( NMODE ) ) ) )
c               ALPHAMAX = MAX( ALPHAMAX, PIM )
c            write(*,*) 'Mode:', CDSQRT( SOL( NMODE ) )       
         END DO
           
      END DO
      RETURN
      END
        
       
