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

    subroutine hs94forcing( day, kappa, im, jm, km, sinlat, coslat, &
         gup, gvp, gtp, gphp, gpp, gduphy, gdvphy, gdtphy )

      use vtype

      implicit none

      real(dp), intent(in) :: day, kappa
      integer(i4b), intent(in) :: im, jm, km
      real(dp), intent(in) :: sinlat( jm ), coslat( jm )
      real(dp), intent(in) :: gup ( im, jm, km   ), gvp( im, jm, km ), &
           gtp( im, jm, km )
      real(dp), intent(in) :: gphp( im, jm, km+1 ), gpp( im, jm, km )
      real(dp), intent(out) :: gduphy( im, jm, km ), gdvphy( im, jm, km ), &
           gdtphy( im, jm, km )


      ! Local variables
      !
      real(dp) :: gps( im, jm )

      integer(i4b) :: i, j, k


      gps( :, : ) = gphp( :, :, km+1 )

      call hs94tforcing( day, kappa, im, jm, km, sinlat, coslat, gpp, gps, &
           gtp, gdtphy )

      call hs94uvforcing( day, im, jm, km, gpp, gps, gup, gvp, gduphy, gdvphy )

    end subroutine hs94forcing

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

    subroutine hs94tforcing( day, kappa, im, jm, km, sinlat, coslat, &
         gpp, gps, gt, gdtphy )

      use vtype

      implicit none

      real(dp), intent(in) :: day, kappa
      integer(i4b), intent(in) :: im, jm, km
      real(dp), intent(in) :: sinlat( jm ), coslat( jm )
      real(dp), intent(in) :: gpp( im, jm, km ), gps( im, jm )
      real(dp), intent(in) :: gt( im, jm, km )
      real(dp), intent(out) :: gdtphy( im, jm, km )


      ! Local variables
      !
      real(dp) :: gteq( im, jm, km )
      real(dp) :: p0, kt, ka, ks, sig, sigb, dty, dthetaz

      integer(i4b) :: i, j, k


      p0 = 1000.0d0 * 1.0d2

      sigb    = 0.7d0
      ka      = 1.0d0 / ( 40.0d0 * day )
      ks      = 1.0d0 / ( 4.0d0 * day )
      dty     = 60.0d0
      dthetaz = 10.0d0

      do k = 1, km
         do j = 1, jm
            do i = 1, im
               gteq( i, j, k ) = max( 200.0d0, &
                    ( 315.0d0 - dty * sinlat( j ) * sinlat( j ) &
                    - dthetaz * log( gpp( i, j, k ) / p0 ) &
                    * coslat( j ) * coslat( j ) ) &
                    *( gpp( i, j, k ) / p0 )**kappa )
            enddo
         enddo
      enddo

      do k = 1, km
         do j = 1, jm
            do i = 1, im
               sig = gpp( i, j, k ) / gps( i, j )
               kt = ka &
                    + ( ks - ka ) &
                    * max( 0.0d0, &
                    (sig-sigb)/(1.0d0-sigb) ) &
                    *coslat(j)*coslat(j)*coslat(j)*coslat(j)

               gdtphy( i, j, k ) = -kt * ( gt( i, j, k ) - gteq( i, j, k ) )
            enddo
         enddo
      enddo


    end subroutine hs94tforcing

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

    subroutine hs94uvforcing( day, im, jm, km, gpp, gps, gu, gv, &
         gduphy, gdvphy )

      use vtype

      implicit none

      real(dp), intent(in) :: day
      integer(i4b), intent(in) :: im, jm, km
      real(dp), intent(in) :: gpp( im, jm, km ), gps( im, jm )
      real(dp), intent(in) :: gu( im, jm, km ), gv( im, jm, km )
      real(dp), intent(out) :: gduphy( im, jm, km ), gdvphy( im, jm, km )


      ! Local variables
      !
      real(dp) :: p0, kv, kf, sig, sigb

      integer(i4b) :: i, j, k


      p0 = 1000.0d0 * 1.0d2

      sigb = 0.7d0
      kf   = 1.0d0 / day

      do k = 1, km
         do j = 1, jm
            do i = 1, im
               sig = gpp( i, j, k ) / gps( i, j )
               kv = kf * max( 0.0d0, (sig-sigb)/(1.0d0-sigb) )

               gduphy( i, j, k ) = - kv * gu( i, j, k )
               gdvphy( i, j, k ) = - kv * gv( i, j, k )
            enddo
         enddo
      enddo


    end subroutine hs94uvforcing

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