berkho.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE berkho
!                   *****************
!
     &(lf)
!
!***********************************************************************
! ARTEMIS   V7P3                                     Aug 2017
!***********************************************************************
!
!brief    SOLVES THE BERKHOFF EQUATION MODIFIED BY
!+                THE INTRODUCTION OF DISSIPATION TERMS.
!code
!+      DIV (C*CG*GRAD(PHI)) + C*CG*( K**2 + I*K*MU ) * PHI = 0
!+                                           ------
!+
!+ PHI IS A COMPLEX FUNCTION (REAL COMPONENT: PHIR AND IMAGINARY
!+ COMPONENT: PHII)
!+
!+ MU IS A DISSIPATION COEFFICIENT (A PRIORI UNKNOWN)
!+
!+ THE BOUNDARY CONDITIONS COUPLE THE EQUATIONS IN PHIR AND PHII
!+ THEY ARE:
!+
!+ D (PHI) /DN - I*K*PHI = D (F) /DN - I*K*F (N: EXTERNAL NORMAL)
!+ FOR A LIQUID BOUNDARY WITH INCIDENT WAVE CONDITION DEFINED
!+ BY THE POTENTIAL F (F=0 FOR A FREE EXIT)
!+
!+ D (PHI) /DN - I* (1-R*EXP (I*ALFA))/(1 + R*EXP (I*ALFA))*K*COS (TETA) *PHI = 0
!+ FOR A SOLID BOUNDARY, WITH WALL REFLEXION COEFFICIENT: R,
!+ ANGLE OF INCIDENCE OF THE WAVES ON THE WALL: TETA, AND DEPHASING
!+ CAUSED BY THE WALL: ALFA.
!+
!+ THUS GENERALLY :
!+ D(PHIR)/DN = APHIRB*PHII + BPHIRB*PHIR + CPHIRB
!+ D(PHII)/DN =-APHIRB*PHIR + BPHIRB*PHII + DPHIRB
!+
!+
!+ AFTER VARIATIONAL FORMULATION :
!+
!+         (  AM1          BM1     )  ( PHIR )   ( CV1 )
!+         (                       )  (      ) = (     )
!+         (                       )  (      )   (     )
!+         (  -BM1         AM1     )  ( PHII )   ( CV2 )
!+
!+           /
!+ AM1 =    / C*CG * GRAD(PSII)*GRAD(PSIJ) DS
!+         /S
!+
!+           /
!+       -  / OMEGA**2 * CG/C * PSII*PSIJ  DS
!+         /S
!+
!+           /
!+       -  /  BPHIRB * PSII*PSIJ  DB
!+         /B
!+
!+           /                         /
!+ BM1 =  - /  APHIR * PSII*PSIJ DB + /  C*CG* K * MU * PSII * PSIJ DS
!+         /B                        /S
!+
!+          /
!+ CV1 =   /   CPHIR * PSII DB
!+        /B
!+
!+          /
!+ CV2 =   /   CPHII * PSII DB
!+        /B
!+
!+
!+ WHERE S IS THE COMPUTATIONAL DOMAIN AND B ITS BOUNDARY
!+       PSII AND PSIJ ARE THE BASIC FUNCTIONS AT NODES I AND J
!+
!+ GIVEN THAT APHII=-APHIR, BM1 IS ALSO IN THE EQUATION IN PHII.
!
!history  J-M HERVOUET (LNH)
!+
!+
!+   LINKED TO BIEF 5.0
!
!history  D. AELBRECHT (LNH)
!+        04/06/1999
!+        V5P1
!+
!
!history
!+        02/04/2007
!+
!+   INVERSION OF THE SECOND EQUATION BEFORE CALL TO SOLVE IF DIRECT
!+   SOLVEUR IS USED
!
!history  C. DENIS (SINETICS)
!+        18/03/2010
!+        V6P0
!+
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!history  C.PEYRARD (EDF)
!+        2011
!+        V6P1
!+        2ND ORDER BOTTOM EFFECTS ADDED TO BERKHOF EQUATION
!
!history  C.PEYRARD (EDF)
!+        2012
!+        V6P2
!+        NEW TYPE OF BOUNDARY CONDITIONS ADDED
!+        DIFFERENT IMPLEMENTATION OF CV1/CV2
!
!history  C.PEYRARD (EDF)
!+        2013
!+        V6P3
!+        INTEGRATION OF WAVE/CURRENT INTERACTION
!+        AND LOOP ON WAVE VECTOR DIRECTION ADDED
!
!history  C.PEYRARD (EDF)
!+        2014
!+        V7P0
!+        AUTOMATIC ANGLE CALCULATION
!
!history  J,RIEHME (ADJOINTWARE)
!+        November 2016
!+        V7P2
!+   Replaced EXTERNAL statements to parallel functions / subroutines
!+   by the INTERFACE_PARALLEL
!
!history  N.DURAND (HRW)
!+        August 2017
!+        V7P3
!+   DEGRAD and RADDEG now defined in DECLARATIONS_ARTEMIS
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| LF             |-->| INDICATOR OF FIRST RESOLUTION FOR DISSIPATION LOOP
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_telemac
      USE declarations_artemis
      USE interface_artemis, ex_berkho => berkho
!
      USE declarations_special
      USE interface_parallel, ONLY : p_max
      IMPLICIT NONE
!
      INTEGER, INTENT(IN) :: LF
!
      INTEGER I,ITERMU,IG

      DOUBLE PRECISION ECRHMU,MODHMU
!
      DOUBLE PRECISION CBID
