source: codes/icosagcm/trunk/src/sphere/metric.f90

Last change on this file was 899, checked in by adurocher, 5 years ago

trunk : Fixed GCC warnings

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