Context Navigation

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

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

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

Download in other formats: