!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
! physics_ground.f90 
! SST, 及び地表面諸量の読み込み
!
! History
!   2005/09/20 Yamada Yukiko     create
!   2007-05-17 M. Ishiwatari  xy_SurfTemp の受授特性を out から inout
!              に変更. 
!
! TODO
! * 各種パラメータはファイルから読み込むようにもできるようにするべし
!   まず, ファイルからの読みこみ, そうでなければデフォルト値の
!   適用, という優先順位が良いように思われる.
! * 1 回設定すれば良い場合と毎ステップ設定しないといけない場合の処理は?

module physics_ground_mod

  use type_mod,    only : REKIND, DBKIND, INTKIND, TOKEN, STRING
  implicit none
  private
  public :: physics_ground

  logical, save :: physics_ground_initialized = .false.

contains

  subroutine physics_ground( &
    & xy_SurfTemp   ,& ! (inout) 地表温度
    & xy_SurfAlbedo, &        !(out) 地表アルベド
    & xy_SurfHumidCoeff, &    !(out) 地表湿潤度
    & xy_SurfRoughLength, &   !(out) 地表粗度長
    & xy_SurfCondition, &     !(out) 地表状態. 
    & xy_SurfHeatCapacity, &  !(out) 地表熱容量
    & xy_GroundTempFlux, &    !(out) 地中熱フラックス
    &  y_Lat         ) ! (in) 緯度座標
    !
    != 地表面パラメータの設定
    ! 
    != TODO (2006-8-16 石渡)
    ! * 他にも設定できるパラメータを増やすべきだ.
    !   例えば, AGCM5 では
    !       REAL       GRZSD  ( IDIM*JDIM )    !" 地表高度分散 = 0.0
    !   というのが設定できるようになっている.
    ! * 設定方法として何を用意するべきか決めないといけない.
    !   考えられる方法は以下の通り.
    !   * namelist で読み込んだデフォルト値(single value)で全部埋める.
    !   * nc ファイルから分布を読み込む
    !   * モデルの中で分布を計算する
    !     例えば, SST の場合は SurfCond にも依存する, 
    !     アルベドは SST に依存するetc
    ! * パラメータ毎に設定サブルーチンを用意するべきかもしれない.
    !
    use type_mod,      only: REKIND, DBKIND, INTKIND, TOKEN, STRING
    use grid_3d_mod,   only: im, jm, km
    use nmlfile_mod, only : nmlfile_init, nmlfile_open, nmlfile_close
    use dc_trace,      only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump
    use dc_message , only : MessageNotify
 
    implicit none
    real(DBKIND), intent(in) :: y_Lat(:)               ! 緯度座標
    real(DBKIND), intent(inout) :: xy_SurfTemp(im,jm)       ! 地表面温度
    real(DBKIND), intent(out) :: xy_SurfAlbedo(im,jm)        ! 地表アルベド
    real(DBKIND), intent(out) :: xy_SurfHumidCoeff(im,jm)    ! 地表湿潤度
    real(DBKIND), intent(out) :: xy_SurfRoughLength(im,jm)   ! 地表粗度長
    real(DBKIND), intent(out) :: xy_SurfHeatCapacity(im,jm)  ! 地表熱容量
    real(DBKIND), intent(out) :: xy_GroundTempFlux(im,jm) ! 地中熱フラックス
    integer(INTKIND), intent(out) :: xy_SurfCondition(im,jm)    ! 地表状態
                                       ! 1 未満は SST 固定. 1 以上は swamp

    character(STRING),  parameter:: subname = "physics_ground"
    real(DBKIND) :: DefaultSurfAlbedo ! 地表アルベドデフォルト値
    real(DBKIND) :: DefaultSurfHumidCoeff ! 地表湿潤度デフォルト値
    real(DBKIND) :: DefaultSurfRoughLength ! 地表粗度長デフォルト値
    real(DBKIND) :: DefaultSurfHeatCapacity ! 地表熱容量デフォルト値
    real(DBKIND) :: DefaultGroundTempFlux ! 地中熱フラックスデフォルト値
    integer(INTKIND) :: DefaultSurfCondition ! 地表状態デフォルト値
    real(DBKIND) :: xy_SeaSurfaceTemp(im,jm)            ! 海表面温度
    integer(INTKIND)            :: nmlstat, nmlunit
    logical                     :: nmlreadable

    integer(INTKIND) :: i, j

    ! Hosaka et al. (1998) SST 用パラメータ
    !   (2006-8-16 石渡)  あくまで暫定処理
    real(DBKIND), parameter :: TEQ   = 302.0d0
    real(DBKIND), parameter :: ALAT0 =   0.0d0
    real(DBKIND), parameter :: ALAT1 =  30.0d0
    real(DBKIND), parameter :: ALPHA =  60.0d0
    real(DBKIND), parameter :: BETA  =  32.0d0
    real(DBKIND), parameter :: GAMMA =   0.0d0
    real(DBKIND), parameter :: ALACON = 7.0d0

    namelist /ground_nml/ &
      & DefaultSurfAlbedo, & ! 地表アルベドデフォルト値
      & DefaultSurfHumidCoeff, & ! 地表湿潤度デフォルト値
      & DefaultSurfRoughLength, & ! 地表粗度長デフォルト値
      & DefaultSurfHeatCapacity, & ! 地表熱容量デフォルト値
      & DefaultGroundTempFlux, & ! 地中熱フラックスデフォルト値
      & DefaultSurfCondition ! 地表状態デフォルト値

    continue

    !   開始処理
    call BeginSub(subname)

    ! (2007-5-17 石渡) 追加
