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

Last change on this file since 146 was 146, checked in by ymipsl, 11 years ago

Set constant sign for wind way :
ne(ij,right)==ne_right=1
ne(ij,rup)==ne_rup=-1
ne(ij,lup)==ne_lup=1
ne(ij,left)==ne_left=-1
ne(ij,ldown)==ne_ldown=1
ne(ij,rdown)==ne_rdown=-1

Modified transfert function to be compliant for this convention.

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