MODULE remap
  INTEGER, PARAMETER ::                           &
       int_kind  = KIND(1),                       &
       real_kind = SELECTED_REAL_KIND(p=14,r=100),&
       dbl_kind  = selected_real_kind(13)        

  INTEGER :: nc,nhe

!  LOGICAL, PARAMETER:: ldbgr_r = .FALSE.
  LOGICAL :: ldbgr
  LOGICAL :: ldbg_global

  REAL(kind=real_kind), PARAMETER ::              &
       one = 1.0                       ,&
       aa  = 1.0                       ,&
       tiny= 1.0E-9  ,&
       bignum = 1.0E20
  REAL (KIND=dbl_kind), parameter :: fuzzy_width = 10.0*tiny  !CAM-SE add           

  contains


  subroutine compute_weights_cell(xcell_in,ycell_in,jx,jy,nreconstruction,xgno,ygno,&
       jx_min, jx_max, jy_min, jy_max,tmp,&
       ngauss,gauss_weights,abscissae,weights,weights_eul_index,jcollect,jmax_segments,&
       nc_in,nhe_in,nvertex,ldbg)

    implicit none
    integer (kind=int_kind)                  , intent(in):: nreconstruction, jx,jy,ngauss,jmax_segments
    real (kind=real_kind)   ,  dimension(0:nvertex+1)   :: xcell_in,ycell_in
!    real (kind=real_kind)   ,  dimension(0:5), intent(in):: xcell_in,ycell_in
    integer (kind=int_kind), intent(in) :: nc_in,nhe_in,nvertex
    logical, intent(in) :: ldbg
    !
    ! ipanel is just for debugging
    !
    integer (kind=int_kind), intent(in)               :: jx_min, jy_min, jx_max, jy_max
    real (kind=real_kind), dimension(-nhe_in:nc_in+2+nhe_in), intent(in) :: xgno
    real (kind=real_kind), dimension(-nhe_in:nc_in+2+nhe_in), intent(in) :: ygno
    !
    ! for Gaussian quadrature
    !
    real (kind=real_kind), dimension(ngauss), intent(in) :: gauss_weights, abscissae
    !
    ! boundaries of domain
    !
    real (kind=real_kind):: tmp
    !
    ! Number of Eulerian sub-cell integrals for the cell in question
    !
    integer (kind=int_kind), intent(out)       :: jcollect
    !
    ! local workspace
    !
    !
    ! max number of line segments is:
    !
    ! (number of longitudes)*(max average number of crossings per line segment = 3)*ncube*2
    !
    real (kind=real_kind)   ,  &
         dimension(jmax_segments,nreconstruction), intent(out) :: weights
    integer (kind=int_kind),  &
         dimension(jmax_segments,2), intent(out)      :: weights_eul_index
    
    real (kind=real_kind), dimension(0:3) :: x,y
    integer (kind=int_kind),dimension(0:5) :: jx_eul, jy_eul
    integer (kind=int_kind) :: jsegment,i
    !
    ! variables for registering crossings with Eulerian latitudes and longitudes
    !
    integer (kind=int_kind)  :: jcross_lat, iter
    !
    ! max. crossings per side is 2*nhe
    !
    real (kind=real_kind), &
         dimension(jmax_segments,2) :: r_cross_lat
    integer (kind=int_kind), &
         dimension(jmax_segments,2) :: cross_lat_eul_index
    real (kind=real_kind)   ,  dimension(1:nvertex) :: xcell,ycell

    real (kind=real_kind) :: eps

    ldbg_global = ldbg
    ldbgr = ldbg

    nc = nc_in
    nhe = nhe_in

    xcell = xcell_in(1:nvertex)
    ycell = ycell_in(1:nvertex)


    !
    ! this is to avoid ill-conditioning problems
    !
    eps = 1.0E-9

    jsegment = 0
    weights  = 0.0D0
    jcross_lat = 0
    !
    !**********************
    !
    ! Integrate cell sides
    !
    !**********************
       
    
    IF (jx<-nhe.OR.jx>nc+1+nhe.OR.jy<-nhe.OR.jy>nc+1+nhe) THEN
      WRITE(*,*) "jx,jy,-nhe,nc+1+nhe",jx,jy,-nhe,nc+1+nhe
      STOP
    END IF

    
    call side_integral(xcell,ycell,nvertex,jsegment,jmax_segments,&
         weights,weights_eul_index,nreconstruction,jx,jy,xgno,ygno,jx_min, jx_max, jy_min, jy_max,&
         ngauss,gauss_weights,abscissae,&
         jcross_lat,r_cross_lat,cross_lat_eul_index)
    
    !
    !**********************
    ! 
    ! Do inner integrals
    !
    !**********************
    !
    call compute_inner_line_integrals_lat_nonconvex(r_cross_lat,cross_lat_eul_index,&
         jcross_lat,jsegment,jmax_segments,xgno,jx_min, jx_max, jy_min, jy_max,&
         weights,weights_eul_index,&
         nreconstruction,ngauss,gauss_weights,abscissae)
    !
    ! collect line-segment that reside in the same Eulerian cell
    !
    if (jsegment>0) then
      call collect(weights,weights_eul_index,nreconstruction,jcollect,jsegment,jmax_segments)
      !
      ! DBG
      !
      tmp=0.0
      do i=1,jcollect     
        tmp=tmp+weights(i,1)
      enddo

      IF (abs(tmp)>0.01) THEN
        WRITE(*,*) "sum of weights too large",tmp
        stop
      END IF
      IF (tmp<-1.0E-9) THEN
        WRITE(*,*) "sum of weights is negative - negative area?",tmp,jx,jy
        !              ldbgr=.TRUE.
        ! Turn this off as part of physgrid work. This utility emits this
        ! message a few tens of times during the course of an ne30pg2 run, but
        ! there appears to be no issue in the final topo dataset. If we leave
        ! this stop statement, we can't get output from cube_to_target.
        !stop
      END IF
    else
      jcollect = 0
    end if
  end subroutine compute_weights_cell

  
  !
  !****************************************************************************
  !
  ! organize data and store it
  !
  !****************************************************************************
  !
  subroutine collect(weights,weights_eul_index,nreconstruction,jcollect,jsegment,jmax_segments)
    implicit none
    integer (kind=int_kind)                                 , intent(in)    :: nreconstruction
    real (kind=real_kind)   , dimension(jmax_segments,nreconstruction), intent(inout) :: weights
    integer (kind=int_kind), dimension(jmax_segments,2     ), intent(inout) :: weights_eul_index
    integer (kind=int_kind),                                  INTENT(OUT  ) :: jcollect
    integer (kind=int_kind),                                  INTENT(IN   ) :: jsegment,jmax_segments
    !
    ! local workspace
    !
    integer (kind=int_kind) :: imin, imax, jmin, jmax, i,j,k,h
    logical                 :: ltmp

    real (kind=real_kind)   , dimension(jmax_segments,nreconstruction) :: weights_out
    integer (kind=int_kind), dimension(jmax_segments,2     ) :: weights_eul_index_out

    weights_out           = 0.0D0
    weights_eul_index_out = -100

    imin = MINVAL(weights_eul_index(1:jsegment,1))
    imax = MAXVAL(weights_eul_index(1:jsegment,1))
    jmin = MINVAL(weights_eul_index(1:jsegment,2))
    jmax = MAXVAL(weights_eul_index(1:jsegment,2))

    ltmp = .FALSE.

    jcollect = 1

    do j=jmin,jmax
       do i=imin,imax
          do k=1,jsegment
             if (weights_eul_index(k,1)==i.AND.weights_eul_index(k,2)==j) then
                weights_out(jcollect,1:nreconstruction) = &
                     weights_out(jcollect,1:nreconstruction) + weights(k,1:nreconstruction)
                ltmp = .TRUE.
                h = k
             endif
          enddo
          if (ltmp) then
             weights_eul_index_out(jcollect,:) = weights_eul_index(h,:)
             jcollect = jcollect+1
          endif
          ltmp = .FALSE.
       enddo
    enddo
    jcollect = jcollect-1
    weights           = weights_out
    weights_eul_index = weights_eul_index_out
  end subroutine collect
  !
  !*****************************************************************************************
  !
  ! 
  !
  !*****************************************************************************************
  !
  subroutine compute_inner_line_integrals_lat(r_cross_lat,cross_lat_eul_index,&
       jcross_lat,jsegment,jmax_segments,xgno,jx_min,jx_max,jy_min, jy_max,weights,weights_eul_index,&
       nreconstruction,ngauss,gauss_weights,abscissae)!phl add jx_min etc.    
    implicit none
    !
    ! for Gaussian quadrature
    !
    real (kind=real_kind), dimension(ngauss), intent(in) :: gauss_weights, abscissae
    !
    ! variables for registering crossings with Eulerian latitudes and longitudes
    !
    integer (kind=int_kind),         intent(in):: jcross_lat, jmax_segments,nreconstruction,ngauss
    integer (kind=int_kind),         intent(inout):: jsegment
    !
    ! max. crossings per side is 2*nhe
    !
    real (kind=real_kind), &
         dimension(jmax_segments,2), intent(in):: r_cross_lat
    integer (kind=int_kind), &
         dimension(jmax_segments,2), intent(in):: cross_lat_eul_index
    integer (kind=int_kind), intent(in)            ::jx_min, jx_max, jy_min, jy_max
    real (kind=real_kind), dimension(-nhe:nc+2+nhe), intent(in) :: xgno
    real (kind=real_kind)   ,  &
         dimension(jmax_segments,nreconstruction), intent(inout) :: weights
    integer (kind=int_kind),  &
         dimension(jmax_segments,2), intent(inout) :: weights_eul_index
    real (kind=real_kind)   , dimension(nreconstruction) :: weights_tmp
    
    integer (kind=int_kind) :: imin, imax, jmin, jmax, i,j,k, isgn, h, eul_jx, eul_jy
    integer (kind=int_kind) :: idx_start_y,idx_end_y
    logical                 :: ltmp,lcontinue
    real (kind=real_kind), dimension(2)  :: rstart,rend,rend_tmp
    real (kind=real_kind), dimension(2)  :: xseg, yseg
