qnlin3.f

Go to the documentation of this file.

!                       *****************
                        SUBROUTINE qnlin3
!                       *****************
!
     &( t1tot,t1der, f2, n1poin2, n1plan, n1f )
!
!***********************************************************************
! TOMAWAC   V6P1                                   24/06/2011
!***********************************************************************
!
!brief    COMPUTES THE CONTRIBUTION OF THE NON-LINEAR INTERACTIONS
!+                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 USES THE OUTPUT FROM 'PRENL3' TO OPTIMISE
!+          THE COMPUTATIONS FOR DIA.
!
!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
!| F              |-->| DIRECTIONAL SPECTRUM
!| F_COEF         |-->| WORK TABLE FOR SPECTRUM INTERPOLATION
!| F_POIN         |-->| WORK TABLE FOR SPECTRUM INTERPOLATION
!| F_PROJ         |-->| WORK TABLE FOR SPECTRUM INTERPOLATION
!| FSEUIL         |-->| WORK TABLE
!| 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
!| IDCONF         |-->| WORK TABLE
!| LBUF           |-->| VARIABLE FOR SPECTRUM INTERPOLATION
!| NB_NOD         |-->| NUMBER OF POINTS IN 2D MESH
!| 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
!| NT             |-->| NUMBER OF DIRECTIONS
!| NT1            |-->| NUMBER OF INTEGRATION POINT ON TETA1
!| RAISF          |-->| FREQUENTIAL RATIO
!| SEUIL          |-->| THRESHOLD0 FOR CONFIGURATIONS ELIMINATION (GQM)
!| T_POIN         |-->| WORK TABLE FOR SPECTRUM INTERPOLATION
!| TB_FAC         |-->| WORK TABLE
!| TB_SCA         |-->| SCALE COEFFICIENT
!| TB_TMP         |-->| WORK TABLE
!| TB_TPM         |-->| WORK TABLE
!| TB_V14         |-->| WORK TABLE
!| TB_V24         |-->| WORK TABLE
!| TB_V34         |-->| WORK TABLE
!| TSDER          |<->| DERIVED PART OF THE SOURCE TERM CONTRIBUTION
!| TSTOT          |<->| TOTAL PART OF THE SOURCE TERM CONTRIBUTION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE interface_tomawac, ex_qnlin3 => qnlin3
      USE declarations_tomawac
!     
      IMPLICIT NONE
!
!.....VARIABLES IN ARGUMENT
!     """"""""""""""""""""
!     A WAY TO CHANGE DIMENSION OF ARRAY
!     TSTOT,TSDER AND F HAVE ONLY ONE DIMENSION
      
      INTEGER, INTENT(IN)    :: N1POIN2,N1PLAN,N1F
      DOUBLE PRECISION, INTENT(INOUT) :: T1TOT( n1poin2,n1plan,n1f)
      DOUBLE PRECISION, INTENT(INOUT) :: T1DER( n1poin2,n1plan,n1f)
      DOUBLE PRECISION, INTENT(INOUT) :: F2(n1poin2,n1plan,n1f)
