!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
! physics_surface_coeff.f90 
!
! History
!   2005/09/21 Yamada Yukiko     create
!

module physics_surface_coeff_mod

  implicit none

  private
  public :: physics_surface_coeff

contains

  subroutine physics_surface_coeff( &
       & xy_SurfVelBulkCoeff , & ! (out) 
       & xy_SurfTempBulkCoeff, & ! (out) 
       & xy_SurfQvapBulkCoeff, & ! (out) 
       & xy_SurfBulkRiNum    , & ! (in) 
       & xy_SurfVelAbs       , & ! (in) 
       & xy_SurfVelRoughLength   , & ! (in) 
       & xy_SurfTempRoughLength  , & ! (in) 
       & xy_SurfGeoPot            ) ! (in) 

    !==== Dependency
    use type_mod,    only: REKIND, DBKIND, INTKIND, TOKEN, STRING
    use grid_3d_mod, only: im, jm, km
    use constants_mod, only: FKarm 
    use dc_trace,    only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump

    implicit none

    !==== Output
    !
    real(DBKIND), intent(out) :: &
         & xy_SurfVelBulkCoeff(im,jm)      , & ! (out) バルク係数：運動量
         & xy_SurfTempBulkCoeff(im,jm)     , & ! (out) バルク係数：温度
         & xy_SurfQvapBulkCoeff(im,jm)         ! (out) バルク係数：比湿

    !==== Input
    !
    real(DBKIND), intent(in) :: &
         & xy_SurfBulkRiNum(im,jm)         , & ! (in) バルクＲｉ数
         & xy_SurfVelAbs(im,jm)            , & ! (in)  風速絶対値
         & xy_SurfVelRoughLength(im,jm)    , & ! (in) 地表粗度長 (運動量)
         & xy_SurfTempRoughLength(im,jm)   , & ! (in) 地表粗度長 (熱)
         & xy_SurfGeoPot(im,jm)                ! (in) 最下層温度

    !----- 作業用内部変数 -----
    character(STRING),  parameter:: subname = "physics_surface_coeff"

    ! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用)
    integer(INTKIND)    :: i, j

    logical            :: Newtral = .FALSE.                 ! 中立であるか否か
    real(DBKIND), parameter :: ConstBulkCoeff  = -1.        ! バルク定数一定値
    real(DBKIND), parameter :: VelBulkCoeffMin  = 0.        ! ｕバルク係数最小値
    real(DBKIND), parameter :: TempBulkCoeffMin = 0.        ! Ｔバルク係数最小値
    real(DBKIND), parameter :: QvapBulkCoeffMin = 0.        ! ｑバルク係数最小値
    real(DBKIND), parameter :: VelBulkCoeffMax  = 1.        ! ｕバルク係数最大値
    real(DBKIND), parameter :: TempBulkCoeffMax = 1.        ! Ｔバルク係数最大値
    real(DBKIND), parameter :: QvapBulkCoeffMax = 1.        ! ｑバルク係数最大値
    
    continue

    !----------------------------------------------------------------
    !   開始処理
    !----------------------------------------------------------------
    call BeginSub(subname)

    !----------------------------------------------------------------
    !   地表面バルク係数計算
    !----------------------------------------------------------------

    ! ---- 1. 中立バルク係数 ----

      if ( ConstBulkCoeff .LT. 0.   ) then

         xy_SurfVelBulkCoeff  = ( FKarm &
              &                  / LOG ( xy_SurfGeoPot / xy_SurfVelRoughLength ) )**2
         xy_SurfTempBulkCoeff = ( FKarm &
              &                  / LOG ( xy_SurfGeoPot / xy_SurfTempRoughLength ) )**2
         xy_SurfQvapBulkCoeff = xy_SurfTempBulkCoeff

      else
         xy_SurfVelBulkCoeff  = ConstBulkCoeff
         xy_SurfTempBulkCoeff = ConstBulkCoeff
         xy_SurfQvapBulkCoeff = ConstBulkCoeff
      end if

    ! ---- 2. 非中立バルク係数 ----

      if ( .NOT. Newtral ) then

         do i = 1, im
            do j = 1, jm
               if ( xy_SurfBulkRiNum(i,j) .GT. 0. ) then 
                  xy_SurfVelBulkCoeff(i,j) = xy_SurfVelBulkCoeff(i,j) &
                       &     / ( 1. + 10.* xy_SurfBulkRiNum(i,j) &
                       &              / SQRT( 1. + 5.* xy_SurfBulkRiNum(i,j) ) )
                  xy_SurfTempBulkCoeff(i,j) = xy_SurfTempBulkCoeff(i,j) &
                       &     / ( 1. + 15.* xy_SurfBulkRiNum(i,j) &
                       &              / SQRT( 1. + 5.* xy_SurfBulkRiNum(i,j) ) )
                  xy_SurfQvapBulkCoeff(i,j) = xy_SurfTempBulkCoeff(i,j)
               else
                  xy_SurfVelBulkCoeff(i,j) = xy_SurfVelBulkCoeff(i,j) &
                       &     * ( 1. - 10.* xy_SurfBulkRiNum(i,j) &
                       &              / ( 1. + 75. * xy_SurfVelBulkCoeff(i,j) &
                       &                      * SQRT( - xy_SurfGeoPot(i,j) &
                       &                               / xy_SurfVelRoughLength(i,j) &
                       &                               * xy_SurfBulkRiNum(i,j) ) &
                       &                ) &
                       &      )

                  xy_SurfTempBulkCoeff(i,j) = xy_SurfTempBulkCoeff(i,j) &
                       &     * ( 1. - 15.* xy_SurfBulkRiNum(i,j) &
                       &              / ( 1. + 75. * xy_SurfTempBulkCoeff(i,j) &
                       &                      * SQRT( - xy_SurfGeoPot(i,j) &
                       &                               / xy_SurfTempRoughLength(i,j) &
                       &                               * xy_SurfBulkRiNum(i,j) ) &
                       &                ) &
                       &      )
                  xy_SurfQvapBulkCoeff(i,j) = xy_SurfTempBulkCoeff(i,j)
               end if               
            end do
         end do

      end if

    ! ---- 3. 最大/最小 判定 ----

      do i = 1, im
         do j = 1, jm
         xy_SurfVelBulkCoeff(i,j)  = &
              & MAX( MIN( xy_SurfVelBulkCoeff(i,j) , VelBulkCoeffMax ), & 
              &      VelBulkCoeffMin )
         xy_SurfTempBulkCoeff(i,j) = &
              & MAX( MIN( xy_SurfTempBulkCoeff(i,j) , TempBulkCoeffMax ), & 
              &      TempBulkCoeffMin )
         xy_SurfQvapBulkCoeff(i,j) = &
              & MAX( MIN( xy_SurfQvapBulkCoeff(i,j) , QvapBulkCoeffMax ), & 
              &      QvapBulkCoeffMin )
         end do
      end do

    !----------------------------------------------------------------
    !   終了処理
    !----------------------------------------------------------------
    call EndSub(subname)

  end subroutine physics_surface_coeff


end module physics_surface_coeff_mod