5   FORMAT(10e14.6)
    
    
    if (jcross_lat>0) then
      do i=MINVAL(cross_lat_eul_index(1:jcross_lat,2)),MAXVAL(cross_lat_eul_index(1:jcross_lat,2))
        !
        ! find "first" crossing with Eulerian cell i
        !
        do k=1,jcross_lat
          if (cross_lat_eul_index(k,2)==i) exit
        enddo
        do j=k+1,jcross_lat
          !
          ! find "second" crossing with Eulerian cell i
          !
          if (cross_lat_eul_index(j,2)==i) then
            if (r_cross_lat(k,1)<r_cross_lat(j,1)) then
              rstart = r_cross_lat(k,1:2)
              rend   = r_cross_lat(j,1:2)
              imin   = cross_lat_eul_index(k,1)
              imax   = cross_lat_eul_index(j,1)
            else
              rstart = r_cross_lat(j,1:2)
              rend   = r_cross_lat(k,1:2)
              imin   = cross_lat_eul_index(j,1)
              imax   = cross_lat_eul_index(k,1)
            endif
            do h=imin,imax
              if (h==imax) then
                rend_tmp = rend
              else
                rend_tmp(1) = xgno(h+1)
                rend_tmp(2) = r_cross_lat(k,2)
              endif
              xseg(1) = rstart(1)
              xseg(2) = rend_tmp(1)
              yseg(1) = rstart(2)
              yseg(2) = rend_tmp(2)
              !                  call get_weights_exact(weights_tmp,xseg,yseg,nreconstruction)
              call get_weights_gauss(weights_tmp,&
                   xseg,yseg,nreconstruction,ngauss,gauss_weights,abscissae)
              
              
              if (i.LE.jy_max-1.AND.i.GE.jy_min.AND.h.LE.jx_max-1.AND.h.GE.jx_min) then
                jsegment=jsegment+1
                weights_eul_index(jsegment,1) = h 
                weights_eul_index(jsegment,2) = i
                weights(jsegment,1:nreconstruction) = -weights_tmp
                if (ldbg_global) then
                  OPEN(unit=40, file='inner_integral.dat',status='old',access='append')
                  WRITE(40,*) xseg(1),yseg(1)
                  WRITE(40,*) xseg(2),yseg(2)
                  WRITE(40,*) "  "
                  CLOSE(40)              
                end if                                
              endif
              
              !
              ! subtract the same weights on the "south" side of the line
              !
              if (i.LE.jy_max.AND.i.GE.jy_min+1.AND.h.LE.jx_max-1.AND.h.GE.jx_min) then
                !phl                   if (i.GE.2.AND.i.LE.nc+1.AND.h.LE.nc.AND.h.GE.1) then
                jsegment = jsegment+1
                weights_eul_index(jsegment,1) = h 
                weights_eul_index(jsegment,2) = i-1
                weights(jsegment,1:nreconstruction) = weights_tmp
              endif
              !
              ! prepare for next iteration
              !
              !                   if (abs(rend_tmp(1)-rend(1))<tiny) then
              !                      EXIT !are we done already?
              !                   else
              rstart = rend_tmp
              !                   endif
            enddo
          endif
        enddo
      enddo
    endif
  end subroutine compute_inner_line_integrals_lat

  subroutine compute_inner_line_integrals_lat_nonconvex(r_cross_lat,cross_lat_eul_index,&
       jcross_lat,jsegment,jmax_segments,xgno,jx_min,jx_max,jy_min, jy_max,weights,weights_eul_index,&
       nreconstruction,ngauss,gauss_weights,abscissae)!phl add jx_min etc.
    
    implicit none
    !
    ! for Gaussian quadrature
    !
    real (kind=real_kind), dimension(ngauss), intent(in) :: gauss_weights, abscissae
    !
    ! variables for registering crossings with Eulerian latitudes and longitudes
    !
    integer (kind=int_kind),         intent(in):: jcross_lat, jmax_segments,nreconstruction,ngauss
    integer (kind=int_kind),         intent(inout):: jsegment
    !
    ! max. crossings per side is 2*nhe
    !
    real (kind=real_kind), &
         dimension(jmax_segments,2), intent(in):: r_cross_lat
    integer (kind=int_kind), &
         dimension(jmax_segments,2), intent(in):: cross_lat_eul_index
    integer (kind=int_kind), intent(in)            ::jx_min, jx_max, jy_min, jy_max
    real (kind=real_kind), dimension(-nhe:nc+2+nhe), intent(in) :: xgno
    real (kind=real_kind)   ,  &
         dimension(jmax_segments,nreconstruction), intent(inout) :: weights
    integer (kind=int_kind),  &
         dimension(jmax_segments,2), intent(inout) :: weights_eul_index
    real (kind=real_kind)   , dimension(nreconstruction) :: weights_tmp
    
    integer (kind=int_kind) :: i,j,k, isgn, h
    logical                 :: ltmp,lcontinue,lclockwise
    
    real (kind=real_kind), dimension(jmax_segments,2)  :: r_cross_lat_seg
    integer (kind=int_kind), dimension(jmax_segments,2):: cross_lat_eul_index_seg
    
    real (kind=real_kind), dimension(jmax_segments,2)  :: r_cross_lat_seg2
    integer (kind=int_kind), dimension(jmax_segments,2):: cross_lat_eul_index_seg2
    
    integer (kind=int_kind) :: count,js,is
    real (kind=real_kind) :: a,a2,b,b2
    
    if (ldbg_global) then
      WRITE(*,*) "from non_convex"
    end if
    
    if (ldbg_global) then
      OPEN(unit=40, file='inner_integral.dat',status='replace')
      WRITE(40,*) "  "
      CLOSE(40)              
      OPEN(unit=41, file='inner_nonconvex.dat',status='replace')
      WRITE(41,*) "  "
      CLOSE(41)              
    end if
    
    
    
    
    if (jcross_lat>0) then      
      do i=MINVAL(cross_lat_eul_index(1:jcross_lat,2)),MAXVAL(cross_lat_eul_index(1:jcross_lat,2))
        !        WRITE(*,*) "looking at latitude ",i !xxxx
        count = 1
        !
        ! find all crossings with Eulerian latitude i
        !
        do k=1,jcross_lat
          if (cross_lat_eul_index(k,2)==i) then
            !            WRITE(*,*) "other crossings with latitude",i ," is ",k!xxxx
            r_cross_lat_seg        (count,:) = r_cross_lat        (k,:)
            cross_lat_eul_index_seg(count,:) = cross_lat_eul_index(k,:)
            
            IF (ldbg_global) then
              WRITE(*,*)  r_cross_lat_seg(count,1),r_cross_lat_seg(count,2)
              WRITE(*,*) "  "
            END IF
            count = count+1
          end if
        enddo
        count = count-1
        IF (ABS((count/2)-DBLE(count)/2.0)<tiny) then
          IF (count.NE.2) THEN
            WRITE(*,*) "non-convex cell", count
            !
            ! sort array from min to max
            !
            !        WRITE(*,*) "before ordering",r_cross_lat_seg(1:count,1)
            do js=2, count
              a =r_cross_lat_seg(js,1)
              a2=r_cross_lat_seg(js,2)
              b =cross_lat_eul_index_seg(js,1) 
              b2=cross_lat_eul_index_seg(js,2) 
              do is=js-1,1,-1
                if (r_cross_lat_seg(is,1)<=a) goto 10
                r_cross_lat_seg(is+1,:)=r_cross_lat_seg(is,:)
                cross_lat_eul_index_seg(is+1,:) = cross_lat_eul_index_seg(is,:)
              end do
              is=0
