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

Last change on this file since 294 was 294, checked in by dubos, 10 years ago

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