!= 
!
!= Slab ocean sea ice utility module : 1D mass and theraml balance calculation
!
! Authors::   Yoshiyuki O. Takahashi
! Version::   $Id: phy_implicit_sdh_V3.f90,v 1.1 2026/06/05 10:00:00 takepiro Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008-2026. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module sosi_utils
  !
  != 
  !
  != Slab ocean sea ice utility module : 1D mass and theraml balance calculation
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  !== Procedures List
  !
  ! SOSIUtilsChkSOSeaIce         :: 海氷インデックスをチェックする
  ! SOSIUtilsSetSOSeaIceLevels   :: 海氷情報のインデックスを設定する
  ! SOSIUtilsSetMissingValue     :: 海氷がないところの海氷温度に欠損値を代入
  
  ! SOSIUtilsAddPhysics          :: 物理過程計算(鉛直 1 次元)
  ! SOSIUtilsInit                :: 初期化
  ! -------------------------------
  ! SOSIUtilsChkSOSeaIce         :: Check index for slab ocean sea ice information
  ! SOSIUtilsSetSOSeaIceLevels   :: Set index for Slab Ocean Sea Ice information
  ! SOSIUtilsSetMissingValue     :: Set missing values at the levels without sea ice
  ! SOSIUtilsAddPhysics          :: Calculate physical processes (vertical 1D)
  ! SOSIUtilsInit                :: Initialization
  !--
  !
  !== NAMELIST
  !
  ! NAMELIST#sosi_utils_nml
  !++

  ! モジュール引用 ; USE statements
  !

  ! 格子点設定
  ! Grid points settings
  !
  use gridset, only:   imax, & ! 経度格子点数. 
                               ! Number of grid points in longitude
    &                  jmax, & ! 緯度格子点数. 
                               ! Number of grid points in latitude
    &                  kmax, & ! 鉛直層数. 
                               ! Number of vertical level
    &                  kslmax, & ! 地下の鉛直層数. 
                               ! Number of subsurface vertical level
    &                  ksimax  ! 海氷の鉛直層数. 
                               ! Number of sea ice vertical level

  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! 倍精度実数型. Double precision. 
    &                 STRING     ! 文字列.       Strings. 

  ! メッセージ出力
  ! Message output
  !
  use dc_message, only: MessageNotify

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  real(DP), parameter :: SIDepthMergin = 1.0e-3_DP
                               ! Mergin for sea ice depth
                               ! This is used to avoid very thin sea ice layer.

  ! 公開手続き
  ! Public procedure
  !
  public :: SOSIUtilsChkSOSeaIce
  public :: SOSIUtilsSetSOSeaIceLevels
  public :: SOSIUtilsSetMissingValue
  public :: SOSIUtilsAddPhysics
  public :: SOSIUtilsInit


  ! 公開変数
  ! Public variables
  !
  real(DP), parameter, public :: SOSeaIceTempMissingValue = -99999.0_DP


  real(DP), save, public :: SOSeaIceMassNegativeThreshold


  logical, save :: sosi_utils_inited = .false.
                              ! 初期設定フラグ. 
                              ! Initialization flag

  character(*), parameter:: module_name = 'sosi_utils'
                              ! モジュールの名称. 
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: sosi_utils.f90,v 2.0 2026/06/05 10:00:00 takepiro Exp $'
                              ! モジュールのバージョン
                              ! Module version