10            r_cross_lat_seg(is+1,1)=a
              r_cross_lat_seg(is+1,2)=a2
              cross_lat_eul_index_seg(is+1,1) = b
              cross_lat_eul_index_seg(is+1,2) = b2
            end do
            r_cross_lat_seg2        (1:count,:) = r_cross_lat_seg        (1:count,:)
            cross_lat_eul_index_seg2(1:count,:) = cross_lat_eul_index_seg(1:count,:)
          end if
        else
          WRITE(*,*) "INCONSISTENCY in number of crossings!", count
          STOP
        END IF
        !
        ! only do every other segment
        !
        IF (ldbg_global) THEN
          WRITE(*,*) "segments send to compute_inner_line_integrals_lat"
        END IF
        do h=1,count-1,2
          r_cross_lat_seg2        (1:2,:) = r_cross_lat_seg        (h:h+1,:)
          cross_lat_eul_index_seg2(1:2,:) = cross_lat_eul_index_seg(h:h+1,:)
          
          IF (ldbg_global) THEN
            OPEN(unit=41, file='inner_nonconvex.dat',status='old',access='append')
            WRITE(41,*) r_cross_lat_seg2(1,1),r_cross_lat_seg2(1,2)
            WRITE(41,*) r_cross_lat_seg2(2,1),r_cross_lat_seg2(2,2)
            WRITE(41,*) "  "
            CLOSE(41)                          
            
            WRITE(*,*) "h=",h
            WRITE(*,*) "from ",r_cross_lat_seg(h,1),r_cross_lat_seg(h,2)
            WRITE(*,*) "to ",r_cross_lat_seg(h+1,1),r_cross_lat_seg(h+1,2)
            WRITE(*,*) "jumping over"
            WRITE(*,*) "from ",r_cross_lat_seg(h+1,1),r_cross_lat_seg(h+1,2)
            WRITE(*,*) "to ",r_cross_lat_seg(h+2,1),r_cross_lat_seg(h+2,2)
          END IF
          
          call compute_inner_line_integrals_lat(r_cross_lat_seg2,cross_lat_eul_index_seg2,&
               2,jsegment,jmax_segments,xgno,jx_min,jx_max,jy_min, jy_max,weights,weights_eul_index,&
               nreconstruction,ngauss,gauss_weights,abscissae)!phl add jx_min etc.
        end do
        
      enddo
    endif
  end subroutine compute_inner_line_integrals_lat_nonconvex

  
  
  !
  ! line integral from (a1_in,a2_in) to (b1_in,b2_in)
  ! If line is coniciding with an Eulerian longitude or latitude the routine
  ! needs to know where an adjacent side is located to determine which
  ! reconstruction must be used. therefore (c1,c2) is passed to the routine
  !
  !   
  
  subroutine side_integral(&
       x_in,y_in,nvertex,jsegment,jmax_segments,&
       weights,weights_eul_index,nreconstruction,jx,jy,xgno,ygno,jx_min,jx_max,jy_min,jy_max,&
       ngauss,gauss_weights,abscissae,&!)!phl add jx_min etc.
       jcross_lat,r_cross_lat,cross_lat_eul_index)
    implicit none
    !
    ! for Gaussian quadrature
    !
    real (kind=real_kind), dimension(ngauss), intent(in) :: gauss_weights, abscissae
    real (kind=real_kind), dimension(1:nvertex)        , intent(in)    :: x_in,y_in
    
    integer (kind=int_kind), intent(in)               :: jx_min, jy_min, jx_max, jy_max
    integer (kind=int_kind), intent(in)               :: nvertex
    real (kind=real_kind), dimension(-nhe:nc+2+nhe), intent(in) :: xgno
    real (kind=real_kind), dimension(-nhe:nc+2+nhe), intent(in) :: ygno
    integer (kind=int_kind),            intent(inout) :: jsegment
    integer (kind=int_kind),            intent(in)    :: nreconstruction,jx,jy,jmax_segments,ngauss
    real (kind=real_kind)   ,  &
         dimension(jmax_segments,nreconstruction), intent(out) :: weights
    integer (kind=int_kind),  &
         dimension(jmax_segments,2), intent(out) :: weights_eul_index
    !
    ! variables for registering crossings with Eulerian latitudes and longitudes
    !
    integer (kind=int_kind),         intent(inout):: jcross_lat
    !
    ! max. crossings per side is 2*nhe
    !
    real (kind=real_kind), &
         dimension(jmax_segments,2), intent(inout):: r_cross_lat
    integer (kind=int_kind), &
         dimension(jmax_segments,2), intent(inout):: cross_lat_eul_index
    !
    ! local variables
    !
    real (kind=real_kind) :: dist_lon,dist_lat, tmp_a1, tmp_a2, tmp_x(1), tmp_b2, a1, a2, b2
    real (kind=real_kind) :: dist
    real (kind=real_kind), dimension(2) :: xseg,yseg 
    real (kind=real_kind), dimension(0:3) :: x,y
    real (kind=real_kind)               :: lon_x,lat_y,lon_y,lat_x
    real (kind=real_kind)               :: xeul,yeul,xcross,ycross,slope
    integer (kind=int_kind) ::    jx_eul_tmp,jy_eul_tmp
    integer (kind=int_kind)            :: xsgn1,ysgn1,xsgn2,ysgn2
    integer (kind=int_kind) :: ifrom_left, iter,previous_jy_eul_cross
    logical :: lcontinue, lregister_cross, lsame_cell_x, lsame_cell_y
    
    integer (kind=int_kind) :: jx_eul, jy_eul, side_count,jdbg
    real (kind=real_kind), dimension(0:nvertex+2)  :: xcell,ycell
    real (kind=real_kind), dimension(0:nvertex+2)  :: xcell_tmp,ycell_tmp
    
    if (ldbg_global) then
      OPEN(unit=40, file='side_integral.dat',status='replace')
      WRITE(40,*) "  "
      CLOSE(40)              
    end if
    
5   FORMAT(10e14.6)
    !
    !***********************************************
    !
    ! find jx_eul and jy_eul for (x(1),y(1))
    !
    !***********************************************
    !
    jx_eul = jx; jy_eul = jy    
    xcell(1:nvertex)=x_in; ycell(1:nvertex)=y_in
    DO iter=1,nvertex
      CALL truncate_vertex(xcell(iter),jx_eul,xgno)
      CALL truncate_vertex(ycell(iter),jy_eul,ygno)
    END DO
    xcell(0) = xcell(nvertex); xcell(nvertex+1)=xcell(1); xcell(nvertex+2)=xcell(2);
    ycell(0) = ycell(nvertex); ycell(nvertex+1)=ycell(1); ycell(nvertex+2)=ycell(2);
    
!    IF (ldbgr) THEN
!      WRITE(*,*) "from side_integral: cell vertices"
!      DO iter=1,nvertex
!        WRITE(*,*) "x(iter),y(iter)",iter, xcell(iter),ycell(iter)
!      END DO
!    END IF
    
    IF (MAXVAL(xcell).LE.xgno(jx_min).OR.MINVAL(xcell).GE.xgno(jx_max).OR.&
         MAXVAL(ycell).LE.ygno(jy_min).OR.MINVAL(ycell).GE.ygno(jy_max)) THEN
      
