module interp_tools

use chimere_params, only : iprec 

!--------------------------------------------------------------------------------------
! Subroutines for the interpolation management1
!--------------------------------------------------------------------------------------
! calculate_corners
! calculate_centers
! findcell
! wcal
! is_inside
! test_rectilon
! test_rectilat
! calcorners
! det
! do_uv_rot_sm: rotate WRF winds (from grid to earth system north-west)
! itpconc : Linear interpolation for concentrations (ensure mass conservation)
! intersec : returns the size of intersection of [mina;maxa] with [minb;maxb]
! cuminbins: Cumul of aerosols concentrations
! distrib_lognorm
! horint
! modisec
! lintp: Linear interpolation in pressure coordinates
! lint1
!--------------------------------------------------------------------------------------

implicit none

contains

subroutine calculate_bound(nxg,nyg,corners_3d,xmin,xmax,ymin,ymax)
! 
! Given a 3d corner array, this subroutine does some trigonometry to calculate the extreme latitudes and longitudes  
! of the mesh, with latitude in degrees N between -90. and +90. (no suprise here)
! and longitude in degrees E, by convention between -180. and 180. (the same convention as in most datasets)
! 
! INPUT variables
integer,intent(in)    :: nxg,nyg                    ! Dimensions of the input grid
real,intent(in)       :: corners_3d(nxg,nyg,4,3)   ! Grid corners of the input grid, in 3d coordinates
real,parameter        :: localpi   = 3.1415926535898d0  ! Redefined not to depend on chimere_consts here
real(kind=8),intent(out) :: xmin,xmax,ymin,ymax   ! ymin and ymax are bw -90. and 90, xmin and xmax bw -180. and 180.
integer :: ix,iy,ic
real(kind=8)          :: x,y

xmin=360d0
xmax=-360d0
ymin=90d0
ymax=-90d0
do ic=1,4
   do iy=1,nyg
      do ix=1,nxg
         y=asin(corners_3d(ix,iy,ic,3))*180d0/localpi
         if(corners_3d(ix,iy,ic,1).eq.0d0) then
           if(corners_3d(ix,iy,ic,2).lt.0) then
             x=-90.
           else
             x=90.
           endif
         else
	   ! This gives a value between -90. and 90. by definition of atan
           x=atan(corners_3d(ix,iy,ic,2)/corners_3d(ix,iy,ic,1))*180d0/localpi          
           if(corners_3d(ix,iy,ic,1).lt.0d0)then   
               if(x.lt.0.)then
                   x=x+180.
               else
                   x=x-180.
               endif
           endif
         endif
         ymin=min(ymin,y)
         ymax=max(ymax,y)
         xmin=min(xmin,x)
         xmax=max(xmax,x)
      enddo
   enddo
enddo

end subroutine calculate_bound

!--------------------------------------------------------------------------------------

subroutine calculate_corners(xg,yg,corners_3d,nxg,nyg,is_global,is_periodic)

! Given grid centers (arrays xg, yg) this subroutine outputs grid corners in corners_3d. 
! The coordinates of the corners are given as unit vectors in cartesian coordinates with axis 1 at 
! the Equator and going through Greenwich Meridian, Axis 2 at the Equator 
! and 90E, and Axis 3 is the polar axis
! These grid corners are checked at the end of the routine by calling findcell on the grid centers and checking 
! that each cell center falls within the respective 4 corners. En error message is produce if this check fails.

! INPUT variables
integer,intent(in)    :: nxg,nyg                    ! Dimensions of the input grid
real,intent(in)       :: xg(nxg,nyg),yg(nxg,nyg)    ! Lat and lon of the grid centers of the input grid, in degrees
logical,optional,intent(in)            :: is_global
logical,optional,intent(in)            :: is_periodic


! OUTPUT variables
real,intent(out)       :: corners_3d(nxg,nyg,4,3)   ! Grid corners of the input grid, in 3d coordinates

! LOCAL variables
real,parameter        :: localpi   = 3.1415926535898d0  ! Redefined not to depend on chimere_consts here
real                  :: xg_3d(nxg,nyg,3)
real,dimension(3)     :: v 
real                  :: clat,clon,slat,slon
integer               :: ix,iy,ixf,iyf,lstix,lstiy
logical               :: found
logical               :: globaldom,periodicdom
integer,parameter     :: idebug=0  ! perform verification. 
                                   ! Very time-consuming for big domains.

if(present(is_global))then
  globaldom=is_global
else
  globaldom=.false.
endif

if(present(is_periodic))then
  periodicdom=is_periodic
else
  periodicdom=.false.
endif

! First calculate the unit vectors corresponding to the input grid
do ix=1,nxg
do iy=1,nyg
   clat=cos(yg(ix,iy)/180.*localpi)
   slat=sin(yg(ix,iy)/180.*localpi)
   clon=cos(xg(ix,iy)/180.*localpi)
   slon=sin(xg(ix,iy)/180.*localpi)
   xg_3d(ix,iy,1)=clat*clon
   xg_3d(ix,iy,2)=clat*slon
   xg_3d(ix,iy,3)=slat
enddo
enddo

! 1. Treat lower left corner for all grid cells

do ix=2,nxg
  do iy=2,nyg
       call normalize(0.5*(xg_3d(ix,iy,:)+xg_3d(ix-1,iy-1,:)),v)
       corners_3d(ix,iy,1,:)=v
  enddo
enddo

! left edge
do iy=2,nyg
       call normalize(corners_3d(2,iy,1,:)+xg_3d(1,iy,:)-xg_3d(2,iy,:),v)
       corners_3d(1,iy,1,:)=v
enddo

!bottom edge
do ix=1,nxg
       call normalize(corners_3d(ix,2,1,:)+xg_3d(ix,1,:)-xg_3d(ix,2,:),v)
       corners_3d(ix,1,1,:)=v
enddo

!2. Treat upper right corner for all grid cells
do ix=1,nxg-1
   do iy=1,nyg-1
      corners_3d(ix,iy,3,:)=corners_3d(ix+1,iy+1,1,:)
   enddo
enddo

! right edge
do iy=1,nyg-1
      call normalize(corners_3d(nxg-1,iy,3,:)+xg_3d(nxg,iy,:)-xg_3d(nxg-1,iy,:),v)
      corners_3d(nxg,iy,3,:)=v
enddo    


! upper edge
do ix=1,nxg
      call normalize(corners_3d(ix,nyg-1,3,:)+xg_3d(ix,nyg,:)-xg_3d(ix,nyg-1,:),v)
      corners_3d(ix,nyg,3,:)=v
enddo    

!3. Treat upper left corner for all grid cells
do ix=1,nxg
   do iy=1,nyg-1
      corners_3d(ix,iy,4,:)=corners_3d(ix,iy+1,1,:)
   enddo
enddo

! upper edge
do ix=2,nxg
      corners_3d(ix,nyg,4,:)=corners_3d(ix-1,nyg,3,:)
enddo

!ul corner
      call normalize(corners_3d(1,nyg-1,4,:)+xg_3d(1,nyg,:)-xg_3d(1,nyg-1,:),v)
      corners_3d(1,nyg,4,:)=v

!4. Treat lower-right corner for all grid cells
do ix=1,nxg-1
  do iy=1,nyg
      corners_3d(ix,iy,2,:)=corners_3d(ix+1,iy,1,:)
  enddo
enddo

! right edge
do iy=2,nyg
      corners_3d(nxg,iy,2,:)=corners_3d(nxg,iy-1,3,:)
enddo

!lr corner
      call normalize(corners_3d(nxg-1,1,2,:)+xg_3d(nxg,1,:)-xg_3d(nxg-1,1,:),v)
      corners_3d(nxg,1,2,:)=v

