      PROGRAM QBO_1

      IMPLICIT NONE

C  DIMENSIONS VERTICALE ET TEMPORELLES
      INTEGER I,N,IT,NITT,NMOIS
      REAL DT,DTVIS

C  DIMENSION DE LA GRIILE VERTICALE
      PARAMETER (N=101)

C***********************************************
C
C 0.  DECLARATIONS ET CONSTANTES PHYSIQUES
C
C***********************************************

C  VITESSE HORIZONTALE A TROIS INSTANTS:

      REAL UB(N),UN(N),UA(N)

C  TENDANCES DUES A LA DIFFUSION< ET AU FORCAGE PAR
C  LES ONDES.

      REAL ZDUDT_VIS(N),ZDUDT_OND(N)

C  TEMPERATURE, DENSITE, TEMPERATURE POTENTIELLE
C  ET VISCOSITE CINEMATIQUE
 
      REAL T_0(N),RHO_0(N),TH_0(N),RNU(N)

C  PRESSION

      REAL ZPBAS,ZPRE_0(N)       

C  DIMENSIONS, ET GRILLE VERTICALE  DU DOMAINE:

      REAL ZBAS,ZHAUT,ZZ(N)       

C  QUELQUES CONSTANTES

      REAL PI,RG,RGA,RCP,RKA,RVIS,ALPHA

C  FICHIERS DE SORTIE

      OPEN(20,file='../idl/qbo_dat',form='formatted')

C  INITIALISATIONS:
      PI=ACOS(-1.)
      RG=9.81
      RGA=287.
      RCP=1004.
      RKA=287./1004.

      ZBAS =10000.
      ZHAUT=100000.
C PRESSION AU PLUS BAS NIVEAU
      ZPBAS=1.E+3
      RVIS=0.1
C
C  TEMPS D'INTEGRATION PAS DE TEMPS FILTRE DE ASSELIN
C
      DT=3600.
      DTVIS=2.*DT
      NMOIS=30*24*3600/DT
      NITT=144*NMOIS

      PRINT *,'NITT=',NITT

      ALPHA=0.1

C****************************************
C
C 1.  INITIALISATIONS DES PROFILS MOYENS
C
C****************************************

      DO 10 I=1,N
        ZZ(I)=ZBAS+(ZHAUT-ZBAS)*FLOAT(I-1)/FLOAT(N-1)
        T_0(I)=273.
 10   CONTINUE

      ZPRE_0(1)=0.
      DO 11 I=2,N
        ZPRE_0(I)=ZPRE_0(I-1)-RG/RGA*(1./T_0(I-1)+1./T_0(I))/2.
     *                            *(ZZ(I)-ZZ(I-1))
 11   CONTINUE

      DO 12 I=1,N
        ZPRE_0(I)=ZPBAS*EXP(ZPRE_0(I))
        RHO_0(I)=ZPRE_0(I)/RG/T_0(I)
        TH_0(I)=T_0(I)*(ZPBAS/ZPRE_0(I))**RKA 
        UB(I)=-20.*(ZZ(I)-ZZ(1))/(ZZ(N)-ZZ(1))
        UN(I)= UB(I)                           
        UA(I)= UB(I)                           
        RNU(I)=RVIS
C       PRINT *,Z(I)/1000.,ZPRE_0(I),RHO_0(I),T_0(I),
C    *          RGA*T_0(I)/RG*ALOG(1.E+5/ZPRE_0(I))/1000.
 12   CONTINUE

C********************************************
C
C 2.  BOUCLE SUR LE TEMPS
C
C********************************************
      DO 2 IT=1,NITT

        CALL ONDES(IT,DT,N,UN,ZZ,RHO_0,TH_0,RG,ZDUDT_OND)
        CALL DIFFU(DTVIS,N,UB,ZZ,RHO_0,RNU,ZDUDT_VIS)

C  PAS DE LEAPFROG

      DO I=2,N-1
         UA(I)=UN(I)+2.*DT*(ZDUDT_OND(I)+ZDUDT_VIS(I))
      ENDDO
         UA(N)=UA(N-1)

C  FILTRE DE ASSELIN

      DO I=1,N
         UB(I)=UN(I)+ALPHA*(UA(I)-2.*UN(I)+UB(I)) 
         UN(I)=UA(I)                               
      ENDDO

C      IF(MOD(IT,25).EQ.1)
C    C PRINT *,IT*DT/3600./24.,UN(N/2)

