source: codes/icosagcm/trunk/src/geometry.f90 @ 195

Last change on this file since 195 was 186, checked in by ymipsl, 10 years ago

Add new openMP parallelism based on distribution of domains on threads. There is no more limitation of number of threads by MPI process.

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)
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.