!-----------------------------------------------------------------------
!
      MODULE module_INTERP2D
!
!-----------------------------------------------------------------------
!
      USE module_KINDS
      USE MPI_MORE,ONLY: GLOB_ABORT
!
!-----------------------------------------------------------------------

      IMPLICIT NONE

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

      CONTAINS

!-----------------------------------------------------------------------
!
      SUBROUTINE INTERP_2D(ARRAY_IN,MASK_IN,IM_IN,JM_IN                 & 
                          ,I_SW,J_SW,SPACE_RATIO                        &
                          ,ARRAY_OUT,MASK_OUT,IM_OUT,JM_OUT)
!
!-----------------------------------------------------------------------
!***  Interpolate a 2-D H-pt array spanning the input (parent) B-grid
!***  to the output (child) B-grid.
!-----------------------------------------------------------------------
      IMPLICIT NONE
!-----------------------------------------------------------------------
!
!------------------------
!***  Argument variables
!------------------------
!
      INTEGER(kind=KINT),INTENT(IN) :: IM_IN ,JM_IN                     &  !<-- I,J extent of the full domain input array
                                      ,IM_OUT,JM_OUT                    &  !<-- I,J extent of the full domain output array 
                                      ,I_SW,J_SW                        &  !<-- Parent I,J of SW corner of child grid
                                      ,SPACE_RATIO                         !<-- Ratio of parent grid increment to child's
!
      REAL(kind=KFPT),DIMENSION(1:IM_IN,1:JM_IN),INTENT(IN):: ARRAY_IN     !<-- Full input (parent) array to be interpolated to child
      REAL(kind=KFPT),DIMENSION(1:IM_IN,1:JM_IN),INTENT(IN):: MASK_IN      !<-- 1.0 -> input array valid; 0.0 -> input array invalid
      REAL(kind=KFPT),DIMENSION(1:IM_OUT,1:JM_OUT),INTENT(IN):: MASK_OUT   !<-- 1.0 -> output array valid; 0.0 -> output array invalid
!
      REAL(kind=KFPT),DIMENSION(1:IM_OUT,1:JM_OUT),INTENT(OUT):: ARRAY_OUT !<-- Full output (child) array 
!
!---------------------
!***  Local variables
!---------------------
!
      integer I_NE,J_NE,i,j,k,ic,jc,IM1,JM1
      real(kind=KFPT) z,a1,a2,x0,x1,x2,x3,y0,y1,y2,y3
      real(kind=KFPT),parameter:: too_low=-1.E20,too_low10=-1.E30
      real,dimension(:,:),allocatable:: f_out
      integer,dimension(:,:),allocatable:: mask_utr
      real(kind=KFPT),dimension(:,:),allocatable:: f_in
!
!-----------------------------------------------------------------------
!***********************************************************************
!-----------------------------------------------------------------------
!
      if (SPACE_RATIO.ne.3) then
        print*,'SPACE_RATIO other than 3 not allowed',SPACE_RATIO
        call GLOB_ABORT(1,'SPACE RATIO.ne.3',1)
      end if

      i=(IM_OUT+1)/SPACE_RATIO
      j=(JM_OUT+1)/SPACE_RATIO
      I_NE=I_SW+i
      J_NE=J_SW+j
      if (I_NE.gt.IM_in .or. J_NE.gt.JM_in) then
        print*,'AM: INTERP_2D: child not inside parent;'// &
        'IM_in,JM_in,IM_out,JM_out,I_SW,J_SW,SPACE_RATIO: ', &
        IM_in,JM_in,IM_out,JM_out,I_SW,J_SW,SPACE_RATIO
        call GLOB_ABORT(1,'AM: INTERP_2D: child not inside parent',1)
      end if
      IM1=i*SPACE_RATIO+1
      JM1=j*SPACE_RATIO+1
      allocate(f_out(IM1,JM1))
      allocate(mask_utr(I_SW:I_NE,J_SW:J_NE))
      allocate(f_in(I_SW:I_NE,J_SW:J_NE))

      f_out=too_low10
      f_in=ARRAY_in(I_SW:I_NE,J_SW:J_NE)
      mask_utr=MASK_IN(I_SW:I_NE,J_SW:J_NE)

      a1=1./3.
      a2=2./3.

!loop
      do j=J_SW,J_NE-1
      do i=I_SW,I_NE-1
!loop
        k=mask_utr(i,j)+mask_utr(i+1,j)+ &
          mask_utr(i+1,j+1)+mask_utr(i,j+1)
!
        if (k.eq.0) CYCLE
!
        if (k.lt.4) then
          z=MASK_IN(i,j)*ARRAY_IN(i,j)+ &
            MASK_IN(i+1,j)*ARRAY_IN(i+1,j)+ &
            MASK_IN(i+1,j+1)*ARRAY_IN(i+1,j+1)+ &
            MASK_IN(i,j+1)*ARRAY_IN(i,j+1)
          z=z/k
          if (mask_utr(i,j).eq.0) f_in(i,j)=z
          if (mask_utr(i+1,j).eq.0) f_in(i+1,j)=z
          if (mask_utr(i+1,j+1).eq.0) f_in(i+1,j+1)=z
          if (mask_utr(i,j+1).eq.0) f_in(i,j+1)=z
        end if

        x0=f_in(i,j)
        x1=a2*f_in(i,j)+a1*f_in(i+1,j)
        x2=a1*f_in(i,j)+a2*f_in(i+1,j)
        x3=f_in(i+1,j)
        y0=f_in(i,j+1)
        y1=a2*f_in(i,j+1)+a1*f_in(i+1,j+1)
        y2=a1*f_in(i,j+1)+a2*f_in(i+1,j+1)
        y3=f_in(i+1,j+1)

        ic=3*(i-I_SW)+1
        jc=3*(j-J_SW)+1
        f_out(ic,jc)=x0
        f_out(ic+1,jc)=x1
        f_out(ic+2,jc)=x2
        f_out(ic+3,jc)=x3
        f_out(ic,jc+1)=a2*x0+a1*y0
        f_out(ic+1,jc+1)=a2*x1+a1*y1
        f_out(ic+2,jc+1)=a2*x2+a1*y2
        f_out(ic+3,jc+1)=a2*x3+a1*y3
        f_out(ic,jc+2)=a1*x0+a2*y0
        f_out(ic+1,jc+2)=a1*x1+a2*y1
        f_out(ic+2,jc+2)=a1*x2+a2*y2
        f_out(ic+3,jc+2)=a1*x3+a2*y3
        f_out(ic,jc+3)=y0
        f_out(ic+1,jc+3)=y1
        f_out(ic+2,jc+3)=y2
        f_out(ic+3,jc+3)=y3
!loop
      end do
      end do
!loop

      where (abs(MASK_OUT-1.).lt.0.1 .and. &
             f_out(1:IM_OUT,1:JM_OUT).gt.too_low)
        ARRAY_OUT=f_out(1:IM_OUT,1:JM_OUT)
      end where

      deallocate(f_out)
      deallocate(mask_utr)
      deallocate(f_in)

!-----------------------------------------------------------------------
!
      END SUBROUTINE INTERP_2D
!
!-----------------------------------------------------------------------
!
      END MODULE module_INTERP2D
!
!-----------------------------------------------------------------------