    !**************************************************************************

    subroutine calcgp( im, jm, kup, km, pint, pt, sigma, gps, gp )

      use vtype

      implicit none

      integer(i4b), intent(in ) :: im, jm, kup, km
      real(dp)    , intent(in ) :: pint, pt, sigma( km+1 ), gps( im, jm )
      real(dp)    , intent(out) :: gp( im, jm, km )


      !
      ! local variables
      !
      real(dp)     :: gpi( im, jm, km+1 ), gph( im, jm, km+1 )
      real(dp)     :: gph1d( km+1 ), gp1d( km )
      integer(i4b) :: i, j, k


      do k = 1, kup
         do j = 1, jm
            do i = 1, im
               gpi( i, j, k ) = pint - pt
            enddo
         enddo
      enddo
      do k = kup+1, km+1
         do j = 1, jm
            do i = 1, im
               gpi( i, j, k ) = gps( i, j ) - pint
            enddo
         enddo
      enddo

      do k = 1, kup
         gph1d( k ) = gpi( 1, 1, k ) * sigma( k ) + pint
      enddo
      do k = 1, kup
         do j = 1, jm
            do i = 1, im
               gph( i, j, k ) = gph1d( k )
            enddo
         enddo
      enddo
      do k = kup+1, km+1
         do j = 1, jm
            do i = 1, im
               gph( i, j, k ) = gpi( i, j, k ) * sigma( k ) + pint
            enddo
         enddo
      enddo

      do k = 1, kup
         gp1d( k ) = sqrt( gph( 1, 1, k ) * gph( 1, 1, k+1 ) )
      enddo
      do k = 1, kup
         do j = 1, jm
            do i = 1, im
               gp( i, j, k ) = gp1d( k )
            enddo
         enddo
      enddo
      do k = kup+1, km
         do j = 1, jm
            do i = 1, im
               gp( i, j, k ) = sqrt( gph( i, j, k ) * gph( i, j, k+1 ) )
            enddo
         enddo
      enddo


    end subroutine calcgp

    !**************************************************************************
