source: codes/icosagcm/devel/src/sphere/compute_geometry.f90 @ 863

Last change on this file since 863 was 863, checked in by dubos, 5 years ago

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