Context Navigation

source: codes/icosagcm/trunk/src/metric.f90 @ 13

Last change on this file since 13 was 13, checked in by ymipsl, 12 years ago

some update

YM

File size: 25.2 KB

Line
1	MODULE metric
2	USE genmod
3	USE grid_param
4	USE domain_param
5
6	TYPE t_cell_glo
7	INTEGER :: neighbour(0:5)
8	INTEGER :: edge(0:5)
9	INTEGER :: assign_face
10	INTEGER :: assign_i
11	INTEGER :: assign_j
12	INTEGER :: assign_domain
13	REAL(rstd) :: xyz(3)
14	END TYPE t_cell_glo
15
16	TYPE t_vertex_glo
17	REAL(rstd) :: xyz(3)
18	INTEGER :: neighbour(0:5)
19	INTEGER :: ne(0:5)
20	INTEGER :: ind
21	INTEGER :: delta
22	END TYPE t_vertex_glo
23
24	TYPE t_edge_glo
25	INTEGER :: e1
26	INTEGER :: e2
27	END TYPE t_edge_glo
28
29
30	TYPE(t_vertex_glo),ALLOCATABLE,SAVE :: vertex_glo(:,:,:)
31	TYPE(t_cell_glo),ALLOCATABLE,SAVE :: cell_glo(:)
32	TYPE(t_edge_glo),ALLOCATABLE,SAVE :: edge_glo(:)
33	INTEGER :: ncell_glo
34
35	INTEGER,ALLOCATABLE,SAVE :: Tab_index(:,:,:)
36	INTEGER,PARAMETER :: right=1
37	INTEGER,PARAMETER :: rup=2
38	INTEGER,PARAMETER :: lup=3
39	INTEGER,PARAMETER :: left=4
40	INTEGER,PARAMETER :: ldown=5
41	INTEGER,PARAMETER :: rdown=6
42
43	INTEGER,PARAMETER :: vrup=1
44	INTEGER,PARAMETER :: vup=2
45	INTEGER,PARAMETER :: vlup=3
46	INTEGER,PARAMETER :: vldown=4
47	INTEGER,PARAMETER :: vdown=5
48	INTEGER,PARAMETER :: vrdown=6
49
50
51	CONTAINS
52
53	SUBROUTINE allocate_metric
54	IMPLICIT NONE
55	INTEGER :: ind
56
57	ALLOCATE(vertex_glo(1-halo:iim_glo+halo,1-halo:jjm_glo+halo,nb_face))
58	ncell_glo=10(iim_glo-2)(jjm_glo-2)+(iim_glo-2)*20+12
59	ALLOCATE(cell_glo(ncell_glo))
60	ALLOCATE(tab_index(nb_face,nb_face,0:5))
61	ALLOCATE(edge_glo(ncell_glo*3))
62
63	DO ind=1,ncell_glo
64	cell_glo(ind)%neighbour(:)=0
65	ENDDO
66
67	END SUBROUTINE allocate_metric
68
69	SUBROUTINE Compute_face
70	USE spherical_geom_mod
71	IMPLICIT NONE
72
73	REAL(rstd) :: len_edge
74	REAL(rstd) :: rot=0.
75	INTEGER :: nf,i,j
76	REAL(rstd),DIMENSION(3) :: p1,p2,p3
77	REAL(rstd) :: d1,d2,d3
78
79
80	len_edge=acos(cos(Pi/5)cos(2Pi/5)/(sin(Pi/5)sin(2Pi/5)))
81
82	DO nf=1,nb_face/2
83	CALL lonlat2xyz(0.,Pi/2,vertex_glo(1,1,nf)%xyz)
84	CALL lonlat2xyz(rot+(nf-1)2Pi/5,Pi/2-len_edge,vertex_glo(iim_glo,1,nf)%xyz)
85	CALL lonlat2xyz(rot+2Pi/5+(nf-1)2*Pi/5,Pi/2-len_edge,vertex_glo(1,jjm_glo,nf)%xyz)
86	CALL lonlat2xyz(rot+Pi/5+(nf-1)2Pi/5,-Pi/2+len_edge,vertex_glo(iim_glo,jjm_glo,nf)%xyz)
87
88	CALL lonlat2xyz(0.,-Pi/2,vertex_glo(1,1,nf+nb_face/2)%xyz)
89	CALL lonlat2xyz(rot+Pi/5+(nf-1)2Pi/5,-(Pi/2-len_edge),vertex_glo(1,jjm_glo,nf+nb_face/2)%xyz)
90	CALL lonlat2xyz(rot+Pi/5+2Pi/5+(nf-1)2*Pi/5,-(Pi/2-len_edge),vertex_glo(iim_glo,1,nf+nb_face/2)%xyz)
91	CALL lonlat2xyz(rot+Pi/5+Pi/5+(nf-1)2Pi/5,-(-Pi/2+len_edge),vertex_glo(iim_glo,jjm_glo,nf+nb_face/2)%xyz)
92	ENDDO
93
94	DO nf=1,nb_face
95	DO i=2,iim_glo-1
96	CALL div_arc(vertex_glo(1,1,nf)%xyz,vertex_glo(iim_glo,1,nf)%xyz,(i-1)*1./(iim_glo-1),vertex_glo(i,1,nf)%xyz)
97	CALL div_arc(vertex_glo(1,jjm_glo,nf)%xyz,vertex_glo(iim_glo,jjm_glo,nf)%xyz, &
98	(i-1)*1./(iim_glo-1),vertex_glo(i,jjm_glo,nf)%xyz)
99	ENDDO
100	DO j=2,jjm_glo-1
101	CALL div_arc(vertex_glo(1,1,nf)%xyz,vertex_glo(1,jjm_glo,nf)%xyz,(j-1)*1./(jjm_glo-1),vertex_glo(1,j,nf)%xyz)
102	CALL div_arc(vertex_glo(iim_glo,1,nf)%xyz,vertex_glo(iim_glo,jjm_glo,nf)%xyz,(j-1)*1./(jjm_glo-1), &
103	vertex_glo(iim_glo,j,nf)%xyz)
104	ENDDO
105
106	DO j=2,jjm_glo-1
107	DO i=2,iim_glo-1
108	IF (i+j-1 > iim_glo) THEN
109	CALL div_arc(vertex_glo(iim_glo,i+j-iim_glo,nf)%xyz,vertex_glo(i+j-jjm_glo,jjm_glo,nf)%xyz, &
110	(iim_glo-i)1./(2iim_glo-(i+j)),vertex_glo(i,j,nf)%xyz)
111	ELSE
112	CALL div_arc(vertex_glo(i+j-1,1,nf)%xyz,vertex_glo(1,i+j-1,nf)%xyz,(j-1)*1./(i+j-2),vertex_glo(i,j,nf)%xyz)
113	ENDIF
114	ENDDO
115	ENDDO
116	ENDDO
117
118	CALL dist_cart(vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,d1)
119	CALL dist_cart(vertex_glo(2,1,1)%xyz,vertex_glo(3,1,1)%xyz,d2)
120	CALL dist_cart(vertex_glo(3,1,1)%xyz,vertex_glo(4,1,1)%xyz,d3)
121	CALL div_arc(vertex_glo(1,1,1)%xyz,vertex_glo(3,1,1)%xyz,0.5,p1)
122	! CALL div_arc(vertex_glo(2,1,1)%xyz,vertex_glo(1,3,1)%xyz,1./3,p1)
123	! CALL div_arc(vertex_glo(1,2,1)%xyz,vertex_glo(3,1,1)%xyz,1./3,p2)
124	! CALL div_arc(vertex_glo(2,2,1)%xyz,vertex_glo(1,1,1)%xyz,1./3,p3)
125	! PRINT *, "dist",d1
126	! PRINT *, "dist",d2
127	! PRINT *, "dist",d3
128	! PRINT *,"dist",vertex_glo(2,1,1)%xyz
129	! PRINT ,"dist",p1/sqrt(sum(p1*2))
130	CALL circonscrit(vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,vertex_glo(1,2,1)%xyz,p1)
131	! CALL Centroide(vertex_glo(1,2,1)%xyz,vertex_glo(2,1,1)%xyz,vertex_glo(1,1,1)%xyz,p1)
132	! CALL Centroide(vertex_glo(1,2,1)%xyz,vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,p1)
133	CALL dist_cart(vertex_glo(1,1,1)%xyz,p1,d1)
134	CALL dist_cart(vertex_glo(2,1,1)%xyz,p1,d2)
135	CALL dist_cart(vertex_glo(1,2,1)%xyz,p1,d3)
136	! PRINT *, "dist",d1
137	! PRINT *, "dist",d2
138	! PRINT *, "dist",d3
139
140	END SUBROUTINE compute_face
141
142	SUBROUTINE Compute_face_projection
143	USE spherical_geom_mod
144	IMPLICIT NONE
145
146	REAL(rstd) :: len_edge
147	REAL(rstd) :: rot=0.
148	INTEGER :: nf,i,j
149	REAL(rstd),DIMENSION(3) :: p1,p2,p3
150	REAL(rstd) :: d1,d2,d3
151
152	len_edge=acos(cos(Pi/5)cos(2Pi/5)/(sin(Pi/5)sin(2Pi/5)))
153
154	DO nf=1,nb_face/2
155	CALL lonlat2xyz(0.,Pi/2,vertex_glo(1,1,nf)%xyz)
156	CALL lonlat2xyz(rot+(nf-1)2Pi/5,Pi/2-len_edge,vertex_glo(iim_glo,1,nf)%xyz)
157	CALL lonlat2xyz(rot+2Pi/5+(nf-1)2*Pi/5,Pi/2-len_edge,vertex_glo(1,jjm_glo,nf)%xyz)
158	CALL lonlat2xyz(rot+Pi/5+(nf-1)2Pi/5,-Pi/2+len_edge,vertex_glo(iim_glo,jjm_glo,nf)%xyz)
159
160	CALL lonlat2xyz(0.,-Pi/2,vertex_glo(1,1,nf+nb_face/2)%xyz)
161	CALL lonlat2xyz(rot+Pi/5+(nf-1)2Pi/5,-(Pi/2-len_edge),vertex_glo(1,jjm_glo,nf+nb_face/2)%xyz)
162	CALL lonlat2xyz(rot+Pi/5+2Pi/5+(nf-1)2*Pi/5,-(Pi/2-len_edge),vertex_glo(iim_glo,1,nf+nb_face/2)%xyz)
163	CALL lonlat2xyz(rot+Pi/5+Pi/5+(nf-1)2Pi/5,-(-Pi/2+len_edge),vertex_glo(iim_glo,jjm_glo,nf+nb_face/2)%xyz)
164	ENDDO
165
166	DO nf=1,nb_face
167	DO i=2,iim_glo-1
168	CALL div_arc_bis(vertex_glo(1,1,nf)%xyz,vertex_glo(iim_glo,1,nf)%xyz,(i-1)*1./(iim_glo-1),vertex_glo(i,1,nf)%xyz)
169	CALL div_arc_bis(vertex_glo(1,jjm_glo,nf)%xyz,vertex_glo(iim_glo,jjm_glo,nf)%xyz, &
170	(i-1)*1./(iim_glo-1),vertex_glo(i,jjm_glo,nf)%xyz)
171	ENDDO
172	DO j=2,jjm_glo-1
173	CALL div_arc_bis(vertex_glo(1,1,nf)%xyz,vertex_glo(1,jjm_glo,nf)%xyz,(j-1)*1./(jjm_glo-1),vertex_glo(1,j,nf)%xyz)
174	CALL div_arc_bis(vertex_glo(iim_glo,1,nf)%xyz,vertex_glo(iim_glo,jjm_glo,nf)%xyz,(j-1)*1./(jjm_glo-1), &
175	vertex_glo(iim_glo,j,nf)%xyz)
176	ENDDO
177
178	DO j=2,jjm_glo-1
179	DO i=2,iim_glo-1
180	IF (i+j-1 > iim_glo) THEN
181	CALL div_arc_bis(vertex_glo(iim_glo,i+j-iim_glo,nf)%xyz,vertex_glo(i+j-jjm_glo,jjm_glo,nf)%xyz, &
182	(iim_glo-i)1./(2iim_glo-(i+j)),vertex_glo(i,j,nf)%xyz)
183	ELSE
184	CALL div_arc_bis(vertex_glo(i+j-1,1,nf)%xyz,vertex_glo(1,i+j-1,nf)%xyz,(j-1)*1./(i+j-2),vertex_glo(i,j,nf)%xyz)
185	ENDIF
186	ENDDO
187	ENDDO
188
189	! DO j=1,jjm_glo
190	! DO i=1,iim_glo
191	! vertex_glo(i,j,nf)%xyz=vertex_glo(i,j,nf)%xyz/sqrt(sum(vertex_glo(i,j,nf)%xyz**2))
192	! ENDDO
193	! ENDDO
194	ENDDO
195
196	! CALL dist_cart(vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,d1)
197	! CALL dist_cart(vertex_glo(2,1,1)%xyz,vertex_glo(3,1,1)%xyz,d2)
198	! CALL dist_cart(vertex_glo(3,1,1)%xyz,vertex_glo(4,1,1)%xyz,d3)
199	! CALL div_arc(vertex_glo(1,1,1)%xyz,vertex_glo(3,1,1)%xyz,0.5,p1)
200	! CALL div_arc(vertex_glo(2,1,1)%xyz,vertex_glo(1,3,1)%xyz,1./3,p1)
201	! CALL div_arc(vertex_glo(1,2,1)%xyz,vertex_glo(3,1,1)%xyz,1./3,p2)
202	! CALL div_arc(vertex_glo(2,2,1)%xyz,vertex_glo(1,1,1)%xyz,1./3,p3)
203	! PRINT *, "dist",d1
204	! PRINT *, "dist",d2
205	! PRINT *, "dist",d3
206	! PRINT *,"dist",vertex_glo(2,1,1)%xyz
207	! PRINT ,"dist",p1/sqrt(sum(p1*2))
208	! CALL Centroide(vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,vertex_glo(1,2,1)%xyz,p1)
209	! CALL Centroide(vertex_glo(1,2,1)%xyz,vertex_glo(2,1,1)%xyz,vertex_glo(1,1,1)%xyz,p1)
210	! CALL Centroide(vertex_glo(1,2,1)%xyz,vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,p1)
211	! CALL dist_cart(vertex_glo(1,1,1)%xyz,p1,d1)
212	! CALL dist_cart(vertex_glo(2,1,1)%xyz,p1,d2)
213	! CALL dist_cart(vertex_glo(1,2,1)%xyz,p1,d3)
214	! PRINT *, "dist",d1
215	! PRINT *, "dist",d2
216	! PRINT *, "dist",d3
217
218	END SUBROUTINE compute_face_projection
219
220
221	SUBROUTINE compute_extended_face_bis
222	IMPLICIT NONE
223
224	INTEGER :: nf
225	INTEGER :: i,j
226	INTEGER :: delta
227	INTEGER :: ind
228	INTEGER :: k,ng
229
230	DO nf=1,nb_face
231	DO i=2,iim_glo-1
232
233	CALL set_neighbour(i,1,nf,ldown-1)
234	CALL set_neighbour(i,1,nf,rdown-1)
235
236	CALL set_neighbour(i,jjm_glo,nf,lup-1)
237	CALL set_neighbour(i,jjm_glo,nf,rup-1)
238
239	ENDDO
240
241	DO j=2,jjm_glo-1
242
243	CALL set_neighbour(1,j,nf,left-1)
244	CALL set_neighbour(1,j,nf,lup-1)
245
246	CALL set_neighbour(iim_glo,j,nf,rdown-1)
247	CALL set_neighbour(iim_glo,j,nf,right-1)
248
249	ENDDO
250
251	CALL set_neighbour(2,0,nf,left-1)
252	CALL set_neighbour(iim_glo,0,nf,right-1)
253	CALL set_neighbour(iim_glo+1,jjm_glo-1,nf,rup-1)
254	CALL set_neighbour(iim_glo-1,jjm_glo+1,nf,right-1)
255	CALL set_neighbour(1,jjm_glo+1,nf,left-1)
256	CALL set_neighbour(0,2,nf,ldown-1)
257
258	CALL set_neighbour(1,0,nf,left-1)
259	CALL set_neighbour(iim_glo,0,nf,right-1)
260	CALL set_neighbour(0,jjm_glo,nf,rup-1)
261	CALL set_neighbour(iim_glo+1,jjm_glo,nf,rup-1)
262	ENDDO
263
264	DO nf=1,nb_face
265	DO j=0,jjm_glo+1
266	DO i=0,iim_glo+1
267	ind=vertex_glo(i,j,nf)%ind
268	delta=vertex_glo(i,j,nf)%delta
269	DO k=0,5
270	ng=MOD(delta+k+6,6)
271	vertex_glo(i,j,nf)%neighbour(k)=cell_glo(ind)%neighbour(ng)
272	ENDDO
273	ENDDO
274	ENDDO
275
276	i=1 ; j=1
277	vertex_glo(i,j,nf)%neighbour(0)=vertex_glo(i+1,j,nf)%ind
278	vertex_glo(i,j,nf)%neighbour(1)=vertex_glo(i,j+1,nf)%ind
279	vertex_glo(i,j,nf)%neighbour(2)=vertex_glo(i-1,j+1,nf)%ind
280	vertex_glo(i,j,nf)%neighbour(3)=vertex_glo(i-1,j,nf)%ind
281	vertex_glo(i,j,nf)%neighbour(4)=vertex_glo(i,j-1,nf)%ind
282	vertex_glo(i,j,nf)%neighbour(5)=vertex_glo(i+1,j-1,nf)%ind
283
284	i=iim_glo ; j=1
285	vertex_glo(i,j,nf)%neighbour(0)=vertex_glo(i+1,j,nf)%ind
286	vertex_glo(i,j,nf)%neighbour(1)=vertex_glo(i,j+1,nf)%ind
287	vertex_glo(i,j,nf)%neighbour(2)=vertex_glo(i-1,j+1,nf)%ind
288	vertex_glo(i,j,nf)%neighbour(3)=vertex_glo(i-1,j,nf)%ind
289	vertex_glo(i,j,nf)%neighbour(4)=vertex_glo(i,j-1,nf)%ind
290	vertex_glo(i,j,nf)%neighbour(5)=vertex_glo(i+1,j-1,nf)%ind
291
292	i=1 ; j=jjm_glo
293	vertex_glo(i,j,nf)%neighbour(0)=vertex_glo(i+1,j,nf)%ind
294	vertex_glo(i,j,nf)%neighbour(1)=vertex_glo(i,j+1,nf)%ind
295	vertex_glo(i,j,nf)%neighbour(2)=vertex_glo(i-1,j+1,nf)%ind
296	vertex_glo(i,j,nf)%neighbour(3)=vertex_glo(i-1,j,nf)%ind
297	vertex_glo(i,j,nf)%neighbour(4)=vertex_glo(i,j-1,nf)%ind
298	vertex_glo(i,j,nf)%neighbour(5)=vertex_glo(i+1,j-1,nf)%ind
299
300	i=iim_glo ; j=jjm_glo
301	vertex_glo(i,j,nf)%neighbour(0)=vertex_glo(i+1,j,nf)%ind
302	vertex_glo(i,j,nf)%neighbour(1)=vertex_glo(i,j+1,nf)%ind
303	vertex_glo(i,j,nf)%neighbour(2)=vertex_glo(i-1,j+1,nf)%ind
304	vertex_glo(i,j,nf)%neighbour(3)=vertex_glo(i-1,j,nf)%ind
305	vertex_glo(i,j,nf)%neighbour(4)=vertex_glo(i,j-1,nf)%ind
306	vertex_glo(i,j,nf)%neighbour(5)=vertex_glo(i+1,j-1,nf)%ind
307	ENDDO
308
309	CONTAINS
310
311	SUBROUTINE set_neighbour(i,j,nf,neighbour)
312	IMPLICIT NONE
313	INTEGER :: i,j,nf,neighbour
314	INTEGER :: in,jn
315	INTEGER :: k,ng
316	INTEGER :: delta
317	INTEGER :: ind,ind2
318
319	SELECT CASE (neighbour)
320	CASE(right-1)
321	in=i+1 ; jn=j
322	CASE(rup-1)
323	in=i ; jn=j+1
324	CASE(lup-1)
325	in=i-1 ; jn=j+1
326	CASE(left-1)
327	in=i-1 ; jn=j
328	CASE(ldown-1)
329	in=i ; jn=j-1
330	CASE(rdown-1)
331	in=i+1; jn=j-1
332	END SELECT
333
334	ind=vertex_glo(i,j,nf)%ind
335	IF (ind==0) THEN
336	PRINT *,"ind=0",i,j,nf
337	ENDIF
338	delta=MOD(vertex_glo(i,j,nf)%delta+neighbour+6,6)
339	ind2=cell_glo(ind)%neighbour(delta)
340	DO k=0,5
341	IF (cell_glo(ind2)%neighbour(k)==ind) ng=k
342	ENDDO
343	vertex_glo(in,jn,nf)%ind=ind2
344	vertex_glo(in,jn,nf)%delta=MOD(ng-neighbour+3+6,6)
345	END SUBROUTINE set_neighbour
346
347	END SUBROUTINE compute_extended_face_bis
348
349
350	SUBROUTINE compute_extended_face
351	IMPLICIT NONE
352
353	INTEGER :: nf,nf2
354	INTEGER :: ind,ind2
355	INTEGER :: i,j,h
356
357
358	DO h=1,halo
359	DO nf=1,nb_face
360	DO i=1-h+1,iim_glo+h-1
361	j=1-h+1
362	ind=vertex_glo(i,j,nf)%ind
363	nf2=cell_glo(ind)%assign_face
364	ind2=cell_glo(ind)%neighbour(tab_index(nf,nf2,ldown-1))
365	vertex_glo(i,j-1,nf)%ind=ind2
366	vertex_glo(i,j-1,nf)%xyz=cell_glo(ind2)%xyz
367	ENDDO
368
369	DO j=1-h,jjm_glo+h-1
370	i=iim_glo+h-1
371	ind=vertex_glo(i,j,nf)%ind
372	nf2=cell_glo(ind)%assign_face
373	ind2=cell_glo(ind)%neighbour(tab_index(nf,nf2,right-1))
374	vertex_glo(i+1,j,nf)%ind=ind2
375	vertex_glo(i+1,j,nf)%xyz=cell_glo(ind2)%xyz
376	ENDDO
377
378	DO i=iim_glo+h,1-h+1,-1
379	j=jjm_glo+h-1
380	ind=vertex_glo(i,j,nf)%ind
381	nf2=cell_glo(ind)%assign_face
382	ind2=cell_glo(ind)%neighbour(tab_index(nf,nf2,rup-1))
383	vertex_glo(i,j+1,nf)%ind=ind2
384	vertex_glo(i,j+1,nf)%xyz=cell_glo(ind2)%xyz
385	ENDDO
386
387	DO j=jjm_glo+h,1-h,-1
388	i=1-h+1
389	ind=vertex_glo(i,j,nf)%ind
390	nf2=cell_glo(ind)%assign_face
391	ind2=cell_glo(ind)%neighbour(tab_index(nf,nf2,left-1))
392	vertex_glo(i-1,j,nf)%ind=ind2
393	vertex_glo(i-1,j,nf)%xyz=cell_glo(ind2)%xyz
394	ENDDO
395	ENDDO
396	ENDDO
397
398
399	END SUBROUTINE compute_extended_face
400
401	SUBROUTINE Set_index
402	IMPLICIT NONE
403	INTEGER :: delta
404	INTEGER :: n1,n2,i
405
406	DO n1=1,5
407	DO delta=-2,2
408	DO i=0,5
409	n2=MOD(n1-1-delta+5,5)+1
410	tab_index(n1,n2,i)=MOD(delta+i+6,6)
411	n2=MOD(n1-1+delta+5,5)+1
412	tab_index(n1+5,n2+5,i)=MOD(delta+i+6,6)
413	ENDDO
414	ENDDO
415	ENDDO
416
417	DO n1=1,5
418	DO i=0,5
419	delta=3
420
421	n2=5+n1
422	Tab_index(n1,n2,i)=MOD(delta+i+6,6)
423	Tab_index(n2,n1,i)=MOD(delta+i+6,6)
424
425	n2=5+n1-1
426	IF (n2==5) n2=10
427	Tab_index(n1,n2,i)=MOD(delta+i+6,6)
428	Tab_index(n2,n1,i)=MOD(delta+i+6,6)
429	ENDDO
430	ENDDO
431
432	END SUBROUTINE set_index
433
434
435	SUBROUTINE set_cell
436	IMPLICIT NONE
437	INTEGER :: ind,ind2
438	INTEGER :: nf,nf2,nfm1,nfp1
439	INTEGER :: i,j,k
440	INTEGER :: delta
441
442	ind=0
443
444	!!!!!!!!!!!!!!!!!!!!!!!!!!!!
445	! Association des cellules !
446	!!!!!!!!!!!!!!!!!!!!!!!!!!!!
447
448	! interieur des faces
449	DO nf=1,nb_face
450	DO i=2,iim_glo-1
451	DO j=2,jjm_glo-1
452	ind=ind+1
453	vertex_glo(i,j,nf)%ind=ind
454	cell_glo(ind)%assign_face=nf
455	cell_glo(ind)%assign_i=i
456	cell_glo(ind)%assign_j=j
457	ENDDO
458	ENDDO
459	ENDDO
460
461	! frontieres faces nord
462	DO nf=1,nb_face/2
463
464	nfm1=nf-1
465	IF (nfm1==0) nfm1=nb_face/2
466
467	nf2=nf+nb_face/2
468
469	DO i=2,iim_glo-1
470	ind=ind+1
471	vertex_glo(i,1,nf)%ind=ind
472	vertex_glo(1,i,nfm1)%ind=ind
473	cell_glo(ind)%assign_face=nf
474	cell_glo(ind)%assign_i=i
475	cell_glo(ind)%assign_j=1
476	ENDDO
477
478	DO i=2,iim_glo-1
479	ind=ind+1
480	vertex_glo(i,jjm_glo,nf)%ind=ind
481	vertex_glo(iim_glo-i+1,jjm_glo,nf2)%ind=ind
482
483	cell_glo(ind)%assign_face=nf
484	cell_glo(ind)%assign_i=i
485	cell_glo(ind)%assign_j=jjm_glo
486	ENDDO
487	ENDDO
488
489	! frontieres faces sud
490
491	DO nf=nb_face/2+1,nb_face
492
493	nfp1=nf+1
494	IF (nfp1==nb_face+1) nfp1=nb_face/2+1
495
496	nf2=nf-nb_face/2+1
497	IF (nf2==nb_face/2+1) nf2=1
498
499	DO i=2,iim_glo-1
500	ind=ind+1
501	vertex_glo(i,1,nf)%ind=ind
502	vertex_glo(1,i,nfp1)%ind=ind
503	cell_glo(ind)%assign_face=nf
504	cell_glo(ind)%assign_i=i
505	cell_glo(ind)%assign_j=1
506	ENDDO
507
508	DO j=2,jjm_glo-1
509	ind=ind+1
510	vertex_glo(iim_glo,j,nf)%ind=ind
511	vertex_glo(iim_glo,jjm_glo-j+1,nf2)%ind=ind
512
513	cell_glo(ind)%assign_face=nf
514	cell_glo(ind)%assign_i=iim_glo
515	cell_glo(ind)%assign_j=j
516	ENDDO
517	ENDDO
518
519	! vertex nord (sur 3 faces)
520
521	DO nf=1,nb_face/2
522
523	nfp1=nf+1
524	IF (nfp1==nb_face/2+1) nfp1=1
525
526	nf2=nf+nb_face/2
527
528	ind=ind+1
529	vertex_glo(1,jjm_glo,nf)%ind=ind
530	vertex_glo(iim_glo,1,nfp1)%ind=ind
531	vertex_glo(iim_glo,jjm_glo,nf2)%ind=ind
532	cell_glo(ind)%assign_face=nf
533	cell_glo(ind)%assign_i=1
534	cell_glo(ind)%assign_j=jjm_glo
535	ENDDO
536
537	! vertex sud (sur 3 faces)
538
539	DO nf=nb_face/2+1,nb_face
540
541	nfp1=nf+1
542	IF (nfp1==nb_face+1) nfp1=nb_face/2+1
543
544	nf2=nf-nb_face/2+1
545	IF (nf2==nb_face/2+1) nf2=1
546
547	ind=ind+1
548	vertex_glo(iim_glo,1,nf)%ind=ind
549	vertex_glo(1,jjm_glo,nfp1)%ind=ind
550	vertex_glo(iim_glo,jjm_glo,nf2)%ind=ind
551	cell_glo(ind)%assign_face=nf
552	cell_glo(ind)%assign_i=iim_glo
553	cell_glo(ind)%assign_j=1
554	ENDDO
555
556
557	! vertex nord (sur 5 faces)
558
559	ind=ind+1
560	vertex_glo(1,1,1)%ind=ind
561	vertex_glo(1,1,2)%ind=ind
562	vertex_glo(1,1,3)%ind=ind
563	vertex_glo(1,1,4)%ind=ind
564	vertex_glo(1,1,5)%ind=ind
565	cell_glo(ind)%assign_face=1
566	cell_glo(ind)%assign_i=1
567	cell_glo(ind)%assign_j=1
568
569	! vertex sud (sur 5 faces)
570
571	ind=ind+1
572	vertex_glo(1,1,6)%ind=ind
573	vertex_glo(1,1,7)%ind=ind
574	vertex_glo(1,1,8)%ind=ind
575	vertex_glo(1,1,9)%ind=ind
576	vertex_glo(1,1,10)%ind=ind
577	cell_glo(ind)%assign_face=6
578	cell_glo(ind)%assign_i=1
579	cell_glo(ind)%assign_j=1
580
581	! assignation des coordonnees xyz
582
583	DO ind=1,ncell_glo
584	nf=cell_glo(ind)%assign_face
585	i=cell_glo(ind)%assign_i
586	j=cell_glo(ind)%assign_j
587	cell_glo(ind)%xyz=vertex_glo(i,j,nf)%xyz
588	ENDDO
589
590	!!!!!!!!!!!!!!!!!!!!!!!!!!!
591	! Association des voisins !
592	!!!!!!!!!!!!!!!!!!!!!!!!!!!
593
594	! interieur des faces
595
596	DO nf=1,nb_face
597	DO j=2,jjm_glo-1
598	DO i=2,iim_glo-1
599	ind=vertex_glo(i,j,nf)%ind
600	! right
601	ind2=vertex_glo(i+1,j,nf)%ind
602	nf2=cell_glo(ind2)%assign_face
603	cell_glo(ind)%neighbour(right-1)=ind2
604	cell_glo(ind2)%neighbour(tab_index(nf,nf2,left-1))=ind
605
606	! rup
607	ind2=vertex_glo(i,j+1,nf)%ind
608	nf2=cell_glo(ind2)%assign_face
609	cell_glo(ind)%neighbour(rup-1)=ind2
610	cell_glo(ind2)%neighbour(tab_index(nf,nf2,ldown-1))=ind
611
612	! lup
613	ind2=vertex_glo(i-1,j+1,nf)%ind
614	nf2=cell_glo(ind2)%assign_face
615	cell_glo(ind)%neighbour(lup-1)=ind2
616	cell_glo(ind2)%neighbour(tab_index(nf,nf2,rdown-1))=ind
617
618	! left
619	ind2=vertex_glo(i-1,j,nf)%ind
620	nf2=cell_glo(ind2)%assign_face
621	cell_glo(ind)%neighbour(left-1)=ind2
622	cell_glo(ind2)%neighbour(tab_index(nf,nf2,right-1))=ind
623	! ldown
624	ind2=vertex_glo(i,j-1,nf)%ind
625	nf2=cell_glo(ind2)%assign_face
626	cell_glo(ind)%neighbour(ldown-1)=ind2
627	cell_glo(ind2)%neighbour(tab_index(nf,nf2,rup-1))=ind
628
629	! rdown
630	ind2=vertex_glo(i+1,j-1,nf)%ind
631	nf2=cell_glo(ind2)%assign_face
632	cell_glo(ind)%neighbour(rdown-1)=ind2
633	cell_glo(ind2)%neighbour(tab_index(nf,nf2,lup-1))=ind
634	ENDDO
635	ENDDO
636	ENDDO
637
638	! frontiere face nord
639
640	DO nf=1,nb_face/2
641	DO i=2,iim_glo-1
642	j=1
643	ind=vertex_glo(i,j,nf)%ind
644	! right
645	ind2=vertex_glo(i+1,j,nf)%ind
646	nf2=cell_glo(ind2)%assign_face
647	cell_glo(ind)%neighbour(right-1)=ind2
648	cell_glo(ind2)%neighbour(tab_index(nf,nf2,left-1))=ind
649	! left
650	ind2=vertex_glo(i-1,j,nf)%ind
651	nf2=cell_glo(ind2)%assign_face
652	cell_glo(ind)%neighbour(left-1)=ind2
653	cell_glo(ind2)%neighbour(tab_index(nf,nf2,right-1))=ind
654	! lup
655	ind2=vertex_glo(i-1,j+1,nf)%ind
656	nf2=cell_glo(ind2)%assign_face
657	cell_glo(ind)%neighbour(lup-1)=ind2
658	cell_glo(ind2)%neighbour(tab_index(nf,nf2,rdown-1))=ind
659
660	ENDDO
661
662	DO i=2,iim_glo-1
663	j=jjm_glo
664	ind=vertex_glo(i,j,nf)%ind
665	! right
666	ind2=vertex_glo(i+1,j,nf)%ind
667	nf2=cell_glo(ind2)%assign_face
668	cell_glo(ind)%neighbour(right-1)=ind2
669	cell_glo(ind2)%neighbour(tab_index(nf,nf2,left-1))=ind
670	! left
671	ind2=vertex_glo(i-1,j,nf)%ind
672	nf2=cell_glo(ind2)%assign_face
673	cell_glo(ind)%neighbour(left-1)=ind2
674	cell_glo(ind2)%neighbour(tab_index(nf,nf2,right-1))=ind
675	! rdown
676	ind2=vertex_glo(i+1,j-1,nf)%ind
677	nf2=cell_glo(ind2)%assign_face
678	cell_glo(ind)%neighbour(rdown-1)=ind2
679	cell_glo(ind2)%neighbour(tab_index(nf,nf2,lup-1))=ind
680	ENDDO
681	ENDDO
682
683
684	! frontiere face sud
685
686	DO nf=nb_face/2+1,nb_face
687	DO i=2,iim_glo-1
688	j=1
689	ind=vertex_glo(i,j,nf)%ind
690	! right
691	ind2=vertex_glo(i+1,j,nf)%ind
692	nf2=cell_glo(ind2)%assign_face
693	cell_glo(ind)%neighbour(right-1)=ind2
694	cell_glo(ind2)%neighbour(tab_index(nf,nf2,left-1))=ind
695	! left
696	ind2=vertex_glo(i-1,j,nf)%ind
697	nf2=cell_glo(ind2)%assign_face
698	cell_glo(ind)%neighbour(left-1)=ind2
699	cell_glo(ind2)%neighbour(tab_index(nf,nf2,right-1))=ind
700	! lup
701	ind2=vertex_glo(i-1,j+1,nf)%ind
702	nf2=cell_glo(ind2)%assign_face
703	cell_glo(ind)%neighbour(lup-1)=ind2
704	cell_glo(ind2)%neighbour(tab_index(nf,nf2,rdown-1))=ind
705	ENDDO
706
707	DO j=2,jjm_glo-1
708	i=iim_glo
709	ind=vertex_glo(i,j,nf)%ind
710	! rup
711	ind2=vertex_glo(i,j+1,nf)%ind
712	nf2=cell_glo(ind2)%assign_face
713	cell_glo(ind)%neighbour(rup-1)=ind2
714	cell_glo(ind2)%neighbour(tab_index(nf,nf2,ldown-1))=ind
715	! ldown
716	ind2=vertex_glo(i,j-1,nf)%ind
717	nf2=cell_glo(ind2)%assign_face
718	cell_glo(ind)%neighbour(ldown-1)=ind2
719	cell_glo(ind2)%neighbour(tab_index(nf,nf2,rup-1))=ind
720	! lup
721	ind2=vertex_glo(i-1,j+1,nf)%ind
722	nf2=cell_glo(ind2)%assign_face
723	cell_glo(ind)%neighbour(lup-1)=ind2
724	cell_glo(ind2)%neighbour(tab_index(nf,nf2,rdown-1))=ind
725
726	ENDDO
727	ENDDO
728
729
730	! vertex nord (sur 3 faces)
731
732	DO nf=1,nb_face/2
733	ind=vertex_glo(1,jjm_glo,nf)%ind
734	cell_glo(ind)%neighbour(2)=cell_glo(ind)%neighbour(1)
735	ENDDO
736
737	! vertex sud (sur 3 faces)
738
739	DO nf=nb_face/2+1,nb_face
740
741	ind=vertex_glo(iim_glo,1,nf)%ind
742	cell_glo(ind)%neighbour(5)=cell_glo(ind)%neighbour(0)
743	ENDDO
744
745
746	! vertex nord (sur 5 faces)
747
748	ind=vertex_glo(1,1,1)%ind
749	cell_glo(ind)%neighbour(3)=cell_glo(ind)%neighbour(4)
750
751	! vertex sud (sur 5 faces)
752
753	ind=vertex_glo(1,1,6)%ind
754	cell_glo(ind)%neighbour(3)=cell_glo(ind)%neighbour(2)
755
756
757	ind=vertex_glo(1,jjm_glo,1)%ind
758	DO i=0,5
759	ind2=cell_glo(ind)%neighbour(i)
760	IF (ind2>0) PRINT *,"neighbour",i,cell_glo(ind2)%assign_face,cell_glo(ind2)%assign_i,cell_glo(ind2)%assign_j
761	ENDDO
762
763
764	!! assignation des delta
765	DO nf=1,nb_face
766	DO j=1,jjm_glo
767	DO i=1,iim_glo
768	ind=vertex_glo(i,j,nf)%ind
769	nf2=cell_glo(ind)%assign_face
770	vertex_glo(i,j,nf)%delta=tab_index(nf,nf2,0)
771	ENDDO
772	ENDDO
773	ENDDO
774
775
776	END SUBROUTINE set_cell
777
778
779	SUBROUTINE set_cell_edge
780	IMPLICIT NONE
781	INTEGER :: ind,k,ind2,ng,ne
782
783	DO ind=1,ncell_glo
784	cell_glo(ind)%edge(:)=0
785	ENDDO
786
787	ne=0
788	DO ind=1,ncell_glo
789	DO k=0,5
790	IF (cell_glo(ind)%edge(k)==0) THEN
791	ind2=cell_glo(ind)%neighbour(k)
792	DO ng=0,5
793	IF (cell_glo(ind2)%neighbour(ng)==ind) THEN
794	IF (cell_glo(ind2)%edge(ng)==0) THEN
795	ne=ne+1
796	cell_glo(ind)%edge(k)=ne
797	edge_glo(ne)%e1=ind
798	edge_glo(ne)%e2=ind2
799	cell_glo(ind2)%edge(ng)=ne
800	ELSE
801	ne=cell_glo(ind2)%edge(ng)
802	cell_glo(ind)%edge(k)=ne
803	ENDIF
804	EXIT
805	ENDIF
806	ENDDO
807	ENDIF
808	ENDDO
809	ENDDO
810
811	END SUBROUTINE set_cell_edge
812
813	SUBROUTINE set_vertex_edge
814	IMPLICIT NONE
815	INTEGER :: nf
816	INTEGER :: i,j
817	INTEGER :: ind,ind2
818	INTEGER :: ne
819	INTEGER :: k,l
820
821	DO nf=1,nb_face
822	DO j=0,jjm_glo+1
823	DO i=0,iim_glo+1
824	ind=vertex_glo(i,j,nf)%ind
825	DO k=0,5
826	ind2=vertex_glo(i,j,nf)%neighbour(k)
827	DO l=0,5
828	ne=cell_glo(ind)%edge(l)
829	IF (edge_glo(ne)%e1==ind .AND. edge_glo(ne)%e2==ind2) THEN
830	vertex_glo(i,j,nf)%ne(k)=1
831	ELSE IF (edge_glo(ne)%e1==ind2 .AND. edge_glo(ne)%e2==ind) THEN
832	vertex_glo(i,j,nf)%ne(k)=-1
833	ENDIF
834	ENDDO
835	ENDDO
836	ENDDO
837	ENDDO
838	ENDDO
839	END SUBROUTINE set_vertex_edge
840
841	SUBROUTINE compute_metric
842	IMPLICIT NONE
843
844	CALL allocate_metric
845	! CALL compute_face
846	CALL compute_face_projection
847	CALL set_index
848	CALL set_cell
849	CALL compute_extended_face_bis
850	CALL set_cell_edge
851	CALL set_vertex_edge
852
853	END SUBROUTINE compute_metric
854
855	END MODULE metric

Note: See TracBrowser for help on using the repository browser.

Download in other formats: