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

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

Merging OpenMP parallisme mode : by subdomain and on vertical level.
This feature is actually experimental but may be retro-compatible with the last method based only on subdomain

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