! Special treatment of left and right edge for periodic domains
if(periodicdom)then
  do iy=2,nyg
     call normalize(0.5*(xg_3d(1,iy,:)+xg_3d(nxg,iy-1,:)),v)
     corners_3d(1,iy,1,:)=v
  enddo
  do iy=1,nyg-1
     call normalize(0.5*(xg_3d(1,iy,:)+xg_3d(nxg,iy+1,:)),v)
     corners_3d(1,iy,4,:)=v
  enddo
  call normalize(corners_3d(1,2,1,:)+xg_3d(1,1,:)-xg_3d(1,2,:),v)
  corners_3d(1,1,1,:)=v

  call normalize(corners_3d(1,nyg-1,4,:)+xg_3d(1,nyg,:)-xg_3d(1,nyg-1,:),v)
  corners_3d(1,nyg,4,:)=v

  do iy=1,nyg
      corners_3d(nxg,iy,3,:)=corners_3d(1,iy,4,:)
      corners_3d(nxg,iy,2,:)=corners_3d(1,iy,1,:)
  enddo
endif

if(idebug==1)then

!Check that findcell works accurately with the current corners
ixf=0
iyf=0
do ix=1,nxg
! print*,ix,'/',nxg
do iy=1,nyg
    found=.false.
    call findcell(xg(ix,iy),yg(ix,iy),ixf,iyf,nxg,nyg,corners_3d,lstix,lstiy,found,globaldom,periodicdom)
    if(ix.ne.ixf.or.iy.ne.iyf)then
     if(.not.(globaldom))then
      print*,'WARNING, calculate_corners_failure, SORRY'
      print*,'Problem with regional domain - ix,iy,ixf,iyf,found,nxg,nyg:',ix,iy,ixf,iyf,found,nxg,nyg
       print*,'xg(ix,iy),yg(ix,iy):',xg(ix,iy),yg(ix,iy)
       print*,'corners :'
       print*,'ll',corners_3d(ix,iy,1,:)
       print*,'lr',corners_3d(ix,iy,2,:)
       print*,'ur',corners_3d(ix,iy,3,:)
       print*,'ul',corners_3d(ix,iy,4,:)
       print*,'expected values',xg_3d(ix,iy,:)
       print*,'periodic domain ?',periodicdom
       print*,'global domain ?',globaldom
       print*,xg(ixf,iyf),yg(ixf,iyf)       
       call exit1('Error 1 in calculatecorners')
     else
      if((ix.ne.1).and.(ix.ne.nxg).and.(iy.ne.1).and.(iy.ne.nyg))then
       print*,'WARNING, calculate_corners_failure, SORRY'
       print*,'Problem with global domain',ix,iy,ixf,iyf,found
       print*,xg(ix,iy),yg(ix,iy),xg(ixf,iyf),yg(ixf,iyf)
       call exit1('Error 2 in calculatecorners')
      endif
     endif
    endif
enddo
enddo
endif

end subroutine calculate_corners

!===========================================================================================

subroutine calculate_centers(xg,yg,centers_3d,ulambda,uphi,nxg,nyg)

! Given grid centers (arrays xg, yg) this subroutine outputs grid centers in corners_3d.
! The coordinates of the corners are given as unit vectors in cartesian coordinates with axis 1 at
! the Equator and going through Greenwich Meridian, Axis 2 at the Equator
! and 90E, and Axis 3 is the polar axis
! ulambda and uphi are the standard unit vectors


! INPUT variables
! Dimensions of the input grid
integer,intent(in)    :: nxg,nyg
! Lat and lon of the grid centers of the input grid, in degrees
real(kind=iprec),intent(in)       :: xg(nxg,nyg),yg(nxg,nyg)

! OUTPUT variables
! Grid centers of the input grid, in 3d coordinates
real(kind=iprec),intent(out)       :: centers_3d(nxg,nyg,3)
real(kind=iprec),intent(out)       :: ulambda(nxg,nyg,3)
real(kind=iprec),intent(out)       :: uphi(nxg,nyg,3)


! LOCAL variables
real(kind=iprec),parameter        :: localpi   = 3.1415926535898d0  ! Redefined not to depend on chimere_consts here
real(kind=iprec)                  :: clat,clon,slat,slon
integer               :: ix,iy

! calculate the unit vectors corresponding to the input grid
do ix=1,nxg
do iy=1,nyg
   clat=cos(yg(ix,iy)/180.*localpi)
   slat=sin(yg(ix,iy)/180.*localpi)
   clon=cos(xg(ix,iy)/180.*localpi)
   slon=sin(xg(ix,iy)/180.*localpi)
   centers_3d(ix,iy,1)=clat*clon
   centers_3d(ix,iy,2)=clat*slon
   centers_3d(ix,iy,3)=slat
   ulambda(ix,iy,1)=-slon
   ulambda(ix,iy,2)= clon
   ulambda(ix,iy,3)=0
   uphi(ix,iy,1)=-slat * clon
   uphi(ix,iy,2)=-slat * slon
   uphi(ix,iy,3)= clat
enddo
enddo

end subroutine calculate_centers

!==============================================================================================
  subroutine findcell(x,y,ix,iy,nxg,nyg,xg3d,lstix,lstiy,found,is_global,is_periodic)
  ! The 3d coordinates of the 4 corners are stored in xg3d
  !                      A :        xg3d(ix,iy,1,:)
  !                      B :        xg3d(ix,iy,2,:)
  !                      C :        xg3d(ix,iy,3,:)
  !                      D :        xg3d(ix,iy,4,:)
  ! The corner coordinates shall be produced by calculate_corners, 
  ! from the target grid center. 3d coordinates are used to avoid problems at the poles.
  ! 
  ! It is assumed that the (ABCD) polygon is numbered in the trigonometric
  ! direction, when looking to the polygon from oustide of the Earth, and towards the Earth center.
  ! 
  ! This is ensured if calculate_corners is called on a grid where grid cell centers are ordered as follows,
  ! which is the historical order in CHIMERE for lat-lon grids :
  !  (1,nyg)----------------->(nxg,nyg)
  !           .....
  !  (1,2  )------------------->(nxg,2  )
  !  (1,1  )------------------->(nxg,1  )
  !
  ! x,y    (input)  : long and lat ot target point M  (input)
  ! ix,iy  (output) : Indices of the grid cell containing (x,y) - i.e., ix and iy such as the point with 
  !          geographical coordinates (x,y) is inside the polygone defined by xg3d(ix,iy,:,:,) 
  ! nxg,nyg (input) : Size of xg3d
  ! xg3d (input)    : Grid corners (see above)
  ! lstix, lstiy(inout) : Initial guess for the location of the target point. 
  ! If target point is found, lstix and lstiy are updated at the end of the routine.
  ! found (output)      : Logical, true of findcell is successfull, false otherwise

  ! The algorithm implemented here is known as "stencil jumping" adapted to 
  ! spherical geometry (determinant instead of cross products).

  

    ! subroutine arguments
! INPUT
    real,intent(in)                        :: x,y
    real,dimension(nxg,nyg,4,3),intent(in) :: xg3d
    integer,intent(in)                     :: nxg,nyg
    logical,optional,intent(in)            :: is_global
    logical,optional,intent(in)            :: is_periodic

! OUTPUT
    integer,intent(out)                    :: ix,iy
    logical,intent(out)   :: found

! INOUT
    integer,intent(inout) :: lstix
    integer,intent(inout) :: lstiy


