LCOV - code coverage report
Current view: top level - dynamics/se/dycore - comp_ctr_vol_around_gll_pts.F90 (source / functions) Hit Total Coverage
Test: coverage.info Lines: 0 1188 0.0 %
Date: 2024-12-17 17:57:11 Functions: 0 21 0.0 %

          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

Generated by: LCOV version 1.14