The TELEMAC-MASCARET system  trunk
fric3d.f
Go to the documentation of this file.
1 !#######################################################################
2  SUBROUTINE fric3d
3  & (cfwc, npoin2, dirhou, u_tel, v_tel, uwbm)
4 ! MODIFIED FRICTION COEFFICIENT DUE TO WAVES+ CURRENTS
5 ! (CHRISTOFFERSEN AND JONSSON THEORY)
6 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
7 !| CFWC |<--| FRICTION COEFFICIENT DUE TO WAV CURRENT
8 !| DIRHOU |-->| WAVE DIRECTION ?????
9 !| NPOIN2 |-->| NUMBER OF POINTS IN 2D MESH
10 !| U_TEL |-->| VELOCITY ALONG X
11 !| V_TEL |-->| VELOCITY ALONG Y
12 !| UWBM |-->| VELOCITY ON BOTTOM
13 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
14  USE bief
15  USE declarations_tomawac, ONLY : depth, deupi , gravit, nf, freq
16  USE interface_tomawac, ex_fric3d => fric3d
17 
18  IMPLICIT NONE
19 
20 !.....VARIABLES IN ARGUMENT
21 ! """"""""""""""""""""
22  INTEGER, INTENT(IN) :: NPOIN2
23  TYPE(bief_obj), INTENT(IN) :: U_TEL,V_TEL
24  TYPE(bief_obj), INTENT(INOUT) :: CFWC
25  DOUBLE PRECISION, INTENT(IN) :: UWBM(npoin2)
26  DOUBLE PRECISION, INTENT(IN) :: DIRHOU(npoin2)
27 !.....LOCAL VARIABLES
28 ! """""""""""""""""
29  DOUBLE PRECISION VITCOU,DIRCOU
30  DOUBLE PRECISION OMEGAA
31  DOUBLE PRECISION XKN , FW, FC, XKAPPA
32  INTEGER MODEL , LUMES
33  INTEGER JF , I
34  INTEGER ITERM , ITER , KONVER
35  DOUBLE PRECISION TOLF , TOLX , ERRF , ERRX , F1MJ, F2MJ,
36  & df1dfc, df1dfw, df2dfc, df2dfw, dfc, dfw,
37  & coef1 , fcinit, fwinit, det, xj,
38  & betamj , zeta , vitnul
39 !
40 !.....VARIABLES LOCALES
41 ! """""""""""""""""
42  DOUBLE PRECISION R, COSDIR, COEFJ ,
43  & mmm , mmmdfc, mmmdfw, jjj , jjjdfc, jjjdfw,
44  & sig , sigdfc, sigdfw, xka , xkadfc, xkadfw,
45  & taw , tawdfc, tawdfw, coefta, coefka, coef2 ,
46  & auxi , coefn
47 
48 !.....VARIABLES LOCALES-FWMOD2
49 
50  DOUBLE PRECISION USKA, C1, YO
51  DOUBLE PRECISION YN , RAC2, C2,G
52 
53  g = sqrt(gravit)
54  lumes = 0
55 !.....PARAMETRES DE L'ALGORITHME ITERATIF
56 ! """""""""""""""""""""""""""""""""""
57 !!!!!!!!!!!!!!!!!!!!!!!!!!comeca da subrotina chris2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
58  vitnul= 1.d-6
59  iterm = 100
60  tolf = 1.d-8
61  tolx = 1.d-6
62  xkappa = 0.40d0
63  betamj = 0.0747d0
64  zeta = 0.5013d0
65  rac2=dsqrt(2.d0)
66  c1 =1.d0/(xkappa*rac2)
67 !""""""""""""""""""""""""""""""""""""""""""""
68 ! 0 C AFFECTATION DES CONSTANTES DU MODELE.
69 !=====C======================================
70 
71  r = 0.45d0
72  coefn = 0.367d0
73  DO i=1,npoin2
74  vitcou=sqrt(u_tel%R(i)**2.d0+v_tel%R(i)**2.d0)
75  IF(vitcou.GE.vitnul) THEN
76  dircou=atan2(u_tel%R(i),v_tel%R(i))
77  ELSE
78  dircou=0
79  ENDIF
80  xkn = 11.d0*max(depth(i),1.d-9) /exp(0.41d0*73.d0/g)
81  DO jf = 1,nf
82 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!ATENÇAO!TOU EU A DAR XKN
83  omegaa=deupi*freq(jf)
84 !.....DETERMINATION DES SOLUTIONS INITIALES (SANS EFFET MUTUEL)
85 ! """""""""""""""""""""""""""""""""""""""""""""""""""""""""
86  xj=(dsqrt(betamj)*uwbm(i)/(xkn*omegaa))**(2.d0/3.d0)
87  coef1 = 30.d0/dexp(1.d0)
88  fcinit = 2.d0*(xkappa/dlog(coef1*depth(i)/xkn))**2
89  IF (uwbm(i).LT.vitnul.OR.xj.LT.vitnul) THEN
90  fwinit=1.d-20
91  ELSEIF (xj.LT.3.47d0) THEN
92  fwinit=2.d0*betamj/xj
93  ELSE
94 ! FUNCAO FWMOD2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
95  uska=uwbm(i)/(xkn*omegaa)
96 !
97  c2=c1*dlog(30.d0*xkappa*dexp(-2.d0*zeta)*uska/rac2)
98  yo=(uska/betamj)**(1.d0/3.d0)/rac2
99  DO iter=1,iterm
100  yn=c2-c1*dlog(yo)
101  IF (dabs(yn-yo).LT.1.d-5) GOTO 205
102  yo=yn
103  ENDDO
104  WRITE(6,*) '/!/ STOP IN FWSEUL - NO CONVERGENCE'
105  WRITE(6,*) 'Uwbm/(Kn.Omega) = ', uska,'XKN,OMEGAA',
106  & xkn,omegaa,'I',i
107  stop
108  205 CONTINUE
109  fwinit=1.d0/(yn*yn)
110 !FUNCAO final FWMOD2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
111  ENDIF
112 !
113 !.....AFFECTATION DE LA SOLUTION DE DEPART
114 ! """"""""""""""""""""""""""""""""""""
115  fc=fcinit
116  fw=fwinit
117 !.....ON S'ARRETE A CES SOLUTIONS SI UNE DES DEUX VITESSES EST
118 ! NULLE OU SI LE RAPPORT DES VITESSES EST TRES GRAND OU PETIT.
119 ! """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""
120  IF ((vitcou.LT.vitnul).OR.(uwbm(i).LT.vitnul)
121  & .OR.(vitcou/uwbm(i).GT.1.d3).OR.
122  & (uwbm(i)/vitcou.GT.1.d3)) THEN
123  cfwc%R(i)=fc
124  ELSE
125 !
126 !.....ALGORITHME ITERATIF DE NEWTON-RAPHSON.
127 ! """"""""""""""""""""""""""""""""""""""
128  iter=0
129  500 CONTINUE
130  iter=iter+1
131 !!!!!!!!!!!!!!!!!!!!!!!!!!final da subrotina chris2!!!!!!!!!!!!!!!
132 !.......CALCUL DES MATRICES ALFA ET BETA POUR X DONNE.
133 !CJSUBR!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
134 !
135 !=====C
136 ! 1 C CALCUL DE SIGMA ET DE SES DERIVEES (COMMUN AUX 2 MODELES)
137 !=====C===========================================================
138  sig = (fc/fw)*(vitcou/uwbm(i))**2
139  sigdfc = sig/fc
140  sigdfw = -sig/fw
141 !=====C
142 ! 2 C CALCUL DE M ET DE SES DERIVEES (COMMUN AUX 2 MODELES)
143 !=====C===========================================================
144  cosdir = dabs(cos(dirhou(i)-dircou))
145  mmm = dsqrt(1.d0+sig*sig+2.d0*sig*cosdir)
146  mmmdfc = (sig+cosdir)/mmm*sigdfc
147  mmmdfw = (sig+cosdir)/mmm*sigdfw
148 !=====C
149 ! 3 C CALCUL DE J ET DE SES DERIVEES (COMMUN AUX 2 MODELES)
150 !=====C===========================================================
151  coefj = uwbm(i)/(xkn*omegaa*dsqrt(2.d0))
152  jjj = coefj*dsqrt(mmm*fw)
153  jjjdfc = 0.5d0*jjj/mmm*mmmdfc
154  jjjdfw = 0.5d0*jjj*(1.d0/fw+1.d0/mmm*mmmdfw)
155 !=====C
156 ! 4 C SELECTION DU MODELE EN FONCTION DE LA VALEUR DE J
157 !=====C==================================================
158  IF (jjj.GT.3.47d0) THEN
159  model=2
160  ELSE
161  model=1
162  ENDIF
163 !=====C
164 ! 5 C CALCUL DE DELTAW ET DE SES DERIVEES (SELON LE MODELE CHOISI)
165 !=====C=============================================================
166  IF (model.EQ.1) THEN
167  coefta = xkn*r*deupi/2.d0*dsqrt(betamj*0.5d0)
168  taw = coefta*dsqrt(jjj)
169  tawdfc = 0.5d0*taw/jjj*jjjdfc
170  tawdfw = 0.5d0*taw/jjj*jjjdfw
171  ELSE
172  coefta = xkn*coefn*xkappa
173  taw = coefta*jjj
174  tawdfc = coefta*jjjdfc
175  tawdfw = coefta*jjjdfw
176  ENDIF
177 !=====C
178 ! 6 C CALCUL DE KA ET DE SES DERIVEES (SELON LE MODELE CHOISI)
179 !=====C=============================================================
180  IF (model.EQ.1) THEN
181  coefka = xkappa/(betamj*xkn)
182  xka = 30.d0*taw*dexp(-coefka*taw*dsqrt(sig/mmm))
183  xkadfc = xka*(tawdfc/taw-coefka*taw*dsqrt(sig/mmm)
184  & *(tawdfc/taw+0.5d0*sigdfc/sig-0.5d0*mmmdfc/mmm))
185  xkadfw = xka*(tawdfw/taw-coefka*taw*dsqrt(sig/mmm)
186  & *(tawdfw/taw+0.5d0*sigdfw/sig-0.5d0*mmmdfw/mmm))
187  ELSE
188  auxi = dsqrt(sig/mmm)
189  xka = xkn*(30.d0*taw/xkn)**(1.d0-auxi)
190  xkadfc = xka*( -0.5d0/auxi*(sigdfc-sig*mmmdfc/mmm)/mmm
191  & *dlog(30.d0*taw/xkn) + (1.d0-auxi)*tawdfc/taw )
192  xkadfw = xka*( -0.5d0/auxi*(sigdfw-sig*mmmdfw/mmm)/mmm
193  & *dlog(30.d0*taw/xkn) + (1.d0-auxi)*tawdfw/taw )
194  ENDIF
195 !=====C
196 ! 7 C CALCUL DE F1 ET DE SES DERIVEES (SELON LE MODELE CHOISI)
197 !=====C=============================================================
198  IF (model.EQ.1) THEN
199  f1mj = fw -2.d0*betamj*mmm/jjj
200  df1dfc = -2.d0*betamj*(mmmdfc/jjj-jjjdfc*mmm/jjj**2)
201  df1dfw = 1.d0 -2.d0*betamj*(mmmdfw/jjj-jjjdfw*mmm/jjj**2)
202  ELSE
203  coef1 = 30.d0*xkappa*dexp(-2.d0*zeta)
204  auxi = dlog(coef1*jjj)
205  f1mj = fw - 2.d0*mmm*(xkappa/auxi)**2
206  df1dfc = - 2.d0*xkappa**2*(mmmdfc/auxi**2
207  & - 2.d0*mmm/jjj*jjjdfc/auxi**3)
208  df1dfw = 1.d0 - 2.d0*xkappa**2*(mmmdfw/auxi**2
209  & - 2.d0*mmm/jjj*jjjdfw/auxi**3)
210  ENDIF
211 !=====C
212 ! 8 C CALCUL DE F2 ET DE SES DERIVEES (COMMUN AUX 2 MODELES)
213 !=====C===========================================================
214  coef2 = 30.d0*depth(i)/dexp(1.d0)
215  auxi = dlog(coef2/xka)
216  f2mj = dsqrt(2.d0/fc) - auxi/xkappa
217  df2dfc = xkadfc/(xka*xkappa)-1.d0/dsqrt(2.d0*fc**3)
218  df2dfw = xkadfw/(xka*xkappa)
219 
220 !CJSUBR!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
221 !
222 !.......CALCUL DE L'ERREUR SUR LA FONCTION
223 ! """"""""""""""""""""""""""""""""""
224  errf=abs(f1mj)+abs(f2mj)
225  IF (errf.LT.tolf) THEN
226  konver=1
227  GOTO 201
228  ENDIF
229 !
230 !.......CALCUL DE DX
231 ! """"""""""""
232  det=df1dfc*df2dfw-df1dfw*df2dfc
233  IF (abs(det).LT.1.d-10) THEN
234  WRITE(6,*) '/!/ ARRET DANS CJSUBR : DETERMINAT NUL DET='
235  stop
236  ENDIF
237  dfc=(-f1mj*df2dfw+f2mj*df1dfw)/det
238  dfw=(-f2mj*df1dfc+f1mj*df2dfc)/det
239 !
240 !.......MISE A JOUR DE LA SOLUTION
241 ! """"""""""""""""""""""""""
242  fc=fc+dfc
243  fw=fw+dfw
244 !
245 !.......TEST SUR LA VALEUR DE DX
246 ! """"""""""""""""""""""""
247  errx=abs(dfc)+abs(dfw)
248  IF (errx.LT.tolx) THEN
249  konver=2
250  GOTO 201
251  ENDIF
252 !
253  IF (iter.LT.iterm) GOTO 500
254 !
255  konver=0
256 !
257  201 CONTINUE
258  IF (lumes.GT.0) WRITE(lumes,2000) konver,iter,fc,fw
259  2000 FORMAT(
260  & 'CONVERGENCE ',i1,' APRES',i5,' ITERATIONS => FC =',
261  & e11.4,' ET FW =',e11.4)
262  ENDIF
263  ENDDO
264  ENDDO
265  RETURN
266  END
267 
subroutine fric3d(CFWC, NPOIN2, DIRHOU, U_TEL, V_TEL, UWBM)
Definition: fric3d.f:5
double precision, dimension(:), pointer depth
double precision, dimension(:), pointer freq
Definition: bief.f:3