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

Last change on this file since 155 was 155, checked in by dubos, 11 years ago

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