      PROGRAM EVOL
      PARAMETER(N=2,NITT=2000000)
CFIN  PARAMETER(N=2,NITT=1000000,NTEST=100)
CPAP     PARAMETER(N=2,NITT=4000000)

      REAL*8 PSA(3),PSK(3),PSL(3),PSC(3),PSM(3),PSN(3),
     *     PSADT,PSKDT,PSLDT,PSCDT,PSMDT,PSNDT

      COMMON/VAR/PSA  ,PSK  ,PSL  ,PSC  ,PSM  ,PSN  ,
     *           PSADT,PSKDT,PSLDT,PSCDT,PSMDT,PSNDT

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

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


      open(10,file='../xmgr/papi_6dg_n2_4v5.dat',form='formatted')
      open(11,file='../idl/evol_6dg_n2_4v5.dat',form='formatted')
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  PARAMETRE NUMERIQUES:

      DT=0.005
      ALPHA=0.001 

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)

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


C  INITIALISATION

         PSAS=4.5/ALPN1*BETN1
C  EVOLUTION VERS LA BRANCHE DU HAUT
         PSA(1)=(1.+1.5)/ALPN1*BETN1
         PSK(1)=1.2*K01/H01    
         PSL(1)=-0.25*K01/H01         
C  EVOLUTION VERS LA BRANCHE DU BAS
         PSA(1)=(1.+0.5)/ALPN1*BETN1
         PSK(1)=1.2*K01/H01    
         PSL(1)=-0.25*K01/H01         
C  PAS D'HARMONIQUES C, M et N
         PSC(1)=0.001  
         PSM(1)=0.001          
         PSN(1)=0.001    

C FORWARD EULER STEP

      CALL DXDT(1)

        PSA(2)=PSA(1)+DT*PSADT
        PSK(2)=PSK(1)+DT*PSKDT
        PSL(2)=PSL(1)+DT*PSLDT
        PSC(2)=PSC(1)+DT*PSCDT
        PSM(2)=PSM(1)+DT*PSMDT
        PSN(2)=PSN(1)+DT*PSNDT


      DO 1 IT=2,NITT

      CALL DXDT(2)

      PSA(3)=PSA(1)+2.*DT*PSADT
      PSK(3)=PSK(1)+2.*DT*PSKDT
      PSL(3)=PSL(1)+2.*DT*PSLDT
      PSC(3)=PSC(1)+2.*DT*PSCDT
      PSM(3)=PSM(1)+2.*DT*PSMDT
      PSN(3)=PSN(1)+2.*DT*PSNDT

C  FILTRE DE ASSELIN:

      PSA(1)=PSA(2)+ALPHA/2*(PSA(3)-2.*PSA(2)+PSA(1))
      PSK(1)=PSK(2)+ALPHA/2*(PSK(3)-2.*PSK(2)+PSK(1))
      PSL(1)=PSL(2)+ALPHA/2*(PSL(3)-2.*PSL(2)+PSL(1))
      PSC(1)=PSC(2)+ALPHA/2*(PSC(3)-2.*PSC(2)+PSC(1))
      PSM(1)=PSM(2)+ALPHA/2*(PSM(3)-2.*PSM(2)+PSM(1))
      PSN(1)=PSN(2)+ALPHA/2*(PSN(3)-2.*PSN(2)+PSN(1))

C  STOCKAGE DE T+1 DANS T

      PSA(2)=PSA(3)
      PSK(2)=PSK(3)
      PSL(2)=PSL(3)
      PSC(2)=PSC(3)
      PSM(2)=PSM(3)
      PSN(2)=PSN(3)

c     IF(MOD(IT,100).EQ.0.AND.IT*DT.GE.35000)PRINT *,PSA(1),PSC(1)
c     IF(MOD(IT,10).EQ.0.AND.IT.GE.NITT/4)
C     IF(MOD(IT,500).EQ.0.AND.IT*DT.GE.9500.)
c    *             PRINT *,PSA(1),PSC(1),PSL(1)
c    *             PRINT 100,FLOAT(IT)*DT-9500.,PSA(1),PSK(1),PSL(1)
C     IF(MOD(IT,500).EQ.0.AND.IT*DT.GE.2000.)
      IF(MOD(IT,500).EQ.0.)                  
     *             WRITE(10,101)PSA(1)/BETN1*ALPN1
     *             ,PSL(1)*H01/K01,PSK(1)*H01/K01
     *             ,PSC(1)/BETN1*ALPN1
     *             ,PSM(1)*H01/K01,PSN(1)*H01/K01
      IF(MOD(IT,500).EQ.0.)                  
     *             WRITE(11,100) FLOAT(IT)*DT
     *             ,PSA(1)             
     *             ,PSL(1)                        
     *             ,PSK(1)                        
     *             ,PSC(1)                        
     *             ,PSM(1)
     *             ,PSN(1)            
CBIF  IF(MOD(IT,100).EQ.0.AND.IT*DT.GE.3500.)
CBIF *             PRINT *,PSAS,PSA(1)

100   FORMAT(1x,f8.2,6(1x,E12.4))
101   FORMAT(7(1x,E12.4))

1     CONTINUE
2     CONTINUE


      STOP
      END

      SUBROUTINE DXDT(NT)

      REAL*8 PSA(3),PSK(3),PSL(3),PSC(3),PSM(3),PSN(3),
     *           PSADT,PSKDT,PSLDT,PSCDT,PSMDT,PSNDT

      COMMON/VAR/PSA  ,PSK  ,PSL  ,PSC  ,PSM  ,PSN  ,
     *           PSADT,PSKDT,PSLDT,PSCDT,PSMDT,PSNDT

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

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

      PSADT=K01*(PSAS-PSA(NT))+H01*PSL(NT)

      PSKDT=     -KN1*PSK(NT)             +BETN1*PSL(NT)
     *           -ALPN1*PSA(NT)*PSL(NT)
     *           -DELN1*PSC(NT)*PSN(NT)

      PSLDT=     -KN1*PSL(NT) -HN1*PSA(NT)-BETN1*PSK(NT)
     *           +ALPN1*PSA(NT)*PSK(NT)
     *           +DELN1*PSC(NT)*PSM(NT)

      PSCDT=     -K02*PSC(NT)+H02*PSN(NT)
     *           +EPSN*(PSK(NT)*PSN(NT)-PSL(NT)*PSM(NT))

      PSMDT=     -KN2*PSM(NT)            +BETN2*PSN(NT)
     *           -ALPN2*PSA(NT)*PSN(NT)
     *           -DELTN2*PSC(NT)*PSL(NT)

      PSNDT=     -KN2*PSN(NT)-HN2*PSC(NT)-BETN2*PSM(NT)
     *           +ALPN2*PSA(NT)*PSM(NT)
     *           +DELTN2*PSC(NT)*PSK(NT)

      RETURN
      END