C  STOCKAGE ECRITURES

      IF(MOD(IT,NMOIS).EQ.1)THEN
        print *,'SAUVEGARDE EN',it/nmois
        IF(IT.EQ.1)THEN
           WRITE(20,*)N,NITT/NMOIS
           DO I=1,N
              WRITE(20,*)ZZ(I)
           ENDDO
         ENDIF
         DO I=1,N
           WRITE(20,*)UN(I)
         ENDDO
      ENDIF

c     DO I=1,N
c        PRINT *,ZZ(I)/1000.,UN(I),ZDUDT_OND(I),ZDUDT_VIS(I)
c     ENDDO
      
  2   CONTINUE

      STOP
      END
      SUBROUTINE ONDES(IT,PDT,KN,PU,PZ,PRHO,PTH,PG,PDUDT)

      IMPLICIT NONE

C*******************************************************
C
C   SOUS PROGRAMME CALCULANT L'ACCELERATION
C   ET LE FREINAGE DUS AUX ONDES DE GRAVITE.
C
C*******************************************************

C***********************************************
C
C 0.  DECLARATIONS ET CONSTANTES PHYSIQUES
C
C***********************************************   
C
C 0.1 TABLEAUX EXTERIEURS.

C INPUT:
C VITESSE, DENSITE, TEMPERATURE POTENTIELLE, PAS DE TEMPS

      INTEGER KN,IT
      REAL PDT,PU(KN),PZ(KN),PRHO(KN),PTH(KN),PG 

C OUTPUT:
C TENDANCE SUR LA VITESSE

      REAL PDUDT(KN)

C 0.2 TABLEAUX LOCAUX

      INTEGER I,NDAT
      PARAMETER (NDAT=201)

      REAL ZN2(NDAT)         

C 0.3 CARACTERISTIQUE DES ONDES,
C
C   QUANTITE D'HARMONIQUES, NIVEAU CRITIQUE
    
      INTEGER MWA,MW,ICRIT
      PARAMETER (MWA=2)

C     NOMBRE D'ONDE, FREQUENCE ABSOLUE, RMS DU VENT HORIZONTAL:

      REAL ZK(MWA),ZW(MWA),ZMS(MWA)

      DATA ZK/6.25E-07, 6.25E-07/
      DATA ZW/6.25E-06,-6.25E-06/
      DATA ZMS/.5,.5/
C
C     PROFILS VERTICAUX DE LA FREQUENCE INTRINSEQUE
C     DU NOMBRE D'ONDE VERTICAL ET DU FLUX DE PSEUDO MOMENT
C
      REAL OME(NDAT),ZM(NDAT),ZPF(NDAT),AMP(NDAT)
C
C***********************************************
C
C 1.  INITIALISATIONS

      IF(KN.GE.NDAT)STOP 'NDAT TROP PETIT'

      DO 10 I=2,KN-1
         ZN2(I)=PG/PTH(I)*(PTH(I+1)-PTH(I-1))/(PZ(I+1)-PZ(I-1))
         PDUDT(I)=0.
         ZPF(I) =0.
10    CONTINUE

      ZN2(1)=ZN2(2)
      ZN2(KN)=ZN2(KN-1)

C**********************************************
C
C 2.  BOUCLE SUR LES ONDES...................
C
C**********************************************

      DO 2 MW=1,MWA
        
       DO 200 I=1,KN
         OME(I)=ZW(MW)-ZK(MW)*PU(I)
         IF(ABS(OME(I)).GT.1.E-10)THEN
           ZM(I) =-ZK(MW)*SQRT(ZN2(I))/OME(I)
           ZPF(I)=-PRHO(1)*ZK(MW)/ZM(1)*ZMS(MW)**2
         ELSE
           ZM(I) = 0.
           ZPF(I)= 0.
         ENDIF
 200   CONTINUE

C   DETECTION DE NIVEAUX CRITIQUES:

      ICRIT=KN+1
      DO 211 I=2,KN
      IF(OME(I)*OME(I-1).LE.0.AND.ICRIT.EQ.KN+1)ICRIT=I-1
 211  CONTINUE

      DO 212 I=ICRIT,KN
      ZPF(I)=0.
 212  CONTINUE


