fric3d.f

Go to the documentation of this file.

!#######################################################################
                        SUBROUTINE fric3d
     &  (cfwc, npoin2, dirhou, u_tel, v_tel, uwbm)
!     MODIFIED FRICTION COEFFICIENT DUE TO WAVES+ CURRENTS
!      (CHRISTOFFERSEN AND JONSSON THEORY)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| CFWC           |<--| FRICTION COEFFICIENT DUE TO WAV CURRENT
!| DIRHOU         |-->| WAVE DIRECTION ?????                        
!| NPOIN2         |-->| NUMBER OF POINTS IN 2D MESH
!| U_TEL          |-->| VELOCITY ALONG X 
!| V_TEL          |-->| VELOCITY ALONG Y 
!| UWBM           |-->| VELOCITY ON BOTTOM
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      USE bief
      USE declarations_tomawac, ONLY : depth, deupi , gravit, nf, freq 
      USE interface_tomawac, ex_fric3d => fric3d

      IMPLICIT NONE
      
!.....VARIABLES IN ARGUMENT
!     """"""""""""""""""""
      INTEGER, INTENT(IN) :: NPOIN2
      TYPE(bief_obj),   INTENT(IN) :: U_TEL,V_TEL
      TYPE(bief_obj),   INTENT(INOUT)   :: CFWC
      DOUBLE PRECISION, INTENT(IN) :: UWBM(npoin2)
      DOUBLE PRECISION, INTENT(IN) :: DIRHOU(npoin2)
!.....LOCAL VARIABLES
!     """""""""""""""""
      DOUBLE PRECISION VITCOU,DIRCOU
      DOUBLE PRECISION OMEGAA
      DOUBLE PRECISION XKN   ,  FW, FC, XKAPPA
      INTEGER          MODEL , LUMES
      INTEGER          JF    , I
      INTEGER          ITERM , ITER  , KONVER
      DOUBLE PRECISION TOLF  , TOLX  , ERRF  , ERRX  , F1MJ, F2MJ,
     &                 df1dfc, df1dfw, df2dfc, df2dfw, dfc, dfw,
     &                 coef1 , fcinit, fwinit, det, xj,
     &                 betamj  , zeta  , vitnul
!
!.....VARIABLES LOCALES
!     """""""""""""""""
      DOUBLE PRECISION R, COSDIR, COEFJ ,
     &                 mmm   , mmmdfc, mmmdfw, jjj   , jjjdfc, jjjdfw,
     &                 sig   , sigdfc, sigdfw, xka   , xkadfc, xkadfw,
     &                 taw   , tawdfc, tawdfw, coefta, coefka, coef2 ,
     &                 auxi  , coefn 

!.....VARIABLES LOCALES-FWMOD2

      DOUBLE PRECISION USKA, C1, YO
      DOUBLE PRECISION YN , RAC2, C2,G

      g = sqrt(gravit)
      lumes = 0
!.....PARAMETRES DE L'ALGORITHME ITERATIF
!     """""""""""""""""""""""""""""""""""
!!!!!!!!!!!!!!!!!!!!!!!!!!comeca da subrotina chris2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      vitnul= 1.d-6
      iterm = 100
      tolf  = 1.d-8
      tolx  = 1.d-6
      xkappa = 0.40d0
      betamj   = 0.0747d0
      zeta   = 0.5013d0
      rac2=dsqrt(2.d0)
      c1  =1.d0/(xkappa*rac2)
!""""""""""""""""""""""""""""""""""""""""""""
!  0  C AFFECTATION DES CONSTANTES DU MODELE.
!=====C======================================

      r      = 0.45d0
      coefn  = 0.367d0
      DO i=1,npoin2
        vitcou=sqrt(u_tel%R(i)**2.d0+v_tel%R(i)**2.d0)
        IF(vitcou.GE.vitnul) THEN 
          dircou=atan2(u_tel%R(i),v_tel%R(i))
        ELSE 
          dircou=0
        ENDIF
        xkn = 11.d0*max(depth(i),1.d-9) /exp(0.41d0*73.d0/g)
        DO jf = 1,nf
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!ATENÇAO!TOU EU A DAR XKN
          omegaa=deupi*freq(jf)
!.....DETERMINATION DES SOLUTIONS INITIALES (SANS EFFET MUTUEL)
!     """""""""""""""""""""""""""""""""""""""""""""""""""""""""
          xj=(dsqrt(betamj)*uwbm(i)/(xkn*omegaa))**(2.d0/3.d0)
          coef1  = 30.d0/dexp(1.d0)
          fcinit = 2.d0*(xkappa/dlog(coef1*depth(i)/xkn))**2
          IF (uwbm(i).LT.vitnul.OR.xj.LT.vitnul) THEN
            fwinit=1.d-20
          ELSEIF (xj.LT.3.47d0) THEN
            fwinit=2.d0*betamj/xj
          ELSE
