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

Last change on this file since 810 was 711, checked in by ymipsl, 6 years ago

Adding halo transfer for scalar field ad vorticity point.
=> Use request "req_z1_scal" with transfer function

YM

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