! LOCAL variables
    real,dimension(3)     :: targ,center
    real                  :: dist,distmin
    real                  :: cp,clat,clon,slat,slon
    integer               :: count ! Counter to avoid infinite loop
    integer               :: iymin,ixmin
    logical               :: inside
    character(len=1)      :: lstmove
    real,parameter        :: localpi   = 3.1415926535898d0  ! Redefined not to importe chimere_consts
    logical               :: globaldom,periodicdom

   cp = -1.
   count = 0

   if(present(is_global))then
     globaldom=is_global
   else
     globaldom=.false.
   endif

   if(present(is_periodic))then
     periodicdom=is_periodic
   else
     periodicdom=.false.
   endif

! Compute the 2d coordinates of target point
   clat=cos(y/180.*localpi)
   slat=sin(y/180.*localpi)
   clon=cos(x/180.*localpi)
   slon=sin(x/180.*localpi)
   targ(1)=clat*clon
   targ(2)=clat*slon
   targ(3)=slat



   inside=.false.

! Initialisation
   if(((lstix.gt.0).and.(lstix.le.nxg)).and.((lstiy.gt.0).and.(lstiy.le.nyg)))then
     ix = lstix
     iy = lstiy
   else
     ix=1
     iy=1
   endif
   lstmove='X'


   do while(.not.(inside))
   if(.not.(((ix.gt.0).and.(ix.le.nxg)).and.((iy.gt.0).and.(iy.le.nyg))))then
     print*,'ix iy',ix,iy
     print*,'nxg nyg',nxg,nyg
   endif
   inside=.true.
   count = count +1
   if (count.gt.(2 * (nxg + nyg + 1)))then
      if(globaldom)then
        exit    !If the domain is a global domain, we need to force success of the procedure
      else
!        print*,ix,iy,x,y,'WARNING INFINITE LOOP ACCURED IN FINDCELL'
        found=.false.
        ix=0
        iy=0
        return
      endif
   endif

! Compute the determinant of (OA,OB,OM)

   cp = det(xg3d(ix,iy,1,:),xg3d(ix,iy,2,:),targ)
   inside=inside.and.(cp.ge.0.)
   if((cp.lt.0.).and.(iy.ge.2.or.globaldom).and.(lstmove.ne.'n'))then
     iy = iy - 1 ! Go one step further south
     if(iy.eq.0.and.globaldom)then
        iy=1
        ix=ix+nxg/2
        if(ix.gt.nxg)ix=ix-nxg
     endif
    lstmove = 's'
   else
     ! compute MB x BC
     cp = det(xg3d(ix,iy,2,:),xg3d(ix,iy,3,:),targ)
     inside=inside.and.(cp.ge.0.)
     if((cp.lt.0.).and.(ix.le.(nxg-1).or.globaldom).and.(lstmove.ne.'w'))then
       ix = ix + 1 ! Go one step further east
       if ((globaldom.or.periodicdom).and.(ix.gt.nxg))ix=1
       lstmove = 'e'
     else
       ! compute MC x CD
       cp = det(xg3d(ix,iy,3,:),xg3d(ix,iy,4,:),targ)
       inside=inside.and.(cp.ge.0.)
       if((cp.lt.0.).and.(iy.le.(nyg-1).or.globaldom).and.(lstmove.ne.'s'))then
         iy = iy + 1 ! Go one step further north
         if ((globaldom).and.(iy.gt.nyg))then 
               iy=nyg
               ix=ix+nxg/2
               if(ix.gt.nxg)ix=ix-nxg
         endif
         lstmove = 'n'
       else
         ! compute MD x DA
         cp = det(xg3d(ix,iy,4,:),xg3d(ix,iy,1,:),targ)
         inside=inside.and.(cp.ge.0.)
         if((cp.lt.0.).and.(lstmove.ne.'e'))then
           ix = ix - 1 ! Go one step further west
           lstmove = 'w'
           if(ix.le.0)then
            if(.not.(globaldom.or.periodicdom))then
             ix = 0
             iy = 0
             found=.false.
             return
            else
             ix=nxg
            endif
           endif
         else
           if(inside)then
            if(globaldom)then
             if(scalar(xg3d(ix,iy,4,:),targ).lt.0.)then !Problem with antipodes
               ix=ix+nxg/2
               if(ix.gt.nxg)ix=ix-nxg
               iy=nyg+1-iy
               lstmove='X'
               cycle
             endif
            endif
            found=.true.
            lstix=ix
            lstiy=iy
            return
           else
            if(globaldom)then
             exit    !Force emergency find
            else
             ix=0
             iy=0
             found=.false.
             return
            endif
           endif
         endif ! end test on MD x DA
       endif ! end test on MC x CD
     endif ! end test on MBxBC
   endif ! end test on MA x AB
   end do ! end research loop

! Emergency find for global domains
   distmin=2.d0 ! Max distance bw two points on a unit sphere
   do ix=1,nxg
   do iy=1,nyg
      center=0.25*xg3d(ix,iy,1,:) &
&           +0.25*xg3d(ix,iy,2,:) &
&           +0.25*xg3d(ix,iy,3,:) &
&           +0.25*xg3d(ix,iy,4,:) 
      dist=norme(targ-center)
      if(dist.lt.distmin)then
        distmin=dist
        ixmin=ix
        iymin=iy
      endif
   enddo
   enddo 
   found=.true.
   ix=ixmin
   iy=iymin
   lstix=ix
   lstiy=iy
! End of emergency find for global domains
  end subroutine findcell

!==============================================================================================

subroutine wcal(x,y,ix,iy,nxg,nyg,xg_corners,wx,wy)

! This routine computes the coefficients wx, wy of the bilinear
! interpolation used later on to interpolate the value at point M of a field defined
! at points 1,2,3,4.

! inputs :
!     x, y   : longitude and latitude of M point
!    ix,iy   : the coordinates (indexes) of the cell which contains M
!    nxg,nyg : grid dimensions (number of corners, not number of cells)
!    xg,yg   : longitudes and latitudes of the grid points
! outputs :
!    wx,wy   : the weights of the bilinear interpolation

! Remark :
! The interpolation formula used later on is equivalent to
! xout=(1-wx)*xin(i) + wx*xin(i+1)
! it expects wx=0 for two coincident points
! Thus the "real" weight is not wx, but (1-wx)
! but for historical reasons, we keep calling wx what should be called (1-wx)

! subroutine arguments
integer,intent(in) :: ix, iy, nxg, nyg
real(kind=iprec),intent(in) :: x,y
real(kind=iprec),dimension(nxg,nyg,4,3),intent(in) :: xg_corners
real(kind=iprec),intent(out) :: wx, wy

! parameters
integer,parameter :: ntry=100

! local variables
real(kind=iprec),dimension(3) :: xx
real(kind=iprec),dimension(3) :: x3d
real(kind=iprec)	      :: clat,slat,clon,slon
real(kind=iprec) :: d2min,d2,b1,b0,b,a1,a0,a
integer :: ntx,nty

! Redefined not to depend on chimere_consts here
real(kind=iprec),parameter :: localpi   = 3.1415926535898d0  

clat=cos(y/180.*localpi)
slat=sin(y/180.*localpi)
clon=cos(x/180.*localpi)
slon=sin(x/180.*localpi)
x3d(1)=clat*clon
x3d(2)=clat*slon
x3d(3)=slat


! Search for weights
d2min = 1e20
do ntx=0,ntry
   a0 = float(ntx)/ntry
   a1 = 1. - a0
   do nty=0,ntry
      b0 = float(nty)/ntry
      b1 = 1. - b0
      xx =     a0*b0*xg_corners(ix,iy,3,:) &
           & + a1*b0*xg_corners(ix,iy,4,:) &
           & + a0*b1*xg_corners(ix,iy,2,:) &
           & + a1*b1*xg_corners(ix,iy,1,:)
      d2 = scalar(x3d-xx,x3d-xx)
      if(d2.lt.d2min) then
         d2min = d2
         a  = a0
         b  = b0
      endif
   enddo