C   DEFERLEMENT

      DO 202 I=2,ICRIT-1
         AMP(I)=-ZK(MW)*ZPF(I-1)/PRHO(I)/ZM(I) 
         IF(AMP(I).LE.OME(I)**2/ZM(I)**2)THEN
           ZPF(I)=ZPF(I-1)
         ELSE
           ZPF(I)=-PRHO(I)/ZK(MW)*OME(I)**2/ZM(I)
         ENDIF
 202  CONTINUE

C  AJOUT DE LA TENDANCE

      DO 203 I=2,KN-1 
        PDUDT(I)=PDUDT(I)
     C           -1./PRHO(I)*(ZPF(I+1)-ZPF(I-1))/(PZ(I+1)-PZ(I-1))
 203  CONTINUE

C  IMPRESSIONS

      IF(IT.EQ.1.OR.IT.EQ.3456)THEN
      IF(MW.EQ.1)THEN
      PRINT *,'# VITESSE MOYENNE:'                      
      DO I=1,KN
          PRINT *,PZ(I)/1000.,PU(I)
        ENDDO
      ENDIF
      PRINT *,'# FLUX DE PM (mPa):',ZW(MW)/ZK(MW),ICRIT
        DO I=1,KN
          PRINT *,PZ(I)/1000.,ZPF(I)*1000.
        ENDDO
      ENDIF

 2    CONTINUE
     

      RETURN
      END 

      SUBROUTINE DIFFU(PDT,KN,PU,PZ,PRHO,PRNU,PDUDT)

      IMPLICIT NONE

C*******************************************************
C
C  SOUS_PROGRAMMES CALCULANT LA DIFFUSION DE LA 
C         VITESSE SOUS FORME IMPLICITE.
C
C  L'INVERSION SOUS JASCENTE D'UNE MATRICE
C  TRI-DIAGONALE SE FAIT SOUS FORME D'UNE 
C           ELIMINATION GAUSSIANNE.
C
C*******************************************************

C***********************************************
C
C 0.  DECLARATIONS ET CONSTANTES PHYSIQUES
C
C***********************************************   
C
C 0.1 TABLEAUX EXTERIEURS.

C INPUT:
C VITESSE AVANT, APRES, DENSITE, VISCOSITE, PAS DE TEMPS

      INTEGER KN
      REAL PDT,PU(KN),PZ(KN),PRHO(KN),PRNU(KN) 
C OUTPUT:
C TENDANCE SUR LA VITESSE

      REAL PDUDT(KN)

C 0.2 TABLEAUX LOCAUX

      INTEGER I,NDAT
      PARAMETER (NDAT=200)

      REAL ZA(NDAT),ZB(NDAT),ZC(NDAT)
      REAL ALPHA(NDAT),BETA(NDAT),ZU2(NDAT)

C***********************************************
C
C 1.  INITIALISATIONS

      IF(KN.GT.NDAT)STOP 'NDAT TROP PETIT'

      DO 10 I=2,KN-1
         ZA(I)=   - 2.*PDT*(PRNU(I+1)+PRNU(I  ))
     C            /(PZ(I+1)-PZ(I  ))/(PZ(I+1)-PZ(I-1)) 
         ZB(I)=1. + 2.*PDT*(PRNU(I  )+PRNU(I-1))
     C            /(PZ(I  )-PZ(I-1))/(PZ(I+1)-PZ(I-1)) 
     C            +2.*PDT*(PRNU(I  )+PRNU(I-1))
     C            /(PZ(I  )-PZ(I-1))/(PZ(I+1)-PZ(I-1)) 
         ZC(I)=-2.*PDT*(PRNU(I  )+PRNU(I-1))
     C            /(PZ(I  )-PZ(I-1))/(PZ(I+1)-PZ(I-1)) 
10    CONTINUE

C
C MONTEE
C

      ALPHA(1)=0.
      BETA(1)=0.

      DO 11 I=2,KN-1
        ALPHA(I)=-ZA(I)/(ZB(I)+ZC(I)*ALPHA(I-1))
        BETA (I)=-(ZC(I)*BETA (I-1)-PU(I))/(ZB(I)+ZC(I)*ALPHA(I-1))
11    CONTINUE

C
C DESCENTE (RESULTAT SOUS FORME DE TENDANCE)
      ZU2(KN)=0.
      DO 12 I=KN-1,2,-1
        ZU2(I)=ALPHA(I)*ZU2(I+1)+BETA(I)
        PDUDT(I)=(ZU2(I)-PU(I))/PDT
 12    CONTINUE

      RETURN
      END 

