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

Last change on this file since 1008 was 993, checked in by rpennel, 5 years ago

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