      PROGRAM EVOL
      IMPLICIT NONE
     
      INTEGER N,NTEST,NPS
      PARAMETER(N=2,NTEST=20000)


      REAL*8       GAMN1,GAMN2,GAMN3,ALPN1,ALPN2,
     *             DELN1,DELN2,EPSN,K01,K02,KN1,KN2,
     *             BETN1,BETN2,H01,H02,HN1,HN2,PSAS

      REAL*8       PI,REA,RG,RH0,PH0,F0

      REAL*8       RK,RL,RLAM2,RH       

      REAL*8      A,B,C,D,P,Q,DET,R

      REAL*8      PHI,PSA,PSL,PSK


C  PARAMETRE TERRESTRES

      PI=ACOS(-1.)
      REA=6400000.
      RG=3.6
      RH0=1.E+4
      PH0=PI/4.
      F0=1.0284E-4

C  PARAMETRE DE CONFIGURATION:

      RK=0.01
      RL=REA/4
      RLAM2=F0*F0/RG/RH0*RL**2
cCD:
      RH=0.2*RH0
cYO:
c     RH=0.1*RH0
cLO:
c     RH=0.5*RH0

C CONSTANTES D'INTEGRATION (EQ.17, Charney De Vore 1979):

      GAMN1=5.* 8.*sqrt(2.)/15./pi*float(n)
      GAMN2=4.* 8.*sqrt(2.)/15./pi*float(n)
      GAMN3=8.* 8.*sqrt(2.)/15./pi*float(n)

      ALPN1=GAMN1*N*N/(N*N+1+RLAM2)
      ALPN2=GAMN2*(N*N+3)/(N*N+4+RLAM2)

      DELN1=GAMN3*N*N/(N*N+1+RLAM2)
      DELN2=GAMN3*(N*N-3)/(N*N+4+RLAM2)

      EPSN=3.*GAMN3/(4.+RLAM2)

      K01=RK/(1+RLAM2)
      K02=RK/(4+RLAM2)
      KN1=RK*(N*N+1)/(N*N+1.+RLAM2)
      KN2=RK*(N*N+4)/(N*N+4.+RLAM2)


      BETN1=RL/REA*COS(PH0)/SIN(PH0)*N/(N*N+1.+RLAM2)
      BETN2=RL/REA*COS(PH0)/SIN(PH0)*N/(N*N+4.+RLAM2)

      H01=GAMN1/(1.+RLAM2)*RH/(2.*RH0)
      H02=GAMN3/(4.+RLAM2)*RH/(2.*RH0)
      HN1=GAMN1/(N*N+1.+RLAM2)*RH/(2.*RH0)
      HN2=GAMN3/(N*N+4.+RLAM2)*RH/(2.*RH0)

      PRINT *,'GAMN1 GAMN2 GAMN3 ALPN1 ALPN2:'
      PRINT *,GAMN1,GAMN2,GAMN3,ALPN1,ALPN2
      PRINT *,'DELN1 DELN2 EPSN:'
      PRINT *,DELN1,DELN2,EPSN
      PRINT *,'K01 K02 KN1 KN2:'
      PRINT *,K01,K02,KN1,KN2
      PRINT *,'BETN1 BETN2:'
      PRINT *,BETN1,BETN2
      PRINT *,'H01 H02 HN1 HN2 PSAS:'
      PRINT *,H01,H02,HN1,HN2,PSAS
      

C  FICHIERS DE SORTIES POUR TRACER AVEC XMGRACE:

      OPEN(11,FILE='../xmgr/bifu_3d_all.dat',form='formatted')

      DO 2 NPS=1,NTEST
      PSAS=FLOAT(NPS)/FLOAT(NTEST)*10.

C  CONSTANTES DE L'EQUATION DU 3eme DEGRE:

      A=+ALPN1**2*K01/KN1/H01
      B=-ALPN1**2*K01/KN1/H01*PSAS-2.*BETN1*K01*ALPN1/KN1/H01
      C=+KN1*K01/H01+HN1+BETN1**2*K01/KN1/H01
     *  +2.*BETN1*K01*ALPN1/KN1/H01*PSAS
      D=-KN1*K01/H01*PSAS-BETN1**2*K01/KN1/H01*PSAS

      P=(3.*A*C-B**2)/9./A**2
      Q=B*B*B/27./A/A/A-B*C/6./A/A+D/2./A
      DET= P*P*P+Q*Q

      R=SQRT(ABS(P))
      IF(Q.LT.0)R=-R

      IF (P.LT.0) THEN
         IF(DET.LE.0)THEN
           PHI=ACOS(Q/R**3)
           PSA=-2.*R*COS(PHI/3.)-B/3./A
           PSL=K01/H01*(PSA-PSAS)
           PSK=1./KN1*(BETN1*PSL-ALPN1*PSA*PSL)
           WRITE(*,*)   PSAS,PSA,PSL,PSK
           WRITE(11,*)   PSAS/BETN1*ALPN1,PSA/BETN1*ALPN1
     C                   ,PSL*H01/K01,PSK*H01/K01
           PSA=2.*R*COS(PI/3.-PHI/3.)-B/3./A
           PSL=K01/H01*(PSA-PSAS)
           PSK=1./KN1*(BETN1*PSL-ALPN1*PSA*PSL)
           WRITE(*,*)   PSAS,PSA,PSL,PSK
           WRITE(11,*)   PSAS/BETN1*ALPN1,PSA/BETN1*ALPN1
     C                   ,PSL*H01/K01,PSK*H01/K01
           PSA=2.*R*COS(PI/3.+PHI/3.)-B/3./A
           PSL=K01/H01*(PSA-PSAS)
           PSK=1./KN1*(BETN1*PSL-ALPN1*PSA*PSL)
           WRITE(*,*)   PSAS,PSA,PSL,PSK
           WRITE(11,*)   PSAS/BETN1*ALPN1,PSA/BETN1*ALPN1
     C                   ,PSL*H01/K01,PSK*H01/K01
         ELSE
           PHI=DLOG(Q/R**3+SQRT(Q**2/R**6-1))
           PSA=-2.*R*COSH(PHI/3)-B/3./A
           PSL=K01/H01*(PSA-PSAS)
           PSK=1./KN1*(BETN1*PSL-ALPN1*PSA*PSL)
           WRITE(*,*)   PSAS,PSA,PSL,PSK
           WRITE(11,*)   PSAS/BETN1*ALPN1,PSA/BETN1*ALPN1
     C                   ,PSL*H01/K01,PSK*H01/K01
         ENDIF   
      ELSE
         PHI=DLOG(Q/R**3+SQRT(Q**2/R**6+1))
         PSA=-2*R*SINH(PHI/3.)-B/3./A
           PSL=K01/H01*(PSA-PSAS)
           PSK=1./KN1*(BETN1*PSL-ALPN1*PSA*PSL)
           WRITE(*,*)   PSAS,PSA,PSL,PSK
           WRITE(11,*)   PSAS/BETN1*ALPN1,PSA/BETN1*ALPN1
     C                   ,PSL*H01/K01,PSK*H01/K01
      ENDIF
 2    CONTINUE
      
      STOP
      END