enddo
wx=real(a)
wy=real(b)

end subroutine wcal


!//////////////////////////////////////////////////////////////////////////////////////////

  function is_inside(x,y,xg,yg)
    ! checks if the point M(x,y) is inside the polygon (ix,iy),(ix+1,iy),(ix,iy+1),(ix+1,iy+1)

    ! definitions
    ! ===========
    !
    ! A cell is defined by its 4 corners 1,2,3,4 which correspond to the index
    ! of xg and yg
    ! if M is on the 1-2 segment, it belongs to the cell
    ! if M is on the 2-3 segment, it belongs to the cell
    ! if M is on the 3-4 segment, it does not belong to the cell
    ! if M is on the 4-1 segment, it does not belong to the cell
    !
    ! if M is at corner 2,     it belongs to the cell
    ! if M is at corners 1,3,4 it does not belong to the cell

    ! arguments
    ! =========
    !
    ! input:
    !    x,y     : longitude and latitude of M
    !    xg,yg   : longitudes and latitudes of the four corners of the grid cell
    !              ordering of corners : walking along the perimeter, direction is not
    !              significant

    ! properties
    ! ==========
    !
    ! if you pave and area with such polygons:
    !     * any point of this area will belong
    !       to one polygon and only one

    ! returns .true. if M is inside


    real(kind=iprec),intent(in),dimension(4) :: xg,yg
    real(kind=iprec),intent(in) :: x,y
    logical :: is_inside

    real(kind=iprec) :: a12,a23,a34,a41

    a12=(x-xg(1))*(y-yg(2))-(x-xg(2))*(y-yg(1))
    a23=(x-xg(2))*(y-yg(3))-(x-xg(3))*(y-yg(2))
     is_inside=.false.
    if((a12*a23).ge.0.)then
      a34=(x-xg(3))*(y-yg(4))-(x-xg(4))*(y-yg(3))
      if((a23*a34).ge.0.)then
       a41=(x-xg(4))*(y-yg(1))-(x-xg(1))*(y-yg(4))
       if((a34*a41).ge.0.)then
          if((a41*a12).ge.0.)then
            if((a41.ne.0.).and.(a34.ne.0.))is_inside=.true.
          endif
        endif
      endif
    endif

  end function is_inside

!//////////////////////////////////////////////////////////////////////////////////////////

function test_rectilon(xlong,nwe,nsn)
!-------------------------------------------------------
real(kind=iprec),dimension(nwe,nsn),intent(in) :: xlong
integer,intent(in) :: nwe,nsn
logical :: test_rectilon
real(kind=iprec) :: deltalon
real(kind=iprec) :: epsx,eps1,eps2,eps3,eps4,epslin
integer :: iwe
!-------------------------------------------------------

deltalon = (xlong(nwe,nsn/2) - xlong(1,nsn/2))/(nwe-1)
epsx = abs(deltalon)/10.

! testing parallelism of Y sides with LAT axis
eps1=abs(xlong(1,nsn)-xlong(1,1))
eps2=abs(xlong(1,nsn/2)-xlong(1,1))
eps3=abs(xlong(nwe,nsn)-xlong(nwe,1))
eps4=abs(xlong(nwe,nsn/2)-xlong(nwe,1))
test_rectilon = (eps1<epsx).and.(eps2<epsx).and.(eps3<epsx).and.(eps4<epsx)

! testing linearity
do iwe=1,nwe
   epslin = abs(xlong(iwe,nsn/2) -  xlong(1,nsn/2) - (iwe-1)*deltalon)
   test_rectilon = test_rectilon .and. (epslin<epsx)
enddo

!-------------------------------------------------------
end function test_rectilon


!//////////////////////////////////////////////////////////////////////////////////////////

function test_rectilat(xlat,nwe,nsn)
!-------------------------------------------------------
real(kind=iprec),dimension(nwe,nsn),intent(in) :: xlat
integer,intent(in) :: nwe,nsn
logical :: test_rectilat
real(kind=iprec) :: deltalat
real(kind=iprec) :: epsx,eps1,eps2,eps3,eps4,epslin
integer :: isn
!-------------------------------------------------------

deltalat = (xlat(nwe/2,nsn)-xlat(nwe/2,1))/(nsn-1)
epsx = abs(deltalat)/10.

! testing parallelism of X sides with LON axis
eps1=abs(xlat(nwe,1)-xlat(1,1))
eps2=abs(xlat(nwe/2,1)-xlat(1,1))
eps3=abs(xlat(nwe,nsn)-xlat(1,nsn))
eps4=abs(xlat(nwe/2,nsn)-xlat(1,nsn))
test_rectilat = (eps1<epsx).and.(eps2<epsx).and.(eps3<epsx).and.(eps4<epsx)

! testing linearity
do isn=1,nsn
   epslin = abs(xlat(nwe/2,isn) -  xlat(nwe/2,1) - (isn-1)*deltalat)
   test_rectilat = test_rectilat .and. (epslin<epsx)
enddo

!-------------------------------------------------------
end function test_rectilat

!//////////////////////////////////////////////////////////////////////////////////////////

  subroutine calcorners(nx,ny,xc,yc,xe,ye)

! Converts grid centers into grid corners. This routine in now almost obsolete in CHIMERE,
! It is used only for regular lat-lon grids. No new calls to this routine should be added in CHIMERE
! Please use calculate_corners, which is now the tool producing the appropriate input for findcell.

    ! Arguments
    integer,intent(in) :: nx,ny
    real(kind=iprec),dimension(nx,ny),intent(in) :: xc,yc
    real(kind=iprec),dimension(nx+1,ny+1),intent(out) :: xe,ye

    ! Local variables
    real(kind=iprec),dimension(0:nx+1,0:ny+1) :: x,y
    integer :: i,j

    ! Build a grid with grid centers plus
    ! a perimetral halo of the same resolution
    x(1:nx,1:ny)=xc
    y(1:nx,1:ny)=yc
    do i=1,nx
       x(i,0) = 2*xc(i,1) - xc(i,2)
       y(i,0) = 2*yc(i,1) - yc(i,2)
       x(i,ny+1) = 2*xc(i,ny) - xc(i,ny-1)
       y(i,ny+1) = 2*yc(i,ny) - yc(i,ny-1)
    enddo
    do j=1,ny
       x(0,j) = 2*xc(1,j) - xc(2,j)
       y(0,j) = 2*yc(1,j) - yc(2,j)
       x(nx+1,j) = 2*xc(nx,j) - xc(nx-1,j)
       y(nx+1,j) = 2*yc(nx,j) - yc(nx-1,j)
    enddo
    x(0,0) = 2*xc(1,1) - xc(2,2)
    y(0,0) = 2*yc(1,1) - yc(2,2)
    x(0,ny+1) = 2*xc(1,ny) - xc(2,ny-1)
    y(0,ny+1) = 2*yc(1,ny) - yc(2,ny-1)
    x(nx+1,0) = 2*xc(nx,1) - xc(nx-1,2)
    y(nx+1,0) = 2*yc(nx,1) - yc(nx-1,2)
    x(nx+1,ny+1) = 2*xc(nx,ny) - xc(nx-1,ny-1)
    y(nx+1,ny+1) = 2*yc(nx,ny) - yc(nx-1,ny-1)

    ! Calculate grid corners
    do j=1,ny+1
       do i=1,nx+1
          xe(i,j) = 0.5*(x(i,j)+x(i-1,j-1))
          ye(i,j) = 0.5*(y(i,j)+y(i-1,j-1))
       enddo
    enddo

  end subroutine calcorners

