module camse_utils

   implicit none

   private
   public :: debugLevel
   public :: camse_init, area_interp_init_camse

   real, parameter :: ZERO    = 0.
   real, parameter :: ZERO_8  = 0._8
   real, parameter :: ONEHALF = .5
   real, parameter :: ONEHALF_8 = .5_8
   real, parameter :: ONE     = 1.
   real, parameter :: ONE_8   = 1._8
   real, parameter :: FOUR    = 4.
   real, parameter :: FOUR_8  = 4._8
   real, parameter :: TEN_8   = 10._8
   real, parameter :: NINETY  = 90.
   real, parameter :: ONE80   = 180.
   real, parameter :: ONE80_8 = 180._8
   real, parameter :: THREE60 = 360.
   real, parameter :: THREE60_8 = 360._8
   real, parameter :: ROUNDOFF = 10.*epsilon(ROUNDOFF)

   integer :: maxVtx
   integer :: debugLevel = 0
   real(8) :: PI, D2R, R2D

   include 'netcdf.inc'

   CONTAINS

   subroutine camse_init( ncid, mdl_grd )

   use anthro_types, only : mdl_poly_type, model_grid_type
   use mdl_grd_def, only  : model_grid_definition
   use netcdf_utils

!-----------------------------------------------------------------------
!  dummy arguments
!-----------------------------------------------------------------------
   integer, intent(in) :: ncid         ! model netcdf file id
   type(model_grid_type), intent(inout) :: mdl_grd

!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------
   integer :: polyNdx, vtxNdx
   integer :: astat
   integer :: dimid, varid
   integer :: nMatchLon, nMatchLat, nMatchVtx
   integer :: m, mp1
   integer :: nVtx, nPolygons
   integer :: CapNdx
   integer :: minmaxLat(2)
   integer :: Quadcnt(4)
   integer, allocatable :: Quad(:)
   real    :: MatchVtxLon, MatchVtxLat
   real, allocatable :: wrk(:,:)
   real, allocatable :: wrkLon(:), wrkLat(:), xprod(:)
   real, allocatable :: CenterLon(:), CenterLat(:), polyArea(:)
   real, allocatable :: VertexLon(:,:), VertexLat(:,:)
   character(len=64) :: mess
   character(len=16) :: mdlType

   mdl_grd%sinceDate = '1750-01-01'
   if( mdl_grd%mdl_is_MPAS ) then
     mdlType = 'MPAS'
   elseif( mdl_grd%mdl_is_CAMSE ) then
     mdlType = 'CAMSE'
   endif

   call model_grid_definition( ncid, mdlType, nPolygons, maxVtx, &
                               CenterLon, CenterLat, polyArea, &
                               VertexLon, VertexLat )

   mdl_grd%nPolygons  = nPolygons
   mdl_grd%maxPolyvtx = maxVtx

!-----------------------------------------------------------------------
!  allocate mdl_poly_type
!-----------------------------------------------------------------------
   allocate( mdl_grd%mdl_poly(nPolygons),stat=astat )
   if( astat /= 0 ) then
     write(*,*) 'camse_init: Failed to allocate model polygons;  error = ',astat
     stop 'AllocErr'
   endif
   do polyNdx = 1,nPolygons
     mdl_grd%mdl_poly(polyNdx)%cntr_lon = CenterLon(polyNdx)
     mdl_grd%mdl_poly(polyNdx)%cntr_lat = CenterLat(polyNdx)
     mdl_grd%mdl_poly(polyNdx)%area     = polyArea(polyNdx)
     allocate( mdl_grd%mdl_poly(polyNdx)%vtx_lon(maxVtx), &
               mdl_grd%mdl_poly(polyNdx)%vtx_lat(maxVtx),stat=astat )
     if( astat /= 0 ) then
       write(*,*) 'camse_init: Failed to allocate polygon vtx lon,lat;  error = ',astat
       stop 'AllocErr'
     endif
     mdl_grd%mdl_poly(polyNdx)%vtx_lon(:) = VertexLon(:,polyNdx)
     mdl_grd%mdl_poly(polyNdx)%vtx_lat(:) = VertexLat(:,polyNdx)
   enddo

   deallocate( CenterLon,CenterLat,polyArea,VertexLon,VertexLat )

