Line data Source code
1 : ! Code that computes control volumes with the same area as the GLL weights
2 : ! (for SCRIP) uses analytic area formula.
3 :
4 : module comp_gll_ctr_vol
5 : use shr_kind_mod, only: r8=>shr_kind_r8, shr_kind_cl
6 : use cam_abortutils, only: endrun
7 : use cam_logfile, only: iulog
8 : use shr_sys_mod, only: shr_sys_flush
9 : use global_norms_mod, only: wrap_repro_sum
10 : use physconst, only: pi
11 : use infnan, only: isnan
12 :
13 : use coordinate_systems_mod, only: cartesian3d_t, cartesian2d_t
14 : use coordinate_systems_mod, only: spherical_polar_t, change_coordinates
15 : use coordinate_systems_mod, only: cubedsphere2cart, cube_face_number_from_cart
16 : use coordinate_systems_mod, only: distance, sphere_tri_area
17 : use parallel_mod, only: global_shared_buf, global_shared_sum
18 : use edgetype_mod, only: EdgeBuffer_t, Ghostbuffer3D_t
19 : use dimensions_mod, only: np, ne
20 : use control_mod, only: fine_ne
21 : use reduction_mod, only: red_sum, parallelmin, parallelmax
22 :
23 : implicit none
24 : private
25 : save
26 :
27 : character(len=16), public :: se_write_gll_grid = "no"
28 :
29 : ! nv_max will be set to 2*max_elements_attached_to_node
30 : ! This works out to 6 for the regular case and 14 for refined meshes
31 : integer :: nv_max=-99
32 :
33 : type, public :: ctrlvol_t
34 : real(r8) :: vol(np,np) ! area of the unit sphere covered (local)
35 : real(r8) :: totvol(np,np) ! area of the unit sphere covered (local)
36 : real(r8) :: invvol(np,np) ! inverse area (includes neigbors)
37 : type(cartesian2d_t) :: cartp_dual(0:np,0:np)
38 : type(cartesian3d_t) :: cart3d_dual(0:np,0:np)
39 : type(cartesian3D_t), allocatable :: vert(:,:,:) ! bounding box for the polygon
40 : type(spherical_polar_t), allocatable :: vert_latlon(:,:,:) ! bounding box for the polygon
41 : integer, allocatable :: face_no(:,:,:) ! face_no of cv vertex. 0 if on cube edge
42 : integer :: nvert(np,np) ! number of vertex per polygon
43 : end type ctrlvol_t
44 :
45 : ! tho options:
46 : ! (1) for NP<>4 or Refined Meshes (this is less accurate)
47 : ! build control volumes out of lines which are
48 : ! always gnomonic coordinate lines. results in hexagon control volumes
49 : ! at cube corners and edges. control volumes at cube-sphere edges are
50 : ! non-convex, which breaks SCRIP.
51 : ! iterative option for NP=4 only:
52 : ! (2) (USE_PENTAGONS)
53 : ! iterate to minimize difference between spherical area and GLL weight
54 : ! introduce pentagons in the center of each element to make areas agree
55 : ! control volumes are triangles, squares or pentagons
56 :
57 : type(ctrlvol_t), allocatable, target :: cvlist(:)
58 : type(EdgeBuffer_t) :: edge1
59 : type(GhostBuffer3D_t) :: ghost_buf
60 :
61 : ! User interface
62 : public :: gll_grid_write ! Write the grid in SCRIP format and exit
63 :
64 : ! Private interfaces
65 : private:: InitControlVolumesData ! allocates internal data structure
66 : private:: InitControlVolumes ! Inits all surfaces: vol,totvol, invvol
67 :
68 : private:: GetVolumeLocal
69 : private:: GetVolume
70 :
71 : logical, private :: initialized = .false.
72 :
73 : CONTAINS
74 :
75 0 : subroutine gll_grid_write(elem, grid_format, filename_in)
76 : use netcdf, only: nf90_strerror
77 : use spmd_utils, only: masterproc, mpicom
78 : use pio, only: var_desc_t, file_desc_t
79 : use pio, only: pio_int, pio_double, PIO_NOERR
80 : use pio, only: pio_put_att, pio_put_var, pio_enddef
81 : use cam_pio_utils, only: cam_pio_createfile, cam_pio_closefile
82 : use cam_pio_utils, only: cam_pio_def_dim, cam_pio_def_var
83 : use cam_grid_support, only: cam_grid_id, cam_grid_dimensions
84 : use cam_grid_support, only: cam_grid_write_dist_array
85 : !!XXgoldyXX: v debug only
86 : #define USE_PIO3D
87 : #ifdef USE_PIO3D
88 : use pio, only: io_desc_t, pio_write_darray, PIO_OFFSET_KIND
89 : use cam_pio_utils, only: cam_pio_newdecomp
90 : use spmd_utils, only: iam
91 : #endif
92 : !!XXgoldyXX: ^ debug only
93 :
94 : use hybrid_mod, only: hybrid_t, config_thread_region
95 : use parallel_mod, only: par
96 : use dimensions_mod, only: nelem, nelemd
97 : use control_mod, only: refined_mesh, fine_ne
98 : use element_mod, only: element_t
99 : use dof_mod, only: UniquePoints
100 : use coordinate_systems_mod, only: cart2spherical
101 :
102 : ! Inputs
103 : type(element_t), intent(in) :: elem(:)
104 : character(len=*), intent(in) :: grid_format
105 : character(len=*), intent(in) :: filename_in
106 :
107 : real(r8), parameter :: rad2deg = 180._r8/pi
108 :
109 : ! Local variables
110 : !!XXgoldyXX: v debug only
111 : #ifdef USE_PIO3D
112 0 : integer(PIO_OFFSET_KIND), allocatable :: ldof(:)
113 : integer :: ii, jj
114 : type(io_desc_t), pointer :: iodesc
115 : #endif
116 : !!XXgoldyXX: ^ debug only
117 : integer :: i, j, ie, ierror, status, ivtx, index
118 : integer :: grid_corners_id, grid_rank_id, grid_size_id
119 : type(var_desc_t) :: grid_dims_id, grid_area_id, grid_center_lat_id
120 : type(var_desc_t) :: grid_center_lon_id, grid_corner_lat_id
121 : type(var_desc_t) :: grid_corner_lon_id, grid_imask_id
122 :
123 : type(file_desc_t) :: file
124 : integer :: gll_grid ! Grid ID
125 : integer :: gridsize ! Total number of unique grid columns
126 : integer :: arr_dims2d(2) ! (/ np*np, nelemed)
127 : integer :: file_dims2d(1) ! (/ gridsize /)
128 : integer :: arr_dims3d(3) ! (/ np*np, nv_max, nelemed)
129 : integer :: file_dims3d(2) ! (/ nv_max, gridsize /)
130 :
131 0 : real(r8), allocatable :: gwork(:,:,:) ! np*np, nv_max, nelemd
132 : type(hybrid_t) :: hybrid
133 : character(len=256) :: errormsg
134 : character(len=shr_kind_cl) :: filename
135 : type(spherical_polar_t) :: sphere
136 : character(len=*), parameter :: subname = 'gll_grid_write'
137 :
138 : !! Check to see if we are doing grid output
139 0 : if (trim(grid_format) == "no") then
140 0 : if (masterproc) then
141 0 : write(iulog, *) subname, ': Not writing phys_grid file.'
142 : end if
143 0 : return
144 0 : else if (trim(grid_format) /= 'SCRIP') then
145 0 : write(errormsg, *) subname, ': ERROR, bad value for se_write_grid, ', &
146 0 : trim(grid_format)
147 0 : call endrun(errormsg)
148 : end if
149 :
150 : ! Set up the control volumes
151 0 : if (refined_mesh) then
152 0 : nv_max = 14
153 : else
154 0 : nv_max = 5
155 : end if
156 0 : if (masterproc) then
157 0 : write(iulog, *) subname, ': computing GLL dual grid for control volumes:'
158 : end if
159 0 : call InitControlVolumesData(par,elem,nelemd)
160 : ! single thread
161 0 : hybrid = config_thread_region(par,'serial')
162 0 : call InitControlVolumes(elem,hybrid,1,nelemd)
163 0 : if (masterproc) then
164 0 : write(6, *) subname, ': done computing GLL dual grid for control volumes.'
165 : end if
166 :
167 : ! Create the NetCDF file
168 0 : if (len_trim(filename_in) == 0) then
169 0 : if (refined_mesh) then
170 0 : if (fine_ne <= 0) then
171 0 : call endrun('gll_grid_write: refined_mesh selected but fine_ne not set')
172 : else
173 0 : write(filename,'(a,i0,a,a,3a)') "ne0np", np, "_refined_", trim(grid_format), ".nc"
174 : end if
175 : else
176 0 : write(filename, '(a,i0,a,i0,a,a,3a)') "ne", ne, "np", np, &
177 0 : "_", trim(grid_format), ".nc"
178 : end if
179 : else
180 0 : filename = trim(filename_in)
181 : end if
182 0 : if (masterproc) then
183 0 : write(iulog, *) 'Writing gll SCRIP grid file: ', trim(filename)
184 0 : call shr_sys_flush(iulog)
185 : end if
186 :
187 0 : call cam_pio_createfile(file, trim(filename))
188 0 : gll_grid = cam_grid_id('GLL')
189 0 : call cam_grid_dimensions(gll_grid, file_dims3d)
190 0 : gridsize = file_dims3d(1)
191 0 : file_dims2d(1) = gridsize
192 0 : file_dims3d(1) = nv_max
193 0 : file_dims3d(2) = gridsize
194 0 : arr_dims2d(1) = np*np
195 0 : arr_dims2d(2) = nelemd
196 0 : arr_dims3d(1) = np*np
197 0 : arr_dims3d(2) = nv_max
198 0 : arr_dims3d(3) = nelemd
199 0 : call cam_pio_def_dim(file, "grid_corners", nv_max, grid_corners_id)
200 0 : call cam_pio_def_dim(file, "grid_rank", 1, grid_rank_id)
201 0 : call cam_pio_def_dim(file, "grid_size", gridsize, grid_size_id)
202 : ! Define the coordinate variables
203 : call cam_pio_def_var(file, "grid_dims", pio_int, (/ grid_rank_id /), &
204 0 : grid_dims_id)
205 :
206 : ! Define grid area
207 : call cam_pio_def_var(file, "grid_area", pio_double, &
208 0 : (/grid_size_id/), grid_area_id)
209 0 : status = pio_put_att(file, grid_area_id, "units", "radians^2")
210 0 : if (status /= pio_noerr) then
211 0 : write(iulog, *) subname,': Error defining units attribute for grid_area'
212 0 : call shr_sys_flush(iulog)
213 0 : call endrun(subname//": "//trim(nf90_strerror(status)))
214 : end if
215 0 : status = pio_put_att(file, grid_area_id, "long_name", "area weights")
216 0 : if (status /= pio_noerr) then
217 0 : write(iulog, *) subname,': Error defining long_name attribute for grid_area'
218 0 : call shr_sys_flush(iulog)
219 0 : call endrun(subname//": "//trim(nf90_strerror(status)))
220 : end if
221 :
222 : ! Define grid center latitudes
223 : call cam_pio_def_var(file, "grid_center_lat", pio_double, &
224 0 : (/grid_size_id/), grid_center_lat_id)
225 0 : status = pio_put_att(file, grid_center_lat_id, "units", "degrees")
226 0 : if (status /= pio_noerr) then
227 0 : write(iulog, *) subname,': Error defining units attribute for grid_center_lat'
228 0 : call shr_sys_flush(iulog)
229 0 : call endrun(subname//": "//trim(nf90_strerror(status)))
230 : end if
231 :
232 : ! Define grid center longitudes
233 : call cam_pio_def_var(file, "grid_center_lon", pio_double, &
234 0 : (/grid_size_id/), grid_center_lon_id)
235 0 : status = pio_put_att(file, grid_center_lon_id, "units", "degrees")
236 0 : if (status /= pio_noerr) then
237 0 : write(iulog, *) subname,': Error defining units attribute for grid_center_lon'
238 0 : call shr_sys_flush(iulog)
239 0 : call endrun(subname//": "//trim(nf90_strerror(status)))
240 : end if
241 :
242 : ! Define grid corner latitudes
243 : call cam_pio_def_var(file, "grid_corner_lat", pio_double, &
244 0 : (/grid_corners_id, grid_size_id/), grid_corner_lat_id)
245 0 : status = pio_put_att(file, grid_corner_lat_id, "units", "degrees")
246 0 : if (status /= pio_noerr) then
247 0 : write(iulog, *) subname,': Error defining units attribute for grid_corner_lon'
248 0 : call shr_sys_flush(iulog)
249 0 : call endrun(subname//": "//trim(nf90_strerror(status)))
250 : end if
251 :
252 : ! Grid corner longitudes
253 : call cam_pio_def_var(file, "grid_corner_lon", pio_double, &
254 0 : (/grid_corners_id, grid_size_id/), grid_corner_lon_id)
255 0 : status = pio_put_att(file, grid_corner_lon_id, "units", "degrees")
256 0 : if (status /= pio_noerr) then
257 0 : write(iulog, *) subname,': Error defining units attribute for grid_corner_lon'
258 0 : call shr_sys_flush(iulog)
259 0 : call endrun(subname//": "//trim(nf90_strerror(status)))
260 : end if
261 :
262 : ! Grid mask
263 : call cam_pio_def_var(file, "grid_imask", pio_double, &
264 0 : (/grid_size_id/), grid_imask_id)
265 :
266 : ! End of NetCDF definitions
267 0 : status = PIO_enddef(file)
268 0 : if (status /= pio_noerr) then
269 0 : write(iulog, *) subname, ': Error calling enddef'
270 0 : call shr_sys_flush(iulog)
271 0 : call endrun(subname//": "//trim(nf90_strerror(status)))
272 : end if
273 :
274 : ! Work array to gather info before writing
275 0 : allocate(gwork(np*np, nv_max, nelemd))
276 :
277 : ! Write grid size
278 0 : status = pio_put_var(file, grid_dims_id, (/ gridsize /))
279 0 : if (status /= pio_noerr) then
280 0 : write(iulog, *) subname, ': Error writing variable grid_dims'
281 0 : call shr_sys_flush(iulog)
282 0 : call endrun(subname//": "//trim(nf90_strerror(status)))
283 : end if
284 : ! Write GLL grid areas
285 0 : do ie = 1, nelemd
286 : index = 1
287 0 : do j = 1, np
288 0 : do i = 1, np
289 0 : gwork(index, 1, ie) = cvlist(ie)%vol(i,j)
290 0 : index = index + 1
291 : end do
292 : end do
293 : end do
294 : call cam_grid_write_dist_array(file, gll_grid, arr_dims2d, file_dims2d, &
295 0 : gwork(:,1,:), grid_area_id)
296 : ! Write GLL grid cell center latitude
297 0 : do ie = 1, nelemd
298 : index = 1
299 0 : do j = 1, np
300 0 : do i = 1, np
301 0 : gwork(index, 1, ie) = elem(ie)%spherep(i,j)%lat * rad2deg
302 0 : index = index + 1
303 : end do
304 : end do
305 : end do
306 : call cam_grid_write_dist_array(file, gll_grid, arr_dims2d, file_dims2d, &
307 0 : gwork(:,1,:), grid_center_lat_id)
308 : ! Write GLL grid cell center longitude
309 0 : do ie = 1, nelemd
310 : index = 1
311 0 : do j = 1, np
312 0 : do i = 1, np
313 0 : gwork(index, 1, ie) = elem(ie)%spherep(i,j)%lon * rad2deg
314 0 : index = index + 1
315 : end do
316 : end do
317 : end do
318 : call cam_grid_write_dist_array(file, gll_grid, arr_dims2d, file_dims2d, &
319 0 : gwork(:,1,:), grid_center_lon_id)
320 :
321 : ! GLL grid corners
322 : ! Collect all information for the grid corner latitude (counter-clockwise)
323 0 : do ie = 1, nelemd
324 0 : do ivtx = 1, nv_max
325 : index = 1
326 0 : do j = 1, np
327 0 : do i = 1, np
328 0 : gwork(index, ivtx, ie) = cvlist(ie)%vert_latlon(ivtx,i,j)%lat * rad2deg
329 0 : index = index + 1
330 : end do
331 : end do
332 : end do
333 : end do
334 : !!XXgoldyXX: v debug only
335 : #ifdef USE_PIO3D
336 0 : allocate(ldof(np*np*nelemd*nv_max))
337 0 : ldof = 0
338 0 : do ie = 1, nelemd
339 0 : do index = 1, elem(ie)%idxP%NumUniquePts
340 0 : i = elem(ie)%idxP%ia(index)
341 0 : j = elem(ie)%idxP%ja(index)
342 0 : ii = (i - 1) + ((j - 1) * np) + ((ie - 1) * np * np * nv_max) + 1
343 0 : jj = (elem(ie)%idxP%UniquePtOffset + index - 2) * nv_max
344 0 : do ivtx = 1, nv_max
345 0 : ldof(ii) = jj + ivtx
346 0 : if ((jj+ivtx < 1) .or. (jj+ivtx > gridsize*nv_max)) then
347 0 : write(errormsg, '(4(a,i0))') ' ERROR (',iam,'): ldof(',ii,') = ',jj + ivtx,' > ',gridsize*nv_max
348 0 : call endrun(subname//trim(errormsg))
349 : end if
350 0 : ii = ii + np*np
351 : end do
352 : end do
353 : end do
354 0 : allocate(iodesc)
355 0 : call cam_pio_newdecomp(iodesc, (/ nv_max, gridsize /), ldof, PIO_double)
356 0 : call pio_write_darray(file, grid_corner_lat_id, iodesc, gwork, status)
357 : #else
358 : !!XXgoldyXX: ^ debug only
359 : call cam_grid_write_dist_array(file, gll_grid, arr_dims3d, file_dims3d, &
360 : gwork, grid_corner_lat_id)
361 : !!XXgoldyXX: v debug only
362 : #endif
363 : !!XXgoldyXX: ^ debug only
364 : ! Collect all information for the grid corner longitude (counter-clockwise)
365 0 : do ie = 1, nelemd
366 0 : do ivtx = 1, nv_max
367 : index = 1
368 0 : do j = 1, np
369 0 : do i = 1, np
370 0 : gwork(index, ivtx, ie) = cvlist(ie)%vert_latlon(ivtx,i,j)%lon * rad2deg
371 0 : index = index + 1
372 : end do
373 : end do
374 : end do
375 : end do
376 : !!XXgoldyXX: v debug only
377 : #ifdef USE_PIO3D
378 0 : call pio_write_darray(file, grid_corner_lon_id, iodesc, gwork, status)
379 : #else
380 : !!XXgoldyXX: ^ debug only
381 : call cam_grid_write_dist_array(file, gll_grid, arr_dims3d, file_dims3d, &
382 : gwork, grid_corner_lon_id)
383 : !!XXgoldyXX: v debug only
384 : #endif
385 : !!XXgoldyXX: ^ debug only
386 : ! Grid imask
387 0 : gwork(:,1,:) = 1.0_r8
388 : call cam_grid_write_dist_array(file, gll_grid, arr_dims2d, file_dims2d, &
389 0 : gwork(:,1,:), grid_imask_id)
390 :
391 0 : call mpi_barrier(mpicom, ierror)
392 : ! Close the file
393 0 : call cam_pio_closefile(file)
394 0 : if(masterproc) then
395 0 : write(iulog, *) 'Finished writing physics grid file: ', trim(filename)
396 0 : call shr_sys_flush(iulog)
397 : end if
398 :
399 0 : end subroutine gll_grid_write
400 :
401 : ! elemid is the local element id (in nets:nete)
402 0 : function GetVolume(elemid) result(vol)
403 :
404 : integer, intent(in) :: elemid
405 : real(kind=r8), pointer :: vol(:,:)
406 :
407 0 : if(.not. initialized) then
408 0 : call endrun('Attempt to use volumes prior to initializing')
409 : end if
410 0 : vol => cvlist(elemid)%totvol
411 :
412 0 : end function GetVolume
413 :
414 0 : function GetVolumeLocal(elemid) result(vol)
415 :
416 : integer, intent(in) :: elemid
417 : real(r8), pointer :: vol(:,:)
418 :
419 0 : if(.not. initialized) then
420 0 : call endrun('Attempt to use volumes prior to initializing')
421 : end if
422 0 : vol => cvlist(elemid)%vol
423 :
424 0 : end function GetVolumeLocal
425 :
426 0 : subroutine InitControlVolumesData(par, elem, nelemd)
427 : use edge_mod, only: initedgebuffer, initGhostBuffer3D
428 : use parallel_mod, only: parallel_t, HME_BNDRY_P2P
429 : use element_mod, only: element_t
430 : use thread_mod, only: horz_num_threads
431 :
432 : type(parallel_t), intent(in) :: par
433 : type(element_t), intent(in) :: elem(:)
434 : integer, intent(in) :: nelemd
435 :
436 : integer :: ie
437 :
438 : ! Cannot be done in a threaded region
439 0 : allocate(cvlist(nelemd))
440 0 : do ie = 1, nelemd
441 0 : allocate(cvlist(ie)%vert(nv_max, np,np))
442 0 : allocate(cvlist(ie)%vert_latlon(nv_max,np,np))
443 0 : allocate(cvlist(ie)%face_no(nv_max,np,np))
444 : end do
445 :
446 0 : call initedgebuffer(par,edge1,elem,3,bndry_type=HME_BNDRY_P2P, nthreads=1)
447 0 : call initGhostBuffer3D(ghost_buf,3,np+1,1)
448 0 : end subroutine InitControlVolumesData
449 :
450 0 : subroutine VerifyAreas(elem,hybrid,nets,nete)
451 :
452 0 : use element_mod, only: element_t
453 : use hybrid_mod, only: hybrid_t
454 :
455 : integer, intent(in) :: nets,nete
456 : type(element_t), intent(in), target :: elem(:)
457 : type(hybrid_t), intent(in) :: hybrid
458 :
459 : integer :: i, j, ie, k, kptr, kmax
460 : real(r8) :: rspheremp(np,np)
461 : real(r8) :: invvol(np,np)
462 : real(r8) :: error, max_error, max_invvol, maxrsphere
463 :
464 0 : error = 0
465 0 : max_error = 0
466 0 : do ie=nets,nete
467 0 : rspheremp = elem(ie)%rspheremp
468 0 : invvol = cvlist(ie)%invvol
469 0 : do j=1,np
470 0 : do i=1,np
471 0 : error = 100*ABS(rspheremp(i,j)-invvol(i,j))/invvol(i,j)
472 0 : if (max_error.lt.error) then
473 0 : max_error = error
474 : max_invvol = invvol(i,j)
475 : maxrsphere = rspheremp(i,j)
476 : end if
477 : end do
478 : end do
479 : end do
480 0 : print '(A,F16.4 )',"Control Volume Stats: Max error percent:", max_error
481 0 : print '(A,F16.12)'," Value From Element:",1/maxrsphere
482 0 : print '(A,F16.12)'," Value From Control Volume:",1/max_invvol
483 0 : max_error = parallelmax(max_error,hybrid)
484 0 : if (hybrid%masterthread) then
485 0 : write(6, '(a,f16.4)') "Control volume area vs. gll area: max error (percent):", max_error
486 : end if
487 :
488 0 : end subroutine VerifyAreas
489 :
490 :
491 0 : subroutine InitControlVolumes(elem, hybrid,nets,nete)
492 : use element_mod, only: element_t
493 : use hybrid_mod, only: hybrid_t
494 : use control_mod, only: refined_mesh
495 :
496 : integer, intent(in) :: nets,nete
497 : type(element_t), intent(in), target :: elem(:)
498 : type(hybrid_t), intent(in) :: hybrid
499 :
500 0 : if (refined_mesh .or. (np /= 4)) then
501 0 : call InitControlVolumes_duel(elem, hybrid,nets,nete)
502 : else
503 0 : call InitControlVolumes_gll(elem, hybrid,nets,nete)
504 0 : call VerifVolumes(elem, hybrid,nets,nete)
505 : end if
506 0 : end subroutine InitControlVolumes
507 :
508 0 : subroutine InitControlVolumes_duel(elem, hybrid,nets,nete)
509 : use bndry_mod, only: bndry_exchange
510 : use edge_mod, only: edgeVpack, edgeVunpack, freeedgebuffer, freeghostbuffer3D
511 : use element_mod, only: element_t, element_var_coordinates, element_var_coordinates3d
512 : use hybrid_mod, only: hybrid_t
513 :
514 : use quadrature_mod, only: quadrature_t, gausslobatto
515 : use coordinate_systems_mod, only: cube_face_number_from_sphere
516 :
517 : integer, intent(in) :: nets,nete
518 : type(element_t), intent(in), target :: elem(:)
519 : type(hybrid_t), intent(in) :: hybrid
520 :
521 : type(quadrature_t) :: gll_pts
522 : type(cartesian3d_t) :: quad(4),corners3d(4)
523 : real(r8) :: cv_pts(0:np) !was kind=longdouble_kind in HOMME
524 : real(r8) :: test(np,np,1)
525 :
526 : integer :: i, j, ie, k, kmax2, kk
527 :
528 0 : gll_pts = gausslobatto(np)
529 : ! gll points
530 0 : cv_pts(0)=-1
531 0 : do i=1,np
532 0 : cv_pts(i) = cv_pts(i-1) + gll_pts%weights(i)
533 : end do
534 0 : cv_pts(np)=1
535 0 : do i=1,np-1
536 0 : if (gll_pts%points(i) > cv_pts(i) .or. cv_pts(i) > gll_pts%points(i+1)) then
537 0 : call endrun("Error: CV and GLL points not interleaved")
538 : end if
539 : end do
540 :
541 :
542 : ! intialize local element areas
543 0 : test = 0
544 0 : do ie=nets,nete
545 0 : cvlist(ie)%cart3d_dual(0:np,0:np) = element_var_coordinates3D(elem(ie)%corners3D, cv_pts)
546 :
547 : ! compute true area of element and SEM area
548 0 : cvlist(ie)%vol=0
549 0 : do i=1,np
550 0 : do j=1,np
551 : ! (gnomonic coordinate lines only), more accurate
552 0 : quad(1) = cvlist(ie)%cart3d_dual(i-1,j-1)
553 0 : quad(2) = cvlist(ie)%cart3d_dual(i,j-1)
554 0 : quad(3) = cvlist(ie)%cart3d_dual(i,j)
555 0 : quad(4) = cvlist(ie)%cart3d_dual(i-1,j)
556 0 : cvlist(ie)%vol(i,j) = surfarea(quad,4)
557 : end do
558 : end do
559 0 : test(:,:,1) = cvlist(ie)%vol(:,:)
560 0 : call edgeVpack(edge1,test,1,0,ie)
561 : end do
562 :
563 0 : call bndry_exchange(hybrid, edge1,location='InitControlVolumes_duel')
564 :
565 0 : test = 0
566 0 : do ie=nets,nete
567 0 : test(:,:,1) = cvlist(ie)%vol(:,:)
568 0 : call edgeVunpack(edge1, test, 1, 0, ie)
569 0 : cvlist(ie)%totvol(:,:) = test(:,:,1)
570 0 : cvlist(ie)%invvol(:,:)=1.0_r8/cvlist(ie)%totvol(:,:)
571 : end do
572 :
573 0 : call VerifyAreas(elem, hybrid, nets, nete)
574 :
575 : ! construct the global CV grid and global CV areas from the
576 : ! local dual grid (cvlist()%cart_dual) and local areas (cvlist()%vol)
577 0 : call construct_cv_duel(elem, hybrid, nets, nete)
578 : ! compute output needed for SCRIP: lat/lon coordinates, and for the
579 : ! control volume with only 3 corners, repeat the last point to make a
580 : ! degenerate quad.
581 0 : kmax2 = 0
582 0 : do ie = nets, nete
583 0 : kmax2 = MAX(kmax2, MAXVAL(cvlist(ie)%nvert))
584 : end do
585 0 : do ie = nets, nete
586 0 : do j = 1, np
587 0 : do i = 1, np
588 0 : cvlist(ie)%vert_latlon(:,i,j)%lat = 0.0_r8
589 0 : cvlist(ie)%vert_latlon(:,i,j)%lon = 0.0_r8
590 0 : k = cvlist(ie)%nvert(i,j)
591 : !
592 : ! follow SCRIP protocol - of kk>k then repeat last vertex
593 : !
594 0 : do kk = k+1, nv_max
595 0 : cvlist(ie)%vert(kk, i, j) = cvlist(ie)%vert(k,i,j)
596 : end do
597 0 : do kk = 1, nv_max
598 0 : cvlist(ie)%vert_latlon(kk, i, j) = change_coordinates(cvlist(ie)%vert(kk, i, j))
599 0 : cvlist(ie)%face_no(kk, i, j) = cube_face_number_from_sphere(cvlist(ie)%vert_latlon(kk, i, j))
600 : end do
601 :
602 : end do
603 : end do
604 : end do
605 : ! Release memory
606 0 : if(hybrid%masterthread) then
607 0 : call freeedgebuffer(edge1)
608 0 : call FreeGhostBuffer3D(ghost_buf)
609 : end if
610 :
611 0 : initialized=.true.
612 0 : end subroutine InitControlVolumes_duel
613 :
614 0 : function average(t, n) result(a)
615 :
616 : integer, intent(in) :: n
617 : type(cartesian3d_t), intent(in) :: t(n)
618 : type(cartesian3d_t) :: a
619 : integer :: i
620 :
621 0 : a%x = 0._r8
622 0 : a%y = 0._r8
623 0 : a%z = 0._r8
624 0 : do i = 1, n
625 0 : a%x = a%x + t(i)%x
626 0 : a%y = a%y + t(i)%y
627 0 : a%z = a%z + t(i)%z
628 : end do
629 0 : a%x = a%x / n
630 0 : a%y = a%y / n
631 0 : a%z = a%z / n
632 0 : return
633 0 : end function average
634 :
635 0 : function make_unique(a, n) result(m)
636 :
637 : integer, intent(in) :: n
638 : real(r8), intent(inout) :: a(n)
639 : integer :: m
640 : integer :: i,j
641 : real(r8) :: delta
642 :
643 0 : do i=1,n-1
644 0 : do j=i+1,n
645 : ! if (ABS(a(j)-a(i)).lt. 1e-6) a(j) = 9999
646 0 : delta = abs(a(j)-a(i))
647 0 : if (delta < 1.e-6_r8) a(j) = 9999.0_r8
648 0 : if (abs((2.0_r8*pi) - delta) < 1.0e-6_r8) a(j) = 9999.0_r8
649 : end do
650 : end do
651 0 : m = 0
652 0 : do i=1,n
653 0 : if (a(i) < 9000.0_r8) m = m + 1
654 : end do
655 0 : if (mod(m,2).ne.0) then
656 0 : do i=1,n
657 0 : print *,'angle with centroid: ',i,a(i),mod(a(i),2*pi)
658 : end do
659 0 : call endrun("Error: Found an odd number or nodes for cv element. Should be even.")
660 : end if
661 : return
662 : end function make_unique
663 :
664 0 : function SortNodes(t3, n) result(m)
665 : use coordinate_systems_mod, only: cube_face_number_from_cart, cart2cubedsphere, change_coordinates
666 :
667 :
668 : integer, intent(in) :: n
669 : type(cartesian3d_t), intent(inout) :: t3(n)
670 :
671 0 : type(cartesian3d_t) :: c3, t(n)
672 : type(cartesian2d_t) :: c2, t2
673 0 : real(r8) :: angle(n)
674 : integer :: i,j,k,m,f
675 0 : integer :: ip(n)
676 :
677 0 : c3 = average(t3, n)
678 0 : f = cube_face_number_from_cart(c3)
679 0 : c2 = cart2cubedsphere(c3, f)
680 :
681 0 : do i=1,n
682 0 : t2 = cart2cubedsphere(t3(i), f)
683 0 : t2%x = t2%x - c2%x
684 0 : t2%y = t2%y - c2%y
685 0 : angle(i) = atan2(t2%y, t2%x)
686 : end do
687 0 : m = make_unique(angle,n)
688 0 : do i=1,m
689 : k = 1
690 0 : do j=2,n
691 0 : if (angle(j)<angle(k)) k=j
692 : end do
693 0 : angle(k) = 9999 ! greater than pi
694 0 : ip(i)=k
695 : end do
696 0 : t(1:m) = t3(ip(1:m))
697 0 : t3(1:m) = t(1:m)
698 : return
699 : end function SortNodes
700 :
701 0 : subroutine construct_cv_duel(elem,hybrid,nets,nete)
702 : !
703 : ! construct global dual grid from local element dual grid cvlist(ie)%cart3d_dual(:,:)
704 : ! all control volumes will be squares or triangles (at cube corners)
705 : !
706 : ! 10/2009: added option to make hexagon control volumes at cube edges and corners
707 : !
708 : use element_mod, only: element_t
709 : use hybrid_mod, only: hybrid_t
710 : use edge_mod, only: ghostVpack3D,ghostVunpack3d
711 : use bndry_mod , only: ghost_exchangeVfull
712 : use dimensions_mod, only: max_corner_elem
713 : use control_mod, only: north, south, east, west, neast, nwest, seast, swest
714 :
715 : type(element_t), intent(in), target :: elem(:)
716 : type(hybrid_t), intent(in) :: hybrid
717 : integer :: nets,nete
718 : ! local
719 : integer :: i,j,k,m,n,o,p,ie,m2
720 : real(r8) :: vertpack ( 0:np, 0:np, 3)
721 : real(r8) :: vertpack2 ( 1:np+1, 1:np+1, 3)
722 : real(r8) :: vertunpack( -1:np+1, -1:np+1, 3)
723 : real(r8) :: vertunpack2( 0:np+2, 0:np+2, 3)
724 0 : real(r8) :: sw( -1: 0, -1: 0, 3, max_corner_elem-1)
725 0 : real(r8) :: se( np:np+1, -1: 0, 3, max_corner_elem-1)
726 0 : real(r8) :: ne( np:np+1, np:np+1, 3, max_corner_elem-1)
727 0 : real(r8) :: nw( -1: 0, np:np+1, 3, max_corner_elem-1)
728 0 : real(r8) :: sw2( 0: 1, 0: 1, 3, max_corner_elem-1)
729 0 : real(r8) :: se2( np+1:np+2, 0: 1, 3, max_corner_elem-1)
730 0 : real(r8) :: ne2( np+1:np+2, np+1:np+2, 3, max_corner_elem-1)
731 0 : real(r8) :: nw2( 0: 1, np+1:np+2, 3, max_corner_elem-1)
732 :
733 0 : type(cartesian3d_t) :: vert(2*nv_max)
734 : type(cartesian3d_t) :: pt_3d
735 : type(cartesian3d_t) :: cv (-1:np+1, -1:np+1)
736 0 : type(cartesian3d_t) :: cv_sw(-1: 0, -1: 0, max_corner_elem-1)
737 0 : type(cartesian3d_t) :: cv_se(np:np+1, -1: 0, max_corner_elem-1)
738 0 : type(cartesian3d_t) :: cv_ne(np:np+1, np:np+1, max_corner_elem-1)
739 0 : type(cartesian3d_t) :: cv_nw(-1: 0, np:np+1, max_corner_elem-1)
740 : integer :: mlt(5:8)
741 :
742 :
743 0 : vertpack = 0
744 0 : vertunpack = 0
745 0 : vertpack2 = 0
746 0 : vertunpack2= 0
747 0 : sw = 0
748 0 : se = 0
749 0 : ne = 0
750 0 : nw = 0
751 0 : do ie=nets,nete
752 0 : do j=0,np
753 0 : do i=0,np
754 0 : pt_3d = cvlist(ie)%cart3d_dual(i,j)
755 0 : vertpack(i,j,1) = pt_3d%x
756 0 : vertpack(i,j,2) = pt_3d%y
757 0 : vertpack(i,j,3) = pt_3d%z
758 : end do
759 : end do
760 0 : do j=0,np
761 0 : do i=0,np
762 0 : do k=1,3
763 0 : vertpack2(i+1,j+1,k) = vertpack(i,j,k)
764 : end do
765 : end do
766 : end do
767 0 : call ghostVpack3D(ghost_buf, vertpack2, 3, 0, elem(ie)%desc)
768 : end do
769 :
770 0 : call ghost_exchangeVfull(hybrid%par,hybrid%ithr,ghost_buf)
771 0 : do ie=nets,nete
772 0 : do j=0,np
773 0 : do i=0,np
774 0 : pt_3d = cvlist(ie)%cart3d_dual(i,j)
775 0 : vertunpack(i,j,1) = pt_3d%x
776 0 : vertunpack(i,j,2) = pt_3d%y
777 0 : vertunpack(i,j,3) = pt_3d%z
778 : end do
779 : end do
780 0 : do j=0,np
781 0 : do i=0,np
782 0 : do k=1,3
783 0 : vertunpack2(i+1,j+1,k) = vertunpack(i,j,k)
784 : end do
785 : end do
786 : end do
787 0 : sw2=0
788 0 : se2=0
789 0 : nw2=0
790 0 : ne2=0
791 0 : call ghostVunpack3d(ghost_buf, vertunpack2, 3, 0, elem(ie)%desc, sw2,se2,nw2,ne2,mlt)
792 0 : do j=0,np+2
793 0 : do i=0,np+2
794 0 : do k=1,3
795 0 : vertunpack(i-1,j-1,k) = vertunpack2(i,j,k)
796 : end do
797 : end do
798 : end do
799 0 : sw=0
800 0 : se=0
801 0 : nw=0
802 0 : ne=0
803 0 : do m=1,mlt(swest)-1
804 0 : do k=1,3
805 0 : do j=0,1
806 0 : do i=0,1
807 0 : sw(i-1,j-1,k,m) = sw2(i,j,k,m)
808 : end do
809 : end do
810 : end do
811 : end do
812 0 : do m=1,mlt(seast)-1
813 0 : do k=1,3
814 0 : do j=0,1
815 0 : do i=np+1,np+2
816 0 : se(i-1,j-1,k,m) = se2(i,j,k,m)
817 : end do
818 : end do
819 : end do
820 : end do
821 0 : do m=1,mlt(nwest)-1
822 0 : do k=1,3
823 0 : do j=np+1,np+2
824 0 : do i=0,1
825 0 : nw(i-1,j-1,k,m) = nw2(i,j,k,m)
826 : end do
827 : end do
828 : end do
829 : end do
830 0 : do m=1,mlt(neast)-1
831 0 : do k=1,3
832 0 : do j=np+1,np+2
833 0 : do i=np+1,np+2
834 0 : ne(i-1,j-1,k,m) = ne2(i,j,k,m)
835 : end do
836 : end do
837 : end do
838 : end do
839 : ! Count and orient vert array
840 : ! Positive: 1->2->3->4 is counter clockwise on the sphere
841 : ! Negative: clockwise orientation
842 :
843 0 : do j=1,np
844 0 : do i=1,np
845 0 : cvlist(ie)%vert(:,i,j)%x = 0.0_r8
846 0 : cvlist(ie)%vert(:,i,j)%y = 0.0_r8
847 0 : cvlist(ie)%vert(:,i,j)%z = 0.0_r8
848 : end do
849 : end do
850 :
851 0 : do j=-1,np+1
852 0 : do i=-1,np+1
853 0 : cv(i,j)%x = vertunpack(i,j,1)
854 0 : cv(i,j)%y = vertunpack(i,j,2)
855 0 : cv(i,j)%z = vertunpack(i,j,3)
856 : end do
857 : end do
858 :
859 0 : do j=-1,0
860 0 : do i=-1,0
861 0 : do k=1,mlt(swest)-1
862 0 : cv_sw(i,j,k)%x = sw(i,j,1,k)
863 0 : cv_sw(i,j,k)%y = sw(i,j,2,k)
864 0 : cv_sw(i,j,k)%z = sw(i,j,3,k)
865 : end do
866 : end do
867 : end do
868 0 : do j=-1,0
869 0 : do i=np,np+1
870 0 : do k=1,mlt(seast)-1
871 0 : cv_se(i,j,k)%x = se(i,j,1,k)
872 0 : cv_se(i,j,k)%y = se(i,j,2,k)
873 0 : cv_se(i,j,k)%z = se(i,j,3,k)
874 : end do
875 : end do
876 : end do
877 0 : do j=np,np+1
878 0 : do i=-1,0
879 0 : do k=1,mlt(nwest)-1
880 0 : cv_nw(i,j,k)%x = nw(i,j,1,k)
881 0 : cv_nw(i,j,k)%y = nw(i,j,2,k)
882 0 : cv_nw(i,j,k)%z = nw(i,j,3,k)
883 : end do
884 : end do
885 : end do
886 0 : do j=np,np+1
887 0 : do i=np,np+1
888 0 : do k=1,mlt(neast)-1
889 0 : cv_ne(i,j,k)%x = ne(i,j,1,k)
890 0 : cv_ne(i,j,k)%y = ne(i,j,2,k)
891 0 : cv_ne(i,j,k)%z = ne(i,j,3,k)
892 : end do
893 : end do
894 : end do
895 :
896 0 : do j=2,np-1
897 0 : do i=2,np-1
898 : ! internal vertex on Cubed sphere
899 : ! Here is the order:
900 : !
901 : ! 4NW <- 3NE
902 : ! | ^
903 : ! v |
904 : ! 1SW -> 2SE
905 0 : vert(1) = cv(i-1, j-1)
906 0 : vert(2) = cv(i , j-1)
907 0 : vert(3) = cv(i , j )
908 0 : vert(4) = cv(i-1, j )
909 0 : cvlist(ie)%vert(1:4,i,j) = vert(1:4)
910 0 : cvlist(ie)%nvert(i,j) = 4
911 0 : m=4
912 : end do
913 : end do
914 :
915 0 : do j=0,np,np
916 0 : do i=2,np-1
917 0 : vert(1) = cv(i-1, j-1)
918 0 : vert(2) = cv(i , j-1)
919 0 : vert(3) = cv(i , j )
920 0 : vert(4) = cv(i , j+1)
921 0 : vert(5) = cv(i-1, j+1)
922 0 : vert(6) = cv(i-1, j )
923 0 : p = j
924 0 : if (p.eq.0) p=1
925 0 : cvlist(ie)%vert(1:6,i,p) = vert(1:6)
926 0 : cvlist(ie)%nvert(i,p) = 6
927 0 : m=6
928 : end do
929 : end do
930 :
931 0 : do j=2,np-1
932 0 : do i=0,np,np
933 0 : vert(1) = cv(i-1, j-1)
934 0 : vert(2) = cv(i , j-1)
935 0 : vert(3) = cv(i+1, j-1)
936 0 : vert(4) = cv(i+1, j )
937 0 : vert(5) = cv(i , j )
938 0 : vert(6) = cv(i-1, j )
939 0 : o = i
940 0 : if (o.eq.0) o=1
941 0 : cvlist(ie)%vert(1:6,o,j) = vert(1:6)
942 0 : cvlist(ie)%nvert(o,j) = 6
943 0 : m=6
944 : end do
945 : end do
946 0 : do j=0,np,np
947 0 : do i=0,np,np
948 0 : m = 0
949 0 : vert(:)%x = 0
950 0 : vert(:)%y = 0
951 0 : vert(:)%z = 0
952 0 : if (i.eq.0.and.j.eq.0) then
953 : ! counterclockwise from lower right
954 0 : vert(m+1) = cv(i+1, j-1) ! 5 4
955 0 : vert(m+2) = cv(i+1, j ) ! (-1,+1) (0,+1) (+1,+1) 3
956 0 : vert(m+3) = cv(i+1, j+1) !
957 0 : vert(m+4) = cv(i , j+1) ! (-1, 0) (i, j) (+1, 0) 2
958 0 : vert(m+5) = cv(i-1, j+1) !
959 0 : vert(m+6) = cv(i-1, j ) ! X X (+1,-1) 1
960 0 : m = m + 6
961 0 : if (mlt(swest).ne.0) then
962 0 : vert(m+1) = cv(i-1, j-1)
963 0 : vert(m+2) = cv(i , j-1)
964 0 : m = m+2
965 0 : do k=1,mlt(swest)-1 ! Bummer, toss in (-1,0) because transpose is undetectable
966 0 : vert(m+1) = cv_sw(i-1, j , k)
967 0 : vert(m+2) = cv_sw(i-1, j-1, k)
968 0 : vert(m+3) = cv_sw(i , j-1, k)
969 0 : m=m+3
970 : end do
971 : end if
972 : end if
973 0 : if (i.eq.np.and.j.eq.0) then
974 0 : if (mlt(seast).ne.0) then
975 0 : vert(m+1) = cv(i+1, j-1)
976 0 : vert(m+2) = cv(i+1, j )
977 0 : m = m+2
978 0 : do k=1,mlt(seast)-1
979 0 : vert(m+1) = cv_se(i , j-1, k)
980 0 : vert(m+2) = cv_se(i+1, j-1, k)
981 0 : vert(m+3) = cv_se(i+1, j , k)
982 0 : m=m+3
983 : end do
984 : end if
985 0 : vert(m+1) = cv(i+1, j+1)
986 0 : vert(m+2) = cv(i , j+1)
987 0 : vert(m+3) = cv(i-1, j+1)
988 0 : vert(m+4) = cv(i-1, j )
989 0 : vert(m+5) = cv(i-1, j-1)
990 0 : vert(m+6) = cv(i , j-1)
991 0 : m = m + 6
992 : end if
993 0 : if (i.eq.np.and.j.eq.np) then
994 0 : vert(1) = cv(i+1, j-1)
995 0 : vert(2) = cv(i+1, j )
996 0 : m = m + 2
997 0 : if (mlt(neast).ne.0) then
998 0 : vert(m+1) = cv(i+1, j+1)
999 0 : vert(m+2) = cv(i , j+1)
1000 0 : m = m+2
1001 0 : do k=1,mlt(neast)-1
1002 0 : vert(m+1) = cv_ne(i+1, j , k)
1003 0 : vert(m+2) = cv_ne(i+1, j+1, k)
1004 0 : vert(m+3) = cv_ne(i , j+1, k)
1005 0 : m=m+3
1006 : end do
1007 : end if
1008 0 : vert(m+1) = cv(i-1, j+1)
1009 0 : vert(m+2) = cv(i-1, j )
1010 0 : vert(m+3) = cv(i-1, j-1)
1011 0 : vert(m+4) = cv(i , j-1)
1012 0 : m = m + 4
1013 : end if
1014 0 : if (i.eq.0.and.j.eq.np) then
1015 0 : vert(m+1) = cv(i+1, j-1)
1016 0 : vert(m+2) = cv(i+1, j )
1017 0 : vert(m+3) = cv(i+1, j+1)
1018 0 : vert(m+4) = cv(i , j+1)
1019 0 : m = m + 4
1020 0 : if (mlt(nwest).ne.0) then
1021 0 : vert(m+1) = cv(i-1, j+1)
1022 0 : vert(m+2) = cv(i-1, j )
1023 0 : m = m+2
1024 0 : do k=1,mlt(nwest)-1
1025 0 : vert(m+1) = cv_nw(i , j+1, k)
1026 0 : vert(m+2) = cv_nw(i-1, j+1, k)
1027 0 : vert(m+3) = cv_nw(i-1, j , k)
1028 0 : m=m+3
1029 : end do
1030 : end if
1031 0 : vert(m+1) = cv(i-1, j-1)
1032 0 : vert(m+2) = cv(i , j-1)
1033 0 : m = m + 2
1034 : end if
1035 0 : o = i
1036 0 : p = j
1037 0 : if (o.eq.0) o=1
1038 0 : if (p.eq.0) p=1
1039 0 : m2=m
1040 0 : if (8 < m) then
1041 0 : m = SortNodes(vert, m2)
1042 : end if
1043 0 : if (m > nv_max) then
1044 0 : call endrun("error: vert dimensioned too small")
1045 : end if
1046 0 : cvlist(ie)%vert(1:m,o,p) = vert(1:m)
1047 0 : cvlist(ie)%nvert(o,p) = m
1048 : end do
1049 : end do
1050 : end do
1051 0 : end subroutine construct_cv_duel
1052 :
1053 0 : function SurfArea( cv, nvert ) result(area)
1054 :
1055 : type(cartesian3D_t), intent(in) :: cv(:)
1056 : integer, intent(in) :: nvert
1057 :
1058 : real(kind=r8) :: area, area1, area2, area3
1059 :
1060 0 : if (abs(nvert) == 3 ) then
1061 0 : area2 = 0.0_r8
1062 0 : area3 = 0.0_r8
1063 0 : if (cv(1)%x == 0) then
1064 0 : call sphere_tri_area(cv(2), cv(3), cv(4), area1)
1065 0 : else if (cv(2)%x == 0) then
1066 0 : call sphere_tri_area(cv(1), cv(3), cv(4), area1)
1067 0 : else if (cv(3)%x == 0) then
1068 0 : call sphere_tri_area(cv(1), cv(2), cv(4), area1)
1069 0 : else if (cv(4)%x == 0) then
1070 : call sphere_tri_area(cv(1), cv(2), cv(3), area1)
1071 : else
1072 0 : write(iulog, *) cv(1)%x, cv(1)%y
1073 0 : write(iulog, *) cv(2)%x, cv(2)%y
1074 0 : write(iulog, *) cv(3)%x, cv(3)%y
1075 0 : write(iulog, *) cv(4)%x, cv(4)%y
1076 0 : write(iulog, *) 'SurfArea error: should never happen'
1077 0 : call shr_sys_flush(iulog)
1078 0 : call endrun('SurfArea: invalid cv coordinates')
1079 : end if
1080 0 : else if (abs(nvert)==4) then
1081 0 : call sphere_tri_area(cv(1), cv(2), cv(3), area1)
1082 0 : call sphere_tri_area(cv(1), cv(3), cv(4), area2)
1083 0 : area3 = 0.0_r8
1084 :
1085 0 : else if (abs(nvert)==5) then
1086 0 : call sphere_tri_area(cv(1),cv(2),cv(3),area1)
1087 0 : call sphere_tri_area(cv(1),cv(3),cv(4),area2)
1088 0 : call sphere_tri_area(cv(1),cv(4),cv(5),area3)
1089 : else
1090 0 : call endrun('SurfArea: nvert > 5 not yet supported')
1091 : end if
1092 0 : area = area1 + area2 + area3
1093 0 : end function SurfArea
1094 :
1095 : ! ^
1096 : ! |dy o
1097 : ! |
1098 : ! (x,y) ---->dx
1099 0 : function SurfArea_dxdy(dx, dy, corner) result(integral)
1100 : use quadrature_mod, only: quadrature_t
1101 :
1102 : real(r8), intent(in) :: dx, dy
1103 : type(cartesian2d_t), intent(in) :: corner
1104 : real(r8) :: integral
1105 :
1106 : real(r8) :: alpha, beta, a1, a2, a3, a4
1107 :
1108 : ! cubed-sphere cell area, from Lauritzen & Nair MWR 2008
1109 : ! central angles:
1110 : ! cube face: -pi/4,-pi/4 -> pi/4,pi/4
1111 : ! this formula gives 2 so normalize by 4pi/6 / 2 = pi/3
1112 0 : alpha = corner%x
1113 0 : beta = corner%y
1114 0 : a1 = acos(-sin(alpha)*sin(beta)) ! 2.094
1115 0 : a2 = -acos(-sin(alpha+dx)*sin(beta) ) ! -1.047
1116 0 : a3 =- acos(-sin(alpha)*sin(beta+dy) ) ! -1.047
1117 0 : a4 = acos(-sin(alpha+dx)*sin(beta+dy) ) ! 2.094
1118 0 : integral = (a1+a2+a3+a4)
1119 : return
1120 : end function SurfArea_dxdy
1121 :
1122 0 : function find_intersect(x1in, x2in, y1in, y2in) result(sect)
1123 :
1124 : type(cartesian2D_t), intent(in) :: x1in, x2in, y1in, y2in
1125 : type(cartesian2D_t) :: sect
1126 :
1127 : type(cartesian2D_t) :: x, y, b, x1, x2, y1, y2
1128 : real(kind=r8) :: s1, s2, detA
1129 :
1130 : ! x1 + (x2-x1)*s1 = y1 + (y2-y1)*s2
1131 : ! b = y1-x1
1132 : ! x=x2-x1
1133 : ! y=y2-y1
1134 : ! x s1 - y s2 = b
1135 : ! x(1) s1 - y(1) s2 = b(1)
1136 : ! x(2) s1 - y(2) s2 = b(2)
1137 : !
1138 : ! x(1) -y(1) s1 = b(1) A s = b
1139 : ! x(2) -y(2) s2 = b(2)
1140 : !
1141 : ! A2= -y(2) y(1)
1142 : ! -x(2) x(1) s = A2 * b /detA
1143 :
1144 : ! convert to gnomonic
1145 0 : x1%x = tan(x1in%x)
1146 0 : x2%x = tan(x2in%x)
1147 0 : y1%x = tan(y1in%x)
1148 0 : y2%x = tan(y2in%x)
1149 0 : x1%y = tan(x1in%y)
1150 0 : x2%y = tan(x2in%y)
1151 0 : y1%y = tan(y1in%y)
1152 0 : y2%y = tan(y2in%y)
1153 :
1154 0 : x%x = x2%x-x1%x
1155 0 : x%y = x2%y-x1%y
1156 0 : y%x = y2%x-y1%x
1157 0 : y%y = y2%y-y1%y
1158 0 : b%x = y1%x-x1%x
1159 0 : b%y = y1%y-x1%y
1160 :
1161 0 : detA = -x%x*y%y + x%y*y%x
1162 :
1163 0 : s1 = (-y%y*b%x + y%x*b%y )/detA
1164 0 : s2 = (-x%y*b%x + x%x*b%y )/detA
1165 :
1166 0 : sect%x = x1%x + (x2%x-x1%x)*s1
1167 0 : sect%y = x1%y + (x2%y-x1%y)*s1
1168 :
1169 0 : sect%x = (sect%x + y1%x + (y2%x-y1%x)*s2)/2
1170 0 : sect%y = (sect%y + y1%y + (y2%y-y1%y)*s2)/2
1171 :
1172 0 : if (s1<0 .or. s1>1) then
1173 0 : write(iulog, *) 'failed: intersection: ',s1,s2
1174 0 : call shr_sys_flush(iulog)
1175 0 : call endrun('find_intersect: intersection failure')
1176 : end if
1177 :
1178 : ! convert back to equal angle:
1179 0 : sect%x = atan(sect%x)
1180 0 : sect%y = atan(sect%y)
1181 0 : end function find_intersect
1182 :
1183 0 : subroutine pentagon_iteration(sq1,sq2,pent,asq1,asq2,apent,faceno,anorm)
1184 : ! sq2
1185 : ! 4 3
1186 : ! 1 2
1187 : !
1188 : ! sq1 4 3
1189 : ! 2 1 5 pent
1190 : ! 3 4 1 2
1191 : !
1192 : !
1193 : ! d/dt sq1(1) = (area(sq1)-asq1) * [ com(sq1)-sq1(1) ]
1194 : ! +(area(sq2)-asq2) * [ com(sq2)-sq1(1) ]
1195 : ! +(area(pent)-apent) * [ com(pent)-sq1(1) ]
1196 : !
1197 : !
1198 : !
1199 : type(cartesian2d_t), intent(inout) :: sq1(4), sq2(4), pent(5)
1200 : real(r8), intent(in) :: asq1, asq2, apent, anorm
1201 : integer, intent(in) :: faceno
1202 :
1203 : type(cartesian3D_t) :: sq1_3d(4), sq2_3d(4), pent_3d(5)
1204 : real(r8) :: isq1, isq2, ipent, diff1, diff2, diffp, err
1205 : real(r8), parameter :: dt = .5_r8
1206 : real(r8), parameter :: tol_pentagon_iteration = 1.0e-10_r8
1207 : type(cartesian2d_t) :: sq1com, sq2com, pentcom, ds1, ds2
1208 : integer :: i, iter
1209 : integer, parameter :: iter_max = 10000
1210 :
1211 : ! compute center of mass:
1212 0 : sq1com%x = sum(sq1(:)%x)/4
1213 0 : sq1com%y = sum(sq1(:)%y)/4
1214 0 : sq2com%x = sum(sq2(:)%x)/4
1215 0 : sq2com%y = sum(sq2(:)%y)/4
1216 0 : pentcom%x = sum(pent(:)%x)/5
1217 0 : pentcom%y = sum(pent(:)%y)/5
1218 :
1219 0 : do i = 1, 4
1220 0 : sq1_3d(i)=cubedsphere2cart(sq1(i),faceno )
1221 0 : sq2_3d(i)=cubedsphere2cart(sq2(i),faceno )
1222 0 : pent_3d(i)=cubedsphere2cart(pent(i),faceno )
1223 : end do
1224 0 : pent_3d(5)=cubedsphere2cart(pent(5),faceno )
1225 :
1226 0 : do iter = 1, iter_max
1227 0 : isq1 = SurfArea(sq1_3d,4)
1228 0 : isq2 = SurfArea(sq2_3d,4)
1229 0 : ipent = SurfArea(pent_3d,5)
1230 :
1231 : ! d/dt sq1(1) = (area(sq1)-asq1) * [ com(sq1)-sq1(1) ]
1232 : ! +(area(sq2)-asq2) * [ com(sq2)-sq1(1) ]
1233 : ! +(area(pent)-apent) * [ com(pent)-sq1(1) ]
1234 : !
1235 0 : diff1 = (isq1-asq1)/anorm
1236 0 : diff2 = (isq2-asq2)/anorm
1237 0 : diffp = (ipent-apent)/anorm
1238 :
1239 0 : err = abs(diff1) + abs(diff2) + abs(diffp)
1240 0 : if (err < tol_pentagon_iteration) exit
1241 0 : if (mod(iter,1000) == 0) then
1242 0 : write(iulog, '(i5,3e18.5)') iter, err
1243 0 : call shr_sys_flush(iulog)
1244 : end if
1245 :
1246 0 : ds1%x = diff1* ( sq1com%x - sq1(1)%x )
1247 0 : ds1%y = diff1* ( sq1com%y - sq1(1)%y )
1248 0 : ds1%x = ds1%x + diffp* ( pentcom%x - sq1(1)%x )
1249 0 : ds1%y = ds1%y + diffp* ( pentcom%y - sq1(1)%y )
1250 :
1251 0 : ds2%x = diff2* ( sq2com%x - sq2(1)%x )
1252 0 : ds2%y = diff2* ( sq2com%y - sq2(1)%y )
1253 0 : ds2%x = ds2%x + diffp* ( pentcom%x - sq2(1)%x )
1254 0 : ds2%y = ds2%y + diffp* ( pentcom%y - sq2(1)%y )
1255 :
1256 0 : sq1(1)%x = sq1(1)%x + dt*ds1%x
1257 0 : sq1(1)%y = sq1(1)%y + dt*ds1%y
1258 0 : sq2(1)%x = sq2(1)%x + dt*ds2%x
1259 0 : sq2(1)%y = sq2(1)%y + dt*ds2%y
1260 0 : pent(4)=sq2(1)
1261 0 : pent(5)=sq1(1)
1262 0 : sq1_3d(1)=cubedsphere2cart(sq1(1),faceno )
1263 0 : sq2_3d(1)=cubedsphere2cart(sq2(1),faceno )
1264 0 : pent_3d(4)=sq2_3d(1)
1265 0 : pent_3d(5)=sq1_3d(1)
1266 : end do
1267 0 : if (iter >= iter_max) then
1268 0 : write(iulog, *) 'pentagon iteration did not converge err=', err
1269 0 : call shr_sys_flush(iulog)
1270 : end if
1271 0 : end subroutine pentagon_iteration
1272 :
1273 0 : subroutine InitControlVolumes_gll(elem, hybrid,nets,nete)
1274 : use edge_mod, only: freeedgebuffer
1275 : use element_mod, only: element_t,element_coordinates
1276 : use hybrid_mod, only: hybrid_t
1277 :
1278 : use quadrature_mod, only: quadrature_t, gausslobatto
1279 : use dimensions_mod, only: nlev
1280 : use cube_mod, only: convert_gbl_index
1281 : use coordinate_systems_mod, only: cart2cubedsphere_failsafe, cart2cubedsphere
1282 : use coordinate_systems_mod, only: cube_face_number_from_sphere
1283 :
1284 : integer, intent(in) :: nets,nete
1285 : type(element_t), intent(in), target :: elem(:)
1286 : type(hybrid_t), intent(in) :: hybrid
1287 :
1288 : type(cartesian2d_t) :: cartp_com(np,np) ! center of mass
1289 : type(cartesian2d_t) :: cartp_nm1(0:np,0:np)
1290 : real(r8) :: delx_k,dely_k,sum_dbg,r
1291 : integer :: i,j,ie,k,kptr,gllpts,nvert,k2,ie1,je1,face_no,kinsert
1292 : integer :: iter,iter_max,i1,j1
1293 : real(r8) :: diff(np,np),diffy(np-1,np-1),diffx(np-1,np-1)
1294 0 : real(r8) :: dx,dy,a1(nets:nete),a2(nets:nete),d1(nets:nete),d1mid(nets:nete)
1295 : real(r8) :: d2,d1_global,d1_global_mid,sphere1,sphere2,diff2,diff3
1296 : real(r8) :: diff23,diff32,diff33,diff22
1297 : real(r8) :: gllnm1(0:np) !was longdouble_kind in HOMME
1298 : type(cartesian2d_t) :: corner,start,endd,cv_loc_2d(4,np,np),cvnew_loc_2d(4,np,np)
1299 0 : type(cartesian3D_t) :: cart,cv_loc_3d(nv_max,np,np)
1300 0 : type(cartesian3D_t) :: temp3d(nv_max)
1301 : type(cartesian2d_t) :: cartp2d(np,np)
1302 : type(cartesian2d_t) :: x1,x2,x3,x
1303 : type(cartesian2d_t) :: sq1(4),sq2(4),pent(5)
1304 : type(cartesian3D_t) :: x1_3d,x2_3d,x3_3d
1305 : type(quadrature_t) :: gll
1306 : type(cartesian2d_t) :: dir,dirsum
1307 : type(spherical_polar_t) :: polar_tmp(0:np,0:np)
1308 : real(r8) :: rvert,area1,area2,ave,lat(4),lon(4)
1309 : real(r8) :: s,ds,triarea,triarea_target
1310 : real(r8) :: xp1,xm1,yp1,ym1,sumdiff
1311 : real(r8) :: tiny = 1e-11_r8,norm
1312 : real(r8) :: tol = 2.e-11_r8 ! convergece outer iteration
1313 : real(r8) :: tol_pentagons = 1.e-13_r8 ! convergece pentagon iteration
1314 :
1315 : ! area difference to trigger pentagons.
1316 : ! if it is too small, we will have pentagons with 1 very short edges
1317 : ! accuracy of surfarea() with very thin triangles seems poor (1e-11)
1318 : ! ne=30 1e-3: add 648 pentagons. area ratio: 1.003
1319 : ! ne=30 1e-4: add 696 pentagons. area ratio: 1.000004102
1320 : ! ne=30 1e-5: add 696 pentagons. area ratio: 1.000004102
1321 : ! ne=240 1e-4: add 5688/ 345600 pentagons, area ratio: 1.0004
1322 : ! ne=240 1e-5: add 5736/ 345600 pentagons, area ratio: 1.000000078
1323 : real(r8) :: tol_use_pentagons=1.0e-5_r8
1324 : logical :: Debug=.FALSE.,keep
1325 :
1326 : integer :: face1,face2,found,ie_max,movex,movey,moved,ii,kmax,kk
1327 : integer :: nskip,npent
1328 0 : integer :: nskipie(nets:nete), npentie(nets:nete)
1329 : type(cartesian2d_t) :: vert1_2d, vert_2d,vert2_2d
1330 : type(cartesian3D_t) :: vert1,vert2,vert_inserted(7)
1331 :
1332 0 : kmax=4
1333 :
1334 0 : gll = gausslobatto(np)
1335 : ! mid point rule:
1336 0 : do i=1,np-1
1337 0 : gllnm1(i) = ( gll%points(i) + gll%points(i+1) ) /2
1338 : end do
1339 : ! check that gll(i) < gllnm1(i) < gll(i+1)
1340 0 : do i=1,np-1
1341 0 : if (gll%points(i) > gllnm1(i) .or. gllnm1(i) > gll%points(i+1)) then
1342 0 : call endrun("InitControlVolumes_gll: CV and GLL points not interleaved")
1343 : end if
1344 : end do
1345 0 : gllnm1(0)=-1
1346 0 : gllnm1(np)=1
1347 :
1348 : ! MNL: dx and dy are no longer part of element_t
1349 : ! but they are easily computed for the
1350 : ! uniform case
1351 0 : dx = pi/(2.0d0*dble(ne))
1352 0 : dy = dx
1353 :
1354 : ! intialize local element dual grid, local element areas
1355 :
1356 0 : do ie=nets,nete
1357 :
1358 0 : call convert_gbl_index(elem(ie)%vertex%number,ie1,je1,face_no)
1359 0 : start%x=-pi/4 + ie1*dx
1360 0 : start%y=-pi/4 + je1*dy
1361 0 : endd%x =start%x + dx
1362 0 : endd%y =start%y + dy
1363 0 : cartp_nm1(0:np,0:np) = element_coordinates(start,endd,gllnm1)
1364 0 : cvlist(ie)%cartp_dual = cartp_nm1
1365 :
1366 : ! compute true area of element and SEM area
1367 0 : a1(ie) = SurfArea_dxdy(dx,dy,elem(ie)%cartp(1,1))
1368 0 : a2(ie) = sum(elem(ie)%spheremp(:,:))
1369 0 : do i=1,np
1370 0 : do j=1,np
1371 : ! (gnomonic coordinate lines only), more accurate
1372 0 : delx_k = cartp_nm1(i,j-1)%x - cartp_nm1(i-1,j-1)%x
1373 0 : dely_k = cartp_nm1(i-1,j)%y - cartp_nm1(i-1,j-1)%y
1374 0 : cvlist(ie)%vol(i,j) = SurfArea_dxdy(delx_k,dely_k,cartp_nm1(i-1,j-1))
1375 : end do
1376 : end do
1377 0 : global_shared_buf(ie,1) = a1(ie)
1378 0 : global_shared_buf(ie,2) = a2(ie)
1379 : end do
1380 0 : call wrap_repro_sum(nvars=2, comm=hybrid%par%comm)
1381 0 : sphere1 = global_shared_sum(1)
1382 0 : sphere2 = global_shared_sum(2)
1383 :
1384 : ! construct the global CV grid and global CV areas from the
1385 : ! local dual grid (cvlist()%cart_dual) and local areas (cvlist()%vol)
1386 0 : call construct_cv_gll(elem,hybrid,nets,nete)
1387 :
1388 0 : iter_max=2000
1389 : if (iter_max>0) then
1390 : ! areas computed from eleemnts on boundaries are from hexagons and pentagons
1391 : ! compute new areas where all CVs are squares or triangles
1392 0 : do ie=nets,nete
1393 0 : do i=1,np
1394 0 : do j=1,np
1395 : ! ifort bug if we try this:
1396 : ! area2 = surfarea(cvlist(ie)%vert(1:4,i,j),cvlist(ie)%nvert(i,j))
1397 0 : cv_loc_3d(:,i,j)=cvlist(ie)%vert(:,i,j)
1398 0 : area2 = surfarea(cv_loc_3d(:,i,j),cvlist(ie)%nvert(i,j))
1399 0 : cvlist(ie)%totvol(i,j)=area2
1400 : end do
1401 : end do
1402 : end do
1403 : end if
1404 : ! iteration over cvlist(ie)%totvol
1405 0 : d1_global=0
1406 0 : do iter=1,iter_max
1407 0 : ie_max=-1
1408 0 : do ie=nets,nete
1409 : ! we want at each point, the gll_area = true_area
1410 : ! but sum(gll_area) = a2 and sum(true_area)=a1
1411 : ! so normalize so that: gll_area/a2 = true_area/a1, or gll_area = area*a2/a1
1412 :
1413 : ! requires more iterations, but the total volume within an
1414 : ! element is always correct
1415 0 : diff(:,:) = ( cvlist(ie)%vol(:,:) - elem(ie)%spheremp(:,:)*a1(ie)/a2(ie) )
1416 0 : sumdiff=sum( cvlist(ie)%vol(:,:)) - a1(ie)
1417 0 : diff(:,:) = diff(:,:)/(a1(ie)/(np*np))
1418 :
1419 :
1420 :
1421 : ! set boundary values (actually not used)
1422 0 : cartp_nm1 = cvlist(ie)%cartp_dual(0:np,0:np)
1423 : ! convert 9 cv corners in this element into cart_nm1 cubed-sphere coordiantes
1424 0 : do i=1,np-1
1425 0 : do j=1,np-1
1426 0 : cartp_nm1(i,j) = cart2cubedsphere( cvlist(ie)%vert(3,i,j),elem(ie)%FaceNum )
1427 : end do
1428 : end do
1429 : ! compute center of mass of control volumes:
1430 : ! todo: move points towards GLL points, not center of mass
1431 : ! center of mass could send up a feedback with CV points!
1432 0 : do i=1,np
1433 0 : do j=1,np
1434 0 : cart%x = sum( cvlist(ie)%vert(:,i,j)%x )/abs(cvlist(ie)%nvert(i,j))
1435 0 : cart%y = sum( cvlist(ie)%vert(:,i,j)%y )/abs(cvlist(ie)%nvert(i,j))
1436 0 : cart%z = sum( cvlist(ie)%vert(:,i,j)%z )/abs(cvlist(ie)%nvert(i,j))
1437 0 : cartp_com(i,j) = cart2cubedsphere( cart,elem(ie)%FaceNum )
1438 : end do
1439 : end do
1440 0 : d2=0
1441 0 : do i=1,np-1
1442 0 : do j=1,np-1
1443 0 : dirsum%x=0
1444 0 : dirsum%y=0
1445 0 : movex=1
1446 0 : movey=1
1447 0 : moved=0
1448 :
1449 0 : do i1=0,1
1450 0 : do j1=0,1
1451 : ! keep=.true. : .85/1.05
1452 : ! corners only: .93/1.07
1453 : ! corners and edges: .89/1.11
1454 0 : keep=.false.
1455 : ! corner volumes
1456 0 : if (i==1 .and. j==1) then
1457 0 : if (i1==0 .and. j1==0) keep=.true.
1458 : moved=1
1459 0 : else if (i==np-1 .and. j==1) then
1460 0 : if (i1==1 .and. j1==0) keep=.true.
1461 : moved=-1
1462 0 : else if (i==1 .and. j==np-1) then
1463 0 : if (i1==0 .and. j1==1) keep=.true.
1464 : moved=-1
1465 0 : else if (i==np-1 .and. j==np-1) then
1466 0 : if (i1==1 .and. j1==1) keep=.true.
1467 : moved=1
1468 : ! edge volumes
1469 :
1470 :
1471 0 : else if (i==1) then
1472 0 : if (i1==0) keep=.true.
1473 0 : else if (i==np-1) then
1474 0 : if (i1==1) keep=.true.
1475 0 : else if (j==1) then
1476 0 : if (j1==0) keep=.true.
1477 0 : else if (j==np-1) then
1478 0 : if (j1==1) keep=.true.
1479 : else
1480 : keep=.true.
1481 : end if
1482 0 : if (keep) then
1483 : ! error weighted direction towards center of mass of area
1484 : ! move towards grid point
1485 0 : dir%x = (elem(ie)%cartp(i+i1,j+j1)%x - cartp_nm1(i,j)%x )*(abs(diff(i+i1,j+j1)))
1486 0 : dir%y = (elem(ie)%cartp(i+i1,j+j1)%y - cartp_nm1(i,j)%y )*(abs(diff(i+i1,j+j1)))
1487 0 : if (moved==1) then
1488 : ! project onto (1,1)/sqrt(2)
1489 0 : dir%x = dir%x/sqrt(2d0) + dir%y/sqrt(2d0)
1490 0 : dir%y = dir%x
1491 : end if
1492 0 : if (moved==-1) then
1493 : ! project onto (-1,1)/sqrt(2)
1494 0 : dir%y = -dir%x/sqrt(2d0) + dir%y/sqrt(2d0)
1495 0 : dir%x = -dir%y
1496 : end if
1497 :
1498 :
1499 0 : if ( diff(i+i1,j+j1) > 0 ) then
1500 : ! this volume is too big, so move cv point towards grid center
1501 : ! weighted by length error
1502 0 : dirsum%x = dirsum%x + movex*dir%x
1503 0 : dirsum%y = dirsum%y + movey*dir%y
1504 : else
1505 0 : dirsum%x = dirsum%x - movex*dir%x
1506 0 : dirsum%y = dirsum%y - movey*dir%y
1507 : end if
1508 : end if
1509 : end do
1510 : end do
1511 0 : d2 = d2 + dirsum%x**2 + dirsum%y**2
1512 0 : cartp_nm1(i,j)%x = cartp_nm1(i,j)%x + 0.25_r8*dirsum%x
1513 0 : cartp_nm1(i,j)%y = cartp_nm1(i,j)%y + 0.25_r8*dirsum%y
1514 :
1515 : end do
1516 : end do
1517 0 : cvlist(ie)%cartp_dual(0:np,0:np) = cartp_nm1
1518 0 : d2=sqrt(d2)
1519 :
1520 0 : d1(ie)=sqrt(sum(diff**2))
1521 :
1522 0 : d1mid(ie)=d1(ie)
1523 : ! ignore center cv's:
1524 0 : diff(2:3,2:3)=0
1525 0 : d1mid(ie)=sqrt(sum(diff**2))
1526 :
1527 : end do ! ie loop
1528 0 : dx=maxval(d1)
1529 0 : d1_global = ParallelMax(dx,hybrid)
1530 0 : dx=maxval(d1mid)
1531 0 : d1_global_mid = ParallelMax(dx,hybrid)
1532 0 : if (mod(iter-1,250).eq.0) then
1533 0 : if (hybrid%masterthread) write(iulog, *) iter,"max d1=",d1_global,d1_global_mid
1534 : end if
1535 : ! compute new global CV (cvlist(ie)%vert from cvlist(ie)%cartp_dual).
1536 : ! cvlist()%totarea incorrect since local volumes not computed above
1537 0 : call construct_cv_gll(elem,hybrid,nets,nete)
1538 :
1539 : ! update totvol (area of multi-element cv)
1540 0 : do ie=nets,nete
1541 0 : do i=1,np
1542 0 : do j=1,np
1543 : ! ifort bug if we try this:
1544 : ! area2 = surfarea(cvlist(ie)%vert(1:4,i,j),cvlist(ie)%nvert(i,j))
1545 0 : cv_loc_3d(:,i,j)=cvlist(ie)%vert(:,i,j)
1546 0 : area2 = surfarea(cv_loc_3d(:,i,j),cvlist(ie)%nvert(i,j))
1547 0 : cvlist(ie)%totvol(i,j) = area2
1548 0 : if (isnan(area2)) then
1549 0 : write(iulog, *) 'ie,i,j',ie,i,j
1550 0 : write(iulog, *) cvlist(ie)%nvert(i,j)
1551 0 : write(iulog, *) cv_loc_3d(1,i,j)
1552 0 : write(iulog, *) cv_loc_3d(2,i,j)
1553 0 : write(iulog, *) cv_loc_3d(3,i,j)
1554 0 : write(iulog, *) cv_loc_3d(4,i,j)
1555 0 : call shr_sys_flush(iulog)
1556 0 : call endrun('InitControlVolumes_gll: area = NaN')
1557 : end if
1558 : end do
1559 : end do
1560 : end do
1561 :
1562 : ! update %vol (local control volume within each element)
1563 0 : do ie=nets,nete
1564 0 : cartp2d = elem(ie)%cartp
1565 0 : do i=1,np
1566 0 : do j=1,np
1567 : ! ifort bug if we try this:
1568 : ! area2 = surfarea(cvlist(ie)%vert(1:4,i,j),cvlist(ie)%nvert(i,j))
1569 :
1570 0 : do ii=1,4
1571 : !
1572 : ! if we do not use _failsafe version of cart2cubedsphere code will fail with "-debug"
1573 : !
1574 0 : cv_loc_2d(ii,i,j) = cart2cubedsphere_failsafe( cvlist(ie)%vert(ii,i,j),elem(ie)%FaceNum )
1575 : end do
1576 0 : if (i==1 .and. j==1) then
1577 0 : cv_loc_2d(1,i,j)=cartp2d(i,j)
1578 : end if
1579 0 : if (i==np .and. j==1) then
1580 0 : cv_loc_2d(2,i,j)=cartp2d(i,j)
1581 : end if
1582 0 : if (i==1 .and. j==np) then
1583 0 : cv_loc_2d(4,i,j)=cartp2d(i,j)
1584 : end if
1585 0 : if (i==np .and. j==np) then
1586 0 : cv_loc_2d(3,i,j)=cartp2d(i,j)
1587 : end if
1588 :
1589 :
1590 0 : cvnew_loc_2d(:,i,j)=cv_loc_2d(:,i,j)
1591 :
1592 : !
1593 : ! 4NW <- 3NE
1594 : ! | ^
1595 : ! v |
1596 : ! 1SW -> 2SE
1597 0 : if (i==1) then
1598 : ! replace points with x< elem(ie)%vert(i,j)%x
1599 0 : if (cv_loc_2d(1,i,j)%x < cartp2d(i,j)%x) then
1600 : cvnew_loc_2d(1,i,j) = find_intersect(&
1601 : cv_loc_2d(1,i,j), cv_loc_2d(2,i,j),&
1602 0 : elem(ie)%cartp(i,1),elem(ie)%cartp(i,np))
1603 : end if
1604 0 : if (cv_loc_2d(4,i,j)%x < cartp2d(i,j)%x) then
1605 : cvnew_loc_2d(4,i,j) = find_intersect(&
1606 : cv_loc_2d(4,i,j), cv_loc_2d(3,i,j),&
1607 0 : elem(ie)%cartp(i,1),elem(ie)%cartp(i,np))
1608 : end if
1609 : end if
1610 :
1611 0 : if (i==np) then
1612 : ! replace points with x> elem(ie)%vert(i,j)%x
1613 0 : if (cv_loc_2d(2,i,j)%x > cartp2d(i,j)%x) then
1614 : cvnew_loc_2d(2,i,j) = find_intersect(&
1615 : cv_loc_2d(1,i,j), cv_loc_2d(2,i,j),&
1616 0 : elem(ie)%cartp(i,1),elem(ie)%cartp(i,np))
1617 : end if
1618 0 : if (cv_loc_2d(3,i,j)%x > cartp2d(i,j)%x) then
1619 : cvnew_loc_2d(3,i,j) = find_intersect(&
1620 : cv_loc_2d(4,i,j), cv_loc_2d(3,i,j),&
1621 0 : elem(ie)%cartp(i,1),elem(ie)%cartp(i,np))
1622 : end if
1623 : end if
1624 : !
1625 : ! 4NW <- 3NE
1626 : ! | ^
1627 : ! v |
1628 : ! 1SW -> 2SE
1629 0 : if (j==1) then
1630 : ! replace points with y < elem(ie)%vert(i,j)%y
1631 0 : if (cv_loc_2d(1,i,j)%y < cartp2d(i,j)%y) then
1632 : cvnew_loc_2d(1,i,j) = find_intersect(&
1633 : cv_loc_2d(1,i,j), cv_loc_2d(4,i,j),&
1634 0 : elem(ie)%cartp(1,j),elem(ie)%cartp(np,j))
1635 : end if
1636 0 : if (cv_loc_2d(2,i,j)%y < cartp2d(i,j)%y) then
1637 : cvnew_loc_2d(2,i,j) = find_intersect(&
1638 : cv_loc_2d(2,i,j), cv_loc_2d(3,i,j),&
1639 0 : elem(ie)%cartp(1,j),elem(ie)%cartp(np,j))
1640 : end if
1641 : end if
1642 0 : if (j==np) then
1643 : ! replace points with y > elem(ie)%vert(i,j)%y
1644 0 : if (cv_loc_2d(4,i,j)%y > cartp2d(i,j)%y) then
1645 : cvnew_loc_2d(4,i,j) = find_intersect(&
1646 : cv_loc_2d(1,i,j), cv_loc_2d(4,i,j),&
1647 0 : elem(ie)%cartp(1,j),elem(ie)%cartp(np,j))
1648 : end if
1649 0 : if (cv_loc_2d(3,i,j)%y > cartp2d(i,j)%y) then
1650 : cvnew_loc_2d(3,i,j) = find_intersect(&
1651 : cv_loc_2d(2,i,j), cv_loc_2d(3,i,j),&
1652 0 : elem(ie)%cartp(1,j),elem(ie)%cartp(np,j))
1653 : end if
1654 : end if
1655 0 : do ii=1,4
1656 0 : cv_loc_3d(ii,i,j)=cubedsphere2cart(cvnew_loc_2d(ii,i,j),elem(ie)%FaceNum )
1657 : end do
1658 0 : area2 = surfarea(cv_loc_3d(:,i,j),4)
1659 0 : cvlist(ie)%vol(i,j) = area2
1660 0 : if (isnan(area2)) then
1661 0 : write(iulog, *) 'ie,i,j',ie,i,j
1662 0 : write(iulog, *) cvlist(ie)%nvert(i,j)
1663 0 : write(iulog, *) cv_loc_3d(1,i,j)
1664 0 : write(iulog, *) cv_loc_3d(2,i,j)
1665 0 : write(iulog, *) cv_loc_3d(3,i,j)
1666 0 : write(iulog, *) cv_loc_3d(4,i,j)
1667 0 : call shr_sys_flush(iulog)
1668 0 : call endrun('InitControlVolumes_gll: area = NaN')
1669 : end if
1670 : end do
1671 : end do
1672 : end do
1673 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1674 :
1675 0 : if ( d1_global > 10.0_r8 .or. d1_global_mid < tol) then
1676 0 : if (hybrid%masterthread) then
1677 0 : write(iulog, *) 'first iteration stopping:'
1678 0 : write(iulog, *) iter, "max error=", d1_global_mid
1679 0 : call shr_sys_flush(iulog)
1680 : end if
1681 : exit
1682 : end if
1683 : end do ! iteration loop
1684 :
1685 0 : kmax=5
1686 :
1687 0 : nskip=0
1688 0 : npent=0
1689 0 : nskipie(:) = 0
1690 0 : npentie(:) = 0
1691 0 : do ie=nets,nete
1692 0 : diff = ( cvlist(ie)%vol(:,:) - elem(ie)%spheremp(:,:)*a1(ie)/a2(ie) )
1693 0 : if ( maxval(abs(diff(2:3,2:3)))/a1(ie) > tol_use_pentagons ) then
1694 0 : npent=npent+1
1695 0 : npentie(ie) = npentie(ie) + 1
1696 : !
1697 : ! 4NW <- 3NE
1698 : ! | ^
1699 : ! v | 23 33
1700 : ! 1SW -> 2SE 22 32
1701 0 : if (diff(2,2)>0 .and. diff(3,3)>0) then
1702 0 : x1 = cart2cubedsphere( cvlist(ie)%vert(3,2,2),elem(ie)%FaceNum )
1703 0 : x2 = cart2cubedsphere( cvlist(ie)%vert(1,2,2),elem(ie)%FaceNum )
1704 0 : s = .99_r8
1705 0 : x3%x = x2%x + (x1%x-x2%x)*s
1706 0 : x3%y = x2%y + (x1%y-x2%y)*s
1707 :
1708 0 : sq1(1) = x3
1709 0 : sq1(2) = cart2cubedsphere( cvlist(ie)%vert(4,2,2),elem(ie)%FaceNum )
1710 0 : sq1(3) = cart2cubedsphere( cvlist(ie)%vert(1,2,2),elem(ie)%FaceNum )
1711 0 : sq1(4) = cart2cubedsphere( cvlist(ie)%vert(2,2,2),elem(ie)%FaceNum )
1712 :
1713 0 : x2 = cart2cubedsphere( cvlist(ie)%vert(3,3,3),elem(ie)%FaceNum )
1714 0 : s = .99_r8
1715 0 : x3%x = x2%x + (x1%x-x2%x)*s
1716 0 : x3%y = x2%y + (x1%y-x2%y)*s
1717 :
1718 0 : sq2(1) = x3
1719 0 : sq2(2) = cart2cubedsphere( cvlist(ie)%vert(2,3,3),elem(ie)%FaceNum )
1720 0 : sq2(3) = cart2cubedsphere( cvlist(ie)%vert(3,3,3),elem(ie)%FaceNum )
1721 0 : sq2(4) = cart2cubedsphere( cvlist(ie)%vert(4,3,3),elem(ie)%FaceNum )
1722 :
1723 0 : pent(1) = cart2cubedsphere( cvlist(ie)%vert(1,3,2),elem(ie)%FaceNum )
1724 0 : pent(2) = cart2cubedsphere( cvlist(ie)%vert(2,3,2),elem(ie)%FaceNum )
1725 0 : pent(3) = cart2cubedsphere( cvlist(ie)%vert(3,3,2),elem(ie)%FaceNum )
1726 0 : pent(4) = sq2(1)
1727 0 : pent(5) = sq1(1)
1728 :
1729 : call pentagon_iteration(sq1,sq2,pent,&
1730 0 : elem(ie)%spheremp(2,2)*a1(ie)/a2(ie), &
1731 : elem(ie)%spheremp(3,3)*a1(ie)/a2(ie), &
1732 0 : elem(ie)%spheremp(3,2)*a1(ie)/a2(ie),elem(ie)%FaceNum,a1(ie))
1733 :
1734 0 : x2_3d=cubedsphere2cart(sq1(1),elem(ie)%FaceNum )
1735 0 : x3_3d=cubedsphere2cart(sq2(1),elem(ie)%FaceNum )
1736 :
1737 0 : cvlist(ie)%vert(3,2,2)=x2_3d
1738 0 : cvlist(ie)%vert(1,3,3)=x3_3d
1739 :
1740 0 : cvlist(ie)%vert(5,2,3)=cvlist(ie)%vert(4,2,3)
1741 0 : cvlist(ie)%vert(4,2,3)=cvlist(ie)%vert(3,2,3)
1742 0 : cvlist(ie)%vert(2,2,3)=x2_3d
1743 0 : cvlist(ie)%vert(3,2,3)=x3_3d
1744 :
1745 0 : cvlist(ie)%vert(5,3,2)=x2_3d
1746 0 : cvlist(ie)%vert(4,3,2)=x3_3d
1747 :
1748 0 : cvlist(ie)%nvert(2,3)=sign(5,cvlist(ie)%nvert(2,3))
1749 0 : cvlist(ie)%nvert(3,2)=sign(5,cvlist(ie)%nvert(3,2))
1750 0 : else if (diff(2,3) >0 .and. diff(3,2)>0) then
1751 : !
1752 : ! 4NW <- 3NE
1753 : ! | ^
1754 : ! v | 23 33
1755 : ! 1SW -> 2SE 22 32
1756 0 : x1 = cart2cubedsphere( cvlist(ie)%vert(2,2,3),elem(ie)%FaceNum )
1757 0 : x2 = cart2cubedsphere( cvlist(ie)%vert(4,2,3),elem(ie)%FaceNum )
1758 0 : s = .99_r8
1759 0 : x3%x = x2%x + (x1%x-x2%x)*s
1760 0 : x3%y = x2%y + (x1%y-x2%y)*s
1761 :
1762 0 : sq1(1) = x3
1763 0 : sq1(2) = cart2cubedsphere( cvlist(ie)%vert(3,2,3),elem(ie)%FaceNum )
1764 0 : sq1(3) = cart2cubedsphere( cvlist(ie)%vert(4,2,3),elem(ie)%FaceNum )
1765 0 : sq1(4) = cart2cubedsphere( cvlist(ie)%vert(1,2,3),elem(ie)%FaceNum )
1766 :
1767 0 : x2 = cart2cubedsphere( cvlist(ie)%vert(2,3,2),elem(ie)%FaceNum )
1768 0 : s = .99_r8
1769 0 : x3%x = x2%x + (x1%x-x2%x)*s
1770 0 : x3%y = x2%y + (x1%y-x2%y)*s
1771 :
1772 0 : sq2(1) = x3
1773 0 : sq2(2) = cart2cubedsphere( cvlist(ie)%vert(1,3,2),elem(ie)%FaceNum )
1774 0 : sq2(3) = cart2cubedsphere( cvlist(ie)%vert(2,3,2),elem(ie)%FaceNum )
1775 0 : sq2(4) = cart2cubedsphere( cvlist(ie)%vert(3,3,2),elem(ie)%FaceNum )
1776 :
1777 0 : pent(1) = cart2cubedsphere( cvlist(ie)%vert(4,2,2),elem(ie)%FaceNum )
1778 0 : pent(2) = cart2cubedsphere( cvlist(ie)%vert(1,2,2),elem(ie)%FaceNum )
1779 0 : pent(3) = cart2cubedsphere( cvlist(ie)%vert(2,2,2),elem(ie)%FaceNum )
1780 0 : pent(4) = sq2(1)
1781 0 : pent(5) = sq1(1)
1782 :
1783 : call pentagon_iteration(sq1,sq2,pent,&
1784 0 : elem(ie)%spheremp(2,3)*a1(ie)/a2(ie), &
1785 : elem(ie)%spheremp(3,2)*a1(ie)/a2(ie), &
1786 0 : elem(ie)%spheremp(2,2)*a1(ie)/a2(ie),elem(ie)%FaceNum,a1(ie))
1787 :
1788 0 : x2_3d=cubedsphere2cart(sq1(1),elem(ie)%FaceNum )
1789 0 : x3_3d=cubedsphere2cart(sq2(1),elem(ie)%FaceNum )
1790 :
1791 0 : cvlist(ie)%vert(2,2,3)=x2_3d
1792 :
1793 0 : cvlist(ie)%vert(4,3,2)=x3_3d
1794 :
1795 0 : cvlist(ie)%vert(5,2,2)=cvlist(ie)%vert(4,2,2)
1796 0 : cvlist(ie)%vert(3,2,2)=x3_3d
1797 0 : cvlist(ie)%vert(4,2,2)=x2_3d
1798 :
1799 :
1800 0 : cvlist(ie)%vert(1,3,3)=x3_3d
1801 0 : cvlist(ie)%vert(5,3,3)=x2_3d
1802 :
1803 0 : cvlist(ie)%nvert(3,3)=sign(5,cvlist(ie)%nvert(3,3))
1804 0 : cvlist(ie)%nvert(2,2)=sign(5,cvlist(ie)%nvert(2,2))
1805 : else
1806 0 : if (hybrid%masterthread) then
1807 0 : write(iulog, *) ie,'bad type'
1808 0 : call shr_sys_flush(iulog)
1809 : end if
1810 0 : call endrun('InitControlVolumes_gll: bad type')
1811 : end if
1812 : ! recompute areas:
1813 0 : do i=2,3
1814 0 : do j=2,3
1815 0 : nvert=abs(cvlist(ie)%nvert(i,j))
1816 0 : temp3d(1:nvert)=cvlist(ie)%vert(1:nvert,i,j)
1817 0 : cvlist(ie)%vol(i,j)=surfarea(temp3d,nvert)
1818 0 : cvlist(ie)%totvol(i,j)=cvlist(ie)%vol(i,j)
1819 : end do
1820 : end do
1821 : else
1822 : !write(iulog, *) 'skipping pentagon procedure ie=',ie
1823 : !write(iulog, *) 'maxval diff: ',maxval(abs(diff(2:3,2:3)))/a1(ie)
1824 0 : nskip=nskip+1
1825 0 : nskipie(ie) = nskipie(ie) + 1
1826 : end if
1827 0 : global_shared_buf(ie,1) = nskipie(ie)
1828 0 : global_shared_buf(ie,2) = npentie(ie)
1829 : end do
1830 0 : call wrap_repro_sum(nvars=2, comm=hybrid%par%comm)
1831 0 : nskip = global_shared_sum(1)
1832 0 : npent = global_shared_sum(2)
1833 0 : if (hybrid%masterthread) then
1834 0 : write(*,'(a,i7,a,i7)') "no. elements where pentagons were added: ",npent,"/",npent+nskip
1835 : end if
1836 :
1837 : ! compute output needed for SCRIP: lat/lon coordinates, and for the
1838 : ! control volume with only 3 corners, repeat the last point to make a
1839 : ! degenerate quad.
1840 0 : do ie=nets,nete
1841 0 : do j=1,np
1842 0 : do i=1,np
1843 0 : cvlist(ie)%vert_latlon(:,i,j)%lat = 0._r8
1844 0 : cvlist(ie)%vert_latlon(:,i,j)%lon = 0._r8
1845 0 : do k = 1, kmax
1846 0 : rvert = cvlist(ie)%vert(k,i,j)%x**2+cvlist(ie)%vert(k,i,j)%y**2+cvlist(ie)%vert(k,i,j)%z**2
1847 0 : if(rvert > 0.9_r8) then
1848 0 : cvlist(ie)%vert_latlon(k,i,j) = change_coordinates(cvlist(ie)%vert(k,i,j))
1849 : else
1850 : ! coordinates = 0, this corner was not set above because this point
1851 : ! only has 3 cells (corner point) pick either neighbor to make a degenerate quad
1852 0 : k2 = k - 1
1853 0 : if (k2 == 0) then
1854 0 : k2 = 2 ! can only happen for corner point with data in 2,3,4
1855 : end if
1856 0 : cvlist(ie)%vert_latlon(k,i,j) = change_coordinates(cvlist(ie)%vert(k2,i,j))
1857 0 : cvlist(ie)%vert(k,i,j) = cvlist(ie)%vert(k2,i,j)
1858 : end if
1859 : end do
1860 : end do
1861 : end do
1862 : end do
1863 : ! Release memory
1864 0 : if(hybrid%masterthread) then
1865 0 : call freeedgebuffer(edge1)
1866 : end if
1867 :
1868 0 : initialized=.true.
1869 0 : end subroutine InitControlVolumes_gll
1870 :
1871 0 : subroutine construct_cv_gll(elem,hybrid,nets,nete)
1872 : !
1873 : ! construct global dual grid from local element dual grid cvlist(ie)%cartp_dual(:,:)
1874 : ! all control volumes will be squares or triangles (at cube corners)
1875 : !
1876 : ! 10/2009: added option to make hexagon control volumes at cube edges and corners
1877 : !
1878 0 : use bndry_mod, only: bndry_exchange
1879 : use element_mod, only: element_t
1880 : use hybrid_mod, only: hybrid_t
1881 : use edge_mod, only: edgeVpack, edgeVunpack, edgeVunpackVert
1882 :
1883 : type(element_t), intent(in), target :: elem(:)
1884 : type(hybrid_t), intent(in) :: hybrid
1885 : integer, intent(in) :: nets,nete
1886 : ! local
1887 : integer :: i,j,k,ie,kptr,nvert,ie2
1888 : logical :: corner
1889 : real(r8) :: test(np,np,1),vertpack(np,np,3),rvert
1890 : type(cartesian2d_t) :: vert(4)
1891 : type(cartesian2d_t) :: cartp_nm1(0:np,0:np)
1892 :
1893 0 : test(:,:,:) = 0
1894 :
1895 0 : do ie=nets,nete
1896 : ! now construct the dual grid
1897 :
1898 0 : cartp_nm1 = cvlist(ie)%cartp_dual
1899 :
1900 0 : do j=1,np
1901 0 : do i=1,np
1902 0 : cvlist(ie)%vert(:,i,j)%x = 0_r8
1903 0 : cvlist(ie)%vert(:,i,j)%y = 0_r8
1904 0 : cvlist(ie)%vert(:,i,j)%z = 0_r8
1905 : end do
1906 : end do
1907 :
1908 : ! interior
1909 :
1910 0 : do j=2,np-1
1911 0 : do i=2,np-1
1912 :
1913 : ! internal vertex on Cubed sphere
1914 : ! Here is the order:
1915 : !
1916 : ! 4NW <- 3NE
1917 : ! | ^
1918 : ! v |
1919 : ! 1SW -> 2SE
1920 0 : vert(1)%x = cartp_nm1(i-1,j-1)%x
1921 0 : vert(1)%y = cartp_nm1(i-1,j-1)%y
1922 0 : vert(2)%x = cartp_nm1(i ,j-1)%x
1923 0 : vert(2)%y = cartp_nm1(i ,j-1)%y
1924 0 : vert(3)%x = cartp_nm1(i ,j )%x
1925 0 : vert(3)%y = cartp_nm1(i ,j )%y
1926 0 : vert(4)%x = cartp_nm1(i-1,j )%x
1927 0 : vert(4)%y = cartp_nm1(i-1,j )%y
1928 :
1929 0 : do k=1,4
1930 0 : cvlist(ie)%vert(k,i,j) = cubedsphere2cart(vert(k),elem(ie)%FaceNum)
1931 : end do
1932 0 : cvlist(ie)%nvert(i,j) = 4
1933 :
1934 : end do
1935 : end do
1936 :
1937 : ! Compute everything on the edges and then sum
1938 0 : do i=2,np-1
1939 0 : j=1
1940 : !
1941 : ! 4NW <- 3NE
1942 : ! | ^
1943 : ! v |
1944 : ! 1SW -> 2SE
1945 : !
1946 : !
1947 : ! only pack top two nodes.
1948 : ! leave other data zero, filled in by edgeexchange
1949 0 : cvlist(ie)%vert(4,i,j)%x = cvlist(ie)%vert(1,i,j+1)%x
1950 0 : cvlist(ie)%vert(4,i,j)%y = cvlist(ie)%vert(1,i,j+1)%y
1951 0 : cvlist(ie)%vert(4,i,j)%z = cvlist(ie)%vert(1,i,j+1)%z
1952 0 : cvlist(ie)%vert(3,i,j)%x = cvlist(ie)%vert(2,i,j+1)%x
1953 0 : cvlist(ie)%vert(3,i,j)%y = cvlist(ie)%vert(2,i,j+1)%y
1954 0 : cvlist(ie)%vert(3,i,j)%z = cvlist(ie)%vert(2,i,j+1)%z
1955 :
1956 :
1957 0 : j=np
1958 :
1959 0 : cvlist(ie)%vert(1,i,j)%x = cvlist(ie)%vert(4,i,j-1)%x
1960 0 : cvlist(ie)%vert(1,i,j)%y = cvlist(ie)%vert(4,i,j-1)%y
1961 0 : cvlist(ie)%vert(1,i,j)%z = cvlist(ie)%vert(4,i,j-1)%z
1962 0 : cvlist(ie)%vert(2,i,j)%x = cvlist(ie)%vert(3,i,j-1)%x
1963 0 : cvlist(ie)%vert(2,i,j)%y = cvlist(ie)%vert(3,i,j-1)%y
1964 0 : cvlist(ie)%vert(2,i,j)%z = cvlist(ie)%vert(3,i,j-1)%z
1965 :
1966 : end do
1967 :
1968 0 : do j=2,np-1
1969 0 : i=1
1970 :
1971 0 : cvlist(ie)%vert(2,i,j)%x = cvlist(ie)%vert(1,i+1,j)%x
1972 0 : cvlist(ie)%vert(2,i,j)%y = cvlist(ie)%vert(1,i+1,j)%y
1973 0 : cvlist(ie)%vert(2,i,j)%z = cvlist(ie)%vert(1,i+1,j)%z
1974 0 : cvlist(ie)%vert(3,i,j)%x = cvlist(ie)%vert(4,i+1,j)%x
1975 0 : cvlist(ie)%vert(3,i,j)%y = cvlist(ie)%vert(4,i+1,j)%y
1976 0 : cvlist(ie)%vert(3,i,j)%z = cvlist(ie)%vert(4,i+1,j)%z
1977 :
1978 0 : i=np
1979 :
1980 0 : cvlist(ie)%vert(4,i,j)%x = cvlist(ie)%vert(3,i-1,j)%x
1981 0 : cvlist(ie)%vert(4,i,j)%y = cvlist(ie)%vert(3,i-1,j)%y
1982 0 : cvlist(ie)%vert(4,i,j)%z = cvlist(ie)%vert(3,i-1,j)%z
1983 0 : cvlist(ie)%vert(1,i,j)%x = cvlist(ie)%vert(2,i-1,j)%x
1984 0 : cvlist(ie)%vert(1,i,j)%y = cvlist(ie)%vert(2,i-1,j)%y
1985 0 : cvlist(ie)%vert(1,i,j)%z = cvlist(ie)%vert(2,i-1,j)%z
1986 :
1987 : end do
1988 :
1989 : ! Corners
1990 : ! SW
1991 0 : cvlist(ie)%vert(3,1,1)%x = cvlist(ie)%vert(1,2,2)%x
1992 0 : cvlist(ie)%vert(3,1,1)%y = cvlist(ie)%vert(1,2,2)%y
1993 0 : cvlist(ie)%vert(3,1,1)%z = cvlist(ie)%vert(1,2,2)%z
1994 :
1995 : ! SE
1996 0 : cvlist(ie)%vert(4,np,1)%x = cvlist(ie)%vert(2,np-1,2)%x
1997 0 : cvlist(ie)%vert(4,np,1)%y = cvlist(ie)%vert(2,np-1,2)%y
1998 0 : cvlist(ie)%vert(4,np,1)%z = cvlist(ie)%vert(2,np-1,2)%z
1999 :
2000 : ! NE
2001 0 : cvlist(ie)%vert(1,np,np)%x = cvlist(ie)%vert(3,np-1,np-1)%x
2002 0 : cvlist(ie)%vert(1,np,np)%y = cvlist(ie)%vert(3,np-1,np-1)%y
2003 0 : cvlist(ie)%vert(1,np,np)%z = cvlist(ie)%vert(3,np-1,np-1)%z
2004 :
2005 : ! NW
2006 0 : cvlist(ie)%vert(2,1,np)%x = cvlist(ie)%vert(4,2,np-1)%x
2007 0 : cvlist(ie)%vert(2,1,np)%y = cvlist(ie)%vert(4,2,np-1)%y
2008 0 : cvlist(ie)%vert(2,1,np)%z = cvlist(ie)%vert(4,2,np-1)%z
2009 :
2010 0 : kptr=0
2011 0 : test(:,:,1) = cvlist(ie)%vol(:,:)
2012 0 : call edgeVpack(edge1,test(1,1,1),1,kptr,ie)
2013 :
2014 0 : cvlist(ie)%invvol(:,:) = cvlist(ie)%vol(:,:)
2015 :
2016 : end do ! loop over NE
2017 :
2018 0 : call bndry_exchange(hybrid,edge1,location='construct_cv_gll #1')
2019 :
2020 0 : do ie=nets,nete
2021 0 : kptr=0
2022 0 : call edgeVunpack(edge1, cvlist(ie)%invvol(1,1),1, kptr, ie)
2023 0 : cvlist(ie)%totvol(:,:)=cvlist(ie)%invvol(:,:)
2024 0 : cvlist(ie)%invvol(:,:)=1.0_r8/cvlist(ie)%invvol(:,:)
2025 : end do
2026 :
2027 : ! Create the polygon at the edges of the element
2028 :
2029 :
2030 : if(.NOT.(MODULO(np,2)==0)) then
2031 : call endrun("surfaces_mod: NV odd not implemented")
2032 : end if
2033 0 : vertpack = 0
2034 0 : do ie=nets,nete
2035 : ! Special messed up copy
2036 : !
2037 : !ASC should be replaced by a edgepack
2038 : ! S+N
2039 0 : do i=1,np/2
2040 0 : j=1
2041 0 : vertpack(i,j,1) = cvlist(ie)%vert(3,i,j)%x
2042 0 : vertpack(i,j,2) = cvlist(ie)%vert(3,i,j)%y
2043 0 : vertpack(i,j,3) = cvlist(ie)%vert(3,i,j)%z
2044 0 : j=np
2045 0 : vertpack(i,j,1) = cvlist(ie)%vert(2,i,j)%x
2046 0 : vertpack(i,j,2) = cvlist(ie)%vert(2,i,j)%y
2047 0 : vertpack(i,j,3) = cvlist(ie)%vert(2,i,j)%z
2048 : end do
2049 :
2050 0 : do i=np/2+1,np
2051 0 : j=1
2052 0 : vertpack(i,j,1) = cvlist(ie)%vert(4,i,j)%x
2053 0 : vertpack(i,j,2) = cvlist(ie)%vert(4,i,j)%y
2054 0 : vertpack(i,j,3) = cvlist(ie)%vert(4,i,j)%z
2055 0 : j=np
2056 0 : vertpack(i,j,1) = cvlist(ie)%vert(1,i,j)%x
2057 0 : vertpack(i,j,2) = cvlist(ie)%vert(1,i,j)%y
2058 0 : vertpack(i,j,3) = cvlist(ie)%vert(1,i,j)%z
2059 : end do
2060 :
2061 : ! E+W
2062 0 : do j=2,np/2
2063 0 : i=1
2064 0 : vertpack(i,j,1) = cvlist(ie)%vert(3,i,j)%x
2065 0 : vertpack(i,j,2) = cvlist(ie)%vert(3,i,j)%y
2066 0 : vertpack(i,j,3) = cvlist(ie)%vert(3,i,j)%z
2067 0 : i=np
2068 0 : vertpack(i,j,1) = cvlist(ie)%vert(4,i,j)%x
2069 0 : vertpack(i,j,2) = cvlist(ie)%vert(4,i,j)%y
2070 0 : vertpack(i,j,3) = cvlist(ie)%vert(4,i,j)%z
2071 : end do
2072 :
2073 0 : do j=np/2+1,np-1
2074 0 : i=1
2075 0 : vertpack(i,j,1) = cvlist(ie)%vert(2,i,j)%x
2076 0 : vertpack(i,j,2) = cvlist(ie)%vert(2,i,j)%y
2077 0 : vertpack(i,j,3) = cvlist(ie)%vert(2,i,j)%z
2078 0 : i=np
2079 0 : vertpack(i,j,1) = cvlist(ie)%vert(1,i,j)%x
2080 0 : vertpack(i,j,2) = cvlist(ie)%vert(1,i,j)%y
2081 0 : vertpack(i,j,3) = cvlist(ie)%vert(1,i,j)%z
2082 : end do
2083 :
2084 0 : do j=2,np-1
2085 0 : do i=2,np-1
2086 0 : vertpack(i,j,1) =0_r8
2087 0 : vertpack(i,j,2) =0_r8
2088 0 : vertpack(i,j,3) =0_r8
2089 : end do
2090 : end do
2091 :
2092 0 : kptr=0
2093 0 : call edgeVpack(edge1,vertpack,3,kptr,ie)
2094 : end do
2095 :
2096 0 : call bndry_exchange(hybrid,edge1,location='construct_cv_gll #2')
2097 :
2098 0 : do ie=nets,nete
2099 0 : kptr=0
2100 0 : call edgeVunpackVert(edge1, cvlist(ie)%vert,ie)
2101 : ! Count and orient vert array
2102 : ! nvert is an integer: -4,-3,3,4
2103 : ! Positive: 1->2->3->4 is counter clockwise on the sphere
2104 : ! Negative: clockwise orientation
2105 0 : do j=1,np
2106 0 : do i=1,np
2107 0 : nvert=0
2108 0 : do k=1,4
2109 0 : rvert = cvlist(ie)%vert(k,i,j)%x**2+cvlist(ie)%vert(k,i,j)%y**2+cvlist(ie)%vert(k,i,j)%z**2
2110 0 : if(rvert>0.9_r8)nvert=nvert+1
2111 : end do
2112 0 : if(.NOT.Orientation(cvlist(ie)%vert(:,i,j),elem(ie)%FaceNum))nvert=-nvert
2113 0 : cvlist(ie)%nvert(i,j) = nvert
2114 : corner = ( ((i==1) .and. (j==1)) .or. &
2115 : ((i==1) .and. (j==np)) .or. &
2116 : ((i==np) .and. (j==1)) .or. &
2117 0 : ((i==np) .and. (j==np)) )
2118 0 : if (abs(nvert)/=4) then
2119 0 : if (abs(nvert)/=3) then
2120 0 : write(iulog, *) 'i,j,nvert=',i,j,nvert
2121 0 : call shr_sys_flush(iulog)
2122 0 : call endrun('construct_cv_gll: bad value of nvert')
2123 : end if
2124 0 : if (.not. corner) then
2125 0 : write(iulog, *) 'non-corner node with only 3 verticies'
2126 0 : write(iulog, *) 'ie,i,j,nvert,corner=',ie,i,j,nvert,corner
2127 0 : write(iulog, *) cvlist(ie)%vert(1,i,j)%x
2128 0 : write(iulog, *) cvlist(ie)%vert(2,i,j)%x
2129 0 : write(iulog, *) cvlist(ie)%vert(3,i,j)%x
2130 0 : write(iulog, *) cvlist(ie)%vert(4,i,j)%x
2131 : !write(iulog, *) 'dual:'
2132 : !do ie2=nets,nete
2133 : ! write(iulog, *) ie2,maxval(cvlist(ie2)%cartp_dual(:,:)%x)
2134 : ! write(iulog, *) ie2,maxval(cvlist(ie2)%cartp_dual(:,:)%y)
2135 : !end do
2136 0 : call shr_sys_flush(iulog)
2137 0 : call endrun('construct_cv_gll: corner point should have nvert=3')
2138 : end if
2139 : ! nvert=3. we are at a cube corner. One of the control volume
2140 : ! nodes from the 'missing' corner element should be all zeros:
2141 0 : if (cvlist(ie)%vert(1,i,j)%x==0) then
2142 : ! ok
2143 0 : else if (cvlist(ie)%vert(2,i,j)%x==0) then
2144 : ! ok
2145 0 : else if (cvlist(ie)%vert(3,i,j)%x==0) then
2146 : ! ok
2147 0 : else if (cvlist(ie)%vert(4,i,j)%x==0) then
2148 : ! ok
2149 : else
2150 0 : write(iulog, *) 'cube corner node with 4 neighbors'
2151 0 : write(iulog, *) 'ie,i,j,nvert,corner=',ie,i,j,nvert,corner
2152 0 : write(iulog, *) cvlist(ie)%vert(1,i,j)%x
2153 0 : write(iulog, *) cvlist(ie)%vert(2,i,j)%x
2154 0 : write(iulog, *) cvlist(ie)%vert(3,i,j)%x
2155 0 : write(iulog, *) cvlist(ie)%vert(4,i,j)%x
2156 0 : call shr_sys_flush(iulog)
2157 0 : call endrun('construct_cv_gll: control volume at cube corner should be a triangle')
2158 : end if
2159 :
2160 : end if
2161 : end do
2162 : end do
2163 : end do
2164 0 : end subroutine construct_cv_gll
2165 :
2166 0 : logical function Orientation(v, FaceNum) result(orient)
2167 :
2168 : type(cartesian3d_t), intent(in) :: v(3)
2169 : integer, intent(in) :: FaceNum
2170 :
2171 : type(cartesian3D_t) :: v12, v23
2172 : real(r8) :: test, cart(3,3)
2173 :
2174 0 : orient = .FALSE.
2175 :
2176 0 : if ((FaceNum == 5).OR.(FaceNum == 6)) then
2177 :
2178 0 : cart(1,1) = v(1)%x
2179 0 : cart(2,1) = v(1)%y
2180 0 : cart(3,1) = v(1)%z
2181 :
2182 0 : cart(1,2) = v(2)%x
2183 0 : cart(2,2) = v(2)%y
2184 0 : cart(3,2) = v(2)%z
2185 :
2186 0 : cart(1,3) = v(3)%x
2187 0 : cart(2,3) = v(3)%y
2188 0 : cart(3,3) = v(3)%z
2189 :
2190 0 : v12%x = cart(1,2) - cart(1,1)
2191 0 : v12%y = cart(2,2) - cart(2,1)
2192 0 : v12%z = cart(3,2) - cart(3,1)
2193 :
2194 0 : v23%x = cart(1,3) - cart(1,2)
2195 0 : v23%y = cart(2,3) - cart(2,2)
2196 0 : v23%z = cart(3,3) - cart(3,2)
2197 :
2198 : test = (v12%y*v23%z - v12%z*v23%y)*v12%x &
2199 : - (v12%x*v23%z - v12%z*v23%x)*v12%y &
2200 0 : + (v12%x*v23%y - v12%y*v23%x)*v12%z
2201 :
2202 0 : if (test > 0_r8)then
2203 0 : orient=.TRUE.
2204 : end if
2205 :
2206 : else
2207 : orient=.TRUE.
2208 : end if
2209 :
2210 0 : end function Orientation
2211 :
2212 0 : subroutine VerifVolumes(elem, hybrid,nets,nete)
2213 : use hybrid_mod, only: hybrid_t
2214 : use element_mod, only: element_t
2215 :
2216 : type(element_t), intent(in) :: elem(:)
2217 : integer, intent(in) :: nets,nete
2218 : type(hybrid_t), intent(in) :: hybrid
2219 :
2220 : real(r8) :: psum,ptot,Vol_tmp(1),corr,maxelem_variation
2221 : real(r8) :: vol(np,np,nets:nete),r,rmin,rmax,a1,a2,locmin,locmax,emin,emax,dx,dy
2222 : integer :: i,j,ie,kptr,face
2223 :
2224 0 : real(r8), pointer :: locvol(:,:)
2225 :
2226 0 : dx = pi/(2.0d0*dble(ne))
2227 0 : dy = dx
2228 :
2229 0 : if(.not. initialized) then
2230 0 : call endrun('VerifyVolumes: Attempt to use volumes prior to initializing')
2231 : end if
2232 0 : rmin=2
2233 0 : rmax=0
2234 0 : maxelem_variation=0
2235 0 : do ie=nets,nete
2236 0 : locvol => GetVolume(ie)
2237 0 : locmin = minval(locvol(:,:)*elem(ie)%rspheremp(:,:))
2238 0 : locmax = maxval(locvol(:,:)*elem(ie)%rspheremp(:,:))
2239 0 : rmin = min(rmin,locmin)
2240 0 : rmax = max(rmax,locmax)
2241 :
2242 0 : if (locmax > 1.01_r8) then
2243 0 : write(iulog, *) 'locmin(:,i)=',ie,locvol(1,1),1/elem(ie)%rspheremp(1,1)
2244 : end if
2245 :
2246 :
2247 0 : if (locmax-locmin > maxelem_variation) then
2248 0 : maxelem_variation = locmax-locmin
2249 0 : emin=locmin
2250 0 : emax=locmax
2251 : end if
2252 : end do
2253 0 : rmin = ParallelMin(rmin,hybrid)
2254 0 : rmax = ParallelMax(rmax,hybrid)
2255 0 : if(hybrid%masterthread) then
2256 0 : write(iulog,'(a,2e14.7)') "Min/max ratio between spherical and GLL area:",rmin,rmax
2257 : end if
2258 0 : if (maxelem_variation == ParallelMax(maxelem_variation,hybrid) ) then
2259 0 : write(iulog,'(a,2e14.7)') "Min/max ratio element with largest variation:",emin,emax
2260 : end if
2261 0 : call shr_sys_flush(iulog)
2262 :
2263 0 : rmin=2
2264 0 : rmax=0
2265 0 : do ie=nets,nete
2266 0 : a1 = SurfArea_dxdy(dx,dy,elem(ie)%cartp(1,1))
2267 0 : a2 = sum(elem(ie)%spheremp(:,:))
2268 0 : r=a1/a2
2269 0 : rmin = min(r,rmin)
2270 0 : rmax = max(r,rmax)
2271 : end do
2272 0 : rmin = ParallelMin(rmin,hybrid)
2273 0 : rmax = ParallelMax(rmax,hybrid)
2274 0 : if(hybrid%masterthread) then
2275 0 : write(*,'(a,2f12.9)') "Min/max ratio spherical and GLL element area:",rmin,rmax
2276 : end if
2277 :
2278 0 : do ie=nets,nete
2279 0 : global_shared_buf(ie,1:6) = 0.d0
2280 0 : face = elem(ie)%FaceNum
2281 0 : locvol => GetVolumeLocal(ie)
2282 0 : do j=1,np
2283 0 : do i=1,np
2284 0 : global_shared_buf(ie,face) = global_shared_buf(ie,face) + locvol(i,j)
2285 : end do
2286 : end do
2287 : end do
2288 0 : call wrap_repro_sum(nvars=6, comm=hybrid%par%comm)
2289 :
2290 0 : ptot=0_r8
2291 0 : do face=1,6
2292 0 : red_sum%buf(1) = global_shared_sum(face)
2293 0 : psum = red_sum%buf(1)
2294 :
2295 0 : ptot = ptot + psum
2296 :
2297 0 : if(hybrid%masterthread) then
2298 0 : write(*,'(a,i2,a,2e23.15)') "cube face:",face," : SURFACE FV =",&
2299 0 : 6_r8*psum/(4_r8 * pi), &
2300 0 : 6_r8*psum/(4_r8 * pi)-1
2301 : end if
2302 : end do
2303 :
2304 0 : if(hybrid%masterthread) then
2305 0 : write(iulog, *) "SURFACE FV (total)= ", ptot/(4_r8 * pi)
2306 : end if
2307 :
2308 0 : end subroutine VerifVolumes
2309 :
2310 0 : end module comp_gll_ctr_vol
|