!//////////////////////////////////////////////////////////////////////////////////////////
  function det(x,y,z)
! Computes the determinant of the three 3d vectors x,y,z
! If strictly positive, (x,y,z) forms a repere direct
! If strictly negative, (x,y,z) forms a repere indirect
! If null z is included in the (x,y) plane
  real(kind=iprec),dimension(3)  :: x,y,z
  real(kind=iprec)               :: det

  det= x(1) * (y(2)*z(3) - z(2)*y(3)) &
&     -x(2) * (y(1)*z(3) - z(1)*y(3)) &
&     +x(3) * (y(1)*z(2) - z(1)*y(2))

  end function det

!//////////////////////////////////////////////////////////////////////////////////////////
  subroutine do_uv_rot_sm(wx,wy,wz,wm,cosalpha,sinalpha)

    ! https://forum.mmm.ucar.edu/phpBB3/viewtopic.php?f=78&t=179
    ! Below are the equations to convert model u and v to earth-relative u and
    ! v:
    ! u-earth = u-model*cosalpha - v-model*sinalpha
    ! v-earth = v-model*cosalpha + u-model*sinalpha
    ! where cosalpha and sinalpha can be found in model input and output files
    ! Warning, the array dimensions are expected to be as follows:
    ! if wz and wm have dimensions [nx,ny,nz] then the assumed dimensions for
    ! other arrays are as follows:
    ! wx[nx+1,ny,nz]    - staggered in x
    ! wy[nx,ny+1,nz]    - staggered in y
    ! cosalpha[nx,ny] ; sinalpha[nx,ny]

    real(kind=iprec),dimension(:,:,:),intent(in)  :: wx                    ! 3d array with "zonal" wind (relative to WRF coordinate system)
    real(kind=iprec),dimension(:,:,:),intent(in)  :: wy                    ! 3d array with "meridional" wind (relative to WRF coordinate system)
    real(kind=iprec),dimension(:,:,:),intent(out) :: wz,wm                 ! 3d array with zonal and meridional wind 
    real(kind=iprec),dimension(:,:),intent(in)    :: cosalpha,sinalpha

    integer :: i,j,k,nx,ny,nl
    real(kind=iprec)    :: wxcenter,wycenter

    nx=size(wz,1)
    ny=size(wz,2)
    nl=size(wz,3)

    do j=1,ny
       do i=1,nx
          do k=1,nl
             wxcenter=0.5 * ( wx(i,j,k) + wx(i+1,j,k) )
             wycenter=0.5 * ( wy(i,j,k) + wy(i,j+1,k) )
             wz(i,j,k)= wxcenter * cosalpha(i,j) - wycenter*sinalpha(i,j)
             wm(i,j,k)= wxcenter * sinalpha(i,j) + wycenter*cosalpha(i,j)
          enddo
       enddo
    enddo

  end subroutine do_uv_rot_sm

!//////////////////////////////////////////////////////////////////////////////////////////
  subroutine do_uv_rot_2d(wx,wy,wz,wm,cosalpha,sinalpha)

    ! u10m and v10m need to be rotated.
    !    https://www-k12.atmos.washington.edu/~ovens/wrfwinds.html
    !    They are given on the centered (non-staggered) grid.
    !
    ! Equations to convert model u and v to earth-relative u and v:
    ! u-earth = u-model*cosalpha - v-model*sinalpha
    ! v-earth = v-model*cosalpha + u-model*sinalpha
    ! where cosalpha and sinalpha can be found in model input and output files
    ! if wz and wm have dimensions [nx,ny] then the assumed dimensions for other
    ! arrays are as follows:
    ! wx[nx,ny]    - non-staggered
    ! wy[nx,ny]    - non-staggered
    ! cosalpha[nx,ny] ; sinalpha[nx,ny]

    real(kind=iprec),dimension(:,:),intent(in)  :: wx, wy
    real(kind=iprec),dimension(:,:),intent(out) :: wz, wm
    real(kind=iprec),dimension(:,:),intent(in)    :: cosalpha,sinalpha

    integer :: i,j,k,nx,ny

    nx=size(wz,1)
    ny=size(wz,2)

    do j=1,ny
    do i=1,nx
       wz(i,j)= wx(i,j)*cosalpha(i,j) - wy(i,j)*sinalpha(i,j)
       wm(i,j)= wx(i,j)*sinalpha(i,j) + wy(i,j)*cosalpha(i,j)
    enddo
    enddo

  end subroutine do_uv_rot_2d

!//////////////////////////////////////////////////////////////////////////////////////////

!-----------------------------------------------------------------------
! Linear interpolation for concentrations (ensure mass conservation)
! from raw (xin) to fine (xout) resolution
! we must have: xin[1] > xout[1] and xin[nxin+1] < xout[nout]
!----------------------------------------------------------------------
!    xin      : Xin intervals [i,i+1] where yin[i] is known
!    yin      : Yin table to interpolate
!    nxin     : intervals of values in yin and xin
!    xout     : Xout bins where we want to interpolate
!    yout     : Yout output of the interpolation
!    nxout    : the number of values of YOUT
!    noutmax  : the max size of table YOUT
!----------------------------------------------------------------------

subroutine itpconc(iin,yin,nxin,iout,yout,nxout,noutmax)

implicit none

integer,intent(in) :: nxin,nxout,noutmax
real(kind=iprec),INTENT(in)  :: yin(nxin),iin(nxin+1),iout(noutmax+1)
real(kind=iprec),INTENT(out) :: yout(noutmax)

real(kind=iprec) :: yin_norm(nxin),inter
integer          :: i,j
! global init

yout=0.

do i=1,nxin
   yin_norm(i)=yin(i)/(log(iin(i+1)/iin(i)))
enddo

do j=1,nxout
  do i=1,nxin
    inter=intersec(log(iout(j)),log(iout(j+1)),log(iin(i)),log(iin(i+1)))
    yout(j)=yout(j)+inter*yin_norm(i)
  enddo
enddo

end subroutine itpconc

!//////////////////////////////////////////////////////////////////////////////////////////

subroutine layave(nlevels,nlayers,alti,x,xlay,hlay) 

! Gives CHIMERE layer averages of pressure
! nlevels, alti, x : number of vertical levels, altitude, and value in the met model
! hlay: layer top altitudes in the output file
! xlay: output averaged values

! This subroutine is the trapezoidal version. It is mandatory to use this one for pressure because if one 
! assumes constant pressure in each WRF layer this may result in errors in subroutines that uses the interpolated 
! pressure to calculate vertical gradients as in calc_deepconv 

use chimere_params, only: iprec
implicit none

!--------------------------------------------------------------------------------
! subroutine arguments
integer,intent(in) :: nlevels
integer,intent(in) :: nlayers
real(kind=iprec),intent(in),dimension(0:nlevels) :: alti
real(kind=iprec),intent(in),dimension(0:nlevels) :: x
real(kind=iprec),intent(out),dimension(nlayers) :: xlay
real(kind=iprec),intent(in),dimension(nlayers) :: hlay

! local variables
real(kind=iprec) :: inter,suminter,hbot,htop
integer          :: i,j

!--------------------------------------------------------------------------------

xlay=0d0

