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

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

Implementation of mixte parallelism MPI/OpenMP into src directory

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