The TELEMAC-MASCARET system  trunk
coupev.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE coupev
3 ! *****************
4 !
5  &(at,z,u,v,w,
6  & shp,imseg,x2dv,y2dv,distor,ikles,
7  & elem,nc2dv,npoin2,nelem2,ncou,fformat,im,jm,
8  & titcas,nva3,tab,textlu,ienre)
9 !
10 !***********************************************************************
11 ! POSTEL3D VERSION 6.2 01/09/99 T. DENOT (LNH) 01 30 87 74 89
12 ! FORTRAN90
13 !***********************************************************************
14 !
15 ! FONCTION : ECRIT POUR CHAQUE COUPE VERTICALES LES VARIABLES
16 ! D'UN PAS DE TEMPS
17 !
18 !-----------------------------------------------------------------------
19 ! ARGUMENTS
20 ! .________________.____.______________________________________________.
21 ! ! NOM !MODE! ROLE !
22 ! !________________!____!______________________________________________!
23 ! ! IREC ! -->! PAS DE TEMPS TRAITE !
24 ! ! AT ! -->! TEMPS CORRESPONDANT AU PAS TRAITE !
25 ! ! Z ! -->! COTES DES NOEUDS !
26 ! ! U,V,W ! -->! COMPOSANTES 3D DE LA VITESSE !
27 ! ! TA,TP ! -->! CONCENTRATIONS DES TRACEURS !
28 ! ! NUX,NUY,NUZ ! -->! COEFFICIENTS DE VISCOSITE POUR LES VITESSES !
29 ! ! NAX,NAY,NAZ ! -->! COEFFICIENTS DE VISCOSITE POUR LES TR.ACTIFS !
30 ! ! NPX,NPY,NPZ ! -->! COEFFICIENTS DE VISCOSITE POUR LES TR.PASSIFS!
31 ! ! RI ! -->! NOMBRE DE RICHARDSON !
32 ! ! AK,EP ! -->! VARIABLES DU MODELE K-EPSILON !
33 ! ! RHO ! -->! ECARTS RELATIFS DE DENSITE !
34 ! ! SHP ! -->! COORDONNEES BARYCENTRIQUES DES PTS DE COUPE !
35 ! ! TAB1,2,3 !<-- ! TABLEAU DE TRAVAIL POUR PROJETER LES VAR. !
36 ! ! NSEG ! -->! NOMBRE DE SEGMENTS CONSTITUANT CHAQUE COUPE !
37 ! ! IMSEG ! -->! NOMBRE DE POINTS PAR SEGMENTS !
38 ! ! X2DV ! -->! ABSCISSES DES SOMMETS DES COUPES VERTICALES !
39 ! ! Y2DV ! -->! ORDONNEES DES SOMMETS DES COUPES VERTICALES !
40 ! ! DISTOR ! -->! DISTORSION SUIVANT Z DE CHAQUE COUPE VERTICALE
41 ! ! IKLES ! -->! TABLE DE CONNECTIVITE !
42 ! ! ELEM ! -->! NUMERO DES ELEMENTS CONTENANT LES PTS DE COUPE
43 ! ! NC2DV ! -->! NOMBRE DE COUPES VERTICALES !
44 ! ! NPOIN2 ! -->! NOMBRE DE POINTS DU MAILLAGE 2D !
45 ! ! NELEM2 ! -->! NOMBRE D'ELEMENTS DU MAILLAGE 2D !
46 ! ! NCOU ! -->! NUMERO DE CANAL - 1 DE LA PREMIERE COUPE !
47 ! ! IM (LU) ! -->! NOMBRE DE PTS DE COUPE SUIVANT L'HORIZONTALE !
48 ! ! JM (=NPLAN) ! -->! NOMBRE DE PTS DE COUPE SUIVANT LA VERTICALE !
49 ! ! SORG3D ! -->! INDICATEUR DES VARIABLES ENREGISTREES !
50 ! ! TITCAS ! -->! TITRE A PORTER SUR CHAQUE COUPE !
51 ! !________________!____!______________________________________________!
52 ! MODE : -->(DONNEE NON MODIFIEE), <--(RESULTAT), <-->(DONNEE MODIFIEE)
53 !-----------------------------------------------------------------------
54 !
55 ! SOUS-PROGRAMME APPELE PAR : POSTEL3D
56 ! SOUS-PROGRAMME APPELES : ECRDEB , ECRI2
57 !
58 !
59 !history Y AUDOUIN (LNHE)
60 !+ 25/05/2015
61 !+ V7P0
62 !+ Modification to comply with the hermes module
63 
64 ! JUNE 2012 - P.LANG / INGEROP : SERAFIN OUTPUT FORMAT
65 !**********************************************************************
66 !
67  USE bief
69  USE declarations_postel3d, ONLY: prez => z
70  USE interface_postel3d, ex_coupev => coupev
71 !
73  IMPLICIT NONE
74 !
75 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
76 !
77  INTEGER,INTENT(IN) :: NPOIN2,NELEM2,IM,JM,NC2DV
78  INTEGER, INTENT(INOUT) :: NCOU(nc2dv)
79  DOUBLE PRECISION ,INTENT(IN) ::AT
80  DOUBLE PRECISION , INTENT(INOUT) :: SHP(im,3,nc2dv)
81  type(bief_obj), INTENT(INOUT) :: tab
82  INTEGER , INTENT(IN) :: IENRE
83  DOUBLE PRECISION, INTENT(IN) :: U(npoin2,jm),V(npoin2,jm)
84  DOUBLE PRECISION, INTENT(IN) :: Z(npoin2,jm),W(npoin2,jm)
85  INTEGER, INTENT(IN) :: IKLES(3,nelem2)
86  INTEGER, INTENT(IN) :: ELEM(im,nc2dv)
87  INTEGER, INTENT(INOUT) :: NVA3
88  INTEGER, INTENT(IN) :: IMSEG(49,nc2dv)
89  DOUBLE PRECISION, INTENT(IN) :: X2DV(50,nc2dv),Y2DV(50,nc2dv)
90  DOUBLE PRECISION, INTENT(IN) :: DISTOR(nc2dv)
91  CHARACTER(LEN=8), INTENT(INOUT) :: FFORMAT
92  CHARACTER(LEN=32), INTENT(IN) ::TEXTLU(100)
93  CHARACTER(LEN=72), INTENT(IN) ::TITCAS
94 !
95 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
96 !
97 !
98  DOUBLE PRECISION TAB1(im,jm),TAB2(im,jm),TAB3(im,jm)
99  DOUBLE PRECISION LGDEB,LGSEG,ALFA,COST,SINT,A1,A2,A3,U1,V1
100 !
101  INTEGER IB(10),IC,N1,N2,N3,I,J,K,CANAL
102  INTEGER ISEG,IDSEG,IFSEG
103 !
104 ! NEW VARIABLES FOR SERAFIN FORMAT
105 !
106  INTEGER IKLE(((im-1)*(jm-1))*2,3),IPOBO(im*jm),NUMELEM
107 !
108 ! END OF NEW VARIABLES
109 !
110  LOGICAL FLAG
111 !
112 !
113 !
114  CHARACTER(LEN=32) :: VAR_NAME
115  INTEGER :: IERR,IREC
116  INTEGER DATE(3),TIME(3)
117 !
118 !***********************************************************************
119 !
120 ! LISTE DE FUTURS PARAMETRES DEJA PREVUS.(SEUL LES PREMIERS SERVENT)
121 !
122  DO i=1,10
123  ib(i)=0
124  ENDDO
125 ! ECRITURE ECLATEE DES RESULTATS (CONVENTION LEONARD)
126  ib(2)=1
127 !
128 !-----------------------------------------------------------------------
129 !
130 ! POUR CHAQUE COUPE VERTICALE FAIRE :
131 !
132  irec = 0
133  DO ic = 1,nc2dv
134 !
135 ! OUVERTURE DU FICHIER + ENREGISTREMENT DES PREMIERS PARAMETRES
136 ! -------------------------------------------------------------
137 !
138  CALL ecrdeb(ncou(ic),fformat,titcas,nva3,.false.,
139  & textlu,ic,ienre)
140 !
141 ! CALCUL DES AUTRES PARAMETRES DE L'ENTETE
142 ! ----------------------------------------
143 !
144 ! MAILLAGE LEONARD ASSOCIE A LA COUPE IC ET AU PAS DE TEMPS IT
145 !
146  iseg = 0
147  ifseg = 1
148  lgdeb = 0.d0
149  lgseg = 0.d0
150 !
151  DO i = 1,im
152 !
153 ! COORDONNEE HORIZONTALE SUIVANT LE PLAN DE COUPE (X)
154 !
155  IF (i.GT.ifseg.OR.i.EQ.1) THEN
156  iseg = iseg + 1
157  idseg = ifseg
158  ifseg = ifseg + imseg(iseg,ic)
159  lgdeb = lgdeb + lgseg
160  lgseg = sqrt((x2dv(iseg+1,ic)-x2dv(iseg,ic))**2
161  & +(y2dv(iseg+1,ic)-y2dv(iseg,ic))**2)
162  ENDIF
163 !
164  tab1(i,1) = lgdeb + float(i-idseg)*lgseg/float(ifseg-idseg)
165 !
166 ! COORDONNEE VERTICALE (Y)
167 !
168  DO j = 1,jm
169 !
170  tab1(i,j) = tab1(i,1)
171  tab2(i,j) = ( shp(i,1,ic)*z(ikles(1,elem(i,ic)),j)
172  & + shp(i,2,ic)*z(ikles(2,elem(i,ic)),j)
173  & + shp(i,3,ic)*z(ikles(3,elem(i,ic)),j) )
174  & * distor(ic)
175 !
176  ENDDO
177  ENDDO !I
178 !
179 ! ENREGISTREMENT DES AUTRES PARAMETRES DE L'ENTETE
180 ! ------------------------------------------------
181 !
182 ! BUILD OF IKLE ARRAY. EACH QUADRANGLE IS DIVIDED INTO 2 TRIANGLES
183 ! NUMELEM VARIABLE MANAGES THE NUMBER OF ELEMENT
184  numelem = 1
185  DO j = 1,jm-1
186  DO i = 1,im-1
187  ikle(numelem,1) = ((j-1)*im)+i
188  ikle(numelem,2) = ((j-1)*im)+i+1
189  ikle(numelem,3) = ((j)*im)+i+1
190  numelem = numelem+1
191  ikle(numelem,1) = ((j-1)*im)+i
192  ikle(numelem,2) = ((j)*im)+i+1
193  ikle(numelem,3) = ((j)*im)+i
194  numelem = numelem+1
195  ENDDO
196  ENDDO
197  numelem = numelem-1
198 ! RECORD NELEM,NPOIN,NDP,1
199  ib(1) = ((im-1)*(jm-1)) * 2
200  ib(2) = im*jm
201  ib(3) = 3
202  ib(4) = 1
203 ! IKLE STORAGE
204 ! IPOBO ARRAY (WITH DUMMY VALUE)
205  DO i=1,ib(2)
206  ipobo(i) = 0
207  ENDDO
208 ! X AND Y COORDINATES
209  date = (/0,0,0/)
210  time = (/0,0,0/)
211  canal = ncou(ic)
212  CALL set_mesh(fformat,canal,2,triangle_elt_type,3,0,0,numelem,
213  & ib(2),ikle,ipobo,ipobo,tab1,tab2,0,
214  & date,time,0,0,ierr)
215  CALL check_call(ierr,'COUPEV:SET_MESH')
216 !
217 !-----------------------------------------------------------------------
218 !
219 ! SORTIE DES VARIABLES
220 !
221 ! 3 COMPOSANTES DE LA VITESSE
222 ! ---------------------------
223 !
224 !
225  iseg = 1
226  ifseg = 1 + imseg(1,ic)
227  alfa = atan2(y2dv(2,ic)-y2dv(1,ic),x2dv(2,ic)-x2dv(1,ic))
228  flag = .true.
229 !
230  DO i = 1,im
231 !
232  IF (flag) cost = cos(alfa)
233  IF (flag) sint = sin(alfa)
234  flag = .false.
235 !
236  IF (i.EQ.ifseg.AND.i.NE.im) THEN
237  flag = .true.
238  iseg = iseg + 1
239  ifseg = ifseg + imseg(iseg,ic)
240  a1 = alfa
241  alfa = atan2(y2dv(iseg+1,ic)-y2dv(iseg,ic),
242  & x2dv(iseg+1,ic)-x2dv(iseg,ic))
243  cost = cos(0.5d0*(alfa+a1))
244  sint = sin(0.5d0*(alfa+a1))
245  ENDIF
246 !
247  n1 = ikles(1,elem(i,ic))
248  n2 = ikles(2,elem(i,ic))
249  n3 = ikles(3,elem(i,ic))
250  a1 = shp(i,1,ic)
251  a2 = shp(i,2,ic)
252  a3 = shp(i,3,ic)
253 !
254  DO j = 1,jm
255 !
256  u1 = a1*u(n1,j) + a2*u(n2,j) + a3*u(n3,j)
257  v1 = a1*v(n1,j) + a2*v(n2,j) + a3*v(n3,j)
258 !
259 ! COMPOSANTE TANGENTIELLE ET HORIZONTALE DE LA VITESSE (UT)
260 !
261  tab1(i,j) = cost*u1 + sint*v1
262 !
263 ! COMPOSANTE VERTICALE DE LA VITESSE (W)
264 !
265  tab2(i,j) = (a1*w(n1,j)+a2*w(n2,j)+a3*w(n3,j))*distor(ic)
266 !
267 ! COMPOSANTE NORMALE ET HORIZONTALE DE LA VITESSE (UN)
268 !
269  tab3(i,j) = -sint*u1 + cost*v1
270 !
271  ENDDO
272  ENDDO !I
273 !
274  IF (lng.EQ.lng_fr) var_name =
275  & 'VITESSE UT M/S '
276  IF (lng.EQ.lng_en) var_name =
277  & 'VELOCITY UT M/S '
278  CALL add_data(fformat,canal,var_name,at,irec,.true.,tab1,
279  & im*jm,ierr)
280  CALL check_call(ierr,'COUPEV:ADD_DATA:UT')
281  IF (lng.EQ.lng_fr) var_name =
282  & 'VITESSE W M/S '
283  IF (lng.EQ.lng_en) var_name =
284  & 'VELOCITY W M/S '
285  CALL add_data(fformat,canal,var_name,at,irec,.false.,tab2,
286  & im*jm,ierr)
287  CALL check_call(ierr,'COUPEV:ADD_DATA:W')
288  IF (lng.EQ.lng_fr) var_name =
289  & 'VITESSE UN M/S '
290  IF (lng.EQ.lng_en) var_name =
291  & 'VELOCITY UN M/S '
292  CALL add_data(fformat,canal,var_name,at,irec,.false.,tab3,
293  & im*jm,ierr)
294  CALL check_call(ierr,'COUPEV:ADD_DATA:UN')
295 !
296 ! Adding z
297 !
298  DO j = 1,jm
299  DO i = 1,im
300  tab1(i,j) = shp(i,1,ic)
301  & *prez(ikles(1,elem(i,ic))+(j-1)*npoin2)
302  & + shp(i,2,ic)
303  & *prez(ikles(2,elem(i,ic))+(j-1)*npoin2)
304  & + shp(i,3,ic)
305  & *prez(ikles(3,elem(i,ic))+(j-1)*npoin2)
306  ENDDO !I
307  ENDDO !J
308  var_name = textlu(1)
309  CALL add_data(fformat,canal,var_name,at,irec,.false.,tab1,
310  & im*jm,ierr)
311  CALL check_call(ierr,'COUPEV:ADD_DATA:Z')
312 !
313 ! other variables
314 !
315  IF (nva3.GT.4) THEN
316  DO k = 5,nva3
317  DO j = 1,jm
318  DO i = 1,im
319  tab1(i,j) = shp(i,1,ic)
320  & *tab%ADR(k-4)%P%R(ikles(1,elem(i,ic))+(j-1)*npoin2)
321  & + shp(i,2,ic)
322  & *tab%ADR(k-4)%P%R(ikles(2,elem(i,ic))+(j-1)*npoin2)
323  & + shp(i,3,ic)
324  & *tab%ADR(k-4)%P%R(ikles(3,elem(i,ic))+(j-1)*npoin2)
325  ENDDO !I
326  ENDDO !J
327  var_name = textlu(k)
328  CALL add_data(fformat,canal,var_name,at,irec,.false.,tab1,
329  & im*jm,ierr)
330  CALL check_call(ierr,'COUPEV:ADD_DATA:K')
331  ENDDO !K
332 !
333  ENDIF
334  CALL close_mesh(fformat,canal,ierr)
335  CALL check_call(ierr,'COUPEV:CLOSE_MESH')
336 !
337  ENDDO !IC
338 !
339 !-----------------------------------------------------------------------
340 !
341  RETURN
342  END SUBROUTINE
subroutine add_data(FFORMAT, FILE_ID, VAR_NAME, TIME, RECORD, FIRST_VAR, VAR_VALUE, N, IERR)
Definition: add_data.f:8
subroutine close_mesh(FFORMAT, FILE_ID, IERR, MESH_NUMBER)
Definition: close_mesh.f:7
subroutine set_mesh(FFORMAT, FILE_ID, MESH_DIM, TYPELM, NDP, NPTFR, NPTIR, NELEM, NPOIN, IKLE, IPOBO, KNOLG, X, Y, NPLAN, DATE, TIME, X_ORIG, Y_ORIG, IERR, Z, IN_PLACE)
Definition: set_mesh.f:9
integer, parameter lng_en
subroutine ecrdeb(CANAL, FFORMAT, TITCAS, NBVAR, C2DH, TEXTLU, IC, N)
Definition: ecrdeb.f:7
integer, parameter triangle_elt_type
integer, parameter lng_fr
Y. AUDOUIN & J-M HERVOUET (EDF LAB, LNHE) 09/05/2014 V7P0 First version.
double precision, dimension(:), pointer z
subroutine coupev(AT, Z, U, V, W, SHP, IMSEG, X2DV, Y2DV, DISTOR, IKLES, ELEM, NC2DV, NPOIN2, NELEM2, NCOU, FFORMAT, IM, JM, TITCAS, NVA3, TAB, TEXTLU, IENRE)
Definition: coupev.f:10
Definition: bief.f:3