do j=1,nlayers
   suminter=0d0
   do i=1,nlevels
      ! We are already over the target CHIMERE layer. Exit the loop on met levels.
      if(alti(i-1).ge.hlay(j))exit       
      if(j.ge.2)then
         ! We are still below the target CHIMERE layer. Cycle to the next i
         if(alti(i).le.hlay(j-1))cycle       
         hbot=max(hlay(j-1),alti(i-1))
      else
         hbot=max(0d0,alti(i-1))
      endif       
      htop=min(hlay(j),alti(i))
      inter=max(htop-hbot,0d0)
      if(inter.gt.0d0)then
         xlay(j)=xlay(j)+inter*(x(i-1)+0.5*(htop+hbot-2d0*alti(i-1))*(x(i)-x(i-1))/(alti(i)-alti(i-1)))
         suminter=suminter+inter
      endif
   enddo
   if(suminter.gt.0d0) xlay(j)=xlay(j)/suminter
enddo

!--------------------------------------------------------------------------------
end subroutine layave

!//////////////////////////////////////////////////////////////////////////////////////////

subroutine chm2wrf_aod_vert(iin,yin,nxin,iout,yout,nxout,noutmax)

implicit none

integer,intent(in) ::  nxin,nxout,noutmax
real(kind=iprec),intent(in)  :: yin(nxin),iin(nxin+1),iout(noutmax+1)
real(kind=iprec),intent(out) :: yout(noutmax)

real(kind=iprec) :: yin_norm(nxin),inter
integer          :: i,j
! global init

yout=0.

do i=1,nxin
   yin_norm(i)=yin(i)/(iin(i+1)-iin(i))
enddo

do j=1,nxout
  do i=1,nxin
    inter=intersec(iout(j),iout(j+1),iin(i),iin(i+1))
    yout(j)=yout(j)+inter*yin_norm(i)
  enddo
enddo

end subroutine chm2wrf_aod_vert

!//////////////////////////////////////////////////////////////////////////////////////////

function intersec(mina,maxa,minb,maxb)
! returns the size of intersection of [mina;maxa] with [minb;maxb]
implicit none

real(kind=iprec)  :: intersec,mina,maxa,minb,maxb

if((mina.gt.maxb).or.(minb.gt.maxa))then
  intersec=0.
else
  intersec=min(maxa,maxb)-max(mina,minb)
endif

end function intersec

!//////////////////////////////////////////////////////////////////////////////////////////

subroutine itpconcLUT(iin,nxin,iout,nxout,noutmax,itvlog,interlut)

! Same as itpconc but returns only:
! 1. the intersec table
! 2. the log(iin(i+1)/iin(i)) term

implicit none

integer :: i,j
integer,intent(in) :: nxin,nxout,noutmax
real(kind=iprec),INTENT(in)  :: iin(nxin+1),iout(noutmax+1)
real(kind=iprec),INTENT(out) :: itvlog(nxin), interlut(nxin,nxout)

do i=1,nxin
   itvlog(i)=log(iin(i+1)/iin(i))
enddo

do j=1,nxout
   do i=1,nxin
      interlut(i,j)=intersec(log(iout(j)),log(iout(j+1)),log(iin(i)),log(iin(i+1)))
   enddo
enddo

end subroutine itpconcLUT

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

subroutine itpconcAPP(yin,nxin,yout,nxout,noutmax,itvlog,interlut)

! Same as itpconc but only apply the Look Up Table to the concentrations

implicit none

integer :: i,j
integer,intent(in) :: nxin,nxout,noutmax
real(kind=iprec),INTENT(in) :: yin(nxin)
real(kind=iprec),INTENT(in) :: itvlog(nxin), interlut(nxin,noutmax)
real(kind=iprec),INTENT(out) :: yout(noutmax)
real(kind=iprec) :: yin_norm(nxin)

yout=0.

do i=1,nxin
   yin_norm(i)=yin(i)/itvlog(i)
enddo

do j=1,nxout
   do i=1,nxin
      yout(j)=yout(j)+interlut(i,j)*yin_norm(i)
   enddo
enddo

end subroutine itpconcAPP


!//////////////////////////////////////////////////////////////////////////////////////////
!-----------------------------------------------------------------------
! Cumul of concentrations
! from fine (xin) to coarse bins intervals (xout) resolution
!----------------------------------------------------------------------
!  xin: mean mass median diameter of the fine distribution
!  yin: value in the fine distribution
!  nxin: number of values of the fine distribution
!  xout: mean mass median diameter intervals of the coarse distribution
!  yout: value cumulated in the coarse bins distribution
!  nxout: the number of values of yout
!----------------------------------------------------------------------
subroutine cuminbins(xin,yin,nxin,xout,yout,nxout)

implicit none

integer nxin,nxout,i,j

real(kind=iprec) :: xin(nxin),yin(nxin)
real(kind=iprec) :: xout(nxout+1),yout(nxout)
!----------------------------------------------------------------------

yout=0.

do j=1,nxout
   do i=1,nxin
      if(xin(i).ge.xout(j).and.xin(i).lt.xout(j+1))then
  yout(j)=yout(j)+yin(i)
      endif
   enddo
enddo

end subroutine cuminbins

!//////////////////////////////////////////////////////////////////////////////////////////

subroutine distrib_lognorm(dmean,sigma,interv,nbins,dist)

implicit none
integer :: nb
integer,intent(in) :: nbins
real(kind=iprec),intent(in) :: dmean, sigma
real(kind=iprec) :: sumdist, diambin1, diambin2, xl1, xl2
real(kind=iprec), dimension(nbins),intent(out) :: dist
real(kind=iprec), dimension(nbins+1), intent(in) :: interv

do nb=1,nbins
   diambin1=interv(nb)
   diambin2=interv(nb+1)
   xl1=erf((log(diambin1/dmean))/(sqrt(2.d0)*log(sigma)))
   xl2=erf((log(diambin2/dmean))/(sqrt(2.d0)*log(sigma)))
   dist(nb)=0.5d0*abs(xl2-xl1)
enddo

sumdist=sum(dist)

do nb=1,nbins
   dist(nb)=dist(nb)/sumdist
enddo


end subroutine distrib_lognorm

!//////////////////////////////////////////////////////////////////////////////////////////

subroutine horint(xin,ideb,nlev,nhz,wx,wy,ix,iy,xout)

! Horizontal interpolation with the weights given by gridloc()
! ideb - starting vertical level for interpolation (if you want to process only
! a few upper levels)

 real(kind=iprec),dimension(:,:,:), intent(in)  :: xin
 integer, intent(in)      :: nlev, nhz, ideb
 real(kind=iprec), dimension(nhz), intent(in) :: wx,wy
 integer,dimension(nhz), intent(in) :: ix,iy
 real(kind=iprec),dimension(nhz,nlev), intent(out)  :: xout
 integer        :: il,ihz

 integer        :: size_x,size_y

 size_x=size(xin,1)
 size_y=size(xin,2)

 do il=ideb,nlev
    do ihz=1,nhz
       if((ix(ihz).lt.size_x).and.(iy(ihz).lt.size_y))then
       xout(ihz,il) =                 &
      & (1.-wx(ihz))*(1.-wy(ihz))*xin(ix(ihz)  ,iy(ihz),  il  ) &
      & +     wx(ihz) *(1.-wy(ihz))*xin(ix(ihz)+1,iy(ihz),  il  ) &
      & + (1.-wx(ihz))*  wy(ihz) *xin(ix(ihz)  ,iy(ihz)+1,il )  &
      & +     wx(ihz) *  wy(ihz) *xin(ix(ihz)+1,iy(ihz)+1,il )
       else
   xout(ihz,il)=xin(ix(ihz)  ,iy(ihz),  il  )
       endif
    enddo
 enddo

end subroutine horint

!//////////////////////////////////////////////////////////////////////////////////////////
!//////////////////////////////////////////////////////////////////////////////////////////


