The TELEMAC-MASCARET system  trunk
ov.f
Go to the documentation of this file.
1 ! *************
2  SUBROUTINE ov
3 ! *************
4 !
5  & ( op , x , y , z , c , dim1 )
6 !
7 !***********************************************************************
8 ! BIEF V6P1 21/08/2010
9 !***********************************************************************
10 !
11 !brief OPERATIONS ON VECTORS.
12 !code
13 !+ OP IS A STRING OF 8 CHARACTERS, WHICH INDICATES THE OPERATION TO BE
14 !+ PERFORMED ON VECTORS X,Y AND Z AND CONSTANT C.
15 !+
16 !+ THE RESULT IS VECTOR X.
17 !+
18 !+ OP = 'X=C ' : SETS X TO C
19 !+ OP = 'X=Y ' : COPIES Y IN X
20 !+ OP = 'X=+Y ' : IDEM
21 !+ OP = 'X=-Y ' : COPIES -Y IN X
22 !+ OP = 'X=1/Y ' : COPIES INVERSE OF Y IN X
23 !+ OP = 'X=Y+Z ' : ADDS Y AND Z
24 !+ OP = 'X=Y-Z ' : SUBTRACTS Z TO Y
25 !+ OP = 'X=YZ ' : Y.Z
26 !+ OP = 'X=-YZ ' : -Y.Z
27 !+ OP = 'X=XY ' : X.Y
28 !+ OP = 'X=X+YZ ' : ADDS Y.Z TO X
29 !+ OP = 'X=X-YZ ' : SUBTRACTS Y.Z FROM X
30 !+ OP = 'X=CXY ' : C.X.Y
31 !+ OP = 'X=CYZ ' : C.Y.Z
32 !+ OP = 'X=CXYZ ' : C.X.Y.Z
33 !+ OP = 'X=X+CYZ ' : ADDS C.Y.Z TO X
34 !+ OP = 'X=Y/Z ' : DIVIDES Y BY Z
35 !+ OP = 'X=CY/Z ' : DIVIDES C.Y BY Z
36 !+ OP = 'X=CXY/Z ' : DIVIDES C.X.Y BY Z
37 !+ OP = 'X=X+CY/Z' : ADDS C.Y/Z TO X
38 !+ OP = 'X=X+Y ' : ADDS Y TO X
39 !+ OP = 'X=X-Y ' : SUBTRACTS Y FROM X
40 !+ OP = 'X=CX ' : MULTIPLIES X BY C
41 !+ OP = 'X=CY ' : MULTIPLIES Y BY C
42 !+ OP = 'X=Y+CZ ' : ADDS C.Z TO Y
43 !+ OP = 'X=X+CY ' : ADDS C.Y TO X
44 !+ OP = 'X=SQR(Y)' : SQUARE ROOT OF Y
45 !+ OP = 'X=ABS(Y)' : ABSOLUTE VALUE OF Y
46 !+ OP = 'X=N(Y,Z)' : NORM OF THE VECTOR WITH COMPONENTS Y AND Z
47 !+ OP = 'X=Y+C ' : ADDS C TO Y
48 !+ OP = 'X=X+C ' : ADDS C TO X
49 !+ OP = 'X=Y**C ' : Y TO THE POWER C
50 !+ OP = 'X=COS(Y)' : COSINE OF Y
51 !+ OP = 'X=SIN(Y)' : SINE OF Y
52 !+ OP = 'X=ATN(Y)' : ARC TANGENT OF Y
53 !+ OP = 'X=A(Y,Z)' : ATAN2(Y,Z)
54 !+ OP = 'X=+(Y,C)' : MAXIMUM OF Y AND C
55 !+ OP = 'X=-(Y,C)' : MINIMUM OF Y AND C
56 !+ OP = 'X=+(Y,Z)' : MAXIMUM OF Y AND Z
57 !+ OP = 'X=-(Y,Z)' : MINIMUM OF Y AND Z
58 !+ OP = 'X=YIFZ<C' : COPIES Y IN X IF Z < C , FOR EACH POINT
59 !+ OP = 'X=C(Y-Z)' : MULTIPLIES C BY (Y-Z)
60 !
61 !warning DIVIDE OPERATIONS INTERNALLY TAKE CARE OF DIVISIONS BY 0.
62 !+ SUCCESSFUL EXIT OF OV IS THEREFORE NOT A PROOF THAT Y OR Z
63 !+ NEVER ARE 0. FOR SUCH OPERATIONS, PREFERABLY USE OVD
64 !
65 !history J-M HERVOUET (LNHE) ; F LEPEINTRE (LNH)
66 !+ 17/03/2010
67 !+ V6P0
68 !+ ORIGINAL IDEA IN ULYSSE. THANK YOU D. LAURENCE
69 !
70 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
71 !+ 13/07/2010
72 !+ V6P0
73 !+ Translation of French comments within the FORTRAN sources into
74 !+ English comments
75 !
76 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
77 !+ 21/08/2010
78 !+ V6P0
79 !+ Creation of DOXYGEN tags for automated documentation and
80 !+ cross-referencing of the FORTRAN sources
81 !
82 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
83 !| C |-->| A GIVEN CONSTANT
84 !| DIM1 |-->| SIZE OF VECTORS
85 !| OP |-->| STRING INDICATING THE OPERATION TO BE DONE
86 !| X |<--| RESULTING VECTOR
87 !| Y |-->| TO BE USED IN THE OPERATION
88 !| Z |-->| TO BE USED IN THE OPERATION
89 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
90 !
92  IMPLICIT NONE
93 !
94 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
95 !
96  INTEGER, INTENT(IN) :: DIM1
97  DOUBLE PRECISION, OPTIONAL, INTENT(IN) :: Y(*),Z(*),C
98  DOUBLE PRECISION, INTENT(INOUT) :: X(dim1)
99  CHARACTER(LEN=8), INTENT(IN) :: OP
100 !
101 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
102 !
103  INTEGER I
104 !
105  INTRINSIC sqrt,abs,cos,sin,atan,max,min,atan2
106 !
107 !-----------------------------------------------------------------------
108 !
109  SELECT CASE(op(3:8))
110 !
111 !-----------------------------------------------------------------------
112 !
113  CASE('C ')
114 !
115  DO i=1,dim1
116  x(i) = c
117  ENDDO
118 !
119 !-----------------------------------------------------------------------
120 !
121  CASE('0 ')
122 !
123  DO i=1,dim1
124  x(i) = 0.d0
125  ENDDO
126 !
127 !-----------------------------------------------------------------------
128 !
129  CASE('Y ')
130 !
131  DO i=1,dim1
132  x(i) = y(i)
133  ENDDO
134 !
135 !-----------------------------------------------------------------------
136 !
137  CASE('+Y ')
138 !
139  DO i=1,dim1
140  x(i) = y(i)
141  ENDDO
142 !
143 !-----------------------------------------------------------------------
144 !
145  CASE('-Y ')
146 !
147  DO i=1,dim1
148  x(i) = - y(i)
149  ENDDO
150 !
151 !-----------------------------------------------------------------------
152 !
153  CASE('1/Y ')
154 !
155  DO i=1,dim1
156  x(i) = 1.d0/y(i)
157  ENDDO
158 !
159 !-----------------------------------------------------------------------
160 !
161  CASE('Y+Z ')
162 !
163  DO i=1,dim1
164  x(i) = y(i) + z(i)
165  ENDDO
166 !
167 !-----------------------------------------------------------------------
168 !
169  CASE('Y-Z ')
170 !
171  DO i=1,dim1
172  x(i) = y(i) - z(i)
173  ENDDO
174 !
175 !-----------------------------------------------------------------------
176 !
177  CASE('YZ ')
178 !
179  DO i=1,dim1
180  x(i) = y(i) * z(i)
181  ENDDO
182 !
183 !-----------------------------------------------------------------------
184 !
185  CASE('-YZ ')
186 !
187  DO i=1,dim1
188  x(i) = - y(i) * z(i)
189  ENDDO
190 !
191 !-----------------------------------------------------------------------
192 !
193  CASE('XY ')
194 !
195  DO i=1,dim1
196  x(i) = x(i) * y(i)
197  ENDDO
198 !
199 !-----------------------------------------------------------------------
200 !
201  CASE('X+YZ ')
202 !
203  DO i=1,dim1
204  x(i) = x(i) + y(i) * z(i)
205  ENDDO
206 !
207 !-----------------------------------------------------------------------
208 !
209  CASE('X-YZ ')
210 !
211  DO i=1,dim1
212  x(i) = x(i) - y(i) * z(i)
213  ENDDO
214 !
215 !-----------------------------------------------------------------------
216 !
217  CASE('CXY ')
218 !
219  DO i=1,dim1
220  x(i) = c * x(i) * y(i)
221  ENDDO
222 !
223 !-----------------------------------------------------------------------
224 !
225  CASE('CYZ ')
226 !
227  DO i=1,dim1
228  x(i) = c * y(i) * z(i)
229  ENDDO
230 !
231 !-----------------------------------------------------------------------
232 !
233  CASE('CXYZ ')
234 !
235  DO i=1,dim1
236  x(i) = c * x(i) * y(i) * z(i)
237  ENDDO
238 !
239 !-----------------------------------------------------------------------
240 !
241  CASE('X+CYZ ')
242 !
243  DO i=1,dim1
244  x(i) = x(i) + c * y(i) * z(i)
245  ENDDO
246 !
247 !-----------------------------------------------------------------------
248 !
249  CASE('Y/Z ')
250 !
251  DO i=1,dim1
252  x(i) = y(i) / z(i)
253  ENDDO
254 !
255 !-----------------------------------------------------------------------
256 !
257  CASE('CY/Z ')
258 !
259  DO i=1,dim1
260  x(i) = c*y(i) / z(i)
261  ENDDO
262 !
263 !-----------------------------------------------------------------------
264 !
265  CASE('CXY/Z ')
266 !
267  DO i=1,dim1
268  x(i) = c*x(i)*y(i) / z(i)
269  ENDDO
270 !
271 !-----------------------------------------------------------------------
272 !
273  CASE('X+CY/Z')
274 !
275  DO i=1,dim1
276  x(i) = x(i) + c * y(i) / z(i)
277  ENDDO
278 !
279 !-----------------------------------------------------------------------
280 !
281  CASE('X+Y ')
282 !
283  DO i=1,dim1
284  x(i) = x(i) + y(i)
285  ENDDO
286 !
287 !-----------------------------------------------------------------------
288 !
289  CASE('X-Y ')
290 !
291  DO i=1,dim1
292  x(i) = x(i) - y(i)
293  ENDDO
294 !
295 !-----------------------------------------------------------------------
296 !
297  CASE('CX ')
298 !
299  DO i=1,dim1
300  x(i) = c * x(i)
301  ENDDO
302 !
303 !-----------------------------------------------------------------------
304 !
305  CASE('CY ')
306 !
307  DO i=1,dim1
308  x(i) = c * y(i)
309  ENDDO
310 !
311 !-----------------------------------------------------------------------
312 !
313  CASE('Y+CZ ')
314 !
315  DO i=1,dim1
316  x(i) = y(i) + c * z(i)
317  ENDDO
318 !
319 !-----------------------------------------------------------------------
320 !
321  CASE('X+CY ')
322 !
323  DO i=1,dim1
324  x(i) = x(i) + c * y(i)
325  ENDDO
326 !
327 !-----------------------------------------------------------------------
328 !
329  CASE('SQR(Y)')
330 !
331  DO i=1,dim1
332  x(i) = sqrt(y(i))
333  ENDDO
334 !
335 !-----------------------------------------------------------------------
336 !
337  CASE('ABS(Y)')
338 !
339  DO i=1,dim1
340  x(i) = abs(y(i))
341  ENDDO
342 !
343 !-----------------------------------------------------------------------
344 !
345  CASE('N(Y,Z)')
346 !
347  DO i=1,dim1
348  x(i) = sqrt( y(i)**2 + z(i)**2 )
349  ENDDO
350 !
351 !-----------------------------------------------------------------------
352 !
353  CASE('Y+C ')
354 !
355  DO i=1,dim1
356  x(i) = y(i) + c
357  ENDDO
358 !
359 !-----------------------------------------------------------------------
360 !
361  CASE('X+C ')
362 !
363  DO i=1,dim1
364  x(i) = x(i) + c
365  ENDDO
366 !
367 !-----------------------------------------------------------------------
368 !
369  CASE('Y**C ')
370 !
371  DO i=1,dim1
372  IF(y(i).GE.0.d0) THEN
373  x(i) = y(i)**c
374  ELSE
375  WRITE(lu,101)
376 101 FORMAT(1x,'OV (BIEF): Y**C FORBIDDEN IF Y < 0')
377  CALL plante(1)
378  stop
379  ENDIF
380  ENDDO
381 !
382 !-----------------------------------------------------------------------
383 !
384  CASE('COS(Y)')
385 !
386  DO i=1,dim1
387  x(i) = cos(y(i))
388  ENDDO
389 !
390 !-----------------------------------------------------------------------
391 !
392  CASE('SIN(Y)')
393 !
394  DO i=1,dim1
395  x(i) = sin(y(i))
396  ENDDO
397 !
398 !-----------------------------------------------------------------------
399 !
400  CASE('ATN(Y)')
401 !
402  DO i=1,dim1
403  x(i) = atan(y(i))
404  ENDDO
405 !
406 !-----------------------------------------------------------------------
407 !
408  CASE('A(Y,Z)')
409 !
410 #ifdef NAGFOR
411  ! Bug in nag 6.2 it crashes id X=0 and Y=0 when is should be
412  ! PI/2
413  DO i=1,dim1
414  IF(y(i).EQ.0.d0.AND.z(i).EQ.0.d0) THEN
415  ! PI/2
416  x(i) = 2.d0*atan(1.d0)
417  ELSE
418  x(i) = atan2(y(i),z(i))
419  ENDIF
420  ENDDO
421 #else
422  DO i=1,dim1
423  x(i) = atan2(y(i),z(i))
424  ENDDO
425 #endif
426 !
427 !-----------------------------------------------------------------------
428 !
429  CASE('+(Y,C)')
430 !
431  DO i=1,dim1
432  x(i) = max(y(i),c)
433  ENDDO
434 !
435 !-----------------------------------------------------------------------
436 !
437  CASE('-(Y,C)')
438 !
439  DO i=1,dim1
440  x(i) = min(y(i),c)
441  ENDDO
442 !
443 !-----------------------------------------------------------------------
444 !
445  CASE('+(Y,Z)')
446 !
447  DO i=1,dim1
448  x(i) = max(y(i),z(i))
449  ENDDO
450 !
451 !-----------------------------------------------------------------------
452 !
453  CASE('-(Y,Z)')
454 !
455  DO i=1,dim1
456  x(i) = min(y(i),z(i))
457  ENDDO
458 !
459 !-----------------------------------------------------------------------
460 !
461  CASE('YIFZ<C')
462 !
463  DO i=1,dim1
464  IF ( z(i).LT.c ) x(i) = y(i)
465  ENDDO
466 !
467 !-----------------------------------------------------------------------
468 !
469  CASE('C(Y-Z)')
470 !
471  DO i=1,dim1
472  x(i) = c*(y(i)-z(i))
473  ENDDO
474 !
475 !-----------------------------------------------------------------------
476 !
477  CASE DEFAULT
478 !
479  WRITE(lu,1001) op
480 1001 FORMAT(1x,'OV (BIEF) : UNKNOWN OPERATION: ',a8)
481  CALL plante(1)
482  stop
483 !
484 !-----------------------------------------------------------------------
485 !
486  END SELECT
487 !
488 !-----------------------------------------------------------------------
489 !
490  RETURN
491  END
subroutine ov(OP, X, Y, Z, C, DIM1)
Definition: ov.f:7