prenl3.f

Go to the documentation of this file.

!                       *****************
                        SUBROUTINE prenl3
!                       *****************
!
!
!***********************************************************************
! TOMAWAC   V6P1                                   24/06/2011
!***********************************************************************
!
!brief   PREPARES THE COMPUTATION FOR THE NON-LINEAR INTERACTION
!+                SOURCE TERM BETWEEN QUADRUPLETS USING THE GQM METHOD
!+                ("GAUSSIAN QUADRATURE METHOD") PROPOSED BY LAVRENOV
!+                (2001)
!
!+            PROCEDURE SPECIFIC TO THE CASE WHERE THE FREQUENCIES
!+                FOLLOW A GEOMETRICAL PROGRESSION AND THE DIRECTIONS
!+                ARE EVENLY DISTRIBUTED OVER [0;2.PI].
!
!+note    THIS SUBROUTINE IS TO BE USED IN CONJONCT
!+          SUBROUTINE QNLIN3, WHICH IT OPTIMISES.
!
!reference  LAVRENOV, I.V. (2001):
!+           "EFFECT OF WIND WAVE PARAMETER FLUCTUATION ON THE NONLINEAR
!+           SPECTRUM EVOLUTION". J. PHYS. OCEANOGR. 31, 861-873.
!
!history  E. GAGNAIRE-RENOU
!+        04/2011
!+        V6P1
!+   CREATED
!
!history  G.MATTAROLO (EDF - LNHE)
!+        24/06/2011
!+        V6P1
!+   Translation of French names of the variables in argument
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| DIMBUF         |-->| VARIABLE FOR SPECTRUM INTERPOLATION
!| ELIM           |<--| FRACTION OF ELIMINATED CONFIGURATIONS
!| F_COEF         |<--| WORK TABLE FOR SPECTRUM INTERPOLATION
!| F_POIN         |<--| WORK TABLE FOR SPECTRUM INTERPOLATION
!| F_PROJ         |<--| WORK TABLE FOR SPECTRUM INTERPOLATION
!| FREQ           |-->| DISCRETIZED FREQUENCIES
!| IDCONF         |<--| WORK TABLE
!| IQ_OM1         |<--| NUMBER OF INTEGRATION POINT ON OMEGA1
!| K_IF1          |<--| WORK TABLE
!| K_IF2          |<--| WORK TABLE
!| K_IF3          |<--| WORK TABLE
!| K_1M           |<--| WORK TABLE
!| K_1M2M         |<--| WORK TABLE
!| K_1M2P         |<--| WORK TABLE
!| K_1M3M         |<--| WORK TABLE
!| K_1M3P         |<--| WORK TABLE
!| K_1P           |<--| WORK TABLE
!| K_1P2M         |<--| WORK TABLE
!| K_1P2P         |<--| WORK TABLE
!| K_1P3M         |<--| WORK TABLE
!| K_1P3P         |<--| WORK TABLE
!| LBUF           |-->| VARIABLE FOR SPECTRUM INTERPOLATION
!| NCONF          |<--| NUMBER OF RETAINED CONFIGURATIONS
!| NCONFM         |-->| MAXIMUM NUMBER OF CONFIGURATIONS
!| NF             |-->| NUMBER OF FREQUENCIES
!| NF1            |-->| NUMBER OF INTEGRATION POINT ON OMEGA1
!| NQ_OM2         |-->| NUMBER OF INTEGRATION POINT ON OMEGA2
!| NQ_TE1         |-->| SETTING FOR INTEGRATION ON THETA1
!| NT             |-->| NUMBER OF DIRECTIONS
!| NT1            |-->| NUMBER OF INTEGRATION POINT ON THETA1
!| RAISF          |-->| FREQUENTIAL RATIO
!| SEUIL1         |-->| THRESHOLD1 FOR CONFIGURATIONS ELIMINATION
!| SEUIL2         |-->| THRESHOLD2 FOR CONFIGURATIONS ELIMINATION
!| T_POIN         |<--| WORK TABLE FOR SPECTRUM INTERPOLATION
!| TAILF          |-->| SPECTRUM QUEUE FACTOR
!| TB_FAC         |<--| WORK TABLE
!| TB_SCA         |<--| SCALE COEFFICIENT
!| TB_TMP         |<--| WORK TABLE
!| TB_TPM         |<--| WORK TABLE
!| TB_V14         |<--| WORK TABLE for F1
!| TB_V24         |<--| WORK TABLE
!| TB_V34         |<--| WORK TABLE
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!Missing the description of the variables in argument
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_tomawac
!
      USE declarations_special
      USE interface_tomawac, ex_prenl3 => prenl3
      IMPLICIT NONE
