source: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/geometry.f90 @ 267

Last change on this file since 267 was 260, checked in by ymipsl, 10 years ago

Implement restartability for dynamico

YM
File size: 17.2 KB
Line 
1MODULE geometry
2  USE field_mod
3 
4  TYPE t_geometry
5    TYPE(t_field),POINTER :: xyz_i(:)
6    TYPE(t_field),POINTER :: centroid(:)
7    TYPE(t_field),POINTER :: xyz_e(:)
8    TYPE(t_field),POINTER :: xyz_v(:)
9    TYPE(t_field),POINTER :: ep_e(:)
10    TYPE(t_field),POINTER :: et_e(:)
11    TYPE(t_field),POINTER :: elon_i(:)
12    TYPE(t_field),POINTER :: elat_i(:)
13    TYPE(t_field),POINTER :: elon_e(:)
14    TYPE(t_field),POINTER :: elat_e(:)
15    TYPE(t_field),POINTER :: Ai(:)
16    TYPE(t_field),POINTER :: Av(:)
17    TYPE(t_field),POINTER :: de(:)
18    TYPE(t_field),POINTER :: le(:)
19    TYPE(t_field),POINTER :: Riv(:)
20    TYPE(t_field),POINTER :: Riv2(:)
21    TYPE(t_field),POINTER :: ne(:)
22    TYPE(t_field),POINTER :: Wee(:)
23    TYPE(t_field),POINTER :: bi(:)
24    TYPE(t_field),POINTER :: fv(:)
25  END TYPE t_geometry   
26
27  TYPE(t_geometry),SAVE,TARGET :: geom
28
29 
30  REAL(rstd),POINTER :: Ai(:)          ! area of a cell
31!$OMP THREADPRIVATE(Ai)
32  REAL(rstd),POINTER :: xyz_i(:,:)     ! coordinate of the center of the cell (voronoi)
33!$OMP THREADPRIVATE(xyz_i)
34  REAL(rstd),POINTER :: centroid(:,:)  ! coordinate of the centroid of the cell
35!$OMP THREADPRIVATE(centroid)
36  REAL(rstd),POINTER :: xyz_e(:,:)     ! coordinate of a wind point on the cell on a edge
37!$OMP THREADPRIVATE(xyz_e)
38  REAL(rstd),POINTER :: ep_e(:,:)      ! perpendicular unit vector of a edge (outsider)
39!$OMP THREADPRIVATE(ep_e)
40  REAL(rstd),POINTER :: et_e(:,:)      ! tangeantial unit vector of a edge
41!$OMP THREADPRIVATE(et_e)
42  REAL(rstd),POINTER :: elon_i(:,:)    ! unit longitude vector on the center
43!$OMP THREADPRIVATE(elon_i)
44  REAL(rstd),POINTER :: elat_i(:,:)    ! unit latitude vector on the center
45!$OMP THREADPRIVATE(elat_i)
46  REAL(rstd),POINTER :: elon_e(:,:)    ! unit longitude vector on a wind point
47!$OMP THREADPRIVATE(elon_e)
48  REAL(rstd),POINTER :: elat_e(:,:)    ! unit latitude vector on a wind point
49!$OMP THREADPRIVATE(elat_e)
50  REAL(rstd),POINTER :: xyz_v(:,:)     ! coordinate of a vertex (center of the dual mesh)
51!$OMP THREADPRIVATE(xyz_v)
52  REAL(rstd),POINTER :: Av(:)          ! area of dual mesk cell
53!$OMP THREADPRIVATE(Av)
54  REAL(rstd),POINTER :: de(:)          ! distance from a neighbour == lenght of an edge of the dual mesh
55!$OMP THREADPRIVATE(de)
56  REAL(rstd),POINTER :: le(:)          ! lenght of a edge
57!$OMP THREADPRIVATE(le)
58  REAL(rstd),POINTER :: Riv(:,:)       ! weight
59!$OMP THREADPRIVATE(Riv)
60  REAL(rstd),POINTER :: Riv2(:,:)      ! weight
61!$OMP THREADPRIVATE(Riv2)
62  INTEGER,POINTER    :: ne(:,:)        ! convention for the way on the normal wind on an edge
63!$OMP THREADPRIVATE(ne)
64  REAL(rstd),POINTER :: Wee(:,:,:)     ! weight
65!$OMP THREADPRIVATE(Wee)
66  REAL(rstd),POINTER :: bi(:)          ! orographie
67!$OMP THREADPRIVATE(bi)
68  REAL(rstd),POINTER :: fv(:)          ! coriolis (evaluted on a vertex)
69!$OMP THREADPRIVATE(fv)
70     
71
72  INTEGER, PARAMETER :: ne_right=1
73  INTEGER, PARAMETER :: ne_rup=-1
74  INTEGER, PARAMETER :: ne_lup=1
75  INTEGER, PARAMETER :: ne_left=-1
76  INTEGER, PARAMETER :: ne_ldown=1
77  INTEGER, PARAMETER :: ne_rdown=-1
78
79CONTAINS
80
81  SUBROUTINE allocate_geometry
82  USE field_mod
83  IMPLICIT NONE
84 
85    CALL allocate_field(geom%Ai,field_t,type_real,name='Ai')
86    CALL allocate_field(geom%xyz_i,field_t,type_real,3)
87    CALL allocate_field(geom%centroid,field_t,type_real,3)
88    CALL allocate_field(geom%xyz_e,field_u,type_real,3)
89    CALL allocate_field(geom%ep_e,field_u,type_real,3)
90    CALL allocate_field(geom%et_e,field_u,type_real,3)
91    CALL allocate_field(geom%elon_i,field_t,type_real,3)
92    CALL allocate_field(geom%elat_i,field_t,type_real,3)
93    CALL allocate_field(geom%elon_e,field_u,type_real,3)
94    CALL allocate_field(geom%elat_e,field_u,type_real,3)
95    CALL allocate_field(geom%xyz_v,field_z,type_real,3)
96    CALL allocate_field(geom%de,field_u,type_real)
97    CALL allocate_field(geom%le,field_u,type_real)
98    CALL allocate_field(geom%bi,field_t,type_real)
99    CALL allocate_field(geom%Av,field_z,type_real)
100    CALL allocate_field(geom%Riv,field_t,type_real,6)
101    CALL allocate_field(geom%Riv2,field_t,type_real,6)
102    CALL allocate_field(geom%ne,field_t,type_integer,6)
103    CALL allocate_field(geom%Wee,field_u,type_real,5,2)
104    CALL allocate_field(geom%bi,field_t,type_real)
105    CALL allocate_field(geom%fv,field_z,type_real)
106
107  END SUBROUTINE allocate_geometry
108 
109 
110  SUBROUTINE swap_geometry(ind)
111  USE field_mod
112  IMPLICIT NONE
113    INTEGER,INTENT(IN) :: ind
114!!$OMP MASTER
115    Ai=geom%Ai(ind)
116    xyz_i=geom%xyz_i(ind)
117    centroid=geom%centroid(ind)
118    xyz_e=geom%xyz_e(ind)
119    ep_e=geom%ep_e(ind)
120    et_e=geom%et_e(ind)
121    elon_i=geom%elon_i(ind)
122    elat_i=geom%elat_i(ind)
123    elon_e=geom%elon_e(ind)
124    elat_e=geom%elat_e(ind)
125    xyz_v=geom%xyz_v(ind)
126    de=geom%de(ind)
127    le=geom%le(ind)
128    Av=geom%Av(ind)
129    Riv=geom%Riv(ind)
130    Riv2=geom%Riv2(ind)
131    ne=geom%ne(ind)
132    Wee=geom%Wee(ind)
133    bi=geom%bi(ind)
134    fv=geom%fv(ind)
135!!$OMP END MASTER
136!!$OMP BARRIER   
137  END SUBROUTINE swap_geometry
138
139  SUBROUTINE update_circumcenters 
140    USE domain_mod
141    USE dimensions
142    USE spherical_geom_mod
143    USE vector
144    USE transfert_mod
145
146    IMPLICIT NONE
147    REAL(rstd) :: x1(3),x2(3)
148    REAL(rstd) :: vect(3,6)
149    REAL(rstd) :: centr(3)
150    INTEGER :: ind,i,j,n,k
151    TYPE(t_message),SAVE :: message
152    LOGICAL, SAVE :: first=.TRUE.
153!$OMP THREADPRIVATE(first)   
154   
155    IF (first) THEN
156      CALL init_message(geom%xyz_i, req_i1 ,message)
157      first=.FALSE.
158    ENDIF
159   
160    CALL transfert_message(geom%xyz_i,message)
161       
162   
163!    CALL transfert_request(geom%xyz_i,req_i1)
164
165    DO ind=1,ndomain
166      IF (.NOT. assigned_domain(ind)) CYCLE
167      CALL swap_dimensions(ind)
168      CALL swap_geometry(ind)
169      DO j=jj_begin,jj_end
170        DO i=ii_begin,ii_end
171          n=(j-1)*iim+i
172          DO k=0,5
173            x1(:) = xyz_i(n+t_pos(k+1),:)
174            x2(:) = xyz_i(n+t_pos(MOD(k+1,6)+1),:)
175            if (norm(x1-x2)<1e-16) x2(:) = xyz_i(n+t_pos(MOD(k+2,6)+1),:)
176            CALL circumcenter(xyz_i(n,:), x1, x2, xyz_v(n+z_pos(k+1),:))
177          ENDDO
178        ENDDO
179      ENDDO
180    ENDDO
181
182  END SUBROUTINE update_circumcenters
183
184  SUBROUTINE optimize_geometry
185  USE metric
186  USE spherical_geom_mod
187  USE domain_mod
188  USE dimensions
189  USE transfert_mod
190  USE vector
191  USE getin_mod
192  IMPLICIT NONE
193    INTEGER :: nb_it=0
194    TYPE(t_domain),POINTER :: d
195    INTEGER :: ind,it,i,j,n,k
196    REAL(rstd) :: x1(3),x2(3)
197    REAL(rstd) :: vect(3,6)
198    REAL(rstd) :: centr(3)
199    REAL(rstd) :: sum
200    LOGICAL    :: check
201   
202   
203    CALL getin('optim_it',nb_it)
204   
205    DO ind=1,ndomain
206      IF (.NOT. assigned_domain(ind)) CYCLE
207      d=>domain(ind)
208      CALL swap_dimensions(ind)
209      CALL swap_geometry(ind)
210      DO j=jj_begin,jj_end
211        DO i=ii_begin,ii_end
212          n=(j-1)*iim+i
213          xyz_i(n,:)=d%xyz(:,i,j) 
214        ENDDO
215      ENDDO
216    ENDDO
217   
218    CALL update_circumcenters
219
220    DO ind=1,ndomain
221      IF (.NOT. assigned_domain(ind)) CYCLE
222      d=>domain(ind)
223      CALL swap_dimensions(ind)
224      CALL swap_geometry(ind)
225      DO j=jj_begin,jj_end
226        DO i=ii_begin,ii_end
227          n=(j-1)*iim+i
228          DO k=0,5
229            x1(:) = xyz_v(n+z_pos(k+1),:)
230            x2(:) = d%vertex(:,k,i,j) 
231            IF (norm(x1-x2)>1e-10) THEN
232              PRINT*,"vertex diff ",ind,i,j,k
233              PRINT*,x1
234              PRINT*,x2
235            ENDIF
236          ENDDO
237        ENDDO
238      ENDDO
239    ENDDO
240   
241   
242    DO it=1,nb_it
243      IF (MOD(it,100)==0) THEN
244        check=.TRUE.
245      ELSE
246        check=.FALSE.
247      ENDIF
248     
249      sum=0
250      DO ind=1,ndomain
251      IF (.NOT. assigned_domain(ind)) CYCLE
252        CALL swap_dimensions(ind)
253        CALL swap_geometry(ind)
254        DO j=jj_begin,jj_end
255          DO i=ii_begin,ii_end
256            n=(j-1)*iim+i
257            vect(:,1)=xyz_v(n+z_rup,:)
258            vect(:,2)=xyz_v(n+z_up,:)
259            vect(:,3)=xyz_v(n+z_lup,:)
260            vect(:,4)=xyz_v(n+z_ldown,:)
261            vect(:,5)=xyz_v(n+z_down,:)
262            vect(:,6)=xyz_v(n+z_rdown,:)
263            CALL compute_centroid(vect,6,centr)
264            IF (check) THEN
265              sum=MAX(sum,norm(xyz_i(n,:)-centr(:)))
266            ENDIF
267            xyz_i(n,:)=centr(:)
268          ENDDO
269        ENDDO
270       
271      ENDDO
272
273      IF (check) THEN
274        PRINT *,"it = ",it,"  diff centroid circumcenter ",sum
275      ENDIF
276
277     CALL update_circumcenters
278
279    ENDDO
280   
281  END SUBROUTINE optimize_geometry
282
283  SUBROUTINE set_geometry
284  USE metric
285  USE vector
286  USE spherical_geom_mod
287  USE domain_mod
288  USE dimensions
289  USE transfert_mod
290  USE getin_mod
291  IMPLICIT NONE
292
293    REAL(rstd) :: surf(6)
294    REAL(rstd) :: surf_v(6)
295    REAL(rstd) :: vect(3,6)
296    REAL(rstd) :: centr(3)
297    REAL(rstd) :: vet(3),vep(3), vertex(3)
298    INTEGER :: ind,i,j,k,n
299    TYPE(t_domain),POINTER :: d
300    REAL(rstd) ::  S
301    REAL(rstd) :: w(6)
302    REAL(rstd) :: lon,lat
303    INTEGER :: ii_glo,jj_glo
304    REAL(rstd) :: S1,S2
305         
306     
307    CALL optimize_geometry
308 
309    DO ind=1,ndomain
310      IF (.NOT. assigned_domain(ind)) CYCLE
311      d=>domain(ind)
312      CALL swap_dimensions(ind)
313      CALL swap_geometry(ind)
314      DO j=jj_begin-1,jj_end+1
315        DO i=ii_begin-1,ii_end+1
316          n=(j-1)*iim+i
317
318          DO k=0,5
319            ne(n,k+1)=d%ne(k,i,j)
320          ENDDO
321
322          vect(:,1)=xyz_v(n+z_rup,:)
323          vect(:,2)=xyz_v(n+z_up,:)
324          vect(:,3)=xyz_v(n+z_lup,:)
325          vect(:,4)=xyz_v(n+z_ldown,:)
326          vect(:,5)=xyz_v(n+z_down,:)
327          vect(:,6)=xyz_v(n+z_rdown,:)
328          CALL compute_centroid(vect,6,centr)
329          centroid(n,:)=centr(:)
330
331               
332          CALL xyz2lonlat(xyz_v(n+z_up,:),lon,lat)
333          fv(n+z_up)=2*sin(lat)*omega
334          CALL xyz2lonlat(xyz_v(n+z_down,:),lon,lat)
335          fv(n+z_down)=2*sin(lat)*omega
336         
337          bi(n)=0. 
338       
339          CALL dist_cart(xyz_i(n,:),xyz_i(n+t_right,:),de(n+u_right))
340          CALL dist_cart(xyz_i(n,:),xyz_i(n+t_lup,:),de(n+u_lup))
341          CALL dist_cart(xyz_i(n,:),xyz_i(n+t_ldown,:),de(n+u_ldown))
342         
343          CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_right,:),0.5,xyz_e(n+u_right,:))
344          CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_lup,:),0.5,xyz_e(n+u_lup,:))
345          CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_ldown,:),0.5,xyz_e(n+u_ldown,:))
346
347          CALL dist_cart(xyz_v(n+z_rdown,:), xyz_v(n+z_rup,:),le(n+u_right))
348          CALL dist_cart(xyz_v(n+z_up,:), xyz_v(n+z_lup,:),le(n+u_lup))
349          CALL dist_cart(xyz_v(n+z_ldown,:), xyz_v(n+z_down,:),le(n+u_ldown))
350         
351          Ai(n)=0
352          DO k=0,5
353            CALL surf_triangle(xyz_i(n,:),xyz_i(n+t_pos(k+1),:),xyz_i(n+t_pos(MOD((k+1+6),6)+1),:),surf_v(k+1))
354            CALL surf_triangle(xyz_i(n,:),xyz_v(n+z_pos(MOD((k-1+6),6)+1),:),xyz_v(n+z_pos(k+1),:),surf(k+1))
355            Ai(n)=Ai(n)+surf(k+1)
356            IF (i==ii_end .AND. j==jj_begin) THEN
357              IF (Ai(n)<1e20) THEN
358              ELSE
359                PRINT *,"PB !!",Ai(n),k,surf(k+1)
360                PRINT*,xyz_i(n,:),xyz_v(n+z_pos(MOD((k-1+6),6)+1),:),xyz_v(n+z_pos(k+1),:)
361              ENDIF
362            ENDIF
363          ENDDO
364
365          ! Sign convention : Ringler et al., JCP 2010, eq. 21 p. 3071
366          ! Normal component is along outgoing normal vector if ne=1
367
368          CALL cross_product2(xyz_v(n+z_rdown,:),xyz_v(n+z_rup,:),vep)
369          IF (norm(vep)>1e-30) THEN
370            vep(:)=vep(:)/norm(vep)                         ! Inward normal vector
371            CALL cross_product2(vep,xyz_e(n+u_right,:),vet) ! Counter-clockwise tangent vector
372            vet(:)=vet(:)/norm(vet)
373            ep_e(n+u_right,:)=-vep(:)*ne(n,right)
374            et_e(n+u_right,:)=vet(:)*ne(n,right)
375          ENDIF
376
377          CALL cross_product2(xyz_v(n+z_up,:),xyz_v(n+z_lup,:),vep)
378          IF (norm(vep)>1e-30) THEN
379            vep(:)=vep(:)/norm(vep)
380            CALL cross_product2(vep,xyz_e(n+u_lup,:),vet)
381            vet(:)=vet(:)/norm(vet)
382            ep_e(n+u_lup,:)=-vep(:)*ne(n,lup)
383            et_e(n+u_lup,:)=vet(:)*ne(n,lup)
384          ENDIF
385
386          CALL cross_product2(xyz_v(n+z_ldown,:),xyz_v(n+z_down,:),vep)
387          IF (norm(vep)>1e-30) THEN
388            vep(:)=vep(:)/norm(vep)
389            CALL cross_product2(vep,xyz_e(n+u_ldown,:),vet)
390            vet(:)=vet(:)/norm(vet)
391            ep_e(n+u_ldown,:)=-vep(:)*ne(n,ldown)
392            et_e(n+u_ldown,:)=vet(:)*ne(n,ldown)
393          ENDIF
394
395          CALL xyz2lonlat(xyz_i(n,:),lon,lat)
396          elon_i(n,1) = -sin(lon)
397          elon_i(n,2) = cos(lon)
398          elon_i(n,3) = 0
399          elat_i(n,1) = -cos(lon)*sin(lat)
400          elat_i(n,2) = -sin(lon)*sin(lat)
401          elat_i(n,3) = cos(lat)
402
403         
404          CALL xyz2lonlat(xyz_e(n+u_right,:),lon,lat)
405          elon_e(n+u_right,1) = -sin(lon)
406          elon_e(n+u_right,2) = cos(lon)
407          elon_e(n+u_right,3) = 0
408          elat_e(n+u_right,1) = -cos(lon)*sin(lat)
409          elat_e(n+u_right,2) = -sin(lon)*sin(lat)
410          elat_e(n+u_right,3) = cos(lat)
411         
412          CALL xyz2lonlat(xyz_e(n+u_lup,:),lon,lat)
413          elon_e(n+u_lup,1) = -sin(lon)
414          elon_e(n+u_lup,2) = cos(lon)
415          elon_e(n+u_lup,3) = 0
416          elat_e(n+u_lup,1) = -cos(lon)*sin(lat)
417          elat_e(n+u_lup,2) = -sin(lon)*sin(lat)
418          elat_e(n+u_lup,3) = cos(lat)
419 
420          CALL xyz2lonlat(xyz_e(n+u_ldown,:),lon,lat)
421          elon_e(n+u_ldown,1) = -sin(lon)
422          elon_e(n+u_ldown,2) = cos(lon)
423          elon_e(n+u_ldown,3) = 0
424          elat_e(n+u_ldown,1) = -cos(lon)*sin(lat)
425          elat_e(n+u_ldown,2) = -sin(lon)*sin(lat)
426          elat_e(n+u_ldown,3) = cos(lat)
427 
428         
429          DO k=0,5
430            CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(k+1),:),S1)
431            CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(MOD(k+1+6,6)+1),:),S2)
432            Riv(n,k+1)=0.5*(S1+S2)/Ai(n)
433            Riv2(n,k+1)=0.5*(S1+S2)/surf_v(k+1)
434          ENDDO
435
436          DO k=1,6
437            IF (ABS(surf_v(k))<1e-30) THEN
438              Riv(n,k)=0.
439            ENDIF
440          ENDDO
441         
442          Av(n+z_up)=surf_v(vup)+1e-100
443          Av(n+z_down)=surf_v(vdown)+1e-100
444       
445        ENDDO
446      ENDDO
447     
448      DO j=jj_begin,jj_end
449        DO i=ii_begin,ii_end
450          n=(j-1)*iim+i
451   
452          CALL compute_wee(n,right,w)
453          Wee(n+u_right,:,1)=w(1:5)
454
455          CALL compute_wee(n+t_right,left,w)
456          Wee(n+u_right,:,2)=w(1:5)
457
458
459          CALL compute_wee(n,lup,w)
460          Wee(n+u_lup,:,1)=w(1:5)
461
462          CALL compute_wee(n+t_lup,rdown,w)
463          Wee(n+u_lup,:,2)=w(1:5)
464
465
466          CALL compute_wee(n,ldown,w)
467          Wee(n+u_ldown,:,1)=w(1:5)
468
469          CALL compute_wee(n+t_ldown,rup,w)
470          Wee(n+u_ldown,:,2)=w(1:5)
471
472        ENDDO
473      ENDDO
474     
475      DO j=jj_begin,jj_end
476        DO i=ii_begin,ii_end
477          n=(j-1)*iim+i
478          ii_glo=d%ii_begin_glo-d%ii_begin+i
479          jj_glo=d%jj_begin_glo-d%jj_begin+j
480         
481          IF (ii_glo==1 .AND. jj_glo==1) THEN
482            le(n+u_ldown)=0
483            xyz_v(n+z_ldown,:)=xyz_v(n+z_down,:)
484                       
485          ENDIF
486
487          IF (ii_glo==iim_glo .AND. jj_glo==1) THEN
488            le(n+u_right)=0
489            xyz_v(n+z_rdown,:)=xyz_v(n+z_rup,:)
490          ENDIF
491
492          IF (ii_glo==iim_glo .AND. jj_glo==jjm_glo) THEN
493            le(n+u_rup)=0
494            xyz_v(n+z_rup,:)=xyz_v(n+z_up,:)
495          ENDIF
496
497          IF (ii_glo==1 .AND. jj_glo==jjm_glo) THEN
498            le(n+u_lup)=0
499            xyz_v(n+z_up,:)=xyz_v(n+z_lup,:)
500          ENDIF
501         
502        ENDDO
503      ENDDO
504
505      DO j=jj_begin-1,jj_end+1
506        DO i=ii_begin-1,ii_end+1
507          n=(j-1)*iim+i
508          xyz_i(n,:)=xyz_i(n,:) * radius
509          xyz_v(n+z_up,:)=xyz_v(n+z_up,:) * radius
510          xyz_v(n+z_down,:)=xyz_v(n+z_down,:) *radius
511          de(n+u_right)=de(n+u_right) * radius
512          de(n+u_lup)=de(n+u_lup)*radius
513          de(n+u_ldown)=de(n+u_ldown)*radius
514          xyz_e(n+u_right,:)=xyz_e(n+u_right,:)*radius
515          xyz_e(n+u_lup,:)=xyz_e(n+u_lup,:)*radius
516          xyz_e(n+u_ldown,:)=xyz_e(n+u_ldown,:)*radius
517          le(n+u_right)=le(n+u_right)*radius
518          le(n+u_lup)=le(n+u_lup)*radius
519          le(n+u_ldown)=le(n+u_ldown)*radius
520          Ai(n)=Ai(n)*radius**2
521          Av(n+z_up)=Av(n+z_up)*radius**2
522          Av(n+z_down)=Av(n+z_down)*radius**2
523        ENDDO
524      ENDDO
525                 
526    ENDDO
527
528    CALL transfert_request(geom%Ai,req_i1)
529    CALL transfert_request(geom%centroid,req_i1)
530
531    CALL surf_triangle(d%xyz(:,ii_begin,jj_begin),d%xyz(:,ii_begin,jj_end),d%xyz(:,ii_end,jj_begin),S)
532 
533  END SUBROUTINE set_geometry
534 
535  SUBROUTINE compute_wee(n,pos,w)
536  IMPLICIT NONE
537    INTEGER,INTENT(IN) :: n
538    INTEGER,INTENT(IN) :: pos
539    REAL(rstd),INTENT(OUT) ::w(6)
540
541    REAL(rstd) :: ne_(0:5)
542    REAL(rstd) :: Riv_(6)
543    INTEGER    :: k
544   
545   
546    DO k=0,5
547      ne_(k)=ne(n,MOD(pos-1+k+6,6)+1) 
548      Riv_(k+1)=Riv(n,MOD(pos-1+k+6,6)+1)
549    ENDDO
550         
551    w(1)=-ne_(0)*ne_(1)*(Riv_(1)-0.5)
552    w(2)=-ne_(2)*(ne_(0)*Riv_(2)-w(1)*ne_(1))
553    w(3)=-ne_(3)*(ne_(0)*Riv_(3)-w(2)*ne_(2))
554    w(4)=-ne_(4)*(ne_(0)*Riv_(4)-w(3)*ne_(3))
555    w(5)=-ne_(5)*(ne_(0)*Riv_(5)-w(4)*ne_(4))
556    w(6)=ne_(0)*ne_(5)*(Riv_(6)-0.5)
557   
558!    IF ( ABS(w(5)-w(6))>1e-20) PRINT *, "pb pour wee : w(5)!=w(6)",sum(Riv_(:))
559
560   END SUBROUTINE compute_wee
561           
562
563 
564  SUBROUTINE compute_geometry
565  IMPLICIT NONE
566    CALL allocate_geometry
567    CALL set_geometry
568   
569  END SUBROUTINE compute_geometry
570 
571END MODULE geometry
Note: See TracBrowser for help on using the repository browser.