!      IF (ldbgr)  WRITE(*,*) "entire cell off panel"
    ELSE             
      jx_eul = jx
      jy_eul = jy
      CALL which_eul_cell(xcell(1:3),jx_eul,xgno)
      CALL which_eul_cell(ycell(1:3),jy_eul,ygno)
!      IF (ldbgr) WRITE(*,*) "x(1),y(1) in cell",jx_eul,jy_eul
     
      side_count = 1
      DO WHILE (side_count<nvertex+1)
        jdbg = 0
        iter = 0
        lcontinue = .TRUE.
        x(0:3) = xcell(side_count-1:side_count+2); y(0:3) = ycell(side_count-1:side_count+2); 
!        IF (ldbgr) WRITE(*,*) "+++++++++++++++++++++++++++++++++++++++"
!        IF (ldbgr) WRITE(*,*) "side",side_count
        DO while (lcontinue)
!          IF (ldbgr) WRITE(*,*) "iter",iter
!          IF (ldbgr) WRITE(*,*) "x,y(1)",x(1),y(1)
!          IF (ldbgr) WRITE(*,*) "x,y(2)",x(2),y(2)
!          IF (ldbgr) WRITE(*,*) "jx_eul,jy_eul",jx_eul,jy_eul
!          IF (ldbgr) WRITE(*,*) "xgno",xgno(jx_eul),xgno(jx_eul+1)
!          IF (ldbgr) WRITE(*,*) "ygno",ygno(jy_eul),ygno(jy_eul+1)          
          iter = iter+1
          IF (iter>1000) THEN
            WRITE(*,*) "search not converging",iter
            STOP
          END IF
          lsame_cell_x = (x(2).GE.xgno(jx_eul).AND.x(2).LE.xgno(jx_eul+1))
          lsame_cell_y = (y(2).GE.ygno(jy_eul).AND.y(2).LE.ygno(jy_eul+1))
!          IF (ldbgr) WRITE(*,*) "lsame_cell_x,lsame_cell_y=",lsame_cell_x,lsame_cell_y
          IF (lsame_cell_x.AND.lsame_cell_y) THEN
            !
            !****************************
            !
            ! same cell integral
            !
            !****************************
            !
!            IF (ldbgr) WRITE(*,*) "same cell integral",jx_eul,jy_eul
            xseg(1) = x(1); yseg(1) = y(1); xseg(2) = x(2); yseg(2) = y(2)
            jx_eul_tmp = jx_eul; jy_eul_tmp = jy_eul; 
            lcontinue = .FALSE.
            !
            ! prepare for next side if (x(2),y(2)) is on a grid line
            !
            IF (x(2).EQ.xgno(jx_eul+1).AND.x(3)>xgno(jx_eul+1)) THEN
              !
              ! cross longitude jx_eul+1
              !
!              IF (ldbgr) WRITE(*,*) "cross longitude",jx_eul+1
              jx_eul=jx_eul+1
            ELSE IF (x(2).EQ.xgno(jx_eul  ).AND.x(3)<xgno(jx_eul)) THEN
              !
              ! cross longitude jx_eul
              !
!              IF (ldbgr) WRITE(*,*) "cross longitude",jx_eul
              jx_eul=jx_eul-1
            END IF
            IF (y(2).EQ.ygno(jy_eul+1).AND.y(3)>ygno(jy_eul+1)) THEN
              !
              ! register crossing with latitude: line-segments point Northward
              !
              jcross_lat = jcross_lat + 1
              jy_eul     = jy_eul     + 1
!              IF (ldbgr) WRITE(*,*) "cross latitude",jy_eul
              cross_lat_eul_index(jcross_lat,1) = jx_eul
              cross_lat_eul_index(jcross_lat,2) = jy_eul
              r_cross_lat(jcross_lat,1) = x(2)
              r_cross_lat(jcross_lat,2) = y(2)
            ELSE IF (y(2).EQ.ygno(jy_eul  ).AND.y(3)<ygno(jy_eul)) THEN
              !
              ! register crossing with latitude: line-segments point Southward
              !
!              IF (ldbgr) WRITE(*,*) "cross latitude",jy_eul
              jcross_lat = jcross_lat+1
              cross_lat_eul_index(jcross_lat,1) = jx_eul
              cross_lat_eul_index(jcross_lat,2) = jy_eul
              r_cross_lat(jcross_lat,1) = x(2)
              r_cross_lat(jcross_lat,2) = y(2)
              
              jy_eul=jy_eul-1
            END IF
            lcontinue=.FALSE.
          ELSE
            !
            !****************************
            !
            ! not same cell integral
            !
            !****************************
            !
            IF (lsame_cell_x) THEN
!              IF (ldbgr) WRITE(*,*) "same cell x"
              ysgn1 = (1+INT(SIGN(1.0D0,y(2)-y(1))))/2 !"1" if y(2)>y(1) else "0"
              ysgn2 = INT(SIGN(1.0D0,y(2)-y(1)))       !"1" if y(2)>y(1) else "-1"
              !
              !*******************************************************************************
              !
              ! there is at least one crossing with latitudes but no crossing with longitudes
              !
              !*******************************************************************************
              !
              yeul   = ygno(jy_eul+ysgn1)
              IF (x(1).EQ.x(2)) THEN
                !
                ! line segment is parallel to longitude (infinite slope)
                !
!                IF (ldbgr) WRITE(*,*) "line segment parallel to longitude"
                xcross = x(1)
              ELSE
                slope  = (y(2)-y(1))/(x(2)-x(1))
                xcross = x_cross_eul_lat(x(1),y(1),yeul,slope)
                !
                ! constrain crossing to be "physically" possible
                !
                xcross = MIN(MAX(xcross,xgno(jx_eul)),xgno(jx_eul+1))

                
!                IF (ldbgr) WRITE(*,*) "cross latitude"
                !
                ! debugging
                !
                IF (xcross.GT.xgno(jx_eul+1).OR.xcross.LT.xgno(jx_eul)) THEN
                  WRITE(*,*) "xcross is out of range",jx,jy
                  WRITE(*,*) "xcross-xgno(jx_eul+1), xcross-xgno(jx_eul))",&
                       xcross-xgno(jx_eul+1), xcross-ygno(jx_eul)
                  STOP
                END IF
              END IF
              xseg(1) = x(1); yseg(1) = y(1); xseg(2) = xcross; yseg(2) = yeul
              jx_eul_tmp = jx_eul; jy_eul_tmp = jy_eul; 
              !
              ! prepare for next iteration
              !
              x(0) = x(1); y(0) = y(1); x(1) = xcross; y(1) = yeul; jy_eul = jy_eul+ysgn2
              !
              ! register crossing with latitude
              !
              jcross_lat = jcross_lat+1
              cross_lat_eul_index(jcross_lat,1) = jx_eul
              if (ysgn2>0) then                
                cross_lat_eul_index(jcross_lat,2) = jy_eul
              else
                cross_lat_eul_index(jcross_lat,2) = jy_eul+1
              end if
              r_cross_lat(jcross_lat,1) = xcross
              r_cross_lat(jcross_lat,2) = yeul
            ELSE IF (lsame_cell_y) THEN
!              IF (ldbgr) WRITE(*,*) "same cell y"
              !
              !*******************************************************************************
              !
              ! there is at least one crossing with longitudes but no crossing with latitudes
              !
              !*******************************************************************************
              !
              xsgn1 = (1+INT(SIGN(1.0D0,x(2)-x(1))))/2 !"1" if x(2)>x(1) else "0"
              xsgn2 = INT(SIGN(1.0D0,x(2)-x(1))) !"1" if x(2)>x(1) else "-1"
              xeul   = xgno(jx_eul+xsgn1)
!              IF (ldbgr) WRITE(*,*) " crossing longitude",jx_eul+xsgn1
              IF (ABS(x(2)-x(1))<fuzzy_width) THEN
                ycross = 0.5*(y(2)-y(1))
                !                IF (ldbgr) WRITE(*,*) "fuzzy crossing"
              ELSE
                slope  = (y(2)-y(1))/(x(2)-x(1))
                ycross = y_cross_eul_lon(x(1),y(1),xeul,slope)
              END IF
              !
              ! constrain crossing to be "physically" possible
              !
              ycross = MIN(MAX(ycross,ygno(jy_eul)),ygno(jy_eul+1))
              
              !
              ! debugging
              !
              IF (ycross.GT.ygno(jy_eul+1).OR.ycross.LT.ygno(jy_eul)) THEN
                WRITE(*,*) "ycross is out of range"
                WRITE(*,*) "jx,jy,jx_eul,jy_eul",jx,jy,jx_eul,jy_eul
                WRITE(*,*) "ycross-ygno(jy_eul+1), ycross-ygno(jy_eul))",&
                     ycross-ygno(jy_eul+1), ycross-ygno(jy_eul)
                STOP
              END IF
              xseg(1) = x(1); yseg(1) = y(1); xseg(2) = xeul; yseg(2) = ycross
              jx_eul_tmp = jx_eul; jy_eul_tmp = jy_eul; 
              !
              ! prepare for next iteration
              !
              x(0) = x(1); y(0) = y(1); x(1) = xeul; y(1) = ycross; jx_eul = jx_eul+xsgn2
            ELSE