!=========================================================
!  function findindex(name, species, length)
!    ! returns the index of string <name> in array of strings <species>
!    ! returns -1 if not found
!
!    integer :: findindex
!    ! arguments
!    character(len=*) :: name
!    character(len=length),dimension(:) :: species
!    integer :: length
!
!    ! local variables
!    integer :: i
!
!    findindex = -1
!    do i=1,size(species,1)
!       !print *,i, name(1:len_trim(name)), ' ',  species(i)(1:len_trim(species(i)))
!       if (name(1:len_trim(name)) == species(i)(1:len_trim(species(i)) ) ) then
!          findindex = i
!          exit
!       end if
!    end do
!
!  end function findindex



!=========================================================

subroutine modisec(dist1,nsec1,dist2,nsec2,nmax,dis1todis2)
! modisec w/o logarithmic ponderation, which is adapted only for aerosols
  implicit none
  integer nsec1,nsec2,nmax
  integer i,j
  integer loca(nsec2)
  real(kind=iprec) dist1(nmax)
  real(kind=iprec),intent(out) :: dis1todis2(nmax,nsec2-1)
  real(kind=iprec) :: dist2(nsec2)

  do i=1,nsec2-1
     do j=1,nsec1-1
        dis1todis2(j,i)=0.0
     enddo
  enddo

  if(dist2(1).lt.dist1(1)) loca(1)=0
  do i=1,nsec2-1
     if(dist2(i+1).lt.dist1(1)) loca(i+1)=0
     do j=1,nsec1-1
        if(dist2(i).lt.dist1(j+1).and.dist2(i).ge.dist1(j))             &
             &   loca(i)=j
     enddo
     if(dist2(i).ge.dist1(nsec1)) loca(i)=nsec1
  enddo
  if(dist2(nsec2).gt.dist1(nsec1)) loca(nsec2)=nsec1
  if(dist2(nsec2).lt.dist1(nsec1)) loca(nsec2)=nsec1-1
  do j=1,nsec1-1
     do i=1,nsec2-1
        if(loca(i).eq.j.and.loca(i+1).eq.j) then
           dis1todis2(j,i)=real((log(dist2(i+1))-log(dist2(i)))/           &
                &                       (log(dist1(j+1))-log(dist1(j))))
        elseif(loca(i).eq.j.and.loca(i+1).gt.j) then
           dis1todis2(j,i)=real((log(dist1(j+1))-log(dist2(i)))/           &
                &                       (log(dist1(j+1))-log(dist1(j))))
        elseif(loca(i).lt.j.and.loca(i+1).eq.j) then
           dis1todis2(j,i)=real((log(dist2(i+1))-log(dist1(j)))/           &
                &                       (log(dist1(j+1))-log(dist1(j))))
        elseif(loca(i).lt.j.and.loca(i+1).gt.j) then
           dis1todis2(j,i)=1.0
        endif
     enddo
  enddo
END subroutine modisec


!-------------------------------------------------------------
subroutine interp_bilin2( xgi, ygi, xin,  xgo, ygo, xout, calc_weights, dealloc )
   ! 2-dim bilinear interpolation using new findcell and wcal
   ! Takes 2D I and O grids and a 2D var, returns 2D bilinear interpolated 2D var

  implicit none

   ! ----------------------------------------------------------------
   ! subroutine arguments
   real(kind=iprec),dimension(:,:),intent(in) :: xgi,ygi ! coords of input grid
   real(kind=iprec),dimension(:,:),intent(in) :: xgo,ygo ! coords of output grid
   real(kind=iprec),dimension(:,:),intent(in) :: xin     ! input variable
   logical,intent(in)             :: calc_weights
   logical,intent(in),optional    :: dealloc
   real(kind=iprec),dimension(:,:),intent(out):: xout    ! output variable

   ! local variables
   integer :: nxi,nyi                            ! dimensions of the input grid
   integer :: nxo,nyo                            ! dimensions of the output grid

   integer,allocatable, dimension(:), save  :: ix,iy   ! indices of the input grid cell
   !                                               ! which contains the current output point M
   real(kind=iprec), allocatable,   dimension(:), save  :: wx,wy   ! bilinear weights of M point
   real(kind=iprec), allocatable,   dimension(:)    :: xgo1d, ygo1d
   real(kind=iprec),allocatable,dimension(:,:,:,:) :: corners_3d
   integer :: i,j
   real(kind=iprec) :: x,y
   integer :: ixg,iyg,lstix,lstiy, npts, np
   real(kind=iprec) :: lxg,lyg
   logical :: found, in_is_recti
   real(kind=iprec) :: deltax,deltay

   ! ----------------------------------------------------------------
   nxi = size(xin,1)
   nyi = size(xin,2)
   nxo = size(xout,1)
   nyo = size(xout,2)
   npts = nxo*nyo

   if (calc_weights) then

      in_is_recti = test_rectilon(xgi,nxi,nyi) .and. &
          &           test_rectilat(ygi,nxi,nyi)
      if (in_is_recti) then
         print *,'**** interp_bilin2: input domain is rectilinear, fast !'
         deltax = (xgi(nxi, nyi/2) - xgi(1, nyi/2))/(nxi-1)
         deltay = (ygi(nxi/2, nyi) - ygi(nxi/2, 1))/(nyi-1)
      else
         print *,'**** interp_bilin2: input domain is NOT rectilinear.'
         ! Calculate the corners
         call calculate_corners(xgi, ygi, corners_3d, nxi, nyi)
      end if

      allocate( ix(nxo*nyo), iy(nxo*nyo), wx(nxo*nyo), wy(nxo*nyo) )
      allocate( xgo1d(nxo*nyo), ygo1d(nxo*nyo) )
      allocate(corners_3d(nxi, nyi, 4, 3))

      xgo1d = reshape(xgo, (/nxo*nyo/))
      ygo1d = reshape(ygo, (/nxo*nyo/))

      ixg = 1
      iyg = 1
      lstix = 1
      lstiy = 1
      ! Calculate grid correspondance indices and weights
      do i=1,npts
         x = xgo1d(i)
         y = ygo1d(i)
         if (in_is_recti) then
            lxg = 1 + (x - xgi(1,     nyi/2))/deltax
            lyg = 1 + (y - ygi(nxi/2, 1    ))/deltay
            found = (lxg>=1).and.(lxg<=nxi).and.(lyg>=1).and.(lyg<=nyi)
            ixg = floor(lxg)
            iyg = floor(lyg)
            if (.not. found) then
               print *,'interp_bilin2: output grid point (', x, ', ', y, ') outside input domain!'
               print *, '   input x: ', xgi(1,1), ':', xgi(nxi,nyi), ';   y: ', ygi(1,1), ':', ygi(nxi,nyi)
               call exit1('Bye!')
            else
               wx(i) = lxg-ixg
               wy(i) = lyg-iyg
               ix(i) = ixg
               iy(i) = iyg
            end if
         else
            call findcell(x,y,ixg,iyg,nxi,nyi, corners_3d, lstix,lstiy,found)
            if (found) then
              call wcal(x,y,ixg,iyg,nxi,nyi,corners_3d,wx(i),wy(i))
              ix(i) = ixg
              iy(i) = iyg
            else
              print *,x,y,xgi(1,1),xgi(nxi,nyi),ygi(1,1),ygi(nxi,nyi)
              print *,'*** interp_bilin2: output point outside input domain ! Bye'
              call exit1('error in interp_bilin2, sorry')
            end if
         end if
      end do
      deallocate( xgo1d, ygo1d )

   end if ! calc_weights

   ! Compute interpolated variable
   do j=1,nyo
      do i=1,nxo
         np=i+(j-1)*nxo
         xout(i,j) = (1.-wx(np))*(1.-wy(np))*xin(ix(np)  ,iy(np)) &
            +     wx(np) *(1.-wy(np))*xin(ix(np)+1,iy(np)) &
            + (1.-wx(np))*    wy(np) *xin(ix(np)  ,iy(np)+1) &
            +     wx(np) *    wy(np) *xin(ix(np)+1,iy(np)+1)
      end do
   end do

   if (present(dealloc)) then
      if (dealloc) deallocate( ix, iy, wx, wy )
   end if

