! First include the set of model-wide compiler flags 
#include "tm5.inc"

module MGridTools

  ! Tool subroutines related to the horizontal and vertical grid

  implicit none

  private

  !  subroutine GetTMCellIndex( xlon,xlat,ix,iy,dx,dy )
  !  subroutine GetTMCellIndex4(xlon,xlat,ix,ix2,iy,iy2,gdx,gdy)

  public :: getTMCellIndex, getTMCellIndex4

contains



  subroutine GetTMCellIndex( im, jm, xlon, xlat, ix, iy, dx, dy )
    !------------------------------------------------------------------------
    ! Determine TM5 cell indices (ix,iy) and the relative position inside 
    ! the cell (dx,dy) from the cell centre of the location (xlon,xlat).
    ! 
    ! im, jm    : TM array dimensions [in]
    ! xlon,xlat : longitude, latitude in degrees [in]
    ! ix,iy     : indices of the TM5 grid cell containing (xlon,xlat) [out]
    ! dx,dy     : relative position with respect to the cell centre, 
    !             el. (-0.5,0.5)  [out] 
    !------------------------------------------------------------------------
    use physconstants,   only : pi

    implicit none

    integer,intent(in)    :: im, jm
    real,intent(in)       :: xlon, xlat
    integer,intent(out)   :: ix, iy
    real,intent(out)      :: dx, dy
    !
    real :: dlon,dlat,dlondeg,dlatdeg,phim1,tetam1
    real :: rx,ry,todegree,firstcellmid,xlon2
    !
    todegree=180.0/pi
    !
    xlon2=xlon
    if( xlon < -180.0 ) xlon2=xlon+360.0
    if( xlon >= 180.0 ) xlon2=xlon-360.0          ! xlon2 in [-180,180)
    !
    if ( im == 0 ) then
       print*,'ERROR in ass_gridtools.f90; im = 0 '
       stop
    end if
    !
    dlon = 2.0*pi/real(im)
    dlondeg = todegree*dlon                   ! deg
    phim1 = 0.5*dlon - pi                     ! lon first grid cell
    firstcellmid = todegree*phim1             ! deg
    rx = (xlon2-firstcellmid)/dlondeg + 1.0   ! old grid: [+1.0,im+1.0)
                                              ! new grid: [+0.5,im+0.5)
    !old        dlon=360.0/real(im)
    !old        rx = xlon/dlon+0.5*im+1.0
    ix = nint(rx)      !  old: [1,im+1]    new: [1,im]
    dx = rx - real(ix)
    if ( ix == im+1 ) ix = 1
    if ( ix == 0    ) ix = im     ! for safety ...
    !
    if ( jm == 0 ) then
       print*,'ERROR in ass_gridtools.f90; jm = 0 '
       stop
    end if
    !
    dlat = pi/real(jm)
    dlatdeg = todegree*dlat                      	! deg
    tetam1 = 0.5*dlat-0.5*pi                   		! lat first grid cell
    firstcellmid = todegree*tetam1               	! deg
    ry = (xlat-firstcellmid)/dlatdeg+1.0       		! [0.5,jm+0.5]
    iy = nint(ry)         				! [0,jm+1] for rounding error
    if ( iy < 1  ) iy = 1     				! [1,jm+1]
    if ( iy > jm ) iy = jm                              ! [1,jm]
    !old    ry = 0.5+(xlat+90.0)*real(jm)/180.0
    !old    iy = 1+int((xlat+90.0)*real(jm)/180.0)
    dy = ry - real(iy)
    !
  end subroutine GetTMCellIndex



  subroutine GetTMCellIndex4 ( im, jm, xlon, xlat, ix4a, iy4a, w4a )
    !------------------------------------------------------------------------
    !
    ! Find the four grid points that
    ! are nearest to the position (lon,lat) on the globe.
    !
    ! im, jm        : TM array dimensions [in]
    ! (xlon,xlat)   : position on the globe, in degrees (in)
    ! ix4a,iy4a     : east-north indices of the four TM5 grid points (out)
    ! w4a           : weight factors of the four points (out)
    !                            Henk Eskes, KNMI
    !------------------------------------------------------------------------
    use physconstants,   only : pi

    implicit none

    ! in/out
    integer,intent(in)               :: im, jm
    real,intent(in)                  :: xlon, xlat
    integer,dimension(4),intent(out) :: ix4a, iy4a
    real,dimension(4),intent(out)    :: w4a

    ! local
    integer :: ix,ix2,iy,iy2
    real    :: gdx,gdy
    real    :: dlon,dlat,xlon2,phim1,tetam1
    real    :: rx,ry,todegree,firstcellmid,dlondeg,dlatdeg
    !
    ! determine south-west model cell indices (ix,iy) 
    ! and inside-cell cordinate (gdx,gdy), both el. [0,1] 
    ! corresponding to the pixel coordinates (longitude=xx,latitude=yy)
    ! ( lon el. [-180,180] )
    !
    todegree=180.0/pi
    !
    xlon2=xlon
    if( xlon < -180.0 ) xlon2=xlon+360.0
    if( xlon >= 180.0 ) xlon2=xlon-360.0          ! xlon2 in [-180,180)
    !
    dlon = 2.0*pi/real(im)
    dlondeg = todegree*dlon                     ! deg
    phim1 = 0.5*dlon-pi                         ! lon first grid cell, in rad
    firstcellmid = todegree*phim1               ! deg
    rx = (xlon2-firstcellmid)/dlondeg + 1.0     ! old grid: [+1.0,im+1.0)
    !                                           new grid: [+0.5,im+0.5)
    ix = int(rx)
    gdx = rx - real(ix)
    if ( ix == im+1 ) ix = 1
    if ( ix == 0    ) ix = im      ! for safety ...
    !
    dlat = pi/real(jm)
    dlatdeg = todegree*dlat                      ! deg
    tetam1 = 0.5*dlat - 0.5*pi                   ! lat first grid cell
    firstcellmid = todegree*tetam1               ! deg
    ry = (xlat-firstcellmid)/dlatdeg + 1.0       ! [0.5,jm+0.5]
    iy = int(ry)
    gdy = ry - real(iy)
    !
    ! determine coordinates of other corner points (ix,iy2),(ix2,iy),(ix2,iy2)
    ix2 = ix + 1
    if( ix2 == im + 1 ) then ! periodic east-west boundary
       ix2 = 1
    end if
    iy2 = iy + 1
    if( iy >= jm ) then ! north pole				
       iy = jm - 1							
       iy2 = jm
       gdy = 1.0									
    end if									
    if( iy <= 0 ) then ! south pole 						
       iy = 1									
       iy2 = 2									
       gdy = 0.0	 											
    end if									
    !
    ix4a(:) = (/ ix, ix2,  ix, ix2 /)
    iy4a(:) = (/ iy,  iy, iy2, iy2 /)
    !if ( (iy == 1) .or. (iy == jm) ) then  ! polar cap exception		!JDMP
    !   ix4a(1) = 1									!JDMP
    !   ix4a(2) = 1									!JDMP
    !end if										!JDMP
    !if ( (iy2 == 1) .or. (iy2 == jm) ) then   ! polar cap exception	!JDMP
    !   ix4a(3) = 1									!JDMP
    !   ix4a(4) = 1									!JDMP
    !end if										!JDMP
    w4a(1) = (1.0-gdx)*(1.0-gdy)       
    w4a(2) =      gdx *(1.0-gdy)       
    w4a(3) = (1.0-gdx)*     gdy       
    w4a(4) =      gdx *     gdy       
    !
  end subroutine GetTMCellIndex4


end module MGridTools