The TELEMAC-MASCARET system  trunk
qporos.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE qporos
3 ! *****************
4 
5  &( tstot , tsder , f , cg, lt,xk,
6  & nf , ndire , npoin2 , amorp )
7 
8 !***********************************************************************
9 ! TOMAWAC V7P3
10 !***********************************************************************
11 !
12 !brief Takes into account the friction due to porous media
13 !
14 !history THIERRY FOUQUET (EDF - LNHE)
15 !+ 27/03/2018
16 !+ V7P3
17 !
18 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
19 !| AMORP |<->| DUMPING DUE TO POROUS MEDIA
20 !| CG |<--| DISCRETIZED GROUP VELOCITY
21 !| F |-->| DIRECTIONAL SPECTRUM
22 !| LT |-->| TIME STEP
23 !| NF |-->| NUMBER OF FREQUENCIES
24 !| NDIRE |-->| NUMBER OF DIRECTIONS
25 !| NPOIN2 |-->| NUMBER OF POINTS IN 2D MESH
26 !| TSDER |<->| DERIVED PART OF THE SOURCE TERM CONTRIBUTION
27 !| TSTOT |<->| TOTAL PART OF THE SOURCE TERM CONTRIBUTION
28 !| XK |<--| WAVE NUMBER
29 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
30 !
31 
32  USE declarations_tomawac, ONLY: deupi, gravit, freq, depth, x, y
33 !
34  USE interface_tomawac, ex_poros => qporos
36  IMPLICIT NONE
37 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
38  INTEGER, INTENT(IN) :: NF,NDIRE,NPOIN2,LT
39  DOUBLE PRECISION, INTENT(INOUT) :: AMORP(npoin2,nf)
40  DOUBLE PRECISION, INTENT(INOUT) :: TSTOT(npoin2,ndire,nf)
41  DOUBLE PRECISION, INTENT(INOUT) :: TSDER(npoin2,ndire,nf)
42  DOUBLE PRECISION, INTENT(INOUT) :: CG(npoin2,nf)
43  DOUBLE PRECISION, INTENT(INOUT) :: XK(npoin2,nf)
44  DOUBLE PRECISION, INTENT(IN) :: F(npoin2,ndire,nf)
45 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
46 
47  INTEGER N
48 ! THe complexe dispersive function and its derivative
49  COMPLEX(KIND=R8) FONC, DERFONC, SECONDFONC, DERTROIS
50  EXTERNAL fonc, derfonc, secondfonc, dertrois
51 ! ALPHA MIGHT BE READ IN A FILE OR SET BY POLYNOME
52  DOUBLE PRECISION GAMMA, ALPHA, INIT
53  COMPLEX(KIND=R8) X0,X1,X2,X3, PHI, PHI1
54  COMPLEX(KIND=R8) LPRIM,INVLPRIM,LSECOND,LTIERCE
55  COMPLEX(KIND=R8) TEMP
56  COMPLEX(KIND=R8) FPRIM,FSECOND
57  DOUBLE PRECISION TAU,COEFF,CA,S,HW
58  INTEGER IP,JF,JP
59  IF (lt.EQ.1) THEN
60 ! TAUX DE POROSITE
61  tau=0.8d0
62 ! COEFFICIENT DE MASSE AJOUTEE
63  ca=2.5d0
64 ! COEFFICIENT D'AMORTISSEMENT
65  coeff=1.d0
66 ! AUTRE COEFF POROSITE
67  s=1+ca*(1.d0-tau)/tau
68  temp=cmplx(s,1.d0, kind=r8)
69 ! COEFFICIENT DE POROSITE
70  phi=tau*temp/(s**2+coeff**2)
71  alpha=0.8d0
72  phi1=cmplx(1.d0,0.d0,kind=r8)
73  DO ip=1,npoin2
74  IF( x(ip).GE.40.d0.AND.x(ip).LE.50.d0
75  & .AND.y(ip).GE.39.d0 .AND.y(ip).LE.61.d0)
76  & THEN
77  alpha=0.8d0
78  hw=(1-alpha)*depth(ip)
79  DO jf=1,nf
80  gamma=(deupi*freq(jf))**2*depth(ip)/gravit
81  init=1.d0
82  DO n=1,10
83  init=gamma/tanh(init)
84  ENDDO
85  x0=cmplx(init,0.d0, kind=r8)
86 ! CALCULATION OF COMPLEXE WAVE NUMBER ACCORDING TO Liou 2006
87  lprim=derfonc(x0,gamma,phi1,alpha)
88  invlprim=1.d0/lprim
89  fprim=derfonc(x0,gamma,phi,alpha)
90  fsecond=secondfonc(x0,gamma,phi,alpha)
91  lsecond=secondfonc(x0,gamma,phi1,alpha)
92  ltierce=dertrois(x0, gamma, phi1,alpha)
93  x1=-fonc(x0,gamma,phi,alpha)*invlprim
94  x2=-invlprim*(0.5d0*lsecond*x1**2+(fprim-lprim)*x1)
95  x3=-invlprim*(lsecond*x1*x2+(fprim-lprim)*x2
96  & +0.5d0*(fsecond-lsecond)*x1**2+1.d0/6.d0*x1**3*ltierce)
97 ! CALCULATION OF CELERITY FOR THE REAL WATER DEPTH
98  gamma=(deupi*freq(jf))**2*hw/gravit
99  init=deupi*freq(jf)*hw/sqrt(gravit*hw)
100  DO n=1,50
101  init=gamma/tanh(init)
102  ENDDO
103  cg(ip,jf)=0.5*(1+2*init/sinh(2*init))
104  & *deupi*freq(jf)*hw/init
105  xk(ip,jf)=init/depth(ip)
106 
107  amorp(ip,jf)=2.d0*cg(ip,jf)*aimag(x0+x1+x2+x3)
108  ENDDO
109  ELSE
110  DO jf=1,nf
111  amorp(ip,jf)=0.d0
112  ENDDO
113  ENDIF
114  ENDDO
115  ENDIF !LT.LT.1
116  DO ip=1,npoin2
117  DO jf=1,nf
118  DO jp=1,ndire
119  tstot(ip,jp,jf) = tstot(ip,jp,jf)+amorp(ip,jf)*f(ip,jp,jf)
120  tsder(ip,jp,jf) = tsder(ip,jp,jf)+amorp(ip,jf)
121  ENDDO
122  ENDDO
123  ENDDO
124  END
125 
126 !-----------------------------------------------------------------------
127  FUNCTION fonc(X, GAMMA, PHI, ALPHA)
128 !
129  IMPLICIT NONE
130  COMPLEX(KIND(1.D0)) FONC
131  DOUBLE PRECISION, INTENT(IN) :: GAMMA, ALPHA
132  COMPLEX(KIND(1.D0)), INTENT(IN) :: X, PHI
133  COMPLEX(KIND(1.D0)) TANHCPLX
134  EXTERNAL tanhcplx
135 
136  DOUBLE PRECISION BETA
137  beta=1.d0-alpha
138  fonc=gamma-x*tanhcplx(beta*x)
139  & -phi*tanhcplx(alpha*x)*(x-gamma*tanhcplx(beta*x))
140  RETURN
141  END
142 !-----------------------------------------------------------------------
143  FUNCTION derfonc(X, GAMMA, PHI, ALPHA)
144 !
145  IMPLICIT NONE
146  COMPLEX(KIND(1.D0)) DERFONC
147  DOUBLE PRECISION, INTENT(IN) :: GAMMA, ALPHA
148  COMPLEX(KIND(1.D0)), INTENT(IN) :: X, PHI
149  COMPLEX(KIND(1.D0)) TANHCPLX,COSHCPLX
150  EXTERNAL tanhcplx,coshcplx
151 
152  DOUBLE PRECISION BETA
153  beta=1.d0-alpha
154  derfonc=-phi*tanhcplx(alpha*x)
155  & *(1.d0-gamma*beta/coshcplx(beta*x)**2)
156  & -phi*alpha*(x-gamma*tanhcplx(beta*x))/coshcplx(alpha*x)**2
157  & -beta*x/coshcplx(beta*x)**2-tanhcplx(beta*x)
158  RETURN
159  END
160 !-----------------------------------------------------------------------
161  FUNCTION secondfonc(X, GAMMA, PHI, ALPHA)
162 !
163  IMPLICIT NONE
164  COMPLEX(KIND(1.D0)) SECONDFONC
165  DOUBLE PRECISION, INTENT(IN) :: GAMMA, ALPHA
166  COMPLEX(KIND(1.D0)), INTENT(IN) :: X, PHI
167  COMPLEX(KIND(1.D0)) TANHCPLX,COSHCPLX
168  EXTERNAL tanhcplx,coshcplx
169 
170  DOUBLE PRECISION BETA
171  beta=1-alpha
172 
173  secondfonc=
174  & -2.d0*phi*alpha*(1.d0-beta*gamma/coshcplx(beta*x)**2)
175  & /coshcplx(alpha*x)**2
176  & -2.d0*phi*gamma*beta**2*tanhcplx(alpha*x)*tanhcplx(beta*x)
177  & /coshcplx(beta*x)**2
178  & +2.d0*alpha**2*phi*(x-gamma*tanhcplx(beta*x))*tanhcplx(alpha*x)
179  & /coshcplx(alpha*x)**2
180  & +(2.d0*x*beta**2*tanhcplx(beta*x)-2.d0*beta)/coshcplx(beta*x)**2
181 
182  RETURN
183  END
184 !-----------------------------------------------------------------------
185  FUNCTION dertrois(X, GAMMA, PHI, ALPHA)
186 !
187  IMPLICIT NONE
188  COMPLEX(KIND(1.D0)) DERTROIS
189  DOUBLE PRECISION, INTENT(IN) :: GAMMA, ALPHA
190  COMPLEX(KIND(1.D0)), INTENT(IN) :: X, PHI
191  COMPLEX(KIND(1.D0)) TANHCPLX,COSHCPLX
192  EXTERNAL tanhcplx,coshcplx
193 
194  DOUBLE PRECISION BETA
195  beta=1.d0-alpha
196 
197  dertrois=
198  & -4.d0*phi*alpha*beta*gamma*(beta*tanhcplx(beta*x)
199  & +alpha*tanhcplx(alpha*x))
200  & /(coshcplx(alpha*x)**2*coshcplx(beta*x)**2)
201  & +4.d0*phi*alpha**2*tanhcplx(alpha*x)/(coshcplx(alpha*x)**2)
202  & -2.d0*phi*gamma*beta**2/(coshcplx(beta*x)**2)*
203  & ( alpha*tanhcplx(beta*x)/(coshcplx(alpha*x)**2)
204  & +beta*tanhcplx(alpha*x)/(coshcplx(beta*x)**2)
205  & -2.d0*beta*tanhcplx(alpha*x)*tanhcplx(beta*x)**2 )
206  & +2.d0*phi*alpha**2/(coshcplx(alpha*x)**2)*
207  & ( (1.d0-gamma*beta/(coshcplx(beta*x)**2))*tanhcplx(alpha*x)
208  & +(x-gamma*tanhcplx(beta*x))*
209  & (alpha/(coshcplx(alpha*x)**2)
210  & -2.d0*alpha*tanhcplx(alpha*x)**2) )
211  & +(2.d0*x*beta**3/(coshcplx(beta*x)**2)
212  & +6.d0*beta**2*tanhcplx(beta*x)
213  & -4.d0*x*beta**3*tanhcplx(beta*x)**2
214  & )/(coshcplx(beta*x)**2)
215  RETURN
216  END
217 !-----------------------------------------------------------------------
218  FUNCTION tanhcplx(Z)
219 !
221  IMPLICIT NONE
222  COMPLEX(KIND(1.D0)) TANHCPLX,III
223  COMPLEX(KIND(1.D0)), INTENT(IN) :: Z
224 
225  iii = cmplx(0.d0,1.d0,kind=r8)
226 
227  tanhcplx= (exp(iii*z)+exp(-iii*z))/(exp(iii*z)-exp(-iii*z))
228  RETURN
229  END
230 !-----------------------------------------------------------------------
231  FUNCTION coshcplx(Z)
232 !
234  IMPLICIT NONE
235  COMPLEX(KIND(1.D0)) COSHCPLX,III
236  COMPLEX(KIND(1.D0)), INTENT(IN) :: Z
237 
238  iii = cmplx(0.d0,1.d0,kind=r8)
239 
240  coshcplx= (exp(iii*z)+exp(-iii*z))*0.5d0
241  RETURN
242  END
complex(kind(1.d0)) function coshcplx(Z)
Definition: qporos.f:233
integer, parameter r8
complex(kind(1.d0)) function dertrois(X, GAMMA, PHI, ALPHA)
Definition: qporos.f:187
double precision, dimension(:), pointer depth
double precision, dimension(:), pointer freq
double precision, dimension(:), pointer y
complex(kind(1.d0)) function fonc(X, GAMMA, PHI, ALPHA)
Definition: qporos.f:129
subroutine qporos
Definition: qporos.f:4
complex(kind(1.d0)) function tanhcplx(Z)
Definition: qporos.f:220
complex(kind(1.d0)) function derfonc(X, GAMMA, PHI, ALPHA)
Definition: qporos.f:145
double precision, dimension(:), pointer x
complex(kind(1.d0)) function secondfonc(X, GAMMA, PHI, ALPHA)
Definition: qporos.f:163