source: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM_OLD/src/geometry.f90 @ 329

Last change on this file since 329 was 318, checked in by ymipsl, 9 years ago

Less output duringg the grid optimization phase

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