!              IF (ldbgr) WRITE(*,*) "not same cell x; not same cell y"
              !
              !*******************************************************************************
              !
              ! there are crossings with longitude(s) and latitude(s)
              !
              !*******************************************************************************
              ! 
              xsgn1 = (1+INT(SIGN(1.0D0,x(2)-x(1))))/2 !"1" if x(2)>x(1) else "0"
              xsgn2 = (INT(SIGN(1.0D0,x(2)-x(1)))) !"1" if x(2)>x(1) else "0"
              xeul   = xgno(jx_eul+xsgn1) 
              ysgn1 = (1+INT(SIGN(1.0D0,y(2)-y(1))))/2 !"1" if y(2)>y(1) else "0"
              ysgn2 = INT(SIGN(1.0D0,y(2)-y(1)))       !"1" if y(2)>y(1) else "-1"
              yeul   = ygno(jy_eul+ysgn1)
              
              slope  = (y(2)-y(1))/(x(2)-x(1))
              IF (ABS(x(2)-x(1))<fuzzy_width) THEN
                ycross = 0.5*(y(2)-y(1))
              ELSE
                ycross = y_cross_eul_lon(x(1),y(1),xeul,slope)
              END IF
              xcross = x_cross_eul_lat(x(1),y(1),yeul,slope)
              
              IF ((xsgn2>0.AND.xcross.LE.xeul).OR.(xsgn2<0.AND.xcross.GE.xeul)) THEN
                !
                ! cross latitude
                !
!                IF (ldbgr) WRITE(*,*) "crossing latitude",jy_eul+ysgn1
                xseg(1) = x(1); yseg(1) = y(1); xseg(2) = xcross; yseg(2) = yeul
                jx_eul_tmp = jx_eul; jy_eul_tmp = jy_eul; 
                !
                ! prepare for next iteration
                !
                x(0) = x(1); y(0) = y(1); x(1) = xcross; y(1) = yeul; jy_eul = jy_eul+ysgn2
                !
                ! register crossing with latitude
                !
                jcross_lat = jcross_lat+1
                cross_lat_eul_index(jcross_lat,1) = jx_eul
                if (ysgn2>0) then                
                  cross_lat_eul_index(jcross_lat,2) = jy_eul
                else
                  cross_lat_eul_index(jcross_lat,2) = jy_eul+1
                end if
                r_cross_lat(jcross_lat,1) = xcross
                r_cross_lat(jcross_lat,2) = yeul
              ELSE
                !
                ! cross longitude
                !
!                IF (ldbgr) WRITE(*,*) "crossing longitude",jx_eul+xsgn1
                xseg(1) = x(1); yseg(1) = y(1); xseg(2) = xeul; yseg(2) = ycross
                jx_eul_tmp = jx_eul; jy_eul_tmp = jy_eul; 
                !
                ! prepare for next iteration
                !
                x(0) = x(1); y(0) = y(1); x(1) = xeul; y(1) = ycross; jx_eul = jx_eul+xsgn2
              END IF
              
            END IF
          END IF
          !
          ! register line-segment (don't register line-segment if outside of panel)
          !
          if (jx_eul_tmp>=jx_min.AND.jy_eul_tmp>=jy_min.AND.&
               jx_eul_tmp<=jx_max-1.AND.jy_eul_tmp<=jy_max-1) then
            !               jx_eul_tmp<=jx_max-1.AND.jy_eul_tmp<=jy_max-1.AND.side_count<3) then
            jsegment=jsegment+1
            weights_eul_index(jsegment,1) = jx_eul_tmp
            weights_eul_index(jsegment,2) = jy_eul_tmp
            call get_weights_gauss(weights(jsegment,1:nreconstruction),&
                 xseg,yseg,nreconstruction,ngauss,gauss_weights,abscissae)
            
!            if (ldbg_global) then
!              OPEN(unit=40, file='side_integral.dat',status='old',access='append')
!              WRITE(40,*) xseg(1),yseg(1)
!              WRITE(40,*) xseg(2),yseg(2)
!              WRITE(40,*) "  "
!              CLOSE(40)              
!            end if
            
            
            jdbg=jdbg+1

            if (xseg(1).EQ.xseg(2))then
              slope = bignum
            else if (abs(yseg(1) -yseg(2))<fuzzy_width) then
              slope = 0.0
            else
              slope    = (yseg(2)-yseg(1))/(xseg(2)-xseg(1))
            end if
          ELSE
!            IF (ldbgr) WRITE(*,*) "segment outside of panel"
          END IF         
        END DO
        side_count = side_count+1
      END DO
    END IF
  end subroutine side_integral
 

  real (kind=real_kind) function compute_slope(x,y)
    implicit none
    real (kind=real_kind), dimension(2), intent(in) :: x,y
    if (fuzzy(ABS(x(2)-x(1)),fuzzy_width)>0) THEN
      compute_slope = (y(2)-y(1))/(x(2)-x(1))
    else
      compute_slope = bignum
    end if
  end function compute_slope

  real (kind=real_kind) function y_cross_eul_lon(x,y,xeul,slope)
    implicit none
    real (kind=real_kind), intent(in) :: x,y
    real (kind=real_kind)              , intent(in) :: xeul,slope
    ! line: y=a*x+b
    real (kind=real_kind) :: a,b
    b = y-slope*x 
    y_cross_eul_lon = slope*xeul+b
  end function y_cross_eul_lon

  real (kind=real_kind) function x_cross_eul_lat(x,y,yeul,slope)
    implicit none
    real (kind=real_kind), intent(in) :: x,y
    real (kind=real_kind)              , intent(in) :: yeul,slope

    if (fuzzy(ABS(slope),fuzzy_width)>0) THEN
        x_cross_eul_lat = x+(yeul-y)/slope
    ELSE
      !      WRITE(*,*) "WARNING: slope is epsilon - ABORT"
      x_cross_eul_lat = bignum
    END IF
  end function x_cross_eul_lat

  subroutine get_weights_exact(weights,xseg,yseg,nreconstruction)
!    use cslam_analytic_mod, only: I_00, I_10, I_01, I_20, I_02, I_11
    implicit none
    integer (kind=int_kind), intent(in) :: nreconstruction
    real (kind=real_kind), dimension(nreconstruction), intent(out) :: weights
    real (kind=real_kind), dimension(2     ), intent(in) :: xseg,yseg
    !
    ! compute weights
    !
    real (kind=real_kind) :: tmp,slope,b,integral,dx2,xc
    integer (kind=int_kind) :: i
!    weights(:) = -half*(xseg(1)*yseg(2)-xseg(2)*yseg(1)) !dummy for testing

    weights(1) = ((I_00(xseg(2),yseg(2))-I_00(xseg(1),yseg(1))))
    if (ABS(weights(1))>1.0) THEN
      WRITE(*,*) "1 exact weights(jsegment)",weights(1),xseg,yseg
      stop
    end if
    if (nreconstruction>1) then
       weights(2) = ((I_10(xseg(2),yseg(2))-I_10(xseg(1),yseg(1))))
       weights(3) = ((I_01(xseg(2),yseg(2))-I_01(xseg(1),yseg(1))))
    endif
    if (nreconstruction>3) then
       weights(4) = ((I_20(xseg(2),yseg(2))-I_20(xseg(1),yseg(1))))
       weights(5) = ((I_02(xseg(2),yseg(2))-I_02(xseg(1),yseg(1))))
       weights(6) = ((I_11(xseg(2),yseg(2))-I_11(xseg(1),yseg(1))))
    endif

  end subroutine get_weights_exact



  subroutine get_weights_gauss(weights,xseg,yseg,nreconstruction,ngauss,gauss_weights,abscissae)
    implicit none
    integer (kind=int_kind), intent(in) :: nreconstruction,ngauss
    real (kind=real_kind), dimension(nreconstruction), intent(out) :: weights
    real (kind=real_kind), dimension(2     ), intent(in) :: xseg,yseg
    real (kind=real_kind) :: slope
    !
    ! compute weights
    !
    !
    ! for Gaussian quadrature
    !
    real (kind=real_kind), dimension(ngauss), intent(in) :: gauss_weights, abscissae

    ! if line-segment parallel to x or y use exact formulaes else use qudrature
    !
    real (kind=real_kind) :: tmp,b,integral,dx2,xc,x,y
    integer (kind=int_kind) :: i




!    if (fuzzy(abs(xseg(1) -xseg(2)),fuzzy_width)==0)then
    if (xseg(1).EQ.xseg(2))then
      weights = 0.0D0
    else if (abs(yseg(1) -yseg(2))<fuzzy_width) then
      !
      ! line segment parallel to latitude - compute weights exactly
      !
      if (ldbgr) write(*,*) "line segment parallel to latitude - compute weights exactly"
      weights(1) = ((I_00(xseg(2),yseg(2))-I_00(xseg(1),yseg(1))))
      if (nreconstruction>1) then
        weights(2) = ((I_10(xseg(2),yseg(2))-I_10(xseg(1),yseg(1))))
        weights(3) = ((I_01(xseg(2),yseg(2))-I_01(xseg(1),yseg(1))))
      endif
      if (nreconstruction>3) then
        weights(4) = ((I_20(xseg(2),yseg(2))-I_20(xseg(1),yseg(1))))
        weights(5) = ((I_02(xseg(2),yseg(2))-I_02(xseg(1),yseg(1))))
        weights(6) = ((I_11(xseg(2),yseg(2))-I_11(xseg(1),yseg(1))))
      endif
    else
      
      
      slope    = (yseg(2)-yseg(1))/(xseg(2)-xseg(1))
      b        = yseg(1)-slope*xseg(1)
      dx2      = 0.5D0*(xseg(2)-xseg(1))
      if (ldbgr) WRITE(*,*) "dx2 and slope in gauss weight",dx2,slope
      xc       = 0.5D0*(xseg(1)+xseg(2))
      integral = 0.0D0
      do i=1,ngauss
        x        = xc+abscissae(i)*dx2
        y        = slope*x+b
        integral = integral+gauss_weights(i)*F_00(x,y)
      enddo
      weights(1) = integral*dx2  
      if (nreconstruction>1) then
        integral = 0.0D0
        do i=1,ngauss
          x        = xc+abscissae(i)*dx2
          y        = slope*x+b
          integral = integral+gauss_weights(i)*F_10(x,y)
        enddo
        weights(2) = integral*dx2  
        integral = 0.0D0
        do i=1,ngauss
          x        = xc+abscissae(i)*dx2
          y        = slope*x+b
          integral = integral+gauss_weights(i)*F_01(x,y)
        enddo
        weights(3) = integral*dx2  
      endif
      if (nreconstruction>3) then
        integral = 0.0D0
        do i=1,ngauss
          x        = xc+abscissae(i)*dx2
          y        = slope*x+b
          integral = integral+gauss_weights(i)*F_20(x,y)
        enddo
        weights(4) = integral*dx2  
        integral = 0.0D0
        do i=1,ngauss
          x        = xc+abscissae(i)*dx2
          y        = slope*x+b
          integral = integral+gauss_weights(i)*F_02(x,y)
        enddo
        weights(5) = integral*dx2  
        integral = 0.0D0
        do i=1,ngauss
          x        = xc+abscissae(i)*dx2
          y        = slope*x+b
          integral = integral+gauss_weights(i)*F_11(x,y)
        enddo
        weights(6) = integral*dx2  
      endif
    end if
  end subroutine get_weights_gauss

  real (kind=real_kind) function F_00(x_in,y_in)
    implicit none
    real (kind=real_kind), intent(in) :: x_in,y_in
    real (kind=real_kind) :: x,y,tmp

    x = x_in
    y = y_in

    F_00 =y/((1.0D0+x*x)*SQRT(1.0D0+x*x+y*y))
  end function F_00

  real (kind=real_kind) function F_10(x_in,y_in)
    implicit none
    real (kind=real_kind), intent(in) :: x_in,y_in
    real (kind=real_kind) :: x,y,tmp

    x = x_in
    y = y_in

    F_10 =x*y/((1.0D0+x*x)*SQRT(1.0D0+x*x+y*y))
  end function F_10

  real (kind=real_kind) function F_01(x_in,y_in)
    implicit none
    real (kind=real_kind), intent(in) :: x_in,y_in
    real (kind=real_kind) :: x,y,tmp

    x = x_in
    y = y_in

    F_01 =-1.0D0/(SQRT(1.0D0+x*x+y*y))
  end function F_01

  real (kind=real_kind) function F_20(x_in,y_in)
    implicit none
    real (kind=real_kind), intent(in) :: x_in,y_in
    real (kind=real_kind) :: x,y,tmp

    x = x_in
    y = y_in

    F_20 =x*x*y/((1.0D0+x*x)*SQRT(1.0D0+x*x+y*y))
  end function F_20

  real (kind=real_kind) function F_02(x_in,y_in)
    implicit none
    real (kind=real_kind), intent(in) :: x_in,y_in
    real (kind=real_kind) :: x,y,alpha, tmp

    x = x_in
    y = y_in

    alpha = ATAN(x)
    tmp=y*COS(alpha)
    F_02 =-y/SQRT(1.0D0+x*x+y*y)+log(tmp+sqrt(tmp*tmp+1))
    
    !
    ! cos(alpha) = 1/sqrt(1+x*x)
    !
  end function F_02

  real (kind=real_kind) function F_11(x_in,y_in)
    implicit none
    real (kind=real_kind), intent(in) :: x_in,y_in
    real (kind=real_kind) :: x,y,tmp

    x = x_in
    y = y_in

    F_11 =-x/(SQRT(1.0D0+x*x+y*y))
  end function F_11

  subroutine which_eul_cell(x,j_eul,gno)
    implicit none
    integer (kind=int_kind)                               , intent(inout) :: j_eul
    real (kind=real_kind), dimension(3)                    , intent(in)    :: x
    real (kind=real_kind), dimension(-nhe:nc+2+nhe), intent(in)    :: gno !phl
!    real (kind=real_kind), intent(in)    :: eps
    
    real (kind=real_kind) :: d1,d2,d3,d1p1
    logical                 :: lcontinue
    integer :: iter


    !
    !  this is not needed in transport code search
    !
!    IF (x(1)<gno(-nhe)) j_eul=-nhe
!    IF (x(1)>gno(nc+2+nhe)) j_eul=nc+1+nhe
!    RETURN

!    j_eul = MIN(MAX(j_eul,-nhe),nc+1+nhe) !added
    
    lcontinue = .TRUE.
    iter = 0 
    IF (ldbgr) WRITE(*,*) "from which_eul_cell",x(1),x(2),x(3)
    DO WHILE (lcontinue)
      iter = iter+1 
      IF (x(1).GE.gno(j_eul).AND.x(1).LT.gno(j_eul+1)) THEN
        lcontinue = .FALSE.
        !
        ! special case when x(1) is on top of grid line
        !
        IF (x(1).EQ.gno(j_eul)) THEN
!        IF (ABS(x(1)-gno(j_eul))<tiny) THEN
          IF (ldbgr) WRITE(*,*) "x(1) is on top of gno(J_eul)"
          IF (x(2).GT.gno(j_eul)) THEN
            j_eul = j_eul
          ELSE IF (x(2).LT.gno(j_eul)) THEN
            j_eul = j_eul-1
          ELSE
            IF (ldbgr) WRITE(*,*) "x(2) is on top of gno(J_eul)"
            !
            ! x(2) is on gno(j_eul) grid line; need x(3) to determine Eulerian cell 
            !
            IF (x(3).GT.gno(j_eul)) THEN
              IF (ldbgr) WRITE(*,*) "x(3) to the right"
              j_eul = j_eul
            ELSE IF (x(3).LT.gno(j_eul)) THEN
              IF (ldbgr) WRITE(*,*) "x(3) to the left x(3)-x(2),x(3),x(2)",x(3)-x(2),x(3),x(2)
              j_eul = j_eul-1
            ELSE
              WRITE(*,*) "inconsistent cell: x(1)=x(2)=x(3)", x(1),x(2),x(3)
