MODULE module_interp

  CONTAINS
!--------------------------------------------------------

  SUBROUTINE interp( data_in, nx, ny, nz, &
                     data_out, nxout, nyout, nzout, &
                     z_data, z_levs, number_of_zlevs )

  USE module_model_basics

  implicit none

 ! Arguments
  integer                                                       :: nx, ny, nz, number_of_zlevs
  real, dimension(west_east_dim,south_north_dim,bottom_top_dim) :: z_data
  real, dimension(nx,ny,nz)                                     :: data_in
  real, pointer, dimension(:,:,:)                               :: data_out
  real, dimension(number_of_zlevs)                              :: z_levs

  ! Local variables
  integer                                                       :: nxout, nyout, nzout
  real, allocatable, dimension(:,:,:)                           :: SCR2
  real, dimension(bottom_top_dim)                               :: data_in_1d, z_data_1d
  real, dimension(number_of_zlevs)                              :: data_out_1d
  integer                                                       :: i,j,k


    IF ( ALLOCATED(SCR2) ) DEALLOCATE(SCR2)
    IF ( ASSOCIATED(data_out) ) DEALLOCATE(data_out)
    nxout = nx
    nyout = ny
    nzout = nz

    !! We may be dealing with a staggered field
    IF ( nx .gt. west_east_dim ) THEN
      ALLOCATE(SCR2(west_east_dim,south_north_dim,bottom_top_dim))
      SCR2 = 0.5*(data_in(1:west_east_dim,:,:)+data_in(2:west_east_dim+1,:,:))
      nxout = west_east_dim
    ELSE IF ( ny .gt. south_north_dim ) THEN
      ALLOCATE(SCR2(west_east_dim,south_north_dim,bottom_top_dim))
      SCR2 = 0.5*(data_in(:,1:south_north_dim,:)+data_in(:,2:south_north_dim+1,:))
      nyout = south_north_dim
    ELSE IF ( nz .gt. bottom_top_dim ) THEN  
      ALLOCATE(SCR2(west_east_dim,south_north_dim,bottom_top_dim))
      SCR2 = 0.5*(data_in(:,:,1:bottom_top_dim)+data_in(:,:,2:bottom_top_dim+1))
      nzout = bottom_top_dim
    ELSE
      ALLOCATE(SCR2(nx,ny,nz))
      SCR2 = data_in
    ENDIF


    IF ( iprogram .ge. 6 .AND. nzout .gt. 10 .AND.                  &
         (vertical_type == 'p' .or. vertical_type == 'z') ) THEN

      ALLOCATE(data_out(west_east_dim,south_north_dim,number_of_zlevs))
      DO i=1,west_east_dim
      DO j=1,south_north_dim
  
        DO k=1,bottom_top_dim
          data_in_1d(k) = SCR2(i,j,k)
          z_data_1d(k)  = z_data(i,j,k)
        ENDDO
  
        CALL interp_1d( data_in_1d, z_data_1d, bottom_top_dim, &
                        data_out_1d, z_levs, number_of_zlevs,  &
                        vertical_type)
  
        DO k=1,number_of_zlevs
          data_out(i,j,k) = data_out_1d(k)
        ENDDO

      ENDDO
      ENDDO
      
      nzout = number_of_zlevs

    ELSE

      ALLOCATE(data_out(nxout,nyout,nzout))
      data_out = SCR2
      DEALLOCATE(SCR2)

    ENDIF

  END SUBROUTINE interp

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

  SUBROUTINE interp_1d( a, xa, na, b, xb, nb, vertical_type)

  implicit none

 ! Arguments
  integer, intent(in)              :: na, nb
  real, intent(in), dimension(na)  :: a, xa
  real, intent(in), dimension(nb)  :: xb
  real, intent(out), dimension(nb) :: b
  character (len=1)                :: vertical_type

  ! Local variables
  real                             :: missing_value
  integer                          :: n_in, n_out
  real                             :: w1, w2
  logical                          :: interp

!  parameter (MISSING_VALUE=1.0E37)
!!!!
!!!!AYM AYM
!!!!
  parameter (MISSING_VALUE=-9999.)



  IF ( vertical_type == 'p' ) THEN

    DO n_out = 1, nb

      b(n_out) = missing_value
      interp = .false.
      n_in = 1

      DO WHILE ( (.not.interp) .and. (n_in < na) )
        IF( (xa(n_in)   >= xb(n_out)) .and. &
            (xa(n_in+1) <= xb(n_out))        ) THEN
          interp = .true.
          w1 = (xa(n_in+1)-xb(n_out))/(xa(n_in+1)-xa(n_in))
          w2 = 1. - w1
          b(n_out) = w1*a(n_in) + w2*a(n_in+1)
        END IF
        n_in = n_in +1
      ENDDO

    ENDDO

  ELSE

    DO n_out = 1, nb
  
      b(n_out) = missing_value
      interp = .false.
      n_in = 1
  
      DO WHILE ( (.not.interp) .and. (n_in < na) )
        IF( (xa(n_in)   <= xb(n_out)) .and. &
            (xa(n_in+1) >= xb(n_out))        ) THEN
          interp = .true.
          w1 = (xa(n_in+1)-xb(n_out))/(xa(n_in+1)-xa(n_in))
          w2 = 1. - w1
          b(n_out) = w1*a(n_in) + w2*a(n_in+1)
        END IF
        n_in = n_in +1
      ENDDO

    ENDDO

  END IF

  END SUBROUTINE interp_1d

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

END MODULE module_interp