contains

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

  subroutine SOSIUtilsChkSOSeaIce(          &
    & xy_SeaIceThickness, xyz_SOSeaIceTemp, & ! (in)
    & ParentRoutine                         & ! (in)
    & )
    !
    ! 海氷インデックスをチェックする
    ! Check index for slab ocean sea ice information
    !

    ! モジュール引用 ; USE statements
    !

    ! 宣言文 ; Declaration statements
    !
    real(DP), intent(in ) :: xy_SeaIceThickness(0:imax-1, 1:jmax)
    real(DP), intent(in ) :: xyz_SOSeaIceTemp  (0:imax-1, 1:jmax, 1:ksimax)
    character(*), intent(in ), optional :: ParentRoutine


    ! 作業変数
    ! Work variables
    !
    integer  :: xy_SOSILocalKMax(0:imax-1, 1:jmax)
    real(DP) :: xyr_SOSILocalDepth(0:imax-1, 1:jmax, 0:ksimax)
    real(DP) :: xyz_SOSILocalDepth(0:imax-1, 1:jmax, 1:ksimax)

    integer:: i               ! 経度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in longitude
    integer:: j               ! 緯度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in latitude
    integer:: k


    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. sosi_utils_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    call SOSIUtilsSetSOSeaIceLevels(                     &
    & xy_SeaIceThickness,                                       & ! (in   )
    & xy_SOSILocalKMax, xyr_SOSILocalDepth, xyz_SOSILocalDepth  & ! (out)
    & )


    do i = 0, imax-1
      do j = 1, jmax
        do k = 1, xy_SOSILocalKMax(i,j)
          if ( xyz_SOSeaIceTemp(i,j,k) < 0.0_DP ) then
            if ( present( ParentRoutine ) ) then
              call MessageNotify( 'M', module_name, &
                & 'Called from %c:', &
                & c1 = trim( ParentRoutine ) )
            end if
            call MessageNotify( 'M', module_name, &
              & 'xyz_SOSeaIceTemp(%d,%d,%d) = %f.', &
              & i = (/i,j,k/), d = (/xyz_SOSeaIceTemp(i,j,k)/) )
          end if
        end do
      end do
    end do


  end subroutine SOSIUtilsChkSOSeaIce

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

  subroutine SOSIUtilsSetSOSeaIceLevels(                     &
    & xy_SeaIceThickness,                                       & ! (in   )
    & xy_SOSILocalKMax, xyr_SOSILocalDepth, xyz_SOSILocalDepth  & ! (out)
    & )
    !
    ! スラブオーシャンの海氷の深さ情報を設定する. 
    ! Set depth information for Slab Ocean Sea Ice.
    !

    ! モジュール引用 ; USE statements
    !

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: &
      & r_SIDepth         ! sea ice grid on interface of layer


    ! 宣言文 ; Declaration statements
    !
    real(DP), intent(in ) :: xy_SeaIceThickness(0:imax-1, 1:jmax)
                              ! 海氷厚さ
                              ! Sea ice thickness
    integer , intent(out) :: xy_SOSILocalKMax  (0:imax-1, 1:jmax)
                              ! 鉛直最大インデックス
                              ! Max. vertical index
    real(DP), intent(out) :: xyr_SOSILocalDepth(0:imax-1, 1:jmax, 0:ksimax)
                              ! 各レベルでの深さ  (\simga 整数レベル)
                              ! Depth at each level (\simga integer)
    real(DP), intent(out) :: xyz_SOSILocalDepth(0:imax-1, 1:jmax, 1:ksimax)
                              ! 各レベルでの深さ  (\simga 半整数レベル)
                              ! Depth at each level (\simga half-integer)

    ! 作業変数
    ! Work variables
    !

    integer:: i               ! 経度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in longitude
    integer:: j               ! 緯度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in latitude
    integer:: k


    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. sosi_utils_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    do i = 0, imax-1
      do j = 1, jmax
        if ( xy_SeaIceThickness(i,j) == 0.0_DP ) then
          xy_SOSILocalKMax(i,j) = 0
        else if ( - xy_SeaIceThickness(i,j) < r_SIDepth(ksimax) ) then
          xy_SOSILocalKMax(i,j) = ksimax
        else
          xy_SOSILocalKMax(i,j) = 0
          search_ksimax : do k = 0+1, ksimax