!     FUNCAO FWMOD2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
            uska=uwbm(i)/(xkn*omegaa)
!     
            c2=c1*dlog(30.d0*xkappa*dexp(-2.d0*zeta)*uska/rac2)
            yo=(uska/betamj)**(1.d0/3.d0)/rac2
            DO iter=1,iterm
              yn=c2-c1*dlog(yo)
              IF (dabs(yn-yo).LT.1.d-5) GOTO 205
              yo=yn
            ENDDO
            WRITE(6,*) '/!/ STOP IN  FWSEUL - NO  CONVERGENCE'
            WRITE(6,*) 'Uwbm/(Kn.Omega) = ', uska,'XKN,OMEGAA',
     &           xkn,omegaa,'I',i
            stop
 205        CONTINUE
            fwinit=1.d0/(yn*yn)
!FUNCAO final FWMOD2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          ENDIF
!
!.....AFFECTATION DE LA SOLUTION DE DEPART
!     """"""""""""""""""""""""""""""""""""
          fc=fcinit
          fw=fwinit
!.....ON S'ARRETE A CES SOLUTIONS SI UNE DES DEUX VITESSES EST
!     NULLE OU SI LE RAPPORT DES VITESSES EST TRES GRAND OU PETIT.
!     """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""
          IF ((vitcou.LT.vitnul).OR.(uwbm(i).LT.vitnul)
     &         .OR.(vitcou/uwbm(i).GT.1.d3).OR.
     &         (uwbm(i)/vitcou.GT.1.d3)) THEN
            cfwc%R(i)=fc
          ELSE
! 
!.....ALGORITHME ITERATIF DE NEWTON-RAPHSON.
!     """"""""""""""""""""""""""""""""""""""
            iter=0
 500        CONTINUE
            iter=iter+1
!!!!!!!!!!!!!!!!!!!!!!!!!!final da subrotina chris2!!!!!!!!!!!!!!!
!.......CALCUL DES MATRICES ALFA ET BETA POUR X DONNE.
!CJSUBR!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!=====C
!  1  C CALCUL DE SIGMA  ET DE SES DERIVEES (COMMUN AUX 2 MODELES)
!=====C===========================================================
            sig    =  (fc/fw)*(vitcou/uwbm(i))**2
            sigdfc =  sig/fc
            sigdfw = -sig/fw
!=====C
!  2  C CALCUL DE M      ET DE SES DERIVEES (COMMUN AUX 2 MODELES)
!=====C===========================================================
            cosdir = dabs(cos(dirhou(i)-dircou))
            mmm    = dsqrt(1.d0+sig*sig+2.d0*sig*cosdir)
            mmmdfc = (sig+cosdir)/mmm*sigdfc
            mmmdfw = (sig+cosdir)/mmm*sigdfw     
!=====C
!  3  C CALCUL DE J      ET DE SES DERIVEES (COMMUN AUX 2 MODELES)
!=====C===========================================================
            coefj  = uwbm(i)/(xkn*omegaa*dsqrt(2.d0))
            jjj    = coefj*dsqrt(mmm*fw)
            jjjdfc = 0.5d0*jjj/mmm*mmmdfc
            jjjdfw = 0.5d0*jjj*(1.d0/fw+1.d0/mmm*mmmdfw)
!=====C
!  4  C SELECTION DU MODELE EN FONCTION DE LA VALEUR DE J
!=====C==================================================
            IF (jjj.GT.3.47d0) THEN
              model=2
            ELSE
              model=1
            ENDIF
!=====C
!  5  C CALCUL DE DELTAW ET DE SES DERIVEES (SELON LE MODELE CHOISI)
!=====C=============================================================
            IF (model.EQ.1) THEN
              coefta = xkn*r*deupi/2.d0*dsqrt(betamj*0.5d0)
              taw    = coefta*dsqrt(jjj)
              tawdfc = 0.5d0*taw/jjj*jjjdfc
              tawdfw = 0.5d0*taw/jjj*jjjdfw
            ELSE
              coefta = xkn*coefn*xkappa
              taw    = coefta*jjj
              tawdfc = coefta*jjjdfc
              tawdfw = coefta*jjjdfw
            ENDIF
!=====C
!  6  C CALCUL DE KA     ET DE SES DERIVEES (SELON LE MODELE CHOISI)
!=====C=============================================================
            IF (model.EQ.1) THEN
              coefka = xkappa/(betamj*xkn)
              xka    = 30.d0*taw*dexp(-coefka*taw*dsqrt(sig/mmm))
              xkadfc = xka*(tawdfc/taw-coefka*taw*dsqrt(sig/mmm)
     &             *(tawdfc/taw+0.5d0*sigdfc/sig-0.5d0*mmmdfc/mmm))
              xkadfw = xka*(tawdfw/taw-coefka*taw*dsqrt(sig/mmm)
     &             *(tawdfw/taw+0.5d0*sigdfw/sig-0.5d0*mmmdfw/mmm))
            ELSE
              auxi   = dsqrt(sig/mmm)
              xka    = xkn*(30.d0*taw/xkn)**(1.d0-auxi)
              xkadfc = xka*( -0.5d0/auxi*(sigdfc-sig*mmmdfc/mmm)/mmm
     &             *dlog(30.d0*taw/xkn) + (1.d0-auxi)*tawdfc/taw )
              xkadfw = xka*( -0.5d0/auxi*(sigdfw-sig*mmmdfw/mmm)/mmm
     &             *dlog(30.d0*taw/xkn) + (1.d0-auxi)*tawdfw/taw )
            ENDIF
!=====C
!  7  C CALCUL DE F1     ET DE SES DERIVEES (SELON LE MODELE CHOISI)
!=====C=============================================================
            IF (model.EQ.1) THEN
              f1mj   = fw -2.d0*betamj*mmm/jjj
              df1dfc =      -2.d0*betamj*(mmmdfc/jjj-jjjdfc*mmm/jjj**2)
              df1dfw = 1.d0 -2.d0*betamj*(mmmdfw/jjj-jjjdfw*mmm/jjj**2)
            ELSE
              coef1  = 30.d0*xkappa*dexp(-2.d0*zeta)
              auxi   = dlog(coef1*jjj)
              f1mj   = fw  - 2.d0*mmm*(xkappa/auxi)**2
              df1dfc =     - 2.d0*xkappa**2*(mmmdfc/auxi**2 
     &             - 2.d0*mmm/jjj*jjjdfc/auxi**3)
              df1dfw = 1.d0 - 2.d0*xkappa**2*(mmmdfw/auxi**2 
     &             - 2.d0*mmm/jjj*jjjdfw/auxi**3)
            ENDIF
!=====C
!  8  C CALCUL DE F2     ET DE SES DERIVEES (COMMUN AUX 2 MODELES)
!=====C===========================================================
            coef2  = 30.d0*depth(i)/dexp(1.d0)
            auxi   = dlog(coef2/xka)
            f2mj     = dsqrt(2.d0/fc) - auxi/xkappa
            df2dfc = xkadfc/(xka*xkappa)-1.d0/dsqrt(2.d0*fc**3)
            df2dfw = xkadfw/(xka*xkappa)

!CJSUBR!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!.......CALCUL DE L'ERREUR SUR LA FONCTION
!       """"""""""""""""""""""""""""""""""
            errf=abs(f1mj)+abs(f2mj)
            IF (errf.LT.tolf) THEN
              konver=1
              GOTO 201
            ENDIF