!
!
!
!.....LOCAL VARIABLES
!     """""""""""""""""
      INTEGER           JF    , JT    , JF1   , JT1   , NF1P1 , IAUX
      INTEGER           IQ_TE1, IQ_OM2
      DOUBLE PRECISION  EPSI_A, AUX   , CCC   , DENO  , AAA   , DP2SG
      DOUBLE PRECISION  V1    , V1_4  , DV1   , DTETAR
      DOUBLE PRECISION  V2    , V2_4  , V3    , V3_4
      DOUBLE PRECISION  W2    , W2_M  , W2_1  , W_MIL , W_RAD 
      DOUBLE PRECISION  RK0   , XK0   , YK0   , RK1   , XK1   , YK1
      DOUBLE PRECISION  RK2   , XK2P  , YK2P  , XK2M  , YK2M
      DOUBLE PRECISION  RK3   , XK3P  , YK3P  , XK3M  , YK3M
      DOUBLE PRECISION  D01P  , C_D01P, S_D01P, D0AP  , C_D0AP, S_D0AP
      DOUBLE PRECISION  GA2P  , C_GA2P, S_GA2P, GA3P  , C_GA3P, S_GA3P
      DOUBLE PRECISION, ALLOCATABLE :: F1SF(:)
!
!.....Variables related to the Gaussian quadratures
      DOUBLE PRECISION  W_CHE_TE1, W_CHE_OM2, C_LEG_OM2
      DOUBLE PRECISION, ALLOCATABLE :: X_CHE_TE1(:),X_CHE_OM2(:),
     &                                 x_leg_om2(:),w_leg_om2(:)
!
!.....Variables related to the configuration selection
      DOUBLE PRECISION  TEST1 , TEST2
      DOUBLE PRECISION, ALLOCATABLE :: MAXCLA(:)
      INTEGER NT 
!
!.....External functions (USE INTERFACE)
!      DOUBLE PRECISION  COUPLE
!      EXTERNAL          COUPLE
!
      nt=ndire
      dtetar=deupi/dble(nt)
!
!=======================================================================
!     INITIALISATION OF AUXILIAIRY TABLES FOR SPECTRUM INTERPOLATION
!=======================================================================
      DO jf=1,dimbuf
        f_poin(jf)=0
        t_poin(jf)=0
        f_coef(jf)=0.d0
        f_proj(jf)=0.d0
        tb_sca(jf)=0.0d0
      ENDDO
      DO jf=1,lbuf
        f_poin(jf)=1
        f_coef(jf)=0.0d0
        f_proj(jf)=0.0d0
      ENDDO
      DO jf=1,nf
        iaux=lbuf+jf
        f_poin(iaux)=jf
        f_coef(iaux)=1.0d0
        f_proj(iaux)=1.0d0
      ENDDO
      aux=1.d0/raisf**tailf
      DO jf=1,lbuf
        iaux=lbuf+nf+jf
        f_poin(iaux)=nf
        f_coef(iaux)=aux**jf
        f_proj(iaux)=0.0d0
      ENDDO
!
      DO jt=lbuf,1,-1
        t_poin(jt)=nt-mod(lbuf-jt,nt)
      ENDDO
      DO jt=1,nt
        t_poin(lbuf+jt)=jt
      ENDDO
      DO jt=1,lbuf
        t_poin(lbuf+nt+jt)=mod(jt-1,nt)+1
      ENDDO
!======================================================================
!
!=======================================================================
!     COMPUTES SCALE COEFFICIENTS FOR THE COUPLING COEFFICIENT
!=======================================================================
      dp2sg=deupi*deupi/gravit
      DO jf=1,lbuf
        aux=freq(1)/raisf**(lbuf-jf+1)
        tb_sca(jf)=(dp2sg*aux**2)**6/(deupi**3*aux)
      ENDDO
      DO jf=1,nf
        tb_sca(lbuf+jf)=(dp2sg*freq(jf)**2)**6/(deupi**3*freq(jf))
      ENDDO
      DO jf=1,lbuf
        iaux=lbuf+nf+jf
        aux=freq(nf)*raisf**jf
        tb_sca(iaux)=(dp2sg*aux**2)**6/(deupi**3*aux)
      ENDDO
!=======================================================================
!
!=======================================================================
!     COMPUTES VALUES FOR GAUSSIAN QUADRATURES
!=======================================================================
      ALLOCATE(x_che_te1(1:nq_te1),x_che_om2(1:nq_om2))
      ALLOCATE(x_leg_om2(1:nq_om2),w_leg_om2(1:nq_om2))
!
!.....Abscissa and weight (constant) for Gauss-Chebyshev
      DO iq_te1=1,nq_te1
        x_che_te1(iq_te1)=cos(pi*(dble(iq_te1)-0.5d0)/dble(nq_te1))
      ENDDO
      w_che_te1=pi/dble(nq_te1)
      DO iq_om2=1,nq_om2
        x_che_om2(iq_om2)=cos(pi*(dble(iq_om2)-0.5d0)/dble(nq_om2))
      ENDDO
      w_che_om2=pi/dble(nq_om2)
!
!.....Abscissa et weight for Gauss-Legendre
      CALL gauleg
     &( w_leg_om2 , x_leg_om2 , nq_om2 )
      DO iq_om2=1,nq_om2
        x_leg_om2(iq_om2)=0.25d0*(1.d0+x_leg_om2(iq_om2))**2
      ENDDO
!=======================================================================
!
!
!=======================================================================
!     COMPUTES VALUES OF RATIO F1/F AS FUNCTION OF THE IQ_OM1 INDICATOR
!=======================================================================
      nf1p1=nf1+1
      ALLOCATE(f1sf(1:nf1p1))
      CALL f1f1f1 ( f1sf  , nf1   , iq_om1)
!=======================================================================
!
!     ==================================================
!     STARTS LOOP 1 OVER THE RATIOS F1/F0
!     ==================================================
      DO jf1=1,nf1
!       ---------Computes and stores v1=f1/f0 and v1**4
        v1=(f1sf(jf1+1)+f1sf(jf1))/2.d0
          k_if1(jf1)=nint(dble(lbuf)+log(v1)/log(raisf))
        v1_4=v1**4
        tb_v14(jf1)=v1_4
!       ---------Computes and stores dv1=df1/f0
        dv1=f1sf(jf1+1)-f1sf(jf1)
!       ---------Computes the A parameter
        aaa=((1.d0+v1)**4-4.d0*(1.d0+v1_4))/(8.d0*v1**2)
!
!
!
!       =================================================
!       STARTS LOOP 2 OVER THE DELTA_1+ VALUES
!       =================================================
        DO jt1=1,nt1
!
!......Computes the Delta1+ values (=Theta_1-Theta_0) between 0 and Pi.
          IF (jt1.LE.nq_te1) THEN
!           ---------First interval : X from -1 to A
            iq_te1=jt1
            c_d01p=(-1.d0+aaa)/2.d0+(1.d0+aaa)/2.d0*x_che_te1(iq_te1)
            ccc=dv1*sqrt((aaa-c_d01p)/(1.d0-c_d01p))*w_che_te1
          ELSE
!           ---------Second interval : X from A to 1
            iq_te1=jt1-nq_te1
            c_d01p=( 1.d0+aaa)/2.d0+(1.d0-aaa)/2.d0*x_che_te1(iq_te1)
            ccc=dv1*sqrt((c_d01p-aaa)/(1.d0+c_d01p))*w_che_te1
          ENDIF
          s_d01p=sqrt(1.d0-c_d01p*c_d01p)
          d01p  =acos(c_d01p)
          k_1p(jt1,jf1)=lbuf+nint(d01p/dtetar)
          k_1m(jt1,jf1)=lbuf-nint(d01p/dtetar)
!
!         ---------Computes Epsilon_a
          epsi_a=2.d0*sqrt(1.d0+v1_4+2.d0*v1*v1*c_d01p)/(1.d0+v1)**2
!         ---------Computes Delta_A+ and its cosinus
          c_d0ap=(1.d0-v1_4+0.25d0*epsi_a**2*(1.d0+v1)**4)
     &           /(epsi_a*(1.d0+v1)**2)
          s_d0ap=sqrt(1.0d0-c_d0ap*c_d0ap)
          d0ap  = acos(c_d0ap)
!
!.......Integration over OMEGA2 depending on EPS_A
          IF (epsi_a.LT.1.d0) THEN
!        - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!........Case of a single singularity (in OMEGA2-)
!        - - - - - - - - - - - - - - - - - - - - - - - - - - - -
            w2_m=0.5d0*(1.d0-epsi_a/2.d0)
            w2_1=0.5d0
!
            w_rad=w2_1-w2_m
            c_leg_om2=sqrt(w_rad)
!
!        ----------------------------------------------------
!........STARTS LOOP 3 OVER OMEGA_2 (CASE Epsilon_A < 1)
!........Case of a single singularity (in OMEGA2-)
!........Integration over OMEGA2 via GAUSS-LEGENDRE quadrature
!        ----------------------------------------------------
            DO iq_om2=1,nq_om2
!             ---------Computes W2, V2, and V3
              w2=w2_m+w_rad*x_leg_om2(iq_om2)
              v2=w2*(1.d0+v1)
              v2_4=v2**4
              tb_v24(iq_om2,jt1,jf1)=v2_4
              k_if2(iq_om2,jt1,jf1) = nint(dble(lbuf)
     &                               + log(v2)/log(raisf))
              v3=1.d0+v1-v2
              v3_4=v3**4
              tb_v34(iq_om2,jt1,jf1)=v3_4
              k_if3(iq_om2,jt1,jf1) = nint(dble(lbuf)
     &                               + log(v3)/log(raisf))
!             ---------Computes Gamma_2+ et Gamma_3+ angles
              c_ga2p=(epsi_a**2/4.d0+w2**4-(1.d0-w2)**4)/(epsi_a*w2*w2)
              c_ga2p=max(min(c_ga2p,1.d0),-1.d0)
              s_ga2p=sqrt(1.d0-c_ga2p*c_ga2p)
              ga2p  =acos(c_ga2p)
              c_ga3p=(epsi_a**2/4.d0-w2**4+(1.d0-w2)**4)/epsi_a
     &              /(1.d0-w2)**2
              c_ga3p=max(min(c_ga3p,1.d0),-1.d0)
              s_ga3p=sqrt(1.d0-c_ga3p*c_ga3p)
              ga3p  =acos(c_ga3p)
!             Shifting of the direction indexes - Config. +Delta1 (SIG=1)
              k_1p2p(iq_om2,jt1,jf1)=nint(( d0ap+ga2p)/dtetar
     &                              +dble(lbuf))
              k_1p3m(iq_om2,jt1,jf1)=nint(( d0ap-ga3p)/dtetar
     &                              +dble(lbuf))
              k_1p2m(iq_om2,jt1,jf1)=nint(( d0ap-ga2p)/dtetar
     &                              +dble(lbuf))
              k_1p3p(iq_om2,jt1,jf1)=nint(( d0ap+ga3p)/dtetar
     &                              +dble(lbuf))
!             Shifting of the direction indexes - Config. -Delta1 (SIG=-1)
              k_1m2p(iq_om2,jt1,jf1)=nint((-d0ap+ga2p)/dtetar
     &                              +dble(lbuf))
              k_1m3m(iq_om2,jt1,jf1)=nint((-d0ap-ga3p)/dtetar
     &                              +dble(lbuf))
              k_1m2m(iq_om2,jt1,jf1)=nint((-d0ap-ga2p)/dtetar
     &                              +dble(lbuf))
              k_1m3p(iq_om2,jt1,jf1)=nint((-d0ap+ga3p)/dtetar
     &                              +dble(lbuf))
!
!.........Computes the coupling coefficients (only for Delta_1+ )
              rk0=1.d0
              rk1=v1*v1
              rk2=v2*v2
              rk3=(1.d0+v1-v2)**2
              xk0  = rk0
              yk0  = 0.0d0
              xk1  = rk1*c_d01p
              yk1  = rk1*s_d01p
              xk2p = rk2*(c_d0ap*c_ga2p-s_d0ap*s_ga2p)
              yk2p = rk2*(s_d0ap*c_ga2p+c_d0ap*s_ga2p)
              xk2m = rk2*(c_d0ap*c_ga2p+s_d0ap*s_ga2p)
              yk2m = rk2*(s_d0ap*c_ga2p-c_d0ap*s_ga2p)
              xk3p = rk3*(c_d0ap*c_ga3p-s_d0ap*s_ga3p)
              yk3p = rk3*(s_d0ap*c_ga3p+c_d0ap*s_ga3p)
              xk3m = rk3*(c_d0ap*c_ga3p+s_d0ap*s_ga3p)
              yk3m = rk3*(s_d0ap*c_ga3p-c_d0ap*s_ga3p)
              tb_tpm(iq_om2,jt1,jf1)=couple
     &( xk0   , yk0   , xk1   , yk1   , xk2p  , yk2p  , xk3m  , yk3m  )
              tb_tmp(iq_om2,jt1,jf1)=couple
     &( xk0   , yk0   , xk1   , yk1   , xk2m  , yk2m  , xk3p  , yk3p  )
!
!.........Computes the multiplicative coefficient for QNL4
              deno=2.d0*sqrt( (0.5d0*(1.d0+epsi_a/2.d0)-w2)
     &                       *((w2-0.5d0)**2+0.25d0*(1.d0+epsi_a))
     &                       *((w2-0.5d0)**2+0.25d0*(1.d0-epsi_a)) )
              tb_fac(iq_om2,jt1,jf1)=1.d0/(deno*v1*w2*(1.d0-w2))
     &                  /(1.d0+v1)**5 * w_leg_om2(iq_om2)*c_leg_om2* ccc
            ENDDO
!        -----------------------------------------------
!........END OF THE LOOP 3 OVER OMEGA_2 (CASE Epsilon_A < 1)
!        -----------------------------------------------
!
          ELSE
!        - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!........STARTS LOOP 3 OVER OMEGA_2 (CASE Epsilon_A > 1)
!........Case of two singularities (in OMEGA2- and OMEGA2_1)
!........Integration over OMEGA2 via GAUSS-CHEBYSCHEV quadrature
!        - - - - - - - - - - - - - - - - - - - - - - - - - - - -
            w2_m=0.5d0*(1.d0-epsi_a/2.d0)
            w2_1=0.5d0*(1.d0-sqrt(epsi_a-1.d0))
!
            w_mil=(w2_m+w2_1)/2.d0
            w_rad=(w2_1-w2_m)/2.d0
!
            DO iq_om2=1,nq_om2
!             ---------Computes W2, V2, and V3
              w2=w_mil+w_rad*x_che_om2(iq_om2)
              v2=w2*(1.d0+v1)
              v2_4=v2**4
              tb_v24(iq_om2,jt1,jf1)=v2_4
                k_if2(iq_om2,jt1,jf1)=nint(dble(lbuf)
     &                                +log(v2)/log(raisf))
              v3=1.d0+v1-v2
              v3_4=v3**4
              tb_v34(iq_om2,jt1,jf1)=v3_4
              k_if3(iq_om2,jt1,jf1)=nint(dble(lbuf)
     &                              +log(v3)/log(raisf))
!             ---------Computes Gamma_2+ et Gamma_3+ angles
              c_ga2p=(epsi_a**2/4.d0+w2**4-(1.d0-w2)**4)/(epsi_a*w2*w2)
              c_ga2p=max(min(c_ga2p,1.d0),-1.d0)
              s_ga2p=sqrt(1.d0-c_ga2p*c_ga2p)
              ga2p  =acos(c_ga2p)
              c_ga3p=(epsi_a**2/4.d0-w2**4+(1.d0-w2)**4)/epsi_a
     &              /(1.d0-w2)**2
              c_ga3p=max(min(c_ga3p,1.d0),-1.d0)
              s_ga3p=sqrt(1.d0-c_ga3p*c_ga3p)
              ga3p  =acos(c_ga3p)
!             Shifts the direction indexes - Config. +Delta1 (SIG=1)
              k_1p2p(iq_om2,jt1,jf1)=nint(( d0ap+ga2p)/dtetar
     &                              +dble(lbuf))
              k_1p3m(iq_om2,jt1,jf1)=nint(( d0ap-ga3p)/dtetar
     &                              +dble(lbuf))
              k_1p2m(iq_om2,jt1,jf1)=nint(( d0ap-ga2p)/dtetar
     &                              +dble(lbuf))
              k_1p3p(iq_om2,jt1,jf1)=nint(( d0ap+ga3p)/dtetar
     &                              +dble(lbuf))
!             Shifts the direction indexes - Config. -Delta1 (SIG=-1)
              k_1m2p(iq_om2,jt1,jf1)=nint((-d0ap+ga2p)/dtetar
     &                              +dble(lbuf))
              k_1m3m(iq_om2,jt1,jf1)=nint((-d0ap-ga3p)/dtetar
     &                              +dble(lbuf))
              k_1m2m(iq_om2,jt1,jf1)=nint((-d0ap-ga2p)/dtetar
     &                              +dble(lbuf))
              k_1m3p(iq_om2,jt1,jf1)=nint((-d0ap+ga3p)/dtetar
     &                              +dble(lbuf))
!
!.........Computes the coupling coefficients (only for Delta_1+ )
              rk0=1.d0
              rk1=v1*v1
              rk2=v2*v2
              rk3=(1.d0+v1-v2)**2
              xk0  = rk0
              yk0  = 0.0d0
              xk1  = rk1*c_d01p
              yk1  = rk1*s_d01p
              xk2p = rk2*(c_d0ap*c_ga2p-s_d0ap*s_ga2p)
              yk2p = rk2*(s_d0ap*c_ga2p+c_d0ap*s_ga2p)
              xk2m = rk2*(c_d0ap*c_ga2p+s_d0ap*s_ga2p)
              yk2m = rk2*(s_d0ap*c_ga2p-c_d0ap*s_ga2p)
              xk3p = rk3*(c_d0ap*c_ga3p-s_d0ap*s_ga3p)
              yk3p = rk3*(s_d0ap*c_ga3p+c_d0ap*s_ga3p)
              xk3m = rk3*(c_d0ap*c_ga3p+s_d0ap*s_ga3p)
              yk3m = rk3*(s_d0ap*c_ga3p-c_d0ap*s_ga3p)
              tb_tpm(iq_om2,jt1,jf1)=couple
     &( xk0   , yk0   , xk1   , yk1   , xk2p  , yk2p  , xk3m  , yk3m  )
              tb_tmp(iq_om2,jt1,jf1)=couple
     &( xk0   , yk0   , xk1   , yk1   , xk2m  , yk2m  , xk3p  , yk3p  )
!
!.........Computes the multiplicative coefficient for QNL4
              deno=2.d0*sqrt( (0.5d0*(1.d0+epsi_a/2.d0)-w2)
     &                       *((w2-0.5d0)**2+0.25d0*(1.d0+epsi_a))
     &                       *(0.5d0*(1.d0+sqrt(epsi_a-1.d0))-w2) )
              tb_fac(iq_om2,jt1,jf1)=1.d0/(deno*v1*w2*(1.d0-w2))
     &              /(1.d0+v1)**5 * w_che_om2* ccc
!
            ENDDO
!        -----------------------------------------------
!........END OF LOOP 3 OVER OMEGA_2 (CASE Epsilon_A > 1)
!        -----------------------------------------------
!
          ENDIF
!
!
        ENDDO
!       =================================================
!       END OF LOOP 2 OVER THE DELTA_1+ VALUES
!       =================================================
!
      ENDDO
!     ==================================================
!     END OF LOOP 1 OVER THE F1/F0 RATIOS
!     ==================================================
      DEALLOCATE(f1sf)
      DEALLOCATE(x_che_te1)
      DEALLOCATE(x_che_om2)
      DEALLOCATE(x_leg_om2)
      DEALLOCATE(w_leg_om2)
!
!
!
!
!     ===========================================================
!     POST-PROCESSING TO ELIMINATE PART OF THE CONFIGURATIONS
!     ===========================================================
!
!.....It looks, for every value of the ratio V1, for the maximum value
!.....of FACTOR*COUPLING : it is stored in the local table NAXCLA(.)
!     """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
      ALLOCATE(maxcla(1:nf1))
      DO jf1=1,nf1
        aux=0.0d0
        DO jt1=1,nt1
          DO iq_om2=1,nq_om2
            aaa=tb_fac(iq_om2,jt1,jf1)*tb_tpm(iq_om2,jt1,jf1)
            IF (aaa.GT.aux) aux=aaa
            ccc=tb_fac(iq_om2,jt1,jf1)*tb_tmp(iq_om2,jt1,jf1)
            IF (ccc.GT.aux) aux=ccc
          ENDDO
        ENDDO
        maxcla(jf1)=aux
      ENDDO