!-----------------------------------------------------------------------
!  delineate the "unique" model grid polygon vertices
!-----------------------------------------------------------------------
   do polyNdx = 1,nPolygons
     matchVtxLon = mdl_grd%mdl_poly(polyNdx)%vtx_lon(maxVtx)
     matchVtxLat = mdl_grd%mdl_poly(polyNdx)%vtx_lat(maxVtx)
     nMatchLon   = count( mdl_grd%mdl_poly(polyNdx)%vtx_lon(:) == matchVtxLon )
     nMatchLat   = count( mdl_grd%mdl_poly(polyNdx)%vtx_lat(:) == matchVtxLat )
     nMatchVtx   = min( nMatchLon,nMatchLat )
     mdl_grd%mdl_poly(polyNdx)%nVtx = maxVtx - nMatchVtx + 1
   enddo

   allocate( Quad(maxVtx),wrkLon(maxVtx),wrkLat(maxVtx),xprod(maxVtx),stat=astat )
   if( astat /= 0 ) then
     write(*,*) 'camse_init: Failed to allocate Quad,wrkLon;  error = ',astat
     stop 'AllocErr'
   endif

   CapNdx = 0
   mdl_grd%mdl_poly(1:nPolygons)%active = .true.
   mdl_grd%mdl_poly(1:nPolygons)%x180   = .false.
!-----------------------------------------------------------------------
!  mark polygons with edges that "cross" 180 longitude
!-----------------------------------------------------------------------
poly_loop: &
   do polyNdx = 1,nPolygons
     nVtx = mdl_grd%mdl_poly(polyNdx)%nVtx
     if( debugLevel >= 300 ) then
       wrkLon(1:nVtx) = mdl_grd%mdl_poly(polyNdx)%vtx_lon(1:nVtx) - mdl_grd%mdl_poly(polyNdx)%cntr_lon
       wrkLat(1:nVtx) = mdl_grd%mdl_poly(polyNdx)%vtx_lat(1:nVtx) - mdl_grd%mdl_poly(polyNdx)%cntr_lat
       do m = 1,nVtx
         if( m /= nVtx ) then
           mp1 = m + 1
         else
           mp1 = 1
         endif
         xprod(m) = wrkLon(m)*wrkLat(mp1) - wrkLon(mp1)*wrkLat(m)
       enddo
       if( any(xprod(1:nVtx) < ZERO) ) then
         write(*,'(''camse_init: negative cross prod for polygon,neg,pos '',i6,2i3)') &
           polyndx,count(xprod(1:nVtx) < ZERO),count(xprod(1:nVtx) > ZERO)
         write(*,'(''camse_init: poly lons = '',1p10g15.7)') mdl_grd%mdl_poly(polyndx)%vtx_lon(1:nVtx)
       endif
     endif
!-----------------------------------------------------------------------
!  put all mdl polygon vertices and centers into range (-180,180)
!-----------------------------------------------------------------------
     if( mdl_grd%mdl_poly(polyNdx)%cntr_lon > ONE80 ) then
       mdl_grd%mdl_poly(polyNdx)%cntr_lon = &
          mdl_grd%mdl_poly(polyNdx)%cntr_lon - THREE60
     endif
     where( mdl_grd%mdl_poly(polyNdx)%vtx_lon(1:nVtx) > ONE80 )
       mdl_grd%mdl_poly(polyNdx)%vtx_lon(1:nVtx) = &
          mdl_grd%mdl_poly(polyNdx)%vtx_lon(1:nVtx) - THREE60
     endwhere
     wrkLon(1:nVtx) = mdl_grd%mdl_poly(polyNdx)%vtx_lon(1:nVtx)
!-----------------------------------------------------------------------
!  place vtx in "quadrant" by longitude
!-----------------------------------------------------------------------
     Quadcnt(:) = 0
vtx_loop: &
     do vtxNdx = 1,nVtx
       if( ZERO <= wrkLon(vtxNdx) .and. wrkLon(vtxNdx) < NINETY ) then
         Quad(vtxNdx) = 1
         Quadcnt(1) = Quadcnt(1) + 1
       elseif( NINETY <= wrkLon(vtxNdx) .and. wrkLon(vtxNdx) < ONE80 ) then
         Quad(vtxNdx) = 2
         Quadcnt(2) = Quadcnt(2) + 1
       elseif( -ONE80 <= wrkLon(vtxNdx) .and. wrkLon(vtxNdx) < -NINETY ) then
         Quad(vtxNdx) = 3
         Quadcnt(3) = Quadcnt(3) + 1
       elseif( -NINETY <= wrkLon(vtxNdx) .and. wrkLon(vtxNdx) < ZERO ) then
         Quad(vtxNdx) = 4
         Quadcnt(4) = Quadcnt(4) + 1
       endif
     enddo vtx_loop
     if( count( Quadcnt(:) /= 0 ) > 2 ) then
       mdl_grd%mdl_poly(polyNdx)%active = .false.
       CapNdx = CapNdx + 1
       if( CapNdx < 3 ) then
         mdl_grd%PolarCapNdx(CapNdx) = polyNdx
       else
         write(*,'(''camse_init: more than two polar cap model grid cells; terminating'')')
         stop 'Runtime-Err'
       endif
     elseif( minval(wrkLon(1:nvtx)) * maxval(wrkLon(1:nVtx)) < 0. ) then
       if( any(Quad(:nVtx) == 2) .and. any(Quad(:nVtx) == 3) ) then
         mdl_grd%mdl_poly(polyNdx)%x180 = .true.
       endif
     endif
   enddo poly_loop

   deallocate( Quad,wrkLon,wrkLat,xprod )

!-----------------------------------------------------------------------
!  diagnostics
!-----------------------------------------------------------------------
   write(*,*) ' '
   write(*,'(''camse_init: There are '',i6,'' total mdl grid polygons'')') nPolygons
   write(*,'(''camse_init: There are '',i2,'' polar cap polygons'')') &
     count( .not. mdl_grd%mdl_poly(:)%active )
   write(*,'(''camse_init: There are '',i6,'' total mdl polygons crossing 180 longitude'')') &
     count( mdl_grd%mdl_poly(:)%x180 )
   minmaxLat(1:1) = minloc( abs(mdl_grd%mdl_poly(:)%cntr_lat),mask=.not. mdl_grd%mdl_poly(:)%x180 )
   minmaxLat(2:2) = maxloc( abs(mdl_grd%mdl_poly(:)%cntr_lat),mask=.not. mdl_grd%mdl_poly(:)%x180 )
   write(*,*) ' '
   write(*,'(''camse_init: masked polygon ('',i6,'') cntr_lat min = '',1pg15.7)') &
     minmaxLat(1),mdl_grd%mdl_poly(minmaxLat(1))%cntr_lat
   write(*,*) ' '
   write(*,'(''camse_init: masked polygon ('',i6,'') cntr_lat max = '',1pg15.7)') &
     minmaxLat(2),mdl_grd%mdl_poly(minmaxLat(2))%cntr_lat

   PI = FOUR*atan(ONE)
   D2R = PI/ONE80_8
   R2D = ONE80_8/PI

   end subroutine camse_init

   subroutine area_interp_init_camse( proj, mdlGrid, d2mMap, mdl_area_type, &
                                      is_CMIP6, is_EPA, srcCellArea, diag_level )
!-----------------------------------------------------------------------
!  initialize area interpolation for camse model
!-----------------------------------------------------------------------

   use anthro_types, only : mdl_poly_type, model_grid_type
   use mapper_types, only : grid_type, area_type, proj_info
   use area_mapper,  only : poly_area
   use constants_module, only : earth_radius_m, pi
!-----------------------------------------------------------------------
!  dummy arguments
!-----------------------------------------------------------------------
   integer,               intent(in)    :: diag_level
   real,                  intent(inout) :: srcCellArea(:,:)
   logical,               intent(in)    :: is_EPA
   logical,               intent(in)    :: is_CMIP6
   type(model_grid_type), intent(inout) :: mdlGrid
   type(grid_type),       intent(in)    :: d2mMap
   type(area_type),       intent(inout) :: mdl_area_type(:)
   type(proj_info),       intent(in)    :: proj
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------
   type enclosingBox
     real :: minLon, maxLon
     real :: minLat, maxLat
   end type enclosingBox

   real, parameter :: msq2kmsq = 1.e-6

   integer :: astat
   integer :: polyNdx
   integer :: lonNdx, latNdx
   integer :: nVtx
   integer :: xcellCnt
   integer :: nData_n_Mdl, nMdl_n_Data
   integer :: ndVtx
   integer :: nX180
   integer :: divCnt
   integer :: tstCnt
   integer :: delPoly, delCnt
   integer :: m, nCells
   integer :: minmaxNdx(2)
   integer :: xLonNdx(d2mMap%nlons*d2mMap%nlats)
   integer :: xLatNdx(d2mMap%nlons*d2mMap%nlats)
   integer :: minmaxPolyNdx(2)
   real    :: swghts
   real    :: AreaAdj
   real    :: wghts(mdlGrid%nPolygons)
   real    :: dataCellArea(d2mMap%nlons,d2mMap%nlats)
   real(8) :: divergence
   real(8) :: xsectArea
   real(8) :: CellArea
   real(8) :: dtotArea, xtotArea, mtotArea
   real(8) :: mCellArea, mxCellArea, dCellArea
   real(8) :: dCellAreaMin, dCellAreaMax
   real(8) :: mCellAreaMin, mCellAreaMax
   real(8) :: TransLon
   real(8) :: dPoly_x(4), dPoly_y(4)
   real(8) :: dPoly_lam(4), dPoly_phi(4)
   real(8) :: mPoly_x(mdlGrid%maxPolyVtx), mPoly_y(mdlGrid%maxPolyVtx)
   real(8) :: mPoly_lam(mdlGrid%maxPolyVtx), mPoly_phi(mdlGrid%maxPolyVtx)
   logical :: mcross180
   logical :: diverges
   logical :: Mask(mdlGrid%nPolygons)
   logical :: Mask1(mdlGrid%nPolygons)
   logical :: ZeroLonMask(d2mMap%nlons,d2mMap%nlats)
   type(enclosingBox) :: mdlBox
   type(enclosingBox), pointer :: dataBox
   type(enclosingBox), target :: dataBoxes(d2mMap%nlons,d2mMap%nlats)
   type(enclosingBox), target :: xdataBoxes(d2mMap%nlons,d2mMap%nlats)

   real(8) :: polyintersectarea