!
!--> VARIABLES FOR CURRENT AND TETAP LOOP
      DOUBLE PRECISION ERREUR1,ERREURT
      DOUBLE PRECISION XMUL,XK,ZERO
      DOUBLE PRECISION TETA(nptfr)   , ANGDIR(nptfr)
!
      INTEGER ITERKN
      INTEGER TMP
!
!-----------------------------------------------------------------------
!
      INTRINSIC abs,min,max,log,cos,sin,atan
!
!-----------------------------------------------------------------------
!
!--------------------------
!      LANGAUTO=.FALSE.
!----------------------------------------------------------------------
!
! INITIALISES MU AND FW: SET TO 0
!         FOR THE FIRST ITERATION
!
      itermu=0
      IF (lf.EQ.0) THEN
        CALL os('X=C     ', x=mu, c=0.d0)
        CALL os('X=C     ', x=fw, c=0.d0)
      ENDIF
!
!-----------------------------------------------------------------------
      iterkn=0
      IF (courant) THEN
        zero=1e-10
!       INITIALISATION OF WAVE VECTOR COMPONENTS X&Y : T5 , T6
        CALL os('X=0     ',x=kancx)
        CALL os('X=0     ',x=kancy)
      ENDIF
!      WRITE(LU,*) 'AVANT BOUCLE 98'

!
!-----------------------------------------------------------------------
!
!     =========================================
!98 : DISSIPATION : ITERATIVE LOOP : ON THE VARIABLE MU R ON THE WAVE VECTOR
98    CONTINUE
!
!   ITERATIVE LOOP ON TETAP (AUTOMATIC CALCULATION) AND WAVE VECTOR (WAVE-CURRENT)
      IF (    ((courant).AND.(iterkn.GT.0))
     &    .OR.((langauto ).AND.(iterkn.GT.0))  ) THEN

!   ==> CURRENT :
        IF (courant) THEN
!       COMPUTE WAVE VECTOR (FIRST ITERATION U=0, SO NO NEED TO DO THAT)
!       ---------------------------------------------------------------
          IF(debug.GT.0) WRITE(lu,*) ' - COMPUTING WAVE VECTOR (1st)'
          DO i=1,npoin
            xk =k%R(i)
            CALL solvelambda(xk,
     &                   uc%R(i),vc%R(i),kancx%R(i),kancy%R(i),h%R(i))
            k%R(i) =xk

            wr%R(i)=sqrt(grav*k%R(i)*tanh(k%R(i)*h%R(i)))
            c%R(i) =wr%R(i)/k%R(i)
            cg%R(i)=0.5d0*c%R(i)*
     &           (1.d0 + 2.d0*k%R(i)*h%R(i)/sinh(2.d0*k%R(i)*h%R(i)))
          ENDDO
          IF(debug.GT.0) WRITE(lu,*) ' - WAVE VECTOR (1st) COMPUTED'
!
        ENDIF
!       ------------------------------------------------------
!       ==> AUTOMATIC ANGLES
        IF (langauto) THEN
!       COMPUTE TETAP ON THE BOUNDARY (FIRST ITERATION TETAP GIVEN BY USER, SO NO NEED TO DO THAT)
!       ---------------------------------------------------------------
!        ==> RELAXATION ON TETAP IF NECESARY : TETAP=TETAP + C*(TETAP-TETAPM) (default : C=0)
          DO i=1,nptfr
            tetap%R(i)=tetapm%R(i)+reltp*(tetap%R(i)-tetapm%R(i))
            tetap%R(i)=min(tetap%R(i),90d0)
          ENDDO
!         TETAP STORAGE IN TETAPM
!         --------------------------
          CALL os( 'X=Y     ' , x=tetapm , y=tetap )
        ENDIF
!
!       ------------------------------------------------------
!        ACTUALIZATION OF BOUNDARY CONDITIONS
!           1/ CURRENT     : K HAS CHANGED
!           2/ AUTO ANGLES : TETAP HAS CHANGED
!       ----------
        IF(debug.GT.0) WRITE(lu,*) ' - ACTUALIZING BOUNDARY CONDITIONS'
        CALL phbor
        IF(debug.GT.0) WRITE(lu,*) ' - BOUNDARY CONDITIONS ACTUALIZED'
!       ----------
!      ------------------------------------------------------
!     END OF FIRST STEP : NEXT STEP IS COMPUTING AM AND BM
      ENDIF
!
!     =========================================
!                   MATRIX AM
!     =========================================
      IF(debug.GT.0) WRITE(lu,*) ' - PREPARING THE AM MATRIX'
!      WRITE(LU,*) 'MATRIX AM'
!     ---------------------------
!     DIFFUSION MATRIX FOR AM1
!     ---------------------------
!
      CALL os('X=YZ    ', x=t1, y=c, z=cg)
      CALL matrix(am1,'M=N     ','MATDIF          ',ielm,ielm,
     &            1.d0,s,s,s,t1,t1,s,mesh,msk,maskel)