!
!.....LOCAL VARIABLES
!     """""""""""""""""
      INTEGER           IP    , JF    , JT    , JF1   , JT1   , IQ_OM2,
     &                  jfm0  , jfm1  , jfm2  , jfm3  , ixf1  , ixf2  ,
     &                  ixf3  , jfmin , jfmax , iconf
      INTEGER           KT1P  , KT1M  , JT1P  , JT1M  , KT1P2P, KT1P2M,
     &                  kt1p3p, kt1p3m, kt1m2p, kt1m2m, kt1m3p, kt1m3m,
     &                  jt1p2p, jt1p2m, jt1p3p, jt1p3m, jt1m2p, jt1m2m,
     &                  jt1m3p, jt1m3m
      DOUBLE PRECISION  V1_4  , V2_4  , V3_4  , Q_2P3M, Q_2M3P, FACTOR,
     &                  t_2p3m, t_2m3p, s_2p3m, s_2m3p, scal_t, t2p3m ,
     &                  t2m3p , sp0   , sp1p  , sp1m  , sp1p2p, sp1p2m,
     &                  sp1p3p, sp1p3m, sp1m2p, sp1m2m, sp1m3p, sp1m3m,
     &                  cf0   , cp0   , cf1   , cp1   , cf2   , cp2   ,
     &                  cf3   , cp3   , q2pd0 , q2pd1 , q2pd2p, q2pd3m,
     &                  q2md0 , q2md1 , q2md2m, q2md3p,
     &                  aux00 , aux01 , aux02 , aux03 , aux04 , aux05 ,
     &                  aux06 , aux07 , aux08 , aux09 , aux10
!
!=======================================================================
!     COMPUTES THE GENERALIZED MIN AND MAX FREQUENCIES : INSTEAD OF GOING
!     FROM 1 TO NF IN FREQ(JF) FOR THE MAIN FREQUENCY, IT GOES FROM JFMIN
!     TO JFMAX
!     JFMIN IS GIVEN BY Fmin=FREQ(1) /Gamma_min
!     JFMAX IS GIVEN BY Fmax=FREQ(NF)*Gamma_max
!     TESTS HAVE SHOWN THAT IT CAN BE ASSUMED Gamma_min=1. (JFMIN=1) AND
!     Gamma_max=1.3 (JFMAX>NF) TO OBTAIN IMPROVED RESULTS
!=======================================================================
      jfmin= 1-int(log(1.0d0)/log(raisf))
      jfmax=nf+int(log(1.3d0)/log(raisf))
!
!=======================================================================
!     COMPUTES THE SPECTRUM THRESHOLD VALUES (BELOW WHICH QNL4 IS NOT
!     CALCULATED). THE THRESHOLD IS SET WITHIN 0 AND 1.
!=======================================================================
      DO ip=1,npoin2
        aux00=0.0d0
        DO jf=1,nf
          DO jt=1,ndire
            IF (f2(ip,jt,jf).GT.aux00) aux00=f2(ip,jt,jf)
          ENDDO
        ENDDO
        tra37(ip)=aux00*seuil
      ENDDO
!=======================================================================
!
!
!
!
!     ==================================================
!     STARTS LOOP 1 OVER THE SELECTED CONFIGURATIONS
!     ==================================================
      DO iconf=1,nconf
!       ---------selected configuration characteristics
        jf1   =idconf(iconf,1)
        jt1   =idconf(iconf,2)
        iq_om2=idconf(iconf,3)
!
!       ---------Recovers V1**4=(f1/f0)**4
        v1_4  =tb_v14(jf1)
!       ---------Recovers the shift of the frequency index on f1
        ixf1  =k_if1(jf1)
!       ---------Recovers the direction indexes for Delat1
        kt1p  =k_1p(jt1,jf1)
        kt1m  =k_1m(jt1,jf1)
!       ---------Recovers V2**4=(f2/f0)**4 and V3**4=(f3/f0)**4
        v2_4  =tb_v24(iq_om2,jt1,jf1)
        v3_4  =tb_v34(iq_om2,jt1,jf1)
!       ---------Recovers the frequency indexes shift on f2 and f3
        ixf2  =k_if2(iq_om2,jt1,jf1)
        ixf3  =k_if3(iq_om2,jt1,jf1)
!       ---------Recovers the direction indexes shift
        kt1p2p=k_1p2p(iq_om2,jt1,jf1)
        kt1p2m=k_1p2m(iq_om2,jt1,jf1)
        kt1p3p=k_1p3p(iq_om2,jt1,jf1)
        kt1p3m=k_1p3m(iq_om2,jt1,jf1)
        kt1m2p=k_1m2p(iq_om2,jt1,jf1)
        kt1m2m=k_1m2m(iq_om2,jt1,jf1)
        kt1m3p=k_1m3p(iq_om2,jt1,jf1)
        kt1m3m=k_1m3m(iq_om2,jt1,jf1)
!       ---------Recovers the coupling coefficients
        t2p3m =tb_tpm(iq_om2,jt1,jf1)
        t2m3p =tb_tmp(iq_om2,jt1,jf1)
!       ---------Recovers the multiplicative factor of QNL4
        factor=tb_fac(iq_om2,jt1,jf1)
!
!       = = = = = = = = = = = = = = = = = = = = = = = = =
!       STARTS LOOP 2 OVER THE SPECTRUM FREQUENCIES
!       = = = = = = = = = = = = = = = = = = = = = = = = =
        DO jf=jfmin,jfmax
!
!.........Recovers the coefficient for the coupling factor
!.........Computes the coupling coefficients for the case +Delta1 (SIG=1)
          scal_t=tb_sca(lbuf+jf)*factor
          t_2p3m=t2p3m*scal_t
          t_2m3p=t2m3p*scal_t
!
!.........Frequency indexes and coefficients
          jfm0=f_poin(jf+lbuf)
          cf0 =f_coef(jf+lbuf)
          cp0 =f_proj(jf+lbuf)
          jfm1=f_poin(jf+ixf1)
          cf1 =f_coef(jf+ixf1)
          cp1 =f_proj(jf+ixf1)
          jfm2=f_poin(jf+ixf2)
          cf2 =f_coef(jf+ixf2)
          cp2 =f_proj(jf+ixf2)
          jfm3=f_poin(jf+ixf3)
          cf3 =f_coef(jf+ixf3)
          cp3 =f_proj(jf+ixf3)
!
!         -------------------------------------------------
!         STARTS LOOP 3 OVER THE SPECTRUM DIRECTIONS
!         -------------------------------------------------
          DO jt=1,ndire
!
!...........Direction indexes
!           direct config (+delta1) (sig =1)
            jt1p  =t_poin(jt+kt1p)
            jt1p2p=t_poin(jt+kt1p2p)
            jt1p2m=t_poin(jt+kt1p2m)
            jt1p3p=t_poin(jt+kt1p3p)
            jt1p3m=t_poin(jt+kt1p3m)
!           image config (-delta1)
            jt1m  =t_poin(jt+kt1m)
            jt1m2p=t_poin(jt+kt1m2p)
            jt1m2m=t_poin(jt+kt1m2m)
            jt1m3p=t_poin(jt+kt1m3p)
            jt1m3m=t_poin(jt+kt1m3m)
!
!           - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!           STARTS LOOP 4 OVER THE MESH NODES
!           - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
            DO ip=1,npoin2
!
              sp0=f2(ip,jt,jfm0)*cf0
!
              IF (sp0.GT.tra37(ip)) THEN
!
!               Config. +Delta1 (SIG=1)
!               =======================
!...............Computes the spectrum values in 1, 2, 3
                sp1p  =f2(ip,jt1p  ,jfm1)*cf1
                sp1p2p=f2(ip,jt1p2p,jfm2)*cf2
                sp1p3m=f2(ip,jt1p3m,jfm3)*cf3
                sp1p2m=f2(ip,jt1p2m,jfm2)*cf2
                sp1p3p=f2(ip,jt1p3p,jfm3)*cf3
!
!...............Computes auxiliary products and variables
                aux01=sp0*v1_4+sp1p
                aux02=sp0*sp1p
                aux03=sp1p2p*sp1p3m
                aux04=sp1p2p*v3_4+sp1p3m*v2_4
                aux05=sp1p2m*sp1p3p
                aux06=sp1p2m*v3_4+sp1p3p*v2_4
                aux07=aux02*v3_4
                aux08=aux02*v2_4
!
!...............Computes the components of the transfer term
                s_2p3m=aux03*aux01-aux02*aux04
                s_2m3p=aux05*aux01-aux02*aux06
                q_2p3m=t_2p3m*s_2p3m
                q_2m3p=t_2m3p*s_2m3p
                aux00 =q_2p3m+q_2m3p
!
!...............Computes the components of the derived terms (dQ/dF)
                q2pd0 =t_2p3m*(aux03*v1_4   - sp1p*aux04)*cf0
                q2pd1 =t_2p3m*(aux03        - sp0 *aux04)*cf1
                q2pd2p=t_2p3m*(aux01*sp1p3m - aux07     )*cf2
                q2pd3m=t_2p3m*(aux01*sp1p2p - aux08     )*cf3
                  q2md0 =t_2m3p*(aux05*v1_4   - sp1p*aux06)*cf0
                  q2md1 =t_2m3p*(aux03        - sp0 *aux06)*cf1
                  q2md2m=t_2m3p*(aux01*sp1p3p - aux07     )*cf2
                  q2md3p=t_2m3p*(aux01*sp1p2m - aux08     )*cf3
                aux09=q2pd0+q2md0
                aux10=q2pd1+q2md1
!
!...............Sum of Qnl4 term in the table TSTOT
                t1tot(ip,jt    ,jfm0)=t1tot(ip,jt    ,jfm0)+aux00 *cp0
                t1tot(ip,jt1p  ,jfm1)=t1tot(ip,jt1p  ,jfm1)+aux00 *cp1
                t1tot(ip,jt1p2p,jfm2)=t1tot(ip,jt1p2p,jfm2)-q_2p3m*cp2
                t1tot(ip,jt1p2m,jfm2)=t1tot(ip,jt1p2m,jfm2)-q_2m3p*cp2
                t1tot(ip,jt1p3m,jfm3)=t1tot(ip,jt1p3m,jfm3)-q_2p3m*cp3
                t1tot(ip,jt1p3p,jfm3)=t1tot(ip,jt1p3p,jfm3)-q_2m3p*cp3
!
!...............Sum of the term dQnl4/dF in the table T1DER
                t1der(ip,jt    ,jfm0)=t1der(ip,jt    ,jfm0)+aux09 *cp0
                t1der(ip,jt1p  ,jfm1)=t1der(ip,jt1p  ,jfm1)+aux10 *cp1
                t1der(ip,jt1p2p,jfm2)=t1der(ip,jt1p2p,jfm2)-q2pd2p*cp2
                t1der(ip,jt1p2m,jfm2)=t1der(ip,jt1p2m,jfm2)-q2md2m*cp2
                t1der(ip,jt1p3m,jfm3)=t1der(ip,jt1p3m,jfm3)-q2pd3m*cp3
                t1der(ip,jt1p3p,jfm3)=t1der(ip,jt1p3p,jfm3)-q2md3p*cp3
!
!               Config. -Delta1 (SIG=-1)
!               ========================
!...............Computes the spectrum values in 1, 2, 3
                sp1m  =f2(ip,jt1m  ,jfm1)*cf1
                sp1m2p=f2(ip,jt1m2p,jfm2)*cf2
                sp1m3m=f2(ip,jt1m3m,jfm3)*cf3
                sp1m2m=f2(ip,jt1m2m,jfm2)*cf2
                sp1m3p=f2(ip,jt1m3p,jfm3)*cf3
!
!...............Computes auxiliary products and variables
                aux01=sp0*v1_4+sp1m
                aux02=sp0*sp1m
                aux03=sp1m2p*sp1m3m
                aux04=sp1m2p*v3_4+sp1m3m*v2_4
                aux05=sp1m2m*sp1m3p
                aux06=sp1m2m*v3_4+sp1m3p*v2_4
                aux07=aux02*v3_4
                aux08=aux02*v2_4
!
!...............Computes the transfer term components
                s_2p3m=aux03*aux01-aux02*aux04
                s_2m3p=aux05*aux01-aux02*aux06
                q_2p3m=t_2m3p*s_2p3m
                q_2m3p=t_2p3m*s_2m3p
                aux00 =q_2p3m+q_2m3p
!
!...............Computes the derived terms components (dQ/dF)
                q2pd0 =t_2p3m*(aux03*v1_4   - sp1m*aux04)*cf0
                q2pd1 =t_2p3m*(aux03        - sp0 *aux04)*cf1
                q2pd2p=t_2p3m*(aux01*sp1m3m - aux07     )*cf2
                q2pd3m=t_2p3m*(aux01*sp1m2p - aux08     )*cf3
                q2md0 =t_2m3p*(aux05*v1_4   - sp1m*aux06)*cf0
                q2md1 =t_2m3p*(aux03        - sp0 *aux06)*cf1
                q2md2m=t_2m3p*(aux01*sp1m3p - aux07     )*cf2
                q2md3p=t_2m3p*(aux01*sp1m2m - aux08     )*cf3
                aux09=q2pd0+q2md0
                aux10=q2pd1+q2md1
!
!...............Sum of Qnl4 term in the table T1TOT
                t1tot(ip,jt    ,jfm0)=t1tot(ip,jt    ,jfm0)+aux00 *cp0
                t1tot(ip,jt1m  ,jfm1)=t1tot(ip,jt1m  ,jfm1)+aux00 *cp1
                t1tot(ip,jt1m2p,jfm2)=t1tot(ip,jt1m2p,jfm2)-q_2p3m*cp2
                t1tot(ip,jt1m2m,jfm2)=t1tot(ip,jt1m2m,jfm2)-q_2m3p*cp2
                t1tot(ip,jt1m3m,jfm3)=t1tot(ip,jt1m3m,jfm3)-q_2p3m*cp3
                t1tot(ip,jt1m3p,jfm3)=t1tot(ip,jt1m3p,jfm3)-q_2m3p*cp3
!
!...............Sum of the term dQnl4/dF in the table T1DER
                t1der(ip,jt    ,jfm0)=t1der(ip,jt    ,jfm0)+aux09 *cp0
                t1der(ip,jt1m  ,jfm1)=t1der(ip,jt1m  ,jfm1)+aux10 *cp1
                t1der(ip,jt1m2p,jfm2)=t1der(ip,jt1m2p,jfm2)-q2pd2p*cp2
                t1der(ip,jt1m2m,jfm2)=t1der(ip,jt1m2m,jfm2)-q2md2m*cp2
                t1der(ip,jt1m3m,jfm3)=t1der(ip,jt1m3m,jfm3)-q2pd3m*cp3
                t1der(ip,jt1m3p,jfm3)=t1der(ip,jt1m3p,jfm3)-q2md3p*cp3
!
              ENDIF
!
            ENDDO
!           - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!           END OF LOOP 4 OVER THE MESH NODES
!           - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
          ENDDO
!         -------------------------------------------------
!         END OF LOOP 3 OVER THE SPECTRUM DIRECTIONS
!         -------------------------------------------------
!
        ENDDO
!       = = = = = = = = = = = = = = = = = = = = = = = = =
!       END OF LOOP 2 OVER THE SPECTRUM FREQUENCIES
!       = = = = = = = = = = = = = = = = = = = = = = = = =
!
      ENDDO
!     ==================================================
!     END OF LOOP 1 OVER THE SELECTED CONFIGURATIONS
!     ==================================================
!
      RETURN
      END