end subroutine interp_bilin2

!-------------------------------------------------------------
subroutine lintp(p1,c1,n1,ps,p2,c2,n2)
! linear interpolation of 1D array in p-coordinates
! p is supposed to be at cell edges, interpolated quantity in cell centers
! ps - surface pressure
!
integer, intent(in)    :: n1,n2
real(kind=iprec),dimension(n1),intent(in) :: p1,c1
real(kind=iprec),intent(in)               :: ps
real(kind=iprec),dimension(n2),intent(in) :: p2
real(kind=iprec),dimension(n2),intent(out) :: c2
real(kind=iprec),dimension(n1)          :: p1c
real(kind=iprec),dimension(n2)          :: p2c
integer                     :: i

! 1. get pressures at cell centers

p1c(1) = (ps + p1(1))*0.5
do i=2,n1
   p1c(i) = (p1(i-1) + p1(i))*0.5
enddo

p2c(1) = (ps + p2(1))*0.5
do i=2,n2
   p2c(i) = (p2(i-1) + p2(i))*0.5
enddo

! 2. interpolate arrays at cell centers

call lint1(p1c,c1,n1,p2c,c2,n2)

end subroutine lintp

!-------------------------------------------------------------
subroutine lint1(x1,y1,n1,x2,y2,n2)
! linear interpolation of 1D array
! Attention : x1 and x2 must decrease monotonously, since those are pressure arrays!!
!
integer, intent(in) :: n1,n2
real(kind=iprec),dimension(n1),intent(in) :: x1,y1
real(kind=iprec),dimension(n2),intent(in) :: x2
real(kind=iprec),dimension(n2),intent(out) :: y2
integer :: i,j
real(kind=iprec) :: w

do i=1,n2

  ! Find corresponding points of the x2.
  j = 1
  do while ( .not.(x1(j) >= x2(i) .and. x2(i) > x1(j+1)) )
     if (j==n1-1) then
        j = 0; exit
     end if
     j = j + 1
  end do

  ! If out of borders -> take boundary values otherwise interpolate

  if (j > 0) then
     w=(x2(i)-x1(j))/(x1(j+1) - x1(j))
     y2(i) = w*y1(j+1)+(1.-w)*y1(j)
  else
     if ( x2(i) > x1(n1) ) then
        y2(i) = y1(1)
     elseif ( x2(i) < x1(1) ) then
        y2(i) = y1(n1)
     else
        print *, 'linint1 error: something wrong with extrapolation condition'
        print *, '  destination point number',i
        print *, '  its x value',x2(i)
        print *, 'x values:', x1(:)
        call exit1('error in linint1, sorry')

     end if
  end if

end do ! i

end subroutine lint1

!-------------------------------------------------------------
! Some basic geometric functions

subroutine normalize(v,vnorm)
use chimere_params, only: iprec
implicit none

real(kind=iprec),dimension(3) :: v,vnorm
real(kind=iprec)              :: vv
! real(kind=iprec)              :: norme

vv=norme(v)

if (vv.ne.0.)then
  vnorm = 1. / vv * v
else
  print*,'ERROR IN normalize (interp_tools.F90) : cannot normalize null vector. Bye !'
  call exit1('error in normalize, sorry')
endif
 
end subroutine normalize


!-------------------------------------------------------------
function norme(x)
   ! Computes module of 3d vector x
   use chimere_params, only: iprec
   implicit none
   
!   real(kind=iprec) scalar Declaration is not permitted because in the same module
   real(kind=iprec),dimension(3)  :: x
   real(kind=iprec) norme
        
   norme=sqrt(scalar(x,x))

end function norme

!-------------------------------------------------------------
function scalar(x,y)
   ! Computes dot product of two 3d vectors x and y
   use chimere_params, only: iprec
   implicit none
   
   real(kind=iprec) scalar
   real(kind=iprec),dimension(3)  :: x,y
   
   scalar = x(1)*y(1) + x(2)*y(2) + x(3) * y(3)

end function scalar

!-------------------------------------------------------------
function normevec(x,y)
  ! Computes the norm of cross product of two 3d vectors x and y
implicit none

real normevec
real,dimension(3)  :: x,y
real,dimension(3)  :: vec

vec(1)=x(2)*y(3)-x(3)*y(2)
vec(2)=x(3)*y(1)-x(1)*y(3)
vec(3)=x(1)*y(2)-x(2)*y(1)

normevec=norme(vec)

end function normevec
!-------------------------------------------------------------
function surface_from3dcorners(x3d)
! So far this routine is used in emiSURF only

! This routine calculates the area in km2 of a square defined by its x-y-z coordinates (on Earth)

! INPUT
real, INTENT(IN), dimension(4,3) :: x3d

! OUTPUT
real				 :: surface_from3dcorners

! LOCAL 	  
real,parameter  		 :: radius=6371000.,degrad=.0174533
real,dimension(3)		 :: xA, xB, xC, xD

xA=x3d(1,:)
xB=x3d(2,:)
xC=x3d(3,:)
xD=x3d(4,:)
surface_from3dcorners=(0.5 * normevec(xB-xA,xD-xA) &
  + 0.5 * normevec(xD-xC,xB-xC))*radius**2
surface_from3dcorners=surface_from3dcorners/1.e6 ! in km2

end function surface_from3dcorners

!==================================================================================

subroutine make_corners_area(lon,lat,nx,ny,corners_3d,is_periodic,xe,xw,ys,yn,domain_is_recti,area)

! So far this routine is used in emiSURF only


! INPUT variables
  integer,INTENT(IN)	  :: nx,ny
  real,INTENT(IN)	  :: lat(nx,ny),lon(nx,ny)

! OUTPUT variables
  real,INTENT(OUT)	  :: corners_3d(nx,ny,4,3)
  real,INTENT(OUT)	  :: area(nx,ny)
  logical,INTENT(OUT)	  :: domain_is_recti,is_periodic
  real,INTENT(OUT)	  :: xe,xw,ys,yn

! Functions
! No function declaration permitted here since we are in the same module (?!?)
!    real                    :: surface_from3dcorners            
!    logical                 :: test_rectilon,test_rectilat

! local variables
integer                 :: ix,iy

! Calculate grid corners for findcell
is_periodic=.false.
call calculate_corners(lon,lat,corners_3d,nx,ny,.false.,is_periodic)
call calculate_bound(nx,ny,corners_3d,xw,xe,ys,yn)

is_periodic=((xe-xw).gt.355.)


! We have detected that the input domain has global coerage in longitude
! Recalculate the corners accordingly

if(is_periodic)then
   call calculate_corners(lon,lat,corners_3d,nx,ny,.false.,is_periodic)
   call calculate_bound(nx,ny,corners_3d,xw,xe,ys,yn)	
endif

! check if grid is rectilinear
domain_is_recti= test_rectilon(lon,nx,ny) .and.test_rectilat(lat,nx,ny)

! Derive grid cell area
do iy=1,ny
do ix=1,nx
  area(ix,iy) = surface_from3dcorners(corners_3d(ix,iy,:,:))
enddo
enddo

end subroutine make_corners_area

!==================================================================================
end module interp_tools