!
!-----------------------------------------------------------------------
!
! PANCHANG, TO BE REVISITED: 7 IS GMRES
!
! THE DIFFUSION MATRIX USED FOR PRECONDITIONING IS STORED
! IF THE METHOD IS THAT OF PANCHANG ET AL. (ISOLVE(1) =7)
!
!     IF (ISOLVE(1).EQ.7) THEN
!
!        CALL OM('M=CN    ',M=AM3,N=AM1,C=1.D0/(RELAX*(2.D0-RELAX)),MESH=MESH)
!
!     ENDIF
!
!-----------------------------------------------------------------------
!
!     -----------------------
!     MASS MATRIX FOR AM1
!     -----------------------
!     (WARNING : CAN'T USE CURRENT AND SECOND ORDER BOTTOM EFFECTS AT THE SAME TIME)
      IF (courant) THEN
        xmul=1.d0
        CALL os( 'X=YZ    ' , x=t1,y=cg,z=c)
        CALL os( 'X=YZ    ' , x=t2 , y=k , z=k)
        CALL os( 'X=XY    ' , x=t1 , y=t2)
        CALL os( 'X=C     ' , x=t2 , c=omega**2)
        IF(iterkn.EQ.0)THEN
          CALL os( 'X=C     ' , x=t3 , c=omega**2)
        ELSE
          CALL os( 'X=YZ    ' , x=t3 , y=wr , z=wr )
        ENDIF
        CALL os( 'X=Y+Z   ' , x=t1 , y=t1 , z=t2)
        CALL os( 'X=Y-Z   ' , x=t1 , y=t1 , z=t3)
        IF (ipentco.GT.(0.5)) THEN
          WRITE(lu,*) 'IT IS NOT POSSIBLE TO USE '
          WRITE(lu,*) 'CURRENT + BOTTOM EFFECTS AT THE SAME TIME'
          WRITE(lu,*) '- FOR CURRENT ALONE, FIX IPENTO=0'
          WRITE(lu,*) '- FOR BOTTOM EFFECTS ALONE, DON T USE CURRENT'
          WRITE(lu,*) '-------------------------'
          WRITE(lu,*) 'THE CODE IS GOING TO STOP'
          WRITE(lu,*) '-------------------------'
          stop
        ENDIF
      ELSE
        CALL os( 'X=Y/Z   ' , x=t1,y=cg,z=c)
        xmul=omega**2
!       SECOND ORDER BOTTOM EFFECTS ?
!          (IPENTCO > 0 --> T1 = T1*(1+F) )
!           0 : NO EFFECT /  1 : GRADIENT / 2 : CURVATURE /  3 : GRADIENT+CURVATURE
        IF ( (ipentco.GT.(0.5)).AND.(ipentco.LT.(3.5)) ) THEN
          CALL pentco(ipentco)
!         WT USED AND TO BE CONSERVED : T3 = 1+F
!         WT USED : T2 T4 T5 T6 T7 T9 T8 T11 T12
          CALL os('X=YZ    ', x=t1, y=t1, z=t3)
        ENDIF
      ENDIF

      CALL matrix(am2,'M=N     ','FMATMA          ', ielm , ielm ,
     &            xmul , t1,s,s,s,s,s,mesh,msk,maskel)
!

!     --------------------------------------------------
!     COMPUTES DIFFUSION MATRIX - MASS MATRIX
!     --------------------------------------------------
      CALL om('M=M+CN  ', m=am1, n=am2, c=-1.d0, mesh=mesh)
!

!
      IF (courant) THEN
!     --------------------------------------------------
!     ADDS CURRENT - CONVECTION MATRIX
!     --------------------------------------------------
!
        CALL matrix(am2,'M=N     ','MAUGUG          ', ielm , ielm ,
     &            1.d0 , s,s,s,uc,vc,s,mesh,msk,maskel)
!
        CALL om('M=M+CN  ', m=am1, n=am2, c=-1.d0, mesh=mesh)
!
      ENDIF
!     --------------------------------
!     ADDS THE BOUNDARY TERM TO AM1
!     --------------------------------
!
!     AM1 et AM2 --> NON-SYMETRIQUES
      CALL om('M=X(M)  ' , m=am1, mesh=mesh)
      CALL om('M=X(M)  ' , m=am2, mesh=mesh)
!
!     ----------------------------
!     BOUNDARY TERM: INCIDENT WAVE
!     ----------------------------
!
      IF (nptfr .GT. 0) THEN
        CALL matrix(mbor,'M=N     ','FMATMA          ',ielmb,ielmb,
     &        -1.d0,bphi1b,s,s,s,s,s,mesh,.true.,mask1)
        CALL om('M=M+N   ' , m=am1 , n=mbor, mesh=mesh)
!
!
        IF(courant) THEN
          CALL matrix(mbor,'M=N     ','MATFGUG         ',ielmb,ielmb,
     &        1.d0,mesh%XSGBOR,mesh%YSGBOR,s,uc,vc,s,
     &        mesh,.true.,mask1)
          CALL om('M=M+N   ', m=am1, n=mbor, mesh=mesh)
        ENDIF
!
      ENDIF
!     ------------------------------
!     BOUNDARY TERM: FREE EXIT
!     ------------------------------
      IF (nptfr .GT. 0) THEN
        CALL matrix(mbor,'M=N     ','FMATMA          ',ielmb,ielmb,
     &        -1.d0,bphi2b,s,s,s,s,s,mesh,.true.,mask2)
        CALL om('M=M+N   ', m=am1, n=mbor, mesh=mesh)
!
        IF(courant) THEN
          CALL matrix(mbor,'M=N     ','MATFGUG         ',ielmb,ielmb,
     &        1.d0,mesh%XSGBOR,mesh%YSGBOR,s,uc,vc,s,
     &        mesh,.true.,mask2)
          CALL om('M=M+N   ', m=am1, n=mbor, mesh=mesh)
        ENDIF
      ENDIF
!     ------------------------------
!     BOUNDARY TERM : SOLID BOUNDARY
!     ------------------------------
      IF (nptfr .GT. 0) THEN
        CALL matrix(mbor,'M=N     ','FMATMA          ',ielmb,ielmb,
     &        -1.d0,bphi3b,s,s,s,s,s,mesh,.true.,mask3)
        CALL om('M=M+N   ', m=am1, n=mbor, mesh=mesh)
      END IF
!     --------------------- ---------
!     BOUNDARY TERM: IMPOSED WAVE
!     ------------------------------
      IF (nptfr .GT. 0) THEN
        CALL matrix(mbor,'M=N     ','FMATMA          ',ielmb,ielmb,
     &        -1.d0,bphi4b,s,s,s,s,s,mesh,.true.,mask4)
        CALL om('M=M+N   ', m=am1, n=mbor, mesh=mesh)
      END IF
!
!     ------------------------------
!     BOUNDARY TERM: INCIDENT POTENTIAL
!     ------------------------------
!
      IF (nptfr .GT. 0) THEN
        CALL matrix(mbor,'M=N     ','FMATMA          ',ielmb,ielmb,
     &        -1.d0,bphi1b,s,s,s,s,s,mesh,.true.,mask5)
        CALL om('M=M+N   ', m=am1, n=mbor, mesh=mesh)
!
        IF(courant) THEN
          CALL matrix(mbor,'M=N     ','MATFGUG         ',ielmb,ielmb,
     &        1.d0,mesh%XSGBOR,mesh%YSGBOR,s,uc,vc,s,
     &        mesh,.true.,mask5)
          CALL om('M=M+N   ', m=am1, n=mbor, mesh=mesh)
        ENDIF
!
      ENDIF
!
      IF(debug.GT.0) WRITE(lu,*) ' - AM MATRIX PREPARED'
!
!     =========================================
!                   SECOND MEMBERS
!     =========================================
      IF(debug.GT.0) WRITE(lu,*) ' - PREPARING SECOND MEMBERS CV1,CV2'
!     WRITE(LU,*) 'CV1 et CV2'
!     ---------------------
!     SECOND MEMBERS : CV1
!     ---------------------
!
      CALL os('X=C     ', x=cv1, c=0.d0)
!     ------------------------------
!     BOUNDARY TERM: INCIDENT WAVE
!     ------------------------------
!       --- CALCUL DE i COS(TETAP) GAMMA
      CALL os('X=CY    ', x=t1,y=tetap,c=degrad)
      CALL os('X=COS(Y)', x=t2,y=t1)
      CALL os('X=YZ    ', x=t3,y=cphi1b,z=t2)
      CALL os('X=C     ', x=t1, c=0.d0)

      IF (nptfr .GT. 0) THEN
        CALL vector(t1,'=','MASVEC          ',ielmb,
     &           -1.d0,t3,s,s,s,s,s,mesh,.true.,mask1)
      END IF
      CALL osdb( 'X=X+Y   ' , cv1 , t1 , sbid , cbid , mesh )

!     --- CALCUL DE GRAD(Gamma).n : REEL
      CALL os('X=Y     ', x=t2,y=cgrx1b)
      CALL os('X=Y     ', x=t3,y=cgry1b)
      CALL os('X=C     ', x=t1, c=0.d0 )
      CALL os('X=C     ', x=t4, c=1.d0 )
      IF (nptfr .GT. 0) THEN
        CALL vector(t1,'=','FLUBDF          ',ielmb,
     &           1.d0,t4,s,s,t2,t3,s,mesh,.true.,mask1)
      END IF
      CALL osdb( 'X=X+Y   ' , cv1 , t1 , sbid , cbid , mesh )
!     ---------------------------------
!     BOUNDARY TERM: INCIDENT POTENTIAL
!     ---------------------------------
!     --- CALCUL DE i COS(TETAP) GAMMA
      CALL os('X=CY    ', x=t1,y=tetap,c=degrad)
      CALL os('X=COS(Y)', x=t2,y=t1)
      CALL os('X=YZ    ', x=t3,y=cphi1b,z=t2)
      CALL os('X=C     ', x=t1,c=0.d0)
      IF (nptfr .GT. 0) THEN
        CALL vector(t1,'=','MASVEC          ',ielmb,
     &           -1.d0,t3,s,s,s,s,s,mesh,.true.,mask5)
      END IF
      CALL osdb( 'X=X+Y   ' , cv1 , t1 , sbid , cbid , mesh )

!     --- CALCUL DE GRAD(Gamma).n : REEL
      CALL os('X=Y     ', x=t2,y=cgrx1b)
      CALL os('X=Y     ', x=t3,y=cgry1b)
      CALL os('X=C     ', x=t1, c=0.d0 )
      CALL os('X=C     ', x=t4, c=1.d0 )
      IF (nptfr .GT. 0) THEN
        CALL vector(t1,'=','FLUBDF          ',ielmb,
     &           1.d0,t4,s,s,t2,t3,s,mesh,.true.,mask5)
      END IF
      CALL osdb( 'X=X+Y   ' , cv1 , t1 , sbid , cbid , mesh )
!     ------------------------------
!     BOUNDARY TERM: FREE EXIT
!     ------------------------------
      IF (nptfr .GT. 0) THEN
        CALL vector(t1,'=','MASVEC          ',ielmb,
     &           1.d0,cphi2b,s,s,s,s,s,mesh,.true.,mask2)
        CALL osdb( 'X=X+Y   ' , cv1 , t1 , sbid , cbid , mesh )
      END IF
!     ------------------------------
!     BOUNDARY TERM: SOLID BOUNDARY
!     ------------------------------
      IF (nptfr .GT. 0) THEN
        CALL vector(t1,'=','MASVEC          ',ielmb,
     &          1.d0,cphi3b,s,s,s,s,s,mesh,.true.,mask3)
        CALL osdb( 'X=X+Y   ' , cv1 , t1 , sbid , cbid , mesh )
      END IF
!     ------------------------------
!     BOUNDARY TERM: IMPOSED WAVE
!     ------------------------------
      IF (nptfr .GT. 0) THEN
        CALL vector(t1,'=','MASVEC          ',ielmb,
     &           1.d0,cphi4b,s,s,s,s,s,mesh,.true.,mask4)
        CALL osdb( 'X=X+Y   ' , cv1 , t1 , sbid , cbid , mesh )
      END IF
!
!     ---------------------
!     SECOND MEMBERS : CV2
!     ---------------------
!
      CALL os('X=C     ', x=cv2, c=0.d0)
!     ------------------------------
!     BOUNDARY TERM: INCIDENT WAVE
!     ------------------------------
!     --- CALCUL DE i COS(TETAP) GAMMA : IMAGINAIRE
      CALL os('X=CY    ' , x=t1,y=tetap,c=degrad)
      CALL os('X=COS(Y)' , x=t2,y=t1)
      CALL os('X=YZ    ' , x=t3,y=dphi1b,z=t2)
      CALL os('X=C     ' , x=t1,c=0.d0)
      IF (nptfr .GT. 0) THEN
        CALL vector(t1,'=','MASVEC          ',ielmb,
     &           -1.d0,t3,s,s,s,s,s,mesh,.true.,mask1)
      END IF
      CALL osdb( 'X=X+Y   ' , cv2 , t1 , sbid , cbid , mesh )

!     --- CALCUL DE GRAD(Gamma).n : IMAGINAIRE
      CALL os('X=Y     ', x=t2,y=dgrx1b)
      CALL os('X=Y     ', x=t3,y=dgry1b)
      CALL os('X=C     ', x=t1,c=0.d0 )
      CALL os('X=C     ', x=t4,c=1.d0 )

      IF (nptfr .GT. 0) THEN
        CALL vector(t1,'=','FLUBDF          ',ielmb,
     &           1.d0,t4,s,s,t2,t3,s,mesh,.true.,mask1)
      END IF
      CALL osdb( 'X=X+Y   ' , cv2 , t1 , sbid , cbid , mesh )
!     ---------------------------------
!     BOUNDARY TERM: INCIDENT POTENTIAL
!     ---------------------------------
!     --- CALCUL DE i COS(TETAP) GAMMA : IMAGINAIRE
      CALL os('X=CY    ', x=t1,y=tetap,c=degrad)
      CALL os('X=COS(Y)', x=t2,y=t1)
      CALL os('X=YZ    ', x=t3,y=dphi1b,z=t2)
      CALL os('X=C     ', x=t1,c=0.d0)
      IF (nptfr .GT. 0) THEN
        CALL vector(t1,'=','MASVEC          ',ielmb,
     &           -1.d0,t3,s,s,s,s,s,mesh,.true.,mask5)
      END IF
      CALL osdb( 'X=X+Y   ' , cv2 , t1 , sbid , cbid , mesh )

!     --- CALCUL DE GRAD(Gamma).n : IMAGINAIRE
      CALL os('X=Y     ', x=t2,y=dgrx1b)
      CALL os('X=Y     ', x=t3,y=dgry1b)
      CALL os('X=C     ', x=t1,c=0.d0)
      CALL os('X=C     ', x=t4,c=1.d0)

      IF (nptfr .GT. 0) THEN
        CALL vector(t1,'=','FLUBDF          ',ielmb,
     &           1.d0,t4,s,s,t2,t3,s,mesh,.true.,mask5)
      END IF
      CALL osdb( 'X=X+Y   ' , cv2 , t1 , sbid , cbid , mesh )
!
!    ------------------------------
!    BOUNDARY TERM: FREE EXIT
!    ------------------------------
      IF (nptfr .GT. 0) THEN
        CALL vector(t1,'=','MASVEC          ',ielmb,
     &           1.d0,dphi2b,s,s,s,s,s,mesh,.true.,mask2)
      END IF
      CALL osdb( 'X=X+Y   ' , cv2 , t1 , sbid , cbid , mesh )
!     ------------------------------
!     BOUNDARY TERM: SOLID BOUNDARY
!     -----------------------------
      IF (nptfr .GT. 0) THEN
        CALL vector(t1,'=','MASVEC          ',ielmb,
     &           1.d0,dphi3b,s,s,s,s,s,mesh,.true.,mask3)
        CALL osdb( 'X=X+Y   ' , cv2 , t1 , sbid , cbid , mesh )
      END IF
!     ------------------------------
!     BOUNDARY TERM: IMPOSED WAVE
!     ------------------------------
      IF (nptfr .GT. 0) THEN
        CALL vector(t1,'=','MASVEC          ',ielmb,
     &           1.d0,dphi4b,s,s,s,s,s,mesh,.true.,mask4)
        CALL osdb( 'X=X+Y   ' , cv2 , t1 , sbid , cbid , mesh )
      END IF
!CP
!          IF (NCSIZE.GT.1) THEN
!           CALL PARCOM(CV2,2,MESH)
!         ENDIF
!CP
      IF(debug.GT.0) WRITE(lu,*) ' - SECOND MEMBERS CV1,CV2 PREPARED'
!
!     =========================================
!                   MATRIX BM
!     =========================================
      IF(debug.GT.0) WRITE(lu,*) ' - PREPARING THE BM MATRIX'
!      WRITE(LU,*) 'MATRIX BM'
!
!     ----------------------------------------------------------
!     COMPUTES THE MATRIX BM1 FOR THE MU VALUES SPECIFIED
!     FOR THE ITERATION 'ITERMU'
!     ----------------------------------------------------------
!
      CALL os('X=YZ    ', x=t1, y=c , z=cg)
      CALL os('X=YZ    ', x=t2, y=k , z=mu)
      CALL os('X=YZ    ', x=t1, y=t1, z=t2)
      CALL matrix(bm1,'M=N     ','FMATMA          ', ielm , ielm ,
     &            1.d0 , t1,s,s,s,s,s,mesh,msk,maskel)
!
      IF ((courant).AND.(iterkn.GT.0)) THEN
!     ----------------------------------------------------------
!     ADD TERMS TO BM1 FOR THE CURRENT 2*OMEGA.U TERM
!     ----------------------------------------------------------
!          ON DESYMETRISE BM1
        CALL om('M=X(M)  ', m=bm1, mesh=mesh)
!
        CALL matrix(bm2,'M=N     ','MATVGR          ',ielm ,ielm ,
     &            2d0*omega , s,s,s,uc,vc,s,
     &            mesh,msk,maskel)
!
        CALL om('M=M+N   ', m=bm1, n=bm2, mesh=mesh)
!
!     ----------------------------------------------------------
!     ADD TERMS TO THE MATRIX BM1 FOR THE CURRENT : DIV(U) TERM
!     ----------------------------------------------------------
!
        CALL vector(t2 , '=' , 'MASBAS          ' , ielm ,
     &            1.d0 , s , s , s , s , s , s ,
     &            mesh , msk  , maskel )
!
        CALL vector(t13 , '=' , 'GRADF          X' , ielm ,
     &            1.d0 , uc , s , s , s , s , s ,
     &            mesh , msk , maskel)
!
        CALL vector(t14 , '=' , 'GRADF          Y' , ielm ,
     &            1.d0 , vc , s , s , s , s , s ,
     &            mesh , msk , maskel)
!
        CALL os('X=Y+Z   ', x=t15 , y=t14 , z=t13)
        CALL os('X=Y/Z   ', x=t16 , y=t15 , z=t2)
!
        CALL matrix(bm2,'M=N     ','FMATMA          ',ielm ,ielm ,
     &            omega , t16,s,s,s,s,s,
     &            mesh,msk,maskel)
!
        CALL om('M=M+N   ', m=bm1, n=bm2, mesh=mesh)
      ENDIF


!     -------------------------------------------
!     ADDS THE BOUNDARY TERM TO BM1
!     -------------------------------------------
!
      IF (nptfr .GT. 0) THEN
!        ------------------------------
!        BOUNDARY TERM: INCIDENT WAVE
!        ------------------------------
        CALL matrix(mbor,'M=N     ','FMATMA          ',ielmb,ielmb,
     &        -1.d0,aphi1b,s,s,s,s,s,mesh,.true.,mask1)
!         WRITE (*,*)'MBOR = ', MBOR%TYPEXT
        CALL om('M=M+N   ', m=bm1, n=mbor, mesh=mesh)
      END IF
!        ------------------------------
!        BOUNDARY TERM: FREE EXIT
!        ------------------------------
      IF (nptfr .GT. 0) THEN
        CALL matrix(mbor,'M=N     ','FMATMA          ',ielmb,ielmb,
     &           -1.d0,aphi2b,s,s,s,s,s,mesh,.true.,mask2)
        CALL om('M=M+N   ', m=bm1, n=mbor, mesh=mesh)
      END IF
!        ------------------------------
!        BOUNDARY TERM: SOLID BOUNDARY
!        ------------------------------
      IF (nptfr .GT. 0) THEN
        CALL matrix(mbor,'M=N     ','FMATMA          ',ielmb,ielmb,
     &        -1.d0,aphi3b,s,s,s,s,s,mesh,.true.,mask3)
        CALL om('M=M+N   ', m=bm1, n=mbor, mesh=mesh)
      END IF
!        ------------------------------
!        BOUNDARY TERM: IMPOSED WAVE
!        ------------------------------
      IF (nptfr .GT. 0) THEN
        CALL matrix(mbor,'M=N     ','FMATMA          ',ielmb,ielmb,
     &           -1.d0,aphi4b,s,s,s,s,s,mesh,.true.,mask4)
        CALL om('M=M+N   ', m=bm1, n=mbor, mesh=mesh)
      END IF
!        ------------------------------
!        BOUNDARY TERM: INCIDENT POTENTIAL
!        ------------------------------
      IF (nptfr .GT. 0) THEN
        CALL matrix(mbor,'M=N     ','FMATMA          ',ielmb,ielmb,
     &        -1.d0,aphi1b,s,s,s,s,s,mesh,.true.,mask5)
        CALL om('M=M+N   ', m=bm1, n=mbor, mesh=mesh)
      END IF

!     ---------
!     AM2 = AM1
!     ---------
!
      CALL om('M=N     ', m=am2, n=am1, mesh=mesh)
!
!     --------------------------
!     BM1 BECOMES NONSYMMETRICAL
!     --------------------------
!
!     SI ON N'A PAS DE COURANT OU SI ON NE L'A PAS ENCORE PRIS EN COMPTE
!
!      WRITE(LU,*) 'MATRIX BM non sym'
      IF ((.NOT.courant).OR.iterkn.EQ.0) THEN
        CALL om('M=X(M)  ', m=bm1, mesh=mesh)
      ENDIF

!
!     ----------------------------
!     TRIES MASS-LUMPING OF BM1
!     ----------------------------
!
!     MASLU = 1.D0
!     CALL LUMP(T1,BM1,MESH,XMESH,MASLU,MSK,MASKEL)
!     CALL OM('M=CM    ', M=BM1, C=1.D0-MASLU, MESH=MESH)
!     CALL OM('M=M+D   ', M=BM1, D=T1, MESH=MESH)
!
!     ----------
!     BM2 = -BM1
!     ----------
!
      CALL om('M=CN    ', m=bm2, n=bm1, c=-1.d0, mesh=mesh)
!
      IF(debug.GT.0) WRITE(lu,*) ' - BM MATRIX PREPARED'
!
!     =======================================
!
!     TAKES INTO ACCOUNT DIRICHLET POINTS
!
!     =======================================
!
      IF ((.NOT.alemon).AND.(.NOT.alemul)) THEN
        IF (deferl .OR. frotte) THEN
          WRITE(lu,221) itermu+1
        ENDIF
      ENDIF
 221  FORMAT(/,1x,'SUB-ITERATION NUMBER :',1x,i3,/)
!
      IF(debug.GT.0) WRITE(lu,*) ' - CALLING DIRICH'
      CALL dirich(unk,mat,rhs,phib,lidir%I,tb,mesh,kent,msk,maskel)
      IF(debug.GT.0) WRITE(lu,*) ' - DIRICH CALLED'
!
!     ===============================================================
!
!     INHIBITS POSSIBLE DIAGONAL PRECONDITIONING
!     IF AN ELEMENT OF DAM1 IS NEGATIVE OR NULL
!
!     ===============================================================
!
      IF(debug.GT.0) WRITE(lu,*) ' - CALLING CNTPRE'
      tmp = slvart%PRECON
      CALL cntpre(am1%D%R,npoin,tmp,slvart%PRECON)
      IF(debug.GT.0) WRITE(lu,*) ' - CNTPRE CALLED'
!      WRITE(LU,231) SLVART%PRECON
! 231  FORMAT(/,1X,'PRECONDITIONNING AFTER CONTROL :',1X,I3)
!
!     ==========================================================
!
!     PRECONDITIONING BLOCK-DIAGONAL:
!                 THE MATRICES BECOME NONSYMMETRICAL.
!
!     ==========================================================
!
      IF (3*(slvart%PRECON/3).EQ.slvart%PRECON) THEN
        CALL om('M=X(M)  ', m=am1, mesh=mesh)
        CALL om('M=X(M)  ', m=am2, mesh=mesh)
      ENDIF
!
!     ==============================
!
!     SOLVES THE LINEAR SYSTEM
!
!     ==============================
!
!     ----------------------------
!     INITIALISES THE UNKNOWN
!     ----------------------------
!
      IF(itermu.EQ.0.AND.lf.EQ.0) THEN
        CALL lump(t1,am1,mesh,1.d0)
        CALL os('X=Y/Z   ', x=phir, y=cv1, z=t1)
        CALL lump(t1,am2,mesh,1.d0)
        CALL os('X=Y/Z   ', x=phii, y=cv2, z=t1)
      ENDIF
!
      WRITE(lu,241)
 241  FORMAT(/,1x,'LINEAR SYSTEM SOLVING (SOLVE)',/)
!
      IF(debug.GT.0) WRITE(lu,*) ' - SOLVING THE LINEAR SYSTEM'
      CALL solve(unk,mat,rhs,tb,slvart,infogr,mesh,am3)
      IF(debug.GT.0) WRITE(lu,*) ' - LINEAR SYSTEM SOLVED'
!
!     ============================================================
!     DIRECTION LOOP
!      - WAVE-CURRENT :checks convergence on the wave vector
!      - AUTO ANGLES  :checks convergence on TETAP
!     ============================================================
      IF (courant.OR.langauto) THEN
!      ------------------------------------------------------
!      COMPUTE WAVE INCIDENCE USING SPEED AT THE FREE SURFACE
!
!       -----------
        IF(debug.GT.0) WRITE(lu,*) ' - CALLING CALDIR'
        CALL caldir()
        IF(debug.GT.0) WRITE(lu,*) ' - CALDIR CALLED'
!       -----------
!          --> PHIR,PHII
!          --  T1,T2,T3,T4
!         <--  INCI
!        -- DIRECTION OF VECTOR K : INCI
        CALL os('X=COS(Y)', x=t5,y=inci)
        CALL os('X=SIN(Y)', x=t6,y=inci)
!       T5 = K  COS(INCIDENCE)
        CALL os('X=XY    ' , x=t5 , y=k)
!       T6 = K  SIN(INCIDENCE)
        CALL os('X=XY    ' , x=t6 , y=k)
!       ------------------------------------------------------
!
!       Error Initialisation
        erreur1=0.d0
        erreurt=0.d0
!     --------------------------------------------------
!     CONVERGENCE CRITERION FOR WAVE-CURRENT INTERACTION
        IF (courant) THEN
!       MAX ERROR CALCULATION : NORM( Kn - Kn-1 )/NORM(Kn) < EPSDIR
!       ----------------------
          IF (iterkn.GT.0) THEN
            DO i=1,npoin
              erreur1=max(erreur1,
     &          sqrt((kancx%R(i)-t5%R(i))**2+(kancy%R(i)-t6%R(i))**2)/
     &          sqrt(t5%R(i)**2+t6%R(i)**2)
     &               )
            ENDDO
            WRITE(lu,*) '--------------------------------------------'
            WRITE(lu,*) 'WAVE-CURRENT : DIFF. BETWEEN 2 ITER. =',
     &                   erreur1
            WRITE(lu,*) 'LOOP FOR WAVE-CURRENT : TOLERANCE    =',
     &                   epsdir
          ELSE
            WRITE(lu,*) 'INITIAL LOOP FOR WAVE-CURRENT COMPLETED'
          ENDIF
          WRITE(lu,*) '----------------------------------------------'
!       OLD WAVE VECTOR STORAGE
!       -----------------------
          CALL os( 'X=Y     ' , x=kancx , y=t5 )
          CALL os( 'X=Y     ' , x=kancy , y=t6 )
        ENDIF
!     --------------------------------------------------
!
!     -----------------------------------------------------
!     CONVERGENCE CRITERION FOR TETAP AUTOMATIC CALCULATION
        IF (langauto) THEN
!       STORAGE OF INCIDENCE ANGLE ON THE BOUNDARY IN A TABLE
          DO i=1,nptfr
            ig       = mesh%NBOR%I(i)
            angdir(i)=inci%R(ig)
          ENDDO
!         TETAP COMPUTATION
          IF(debug.GT.0) WRITE(lu,*) ' - CALLING CALTETAP'
          CALL caltetap(teta,
     &                  mesh%XSGBOR%R,mesh%YSGBOR%R,angdir,nptfr)
          IF(debug.GT.0) WRITE(lu,*) ' - CALTETAP CALLED'
!       MAX ERROR CALCULATION : MAX(cos(TETAPnew) - cos(TETAPold)) < EPSTP
!       ----------------------
          IF (iterkn.GT.0) THEN
            DO i=1,nptfr
              tetap%R(i)=teta(i)
!             ADD FACTOR 1-R%P EXP(iALFA) ?  ADD PHI?
              erreurt=max(erreurt,
     &             abs(cos(tetap%R(i)*degrad)-cos(tetapm%R(i)*degrad))
     &                )
            ENDDO
            WRITE(lu,*) '-------------------------------------------'
            WRITE(lu,*) 'AUTO-ANGLES : DIFF. BETWEEN 2 ITER. =',
     &                   erreurt
            WRITE(lu,*) 'LOOP FOR AUTO-ANGLE  : TOLERANCE    =',
     &                   epstp
          ELSE
            DO i=1,nptfr
              tetapm%R(i)=tetap%R(i)
              tetap%R(i) =teta(i)
            ENDDO
            WRITE(lu,*) 'INITIAL LOOP FOR AUTOMATIC ANGLES COMPLETED'
          ENDIF
          WRITE(lu,*) '--------------------------------------------'
        ENDIF
!     -----------------------------------------------------
!
!
!      MAX ERROR FOR N PROC
        IF (ncsize.GT.1) THEN
          erreurt = p_max(erreurt)
          erreur1 = p_max(erreur1)
        END IF
!
!      ----------------------------------------------------
!      CHECK CONVERGENCE FOR DIRECTION LOOP
        IF ( (erreur1.GT.epsdir).OR.(erreurt.GT.epstp)
     &                         .OR.(iterkn.EQ.0)     ) THEN
!         NEW ITERATION
          iterkn = iterkn + 1
          IF (iterkn.LE.nittp) THEN
            GOTO 98
          ELSE
            WRITE(lu,101) iterkn
          ENDIF
        ENDIF
        WRITE(lu,203) iterkn
!       REMISE A 1 DU NOMBRE d'ITERATION SUR LE COURRANT ET LA DIRECTION
        iterkn=1
!       =================================================
!       END OF THE LOOP ON DIRECTIONS AND WAVE NUMBER
!       =================================================
      ENDIF
!    ----------------------------------------------------
!
 101  FORMAT(/,1x,'BERKHO (ARTEMIS): YOU REACHED THE MAXIMUM',
     & 1x,'NUMBER OF SUB-ITERATIONS FOR CURRENT OR TETAP :)',1x,i3)

 203  FORMAT(/,1x,'NUMBER OF SUB-ITERATIONS DIRECTION / CURRENT :',
     &   1x,i3)

!     ============================================================
!
!     COMPUTES THE TOTAL DISSIPATION COEFFICIENT MU_DEFERL + MU_FROTTE
!     FOR REGULAR WAVES
!     DISSIPATION FOR IRREGULAR WAVE IS LOOKED AT INTO ARTEMIS.F)
!     ============================================================
!
      IF (.NOT. alemon .AND. .NOT. alemul) THEN
        IF (deferl .OR. frotte) THEN
          ecrhmu=0d0
!         COMPUTES DISSIPATION COEFFICIENT MU2
          IF(debug.GT.0) WRITE(lu,*) ' - CALLING CALCMU'
          CALL calcmu(itermu)
          IF(debug.GT.0) WRITE(lu,*) ' - CALCMU CALLED'
!              WORK TABLE USED                      : T1,T4
!              WORK TABLE USED AND TO BE CONSERVED  : T3 => QB
!
!         USE RELAXATION METHOD FOR DISSPATION COEFFICIENT MU
          IF(debug.GT.0) WRITE(lu,*) ' - CALLING RELAXMU'
          CALL relaxmu(ecrhmu,modhmu,itermu)
          IF(debug.GT.0) WRITE(lu,*) ' - CALLING RELAXMU'

!              WORK TABLE USED                      : NONE

!         ----------------------------------------------------
!         CHECKS CONVERGENCE ON THE DISSIPATION ITERATIVE LOOP
!         ----------------------------------------------------
          WRITE(lu,*) ' '
          WRITE(lu,*) '----------------------------------------------- '
          IF (ecrhmu.GT.epsdis*modhmu) GOTO 98
!
!         QB STORAGE
!         ----------
          CALL os('X=Y     ', x=qb,y=t3)
!
          WRITE(lu,201) itermu
 201      FORMAT(/,1x,'NUMBER OF SUB-ITERATIONS FOR DISSIPATION:',
     &    1x,i3)
!
        ENDIF
      ENDIF
!
!
!-----------------------------------------------------------------------
!
      RETURN
      END