!-----------------------------------------------------------------------
!  set enclosing "box" for overall data domain
!-----------------------------------------------------------------------
   dataBox => dataBoxes(1,1)
   dataBox%minLon = minval(d2mMap%xedge_2d(:,:) )
   dataBox%minLat = minval(d2mMap%yedge_2d(:,:) )
   dataBox%maxLon = maxval(d2mMap%xedge_2d(:,:) )
   dataBox%maxLat = maxval(d2mMap%yedge_2d(:,:) )

   write(*,*) ' '
   write(*,*) 'There are ',count(mdlGrid%mdl_poly(:)%active),' active polygons before box test'
   write(*,*) 'dataBox = ',dataBox

!-----------------------------------------------------------------------
!  weed out mdl polygons that are completely outside data grid
!-----------------------------------------------------------------------
!$OMP PARALLEL DO PRIVATE( polyNdx,nVtx,mdlBox )
   do polyNdx = 1,mdlGrid%nPolygons
     if( mdlGrid%mdl_poly(polyNdx)%active ) then
       nVtx = mdlGrid%mdl_poly(polyNdx)%nVtx
!-----------------------------------------------------------------------
!  set mdl cell enclosing "box"
!-----------------------------------------------------------------------
       mdlBox%minLon = minval(mdlGrid%mdl_poly(polyNdx)%vtx_lon(1:nVtx) )
       mdlBox%minLat = minval(mdlGrid%mdl_poly(polyNdx)%vtx_lat(1:nVtx) )
       mdlBox%maxLon = maxval(mdlGrid%mdl_poly(polyNdx)%vtx_lon(1:nVtx) )
       mdlBox%maxLat = maxval(mdlGrid%mdl_poly(polyNdx)%vtx_lat(1:nVtx) )
       if( .not. Boxes_overlap() ) then
         mdlGrid%mdl_poly(polyNdx)%active = .false.
       else
         mdlGrid%mdl_poly(polyNdx)%active = .true.
       endif
     endif
   enddo
!$OMP END PARALLEL DO

   write(*,*) ' '
   write(*,*) 'There are ',count(mdlGrid%mdl_poly(:)%active),' active polygons after  box test'

   Mask(:) = mdlGrid%mdl_poly(:)%active
   mdl_area_type(1:mdlGrid%nPolygons)%has_data = .false.
   mdl_area_type(1:mdlGrid%nPolygons)%active_dcell_cnt = 0
   mdl_area_type(1:mdlGrid%nPolygons)%interior_dcell_cnt = 0
   mdl_area_type(1:mdlGrid%nPolygons)%partial_dcell_cnt = 0

!-----------------------------------------------------------------------
!  a few more diagnostics
!-----------------------------------------------------------------------
   Mask1(:) = Mask(:) .and. mdlGrid%mdl_poly(:)%nVtx == 4
   write(*,*) ' '
   write(*,*) 'There are ',count(Mask1(:)),' active mdl quad polygons after  box test'
   if( count(Mask1(:)) > 0 ) then
     write(*,'(''Min quad area = '',1pg15.7,'' km^2'')') msq2kmsq*earth_radius_m**2*minval(mdlGrid%mdl_poly(:)%area,mask=Mask1)
     write(*,'(''Max quad area = '',1pg15.7,'' km^2'')') msq2kmsq*earth_radius_m**2*maxval(mdlGrid%mdl_poly(:)%area,mask=Mask1)
   endif
   minmaxNdx(1:1) = minloc(mdlGrid%mdl_poly(:)%area,mask=Mask)
   minmaxNdx(2:2) = maxloc(mdlGrid%mdl_poly(:)%area,mask=Mask)
   write(*,*) ' '
   write(*,'(''Min mdl poly area,ndx,#vtx = '',1pg15.7,i6,i3)') &
         msq2kmsq*earth_radius_m**2*mdlGrid%mdl_poly(minmaxNdx(1))%area, &
                                    minmaxNdx(1),mdlGrid%mdl_poly(minmaxNdx(1))%nVtx
   write(*,'(''Max mdl poly area,ndx,#vtx = '',1pg15.7,i6,i3)') &
         msq2kmsq*earth_radius_m**2*mdlGrid%mdl_poly(minmaxNdx(2))%area, &
                                    minmaxNdx(2),mdlGrid%mdl_poly(minmaxNdx(2))%nVtx

   dtotArea = 0._8
   dCellAreaMin = 1.e-3_8; dcellAreaMax = 0._8
   ndVtx    = 4
   nX180 = count( mdlGrid%mdl_poly(:)%x180 )
   ZeroLonMask(:,:) = .false.