!
!.....It looks for the max V1 value
!     """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""
      aux=0.0d0
      DO jf1=1,nf1
        IF (maxcla(jf1).GT.aux) aux=maxcla(jf1)
      ENDDO
      test1=seuil1*aux
!
!.....Set to zero the coupling coefficients not used
!     """""""""""""""""""""""""""""""""""""""""""""""""""""
      nconf=0
      DO jf1=1,nf1
        test2 =seuil2*maxcla(jf1)
        DO jt1=1,nt1
          DO iq_om2=1,nq_om2
            aaa=tb_fac(iq_om2,jt1,jf1)*tb_tpm(iq_om2,jt1,jf1)
            ccc=tb_fac(iq_om2,jt1,jf1)*tb_tmp(iq_om2,jt1,jf1)
            IF ((aaa.GT.test1.OR.aaa.GT.test2).OR.
     &          (ccc.GT.test1.OR.ccc.GT.test2)) THEN
              nconf=nconf+1
              idconf(nconf,1)=jf1
              idconf(nconf,2)=jt1
              idconf(nconf,3)=iq_om2
            ENDIF
          ENDDO
        ENDDO
      ENDDO
      DEALLOCATE(maxcla)
!
!.....It counts the fraction of the eliminated configurations
!     """"""""""""""""""""""""""""""""""""""""""""""""""""""
      elim=(1.d0-dble(nconf)/dble(nconfm))*100.d0
      WRITE(lu,*) 'PART DE CONFIGURATIONS ANNULEES : ',elim,' %'
!
      RETURN
      END