The TELEMAC-MASCARET system  trunk
conw4d.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE conw4d
3 ! *****************
4 !
5  &(cx,cy,ct,cf,xk,cg, npoin2,ndire,jf,nf)
6 !
7 !***********************************************************************
8 ! TOMAWAC V7P1
9 !***********************************************************************
10 !
11 !brief COMPUTES THE ADVECTION FIELD.
12 !
13 !warning IN THIS CASE THE X AXIS IS VERTICAL ORIENTED UPWARDS AND
14 !+ THE Y AXIS IS HORIZONTAL ORIENTED TOWARDS THE RIGHT;
15 !+ TETA IS THE DIRECTION WRT NORTH, CLOCKWISE
16 !
17 !history F MARCOS (LNH)
18 !+ 01/02/95
19 !+ V1P0
20 !+
21 !
22 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
23 !+ 13/07/2010
24 !+ V6P0
25 !+ Translation of French comments within the FORTRAN sources into
26 !+ English comments
27 !
28 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
29 !+ 21/08/2010
30 !+ V6P0
31 !+ Creation of DOXYGEN tags for automated documentation and
32 !+ cross-referencing of the FORTRAN sources
33 !
34 !history G.MATTAROLO (EDF - LNHE)
35 !+ 14/06/2011
36 !+ V6P1
37 !+ Translation of French names of the variables in argument
38 !
39 !history J-M HERVOUET (EDF-LNHE)
40 !+ 27/11/2012
41 !+ V6P3
42 !+ Optimisation (loops on NPOIN2 and NDIRE swapped to get smaller
43 !+ strides, work array TRA01 differently used, etc.)
44 !
45 !history J-M HERVOUET (EDF-LNHE)
46 !+ 08/01/2014
47 !+ V7P0
48 !+ Was always called with COURAN=.TRUE., with heavy tests on this
49 !+ variable. COURAN now considered=.TRUE. and suppressed.
50 !
51 !history J-M HERVOUET (EDF-LNHE)
52 !+ 16/11/2015
53 !+ V7P1
54 !+ DZHDT always included in formulas, even if depth not varying. This
55 !+ assumes that DZHDT is duly set to 0 in such cases. Testing MAREE
56 !+ was a bug since it is true only with depth given in a file...
57 !
58 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
59 !| CG |-->| DISCRETIZED GROUP VELOCITY
60 !| COSF |-->| COSINE OF THE LATITUDES OF THE POINTS 2D
61 !| COSTET |-->| COSINE OF TETA ANGLE
62 !| CT |<--| ADVECTION FIELD ALONG TETA
63 !| CX |<--| ADVECTION FIELD ALONG X(OR PHI)
64 !| CY |<--| ADVECTION FIELD ALONG Y(OR LAMBDA)
65 !| CF |<--| ADVECTION FIELD ALONG FREQUENCY
66 !| DEPTH |-->| WATER DEPTH
67 !| DVY |-->| DERIVATIVE OF CURRENT SPEED DU/DX
68 !| DVX |-->| DERIVATIVE OF CURRENT SPEED DU/DY
69 !| DUY |-->| DERIVATIVE OF CURRENT SPEED DV/DX
70 !| DUX |-->| DERIVATIVE OF CURRENT SPEED DV/DY
71 !| DZHDT |-->| WATER DEPTH DERIVATIVE WITH RESPECT TO T
72 !| DZY |-->| SEA BOTTOM SLOPE ALONG X
73 !| DZX |-->| SEA BOTTOM SLOPE ALONG Y
74 !| FREQ |-->| DISCRETIZED FREQUENCIES
75 !| JF |-->| INDEX OF THE FREQUENCY
76 !| NF |-->| NUMBER OF FREQUENCIES
77 !| NDIRE |-->| NUMBER OF DIRECTIONS
78 !| NPOIN2 |-->| NUMBER OF POINTS IN 2D MESH
79 !| PROINF |-->| LOGICAL INDICATING INFINITE DEPTH ASSUMPTION
80 !| SINTET |-->| SINE OF TETA ANGLE
81 !| SPHE |-->| LOGICAL INDICATING SPHERICAL COORD ASSUMPTION
82 !| TGF |-->| TANGENT OF THE LATITUDES OF THE POINTS 2D
83 !| U |-->| CURRENT SPEED ALONG X
84 !| V |-->| CURRENT SPEED ALONG Y
85 !| XK |-->| DISCRETIZED WAVE NUMBER
86 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
87 !
89  &proinf, sphe, cosf,tgf,depth,dzhdt,dzy,dzx,freq,costet,sintet,
90  & uc, vc, dux, duy, dvx, dvy, tra01
91 !
93  USE interface_tomawac, ex_conw4d => conw4d
94  IMPLICIT NONE
95 !
96 !
97 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
98 !
99  INTEGER, INTENT(IN) :: NF,NDIRE,NPOIN2,JF
100 !
101  DOUBLE PRECISION, INTENT(INOUT) :: CX(npoin2,ndire,nf)
102  DOUBLE PRECISION, INTENT(INOUT) :: CY(npoin2,ndire,nf)
103  DOUBLE PRECISION, INTENT(INOUT) :: CT(npoin2,ndire,nf)
104  DOUBLE PRECISION, INTENT(INOUT) :: CF(npoin2,ndire,nf)
105  DOUBLE PRECISION, INTENT(IN) :: CG(npoin2,nf),XK(npoin2,nf)
106 !
107 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
108 !
109  INTEGER IP,IPOIN
110  DOUBLE PRECISION GSQP,SRCF,TFSR,DDDN,LSDUDN,LSDUDS
111  DOUBLE PRECISION USGD,DEUKD,TR1,TR2
112 !
113 !***********************************************************************
114 !
115  gsqp=gravit/(2.d0*deupi)
116 !
117 !-----------------------------------------------------------------------
118 ! INFINITE WATER DEPTH ...
119 !-----------------------------------------------------------------------
120 !
121  IF(proinf) THEN
122 !
123 ! ----------------------------------------------------------------
124 ! ... AND IN CARTESIAN COORDINATE SYSTEM
125 ! ----------------------------------------------------------------
126 !
127  IF(.NOT.sphe) THEN
128 !
129  DO ip=1,ndire
130  tr1=gsqp/freq(jf)*costet(ip)
131  tr2=gsqp/freq(jf)*sintet(ip)
132  DO ipoin=1,npoin2
133  lsdudn= sintet(ip)*
134  & (-costet(ip)*dvy(ipoin)-sintet(ip)*duy(ipoin))
135  & + costet(ip)*
136  & ( costet(ip)*dvx(ipoin)+sintet(ip)*dux(ipoin))
137  lsduds= costet(ip)*
138  & (costet(ip)*dvy(ipoin)+sintet(ip)*duy(ipoin))
139  & + sintet(ip)*
140  & (costet(ip)*dvx(ipoin)+sintet(ip)*dux(ipoin))
141  cx(ipoin,ip,jf)=tr2+uc(ipoin)
142  cy(ipoin,ip,jf)=tr1+vc(ipoin)
143  ct(ipoin,ip,jf)=-lsdudn
144  cf(ipoin,ip,jf)=-cg(ipoin,jf)*xk(ipoin,jf)*lsduds*usdpi
145  ENDDO
146  ENDDO
147 !
148 ! ----------------------------------------------------------------
149 ! ... AND IN SPHERICAL COORDINATE SYSTEM
150 ! ----------------------------------------------------------------
151 !
152  ELSE
153 !
154  DO ip=1,ndire
155  tr1=gsqp/freq(jf)*costet(ip)
156  tr2=gsqp/freq(jf)*sintet(ip)
157  DO ipoin=1,npoin2
158  srcf=sr/cosf(ipoin)
159  lsdudn= sintet(ip)*sr*
160  & (-costet(ip)*dvy(ipoin)-sintet(ip)*duy(ipoin))
161  & + costet(ip)*srcf*
162  & ( costet(ip)*dvx(ipoin)+sintet(ip)*dux(ipoin))
163  lsduds= costet(ip)*sr*
164  & (costet(ip)*dvy(ipoin)+sintet(ip)*duy(ipoin))
165  & + sintet(ip)*srcf*
166  & (costet(ip)*dvx(ipoin)+sintet(ip)*dux(ipoin))
167  cx(ipoin,ip,jf)=(tr2+uc(ipoin))*gradeg*srcf
168  cy(ipoin,ip,jf)=(tr1+vc(ipoin))*gradeg*sr
169  ct(ipoin,ip,jf)=tr2*tgf(ipoin)*sr - lsdudn*gradeg
170  cf(ipoin,ip,jf)= - lsduds*gradeg*
171  & cg(ipoin,jf)*xk(ipoin,jf)*usdpi
172  ENDDO
173  ENDDO
174  ENDIF
175 !
176 !-----------------------------------------------------------------------
177 ! FINITE WATER DEPTH ....
178 !-----------------------------------------------------------------------
179 !
180  ELSE
181 !
182 ! ----------------------------------------------------------------
183 ! ... AND IN CARTESIAN COORDINATE SYSTEM
184 ! ----------------------------------------------------------------
185 !
186  IF(.NOT.sphe) THEN
187 !
188  DO ipoin=1,npoin2
189  deukd=2.d0*xk(ipoin,jf)*depth(ipoin)
190  IF(deukd.GT.7.d2) THEN
191  tra01(ipoin) = 0.d0
192  ELSE
193  tra01(ipoin) = deupi*freq(jf)/sinh(deukd)
194  ENDIF
195  ENDDO
196 !
197  DO ip=1,ndire
198  DO ipoin=1,npoin2
199  dddn=-sintet(ip)*dzy(ipoin)+costet(ip)*dzx(ipoin)
200  cy(ipoin,ip,jf)=cg(ipoin,jf)*costet(ip)
201  cx(ipoin,ip,jf)=cg(ipoin,jf)*sintet(ip)
202  ct(ipoin,ip,jf)=-tra01(ipoin)*dddn
203  ENDDO
204  ENDDO
205 !
206  DO ipoin=1,npoin2
207  deukd=2.d0*xk(ipoin,jf)*depth(ipoin)
208  IF(deukd.GT.7.d2) THEN
209  tra01(ipoin)=0.d0
210  ELSE
211  tra01(ipoin)=xk(ipoin,jf)*deupi*freq(jf)/sinh(deukd)
212  ENDIF
213  ENDDO
214 !
215  DO ip=1,ndire
216  DO ipoin=1,npoin2
217  lsdudn= sintet(ip)*
218  & (-costet(ip)*dvy(ipoin)-sintet(ip)*duy(ipoin))
219  & + costet(ip)*
220  & ( costet(ip)*dvx(ipoin)+sintet(ip)*dux(ipoin))
221  lsduds= costet(ip)*
222  & (costet(ip)*dvy(ipoin)+sintet(ip)*duy(ipoin))
223  & + sintet(ip)*
224  & (costet(ip)*dvx(ipoin)+sintet(ip)*dux(ipoin))
225  usgd=vc(ipoin)*dzy(ipoin)+uc(ipoin)*dzx(ipoin)
226  cx(ipoin,ip,jf)=cx(ipoin,ip,jf) + uc(ipoin)
227  cy(ipoin,ip,jf)=cy(ipoin,ip,jf) + vc(ipoin)
228  ct(ipoin,ip,jf)=ct(ipoin,ip,jf) - lsdudn
229  cf(ipoin,ip,jf)= (tra01(ipoin)*(usgd+dzhdt(ipoin))
230  & - lsduds*cg(ipoin,jf)*xk(ipoin,jf))*usdpi
231  ENDDO
232  ENDDO
233 !
234 ! --------------------------------------------------------------
235 ! ... AND IN SPHERICAL COORDINATE SYSTEM
236 ! --------------------------------------------------------------
237 !
238  ELSE
239 !
240  DO ipoin=1,npoin2
241  deukd=2.d0*xk(ipoin,jf)*depth(ipoin)
242  IF(deukd.GT.7.d2) THEN
243  tra01(ipoin) = 0.d0
244  ELSE
245  tra01(ipoin) = deupi*freq(jf)/sinh(deukd)
246  ENDIF
247  ENDDO
248 !
249  DO ip=1,ndire
250  DO ipoin=1,npoin2
251  srcf=sr/cosf(ipoin)
252  tfsr=tgf(ipoin)*sr
253  dddn=-sintet(ip)*dzy(ipoin)*sr+costet(ip)*dzx(ipoin)*srcf
254  cy(ipoin,ip,jf)=(cg(ipoin,jf)*costet(ip))*sr*gradeg
255  cx(ipoin,ip,jf)=(cg(ipoin,jf)*sintet(ip))*srcf*gradeg
256  ct(ipoin,ip,jf)=cg(ipoin,jf)*sintet(ip)*tfsr
257  & -tra01(ipoin)*dddn*gradeg
258  ENDDO
259  ENDDO
260 !
261  DO ipoin=1,npoin2
262  deukd=2.d0*xk(ipoin,jf)*depth(ipoin)
263  IF(deukd.GT.7.d2) THEN
264  tra01(ipoin)=0.d0
265  ELSE
266  tra01(ipoin)=xk(ipoin,jf)*deupi*freq(jf)/sinh(deukd)
267  ENDIF
268  ENDDO
269 !
270  DO ip=1,ndire
271  DO ipoin=1,npoin2
272  srcf=sr/cosf(ipoin)
273  lsdudn= sintet(ip)*sr*
274  & (-costet(ip)*dvy(ipoin)-sintet(ip)*duy(ipoin))
275  & + costet(ip)*srcf*
276  & ( costet(ip)*dvx(ipoin)+sintet(ip)*dux(ipoin))
277  lsduds= costet(ip)*sr*
278  & ( costet(ip)*dvy(ipoin)+sintet(ip)*duy(ipoin))
279  & + sintet(ip)*srcf*
280  & ( costet(ip)*dvx(ipoin)+sintet(ip)*dux(ipoin))
281  usgd=vc(ipoin)*dzy(ipoin)*sr
282  & +uc(ipoin)*dzx(ipoin)*srcf
283  cy(ipoin,ip,jf)=cy(ipoin,ip,jf)+vc(ipoin)*sr*gradeg
284  cx(ipoin,ip,jf)=cx(ipoin,ip,jf)+uc(ipoin)*srcf*gradeg
285  ct(ipoin,ip,jf)=ct(ipoin,ip,jf)-lsdudn*gradeg
286  cf(ipoin,ip,jf)=
287  & (tra01(ipoin)*(usgd*gradeg+dzhdt(ipoin))
288  & -lsduds*gradeg*cg(ipoin,jf)*xk(ipoin,jf))*usdpi
289  ENDDO
290  ENDDO
291 !
292  ENDIF
293 !
294  ENDIF
295 !
296 !-----------------------------------------------------------------------
297 !
298  RETURN
299  END
300 
subroutine conw4d(CX, CY, CT, CF, XK, CG, NPOIN2, NDIRE, JF, NF)
Definition: conw4d.f:7