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

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

Grid optimisation : report correction of rev r294 in Satrun branch

YM
File size: 17.3 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 :: message1,message0
152    LOGICAL, SAVE :: first=.TRUE.
153!$OMP THREADPRIVATE(first)   
154   
155    IF (first) THEN
156      CALL init_message(geom%xyz_i, req_i1 ,message1)
157      CALL init_message(geom%xyz_i,req_i0, message0)
158      first=.FALSE.
159    ENDIF
160   
161    CALL transfert_message(geom%xyz_i,message0)
162    CALL transfert_message(geom%xyz_i,message1)
163       
164   
165!    CALL transfert_request(geom%xyz_i,req_i1)
166
167    DO ind=1,ndomain
168      IF (.NOT. assigned_domain(ind)) CYCLE
169      CALL swap_dimensions(ind)
170      CALL swap_geometry(ind)
171      DO j=jj_begin,jj_end
172        DO i=ii_begin,ii_end
173          n=(j-1)*iim+i
174          DO k=0,5
175            x1(:) = xyz_i(n+t_pos(k+1),:)
176            x2(:) = xyz_i(n+t_pos(MOD(k+1,6)+1),:)
177            if (norm(x1-x2)<1e-16) x2(:) = xyz_i(n+t_pos(MOD(k+2,6)+1),:)
178            CALL circumcenter(xyz_i(n,:), x1, x2, xyz_v(n+z_pos(k+1),:))
179          ENDDO
180        ENDDO
181      ENDDO
182    ENDDO
183
184  END SUBROUTINE update_circumcenters
185
186  SUBROUTINE optimize_geometry
187  USE metric
188  USE spherical_geom_mod
189  USE domain_mod
190  USE dimensions
191  USE transfert_mod
192  USE vector
193  USE getin_mod
194  IMPLICIT NONE
195    INTEGER :: nb_it=0
196    TYPE(t_domain),POINTER :: d
197    INTEGER :: ind,it,i,j,n,k
198    REAL(rstd) :: x1(3),x2(3)
199    REAL(rstd) :: vect(3,6)
200    REAL(rstd) :: centr(3)
201    REAL(rstd) :: sum
202    LOGICAL    :: check
203   
204   
205    CALL getin('optim_it',nb_it)
206   
207    DO ind=1,ndomain
208      IF (.NOT. assigned_domain(ind)) CYCLE
209      d=>domain(ind)
210      CALL swap_dimensions(ind)
211      CALL swap_geometry(ind)
212      DO j=jj_begin,jj_end
213        DO i=ii_begin,ii_end
214          n=(j-1)*iim+i
215          xyz_i(n,:)=d%xyz(:,i,j) 
216        ENDDO
217      ENDDO
218    ENDDO
219   
220    CALL update_circumcenters
221
222    DO ind=1,ndomain
223      IF (.NOT. assigned_domain(ind)) CYCLE
224      d=>domain(ind)
225      CALL swap_dimensions(ind)
226      CALL swap_geometry(ind)
227      DO j=jj_begin,jj_end
228        DO i=ii_begin,ii_end
229          n=(j-1)*iim+i
230          DO k=0,5
231            x1(:) = xyz_v(n+z_pos(k+1),:)
232            x2(:) = d%vertex(:,k,i,j) 
233            IF (norm(x1-x2)>1e-10) THEN
234              PRINT*,"vertex diff ",ind,i,j,k
235              PRINT*,x1
236              PRINT*,x2
237            ENDIF
238          ENDDO
239        ENDDO
240      ENDDO
241    ENDDO
242   
243   
244    DO it=1,nb_it
245      IF (MOD(it,100)==0) THEN
246        check=.TRUE.
247      ELSE
248        check=.FALSE.
249      ENDIF
250     
251      sum=0
252      DO ind=1,ndomain
253      IF (.NOT. assigned_domain(ind)) CYCLE
254        CALL swap_dimensions(ind)
255        CALL swap_geometry(ind)
256        DO j=jj_begin,jj_end
257          DO i=ii_begin,ii_end
258            n=(j-1)*iim+i
259            vect(:,1)=xyz_v(n+z_rup,:)
260            vect(:,2)=xyz_v(n+z_up,:)
261            vect(:,3)=xyz_v(n+z_lup,:)
262            vect(:,4)=xyz_v(n+z_ldown,:)
263            vect(:,5)=xyz_v(n+z_down,:)
264            vect(:,6)=xyz_v(n+z_rdown,:)
265            CALL compute_centroid(vect,6,centr)
266            IF (check) THEN
267              sum=MAX(sum,norm(xyz_i(n,:)-centr(:)))
268            ENDIF
269            xyz_i(n,:)=centr(:)
270          ENDDO
271        ENDDO
272       
273      ENDDO
274
275      IF (check) THEN
276        PRINT *,"it = ",it,"  diff centroid circumcenter ",sum
277      ENDIF
278
279     CALL update_circumcenters
280
281    ENDDO
282   
283  END SUBROUTINE optimize_geometry
284
285  SUBROUTINE set_geometry
286  USE metric
287  USE vector
288  USE spherical_geom_mod
289  USE domain_mod
290  USE dimensions
291  USE transfert_mod
292  USE getin_mod
293  IMPLICIT NONE
294
295    REAL(rstd) :: surf(6)
296    REAL(rstd) :: surf_v(6)
297    REAL(rstd) :: vect(3,6)
298    REAL(rstd) :: centr(3)
299    REAL(rstd) :: vet(3),vep(3), vertex(3)
300    INTEGER :: ind,i,j,k,n
301    TYPE(t_domain),POINTER :: d
302    REAL(rstd) ::  S
303    REAL(rstd) :: w(6)
304    REAL(rstd) :: lon,lat
305    INTEGER :: ii_glo,jj_glo
306    REAL(rstd) :: S1,S2
307         
308     
309    CALL optimize_geometry
310 
311    DO ind=1,ndomain
312      IF (.NOT. assigned_domain(ind)) CYCLE
313      d=>domain(ind)
314      CALL swap_dimensions(ind)
315      CALL swap_geometry(ind)
316      DO j=jj_begin-1,jj_end+1
317        DO i=ii_begin-1,ii_end+1
318          n=(j-1)*iim+i
319
320          DO k=0,5
321            ne(n,k+1)=d%ne(k,i,j)
322          ENDDO
323
324          vect(:,1)=xyz_v(n+z_rup,:)
325          vect(:,2)=xyz_v(n+z_up,:)
326          vect(:,3)=xyz_v(n+z_lup,:)
327          vect(:,4)=xyz_v(n+z_ldown,:)
328          vect(:,5)=xyz_v(n+z_down,:)
329          vect(:,6)=xyz_v(n+z_rdown,:)
330          CALL compute_centroid(vect,6,centr)
331          centroid(n,:)=centr(:)
332
333               
334          CALL xyz2lonlat(xyz_v(n+z_up,:),lon,lat)
335          fv(n+z_up)=2*sin(lat)*omega
336          CALL xyz2lonlat(xyz_v(n+z_down,:),lon,lat)
337          fv(n+z_down)=2*sin(lat)*omega
338         
339          bi(n)=0. 
340       
341          CALL dist_cart(xyz_i(n,:),xyz_i(n+t_right,:),de(n+u_right))
342          CALL dist_cart(xyz_i(n,:),xyz_i(n+t_lup,:),de(n+u_lup))
343          CALL dist_cart(xyz_i(n,:),xyz_i(n+t_ldown,:),de(n+u_ldown))
344         
345          CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_right,:),0.5,xyz_e(n+u_right,:))
346          CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_lup,:),0.5,xyz_e(n+u_lup,:))
347          CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_ldown,:),0.5,xyz_e(n+u_ldown,:))
348
349          CALL dist_cart(xyz_v(n+z_rdown,:), xyz_v(n+z_rup,:),le(n+u_right))
350          CALL dist_cart(xyz_v(n+z_up,:), xyz_v(n+z_lup,:),le(n+u_lup))
351          CALL dist_cart(xyz_v(n+z_ldown,:), xyz_v(n+z_down,:),le(n+u_ldown))
352         
353          Ai(n)=0
354          DO k=0,5
355            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))
356            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))
357            Ai(n)=Ai(n)+surf(k+1)
358            IF (i==ii_end .AND. j==jj_begin) THEN
359              IF (Ai(n)<1e20) THEN
360              ELSE
361                PRINT *,"PB !!",Ai(n),k,surf(k+1)
362                PRINT*,xyz_i(n,:),xyz_v(n+z_pos(MOD((k-1+6),6)+1),:),xyz_v(n+z_pos(k+1),:)
363              ENDIF
364            ENDIF
365          ENDDO
366
367          ! Sign convention : Ringler et al., JCP 2010, eq. 21 p. 3071
368          ! Normal component is along outgoing normal vector if ne=1
369
370          CALL cross_product2(xyz_v(n+z_rdown,:),xyz_v(n+z_rup,:),vep)
371          IF (norm(vep)>1e-30) THEN
372            vep(:)=vep(:)/norm(vep)                         ! Inward normal vector
373            CALL cross_product2(vep,xyz_e(n+u_right,:),vet) ! Counter-clockwise tangent vector
374            vet(:)=vet(:)/norm(vet)
375            ep_e(n+u_right,:)=-vep(:)*ne(n,right)
376            et_e(n+u_right,:)=vet(:)*ne(n,right)
377          ENDIF
378
379          CALL cross_product2(xyz_v(n+z_up,:),xyz_v(n+z_lup,:),vep)
380          IF (norm(vep)>1e-30) THEN
381            vep(:)=vep(:)/norm(vep)
382            CALL cross_product2(vep,xyz_e(n+u_lup,:),vet)
383            vet(:)=vet(:)/norm(vet)
384            ep_e(n+u_lup,:)=-vep(:)*ne(n,lup)
385            et_e(n+u_lup,:)=vet(:)*ne(n,lup)
386          ENDIF
387
388          CALL cross_product2(xyz_v(n+z_ldown,:),xyz_v(n+z_down,:),vep)
389          IF (norm(vep)>1e-30) THEN
390            vep(:)=vep(:)/norm(vep)
391            CALL cross_product2(vep,xyz_e(n+u_ldown,:),vet)
392            vet(:)=vet(:)/norm(vet)
393            ep_e(n+u_ldown,:)=-vep(:)*ne(n,ldown)
394            et_e(n+u_ldown,:)=vet(:)*ne(n,ldown)
395          ENDIF
396
397          CALL xyz2lonlat(xyz_i(n,:),lon,lat)
398          elon_i(n,1) = -sin(lon)
399          elon_i(n,2) = cos(lon)
400          elon_i(n,3) = 0
401          elat_i(n,1) = -cos(lon)*sin(lat)
402          elat_i(n,2) = -sin(lon)*sin(lat)
403          elat_i(n,3) = cos(lat)
404
405         
406          CALL xyz2lonlat(xyz_e(n+u_right,:),lon,lat)
407          elon_e(n+u_right,1) = -sin(lon)
408          elon_e(n+u_right,2) = cos(lon)
409          elon_e(n+u_right,3) = 0
410          elat_e(n+u_right,1) = -cos(lon)*sin(lat)
411          elat_e(n+u_right,2) = -sin(lon)*sin(lat)
412          elat_e(n+u_right,3) = cos(lat)
413         
414          CALL xyz2lonlat(xyz_e(n+u_lup,:),lon,lat)
415          elon_e(n+u_lup,1) = -sin(lon)
416          elon_e(n+u_lup,2) = cos(lon)
417          elon_e(n+u_lup,3) = 0
418          elat_e(n+u_lup,1) = -cos(lon)*sin(lat)
419          elat_e(n+u_lup,2) = -sin(lon)*sin(lat)
420          elat_e(n+u_lup,3) = cos(lat)
421 
422          CALL xyz2lonlat(xyz_e(n+u_ldown,:),lon,lat)
423          elon_e(n+u_ldown,1) = -sin(lon)
424          elon_e(n+u_ldown,2) = cos(lon)
425          elon_e(n+u_ldown,3) = 0
426          elat_e(n+u_ldown,1) = -cos(lon)*sin(lat)
427          elat_e(n+u_ldown,2) = -sin(lon)*sin(lat)
428          elat_e(n+u_ldown,3) = cos(lat)
429 
430         
431          DO k=0,5
432            CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(k+1),:),S1)
433            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)
434            Riv(n,k+1)=0.5*(S1+S2)/Ai(n)
435            Riv2(n,k+1)=0.5*(S1+S2)/surf_v(k+1)
436          ENDDO
437
438          DO k=1,6
439            IF (ABS(surf_v(k))<1e-30) THEN
440              Riv(n,k)=0.
441            ENDIF
442          ENDDO
443         
444          Av(n+z_up)=surf_v(vup)+1e-100
445          Av(n+z_down)=surf_v(vdown)+1e-100
446       
447        ENDDO
448      ENDDO
449     
450      DO j=jj_begin,jj_end
451        DO i=ii_begin,ii_end
452          n=(j-1)*iim+i
453   
454          CALL compute_wee(n,right,w)
455          Wee(n+u_right,:,1)=w(1:5)
456
457          CALL compute_wee(n+t_right,left,w)
458          Wee(n+u_right,:,2)=w(1:5)
459
460
461          CALL compute_wee(n,lup,w)
462          Wee(n+u_lup,:,1)=w(1:5)
463
464          CALL compute_wee(n+t_lup,rdown,w)
465          Wee(n+u_lup,:,2)=w(1:5)
466
467
468          CALL compute_wee(n,ldown,w)
469          Wee(n+u_ldown,:,1)=w(1:5)
470
471          CALL compute_wee(n+t_ldown,rup,w)
472          Wee(n+u_ldown,:,2)=w(1:5)
473
474        ENDDO
475      ENDDO
476     
477      DO j=jj_begin,jj_end
478        DO i=ii_begin,ii_end
479          n=(j-1)*iim+i
480          ii_glo=d%ii_begin_glo-d%ii_begin+i
481          jj_glo=d%jj_begin_glo-d%jj_begin+j
482         
483          IF (ii_glo==1 .AND. jj_glo==1) THEN
484            le(n+u_ldown)=0
485            xyz_v(n+z_ldown,:)=xyz_v(n+z_down,:)
486                       
487          ENDIF
488
489          IF (ii_glo==iim_glo .AND. jj_glo==1) THEN
490            le(n+u_right)=0
491            xyz_v(n+z_rdown,:)=xyz_v(n+z_rup,:)
492          ENDIF
493
494          IF (ii_glo==iim_glo .AND. jj_glo==jjm_glo) THEN
495            le(n+u_rup)=0
496            xyz_v(n+z_rup,:)=xyz_v(n+z_up,:)
497          ENDIF
498
499          IF (ii_glo==1 .AND. jj_glo==jjm_glo) THEN
500            le(n+u_lup)=0
501            xyz_v(n+z_up,:)=xyz_v(n+z_lup,:)
502          ENDIF
503         
504        ENDDO
505      ENDDO
506
507      DO j=jj_begin-1,jj_end+1
508        DO i=ii_begin-1,ii_end+1
509          n=(j-1)*iim+i
510          xyz_i(n,:)=xyz_i(n,:) * radius
511          xyz_v(n+z_up,:)=xyz_v(n+z_up,:) * radius
512          xyz_v(n+z_down,:)=xyz_v(n+z_down,:) *radius
513          de(n+u_right)=de(n+u_right) * radius
514          de(n+u_lup)=de(n+u_lup)*radius
515          de(n+u_ldown)=de(n+u_ldown)*radius
516          xyz_e(n+u_right,:)=xyz_e(n+u_right,:)*radius
517          xyz_e(n+u_lup,:)=xyz_e(n+u_lup,:)*radius
518          xyz_e(n+u_ldown,:)=xyz_e(n+u_ldown,:)*radius
519          le(n+u_right)=le(n+u_right)*radius
520          le(n+u_lup)=le(n+u_lup)*radius
521          le(n+u_ldown)=le(n+u_ldown)*radius
522          Ai(n)=Ai(n)*radius**2
523          Av(n+z_up)=Av(n+z_up)*radius**2
524          Av(n+z_down)=Av(n+z_down)*radius**2
525        ENDDO
526      ENDDO
527                 
528    ENDDO
529
530    CALL transfert_request(geom%Ai,req_i1)
531    CALL transfert_request(geom%centroid,req_i1)
532
533    CALL surf_triangle(d%xyz(:,ii_begin,jj_begin),d%xyz(:,ii_begin,jj_end),d%xyz(:,ii_end,jj_begin),S)
534 
535  END SUBROUTINE set_geometry
536 
537  SUBROUTINE compute_wee(n,pos,w)
538  IMPLICIT NONE
539    INTEGER,INTENT(IN) :: n
540    INTEGER,INTENT(IN) :: pos
541    REAL(rstd),INTENT(OUT) ::w(6)
542
543    REAL(rstd) :: ne_(0:5)
544    REAL(rstd) :: Riv_(6)
545    INTEGER    :: k
546   
547   
548    DO k=0,5
549      ne_(k)=ne(n,MOD(pos-1+k+6,6)+1) 
550      Riv_(k+1)=Riv(n,MOD(pos-1+k+6,6)+1)
551    ENDDO
552         
553    w(1)=-ne_(0)*ne_(1)*(Riv_(1)-0.5)
554    w(2)=-ne_(2)*(ne_(0)*Riv_(2)-w(1)*ne_(1))
555    w(3)=-ne_(3)*(ne_(0)*Riv_(3)-w(2)*ne_(2))
556    w(4)=-ne_(4)*(ne_(0)*Riv_(4)-w(3)*ne_(3))
557    w(5)=-ne_(5)*(ne_(0)*Riv_(5)-w(4)*ne_(4))
558    w(6)=ne_(0)*ne_(5)*(Riv_(6)-0.5)
559   
560!    IF ( ABS(w(5)-w(6))>1e-20) PRINT *, "pb pour wee : w(5)!=w(6)",sum(Riv_(:))
561
562   END SUBROUTINE compute_wee
563           
564
565 
566  SUBROUTINE compute_geometry
567  IMPLICIT NONE
568    CALL allocate_geometry
569    CALL set_geometry
570   
571  END SUBROUTINE compute_geometry
572 
573END MODULE geometry
Note: See TracBrowser for help on using the repository browser.