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

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

devel : store cell bounds once, use them for XIOS later
File size: 16.6 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    DO ind=1,ndomain
126      IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master ) CYCLE
127      d=>domain(ind)
128      CALL swap_dimensions(ind)
129      CALL swap_geometry(ind)
130      DO j=jj_begin,jj_end
131        DO i=ii_begin,ii_end
132          n=(j-1)*iim+i
133          DO k=0,5
134            x1(:) = xyz_v(n+z_pos(k+1),:)
135            x2(:) = d%vertex(:,k,i,j) 
136            IF (norm(x1-x2)>1e-10) THEN
137              PRINT*,"vertex diff ",ind,i,j,k
138              PRINT*,x1
139              PRINT*,x2
140            ENDIF
141          ENDDO
142        ENDDO
143      ENDDO
144    ENDDO
145   
146   
147    DO it=1,nb_it
148      IF (MOD(it,100)==0) THEN
149        check=is_master
150      ELSE
151        check=.FALSE.
152      ENDIF
153     
154      sum=0
155      DO ind=1,ndomain
156      IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master ) CYCLE
157        CALL swap_dimensions(ind)
158        CALL swap_geometry(ind)
159        DO j=jj_begin,jj_end
160          DO i=ii_begin,ii_end
161            n=(j-1)*iim+i
162            vect(:,1)=xyz_v(n+z_rup,:)
163            vect(:,2)=xyz_v(n+z_up,:)
164            vect(:,3)=xyz_v(n+z_lup,:)
165            vect(:,4)=xyz_v(n+z_ldown,:)
166            vect(:,5)=xyz_v(n+z_down,:)
167            vect(:,6)=xyz_v(n+z_rdown,:)
168            CALL compute_centroid(vect,6,centr)
169            IF (check) THEN
170              sum=MAX(sum,norm(xyz_i(n,:)-centr(:)))
171            ENDIF
172            xyz_i(n,:)=centr(:)
173          ENDDO
174        ENDDO
175       
176      ENDDO
177
178      IF (check) PRINT *,"it = ",it,"  diff centroid circumcenter ",sum
179
180     CALL update_circumcenters
181
182    ENDDO
183   
184  END SUBROUTINE optimize_geometry
185
186  SUBROUTINE update_domain 
187    ! copy position of generators and vertices back into domain(:)%xyz/vertex
188    ! so that XIOS/create_header_gen gets the right values
189    USE transfert_mpi_mod, ONLY : gather_field, bcast_field
190
191    INTEGER :: ind,i,j,k,n
192    TYPE(t_domain),POINTER :: d
193    TYPE(t_field),POINTER,SAVE :: xyz_glo(:), xyz_loc(:), vertex_glo(:), vertex_loc(:)
194    REAL(rstd), POINTER :: xyz(:,:), vertex(:,:)
195   
196    CALL allocate_field(xyz_loc, field_t, type_real, 3)
197    CALL allocate_field(vertex_loc, field_z, type_real, 3)
198
199    DO ind=1,ndomain
200       IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master ) CYCLE
201       CALL swap_dimensions(ind)
202       CALL swap_geometry(ind)
203       xyz = xyz_loc(ind)
204       xyz(:,:) = xyz_i(:,:)
205       vertex = vertex_loc(ind)
206       vertex(:,:) = xyz_v(:,:)
207    END DO
208
209!$OMP BARRIER   
210!$OMP MASTER
211    CALL allocate_field_glo(xyz_glo, field_t, type_real, 3)
212    CALL allocate_field_glo(vertex_glo, field_z, type_real, 3)
213
214    CALL gather_field(xyz_loc, xyz_glo)
215    CALL gather_field(vertex_loc, vertex_glo)
216
217    CALL bcast_field(xyz_glo)
218    CALL bcast_field(vertex_glo)
219
220    DO ind=1,ndomain_glo
221       d=>domain_glo(ind)
222       xyz = xyz_glo(ind)
223       vertex = vertex_glo(ind)
224       DO j=d%jj_begin,d%jj_end
225          DO i=d%ii_begin,d%ii_end         
226             n=(j-1)*d%iim+i
227             d%xyz(:,i,j)=xyz(n,:)
228             DO k=0,5
229                d%vertex(:,k,i,j) = vertex(n+d%z_pos(k+1),:)
230             END DO
231          END DO
232       END DO
233    END DO
234
235    CALL deallocate_field_glo(vertex_glo)
236    CALL deallocate_field_glo(xyz_glo)
237!$OMP END MASTER
238!$OMP BARRIER
239    CALL deallocate_field(vertex_loc)
240    CALL deallocate_field(xyz_loc)
241
242  END SUBROUTINE update_domain
243 
244  SUBROUTINE set_geometry
245  USE metric
246
247    REAL(rstd) :: surf(6)
248    REAL(rstd) :: surf_v(6)
249    REAL(rstd) :: vect(3,6)
250    REAL(rstd) :: centr(3)
251    REAL(rstd) :: vet(3),vep(3), vertex(3)
252    INTEGER :: ind,i,j,k,n
253    TYPE(t_domain),POINTER :: d
254    REAL(rstd) ::  S12
255    REAL(rstd) :: w(6)
256    REAL(rstd) :: lon,lat
257    INTEGER :: ii_glo,jj_glo
258    REAL(rstd) :: S
259         
260     
261    CALL optimize_geometry
262    CALL remap_schmidt_loc
263    CALL update_circumcenters
264
265    DO ind=1,ndomain
266      IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master ) CYCLE
267      d=>domain(ind)
268      CALL swap_dimensions(ind)
269      CALL swap_geometry(ind)
270      lon_i(:)=0 ; lat_i(:)=0
271      lon_e(:)=0 ; lat_e(:)=0
272      DO j=jj_begin-1,jj_end+1
273        DO i=ii_begin-1,ii_end+1
274          n=(j-1)*iim+i
275
276          DO k=0,5
277            ne(n,k+1)=d%ne(k,i,j)
278          ENDDO
279
280          vect(:,1)=xyz_v(n+z_rup,:)
281          vect(:,2)=xyz_v(n+z_up,:)
282          vect(:,3)=xyz_v(n+z_lup,:)
283          vect(:,4)=xyz_v(n+z_ldown,:)
284          vect(:,5)=xyz_v(n+z_down,:)
285          vect(:,6)=xyz_v(n+z_rdown,:)
286          CALL compute_centroid(vect,6,centr)
287          centroid(n,:)=centr(:)
288
289               
290          CALL xyz2lonlat(xyz_v(n+z_up,:),lon,lat)
291          fv(n+z_up)=2*sin(lat)*omega
292          CALL xyz2lonlat(xyz_v(n+z_down,:),lon,lat)
293          fv(n+z_down)=2*sin(lat)*omega
294         
295          bi(n)=0. 
296       
297          CALL dist_cart(xyz_i(n,:),xyz_i(n+t_right,:),de(n+u_right))
298          CALL dist_cart(xyz_i(n,:),xyz_i(n+t_lup,:),de(n+u_lup))
299          CALL dist_cart(xyz_i(n,:),xyz_i(n+t_ldown,:),de(n+u_ldown))
300         
301          CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_right,:),0.5,xyz_e(n+u_right,:))
302          CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_lup,:),0.5,xyz_e(n+u_lup,:))
303          CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_ldown,:),0.5,xyz_e(n+u_ldown,:))
304
305          CALL dist_cart(xyz_v(n+z_rdown,:), xyz_v(n+z_rup,:),le(n+u_right))
306          CALL dist_cart(xyz_v(n+z_up,:), xyz_v(n+z_lup,:),le(n+u_lup))
307          CALL dist_cart(xyz_v(n+z_ldown,:), xyz_v(n+z_down,:),le(n+u_ldown))
308
309          le_de(n+u_right)=le(n+u_right)/de(n+u_right) ! NaN possible but should be harmless
310          le_de(n+u_lup)  =le(n+u_lup)  /de(n+u_lup)
311          le_de(n+u_ldown)=le(n+u_ldown)/de(n+u_ldown)
312
313          Ai(n)=0
314          DO k=0,5
315            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))
316            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))
317            Ai(n)=Ai(n)+surf(k+1)
318            IF (i==ii_end .AND. j==jj_begin) THEN
319              IF (Ai(n)<1e20) THEN
320              ELSE
321                PRINT *,"PB !!",Ai(n),k,surf(k+1)
322                PRINT*,xyz_i(n,:),xyz_v(n+z_pos(MOD((k-1+6),6)+1),:),xyz_v(n+z_pos(k+1),:)
323              ENDIF
324            ENDIF
325          ENDDO
326
327          ! Sign convention : Ringler et al., JCP 2010, eq. 21 p. 3071
328          ! Normal component is along outgoing normal vector if ne=1
329
330          CALL cross_product2(xyz_v(n+z_rdown,:),xyz_v(n+z_rup,:),vep)
331          IF (norm(vep)>1e-30) THEN
332            vep(:)=vep(:)/norm(vep)                         ! Inward normal vector
333            CALL cross_product2(vep,xyz_e(n+u_right,:),vet) ! Counter-clockwise tangent vector
334            vet(:)=vet(:)/norm(vet)
335            ep_e(n+u_right,:)=-vep(:)*ne(n,right)
336            et_e(n+u_right,:)=vet(:)*ne(n,right)
337          ENDIF
338
339          CALL cross_product2(xyz_v(n+z_up,:),xyz_v(n+z_lup,:),vep)
340          IF (norm(vep)>1e-30) THEN
341            vep(:)=vep(:)/norm(vep)
342            CALL cross_product2(vep,xyz_e(n+u_lup,:),vet)
343            vet(:)=vet(:)/norm(vet)
344            ep_e(n+u_lup,:)=-vep(:)*ne(n,lup)
345            et_e(n+u_lup,:)=vet(:)*ne(n,lup)
346          ENDIF
347
348          CALL cross_product2(xyz_v(n+z_ldown,:),xyz_v(n+z_down,:),vep)
349          IF (norm(vep)>1e-30) THEN
350            vep(:)=vep(:)/norm(vep)
351            CALL cross_product2(vep,xyz_e(n+u_ldown,:),vet)
352            vet(:)=vet(:)/norm(vet)
353            ep_e(n+u_ldown,:)=-vep(:)*ne(n,ldown)
354            et_e(n+u_ldown,:)=vet(:)*ne(n,ldown)
355          ENDIF
356
357          CALL xyz2lonlat(xyz_i(n,:),lon,lat)
358          lon_i(n)=lon
359          lat_i(n)=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          lon_e(n+u_right)=lon
370          lat_e(n+u_right)=lat
371          elon_e(n+u_right,1) = -sin(lon)
372          elon_e(n+u_right,2) = cos(lon)
373          elon_e(n+u_right,3) = 0
374          elat_e(n+u_right,1) = -cos(lon)*sin(lat)
375          elat_e(n+u_right,2) = -sin(lon)*sin(lat)
376          elat_e(n+u_right,3) = cos(lat)
377         
378          CALL xyz2lonlat(xyz_e(n+u_lup,:),lon,lat)
379          lon_e(n+u_lup)=lon
380          lat_e(n+u_lup)=lat
381          elon_e(n+u_lup,1) = -sin(lon)
382          elon_e(n+u_lup,2) = cos(lon)
383          elon_e(n+u_lup,3) = 0
384          elat_e(n+u_lup,1) = -cos(lon)*sin(lat)
385          elat_e(n+u_lup,2) = -sin(lon)*sin(lat)
386          elat_e(n+u_lup,3) = cos(lat)
387 
388          CALL xyz2lonlat(xyz_e(n+u_ldown,:),lon,lat)
389          lon_e(n+u_ldown)=lon
390          lat_e(n+u_ldown)=lat
391          elon_e(n+u_ldown,1) = -sin(lon)
392          elon_e(n+u_ldown,2) = cos(lon)
393          elon_e(n+u_ldown,3) = 0
394          elat_e(n+u_ldown,1) = -cos(lon)*sin(lat)
395          elat_e(n+u_ldown,2) = -sin(lon)*sin(lat)
396          elat_e(n+u_ldown,3) = cos(lat)
397 
398         
399          DO k=0,5
400            CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(k+1),:),S1(n,k+1) )
401            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) )
402            S12 = .5*(S1(n,k+1)+S2(n,k+1))
403            Riv(n,k+1)=S12/Ai(n)
404            Riv2(n,k+1)=S12/surf_v(k+1)
405          ENDDO
406
407          DO k=1,6
408            IF (ABS(surf_v(k))<1e-30) THEN
409              Riv(n,k)=0.
410            ENDIF
411          ENDDO
412         
413          Av(n+z_up)=surf_v(vup)+1e-100
414          Av(n+z_down)=surf_v(vdown)+1e-100
415       
416        ENDDO
417      ENDDO
418     
419      DO j=jj_begin,jj_end
420        DO i=ii_begin,ii_end
421          n=(j-1)*iim+i
422   
423          CALL compute_wee(n,right,w)
424          Wee(n+u_right,:,1)=w(1:5)
425
426          CALL compute_wee(n+t_right,left,w)
427          Wee(n+u_right,:,2)=w(1:5)
428
429
430          CALL compute_wee(n,lup,w)
431          Wee(n+u_lup,:,1)=w(1:5)
432
433          CALL compute_wee(n+t_lup,rdown,w)
434          Wee(n+u_lup,:,2)=w(1:5)
435
436
437          CALL compute_wee(n,ldown,w)
438          Wee(n+u_ldown,:,1)=w(1:5)
439
440          CALL compute_wee(n+t_ldown,rup,w)
441          Wee(n+u_ldown,:,2)=w(1:5)
442
443        ENDDO
444      ENDDO
445     
446      DO j=jj_begin,jj_end
447        DO i=ii_begin,ii_end
448          n=(j-1)*iim+i
449          ii_glo=d%ii_begin_glo-d%ii_begin+i
450          jj_glo=d%jj_begin_glo-d%jj_begin+j
451         
452          IF (ii_glo==1 .AND. jj_glo==1) THEN
453            le(n+u_ldown)=0
454            le_de(n+u_ldown)=0
455            xyz_v(n+z_ldown,:)=xyz_v(n+z_down,:)
456                       
457          ENDIF
458
459          IF (ii_glo==iim_glo .AND. jj_glo==1) THEN
460            le(n+u_right)=0
461            le_de(n+u_right)=0
462            xyz_v(n+z_rdown,:)=xyz_v(n+z_rup,:)
463          ENDIF
464
465          IF (ii_glo==iim_glo .AND. jj_glo==jjm_glo) THEN
466            le(n+u_rup)=0
467            le_de(n+u_rup)=0
468            xyz_v(n+z_rup,:)=xyz_v(n+z_up,:)
469          ENDIF
470
471          IF (ii_glo==1 .AND. jj_glo==jjm_glo) THEN
472            le(n+u_lup)=0
473            le_de(n+u_lup)=0
474            xyz_v(n+z_up,:)=xyz_v(n+z_lup,:)
475          ENDIF
476         
477        ENDDO
478      ENDDO
479
480      DO j=jj_begin-1,jj_end+1
481        DO i=ii_begin-1,ii_end+1
482          n=(j-1)*iim+i
483          xyz_i(n,:)=xyz_i(n,:) * radius
484          xyz_v(n+z_up,:)=xyz_v(n+z_up,:) * radius
485          xyz_v(n+z_down,:)=xyz_v(n+z_down,:) *radius
486          de(n+u_right)=de(n+u_right) * radius
487          de(n+u_lup)=de(n+u_lup)*radius
488          de(n+u_ldown)=de(n+u_ldown)*radius
489          xyz_e(n+u_right,:)=xyz_e(n+u_right,:)*radius
490          xyz_e(n+u_lup,:)=xyz_e(n+u_lup,:)*radius
491          xyz_e(n+u_ldown,:)=xyz_e(n+u_ldown,:)*radius
492          le(n+u_right)=le(n+u_right)*radius
493          le(n+u_lup)=le(n+u_lup)*radius
494          le(n+u_ldown)=le(n+u_ldown)*radius
495          Ai(n)=Ai(n)*radius**2
496          Av(n+z_up)=Av(n+z_up)*radius**2
497          Av(n+z_down)=Av(n+z_down)*radius**2
498        ENDDO
499      ENDDO
500                 
501    ENDDO
502
503    CALL transfert_request(geom%Ai,req_i1)
504    CALL transfert_request(geom%centroid,req_i1)
505
506!    CALL surf_triangle(d%xyz(:,ii_begin,jj_begin),d%xyz(:,ii_begin,jj_end),d%xyz(:,ii_end,jj_begin),S)
507 
508  END SUBROUTINE set_geometry
509
510  SUBROUTINE compute_wee(n,pos,w)
511    INTEGER,INTENT(IN) :: n
512    INTEGER,INTENT(IN) :: pos
513    REAL(rstd),INTENT(OUT) ::w(6)
514
515    REAL(rstd) :: ne_(0:5)
516    REAL(rstd) :: Riv_(6)
517    INTEGER    :: k
518   
519   
520    DO k=0,5
521      ne_(k)=ne(n,MOD(pos-1+k+6,6)+1) 
522      Riv_(k+1)=Riv(n,MOD(pos-1+k+6,6)+1)
523    ENDDO
524         
525    w(1)=-ne_(0)*ne_(1)*(Riv_(1)-0.5)
526    w(2)=-ne_(2)*(ne_(0)*Riv_(2)-w(1)*ne_(1))
527    w(3)=-ne_(3)*(ne_(0)*Riv_(3)-w(2)*ne_(2))
528    w(4)=-ne_(4)*(ne_(0)*Riv_(4)-w(3)*ne_(3))
529    w(5)=-ne_(5)*(ne_(0)*Riv_(5)-w(4)*ne_(4))
530    w(6)=ne_(0)*ne_(5)*(Riv_(6)-0.5)
531   
532!    IF ( ABS(w(5)-w(6))>1e-20) PRINT *, "pb pour wee : w(5)!=w(6)",sum(Riv_(:))
533
534   END SUBROUTINE compute_wee
535           
536 
537  SUBROUTINE compute_geometry
538    USE grid_param
539    USE init_unstructured_mod, ONLY : read_local_mesh
540    USE set_bounds_mod, ONLY : set_bounds
541    CALL allocate_geometry
542
543    SELECT CASE(grid_type)
544    CASE(grid_ico)
545       CALL set_geometry
546       ! copy position of generators and vertices back into domain_glo(:)%xyz/vertex
547       ! so that write_field gets the right values
548       CALL update_domain
549       CALL set_bounds(domain_glo, .TRUE.)
550       CALL set_bounds(domain, .FALSE.)
551    CASE(grid_unst)
552       swap_needed=.FALSE.
553       CALL read_local_mesh
554    CASE DEFAULT
555       STOP 'Invalid value of grid_type encountered in compute_geometry'
556    END SELECT
557  END SUBROUTINE compute_geometry
558 
559END MODULE compute_geometry_mod
Note: See TracBrowser for help on using the repository browser.