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