!    if (physics_ground_initialized) then
!!      call EndSub( subname, '%c is already called.', c1=trim(subname) )
!      return
!       continue
!    else
!    endif

    ! 地表状態の設定
    DefaultSurfAlbedo = 0.0         ! 地表アルベド
    DefaultSurfHumidCoeff = 1.0      ! 地表湿潤度
    DefaultSurfRoughLength = 0.0001  ! 地表粗度長
    DefaultSurfHeatCapacity = 0.0    ! 地表熱容量
    DefaultGroundTempFlux  = 0.0     ! 地中熱フラックス
    DefaultSurfCondition = 0 ! 地表状態

    call nmlfile_init
    call nmlfile_open(nmlunit, nmlreadable)
    if (nmlreadable) then
       read(nmlunit, nml=ground_nml, iostat=nmlstat)
       call DbgMessage('Stat of NAMELIST ground_nml Input is <%d>', &
            &           i=(/nmlstat/))
       write(0, nml=ground_nml)
    else
       call DbgMessage('Not Read NAMELIST ground_nml')
       call MessageNotify('W', subname, &
            & 'Can not Read NAMELIST ground_nml. Force Use Default Value.')
    end if
    call nmlfile_close
    ! ファイルからも入力できるようにするべし.


    ! (2007-5-20 石渡) 追加
    ! その他の地表面パラメータの設定.
    ! 地表状態と同様に処理するべき.
    xy_SurfAlbedo = DefaultSurfAlbedo        
    xy_SurfHumidCoeff = DefaultSurfHumidCoeff
    xy_SurfRoughLength = DefaultSurfRoughLength
    xy_SurfHeatCapacity = DefaultSurfHeatCapacity
    xy_GroundTempFlux  = DefaultGroundTempFlux
    xy_SurfCondition = DefaultSurfCondition

    !----------------------------------------------------------------
    !   地表面温度の設定
    !----------------------------------------------------------------
    
    !----- 海面温度の設定 -----
    ! Hosaka et al. (1998) の SST 分布を与えてしまう.
    ! これは暫定的処置. 
    ! 本当は
    !   1. SST ファイルが指定されていればそれを入力
    !   2. そうでなければデフォルトの分布を与える
    ! とするべき.
    
    call mksst( &
         & xy_SeaSurfaceTemp                           , &
         & y_Lat                                       , &
         & TEQ , ALAT0 , ALAT1 , ALPHA , BETA  , GAMMA , &
         & ALACON                                       )

    !----- 地表面温度を海面温度で置き換え-----
    !   (2007-08-27 石渡) 
    !   * 次は SST の初期値の設定をしているだけ.
    !   * 本当は xy_SurfCondition を参照して SST 固定グリッドだけ
    !     SST 値を与えるようにしないといけない 
    !     (SST を固定する grid と SST を時間発展させる grid が
    !     混在していても可とするべき).
    !   * 現状では, physics_ground は SST の初期値の設定をしている
    !     だけであるが, いずれは SST 他の値に応じて表面アルベドの
    !     変えることをここでやれるようにしたい.

    if (physics_ground_initialized) then
    else
      xy_SurfTemp = xy_SeaSurfaceTemp
      physics_ground_initialized = .true.
    endif


    ! 終了処理
    call EndSub(subname)

  end subroutine physics_ground


  subroutine mksst( &
    & xy_SeaSurfaceTemp                           , &
    & y_Lat                                       , &
    & TEQ , ALAT0 , ALAT1 , ALPHA , BETA  , GAMMA , &
    & ALACON                                       )
    !
    use type_mod,      only: REKIND, DBKIND, INTKIND, TOKEN, STRING
    use grid_3d_mod,   only: im, jm, km
    use constants_mod, only: PI 
    use dc_trace,      only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump
    implicit none
    real(DBKIND), intent(in) ::   &
         &  y_Lat(:)             ,& ! (in) 経度座標
         &  TEQ , ALAT0 , ALAT1 , ALPHA , BETA  , GAMMA , ALACON
    real(DBKIND), intent(out) ::  &
         & xy_SeaSurfaceTemp(im,jm)       ! (out) 地表面温度
    character(STRING),  parameter:: subname = "mksst"
    integer(INTKIND)    :: i, j
            ! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用)
    real(DBKIND) :: PHI1, AB4, PHI, GSSTP, ALATP, ALATM, GSSTMX
    integer(INTKIND)    :: JP, JMM

    continue

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

    !----------------------------------------------------------------
    !   海表面温度の設定
    !----------------------------------------------------------------

    PHI1 = ABS( ALAT1 * PI/180. )
    AB4  = 2. *( PHI1**3 )*BETA/ALPHA
    
    do j = 1, jm
       PHI   = ABS( y_Lat( J )*PI/180. - ALAT0 *PI/180. )
       GSSTP = TEQ                                                      &
            &         - ALPHA/2.                                        &
            &           * ( PHI - MAX(   SQRT(  PHI1**2+AB4 )      &
            &                     - SQRT( (PHI-PHI1)**2+AB4 ), 0.  )  ) &
            &         + GAMMA *( PHI**3 )
       do i = 1, im
          xy_SeaSurfaceTemp(i,j) = GSSTP
       end do
    end do
        
    ! ---- 中心 ALAT0 +/- ALACON の間を平坦に -----
    ALATP = ( ALAT0 + ALACON )*PI/180.
    ALATM = ( ALAT0 - ALACON )*PI/180.

    JP = 1
    JMM = jm
    do  j = 1, jm
       IF ( (-y_Lat( J )*PI/180.) .GT. ALATP ) THEN
          JMM = J
       ENDIF
       IF ( (-y_Lat( J )*PI/180.) .GE. ALATM ) THEN
          JP = J
       ENDIF
    end do

    GSSTMX = ( xy_SeaSurfaceTemp(1,JP) * ( ALATP - y_Lat(JP+1)*PI/180.)     &
         &   +  xy_SeaSurfaceTemp(1,JP+1) * ( y_Lat(JP)*PI/180. - ALATP ))  &
         &        /( y_Lat(JP)*PI/180. - y_Lat(JP+1)*PI/180. ) 

    do j = JMM+1, JP 
       do i = 1, im 
          xy_SeaSurfaceTemp(i,j) = GSSTMX
       end do
    end do

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

  end subroutine mksst

end module physics_ground_mod