!!$            if ( - xy_SeaIceThickness(i,j) >= r_SIDepth(k) ) then
            ! This SIDepthMergin avoids very thin lowest layer.
            if ( - xy_SeaIceThickness(i,j) >= r_SIDepth(k)-SIDepthMergin ) then
                xy_SOSILocalKMax(i,j) = k
              exit search_ksimax
            end if
          end do search_ksimax
        end if
      end do
    end do
    do j = 1, jmax
      do i = 0, imax-1
        do k = 0, xy_SOSILocalKMax(i,j)-1
          xyr_SOSILocalDepth(i,j,k) = r_SIDepth(k)
        end do
        k = xy_SOSILocalKMax(i,j)
        xyr_SOSILocalDepth(i,j,k) = - xy_SeaIceThickness(i,j)
        do k = xy_SOSILocalKMax(i,j)+1, ksimax
          xyr_SOSILocalDepth(i,j,k) = -1.0e100_DP
        end do
        !
        do k = 1, xy_SOSILocalKMax(i,j)
          xyz_SOSILocalDepth(i,j,k) = &
            & ( xyr_SOSILocalDepth(i,j,k-1) + xyr_SOSILocalDepth(i,j,k) ) / 2.0_DP
        end do
        do k = xy_SOSILocalKMax(i,j)+1, ksimax
          xyz_SOSILocalDepth(i,j,k) = -1.0e100_DP
        end do
      end do
    end do


  end subroutine SOSIUtilsSetSOSeaIceLevels

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

  subroutine SOSIUtilsSetMissingValue( &
    & xy_SOSeaIceMass,                 & !(in )
    & xyz_SOSeaIceTemp,                & !(inout)
    & SOSeaIceValue                    & !(in ) optional
    & )
    !
    ! 海氷がないところの海氷温度に欠損値を代入
    ! Set missing values at the levels without sea ice
    
    ! モジュール引用 ; USE statements
    !

    ! 雪と海氷の定数の設定
    ! Setting constants of snow and sea ice
    !
    use constants_snowseaice, only: &
      & TempBelowSeaIce, &
      & SeaIceDen


    ! 宣言文 ; Declaration statements
    !
    real(DP), intent(in ) :: xy_SOSeaIceMass    (0:imax-1, 1:jmax)
                              ! $ M_si $ . 海氷質量 (kg m-2)
                              ! Slab ocean sea ice mass (kg m-2)
    real(DP), intent(inout) :: xyz_SOSeaIceTemp      (0:imax-1, 1:jmax, 1:ksimax)
    real(DP), intent(in   ), optional :: SOSeaIceValue


    ! 作業変数
    ! Work variables
    !
    real(DP) :: xy_SeaIceThickness(0:imax-1, 1:jmax)
                 !
                 ! Sea ice thickness
    integer  :: xy_SOSILocalKMax  (0:imax-1, 1:jmax)
    real(DP) :: xyr_SOSILocalDepth(0:imax-1, 1:jmax, 0:ksimax)
    real(DP) :: xyz_SOSILocalDepth(0:imax-1, 1:jmax, 1:ksimax)

    real(DP) :: xy_SeaIceThicknessA(0:imax-1, 1:jmax)
                 !
                 ! Sea ice thickness
    integer  :: xy_SOSILocalKMaxA  (0:imax-1, 1:jmax)
    real(DP) :: xyr_SOSILocalDepthA(0:imax-1, 1:jmax, 0:ksimax)
    real(DP) :: xyz_SOSILocalDepthA(0:imax-1, 1:jmax, 1:ksimax)

    real(DP) :: SetValue

    integer:: i               ! 東西方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in zonal direction
    integer:: j               ! 南北方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in meridional direction
    integer:: k


    if ( .not. sosi_utils_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    ! 海氷の厚さ計算
    ! Calcuate sea ice thickness
    !
    xy_SeaIceThickness = xy_SOSeaIceMass / SeaIceDen
    !
    ! 海氷の深さレベルを設定
    ! Set slab ocean sea ice levels
    !
    call SOSIUtilsSetSOSeaIceLevels(                              &
      & xy_SeaIceThickness,                                       & ! (in   )
      & xy_SOSILocalKMax, xyr_SOSILocalDepth, xyz_SOSILocalDepth  & ! (out)
      & )


    if ( present ( SOSeaIceValue ) ) then
      SetValue = SOSeaIceValue
    else
      SetValue = SOSeaIceTempMissingValue
    end if
    !
    ! 海氷がないところの海氷温度に欠損値を代入
    ! Set missing values at the levels without sea ice
    !
    do j = 1, jmax
      do i = 0, imax-1
        do k = xy_SOSILocalKMax(i,j)+1, ksimax
          xyz_SOSeaIceTemp(i,j,k) = SetValue
        end do
      end do
    end do

  end subroutine SOSIUtilsSetMissingValue

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

  subroutine SOSIUtilsAddPhysics(                         &
    & xy_SOSeaIceMass, xyz_SOSeaIceTemp,                  & !(inout)
    & xy_DSOSeaIceMassDtPhyTop, xy_DSOSeaIceMassDtPhyBot, & !(in   )
    & xyz_DSOSeaIceTempDtPhy                              & !(in   )
    & )
    !
    ! スラブオーシャン海氷の物理過程計算(鉛直 1 次元)
    ! Calculate physical processes of slab ocean sea ice.

    ! モジュール引用 ; USE statements
    !

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut

    use timeset    , only : &
      & TimeN, &
      & DelTime
                              ! $\Delta t$

    ! 雪と海氷の定数の設定
    ! Setting constants of snow and sea ice
    !
    use constants_snowseaice, only: &
      & TempBelowSeaIce, &
      & SeaIceDen


    ! 宣言文 ; Declaration statements
    !
    real(DP), intent(inout) :: xy_SOSeaIceMass    (0:imax-1, 1:jmax)
                              ! $ M_si $ . 海氷質量 (kg m-2)
                              ! Slab ocean sea ice mass (kg m-2)
    real(DP), intent(inout) :: xyz_SOSeaIceTemp      (0:imax-1, 1:jmax, 1:ksimax)
    real(DP), intent(in   ) :: xy_DSOSeaIceMassDtPhyTop(0:imax-1, 1:jmax)
    real(DP), intent(in   ) :: xy_DSOSeaIceMassDtPhyBot(0:imax-1, 1:jmax)

    real(DP), intent(in   ) :: xyz_DSOSeaIceTempDtPhy(0:imax-1, 1:jmax, 1:ksimax)


    ! 作業変数
    ! Work variables
    !
    real(DP) :: xy_SeaIceThickness(0:imax-1, 1:jmax)
    integer  :: xy_SOSILocalKMax  (0:imax-1, 1:jmax)
    real(DP) :: xyr_SOSILocalDepth(0:imax-1, 1:jmax, 0:ksimax)
    real(DP) :: xyz_SOSILocalDepth(0:imax-1, 1:jmax, 1:ksimax)

    ! 海氷下面からの変化
    real(DP) :: xy_SOSeaIceMassTent1   (0:imax-1, 1:jmax)
    real(DP) :: xy_SeaIceThicknessTent1(0:imax-1, 1:jmax)
    integer  :: xy_SOSILocalKMaxTent1  (0:imax-1, 1:jmax)
    real(DP) :: xyr_SOSILocalDepthTent1(0:imax-1, 1:jmax, 0:ksimax)
    real(DP) :: xyz_SOSILocalDepthTent1(0:imax-1, 1:jmax, 1:ksimax)

    ! 海氷上面からの変化
    real(DP) :: xy_SOSeaIceMassTent2   (0:imax-1, 1:jmax)
    real(DP) :: xy_SeaIceThicknessTent2(0:imax-1, 1:jmax)
    integer  :: xy_SOSILocalKMaxTent2  (0:imax-1, 1:jmax)
    real(DP) :: xyr_SOSILocalDepthTent2(0:imax-1, 1:jmax, 0:ksimax)
    real(DP) :: xyz_SOSILocalDepthTent2(0:imax-1, 1:jmax, 1:ksimax)

    real(DP) :: xyr_SOSILocalDepthT2MapToT1(0:imax-1, 1:jmax, 0:ksimax)

    real(DP) :: xyz_SOSeaIceTempTent1(0:imax-1, 1:jmax, 1:ksimax)
    real(DP) :: xyz_SOSeaIceTempTent2(0:imax-1, 1:jmax, 1:ksimax)

    integer:: i               ! 東西方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in zonal direction
    integer:: j               ! 南北方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in meridional direction
    integer:: k
    integer:: kk
    integer:: kTop
    integer:: kBot


    if ( .not. sosi_utils_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    ! Calculate sea ice mass at next time step
    !

    !   Add sea ice mass change at the bottom of sea ice
    xy_SOSeaIceMassTent1 = xy_SOSeaIceMass &
      & + xy_DSOSeaIceMassDtPhyBot * DelTime

    !   Add sea ice mass change at the top of sea ice
    xy_SOSeaIceMassTent2 = xy_SOSeaIceMassTent1 &
      & + xy_DSOSeaIceMassDtPhyTop * DelTime


    ! 海氷温度時間積分
    ! Time integration of sea ice temperature
    !
    xyz_SOSeaIceTempTent1 = xyz_SOSeaIceTemp &
      & + xyz_DSOSeaIceTempDtPhy * DelTime


    !
    ! Adjust temperature
    !

    !
    !   Calcuate sea ice thickness
    !
    xy_SeaIceThickness = xy_SOSeaIceMass / SeaIceDen
    !
    !   Set slab ocean sea ice levels
    !
    call SOSIUtilsSetSOSeaIceLevels(                              &
      & xy_SeaIceThickness,                                       & ! (in )
      & xy_SOSILocalKMax, xyr_SOSILocalDepth, xyz_SOSILocalDepth  & ! (out)
      & )

    !
    !   Calcuate sea ice thickness
    !
    xy_SeaIceThicknessTent1 = xy_SOSeaIceMassTent1 / SeaIceDen
    !
    !   Set slab ocean sea ice levels
    !
    call SOSIUtilsSetSOSeaIceLevels(                     &
      & xy_SeaIceThicknessTent1,                                       & ! (in   )
      & xy_SOSILocalKMaxTent1, xyr_SOSILocalDepthTent1, xyz_SOSILocalDepthTent1  & ! (out)
      & )

    !   Adjust levels
    !
    do j = 1, jmax
      do i = 0, imax-1
        if ( xy_SOSILocalKMaxTent1(i,j) > xy_SOSILocalKMax(i,j) ) then
          ! sea ice thickness increases
          if ( xy_SOSILocalKMax(i,j) == 0 ) then
            do k = 1, xy_SOSILocalKMaxTent1(i,j)
              xyz_SOSeaIceTempTent1(i,j,k) = TempBelowSeaIce
            end do
          else
!!$            do k = 1, xy_SOSILocalKMaxB(i,j)
!!$              ! Do nothing
!!$            end do
            do k = xy_SOSILocalKMax(i,j)+1, xy_SOSILocalKMaxTent1(i,j)
              kk = xy_SOSILocalKMax(i,j)
              xyz_SOSeaIceTempTent1(i,j,k) = xyz_SOSeaIceTempTent1(i,j,kk)
            end do
          end if
        else if ( xy_SOSILocalKMaxTent1(i,j) < xy_SOSILocalKMax(i,j) ) then
          ! sea ice thickness decreases
          !   Do nothing
          !   Melted sea ice had freezing temperature
        end if
      end do
    end do



    !
    !   Calcuate sea ice thickness
    !
    xy_SeaIceThicknessTent2 = xy_SOSeaIceMassTent2 / SeaIceDen
    !
    !   Set slab ocean sea ice levels
    !
    call SOSIUtilsSetSOSeaIceLevels(                     &
      & xy_SeaIceThicknessTent2,                                       & ! (in   )
      & xy_SOSILocalKMaxTent2, xyr_SOSILocalDepthTent2, xyz_SOSILocalDepthTent2  & ! (out)
      & )

    !   Adjust levels
    !
    do j = 1, jmax
      do i = 0, imax-1
        if ( xy_SOSILocalKMaxTent2(i,j) == 0 ) then
          do k = 0, ksimax
            xyr_SOSILocalDepthT2MapToT1(i,j,k) = 0.0_DP
          end do
        else
          do k = 0, xy_SOSILocalKMaxTent2(i,j)
            xyr_SOSILocalDepthT2MapToT1(i,j,k) = &
              &   xyr_SOSILocalDepthTent2(i,j,k) &
              & - ( xy_SeaIceThicknessTent1(i,j) - xy_SeaIceThicknessTent2(i,j) )
          end do
          do k = xy_SOSILocalKMaxTent2(i,j)+1, ksimax
            xyr_SOSILocalDepthT2MapToT1(i,j,k) = 0.0_DP
          end do
        end if
      end do
    end do
    !
    do j = 1, jmax
      do i = 0, imax-1

        if ( xy_SOSILocalKMaxTent2(i,j) == 0 ) then
          do k = 1, ksimax
            xyz_SOSeaIceTempTent2(i,j,k) = SOSeaIceTempMissingValue
          end do
        else

          do k = 1, xy_SOSILocalKMaxTent2(i,j)

            search_kTop : do kTop = 1, xy_SOSILocalKMaxTent1(i,j)
              if ( xyr_SOSILocalDepthT2MapToT1(i,j,k-1) >= xyr_SOSILocalDepthTent1(i,j,kTop) ) then
                exit search_kTop
              end if
            end do search_kTop
            search_kBot : do kBot = kTop, xy_SOSILocalKMaxTent1(i,j)
              if ( xyr_SOSILocalDepthT2MapToT1(i,j,k  ) >= xyr_SOSILocalDepthTent1(i,j,kBot) ) then
                exit search_kBot
              end if
            end do search_kBot

            if ( kTop == kBot ) then
              kk = kTop
              xyz_SOSeaIceTempTent2(i,j,k) = &
                & + xyz_SOSeaIceTempTent1(i,j,kk)                &
                &     * (   xyr_SOSILocalDepthT2MapToT1(i,j,k-1) &
                &         - xyr_SOSILocalDepthT2MapToT1(i,j,k  ) )
            else
              xyz_SOSeaIceTempTent2(i,j,k) = 0.0_DP
              kk = kTop
              xyz_SOSeaIceTempTent2(i,j,k) = xyz_SOSeaIceTempTent2(i,j,k) &
                & + xyz_SOSeaIceTempTent1(i,j,kk)                      &
                &     * (   xyr_SOSILocalDepthT2MapToT1(i,j,k-1) &
                &         - xyr_SOSILocalDepthTent1(i,j,kk) )
              do kk = kTop+1, kBot-1
                xyz_SOSeaIceTempTent2(i,j,k) = xyz_SOSeaIceTempTent2(i,j,k) &
                  & + xyz_SOSeaIceTempTent1(i,j,kk)                      &
                  &     * (   xyr_SOSILocalDepthTent1(i,j,kk-1)          &
                  &         - xyr_SOSILocalDepthTent1(i,j,kk  ) )
              end do
              kk = kBot
              xyz_SOSeaIceTempTent2(i,j,k) = xyz_SOSeaIceTempTent2(i,j,k) &
                & + xyz_SOSeaIceTempTent1(i,j,kk)                      &
                &     * (   xyr_SOSILocalDepthTent1(i,j,kk-1)          &
                &         - xyr_SOSILocalDepthT2MapToT1(i,j,k) )
            end if
            !
            xyz_SOSeaIceTempTent2(i,j,k) = &
              & xyz_SOSeaIceTempTent2(i,j,k) &
              & / ( xyr_SOSILocalDepthTent2(i,j,k-1) - xyr_SOSILocalDepthTent2(i,j,k) )

          end do
          do k = xy_SOSILocalKMaxTent2(i,j)+1, ksimax
            xyz_SOSeaIceTempTent2(i,j,k) = SOSeaIceTempMissingValue
          end do

        end if

      end do
    end do


    ! Update sea ice temperature
    !
    xyz_SOSeaIceTemp = xyz_SOSeaIceTempTent2

    ! Update sea ice mass
    !
    xy_SOSeaIceMass = xy_SOSeaIceMassTent2


    ! Check sea ice mass
    !
    do j = 1, jmax
      do i = 0, imax-1
        if ( xy_SOSeaIceMass(i,j) < 0 ) then
          if ( xy_SOSeaIceMass(i,j) < SOSeaIceMassNegativeThreshold ) then
            call MessageNotify( 'M', module_name, &
              & '  Slab sea ice mass is negative after physics, %f, and this is set to zero.', &
              & d = (/ xy_SOSeaIceMass(i,j) /) )
          end if
          xy_SOSeaIceMass(i,j) = 0.0_DP
        end if
      end do
    end do


    ! Check
    !
    xy_SeaIceThickness = xy_SOSeaIceMass / SeaIceDen
    !
    call SOSIUtilsChkSOSeaIce(          &
      & xy_SeaIceThickness, xyz_SOSeaIceTemp, & ! (in)
      & "SOSIUtilsAddPhysics"               & ! (in)
      & )


  end subroutine SOSIUtilsAddPhysics

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

!!$  subroutine BK_SOSIUtilsAddPhysics(                         &
!!$    & xy_SOSeaIceMassB,                                   & !(in )
!!$    & xy_DSOSeaIceMassDtPhyTop, xy_DSOSeaIceMassDtPhyBot, & !(in )
!!$    & xy_SOSeaIceMassA,                                   & !(out)
!!$    & xyz_SOSeaIceTemp,                                   & !(inout)
!!$    & xyz_DSOSeaIceTempDtPhy                              & !(in   )
!!$    & )
!!$    ! 
!!$    ! Calculates slab sea ice horizontal transports by diffusion
!!$
!!$    ! ヒストリデータ出力
!!$    ! History data output
!!$    !
!!$    use gtool_historyauto, only: HistoryAutoPut
!!$
!!$    use timeset    , only : &
!!$      & TimeN, &
!!$      & DelTime
!!$                              ! $\Delta t$
!!$
!!$    ! 雪と海氷の定数の設定
!!$    ! Setting constants of snow and sea ice
!!$    !
!!$    use constants_snowseaice, only: &
!!$      & TempBelowSeaIce, &
!!$      & SeaIceDen
!!$
!!$
!!$    ! 宣言文 ; Declaration statements
!!$    !
!!$    real(DP), intent(in ) :: xy_SOSeaIceMassB    (0:imax-1, 1:jmax)
!!$                              ! $ M_si (t-1) $ . 海氷質量 (kg m-2)
!!$                              ! Slab ocean sea ice mass (kg m-2)
!!$    real(DP), intent(in ) :: xy_DSOSeaIceMassDtPhyTop(0:imax-1, 1:jmax)
!!$    real(DP), intent(in ) :: xy_DSOSeaIceMassDtPhyBot(0:imax-1, 1:jmax)
!!$
!!$    real(DP), intent(out) :: xy_SOSeaIceMassA(0:imax-1, 1:jmax)
!!$                              !
!!$                              ! Slab sea ice mass at next time step
!!$    real(DP), intent(inout) :: xyz_SOSeaIceTemp      (0:imax-1, 1:jmax, 1:ksimax)
!!$    real(DP), intent(in   ) :: xyz_DSOSeaIceTempDtPhy(0:imax-1, 1:jmax, 1:ksimax)
!!$
!!$
!!$    ! 作業変数
!!$    ! Work variables
!!$    !
!!$    real(DP) :: xy_SeaIceThicknessB(0:imax-1, 1:jmax)
!!$    integer  :: xy_SOSILocalKMaxB  (0:imax-1, 1:jmax)
!!$    real(DP) :: xyr_SOSILocalDepthB(0:imax-1, 1:jmax, 0:ksimax)
!!$    real(DP) :: xyz_SOSILocalDepthB(0:imax-1, 1:jmax, 1:ksimax)
!!$
!!$    real(DP) :: xy_SOSeaIceMassA1(0:imax-1, 1:jmax)
!!$    real(DP) :: xy_SeaIceThicknessA1(0:imax-1, 1:jmax)
!!$    integer  :: xy_SOSILocalKMaxA1  (0:imax-1, 1:jmax)
!!$    real(DP) :: xyr_SOSILocalDepthA1(0:imax-1, 1:jmax, 0:ksimax)
!!$    real(DP) :: xyz_SOSILocalDepthA1(0:imax-1, 1:jmax, 1:ksimax)
!!$
!!$    real(DP) :: xy_SOSeaIceMassA2(0:imax-1, 1:jmax)
!!$    real(DP) :: xy_SeaIceThicknessA2(0:imax-1, 1:jmax)
!!$    integer  :: xy_SOSILocalKMaxA2  (0:imax-1, 1:jmax)
!!$    real(DP) :: xyr_SOSILocalDepthA2(0:imax-1, 1:jmax, 0:ksimax)
!!$    real(DP) :: xyz_SOSILocalDepthA2(0:imax-1, 1:jmax, 1:ksimax)
!!$
!!$    real(DP) :: xyr_SOSILocalDepthA2MappedToA1(0:imax-1, 1:jmax, 0:ksimax)
!!$
!!$    real(DP) :: xyz_SOSeaIceTempA1(0:imax-1, 1:jmax, 1:ksimax)
!!$    real(DP) :: xyz_SOSeaIceTempA2(0:imax-1, 1:jmax, 1:ksimax)
!!$
!!$    real(DP) :: xyz_DSOSeaIceTempDtPhyUpdt(0:imax-1, 1:jmax, 1:ksimax)
!!$
!!$    integer:: i               ! 東西方向に回る DO ループ用作業変数
!!$                              ! Work variables for DO loop in zonal direction
!!$    integer:: j               ! 南北方向に回る DO ループ用作業変数
!!$                              ! Work variables for DO loop in meridional direction
!!$    integer:: k
!!$    integer:: kk
!!$    integer:: kTop
!!$    integer:: kBot
!!$
!!$
!!$    if ( .not. sosi_utils_inited ) then
!!$      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
!!$    end if
!!$
!!$
!!$    ! Calculate sea ice mass at next time step
!!$    !
!!$
!!$    !   Add sea ice mass change at the bottom of sea ice
!!$    xy_SOSeaIceMassA1 = xy_SOSeaIceMassB &
!!$      & + xy_DSOSeaIceMassDtPhyBot * ( 2.0_DP * DelTime )
!!$
!!$    !   Add sea ice mass change at the top of sea ice
!!$    xy_SOSeaIceMassA2 = xy_SOSeaIceMassA1 &
!!$      & + xy_DSOSeaIceMassDtPhyTop * ( 2.0_DP * DelTime )
!!$
!!$    !   Update sea ice mass
!!$    xy_SOSeaIceMassA = xy_SOSeaIceMassA2
!!$
!!$
!!$    ! Check sea ice mass
!!$    !
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        if ( xy_SOSeaIceMassA(i,j) < 0 ) then
!!$          if ( xy_SOSeaIceMassA(i,j) < SOSeaIceMassNegativeThreshold ) then
!!$            call MessageNotify( 'M', module_name, &
!!$              & '  Slab sea ice mass is negative after physics, %f, and this is set to zero.', &
!!$              & d = (/ xy_SOSeaIceMassA(i,j) /) )
!!$          end if
!!$          xy_SOSeaIceMassA(i,j) = 0.0_DP
!!$        end if
!!$      end do
!!$    end do
!!$
!!$
!!$    ! 海氷温度時間積分
!!$    ! Time integration of sea ice temperature
!!$    !
!!$    xyz_SOSeaIceTemp = xyz_SOSeaIceTemp &
!!$      & + xyz_DSOSeaIceTempDtPhy * DelTime
!!$
!!$    !
!!$    ! Calcuate sea ice thickness
!!$    !
!!$
!!$    !
!!$    !   Calcuate sea ice thickness
!!$    !
!!$    xy_SeaIceThicknessB = xy_SOSeaIceMassB / SeaIceDen
!!$    !
!!$    !   Set slab ocean sea ice levels
!!$    !
!!$    call SOSIUtilsSetSOSeaIceLevels(                                 &
!!$      & xy_SeaIceThicknessB,                                         & ! (in )
!!$      & xy_SOSILocalKMaxB, xyr_SOSILocalDepthB, xyz_SOSILocalDepthB  & ! (out)
!!$      & )
!!$
!!$    !
!!$    !   Calcuate sea ice thickness
!!$    !
!!$    xy_SeaIceThicknessA1 = xy_SOSeaIceMassA1 / SeaIceDen
!!$    !
!!$    !   Set slab ocean sea ice levels
!!$    !
!!$    call SOSIUtilsSetSOSeaIceLevels(                     &
!!$      & xy_SeaIceThicknessA1,                                       & ! (in   )
!!$      & xy_SOSILocalKMaxA1, xyr_SOSILocalDepthA1, xyz_SOSILocalDepthA1  & ! (out)
!!$      & )
!!$
!!$    !   Adjust levels
!!$    !
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        if ( xy_SOSILocalKMaxA1(i,j) > xy_SOSILocalKMaxB(i,j) ) then
!!$          ! sea ice thickness increases
!!$          if ( xy_SOSILocalKMaxB(i,j) == 0 ) then
!!$            do k = 1, xy_SOSILocalKMaxA1(i,j)
!!$              xyz_SOSeaIceTempA1(i,j,k) = TempBelowSeaIce
!!$            end do
!!$          else
!!$!            do k = 1, xy_SOSILocalKMaxB(i,j)
!!$!              ! Do nothing
!!$!            end do
!!$            do k = xy_SOSILocalKMaxB(i,j)+1, xy_SOSILocalKMaxA1(i,j)
!!$              kk = xy_SOSILocalKMaxB(i,j)
!!$              xyz_SOSeaIceTempA1(i,j,k) = xyz_SOSeaIceTempA1(i,j,kk)
!!$            end do
!!$          end if
!!$        else if ( xy_SOSILocalKMaxA1(i,j) < xy_SOSILocalKMaxB(i,j) ) then
!!$          ! sea ice thickness decreases
!!$          !   Do nothing
!!$          !   Melted sea ice had freezing temperature
!!$        end if
!!$      end do
!!$    end do
!!$
!!$
!!$
!!$    !
!!$    !   Calcuate sea ice thickness
!!$    !
!!$    xy_SeaIceThicknessA2 = xy_SOSeaIceMassA2 / SeaIceDen
!!$    !
!!$    !   Set slab ocean sea ice levels
!!$    !
!!$    call SOSIUtilsSetSOSeaIceLevels(                     &
!!$      & xy_SeaIceThicknessA2,                                       & ! (in   )
!!$      & xy_SOSILocalKMaxA2, xyr_SOSILocalDepthA2, xyz_SOSILocalDepthA2  & ! (out)
!!$      & )
!!$
!!$    !   Adjust levels
!!$    !
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        if ( xy_SOSILocalKMaxA2(i,j) == 0 ) then
!!$          do k = 0, ksimax
!!$            xyr_SOSILocalDepthA2MappedToA1(i,j,k) = 0.0_DP
!!$          end do
!!$        else
!!$          do k = 0, xy_SOSILocalKMaxA2(i,j)
!!$            xyr_SOSILocalDepthA2MappedToA1(i,j,k) = &
!!$              &   xyr_SOSILocalDepthA2(i,j,k) &
!!$              & - ( xy_SeaIceThicknessA1(i,j) - xy_SeaIceThicknessA2(i,j) )
!!$          end do
!!$          do k = xy_SOSILocalKMaxA2(i,j)+1, ksimax
!!$            xyr_SOSILocalDepthA2MappedToA1(i,j,k) = 0.0_DP
!!$          end do
!!$        end if
!!$      end do
!!$    end do
!!$    !
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$
!!$        if ( xy_SOSILocalKMaxA2(i,j) == 0 ) then
!!$          do k = 1, ksimax
!!$            xyz_SOSeaIceTempA2(i,j,k) = SOSeaIceTempMissingValue
!!$          end do
!!$        else
!!$
!!$          do k = 1, xy_SOSILocalKMaxA2(i,j)
!!$
!!$            search_kTop : do kTop = 1, xy_SOSILocalKMaxA1(i,j)
!!$              if ( xyr_SOSILocalDepthA2MappedToA1(i,j,k-1) >= xyr_SOSILocalDepthA1(i,j,kTop) ) then
!!$                exit search_kTop
!!$              end if
!!$            end do search_kTop
!!$            search_kBot : do kBot = kTop, xy_SOSILocalKMaxA1(i,j)
!!$              if ( xyr_SOSILocalDepthA2MappedToA1(i,j,k  ) >= xyr_SOSILocalDepthA1(i,j,kBot) ) then
!!$                exit search_kBot
!!$              end if
!!$            end do search_kBot
!!$
!!$            if ( kTop == kBot ) then
!!$              kk = kTop
!!$              xyz_SOSeaIceTempA2(i,j,k) = &
!!$                & + xyz_SOSeaIceTempA1(i,j,kk)                      &
!!$                &     * (   xyr_SOSILocalDepthA2MappedToA1(i,j,k-1) &
!!$                &         - xyr_SOSILocalDepthA2MappedToA1(i,j,k  ) )
!!$            else
!!$              xyz_SOSeaIceTempA2(i,j,k) = 0.0_DP
!!$              kk = kTop
!!$              xyz_SOSeaIceTempA2(i,j,k) = xyz_SOSeaIceTempA2(i,j,k) &
!!$                & + xyz_SOSeaIceTempA1(i,j,kk)                      &
!!$                &     * (   xyr_SOSILocalDepthA2MappedToA1(i,j,k-1) &
!!$                &         - xyr_SOSILocalDepthA1(i,j,kk) )
!!$              do kk = kTop+1, kBot-1
!!$                xyz_SOSeaIceTempA2(i,j,k) = xyz_SOSeaIceTempA2(i,j,k) &
!!$                  & + xyz_SOSeaIceTempA1(i,j,kk)                      &
!!$                  &     * (   xyr_SOSILocalDepthA1(i,j,kk-1)          &
!!$                  &         - xyr_SOSILocalDepthA1(i,j,kk  ) )
!!$              end do
!!$              kk = kBot
!!$              xyz_SOSeaIceTempA2(i,j,k) = xyz_SOSeaIceTempA2(i,j,k) &
!!$                & + xyz_SOSeaIceTempA1(i,j,kk)                      &
!!$                &     * (   xyr_SOSILocalDepthA1(i,j,kk-1)          &
!!$                &         - xyr_SOSILocalDepthA2MappedToA1(i,j,k) )
!!$            end if
!!$            !
!!$            xyz_SOSeaIceTempA2(i,j,k) = &
!!$              & xyz_SOSeaIceTempA2(i,j,k) &
!!$              & / ( xyr_SOSILocalDepthA2(i,j,k-1) - xyr_SOSILocalDepthA2(i,j,k) )
!!$
!!$          end do
!!$          do k = xy_SOSILocalKMaxA2(i,j)+1, ksimax
!!$            xyz_SOSeaIceTempA2(i,j,k) = SOSeaIceTempMissingValue
!!$          end do
!!$
!!$        end if
!!$
!!$      end do
!!$    end do
!!$
!!$
!!$    !   Calculation of updated tendency
!!$    !
!!$!    do j = 1, jmax
!!$!      do i = 0, imax-1
!!$!    xyz_DSOSeaIceTempDtPhyUpdt = &
!!$!      & ( xyz_SOSeaIceTempA2 - xyz_SOSeaIceTemp ) / ( 2.0_DP * DelTime )
!!$
!!$
!!$    ! Update sea ice temperature
!!$    !
!!$    xyz_SOSeaIceTemp = xyz_SOSeaIceTemp + xyz_DSOSeaIceTempDtPhyUpdt * DelTime
!!$
!!$
!!$
!!$  end subroutine BK_SOSIUtilsAddPhysics

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

  subroutine SOSIUtilsInit

    !
    ! sosi_utils モジュールの初期化を行います. 
    ! NAMELIST#sosi_utils_nml の読み込みはこの手続きで行われます. 
    !
    ! "sosi_utils" module is initialized. 
    ! "NAMELIST#sosi_utils_nml" is loaded in this procedure. 
    !

    ! モジュール引用 ; USE statements
    !

    ! NAMELIST ファイル入力に関するユーティリティ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid

    ! ファイル入出力補助
    ! File I/O support
    !
    use dc_iounit, only: FileOpen

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: STDOUT ! 標準出力の装置番号. Unit number of standard output

    ! 文字列操作
    ! Character handling
    !
    use dc_string, only: StoA


    ! 宣言文 ; Declaration statements
    !



    ! 作業変数
    ! Work variables
    !
    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号. 
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT. 
                              ! IOSTAT of NAMELIST read

    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /sosi_utils_nml/ &
      & SOSeaIceMassNegativeThreshold
!!$      & SOMass,         &          ! Slab ocean heat capacity (J m-2 K-1)
!!$      & NumMaxItr,      &          ! Number of interation
!!$      & TempItrCrit,    &
!!$      & FlagSublimation
          !
          ! デフォルト値については初期化手続 "sosi_utils#SOSOUtilsInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "sosi_utils#SOSIUtilsInit" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    if ( sosi_utils_inited ) return


    ! デフォルト値の設定
    ! Default values settings
    !
    SOSeaIceMassNegativeThreshold = -1.0e-10_DP


    ! NAMELIST の読み込み
    ! NAMELIST is input
    !
    if ( trim(namelist_filename) /= '' ) then
      call FileOpen( unit_nml, &          ! (out)
        & namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml,                    &  ! (in)
        & nml = sosi_utils_nml,          &  ! (out)
        & iostat = iostat_nml )             ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
    end if



    ! Initialization of modules used in this model
    !


    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, '  SOSeaIceMassNegativeThreshold = %f', &
      & d = (/ SOSeaIceMassNegativeThreshold /) )
!!$    call MessageNotify( 'M', module_name, '  SOMass          = %f', d = (/ SOMass /) )
!!$    call MessageNotify( 'M', module_name, '  SOHeatCapacity  = %f', d = (/ SOHeatCapacity /) )
!!$    call MessageNotify( 'M', module_name, '  NumMaxItr       = %d', i = (/ NumMaxItr /) )
!!$    call MessageNotify( 'M', module_name, '  TempItrCrit     = %f', d = (/ TempItrCrit /) )
!!$    call MessageNotify( 'M', module_name, '  FlagSublimation = %b', l = (/ FlagSublimation /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )


    sosi_utils_inited = .true.

  end subroutine SOSIUtilsInit

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

end module sosi_utils