!              WRITE(*,*) "gno(j_eul),j_eul",gno(j_eul),j_eul
              STOP
            END IF
          END IF
        END IF
      ELSE
        ! 
        ! searching - prepare for next iteration
        !
        IF (x(1).GE.gno(j_eul+1)) THEN
          j_eul = j_eul + 1
        ELSE
          !
          ! x(1).LT.gno(j_eul)
          !
          j_eul = j_eul - 1
        END IF
      END IF
      IF (iter>1000.OR.j_eul<-nhe.OR.j_eul>nc+2+nhe) THEN
        WRITE(*,*) "search in which_eul_cell not converging!", iter,j_eul
        WRITE(*,*) "input", x
        WRITE(*,*) "gno", gno(nc),gno(nc+1),gno(nc+2),gno(nc+3)
        STOP
      END IF
    END DO
  END subroutine which_eul_cell


  subroutine truncate_vertex(x,j_eul,gno)
    implicit none
    integer (kind=int_kind)                               , intent(inout) :: j_eul
    real (kind=real_kind)                    , intent(inout)    :: x
    real (kind=real_kind), dimension(-nhe:nc+2+nhe), intent(in)    :: gno !phl
!    real (kind=real_kind), intent(in)    :: eps
    
    logical                 :: lcontinue
    integer :: iter
    real (kind=real_kind) :: xsgn,dist,dist_new,tmp
    
    !
    !  this is not needed in transport code search
    !
!    IF (x<gno(-nhe)) j_eul=-nhe
!    IF (x>gno(nc+2+nhe)) j_eul=nc+1+nhe
!
!    RETURN


    lcontinue = .TRUE.
    iter = 0 
    dist = bignum
!    j_eul = MIN(MAX(j_eul,-nhe),nc+1+nhe) !added
    xsgn     = INT(SIGN(1.0_dbl_kind,x-gno(j_eul)))
    DO WHILE (lcontinue)
      iter     = iter+1 
      tmp      = x-gno(j_eul)
      dist_new = ABS(tmp)
      IF (dist_new>dist) THEN
        lcontinue = .FALSE.
!      ELSE IF (ABS(tmp)<1.0E-11) THEN
      ELSE IF (ABS(tmp)<1.0E-9) THEN
!      ELSE IF (ABS(tmp)<1.0E-4) THEN
        x = gno(j_eul)
        lcontinue = .FALSE.
      ELSE
        j_eul = j_eul+xsgn
        dist = dist_new
      END IF
      IF (iter>10000) THEN
        WRITE(*,*) "truncate vertex not converging"
        STOP
      END IF
    END DO
  END subroutine truncate_vertex




!********************************************************************************
!
! Gauss-Legendre quadrature
!
! Tabulated values
!
!********************************************************************************
subroutine gauss_points(n,weights,points)
  implicit none
  real (kind=real_kind), dimension(n), intent(out) :: weights, points
  integer (kind=int_kind)           , intent(in ) :: n
  
  select case (n)
!    CASE(1)
!       abscissae(1) = 0.0D0
!       weights(1)   = 2.0D0
  case(2)
     points(1)    = -sqrt(1.0D0/3.0D0)
     points(2)    =  sqrt(1.0D0/3.0D0)
     weights(1)   =  1.0D0
     weights(2)   =  1.0D0
  case(3)
     points(1)    = -0.774596669241483377035853079956D0
     points(2)    =  0.0D0
     points(3)    =  0.774596669241483377035853079956D0
     weights(1)   =  0.555555555555555555555555555556D0
     weights(2)   =  0.888888888888888888888888888889D0
     weights(3)   =  0.555555555555555555555555555556D0
  case(4)
     points(1)    = -0.861136311594052575223946488893D0
     points(2)    = -0.339981043584856264802665659103D0
     points(3)    =  0.339981043584856264802665659103D0
     points(4)    =  0.861136311594052575223946488893D0
     weights(1)   =  0.347854845137453857373063949222D0
     weights(2)   =  0.652145154862546142626936050778D0 
     weights(3)   =  0.652145154862546142626936050778D0 
     weights(4)   =  0.347854845137453857373063949222D0      
  case(5)
     points(1)    = -(1.0D0/3.0D0)*sqrt(5.0D0+2.0D0*sqrt(10.0D0/7.0D0))
     points(2)    = -(1.0D0/3.0D0)*sqrt(5.0D0-2.0D0*sqrt(10.0D0/7.0D0))
     points(3)    =  0.0D0
     points(4)    =  (1.0D0/3.0D0)*sqrt(5.0D0-2.0D0*sqrt(10.0D0/7.0D0))
     points(5)    =  (1.0D0/3.0D0)*sqrt(5.0D0+2.0D0*sqrt(10.0D0/7.0D0))
     weights(1)   = (322.0D0-13.0D0*sqrt(70.0D0))/900.0D0
     weights(2)   = (322.0D0+13.0D0*sqrt(70.0D0))/900.0D0
     weights(3)   = 128.0D0/225.0D0
     weights(4)   = (322.0D0+13.0D0*sqrt(70.0D0))/900.0D0
     weights(5)   = (322.0D0-13.0D0*sqrt(70.0D0))/900.0D0
  case default
     write(*,*) 'n out of range in glwp of module gll. n=',n
     write(*,*) '0<n<5'
     stop
  end select

end subroutine gauss_points

!------------------------------------------------------------------------------
! FUNCTION SIGNUM
!
! Description:
!   Gives the sign of the given real number.
!------------------------------------------------------------------------------
  function signum(x)
    implicit none

    real (kind=real_kind) :: signum
    real (kind=real_kind) :: x

    IF (x > 0.0D0) THEN
      signum = 1.0D0
    ELSEIF (x < 0.0D0) THEN
      signum = -1.0D0
    ELSE
      signum = 0.0D0
    ENDIF
  end function
  
!------------------------------------------------------------------------------
! FUNCTION SIGNUM_FUZZY
!
! Description:
!   Gives the sign of the given real number, returning zero if x is within 
!     a small amount from zero.
!------------------------------------------------------------------------------
  function signum_fuzzy(x)
    implicit none

    real (kind=real_kind) :: signum_fuzzy
    real (kind=real_kind) :: x

    IF (x > fuzzy_width) THEN
      signum_fuzzy = 1.0D0
    ELSEIF (x < fuzzy_width) THEN
      signum_fuzzy = -1.0D0
    ELSE
      signum_fuzzy = 0.0D0
    ENDIF
  end function

  function fuzzy(x,epsilon)
    implicit none

    integer (kind=int_kind) :: fuzzy
    real (kind=real_kind), intent(in) :: epsilon
    real (kind=real_kind) :: x

    IF (ABS(x)<epsilon) THEN
      fuzzy = 0
    ELSE IF (x >epsilon) THEN
      fuzzy = 1
    ELSE !IF (x < fuzzy_width) THEN
      fuzzy = -1
    ENDIF
  end function

!
! see, e.g., http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline2d/
!
subroutine check_lines_cross(x1,x2,x3,x4,y1,y2,y3,y4,lcross)
  implicit none
  real (kind=real_kind), INTENT(IN) :: x1,x2,x3,x4,y1,y2,y3,y4
  LOGICAL, INTENT(OUT) :: lcross
  !
  ! local workspace
  !
  real (kind=real_kind)    :: cp,tx,ty

  cp = (y4-y3)*(x2-x1)-(x4-x3)*(y2-y1)
  IF (ABS(cp)<tiny) THEN
    !
    ! lines are parallel
    !
    lcross = .FALSE.
!      WRITE(*,*) "lines parallel"
  ELSE
    tx = ((x4-x3)*(y1-y3)-(y4-y3)*(x1-x3))/cp
    ty = ((x2-x1)*(y1-y3)-(y2-y1)*(x1-x3))/cp
    IF (tx>-tiny.AND.tx<1.0D0+tiny.AND.&
        ty>-tiny.AND.ty<1.0D0+tiny) THEN
      lcross = .TRUE.
    ELSE
      lcross = .FALSE.
!        WRITE(*,*) "not parallel but not crossing,",tx,ty
    ENDIF
  ENDIF
