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

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

correction for compiling with gfortran (line too long)
improvement for splitting domain
Call twice transfert request for field u is no longer necessary

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