!
!.......CALCUL DE DX
!       """"""""""""
            det=df1dfc*df2dfw-df1dfw*df2dfc
            IF (abs(det).LT.1.d-10) THEN
              WRITE(6,*) '/!/ ARRET DANS CJSUBR : DETERMINAT NUL DET='
              stop
            ENDIF
            dfc=(-f1mj*df2dfw+f2mj*df1dfw)/det
            dfw=(-f2mj*df1dfc+f1mj*df2dfc)/det
!     
!.......MISE A JOUR DE LA SOLUTION
!       """"""""""""""""""""""""""
            fc=fc+dfc
            fw=fw+dfw
!     
!.......TEST SUR LA VALEUR DE DX
!       """"""""""""""""""""""""
            errx=abs(dfc)+abs(dfw)
            IF (errx.LT.tolx) THEN
              konver=2
              GOTO 201
            ENDIF
!
            IF (iter.LT.iterm) GOTO 500
!
            konver=0
!
 201        CONTINUE
            IF (lumes.GT.0) WRITE(lumes,2000) konver,iter,fc,fw
 2000       FORMAT(
     &     'CONVERGENCE ',i1,' APRES',i5,' ITERATIONS   =>  FC =',
     &     e11.4,'  ET FW =',e11.4)
          ENDIF
        ENDDO
      ENDDO
      RETURN
      END