end subroutine check_lines_cross


  REAL (KIND=dbl_kind) FUNCTION I_00(x_in,y_in)
    IMPLICIT NONE
    REAL (KIND=dbl_kind), INTENT(IN) :: x_in,y_in
    REAL (KIND=dbl_kind) :: x,y

    x = x_in/aa
    y = y_in/aa
!    x = x_in
!    y = y_in
    I_00 = ATAN(x*y/SQRT(one+x*x+y*y))
  END FUNCTION I_00

  REAL (KIND=dbl_kind) FUNCTION I_10(x_in,y_in)
    IMPLICIT NONE
    REAL (KIND=dbl_kind), INTENT(IN) :: x_in,y_in
    REAL (KIND=dbl_kind) :: x,y,tmp

    x = x_in/aa
    y = y_in/aa
    tmp = ATAN(x)
    I_10 = -ASINH(y*COS(tmp))
    !
    ! = -arcsinh(y/sqrt(1+x^2))
    !
  END FUNCTION I_10

  REAL (KIND=dbl_kind) FUNCTION I_10_ab(alpha,beta)
    IMPLICIT NONE
    REAL (KIND=dbl_kind), INTENT(IN) :: alpha,beta
    I_10_ab = -ASINH(COS(alpha) * TAN(beta))
  END FUNCTION I_10_AB

  REAL (KIND=dbl_kind) FUNCTION I_01(x_in,y_in)
    IMPLICIT NONE
    REAL (KIND=dbl_kind), INTENT(IN) :: x_in,y_in
    REAL (KIND=dbl_kind) :: x,y!,beta

    x = x_in/aa
    y = y_in/aa
!    beta = ATAN(y)
!    I_01 = -ASINH(x*COS(beta))
    I_01 = -ASINH(x/SQRT(1+y*y))
  END FUNCTION I_01

  REAL (KIND=dbl_kind) FUNCTION I_01_ab(alpha,beta)
    IMPLICIT NONE
    REAL (KIND=dbl_kind), INTENT(IN) :: alpha,beta
    I_01_ab = -ASINH(COS(beta) * TAN(alpha))
  END FUNCTION I_01_AB

  REAL (KIND=dbl_kind) FUNCTION I_20(x_in,y_in)
    IMPLICIT NONE
    REAL (KIND=dbl_kind), INTENT(IN) :: x_in,y_in
    REAL (KIND=dbl_kind) :: x,y, tmp!,alpha,beta

    x = x_in/aa
    y = y_in/aa
!    alpha = aa*ATAN(x)
!    beta  = aa*ATAN(y)

    tmp = one+y*y

!    I_20 = y*ASINH(COS(beta)*x)+ACOS(SIN(alpha)*SIN(beta))
    I_20 = y*ASINH(x/SQRT(tmp))+ACOS(x*y/(SQRT((one+x*x)*tmp)))
  END FUNCTION I_20

  REAL (KIND=dbl_kind) FUNCTION I_20_ab(alpha,beta)
    IMPLICIT NONE
    REAL (KIND=dbl_kind), INTENT(IN) :: alpha,beta

    I_20_ab = TAN(beta)*ASINH(COS(beta)*TAN(alpha))+ACOS(SIN(alpha)*SIN(beta))
  END FUNCTION I_20_AB

  REAL (KIND=dbl_kind) FUNCTION I_02(x_in,y_in)
    IMPLICIT NONE
    REAL (KIND=dbl_kind), INTENT(IN) :: x_in,y_in
    REAL (KIND=dbl_kind) :: x,y, tmp!,alpha,beta

    x = x_in/aa
    y = y_in/aa
!    alpha = aa*ATAN(x)
!    beta  = aa*ATAN(y)

    tmp=one+x*x

    I_02 = x*ASINH(y/SQRT(tmp))+ACOS(x*y/SQRT(tmp*(1+y*y)))
  END FUNCTION I_02

  REAL (KIND=dbl_kind) FUNCTION I_02_ab(alpha,beta)
    IMPLICIT NONE
    REAL (KIND=dbl_kind), INTENT(IN) :: alpha,beta

    I_02_ab = TAN(alpha)*ASINH(TAN(beta)*COS(alpha))+ACOS(SIN(alpha)*SIN(beta))
  END FUNCTION I_02_AB


  REAL (KIND=dbl_kind) FUNCTION I_11(x_in,y_in)
    IMPLICIT NONE
    REAL (KIND=dbl_kind), INTENT(IN) :: x_in,y_in
    REAL (KIND=dbl_kind) :: x,y

    x = x_in/aa
    y = y_in/aa

    I_11 = -SQRT(1+x*x+y*y)
  END FUNCTION I_11

  REAL (KIND=dbl_kind) FUNCTION I_11_ab(alpha,beta)
    IMPLICIT NONE
    REAL (KIND=dbl_kind), INTENT(IN) :: alpha,beta

    I_11_ab = -SQRT(one+TAN(alpha)**2+TAN(beta)**2)
  END FUNCTION I_11_AB
!------------------------------------------------------------------------------
! FUNCTION ASINH
!
! Description:
!   Hyperbolic arcsin function
!------------------------------------------------------------------------------
  FUNCTION ASINH(x)
    IMPLICIT NONE

    REAL (KIND=dbl_kind) :: ASINH
    REAL (KIND=dbl_kind) :: x

    ASINH = LOG(x + SQRT(x * x + one))
  END FUNCTION


  !********************************************************************************
  !
  ! Gauss-Legendre quadrature
  !
  ! Tabulated values
  !
  !********************************************************************************
  SUBROUTINE glwp(n,weights,abscissae)
    IMPLICIT NONE
    REAL (KIND=dbl_kind), DIMENSION(n), INTENT(OUT) :: weights, abscissae
    INTEGER (KIND=int_kind)           , INTENT(IN ) :: n

    SELECT CASE (n)
    CASE(1)
       abscissae(1) = 0.0
       weights(1)   = 2.0
    CASE(2)
       abscissae(1) = -SQRT(1.0/3.0)
       abscissae(2) = SQRT(1.0/3.0)
       weights(1)   =  1.0
       weights(2)   =  1.0
    CASE(3)
       abscissae(1) = -0.774596669241483377035853079956_dbl_kind
       abscissae(2) =  0.0
       abscissae(3) =  0.774596669241483377035853079956_dbl_kind
       weights(1)   =  0.555555555555555555555555555556_dbl_kind
       weights(2)   =  0.888888888888888888888888888889_dbl_kind
       weights(3)   =  0.555555555555555555555555555556_dbl_kind
    CASE(4)
       abscissae(1) = -0.861136311594052575223946488893_dbl_kind
       abscissae(2) = -0.339981043584856264802665659103_dbl_kind
       abscissae(3) =  0.339981043584856264802665659103_dbl_kind
       abscissae(4) =  0.861136311594052575223946488893_dbl_kind
       weights(1)   =  0.347854845137453857373063949222_dbl_kind
       weights(2)   =  0.652145154862546142626936050778_dbl_kind
       weights(3)   =  0.652145154862546142626936050778_dbl_kind
       weights(4)   =  0.347854845137453857373063949222_dbl_kind
    CASE(5)
       abscissae(1) = -(1.0/3.0)*SQRT(5.0+2.0*SQRT(10.0/7.0))
       abscissae(2) = -(1.0/3.0)*SQRT(5.0-2.0*SQRT(10.0/7.0))
       abscissae(3) =  0.0
       abscissae(4) =  (1.0/3.0)*SQRT(5.0-2.0*SQRT(10.0/7.0))
       abscissae(5) =  (1.0/3.0)*SQRT(5.0+2.0*SQRT(10.0/7.0))
       weights(1)   = (322.0_dbl_kind-13.0_dbl_kind*SQRT(70.0_dbl_kind))/900.0_dbl_kind
       weights(2)   = (322.0_dbl_kind+13.0_dbl_kind*SQRT(70.0_dbl_kind))/900.0_dbl_kind
       weights(3)   = 128.0_dbl_kind/225.0_dbl_kind
       weights(4)   = (322.0_dbl_kind+13.0_dbl_kind*SQRT(70.0_dbl_kind))/900.0_dbl_kind
       weights(5)   = (322.0_dbl_kind-13.0_dbl_kind*SQRT(70.0_dbl_kind))/900.0_dbl_kind
    CASE DEFAULT
       WRITE(*,*) 'n out of range in glwp of module gll. n=',n
       WRITE(*,*) '0<n<5'
       STOP
    END SELECT

  END SUBROUTINE glwp


END MODULE remap
