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

Last change on this file since 131 was 131, checked in by ymipsl, 11 years ago

Some operations must be only done by the mpi master task.

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