!-----------------------------------------------------------------------
!  set dataset total area
!-----------------------------------------------------------------------
!$OMP PARALLEL PRIVATE( latNdx,lonNdx,dPoly_x,dPoly_lam,dPoly_y,dPoly_phi, &
!$OMP                   dCellArea ) SHARED( dtotArea,dCellAreaMin,dcellAreaMax )
!$OMP DO
   do latNdx = 1,d2mMap%nlats
     do lonNdx = 1,d2mMap%nlons
       dPoly_x(1) = real(d2mMap%xedge_2d(lonNdx,latNdx),8)
       dPoly_x(2) = real(d2mMap%xedge_2d(lonNdx+1,latNdx),8)
       dPoly_x(3) = real(d2mMap%xedge_2d(lonNdx+1,latNdx+1),8)
       dPoly_x(4) = real(d2mMap%xedge_2d(lonNdx,latNdx+1),8)
       dPoly_y(1) = real(d2mMap%yedge_2d(lonNdx,latNdx),8)
       dPoly_y(2) = real(d2mMap%yedge_2d(lonNdx+1,latNdx),8)
       dPoly_y(3) = real(d2mMap%yedge_2d(lonNdx+1,latNdx+1),8)
       dPoly_y(4) = real(d2mMap%yedge_2d(lonNdx,latNdx+1),8)
       dPoly_lam(1:4) = dPoly_x(1:4)*D2R
       dPoly_phi(1:4) = sin(dPoly_y(1:4)*D2R)
       dCellArea = poly_area( 4, dPoly_lam, dPoly_phi )
!$OMP CRITICAL
       dtotArea = dtotArea + dCellArea
       dCellAreaMin = min( dCellAreaMin,dCellArea )
       dCellAreaMax = max( dCellAreaMax,dCellArea )
!$OMP END CRITICAL
       dataCellArea(lonNdx,latNdx) = real( dCellArea,4 )
       if( is_CMIP6 ) then
         srcCellArea(lonNdx,latNdx)  = real( dCellArea,4 )
       endif
!-----------------------------------------------------------------------
!  "box" enclosing data cell
!-----------------------------------------------------------------------
       dataBoxes(lonNdx,latNdx)%minLon = real(minval(dPoly_x(:)),4)
       dataBoxes(lonNdx,latNdx)%minLat = real(minval(dPoly_y(:)),4)
       dataBoxes(lonNdx,latNdx)%maxLon = real(maxval(dPoly_x(:)),4)
       dataBoxes(lonNdx,latNdx)%maxLat = real(maxval(dPoly_y(:)),4)
       if( any(dPoly_x(:) < ZERO_8) .and. any(dPoly_x(:) >= ZERO_8) ) then
         ZeroLonMask(lonNdx,latNdx) = .true.
       endif
       if( nX180 > 0 ) then
         where( dPoly_x(1:4) < ZERO_8 )
           dPoly_x(1:4) = dPoly_x(1:4) + THREE60_8
         endwhere
         xdataBoxes(lonNdx,latNdx)%minLon = real(minval(dPoly_x(:)),4)
         xdataBoxes(lonNdx,latNdx)%maxLon = real(maxval(dPoly_x(:)),4)
         xdataBoxes(lonNdx,latNdx)%minLat = dataBoxes(lonNdx,latNdx)%minLat
         xdataBoxes(lonNdx,latNdx)%maxLat = dataBoxes(lonNdx,latNdx)%maxLat
       endif
     enddo
   enddo
!$OMP END DO
!$OMP END PARALLEL

   write(*,*) ' '
   write(*,'(''ZeroLon count = '',i6)') count(ZeroLonMask(:,:))
   dCellAreaMin = msq2kmsq*earth_radius_m**2*dCellAreaMin
   dCellAreaMax = msq2kmsq*earth_radius_m**2*dCellAreaMax
   write(*,*) 'Min,Max dCell area = ',dCellAreaMin,dCellAreaMax

   stop 'OpenMP TESTING'
!-----------------------------------------------------------------------
!  debug section
!-----------------------------------------------------------------------
   mtotArea = ZERO_8
   do polyNdx = 1,mdlGrid%nPolygons
     if( Mask(polyNdx) ) then
       mcross180 = mdlGrid%mdl_poly(polyNdx)%x180
       nVtx = mdlGrid%mdl_poly(polyNdx)%nVtx
       mPoly_y(1:nVtx) = real(mdlGrid%mdl_poly(polyNdx)%vtx_lat(1:nVtx),8)
       mPoly_phi(1:nVtx) = sin(mPoly_y(1:nVtx)*D2R)
       mPoly_x(1:nVtx) = real(mdlGrid%mdl_poly(polyNdx)%vtx_lon(1:nVtx),8)
       if( mcross180 ) then
         where( mPoly_x(1:nVtx) < ZERO_8 )
           mPoly_x(1:nVtx) = mPoly_x(1:nVtx) + THREE60_8
         endwhere
       endif
       mPoly_lam(1:nVtx) = mPoly_x(1:nVtx)*D2R
       mCellArea = poly_area( nVtx, mPoly_lam, mPoly_phi ) 
       mtotArea = mCellArea + mtotArea
     endif
   enddo

   xtotArea = ZERO_8
   mxCellArea = ZERO_8
   mCellAreaMin = 1.e-3_8; mcellAreaMax = ZERO_8
   nData_n_Mdl = 0 ; nMdl_n_Data = 0
   divCnt = 0
   delPoly = mdlGrid%nPolygons/5 ; delCnt = 0
   write(*,*) ' '
   write(*,'(''area_interp_init_camse: Completed '',f5.1,'' %'')',advance='no') delCnt*20.
!-----------------------------------------------------------------------
!  find intersection between mdl polygon and source polygon
!-----------------------------------------------------------------------
!$OMP PARALLEL PRIVATE( polyNdx,latNdx,lonNdx,dPoly_x,dPoly_lam,dPoly_y,dPoly_phi, &
!$OMP                   CellArea,mcross180,nVtx,mPoly_x,mPoly_lam,mPoly_y,mPoly_phi, &
!$OMP                   mCellArea,xCellCnt,tstCnt,dataBox,xsectArea,nData_n_Mdl, &
!$OMP                   nMdl_n_Data,swghts,divergence,diverges,TransLon,AreaAdj,wghts, &
!$OMP                   xLonNdx,xLatNdx ) &
!$OMP          SHARED( xtotArea,divcnt,mCellAreaMin,mCellAreaMax )
!$OMP DO
poly_loop: &
   do polyNdx = 1,mdlGrid%nPolygons
mdl_cell_is_active: &
     if( Mask(polyNdx) ) then
       if( mod( polyNdx,delPoly ) == 0 ) then
         delCnt = delCnt + 1
         write(*,'(tl5,f5.1,'' %'')',advance='no') delCnt*20.
       endif
       CellArea = ZERO_8
       mcross180 = mdlGrid%mdl_poly(polyNdx)%x180
       nVtx = mdlGrid%mdl_poly(polyNdx)%nVtx
       mPoly_y(1:nVtx) = real(mdlGrid%mdl_poly(polyNdx)%vtx_lat(1:nVtx),8)
       mPoly_phi(1:nVtx) = sin(mPoly_y(1:nVtx)*D2R)
       mPoly_x(1:nVtx) = real(mdlGrid%mdl_poly(polyNdx)%vtx_lon(1:nVtx),8)
       if( mcross180 ) then
         where( mPoly_x(1:nVtx) < ZERO_8 )
           mPoly_x(1:nVtx) = mPoly_x(1:nVtx) + THREE60_8
         endwhere
       endif
!-----------------------------------------------------------------------
!  "box" enclosing mdl cell
!-----------------------------------------------------------------------
       mdlBox%minLon = real(minval(mPoly_x(1:nVtx)),4)
       mdlBox%maxLon = real(maxval(mPoly_x(1:nVtx)),4)
       mdlBox%minLat = real(minval(mPoly_y(1:nVtx)),4)
       mdlBox%maxLat = real(maxval(mPoly_y(1:nVtx)),4)

       mPoly_lam(1:nVtx) = mPoly_x(1:nVtx)*D2R
       mCellArea = poly_area( nVtx, mPoly_lam, mPoly_phi ) 
!$OMP CRITICAL
       mCellAreaMin = min( mCellAreaMin,mCellArea )
       mCellAreaMax = max( mCellAreaMax,mCellArea )
!$OMP END CRITICAL

Tst_loop: &
       do tstCnt = 1,2
       xcellCnt = 0
dlat_loop: &
       do latNdx = 1,d2mMap%nlats
dlon_loop: &
         do lonNdx = 1,d2mMap%nlons
!-----------------------------------------------------------------------
!  "box" enclosing data cell
!-----------------------------------------------------------------------
           if( .not. mcross180 ) then
             dataBox => dataBoxes(lonNdx,latNdx)
           else
             if( .not. ZeroLonMask(lonNdx,latNdx) ) then
               dataBox => xdataBoxes(lonNdx,latNdx)
             else
               cycle dlon_loop
             endif
           endif
!-----------------------------------------------------------------------
!  set xsection parameters if mdl, data enclosing "boxes" overlap
!-----------------------------------------------------------------------
overlap:   if( Boxes_overlap() ) then
             dPoly_x(1) = real(d2mMap%xedge_2d(lonNdx,latNdx),8)
             dPoly_x(2) = real(d2mMap%xedge_2d(lonNdx+1,latNdx),8)
             dPoly_x(3) = real(d2mMap%xedge_2d(lonNdx+1,latNdx+1),8)
             dPoly_x(4) = real(d2mMap%xedge_2d(lonNdx,latNdx+1),8)
             if( tstCnt == 2 ) then
               dPoly_x(:) = dPoly_x(:) + TransLon
             endif
             if( mcross180 ) then
               where( dPoly_x(:) < ZERO_8 )
                 dPoly_x(:) = dPoly_x(:) + THREE60_8
               endwhere
             endif
             dPoly_lam(1:4) = dPoly_x(1:4)*D2R
             dPoly_y(1) = real(d2mMap%yedge_2d(lonNdx,latNdx),8)
             dPoly_y(2) = real(d2mMap%yedge_2d(lonNdx+1,latNdx),8)
             dPoly_y(3) = real(d2mMap%yedge_2d(lonNdx+1,latNdx+1),8)
             dPoly_y(4) = real(d2mMap%yedge_2d(lonNdx,latNdx+1),8)
             dPoly_phi(1:4) = sin(dPoly_y(1:4)*D2R)
             xsectArea = polyintersectarea( nData_n_Mdl, nMdl_n_Data, nVtx, ndVtx, mPoly_lam, mPoly_phi, dPoly_lam, dPoly_phi )
             if( xsectArea /= ZERO_8 ) then
               xcellCnt = xcellCnt + 1
               if( is_EPA ) then
                 wghts(xcellCnt) = real(xsectArea,4)/dataCellArea(lonNdx,latNdx)
               elseif( is_CMIP6 ) then
                 wghts(xcellCnt) = real(xsectArea,4)
               endif
               CellArea = CellArea + xsectArea
!              xtotArea = xtotArea + xsectArea
               xLonNdx(xcellCnt) = lonNdx ; xLatNdx(xcellCnt) = latNdx
               mdl_area_type(polyndx)%partial_dcell_cnt &
                 = mdl_area_type(polyndx)%partial_dcell_cnt + 1
             endif
           endif overlap
         enddo dlon_loop
       enddo dlat_loop

       mdl_area_type(polyNdx)%active_dcell_cnt = xcellCnt
has_xsects: &
       if( xcellCnt > 0 ) then
         if( is_CMIP6 ) then
           swghts     = sum(wghts(1:xcellCnt))
           divergence = 100._8*(real(swghts,8)/mCellArea - ONE_8)
           diverges   = abs(divergence) > 1.e-2_8
           if( diverges ) then
             if( tstCnt == 1 ) then
               TransLon = abs( real(mdlGrid%mdl_poly(polyNdx)%cntr_lon,8) ) + TEN_8
               mPoly_x(1:nVtx)   = mPoly_x(1:nVtx) + TransLon
               mPoly_lam(1:nVtx) = mPoly_x(1:nVtx)*D2R
               write(*,'(''('',i3.3,'') A_I_init_camse: WARNING mapping divergence ('',1pg9.4,''%) exceeds norm @ polyndx,xcellcnt,X180,tstcnt = '',2i6,2x,l1,i2,&
                         '' ; mdl lon,lat_cntr = '',1p2g15.7)') &
               divcnt,divergence,polyNdx,xcellCnt,mcross180,tstCnt,mdlGrid%mdl_poly(polyNdx)%cntr_lon ,mdlGrid%mdl_poly(polyNdx)%cntr_lat
               cycle Tst_loop
             elseif( tstCnt == 2 ) then
               divCnt = divCnt + 1
               if( divCnt == 1 ) then
                 write(*,*) ' '
               endif
               write(*,'(''('',i3.3,'') A_I_init_camse: WARNING mapping divergence ('',1pg9.4,''%) exceeds norm @ polyndx,xcellcnt,X180,tstcnt = '',2i6,2x,l1,i2,&
                         '' ; mdl lon,lat_cntr = '',1p2g15.7)') &
               divcnt,divergence,polyNdx,xcellCnt,mcross180,tstCnt,mdlGrid%mdl_poly(polyNdx)%cntr_lon ,mdlGrid%mdl_poly(polyNdx)%cntr_lat
             endif
!-----------------------------------------------------------------------
!  this, hopefully, is a temp kludge
!-----------------------------------------------------------------------
             write(*,'(''forcing area kludge @ polyNdx = '',i6)') polyNdx
             AreaAdj = real(mCellArea,4)/swghts
             wghts(1:xcellCnt) = AreaAdj*wghts(1:xcellCnt)
             swghts = sum(wghts(1:xcellCnt))
           endif
         endif
         allocate( mdl_area_type(polyNdx)%wght(xcellCnt), &
                   mdl_area_type(polyNdx)%dcell_lon_ndx(xcellCnt), &
                   mdl_area_type(polyNdx)%dcell_lat_ndx(xcellCnt),stat=astat )
         if( astat /= 0 ) then
           write(*,*) 'area_interp_init_camse: Failed to mdl_area_type array',astat
           Stop 'Alloc-ERR'
         endif
         mdl_area_type(polyNdx)%has_data = .true.
         mdl_area_type(polyNdx)%wght(1:xcellCnt) = wghts(1:xcellCnt)
         mdl_area_type(polyNdx)%dcell_lon_ndx(1:xcellCnt) = xLonNdx(1:xcellCnt)
         mdl_area_type(polyNdx)%dcell_lat_ndx(1:xcellCnt) = xLatNdx(1:xcellCnt)
         xtotArea = xtotArea + swghts
         exit Tst_loop
       endif has_xsects
       enddo Tst_loop
     endif mdl_cell_is_active
   enddo poly_loop
!$OMP END DO
!$OMP END PARALLEL

   write(*,*) ' '
   write(*,'(''area_interp_init_camse: '',i6,'' divergent cells; '',1pg15.7,''%'')') &
     divCnt,real(divCnt)/real(mdlGrid%nPolygons)

   write(*,*) ' '
   mCellAreaMin = msq2kmsq*earth_radius_m**2*mCellAreaMin
   mCellAreaMax = msq2kmsq*earth_radius_m**2*mCellAreaMax
   write(*,'(''Min,Max mCell area = '',1p2g15.7,'' km^2'')') mCellAreaMin,mCellAreaMax
   write(*,*) ' '
   write(*,'(''area_interp_init_camse: nData_n_Mdl,nMdl_n_Data = '',2i6)') nData_n_Mdl,nMdl_n_Data
   minmaxPolyNdx(1:1) = minloc( mdl_area_type(1:mdlgrid%nPolygons)%active_dcell_cnt, &
              mask=mdl_area_type(1:mdlgrid%nPolygons)%active_dcell_cnt > 0 )
   minmaxPolyNdx(2:2) = maxloc( mdl_area_type(1:mdlgrid%nPolygons)%active_dcell_cnt, &
              mask=mdl_area_type(1:mdlgrid%nPolygons)%active_dcell_cnt > 0 )
   write(*,'(''area_interp_init_camse: min dCell cnt = '',1pg15.7,'' @ polyNdx = '',i6)') &
      mdl_area_type(minmaxPolyNdx(1))%active_dcell_cnt,minmaxPolyNdx(1)
   write(*,'(''area_interp_init_camse: max dCell cnt = '',1pg15.7,'' @ polyNdx = '',i6)') &
      mdl_area_type(minmaxPolyNdx(2))%active_dcell_cnt,minmaxPolyNdx(2)
   write(*,'(''area_interp_init_camse: dtotArea,xtotArea,mxCellArea = '',1p3g15.7)') &
      dtotArea,xtotArea,mxCellArea

!-----------------------------------------------------------------------
!  a second sanity check on wghts
!-----------------------------------------------------------------------
!  if( debugLevel >= 200 .and. is_CMIP6 ) then
   if( is_CMIP6 ) then
     do latNdx = 1,d2mMap%nlats
       do lonNdx = 1,d2mMap%nlons
         dtotarea = ZERO_8
         nCells = 0
         do polyNdx = 1,mdlGrid%nPolygons
           if( Mask(polyNdx) ) then
             do m = 1,mdl_area_type(polyNdx)%active_dcell_cnt
               if( mdl_area_type(polyNdx)%dcell_lon_ndx(m) == lonNdx .and. &
                   mdl_area_type(polyNdx)%dcell_lat_ndx(m) == latNdx ) then
                 nCells = nCells + 1
                 dtotarea = dtotarea + real(mdl_area_type(polyNdx)%wght(m),8)
                 exit
               endif
             enddo
           endif
         enddo
         write(*,'(''area_interp_init_camse: dtotarea,srcarea,nCells,lon,latndx = '',1p2g15.7,3i4)') &
           dtotarea,srcCellArea(lonNdx,latNdx),nCells,lonndx,latndx
       enddo
     enddo
   endif

   stop 'DEBUGGING'

   CONTAINS

   function Boxes_overlap() result(Ovrlap)
     
     logical :: Ovrlap

     logical :: NoOvrlap

     NoOvrlap = (dataBox%minLon >= mdlBox%maxLon) .or. (dataBox%maxLon <= mdlBox%minLon)
     if( .not. NoOvrlap ) then
       NoOvrlap = (dataBox%minLat >= mdlBox%maxLat) .or. (dataBox%maxLat <= mdlBox%minLat)
     endif
     Ovrlap = .not. NoOvrlap

   end function Boxes_overlap

   end subroutine area_interp_init_camse

   end module camse_utils