!= 
!
!= Slab ocean sea ice horizontal transport
!
! Authors::   Yoshiyuki O. TAKAHASHI
! Version::   $Id: gridset.F90,v 1.4 2025/09/17 22:00:00 takepiro Exp $ 
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2013-2025. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module sosi_dynamics
  !
  != 
  !
  != Slab sea ice horizontal transport
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! Horizontal transport of slab sea ice (sea ice on slab ocean) is calculated
  ! based on diffusion.
  !
  !== Procedures List
  !
!!$  ! SLTTMain     :: 移流計算
!!$  ! SLTTInit     :: 初期化
!!$  ! SLTTTest     :: 移流テスト用
!!$  ! ---------------------     :: ------------
!!$  ! SLTTMain     :: Main subroutine for SLTT
!!$  ! SLTTInit     :: Initialization for SLTT
!!$  ! SLTTTest     :: Generate velocity for SLTT Test 
  !
  !== NAMELIST
  !
  ! NAMELIST#
  !
  !== References
  !
  !
  ! モジュール引用 ; USE statements
  !
  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP,  & ! 倍精度実数型. Double precision.
    &                 TOKEN  ! キーワード.   Keywords. 

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

  !
  ! MPI
  !
  use mpi_wrapper, only : MPIWrapperFindMaxVal

  ! 格子点設定
  ! 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
    &                ksimax  ! 海氷の鉛直層数.
                             ! Number of sea ice vertical level

  ! 組成に関わる配列の設定
  ! Settings of array for atmospheric composition
  !
  use composition, only:                              &
    &                    ncmax
                             ! 成分の数
                             ! Number of composition

  ! 質量の補正
  ! Mass fixer
  !
!!$  use mass_fixer, only: &
!!$    & MassFixerBC02, MassFixerBC02Layer, MassFixerBC02Column, &
!!$    & MassFixer, MassFixerR95, MassFixerWO94, MassFixerColumn!, MassFixerLayer


  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedure
  !
  public :: SOSIDynamics
  public :: SOSIDynamicsInit


  ! 公開変数
  ! Public variables
  !

  ! 非公開変数
  ! Private variables
  !
  logical, save      :: FlagSlabOcean
  ! flag for use of slab ocean

  real(DP)    , save              :: SOSeaIceDiffCoef

  real(DP)    , save, allocatable :: y_CosLat(:)
                              ! $\cos\varphai$

  real(DP)    , save, allocatable :: x_LonS   (:)
                              ! $\lambda_S$ 南半球の経度。
                              ! longitude in SH.
  real(DP)    , save, allocatable :: x_SinLonS(:)
                              ! $\sin\lambda_S$
  real(DP)    , save, allocatable :: x_CosLonS(:)
                              ! $\cos\lambda_S$
  real(DP)    , save, allocatable :: y_LatS   (:)
                              ! $\varphi_S$ 南半球の緯度。
                              ! latitude in SH.
  real(DP)    , save, allocatable :: y_SinLatS(:)
                              ! $\sin\varphai_S$
  real(DP)    , save, allocatable :: y_CosLatS(:)
                              ! $\cos\varphai_S$
  real(DP)    , save, allocatable :: x_ExtLonS(:)
                              ! $ x_LonSの拡張配列。
                              !Extended array of x_LonS.
  real(DP)    , save, allocatable :: y_ExtLatS(:)
                              ! $ x_LatSの拡張配列。
                              !Extended array of x_LatS.
  real(DP)    , save, allocatable :: y_ExtCosLatS(:)
                              ! $ y_CosLatS の拡張配列。
                              !Extended array of y_CosLatS.

  real(DP)    , save, allocatable :: x_LonN   (:)
                              ! $\lambda_N$ 北半球の経度。
                              ! longitude in NH.
  real(DP)    , save, allocatable :: x_SinLonN(:)
                              ! $\sin\lambda_N$
  real(DP)    , save, allocatable :: x_CosLonN(:)
                              ! $\cos\lambda_N$
  real(DP)    , save, allocatable :: y_LatN   (:)
                              ! $\varphi_N$ 北半球の緯度。
                              ! latitude in NH.
  real(DP)    , save, allocatable :: y_SinLatN(:)
                              ! $\sin\varphai_N$
  real(DP)    , save, allocatable :: y_CosLatN(:)
                              ! $\cos\varphai_N$
  real(DP)    , save, allocatable :: x_ExtLonN(:)
                              ! $ x_LonNの拡張配列。
                              !Extended array of x_LonN.
  real(DP)    , save, allocatable :: y_ExtLatN(:)
                              ! $ x_LatNの拡張配列。
                              !Extended array of x_LatN.
  real(DP)    , save, allocatable :: y_ExtCosLatN(:)
                              ! $ y_CosLatNの拡張配列。
                              !Extended array of y_CosLatN.

  real(DP)    , save, allocatable :: p_LonS   (:)
  real(DP)    , save, allocatable :: q_LatS   (:)
  real(DP)    , save, allocatable :: q_CosLatS(:)
  real(DP)    , save, allocatable :: p_LonN   (:)
  real(DP)    , save, allocatable :: q_LatN   (:)
  real(DP)    , save, allocatable :: q_CosLatN(:)


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


  character(*), parameter:: module_name = 'sosi_dynamics'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: sltt.F90,v 1.8 2025/09/17 22:00:00 takepiro Exp $'
                              ! モジュールのバージョン
                              ! Module version


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

contains

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

  subroutine SOSIDynamics(                                &
    & xy_SurfType,                                        & !(in   )
    & xy_SurfTemp, xy_SOSeaIceMass, xyz_SOSeaIceTemp,     & !(inout)
    & xy_DSOSeaIceMassDtPhyTop, xy_DSOSeaIceMassDtPhyBot, & !(in   )
    & xyz_DSOSeaIceTempDtPhy                              & !(in   )
    & )
    ! 
    ! Calculates slab sea ice horizontal transports by diffusion

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

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

    ! 物理・数学定数設定
    ! Physical and mathematical constants settings
    !
    use constants0, only: &
      & WaterHeatCap

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & SOMass
                              ! Slab ocean mass

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

    !
    ! Slab ocean sea ice utility module
    !
    use sosi_utils, only :             &
      & SOSIUtilsChkSOSeaIce,          &
      & SOSIUtilsSetSOSeaIceLevels,    &
      & SOSeaIceMassNegativeThreshold, &
      & SOSIUtilsAddPhysics


    ! 宣言文 ; Declaration statements
    !
    integer , intent(in ) :: xy_SurfType       (0:imax-1, 1:jmax)
                              ! 土地利用.
                              ! Surface index
    real(DP), intent(inout) :: xy_SurfTemp        (0:imax-1, 1:jmax)
    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)
                              !
                              ! Slab sea ice mass at next time step
    real(DP), intent(in   ) :: xyz_DSOSeaIceTempDtPhy(0:imax-1, 1:jmax, 1:ksimax)


    ! 作業変数
    ! Work variables
    !
    real(DP) :: xy_SurfTempSave     (0:imax-1, 1:jmax)
    real(DP) :: xy_SOSeaIceMassSave (0:imax-1, 1:jmax)
    real(DP) :: xyz_SOSeaIceTempSave(0:imax-1, 1:jmax, ksimax)

    real(DP) :: xy_SOSeaIceThickness(0:imax-1, 1:jmax)

    real(DP) :: xyz_SOSeaIceThickness(0:imax-1, 1:jmax, 1:ksimax)
                 !
                 ! 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)

    logical  :: xy_FlagSlabOcean(0:imax-1, 1:jmax)

    real(DP) :: xyz_DSOSeaIceThicknessDt(0:imax-1, 1:jmax, 1:ksimax)
    real(DP) :: xyz_DSOSeaIceTempDt     (0:imax-1, 1:jmax, 1:ksimax)
    real(DP) :: xy_DSOTempDt            (0:imax-1, 1:jmax)

!!$    real(DP) :: xy_DSOSeaIceMassDt(0:imax-1, 1:jmax)
!!$                              !
!!$                              ! Slab sea ice mass tendency

    real(DP) :: xy_SOTemp               (0:imax-1, 1:jmax)

    real(DP) :: xyz_SOSeaIceThicknessAdv(0:imax-1, 1:jmax, 1:ksimax)
    real(DP) :: xyz_SOSeaIceTempAdv     (0:imax-1, 1:jmax, 1:ksimax)
    real(DP) :: xy_SOTempAdv            (0:imax-1, 1:jmax)

!!$    real(DP) :: xy_TempSlabOcean (0:imax-1, 1:jmax)
!!$                              !
!!$                              ! Slab ocean temperature
!!$    real(DP) :: xy_TempSlabSeaIce(0:imax-1, 1:jmax)
!!$                              !
!!$                              ! Slab sea ice temperature

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

    real(DP) :: SOSIMass
    real(DP) :: SOSIMass1L
    real(DP) :: DelSOSIMass

!!$    real(DP) :: SurfTempTent
    real(DP) :: SOTempTent
    real(DP) :: SOTempTent1st

    logical, parameter :: FlagSOSIAdv = .false.

    logical  :: FlagCalc

    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_dynamics_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    ! add physics tendency
    call SOSIUtilsAddPhysics(                         &
      & xy_SOSeaIceMass, xyz_SOSeaIceTemp,            & !(inout)
      & xy_DSOSeaIceMassDtPhyTop, xy_DSOSeaIceMassDtPhyBot, & !(in   )
      & xyz_DSOSeaIceTempDtPhy                           & !(in   )
      & )


!!$    xy_SOSeaIceMassA = xy_SOSeaIceMassB &
!!$      & + xy_DSOSeaIceMassDtPhy * ( 2.0_DP * DelTime )
!!$
!!$
!!$    !
!!$    ! Calcuate sea ice thickness
!!$    !
!!$    xy_SeaIceThickness = xy_SOSeaIceMassB / SeaIceDen
!!$    !
!!$    ! Set slab ocean sea ice levels
!!$    !
!!$    call SOSIUtilsSetSOSeaIceLevels(                     &
!!$      & xy_SeaIceThickness,                                       & ! (in   )
!!$      & xy_SOSILocalKMax, xyr_SOSILocalDepth, xyz_SOSILocalDepth  & ! (out)
!!$      & )
!!$
!!$    !
!!$    ! Calcuate sea ice thickness
!!$    !
!!$    xy_SeaIceThicknessA = xy_SOSeaIceMassA / SeaIceDen
!!$    !
!!$    ! Set slab ocean sea ice levels
!!$    !
!!$    call SOSIUtilsSetSOSeaIceLevels(                     &
!!$      & xy_SeaIceThicknessA,                                       & ! (in   )
!!$      & xy_SOSILocalKMaxA, xyr_SOSILocalDepthA, xyz_SOSILocalDepthA  & ! (out)
!!$      & )
!!$
!!$
!!$    ! 海氷温度時間積分
!!$    ! Time integration of sea ice temperature
!!$    !
!!$    xyz_SOSeaIceTemp = xyz_SOSeaIceTemp + xyz_DSOSeaIceTempDtPhy * DelTime
!!$
!!$
!!$    ! Adjust levels
!!$    !
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        if ( xy_SOSILocalKMaxA(i,j) > xy_SOSILocalKMax(i,j) ) then
!!$          ! sea ice thickness increases
!!$          do k = 1, xy_SOSILocalKMax(i,j)
!!$            xyz_SOSeaIceTemp(i,j,k) = xyz_SOSeaIceTemp(i,j,k) &
!!$              & + xyz_DSOSeaIceTempDtPhy(i,j,k) * DelTime
!!$          end do
!!$          do k = xy_SOSILocalKMax(i,j)+1, xy_SOSILocalKMaxA(i,j)
!!$            kk = xy_SOSILocalKMax(i,j)
!!$            xyz_SOSeaIceTemp(i,j,k) = xyz_SOSeaIceTemp(i,j,kk)
!!$          end do
!!$        else if ( xy_SOSILocalKMaxA(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



    xy_FlagSlabOcean = .false.
    FlagCalc         = .false.
    if ( FlagSlabOcean ) then
      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_SurfType(i,j) <= 0 ) then
            ! slab ocean
            xy_FlagSlabOcean(i,j) = .true.
            FlagCalc              = .true.
          end if
        end do
      end do
    end if

    if ( SOSeaIceDiffCoef <= 0.0_DP ) then
      FlagCalc = .false.
!!$    else
!!$      call MessageNotify( 'E', module_name, &
!!$        & '  Now, SOSeaIceDiffCoef has to be zero.' )
    end if

    if ( .not. FlagCalc ) then
      return
    end if



    ! Save values
    !
    xy_SurfTempSave      = xy_SurfTemp
    xy_SOSeaIceMassSave  = xy_SOSeaIceMass
    xyz_SOSeaIceTempSave = xyz_SOSeaIceTemp


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

    ! Slab ocean temperature
    !
    do j = 1, jmax
      do i = 0, imax-1
!!$        if ( xy_SOSILocalKMax(i,j) <= 0 ) then
        if ( xy_SOSeaIceMass(i,j) < SOSeaIceThresholdMass ) then
          xy_SOTemp(i,j) = xy_SurfTemp(i,j)
        else
          xy_SOTemp(i,j) = TempBelowSeaIce
        end if
      end do
    end do

    do j = 1, jmax
      do i = 0, imax-1
        do k = 1, xy_SOSILocalKMax(i,j)
          xyz_SOSeaIceThicknessAdv(i,j,k) = &
            & xyr_SOSILocalDepth(i,j,k-1) - xyr_SOSILocalDepth(i,j,k)
          xyz_SOSeaIceTempAdv     (i,j,k) = xyz_SOSeaIceTemp(i,j,k)
        end do
        do k = xy_SOSILocalKMax(i,j)+1, ksimax
          xyz_SOSeaIceThicknessAdv(i,j,k) = 0.0_DP
          xyz_SOSeaIceTempAdv     (i,j,k) = TempCondWater
!!$          xyz_SOSeaIceTempAdv     (i,j,k) = TempBelowSeaIce
        end do
      end do
    end do
    !
    xy_SOTempAdv = xy_SOTemp


!!$    call SOSIUtilsSetMissingValue( &
!!$      & xy_SOSeaIceMass,                 & !(in )
!!$      & xyz_SOSeaIceTemp,                & !(inout)
!!$      & SOSeaIceValue                    & !(in ) optional
!!$      & )

    call SOSIHorTransportDiff(                          &
      & xy_FlagSlabOcean,                               & ! (in)
      & xy_SOSeaIceMass,                                & ! (in)
      & xy_SOTempAdv,                                   & ! (in)
      & xyz_SOSeaIceThicknessAdv, xyz_SOSeaIceTempAdv,  & ! (in)
      & xy_DSOTempDt,                                   & ! (out)
      & xyz_DSOSeaIceThicknessDt, xyz_DSOSeaIceTempDt   & ! (out)
      & )


    xyz_SOSeaIceThickness = xyz_SOSeaIceThicknessAdv
    xyz_SOSeaIceTemp      = xyz_SOSeaIceTempAdv
    xy_SOTemp             = xy_SOTempAdv
    !
    if ( FlagSOSIAdv ) then
      xyz_SOSeaIceThickness = xyz_SOSeaIceThickness &
        & + xyz_DSOSeaIceThicknessDt * DelTime
      xyz_SOSeaIceTemp      = xyz_SOSeaIceTemp      &
        & + xyz_DSOSeaIceTempDt      * DelTime
    else
      xy_SOTemp = xy_SOTemp &
        & + xy_DSOTempDt * DelTime
    end if



    ! Adjustment
    !
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        if ( xy_FlagSlabOcean(i,j) ) then
!!$          do k = xy_SOSILocalKMax(i,j), 1, -1
!!$            if ( xyz_SOSeaIceThickness(i,j,k) > 0.0_DP ) then
!!$              if ( xy_SurfTemp(i,j) > TempCondWater ) then
!!$                ! Check whether all sea ice melt
!!$                SOSIMass = SeaIceDen * xyz_SOSeaIceThickness(i,j,k)
!!$                SurfTempTent = &
!!$                  &   (   WaterHeatCap * ( SOMass - SOSIMass ) * xy_SurfTemp(i,j) &
!!$                  &     + SeaIceHeatCap * SOSIMass * xyz_SOSeaIceTemp(i,j,k)      &
!!$                  &     - LatentHeatFusion * SOSIMass )                           &
!!$                  & / ( WaterHeatCap * SOMass )
!!$
!!$                if ( SurfTempTent >= TempCondWater ) then
!!$                  ! Case in which all sea ice melt
!!$                  xy_SurfTemp          (i,j)   = SurfTempTent
!!$                  xyz_SOSeaIceThickness(i,j,k) = 0.0_DP
!!$                  xyz_SOSeaIceTemp     (i,j,k) = TempCondWater
!!$                else
!!$                  ! Case in which a part of sea ice melt
!!$                  DelSOSIMass =                                  &
!!$                    & - WaterHeatCap * ( SOMass - SOSIMass )     &
!!$                    &   * ( xy_SurfTemp(i,j) - TempBelowSeaIce ) &
!!$                    &   / (                                                   &
!!$                    &         LatentHeatFusion                                &
!!$                    &       + SeaIceHeatCap                                   &
!!$                    &         * ( TempBelowSeaIce - xyz_SOSeaIceTemp(i,j,k) ) &
!!$                    &      )
!!$
!!$                  xy_SurfTemp(i,j) = TempCondWater
!!$                  xyz_SOSeaIceThickness(i,j,k) =     &
!!$                    &   xyz_SOSeaIceThickness(i,j,k) &
!!$                    & + DelSOSIMass / SeaIceDen
!!$                end if
!!$              end if
!!$            end if
!!$          end do
!!$        end if
!!$      end do
!!$    end do
!!$    !
!!$    xy_SOSeaIceMass = 0.0_DP
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        do k = 1, xy_SOSILocalKMax(i,j)
!!$          xy_SOSeaIceMass(i,j) = xy_SOSeaIceMass(i,j) &
!!$            & + SeaIceDen * xyz_SOSeaIceThickness(i,j,k)
!!$        end do
!!$      end do
!!$    end do

    do j = 1, jmax
      do i = 0, imax-1
        do k = 1, xy_SOSILocalKMax(i,j)
          xyz_SOSIMassEachLayer(i,j,k) = SeaIceDen * xyz_SOSeaIceThickness(i,j,k)
        end do
        do k = xy_SOSILocalKMax(i,j)+1, ksimax
          xyz_SOSIMassEachLayer(i,j,k) = 0.0_DP
        end do
      end do
    end do

    xy_SOSeaIceMass = 0.0_DP
    do j = 1, jmax
      do i = 0, imax-1
        do k = 1, xy_SOSILocalKMax(i,j)
          xy_SOSeaIceMass(i,j) = xy_SOSeaIceMass(i,j) &
            & + xyz_SOSIMassEachLayer(i,j,k)
        end do
      end do
    end do

    do j = 1, jmax
      do i = 0, imax-1
        if ( xy_FlagSlabOcean(i,j) ) then
          do k = xy_SOSILocalKMax(i,j), 1, -1
            if ( xy_SOTemp(i,j) > TempBelowSeaIce ) then
              ! Check whether all sea ice melt
              SOSIMass    = xy_SOSeaIceMass(i,j)
              SOSIMass1L  = xyz_SOSIMassEachLayer(i,j,k)
              DelSOSIMass = - SOSIMass1L

!!$              SOTempTent =                                                    &
!!$                &   (   WaterHeatCap * ( SOMass - SOSIMass ) * xy_SOTemp(i,j) &
!!$                &     + SeaIceHeatCap * SOSIMass1L * xyz_SOSeaIceTemp(i,j,k)  &
!!$                &     + LatentHeatFusion * DelSOSIMass )                      &
!!$                & / ( WaterHeatCap * ( ( SOMass - SOSIMass ) - DelSOSIMass ) )

              SOTempTent1st = TempBelowSeaIce
              SOTempTent =                                                    &
                &   (   WaterHeatCap * ( SOMass - SOSIMass ) * xy_SOTemp(i,j) &
                &     - WaterHeatCap * DelSOSIMass * SOTempTent1st            &
                &     + SeaIceHeatCap * DelSOSIMass                           &
                &         * ( SOTempTent1st - xyz_SOSeaIceTemp(i,j,k) )       &
                &     + LatentHeatFusionBelowSeaice * DelSOSIMass )           &
                & / ( WaterHeatCap * ( ( SOMass - SOSIMass ) - DelSOSIMass ) )

              if ( SOTempTent >= TempBelowSeaIce ) then
                ! Case in which all sea ice melt
                xy_SOTemp            (i,j)   = SOTempTent
                xyz_SOSeaIceTemp     (i,j,k) = TempBelowSeaIce
              else
                ! Case in which a part of sea ice melt
                SOTempTent = TempBelowSeaIce

!!$                DelSOSIMass =                                  &
!!$                  &   ( &
!!$                  &     - WaterHeatCap * ( SOMass - SOSIMass )     &
!!$                  &         * ( SOTempTent - xy_SOTemp(i,j) )      &
!!$                  &     - SeaIceHeatCap * SOSIMass1L               &
!!$                  &         * ( xyz_SOSeaIceTemp(i,j,k) - xyz_SOSeaIceTemp(i,j,k) )  &
!!$                  &   ) &
!!$                  & / (                                                 &
!!$                  &     - WaterHeatCap * SOTempTent                     &
!!$                  &     + SeaIceHeatCap * xyz_SOSeaIceTemp(i,j,k)       &
!!$                  &     - LatentHeatFusion                              &
!!$                  &   )

                SOTempTent1st = TempBelowSeaIce
                DelSOSIMass =                                  &
                  &   ( &
                  &       WaterHeatCap * ( SOMass - SOSIMass ) &
                  &         * ( SOTempTent - xy_SOTemp(i,j) )  &
                  &   ) &
                  & / (                                                               &
                  &       WaterHeatCap * ( SOTempTent - SOTempTent1st )               &
                  &     + SeaIceHeatCap * ( SOTempTent1st - xyz_SOSeaIceTemp(i,j,k) ) &
                  &     + LatentHeatFusionBelowSeaIce                                 &
                  &   )

                xy_SOTemp(i,j) = SOTempTent
              end if
              xyz_SOSIMassEachLayer(i,j,k) = &
                &   xyz_SOSIMassEachLayer(i,j,k) + DelSOSIMass
              xy_SOSeaIceMass(i,j) = xy_SOSeaIceMass(i,j) + DelSOSIMass
            end if
          end do
        end if
      end do
    end do

    do j = 1, jmax
      do i = 0, imax-1
        if ( xy_SOSeaIceMass(i,j) > SOSeaIceThresholdMass ) then
          ! in case with sea ice remained
          xy_SurfTemp(i,j) = xy_SurfTemp(i,j)
        else if ( xy_SOSeaIceMassSave(i,j) > SOSeaIceThresholdMass ) then
          ! in case with sea ice was present but melts completely
          xy_SurfTemp(i,j) = xy_SOTemp(i,j)
        else
          ! in case with sea ice was/is not present
!!$          xy_SurfTemp(i,j) = xy_SurfTemp(i,j)
          xy_SurfTemp(i,j) = xy_SOTemp(i,j)
        end if
      end do
    end do


    do j = 1, jmax
      do i = 0, imax-1
        if ( xy_FlagSlabOcean(i,j) ) then
          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 diffusion, %f, and this is set to zero.', &
                & d = (/ xy_SOSeaIceMass(i,j) /) )
            end if
            xy_SOSeaIceMass(i,j) = 0.0_DP
          end if
        end if
      end do
    end do


    ! Check
    !
    xy_SOSeaIceThickness = xy_SOSeaIceMass / SeaIceDen
    !
    call SOSIUtilsChkSOSeaIce(                &
      & xy_SOSeaIceThickness, xyz_SOSeaIceTemp, & ! (in)
      & "SOSIDynamics"                        & ! (in)
      & )


  end subroutine SOSIDynamics

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

  subroutine SOSIHorTransportFFSL(                    &
    & xy_FlagSlabOcean,                               & ! (in)
    & xy_SOSeaIceMass,                                & ! (in)
    & xy_SurfTemp,                                    & ! (in)
    & xyz_SOSeaIceThickness, xyz_SOSeaIceTemp,        & ! (in)
    & xy_DSurfTempDt,                                 & ! (out)
    & xyz_DSOSeaIceThicknessDt, xyz_DSOSeaIceTempDt   & ! (out)
    & )
    ! 
    ! Calculates slab sea ice transport by diffusion

    use axesset    , only : x_Lon
                              ! $\lambda, \varphai$ lon and lat
    use sltt_const , only : dtjw, iexmin, iexmax, jexmin, jexmax
    use sosi_dynamics_extarr, only : SLTTExtArrExt2
                              ! 配列拡張ルーチン
                              ! Expansion of arrays
    ! 雪と海氷の定数の設定
    ! Setting constants of snow and sea ice
    !
    use constants_snowseaice, only: &
      & SeaIceDen


    logical , intent(in ) :: xy_FlagSlabOcean(0:imax-1, 1:jmax)
    real(DP), intent(in ) :: xy_SOSeaIceMass      (0:imax-1, 1:jmax)
    real(DP), intent(in ) :: xy_SurfTemp          (0:imax-1, 1:jmax)
    real(DP), intent(in ) :: xyz_SOSeaIceThickness(0:imax-1, 1:jmax, 1:ksimax)
    real(DP), intent(in ) :: xyz_SOSeaIceTemp     (0:imax-1, 1:jmax, 1:ksimax)
    real(DP), intent(out) :: xy_DSurfTempDt          (0:imax-1, 1:jmax)
    real(DP), intent(out) :: xyz_DSOSeaIceThicknessDt(0:imax-1, 1:jmax, 1:ksimax)
    real(DP), intent(out) :: xyz_DSOSeaIceTempDt     (0:imax-1, 1:jmax, 1:ksimax)


    !
    ! local variables
    !
    real(DP) :: xy_SurfTempA          (0:imax-1, 1:jmax)
    real(DP) :: xyz_SOSeaIceThicknessA(0:imax-1, 1:jmax, 1:ksimax)
    real(DP) :: xyz_SOSeaIceTempA     (0:imax-1, 1:jmax, 1:ksimax)

    real(DP) :: xy_SOSeaIceMassA      (0:imax-1, 1:jmax)

    real(DP) :: xyza_TMP4DArray(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)

    real(DP) :: xyza_ExtTMP4DArrayS(iexmin:iexmax, jexmin:jexmax, 1:kmax, 1:ncmax)
                              ! 
                              ! Extended 4D array (SH)
    real(DP) :: xyza_ExtTMP4DArrayN(iexmin:iexmax, jexmin:jexmax, 1:kmax, 1:ncmax)
                              ! 
                              ! Extended 4D array (NH)

    real(DP) :: xyz_ExtSOSeaIceThicknessS(iexmin:iexmax, jexmin:jexmax, 1:ksimax)
    real(DP) :: xyz_ExtSOSeaIceTempS     (iexmin:iexmax, jexmin:jexmax, 1:ksimax)
    real(DP) :: xyz_ExtSurfTempS         (iexmin:iexmax, jexmin:jexmax, 1:ksimax)
    real(DP) :: xy_ExtSOSeaIceMassS      (iexmin:iexmax, jexmin:jexmax)
                              ! 
                              ! Extended 4D array (SH)
    real(DP) :: xyz_ExtSOSeaIceThicknessN(iexmin:iexmax, jexmin:jexmax, 1:ksimax)
    real(DP) :: xyz_ExtSOSeaIceTempN     (iexmin:iexmax, jexmin:jexmax, 1:ksimax)
    real(DP) :: xyz_ExtSurfTempN         (iexmin:iexmax, jexmin:jexmax, 1:ksimax)
    real(DP) :: xy_ExtSOSeaIceMassN      (iexmin:iexmax, jexmin:jexmax)
                              ! 
                              ! Extended 4D array (NH)

    logical  :: xy_ExtFlagSlabOceanS(iexmin:iexmax, jexmin:jexmax)
                              ! 
                              ! Extended 4D array (SH)
    logical  :: xy_ExtFlagSlabOceanN(iexmin:iexmax, jexmin:jexmax)
                              ! 
                              ! Extended 4D array (NH)

    real(DP) :: xyz_DSOSeaIceThicknessDtS(0:imax-1, 1:jmax/2, 1:ksimax)
    real(DP) :: xyz_DSOSeaIceTempDtS     (0:imax-1, 1:jmax/2, 1:ksimax)
    real(DP) :: xya_DSurfTempDtS         (0:imax-1, 1:jmax/2, 1:ksimax)
    real(DP) :: xyz_DSOSeaIceThicknessDtN(0:imax-1, 1:jmax/2, 1:ksimax)
    real(DP) :: xyz_DSOSeaIceTempDtN     (0:imax-1, 1:jmax/2, 1:ksimax)
    real(DP) :: xya_DSurfTempDtN         (0:imax-1, 1:jmax/2, 1:ksimax)



    real(DP) :: PM
    ! 配列拡張する際、極ごえ後に符号が変わる場合は -1.0を与える。
    ! そうでない場合は1.0を与える。
    ! Sign change flag for array extension; -1.0 for sign change
    ! over the pole, 1.0 for no sign change



    integer:: i               ! 東西方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in zonal direction
    integer:: j               ! 南北方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in meridional direction
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
    integer:: n               ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituents
    integer:: kk


    ! 実行文 ; Executable statement
    !

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


    !
    ! Longitudinal diffusion
    !

#if defined(AXISYMMETRY) || defined(AXISYMMETRY_SJPACK)
    xy_SurfTempA           = xy_SurfTemp
    xyz_SOSeaIceThicknessA = xyz_SOSeaIceThickness
    xyz_SOSeaIceTempA      = xyz_SOSeaIceTemp
#else
    call SOSIDiffusionX(                               &
      & x_Lon(1)-x_Lon(0), y_CosLat, xy_FlagSlabOcean, & ! (in)
      & xy_SOSeaIceMass,                                             & ! (in)
      & xy_SurfTemp , xyz_SOSeaIceThickness , xyz_SOSeaIceTemp,      & ! (in)
      & xy_SurfTempA, xyz_SOSeaIceThicknessA, xyz_SOSeaIceTempA      & ! (out)
      & )
#endif

    xy_SOSeaIceMassA = 0.0_DP
    do k = 1, ksimax
      xy_SOSeaIceMassA = xy_SOSeaIceMassA &
        & + SeaIceDen * xyz_SOSeaIceThicknessA(:,:,k)
    end do


    !
    ! Latitudinal diffusion
    !

    ! 配列の分割と拡張
    ! Division and extension of arrays
    !

    if ( ksimax > kmax ) then
      call MessageNotify( 'E', module_name, 'ksimax > kmax.' )
    end if
    if ( kmax < 3 ) then
      call MessageNotify( 'E', module_name, 'kmax < 3.' )
    end if
    if ( ncmax < 3 ) then
      call MessageNotify( 'E', module_name, 'ncmax < 3.' )
    end if


    n = 1
    do k = 1, ksimax
      xyza_TMP4DArray(:,:,k,n) = xyz_SOSeaIceThicknessA(:,:,k)
    end do
    do k = ksimax+1, kmax
      xyza_TMP4DArray(:,:,k,n) = 0.0_DP
    end do
    !
    n = 2
    do k = 1, ksimax
      xyza_TMP4DArray(:,:,k,n) = xyz_SOSeaIceTempA(:,:,k)
    end do
    do k = ksimax+1, kmax
      xyza_TMP4DArray(:,:,k,n) = 0.0_DP
    end do
    !
    n = 3
    k = 1
    do j = 1, jmax
      do i = 0, imax-1
        if ( xy_FlagSlabOcean(i,j) ) then
          xyza_TMP4DArray(i,j,k,n) = 1.0_DP
        else
          xyza_TMP4DArray(i,j,k,n) = 0.0_DP
        end if
      end do
    end do
    k = 2
    xyza_TMP4DArray(:,:,k,n) = xy_SurfTempA
    k = 3
    xyza_TMP4DArray(:,:,k,n) = xy_SOSeaIceMassA
    do k = 3+1, kmax
      xyza_TMP4DArray(:,:,k,n) = 0.0_DP
    end do
    !
    do n = 3+1, ncmax
      xyza_TMP4DArray(:,:,:,n) = 0.0_DP
    end do

    PM = 1.0_DP
    call SLTTExtArrExt2(                            &
      & x_SinLonS, x_CosLonS, x_SinLonN, x_CosLonN, & ! (in)
      & xyza_TMP4DArray, PM,                        & ! (in)
      & xyza_ExtTMP4DArrayS, xyza_ExtTMP4DArrayN    & ! (out)
      & )

    n = 1
    do k = 1, ksimax
      xyz_ExtSOSeaIceThicknessS(:,:,k) = xyza_ExtTMP4DArrayS(:,:,k,n)
      xyz_ExtSOSeaIceThicknessN(:,:,k) = xyza_ExtTMP4DArrayN(:,:,k,n)
    end do
    n = 2
    do k = 1, ksimax
      xyz_ExtSOSeaIceTempS     (:,:,k) = xyza_ExtTMP4DArrayS(:,:,k,n)
      xyz_ExtSOSeaIceTempN     (:,:,k) = xyza_ExtTMP4DArrayN(:,:,k,n)
    end do
    n = 3
    k = 1
    do j = jexmin, jexmax
      do i = iexmin, iexmax
        if ( xyza_ExtTMP4DArrayS(i,j,k,n) == 1.0_DP ) then
          xy_ExtFlagSlabOceanS(i,j) = .true.
        else
          xy_ExtFlagSlabOceanS(i,j) = .false.
        end if
        if ( xyza_ExtTMP4DArrayN(i,j,k,n) == 1.0_DP ) then
          xy_ExtFlagSlabOceanN(i,j) = .true.
        else
          xy_ExtFlagSlabOceanN(i,j) = .false.
        end if
      end do
    end do
    k = 2
    do kk = 1, ksimax
      xyz_ExtSurfTempS(:,:,kk) = xyza_ExtTMP4DArrayS(:,:,k,n)
      xyz_ExtSurfTempN(:,:,kk) = xyza_ExtTMP4DArrayN(:,:,k,n)
    end do
    k = 3
    do kk = 1, ksimax
      xy_ExtSOSeaIceMassS = xyza_ExtTMP4DArrayS(:,:,k,n)
      xy_ExtSOSeaIceMassN = xyza_ExtTMP4DArrayN(:,:,k,n)
    end do


!!$    call SSIDiffusion(                             &
!!$      & x_ExtLonS, y_ExtLatS, y_ExtCosLatS,        & ! (in)
!!$      & p_LonS, q_LatS, q_CosLatS,                 & ! (in)
!!$      & xy_ExtFlagSlabOceanS, xy_ExtSOSeaIceMassS, & ! (in)
!!$      & xy_DSOSeaIceMassDtS                        & ! (out)
!!$      & )
!!$    call SSIDiffusion(                             &
!!$      & x_ExtLonN, y_ExtLatN, y_ExtCosLatN,        & ! (in)
!!$      & p_LonN, q_LatN, q_CosLatN,                 & ! (in)
!!$      & xy_ExtFlagSlabOceanN, xy_ExtSOSeaIceMassN, & ! (in)
!!$      & xy_DSOSeaIceMassDtN                        & ! (out)
!!$      & )

    call SOSIDiffusionY(                                 &
      & y_ExtLatS, y_ExtCosLatS,                         & ! (in)
      & q_LatS, q_CosLatS,                               & ! (in)
      & xy_ExtFlagSlabOceanS, xyz_ExtSOSeaIceThicknessS, & ! (in)
      & xyz_DSOSeaIceThicknessDtS                        & ! (out)
      & )
    call SOSIDiffusionY(                                 &
      & y_ExtLatN, y_ExtCosLatN,                         & ! (in)
      & q_LatN, q_CosLatN,                               & ! (in)
      & xy_ExtFlagSlabOceanN, xyz_ExtSOSeaIceThicknessN, & ! (in)
      & xyz_DSOSeaIceThicknessDtN                        & ! (out)
      & )
    !
    call SOSIDiffusionY(                                 &
      & y_ExtLatS, y_ExtCosLatS,                         & ! (in)
      & q_LatS, q_CosLatS,                               & ! (in)
      & xy_ExtFlagSlabOceanS, xyz_ExtSOSeaIceTempS, & ! (in)
      & xyz_DSOSeaIceTempDtS                        & ! (out)
      & )
    call SOSIDiffusionY(                                 &
      & y_ExtLatN, y_ExtCosLatN,                         & ! (in)
      & q_LatN, q_CosLatN,                               & ! (in)
      & xy_ExtFlagSlabOceanN, xyz_ExtSOSeaIceTempN, & ! (in)
      & xyz_DSOSeaIceTempDtN                        & ! (out)
      & )
    !
    call SOSIDiffusionY(                                 &
      & y_ExtLatS, y_ExtCosLatS,                         & ! (in)
      & q_LatS, q_CosLatS,                               & ! (in)
      & xy_ExtFlagSlabOceanS, xyz_ExtSurfTempS,     & ! (in)
      & xya_DSurfTempDtS,                           & ! (out)
      & xy_ExtSOSeaIceMassS                         & ! (in) optional
      & )
    call SOSIDiffusionY(                                 &
      & y_ExtLatN, y_ExtCosLatN,                         & ! (in)
      & q_LatN, q_CosLatN,                               & ! (in)
      & xy_ExtFlagSlabOceanN, xyz_ExtSurfTempN,     & ! (in)
      & xya_DSurfTempDtN,                           & ! (out)
      & xy_ExtSOSeaIceMassN                         & ! (in) optional
      & )

    xyz_DSOSeaIceThicknessDt(:,1:jmax/2     ,:) = xyz_DSOSeaIceThicknessDtS(:,1:jmax/2,:)
    xyz_DSOSeaIceThicknessDt(:,jmax/2+1:jmax,:) = xyz_DSOSeaIceThicknessDtN(:,1:jmax/2,:)
    !
    xyz_DSOSeaIceTempDt(:,1:jmax/2     ,:) = xyz_DSOSeaIceTempDtS(:,1:jmax/2,:)
    xyz_DSOSeaIceTempDt(:,jmax/2+1:jmax,:) = xyz_DSOSeaIceTempDtN(:,1:jmax/2,:)
    !
    xy_DSurfTempDt(:,1:jmax/2     ) = xya_DSurfTempDtS(:,1:jmax/2,1)
    xy_DSurfTempDt(:,jmax/2+1:jmax) = xya_DSurfTempDtN(:,1:jmax/2,1)


!!$    ! sea ice mass at next time step is calculated temporarily
!!$    xy_SOSeaIceMassA = xy_SOSeaIceMassA &
!!$      & + xy_DSOSeaIceMassDt * ( 2.0_DP * DelTime )
!!$
!!$    ! tendency is calculated
!!$    xy_DSOSeaIceMassDt = &
!!$      & ( xy_SOSeaIceMassA - xy_SOSeaIceMassB ) / ( 2.0_DP * DelTime )


!!$    py_ExtSOSeaIceMassXFlux = 0.0_DP
!!$    xq_ExtSOSeaIceMassYFlux = 0.0_DP
!!$    do j = jexmin, jexmax-1
!!$      do i = iexmin, iexmax-1
!!$        py_ExtSOSeaIceMassXFlux(i,j) = &
!!$          & - SOSeaIceDiffCoef &
!!$          &   * ( xy_ExtSOSeaIceMassS(i+1,j) - xy_ExtSOSeaIceMassS(i,j) ) &
!!$          &   / ( RPlanet * y_CosLatS(j) * ( x_ExtLonS(i+1) - x_ExtLonS(i) ) )
!!$        xq_ExtSOSeaIceMassYFlux(i,j) = &
!!$          & - SOSeaIceDiffCoef &
!!$          &   * ( xy_ExtSOSeaIceMassS(i,j+1) - xy_ExtSOSeaIceMassS(i,j) ) &
!!$          &   / ( RPlanet * ( y_ExtLatS(j+1) - y_ExtLatS(j) ) )
!!$      end do
!!$    end do
!!$    if ( myrank == nprocs-1 ) then
!!$      j = 0
!!$      do i = iexmin, iexmax-1
!!$        xq_ExtSOSeaIceMassYFlux(i,j) = 0.0_DP
!!$      end do
!!$    end if
!!$
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        xy_DSOSeaIceMassDt(i,j) = &
!!$          & - ( py_ExtSOSeaIceMassXFlux(i,j) - py_ExtSOSeaIceMassXFlux(i-1,j) ) &
!!$          &   / ( RPlanet * y_CosLatS(j) * ( p_Lon(i) - p_Lon(i-1) ) ) &
!!$          & - (   xq_ExtSOSeaIceMassYFlux(i,j  ) * q_CosLatS(j  )   &
!!$          &     - xq_ExtSOSeaIceMassYFlux(i,j-1) * q_CosLatS(j-1) ) &
!!$          &   / ( RPlanet * y_CosLatS(j) * ( q_LatS(j) - q_LatS(j-1) ) )
!!$      end do
!!$    end do




  end subroutine SOSIHorTransportFFSL

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

!!$  subroutine SOSIFFSLX(                          &
!!$    & DelLon, y_CosLat, xy_FlagSlabOcean,        & ! (in)
!!$    & xy_SOSeaIceMass,                           & ! (in)
!!$    & xy_SurfTemp , xyz_SOSeaIceThickness , xyz_SOSeaIceTemp,  & ! (in)
!!$    & xy_SurfTempA, xyz_SOSeaIceThicknessA, xyz_SOSeaIceTempA  & ! (out)
!!$    & )
!!$    ! 
!!$    ! Calculates slab sea ice transport by longitudinal diffusion
!!$
!!$
!!$    !
!!$    !
!!$    use ludecomp_module, only : &
!!$      & ludecomp_prep_simple_many, &
!!$      & ludecomp_solve_simple_many
!!$
!!$    use constants, only: &
!!$      & RPlanet, &
!!$      ! $ a $ [m].
!!$      ! 惑星半径.
!!$      ! Radius of planet
!!$      & SOMass
!!$                              ! Slab ocean mass
!!$
!!$    real(DP), intent(in ) :: DelLon
!!$    real(DP), intent(in ) :: y_CosLat(1:jmax)
!!$    logical , intent(in ) :: xy_FlagSlabOcean(0:imax-1, 1:jmax)
!!$    real(DP), intent(in ) :: xy_SOSeaIceMass       (0:imax-1, 1:jmax)
!!$    real(DP), intent(in ) :: xy_SurfTemp           (0:imax-1, 1:jmax)
!!$    real(DP), intent(in ) :: xyz_SOSeaIceThickness (0:imax-1, 1:jmax, 1:ksimax)
!!$    real(DP), intent(in ) :: xyz_SOSeaIceTemp      (0:imax-1, 1:jmax, 1:ksimax)
!!$    real(DP), intent(out) :: xy_SurfTempA          (0:imax-1, 1:jmax)
!!$    real(DP), intent(out) :: xyz_SOSeaIceThicknessA(0:imax-1, 1:jmax, 1:ksimax)
!!$    real(DP), intent(out) :: xyz_SOSeaIceTempA     (0:imax-1, 1:jmax, 1:ksimax)
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    real(DP) :: SOSeaIceU0 = 0.1_DP
!!$
!!$    real(DP) :: y_Factor(1:jmax)
!!$
!!$    real(DP) :: py_SOSeaIceU(0:imax-1, 1:jmax)
!!$                              ! 
!!$                              ! Longitudional Flux of slab sea ice
!!$
!!$    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:: l
!!$    integer:: ii
!!$    integer:: iprev
!!$    integer:: inext
!!$
!!$
!!$    ! 実行文 ; Executable statement
!!$    !
!!$
!!$    ! 初期化確認
!!$    ! Initialization check
!!$    !
!!$    if ( .not. sosi_dynamics_inited ) then
!!$      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
!!$    end if
!!$
!!$
!!$    y_Factor = 1.0_DP / ( RPlanet * y_CosLat * DelLon )**2
!!$
!!$    ! grid arrangement
!!$    !
!!$    !       |   .   |   .   |   .   |   .   |   .   |   .
!!$    !
!!$    ! x         0       1       2       3       4
!!$    ! p             0       1       2       3       4
!!$    !
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        if ( i == imax-1 ) then
!!$          iprev = i
!!$          inext = 0
!!$        else
!!$          iprev = i
!!$          inext = i+1
!!$        end if
!!$        if ( xy_FlagSlabOcean(iprev,j) .and. xy_FlagSlabOcean(inext,j) ) then
!!$          py_SOSeaIceU(i,j) = SOSeaIceU0
!!$        else
!!$          py_SOSeaIceU(i,j) = 0.0_DP
!!$        end if
!!$      end do
!!$    end do
!!$
!!$!    do i = 0, imax-1
!!$!      p_Lon(i) = DelLon / 2.0_DP + DelLon * i
!!$!    end do
!!$    p_Lon = x_Lon + ( x_Lon(1) - x_Lon(0) ) / 2.0_DP
!!$
!!$    ! Departure point
!!$    do j = 1, jmax
!!$      py_DPLon(:,j) = &
!!$        & p_Lon - py_SOSeaIceU(:,j) * DelTime / ( RPlanet * y_CosLat(j) )
!!$    end do
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        if ( py_DPLon(i,j) < 0.0_DP ) then
!!$          py_DPLon(i,j) = py_DPLon(i,j) + 2.0_DP * PI
!!$        else if ( py_DPLon(i,j) > 2.0_DP * PI ) then
!!$          py_DPLon(i,j) = py_DPLon(i,j) - 2.0_DP * PI
!!$        end if
!!$      end do
!!$    end do
!!$
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        dp_search : do ii = 0, imax-1
!!$          if ( py_DPLon(i,j) <= p_Lon(ii) ) then
!!$            exit dp_search
!!$          end if
!!$        end do dp_search
!!$        py_DPIndex
!!$        if (
!!$        py_SOSeaIceDPLon(i,j) = p_Lon
!!$      end do
!!$    end do
!!$
!!$
!!$
!!$
!!$  end subroutine SOSIFFSLX

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

  subroutine SOSIHorTransportDiff(                    &
    & xy_FlagSlabOcean,                               & ! (in)
    & xy_SOSeaIceMass,                                & ! (in)
    & xy_SurfTemp,                                    & ! (in)
    & xyz_SOSeaIceThickness, xyz_SOSeaIceTemp,        & ! (in)
    & xy_DSurfTempDt,                                 & ! (out)
    & xyz_DSOSeaIceThicknessDt, xyz_DSOSeaIceTempDt   & ! (out)
    & )
    ! 
    ! Calculates slab sea ice transport by diffusion

    use axesset    , only : x_Lon
                              ! $\lambda, \varphai$ lon and lat
    use sltt_const , only : dtjw, iexmin, iexmax, jexmin, jexmax
    use sosi_dynamics_extarr, only : SLTTExtArrExt2
                              ! 配列拡張ルーチン
                              ! Expansion of arrays
    ! 雪と海氷の定数の設定
    ! Setting constants of snow and sea ice
    !
    use constants_snowseaice, only: &
      & SeaIceDen


    logical , intent(in ) :: xy_FlagSlabOcean(0:imax-1, 1:jmax)
    real(DP), intent(in ) :: xy_SOSeaIceMass      (0:imax-1, 1:jmax)
    real(DP), intent(in ) :: xy_SurfTemp          (0:imax-1, 1:jmax)
    real(DP), intent(in ) :: xyz_SOSeaIceThickness(0:imax-1, 1:jmax, 1:ksimax)
    real(DP), intent(in ) :: xyz_SOSeaIceTemp     (0:imax-1, 1:jmax, 1:ksimax)
    real(DP), intent(out) :: xy_DSurfTempDt          (0:imax-1, 1:jmax)
    real(DP), intent(out) :: xyz_DSOSeaIceThicknessDt(0:imax-1, 1:jmax, 1:ksimax)
    real(DP), intent(out) :: xyz_DSOSeaIceTempDt     (0:imax-1, 1:jmax, 1:ksimax)


    !
    ! local variables
    !
    real(DP) :: xy_SurfTempA          (0:imax-1, 1:jmax)
    real(DP) :: xyz_SOSeaIceThicknessA(0:imax-1, 1:jmax, 1:ksimax)
    real(DP) :: xyz_SOSeaIceTempA     (0:imax-1, 1:jmax, 1:ksimax)

    real(DP) :: xy_SOSeaIceMassA      (0:imax-1, 1:jmax)

    real(DP) :: xyza_TMP4DArray(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)

    real(DP) :: xyza_ExtTMP4DArrayS(iexmin:iexmax, jexmin:jexmax, 1:kmax, 1:ncmax)
                              ! 
                              ! Extended 4D array (SH)
    real(DP) :: xyza_ExtTMP4DArrayN(iexmin:iexmax, jexmin:jexmax, 1:kmax, 1:ncmax)
                              ! 
                              ! Extended 4D array (NH)

    real(DP) :: xyz_ExtSOSeaIceThicknessS(iexmin:iexmax, jexmin:jexmax, 1:ksimax)
    real(DP) :: xyz_ExtSOSeaIceTempS     (iexmin:iexmax, jexmin:jexmax, 1:ksimax)
    real(DP) :: xyz_ExtSurfTempS         (iexmin:iexmax, jexmin:jexmax, 1:ksimax)
    real(DP) :: xy_ExtSOSeaIceMassS      (iexmin:iexmax, jexmin:jexmax)
                              ! 
                              ! Extended 4D array (SH)
    real(DP) :: xyz_ExtSOSeaIceThicknessN(iexmin:iexmax, jexmin:jexmax, 1:ksimax)
    real(DP) :: xyz_ExtSOSeaIceTempN     (iexmin:iexmax, jexmin:jexmax, 1:ksimax)
    real(DP) :: xyz_ExtSurfTempN         (iexmin:iexmax, jexmin:jexmax, 1:ksimax)
    real(DP) :: xy_ExtSOSeaIceMassN      (iexmin:iexmax, jexmin:jexmax)
                              ! 
                              ! Extended 4D array (NH)

    logical  :: xy_ExtFlagSlabOceanS(iexmin:iexmax, jexmin:jexmax)
                              ! 
                              ! Extended 4D array (SH)
    logical  :: xy_ExtFlagSlabOceanN(iexmin:iexmax, jexmin:jexmax)
                              ! 
                              ! Extended 4D array (NH)

    real(DP) :: xyz_DSOSeaIceThicknessDtS(0:imax-1, 1:jmax/2, 1:ksimax)
    real(DP) :: xyz_DSOSeaIceTempDtS     (0:imax-1, 1:jmax/2, 1:ksimax)
    real(DP) :: xya_DSurfTempDtS         (0:imax-1, 1:jmax/2, 1:ksimax)
    real(DP) :: xyz_DSOSeaIceThicknessDtN(0:imax-1, 1:jmax/2, 1:ksimax)
    real(DP) :: xyz_DSOSeaIceTempDtN     (0:imax-1, 1:jmax/2, 1:ksimax)
    real(DP) :: xya_DSurfTempDtN         (0:imax-1, 1:jmax/2, 1:ksimax)



    real(DP) :: PM
    ! 配列拡張する際、極ごえ後に符号が変わる場合は -1.0を与える。
    ! そうでない場合は1.0を与える。
    ! Sign change flag for array extension; -1.0 for sign change
    ! over the pole, 1.0 for no sign change



    integer:: i               ! 東西方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in zonal direction
    integer:: j               ! 南北方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in meridional direction
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
    integer:: n               ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituents
    integer:: kk


    ! 実行文 ; Executable statement
    !

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


    !
    ! Longitudinal diffusion
    !

#if defined(AXISYMMETRY) || defined(AXISYMMETRY_SJPACK)
    xy_SurfTempA           = xy_SurfTemp
    xyz_SOSeaIceThicknessA = xyz_SOSeaIceThickness
    xyz_SOSeaIceTempA      = xyz_SOSeaIceTemp
#else
    call SOSIDiffusionX(                               &
      & x_Lon(1)-x_Lon(0), y_CosLat, xy_FlagSlabOcean, & ! (in)
      & xy_SOSeaIceMass,                                             & ! (in)
      & xy_SurfTemp , xyz_SOSeaIceThickness , xyz_SOSeaIceTemp,      & ! (in)
      & xy_SurfTempA, xyz_SOSeaIceThicknessA, xyz_SOSeaIceTempA      & ! (out)
      & )
#endif

    xy_SOSeaIceMassA = 0.0_DP
    do k = 1, ksimax
      xy_SOSeaIceMassA = xy_SOSeaIceMassA &
        & + SeaIceDen * xyz_SOSeaIceThicknessA(:,:,k)
    end do


    !
    ! Latitudinal diffusion
    !

    ! 配列の分割と拡張
    ! Division and extension of arrays
    !

    if ( ksimax > kmax ) then
      call MessageNotify( 'E', module_name, 'ksimax > kmax.' )
    end if
    if ( kmax < 3 ) then
      call MessageNotify( 'E', module_name, 'kmax < 3.' )
    end if
    if ( ncmax < 3 ) then
      call MessageNotify( 'E', module_name, 'ncmax < 3.' )
    end if


    n = 1
    do k = 1, ksimax
      xyza_TMP4DArray(:,:,k,n) = xyz_SOSeaIceThicknessA(:,:,k)
    end do
    do k = ksimax+1, kmax
      xyza_TMP4DArray(:,:,k,n) = 0.0_DP
    end do
    !
    n = 2
    do k = 1, ksimax
      xyza_TMP4DArray(:,:,k,n) = xyz_SOSeaIceTempA(:,:,k)
    end do
    do k = ksimax+1, kmax
      xyza_TMP4DArray(:,:,k,n) = 0.0_DP
    end do
    !
    n = 3
    k = 1
    do j = 1, jmax
      do i = 0, imax-1
        if ( xy_FlagSlabOcean(i,j) ) then
          xyza_TMP4DArray(i,j,k,n) = 1.0_DP
        else
          xyza_TMP4DArray(i,j,k,n) = 0.0_DP
        end if
      end do
    end do
    k = 2
    xyza_TMP4DArray(:,:,k,n) = xy_SurfTempA
    k = 3
    xyza_TMP4DArray(:,:,k,n) = xy_SOSeaIceMassA
    do k = 3+1, kmax
      xyza_TMP4DArray(:,:,k,n) = 0.0_DP
    end do
    !
    do n = 3+1, ncmax
      xyza_TMP4DArray(:,:,:,n) = 0.0_DP
    end do

    PM = 1.0_DP
    call SLTTExtArrExt2(                            &
      & x_SinLonS, x_CosLonS, x_SinLonN, x_CosLonN, & ! (in)
      & xyza_TMP4DArray, PM,                        & ! (in)
      & xyza_ExtTMP4DArrayS, xyza_ExtTMP4DArrayN    & ! (out)
      & )

    n = 1
    do k = 1, ksimax
      xyz_ExtSOSeaIceThicknessS(:,:,k) = xyza_ExtTMP4DArrayS(:,:,k,n)
      xyz_ExtSOSeaIceThicknessN(:,:,k) = xyza_ExtTMP4DArrayN(:,:,k,n)
    end do
    n = 2
    do k = 1, ksimax
      xyz_ExtSOSeaIceTempS     (:,:,k) = xyza_ExtTMP4DArrayS(:,:,k,n)
      xyz_ExtSOSeaIceTempN     (:,:,k) = xyza_ExtTMP4DArrayN(:,:,k,n)
    end do
    n = 3
    k = 1
    do j = jexmin, jexmax
      do i = iexmin, iexmax
        if ( xyza_ExtTMP4DArrayS(i,j,k,n) == 1.0_DP ) then
          xy_ExtFlagSlabOceanS(i,j) = .true.
        else
          xy_ExtFlagSlabOceanS(i,j) = .false.
        end if
        if ( xyza_ExtTMP4DArrayN(i,j,k,n) == 1.0_DP ) then
          xy_ExtFlagSlabOceanN(i,j) = .true.
        else
          xy_ExtFlagSlabOceanN(i,j) = .false.
        end if
      end do
    end do
    k = 2
    do kk = 1, ksimax
      xyz_ExtSurfTempS(:,:,kk) = xyza_ExtTMP4DArrayS(:,:,k,n)
      xyz_ExtSurfTempN(:,:,kk) = xyza_ExtTMP4DArrayN(:,:,k,n)
    end do
    k = 3
    do kk = 1, ksimax
      xy_ExtSOSeaIceMassS = xyza_ExtTMP4DArrayS(:,:,k,n)
      xy_ExtSOSeaIceMassN = xyza_ExtTMP4DArrayN(:,:,k,n)
    end do


!!$    call SSIDiffusion(                             &
!!$      & x_ExtLonS, y_ExtLatS, y_ExtCosLatS,        & ! (in)
!!$      & p_LonS, q_LatS, q_CosLatS,                 & ! (in)
!!$      & xy_ExtFlagSlabOceanS, xy_ExtSOSeaIceMassS, & ! (in)
!!$      & xy_DSOSeaIceMassDtS                        & ! (out)
!!$      & )
!!$    call SSIDiffusion(                             &
!!$      & x_ExtLonN, y_ExtLatN, y_ExtCosLatN,        & ! (in)
!!$      & p_LonN, q_LatN, q_CosLatN,                 & ! (in)
!!$      & xy_ExtFlagSlabOceanN, xy_ExtSOSeaIceMassN, & ! (in)
!!$      & xy_DSOSeaIceMassDtN                        & ! (out)
!!$      & )

    call SOSIDiffusionY(                                 &
      & y_ExtLatS, y_ExtCosLatS,                         & ! (in)
      & q_LatS, q_CosLatS,                               & ! (in)
      & xy_ExtFlagSlabOceanS, xyz_ExtSOSeaIceThicknessS, & ! (in)
      & xyz_DSOSeaIceThicknessDtS                        & ! (out)
      & )
    call SOSIDiffusionY(                                 &
      & y_ExtLatN, y_ExtCosLatN,                         & ! (in)
      & q_LatN, q_CosLatN,                               & ! (in)
      & xy_ExtFlagSlabOceanN, xyz_ExtSOSeaIceThicknessN, & ! (in)
      & xyz_DSOSeaIceThicknessDtN                        & ! (out)
      & )
    !
    call SOSIDiffusionY(                                 &
      & y_ExtLatS, y_ExtCosLatS,                         & ! (in)
      & q_LatS, q_CosLatS,                               & ! (in)
      & xy_ExtFlagSlabOceanS, xyz_ExtSOSeaIceTempS, & ! (in)
      & xyz_DSOSeaIceTempDtS                        & ! (out)
      & )
    call SOSIDiffusionY(                                 &
      & y_ExtLatN, y_ExtCosLatN,                         & ! (in)
      & q_LatN, q_CosLatN,                               & ! (in)
      & xy_ExtFlagSlabOceanN, xyz_ExtSOSeaIceTempN, & ! (in)
      & xyz_DSOSeaIceTempDtN                        & ! (out)
      & )
    !
    call SOSIDiffusionY(                                 &
      & y_ExtLatS, y_ExtCosLatS,                         & ! (in)
      & q_LatS, q_CosLatS,                               & ! (in)
      & xy_ExtFlagSlabOceanS, xyz_ExtSurfTempS,     & ! (in)
      & xya_DSurfTempDtS,                           & ! (out)
      & xy_ExtSOSeaIceMassS                         & ! (in) optional
      & )
    call SOSIDiffusionY(                                 &
      & y_ExtLatN, y_ExtCosLatN,                         & ! (in)
      & q_LatN, q_CosLatN,                               & ! (in)
      & xy_ExtFlagSlabOceanN, xyz_ExtSurfTempN,     & ! (in)
      & xya_DSurfTempDtN,                           & ! (out)
      & xy_ExtSOSeaIceMassN                         & ! (in) optional
      & )

    xyz_DSOSeaIceThicknessDt(:,1:jmax/2     ,:) = xyz_DSOSeaIceThicknessDtS(:,1:jmax/2,:)
    xyz_DSOSeaIceThicknessDt(:,jmax/2+1:jmax,:) = xyz_DSOSeaIceThicknessDtN(:,1:jmax/2,:)
    !
    xyz_DSOSeaIceTempDt(:,1:jmax/2     ,:) = xyz_DSOSeaIceTempDtS(:,1:jmax/2,:)
    xyz_DSOSeaIceTempDt(:,jmax/2+1:jmax,:) = xyz_DSOSeaIceTempDtN(:,1:jmax/2,:)
    !
    xy_DSurfTempDt(:,1:jmax/2     ) = xya_DSurfTempDtS(:,1:jmax/2,1)
    xy_DSurfTempDt(:,jmax/2+1:jmax) = xya_DSurfTempDtN(:,1:jmax/2,1)


!!$    ! sea ice mass at next time step is calculated temporarily
!!$    xy_SOSeaIceMassA = xy_SOSeaIceMassA &
!!$      & + xy_DSOSeaIceMassDt * ( 2.0_DP * DelTime )
!!$
!!$    ! tendency is calculated
!!$    xy_DSOSeaIceMassDt = &
!!$      & ( xy_SOSeaIceMassA - xy_SOSeaIceMassB ) / ( 2.0_DP * DelTime )


!!$    py_ExtSOSeaIceMassXFlux = 0.0_DP
!!$    xq_ExtSOSeaIceMassYFlux = 0.0_DP
!!$    do j = jexmin, jexmax-1
!!$      do i = iexmin, iexmax-1
!!$        py_ExtSOSeaIceMassXFlux(i,j) = &
!!$          & - SOSeaIceDiffCoef &
!!$          &   * ( xy_ExtSOSeaIceMassS(i+1,j) - xy_ExtSOSeaIceMassS(i,j) ) &
!!$          &   / ( RPlanet * y_CosLatS(j) * ( x_ExtLonS(i+1) - x_ExtLonS(i) ) )
!!$        xq_ExtSOSeaIceMassYFlux(i,j) = &
!!$          & - SOSeaIceDiffCoef &
!!$          &   * ( xy_ExtSOSeaIceMassS(i,j+1) - xy_ExtSOSeaIceMassS(i,j) ) &
!!$          &   / ( RPlanet * ( y_ExtLatS(j+1) - y_ExtLatS(j) ) )
!!$      end do
!!$    end do
!!$    if ( myrank == nprocs-1 ) then
!!$      j = 0
!!$      do i = iexmin, iexmax-1
!!$        xq_ExtSOSeaIceMassYFlux(i,j) = 0.0_DP
!!$      end do
!!$    end if
!!$
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        xy_DSOSeaIceMassDt(i,j) = &
!!$          & - ( py_ExtSOSeaIceMassXFlux(i,j) - py_ExtSOSeaIceMassXFlux(i-1,j) ) &
!!$          &   / ( RPlanet * y_CosLatS(j) * ( p_Lon(i) - p_Lon(i-1) ) ) &
!!$          & - (   xq_ExtSOSeaIceMassYFlux(i,j  ) * q_CosLatS(j  )   &
!!$          &     - xq_ExtSOSeaIceMassYFlux(i,j-1) * q_CosLatS(j-1) ) &
!!$          &   / ( RPlanet * y_CosLatS(j) * ( q_LatS(j) - q_LatS(j-1) ) )
!!$      end do
!!$    end do




  end subroutine SOSIHorTransportDiff

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

  subroutine SOSIDiffusion(                      &
    & x_ExtLonH, y_ExtLatH, y_ExtCosLatH,        & ! (in)
    & p_LonH, q_LatH, q_CosLatH,                 & ! (in)
    & xy_ExtFlagSlabOceanH, xy_ExtSOSeaIceMassH, & ! (in)
    & xy_DSOSeaIceMassDtH                        & ! (out)
    & )

    ! 
    ! Calculates slab sea ice transport by diffusion

    use sltt_const , only : iexmin, iexmax, jexmin, jexmax
                              ! 配列拡張ルーチン
                              ! Expansion of arrays

    use constants, only: &
      & RPlanet
      ! $ a $ [m].
      ! 惑星半径.
      ! Radius of planet

    real(DP), intent(in ) :: x_ExtLonH   (iexmin:iexmax)
    real(DP), intent(in ) :: y_ExtLatH   (jexmin:jexmax)
    real(DP), intent(in ) :: y_ExtCosLatH(jexmin:jexmax)
    real(DP), intent(in ) :: p_LonH   (0-1:imax-1+1)
    real(DP), intent(in ) :: q_LatH   (1-1:jmax/2+1)
    real(DP), intent(in ) :: q_CosLatH(1-1:jmax/2+1)
    real(DP), intent(in ) :: xy_ExtSOSeaIceMassH(iexmin:iexmax, jexmin:jexmax)
                              ! 
                              ! Extended 4D array
    logical , intent(in ) :: xy_ExtFlagSlabOceanH(iexmin:iexmax, jexmin:jexmax)
                              ! 
                              ! Extended 4D array
    real(DP), intent(out) :: xy_DSOSeaIceMassDtH(0:imax-1, 1:jmax/2)
                              ! 
                              ! Slab sea ice mass tendency

    !
    ! local variables
    !

    real(DP) :: py_ExtSOSeaIceMassXFlux(iexmin-1:iexmax, jexmin:jexmax)
                              ! 
                              ! Longitudional Flux of slab sea ice
    real(DP) :: xq_ExtSOSeaIceMassYFlux(iexmin:iexmax, jexmin-1:jexmax)
                              ! 
                              ! Latitudinal Flux of slab sea ice


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


    ! 実行文 ; Executable statement
    !

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


    py_ExtSOSeaIceMassXFlux = 0.0_DP
    xq_ExtSOSeaIceMassYFlux = 0.0_DP
    do j = jexmin, jexmax-1
      do i = iexmin, iexmax-1
        if ( xy_ExtFlagSlabOceanH(i,j) .and. xy_ExtFlagSlabOceanH(i+1,j) ) then
          py_ExtSOSeaIceMassXFlux(i,j) =                                    &
            & - SOSeaIceDiffCoef                                            &
            &   * ( xy_ExtSOSeaIceMassH(i+1,j) - xy_ExtSOSeaIceMassH(i,j) ) &
            &   / ( RPlanet * y_ExtCosLatH(j) * ( x_ExtLonH(i+1) - x_ExtLonH(i) ) )
        else
          py_ExtSOSeaIceMassXFlux(i,j) = 0.0_DP
        end if
        if ( xy_ExtFlagSlabOceanH(i,j) .and. xy_ExtFlagSlabOceanH(i,j+1) ) then
          xq_ExtSOSeaIceMassYFlux(i,j) =                                    &
            & - SOSeaIceDiffCoef                                            &
            &   * ( xy_ExtSOSeaIceMassH(i,j+1) - xy_ExtSOSeaIceMassH(i,j) ) &
            &   / ( RPlanet * ( y_ExtLatH(j+1) - y_ExtLatH(j) ) )
        else
          xq_ExtSOSeaIceMassYFlux(i,j) = 0.0_DP
        end if
      end do
    end do
!!$    if ( myrank == nprocs-1 ) then
!!$      j = 0
!!$      do i = iexmin, iexmax-1
!!$        xq_ExtSOSeaIceMassYFlux(i,j) = 0.0_DP
!!$      end do
!!$    end if

    do j = 1, jmax/2
      do i = 0, imax-1
        xy_DSOSeaIceMassDtH(i,j) =                 &
          & - (   py_ExtSOSeaIceMassXFlux(i  ,j)   &
          &     - py_ExtSOSeaIceMassXFlux(i-1,j) ) &
          &   / ( RPlanet * y_ExtCosLatH(j) * ( p_LonH(i) - p_LonH(i-1) ) ) &
          & - (   xq_ExtSOSeaIceMassYFlux(i,j  ) * q_CosLatH(j  )           &
          &     - xq_ExtSOSeaIceMassYFlux(i,j-1) * q_CosLatH(j-1) )         &
          &   / ( RPlanet * y_ExtCosLatH(j) * ( q_LatH(j) - q_LatH(j-1) ) )
      end do
    end do


  end subroutine SOSIDiffusion

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

  subroutine SOSIDiffusionY(                      &
    & y_ExtLatH, y_ExtCosLatH,                    & ! (in)
    & q_LatH, q_CosLatH,                          & ! (in)
    & xy_ExtFlagSlabOceanH, xyz_ExtSOSeaIceMassH, & ! (in)
    & xyz_DSOSeaIceMassDtH,                       & ! (out)
    & xy_ExtSOSeaIceMassH                         & ! (in) optional
    & )

    ! 
    ! Calculates slab sea ice transport by diffusion

    use sltt_const , only : iexmin, iexmax, jexmin, jexmax
                              ! 配列拡張ルーチン
                              ! Expansion of arrays

    use constants, only: &
      & RPlanet, &
      ! $ a $ [m].
      ! 惑星半径.
      ! Radius of planet
      & SOMass
                              ! Slab ocean mass

    real(DP), intent(in ) :: y_ExtLatH   (jexmin:jexmax)
    real(DP), intent(in ) :: y_ExtCosLatH(jexmin:jexmax)
    real(DP), intent(in ) :: q_LatH   (1-1:jmax/2+1)
    real(DP), intent(in ) :: q_CosLatH(1-1:jmax/2+1)
    logical , intent(in ) :: xy_ExtFlagSlabOceanH(iexmin:iexmax, jexmin:jexmax)
                              ! 
                              ! Extended 4D array
    real(DP), intent(in ) :: xyz_ExtSOSeaIceMassH(iexmin:iexmax, jexmin:jexmax, 1:ksimax)
                              ! 
                              ! Extended 4D array
    real(DP), intent(out) :: xyz_DSOSeaIceMassDtH(0:imax-1, 1:jmax/2, 1:ksimax)
                              ! 
                              ! Slab sea ice mass tendency
    real(DP), intent(in ), optional :: xy_ExtSOSeaIceMassH(iexmin:iexmax, jexmin:jexmax)


    !
    ! local variables
    !

    real(DP) :: xqz_ExtSODiffCoef       (iexmin:iexmax, jexmin-1:jexmax, 1:ksimax)
    real(DP) :: xqz_ExtSOSeaIceMassYFlux(iexmin:iexmax, jexmin-1:jexmax, 1:ksimax)
                              ! 
                              ! Latitudinal Flux of slab sea ice


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


    ! 実行文 ; Executable statement
    !

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


    do k = 1, ksimax
      do j = jexmin, jexmax-1
        do i = iexmin, iexmax-1
          if ( xy_ExtFlagSlabOceanH(i,j) .and. xy_ExtFlagSlabOceanH(i,j+1) ) then
            xqz_ExtSODiffCoef(i,j,k) = SOSeaIceDiffCoef
            if ( present( xy_ExtSOSeaIceMassH ) ) then
              xqz_ExtSODiffCoef(i,j,k) = xqz_ExtSODiffCoef(i,j,k) &
                & * min( SOMass - xy_ExtSOSeaIceMassH(i,j  ), &
                &        SOMass - xy_ExtSOSeaIceMassH(i,j+1) )
            end if
          else
            xqz_ExtSODiffCoef(i,j,k) = 0.0_DP
          end if
        end do
      end do
    end do

    do k = 1, ksimax
      do j = jexmin, jexmax-1
        do i = iexmin, iexmax-1
          xqz_ExtSOSeaIceMassYFlux(i,j,k) =                                       &
            & - xqz_ExtSODiffCoef(i,j,k)                                          &
            &   * ( xyz_ExtSOSeaIceMassH(i,j+1,k) - xyz_ExtSOSeaIceMassH(i,j,k) ) &
            &   / ( RPlanet * ( y_ExtLatH(j+1) - y_ExtLatH(j) ) )
        end do
      end do
    end do
!!$    if ( myrank == nprocs-1 ) then
!!$      j = 0
!!$      do i = iexmin, iexmax-1
!!$        xq_ExtSOSeaIceMassYFlux(i,j) = 0.0_DP
!!$      end do
!!$    end if

    do k = 1, ksimax
      do j = 1, jmax/2
        do i = 0, imax-1
          xyz_DSOSeaIceMassDtH(i,j,k) =                                      &
            & - (   xqz_ExtSOSeaIceMassYFlux(i,j  ,k) * q_CosLatH(j  )       &
            &     - xqz_ExtSOSeaIceMassYFlux(i,j-1,k) * q_CosLatH(j-1) )     &
            &   / ( RPlanet * y_ExtCosLatH(j) * ( q_LatH(j) - q_LatH(j-1) ) )
          if ( present( xy_ExtSOSeaIceMassH ) ) then
            if ( SOMass - xy_ExtSOSeaIceMassH(i,j) > 0.0_DP ) then
              xyz_DSOSeaIceMassDtH(i,j,k) = xyz_DSOSeaIceMassDtH(i,j,k) &
                & / ( SOMass - xy_ExtSOSeaIceMassH(i,j) )
            else
              xyz_DSOSeaIceMassDtH(i,j,k) = 0.0_DP
            end if
          end if
        end do
      end do
    end do


  end subroutine SOSIDiffusionY

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

  subroutine SOSIDiffusionX(                     &
    & DelLon, y_CosLat, xy_FlagSlabOcean,        & ! (in)
    & xy_SOSeaIceMass,                           & ! (in)
    & xy_SurfTemp , xyz_SOSeaIceThickness , xyz_SOSeaIceTemp,  & ! (in)
    & xy_SurfTempA, xyz_SOSeaIceThicknessA, xyz_SOSeaIceTempA  & ! (out)
    & )
    ! 
    ! Calculates slab sea ice transport by longitudinal diffusion


    !
    !
    use ludecomp_module, only : &
      & ludecomp_prep_simple_many, &
      & ludecomp_solve_simple_many

    use timeset    , only : DelTime
                              ! $\Delta t$
    use constants, only: &
      & RPlanet, &
      ! $ a $ [m].
      ! 惑星半径.
      ! Radius of planet
      & SOMass
                              ! Slab ocean mass

    real(DP), intent(in ) :: DelLon
    real(DP), intent(in ) :: y_CosLat(1:jmax)
    logical , intent(in ) :: xy_FlagSlabOcean(0:imax-1, 1:jmax)
    real(DP), intent(in ) :: xy_SOSeaIceMass       (0:imax-1, 1:jmax)
    real(DP), intent(in ) :: xy_SurfTemp           (0:imax-1, 1:jmax)
    real(DP), intent(in ) :: xyz_SOSeaIceThickness (0:imax-1, 1:jmax, 1:ksimax)
    real(DP), intent(in ) :: xyz_SOSeaIceTemp      (0:imax-1, 1:jmax, 1:ksimax)
    real(DP), intent(out) :: xy_SurfTempA          (0:imax-1, 1:jmax)
    real(DP), intent(out) :: xyz_SOSeaIceThicknessA(0:imax-1, 1:jmax, 1:ksimax)
    real(DP), intent(out) :: xyz_SOSeaIceTempA     (0:imax-1, 1:jmax, 1:ksimax)

    !
    ! local variables
    !

    real(DP) :: aax_LUMat(1:jmax*ksimax, 0:imax-1, 0:imax-1)
    real(DP) :: aa_LUVec (1:jmax*ksimax, 0:imax-1)

    real(DP) :: y_Factor(1:jmax)

    real(DP) :: pyz_SOSeaIceDiffCoef(0:imax-1, 1:jmax, 1:ksimax)
                              ! 
                              ! Longitudional Flux of slab sea ice

    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:: l
    integer:: ii
    integer:: iprev
    integer:: inext


    ! 実行文 ; Executable statement
    !

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


    y_Factor = 1.0_DP / ( RPlanet * y_CosLat * DelLon )**2

    do k = 1, ksimax
      do j = 1, jmax
        do i = 0, imax-1
          if ( i == imax-1 ) then
            iprev = i
            inext = 0
          else
            iprev = i
            inext = i+1
          end if
          if ( xy_FlagSlabOcean(iprev,j) .and. xy_FlagSlabOcean(inext,j) ) then
            pyz_SOSeaIceDiffCoef(i,j,k) = SOSeaIceDiffCoef
          else
            pyz_SOSeaIceDiffCoef(i,j,k) = 0.0_DP
          end if
        end do
      end do
    end do

    aax_LUMat = 0.0_DP
    do k = 1, ksimax
      do j = 1, jmax
        l = (k-1)*jmax+j

        do ii = 0, 0
          i = ii
          aax_LUMat(l,ii,imax-1) = &
            & - pyz_SOSeaIceDiffCoef(imax-1,j,k) &
            &   * y_Factor(j)
          aax_LUMat(l,ii,i  ) = &
            &   1.0_DP / ( 1.0_DP * DelTime ) &
            & + ( pyz_SOSeaIceDiffCoef(imax-1,j,k) + pyz_SOSeaIceDiffCoef(i,j,k) ) &
            &   * y_Factor(j)
          aax_LUMat(l,ii,i+1) = &
            & - pyz_SOSeaIceDiffCoef(i     ,j,k) &
            &   * y_Factor(j)
        end do
        do ii = 0+1, (imax-1)-1
          i = ii
          aax_LUMat(l,ii,i-1) = &
            & - pyz_SOSeaIceDiffCoef(i-1,j,k) &
            &   * y_Factor(j)
          aax_LUMat(l,ii,i  ) = &
            &   1.0_DP / ( 1.0_DP * DelTime ) &
            & + ( pyz_SOSeaIceDiffCoef(i-1,j,k) + pyz_SOSeaIceDiffCoef(i,j,k) ) &
            &   * y_Factor(j)
          aax_LUMat(l,ii,i+1) = &
            & - pyz_SOSeaIceDiffCoef(i  ,j,k) &
            &   * y_Factor(j)
        end do
        do ii = imax-1, imax-1
          i = ii
          aax_LUMat(l,ii,i-1) = &
            & - pyz_SOSeaIceDiffCoef(i-1,j,k) &
            &   * y_Factor(j)
          aax_LUMat(l,ii,i  ) = &
            &   1.0_DP / ( 1.0_DP * DelTime ) &
            & + ( pyz_SOSeaIceDiffCoef(i-1,j,k) + pyz_SOSeaIceDiffCoef(i,j,k) ) &
            &   * y_Factor(j)
          aax_LUMat(l,ii,0  ) = &
            & - pyz_SOSeaIceDiffCoef(imax-1,j,k) &
            &   * y_Factor(j)
        end do


      end do
    end do

    call ludecomp_prep_simple_many ( jmax*ksimax, imax, aax_LUMat, 1, jmax*ksimax )


    do k = 1, ksimax
      do j = 1, jmax
        l = (k-1)*jmax+j
        do ii = 0, imax-1
          i = ii
          aa_LUVec(l,ii) = xyz_SOSeaIceThickness(i,j,k) / ( 1.0_DP * DelTime )
        end do
      end do
    end do
    call ludecomp_solve_simple_many( jmax*ksimax, imax, aax_LUMat, aa_LUVec, 1, jmax*ksimax )
    do k = 1, ksimax
      do j = 1, jmax
        l = (k-1)*jmax+j
        do ii = 0, imax-1
          i = ii
          xyz_SOSeaIceThicknessA(i,j,k) = aa_LUVec(l,ii)
        end do
      end do
    end do

    do k = 1, ksimax
      do j = 1, jmax
        l = (k-1)*jmax+j
        do ii = 0, imax-1
          i = ii
          aa_LUVec(l,ii) = xyz_SOSeaIceTemp(i,j,k) / ( 1.0_DP * DelTime )
        end do
      end do
    end do
    call ludecomp_solve_simple_many( jmax*ksimax, imax, aax_LUMat, aa_LUVec, 1, jmax*ksimax )
    do k = 1, ksimax
      do j = 1, jmax
        l = (k-1)*jmax+j
        do ii = 0, imax-1
          i = ii
          xyz_SOSeaIceTempA(i,j,k) = aa_LUVec(l,ii)
        end do
      end do
    end do


    ! Diffusion of slab ocean temperature
    !
    do k = 1, ksimax
      do j = 1, jmax
        do i = 0, imax-1
          if ( i == imax-1 ) then
            iprev = i
            inext = 0
          else
            iprev = i
            inext = i+1
          end if
          if ( xy_FlagSlabOcean(iprev,j) .and. xy_FlagSlabOcean(inext,j) ) then
            pyz_SOSeaIceDiffCoef(i,j,k) = SOSeaIceDiffCoef &
              & * min( SOMass - xy_SOSeaIceMass(iprev,j), &
              &        SOMass - xy_SOSeaIceMass(inext,j) )
          else
            pyz_SOSeaIceDiffCoef(i,j,k) = 0.0_DP
          end if
        end do
      end do
    end do

    aax_LUMat = 0.0_DP
    do k = 1, ksimax
      do j = 1, jmax
        l = (k-1)*jmax+j

        do ii = 0, 0
          i = ii
          aax_LUMat(l,ii,imax-1) = &
            & - pyz_SOSeaIceDiffCoef(imax-1,j,k) &
            &   * y_Factor(j)
          aax_LUMat(l,ii,i  ) = &
            &   1.0_DP / ( 1.0_DP * DelTime )          &
            &     * ( SOMass - xy_SOSeaIceMass(i,j) )  &
            & + ( pyz_SOSeaIceDiffCoef(imax-1,j,k) + pyz_SOSeaIceDiffCoef(i,j,k) ) &
            &   * y_Factor(j)
          aax_LUMat(l,ii,i+1) = &
            & - pyz_SOSeaIceDiffCoef(i     ,j,k) &
            &   * y_Factor(j)
        end do
        do ii = 0+1, (imax-1)-1
          i = ii
          aax_LUMat(l,ii,i-1) = &
            & - pyz_SOSeaIceDiffCoef(i-1,j,k) &
            &   * y_Factor(j)
          aax_LUMat(l,ii,i  ) = &
            &   1.0_DP / ( 1.0_DP * DelTime )          &
            &     * ( SOMass - xy_SOSeaIceMass(i,j) )  &
            & + ( pyz_SOSeaIceDiffCoef(i-1,j,k) + pyz_SOSeaIceDiffCoef(i,j,k) ) &
            &   * y_Factor(j)
          aax_LUMat(l,ii,i+1) = &
            & - pyz_SOSeaIceDiffCoef(i  ,j,k) &
            &   * y_Factor(j)
        end do
        do ii = imax-1, imax-1
          i = ii
          aax_LUMat(l,ii,i-1) = &
            & - pyz_SOSeaIceDiffCoef(i-1,j,k) &
            &   * y_Factor(j)
          aax_LUMat(l,ii,i  ) = &
            &   1.0_DP / ( 1.0_DP * DelTime )          &
            &     * ( SOMass - xy_SOSeaIceMass(i,j) )  &
            & + ( pyz_SOSeaIceDiffCoef(i-1,j,k) + pyz_SOSeaIceDiffCoef(i,j,k) ) &
            &   * y_Factor(j)
          aax_LUMat(l,ii,0  ) = &
            & - pyz_SOSeaIceDiffCoef(imax-1,j,k) &
            &   * y_Factor(j)
        end do


      end do
    end do
    !
    do k = 1, ksimax
      do j = 1, jmax
        l = (k-1)*jmax+j
        do ii = 0, 0
          i = ii
          if ( SOMass - xy_SOSeaIceMass(i,j) <= 0.0_DP ) then
            aax_LUMat(l,ii,imax-1) = 0.0_DP
            aax_LUMat(l,ii,i     ) = 1.0_DP
            aax_LUMat(l,ii,i+1   ) = 0.0_DP
          end if
        end do
        do ii = 0+1, (imax-1)-1
          i = ii
          if ( SOMass - xy_SOSeaIceMass(i,j) <= 0.0_DP ) then
            aax_LUMat(l,ii,i-1) = 0.0_DP
            aax_LUMat(l,ii,i  ) = 1.0_DP
            aax_LUMat(l,ii,i+1) = 0.0_DP
          end if
        end do
        do ii = imax-1, imax-1
          i = ii
          if ( SOMass - xy_SOSeaIceMass(i,j) <= 0.0_DP ) then
            aax_LUMat(l,ii,i-1) = 0.0_DP
            aax_LUMat(l,ii,i  ) = 1.0_DP
            aax_LUMat(l,ii,0  ) = 0.0_DP
          end if
        end do


      end do
    end do


    call ludecomp_prep_simple_many ( jmax*ksimax, imax, aax_LUMat, 1, jmax*ksimax )

    do k = 1, ksimax
      do j = 1, jmax
        l = (k-1)*jmax+j
        do ii = 0, imax-1
          i = ii
          if ( SOMass - xy_SOSeaIceMass(i,j) <= 0.0_DP ) then
            aa_LUVec(l,ii) = xy_SurfTemp(i,j)
          else
            aa_LUVec(l,ii) = xy_SurfTemp(i,j) &
              & * ( SOMass - xy_SOSeaIceMass(i,j) )  &
              & / ( 1.0_DP * DelTime )
          end if
        end do
      end do
    end do
    call ludecomp_solve_simple_many( jmax*ksimax, imax, aax_LUMat, aa_LUVec, 1, jmax*ksimax )
    do k = 1, 1
      do j = 1, jmax
        l = (k-1)*jmax+j
        do ii = 0, imax-1
          i = ii
          xy_SurfTempA(i,j) = aa_LUVec(l,ii)
        end do
      end do
    end do


  end subroutine SOSIDiffusionX

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

  subroutine SOSIDynamicsInit( &
    & ArgFlagSlabOcean &
    & )
    ! flag for use of slab ocean
    ! 
    ! Initialization of module


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

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

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

    use mpi_wrapper   , only : myrank, nprocs

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

    ! 組成に関わる配列の設定
    ! Settings of array for atmospheric composition
    !
    use composition, only:                              &
      &                    ncmax
                             ! 成分の数
                             ! Number of composition

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: &
      & x_Lon, y_Lat, &
      & AxNameX, AxNameY, AxNameZ, AxNameT

    use constants0, only : PI

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

    !
    ! Slab ocean sea ice utility module
    !
    use sosi_utils, only : SOSIUtilsInit

    use sltt_const , only : SLTTConstInit
    use sosi_dynamics_extarr, only : SLTTExtArrInit


    use sltt_const , only : iexmin, iexmax, jexmin, jexmax

    ! メッセージ制御
    ! Message control
    !
    use mpi_messagecntl, only : DoesOutputMPIMessage


    logical, intent(in ) :: ArgFlagSlabOcean


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

    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_dynamics_nml/ &
      & SOSeaIceDiffCoef


    ! 実行文 ; Executable statement
    !

    if ( sosi_dynamics_inited ) return



    FlagSlabOcean = ArgFlagSlabOcean


    if ( mod( jmax, 2 ) /= 0 ) then
      stop 'jmax cannot be divided by 2.'
    end if


    ! Initialization of modules used in this module
    !

    ! 雪と海氷の定数の設定
    ! Setting constants of snow and sea ice
    !
    call ConstantsSnowSeaIceInit

    !
    ! Slab ocean sea ice utility module
    !
    call SOSIUtilsInit


    call SLTTConstInit


    ! デフォルト値の設定
    ! Default values settings
    !
    SOSeaIceDiffCoef              =  0.0e0_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_dynamics_nml, &  ! (out)
        & iostat = iostat_nml )       ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
      if ( iostat_nml == 0 .AND. DoesOutputMPIMessage() ) write( STDOUT, nml = sosi_dynamics_nml )
    end if


    allocate( y_CosLat(1:jmax) )
    y_CosLat = cos( y_Lat )

    allocate( x_LonS   (0:imax-1) )
    allocate( x_SinLonS(0:imax-1) )
    allocate( x_CosLonS(0:imax-1) )
    allocate( y_latS   (1:jmax/2) )
    allocate( y_SinLatS(1:jmax/2) )
    allocate( y_CosLatS(1:jmax/2) )
    do i = 0, imax-1
      x_LonS   (i) = x_Lon(i)
      x_SinLonS(i) = sin( x_LonS(i) )
      x_CosLonS(i) = cos( x_LonS(i) )
    end do
    do j = 1, jmax/2
      y_LatS   (j) = y_Lat(j)
      y_SinLatS(j) = sin( y_LatS(j) )
      y_CosLatS(j) = cos( y_LatS(j) )
    end do

    allocate( x_LonN   (0:imax-1) )
    allocate( x_SinLonN(0:imax-1) )
    allocate( x_CosLonN(0:imax-1) )
    allocate( y_latN   (1:jmax/2) )
    allocate( y_SinLatN(1:jmax/2) )
    allocate( y_CosLatN(1:jmax/2) )
    do i = 0, imax-1
      x_LonN   (i) = x_Lon(i)
      x_SinLonN(i) = sin( x_LonN(i) )
      x_CosLonN(i) = cos( x_LonN(i) )
    end do
    do j = 1, jmax/2
      y_LatN   (j) = y_Lat(j+jmax/2)
      y_SinLatN(j) = sin( y_LatN(j) )
      y_CosLatN(j) = cos( y_LatN(j) )
    end do

    allocate( x_ExtLonS( iexmin:iexmax ) )
    allocate( x_ExtLonN( iexmin:iexmax ) )

    allocate( y_ExtLatS( jexmin:jexmax ) )
    allocate( y_ExtLatN( jexmin:jexmax ) )

    allocate( y_ExtCosLatS( jexmin:jexmax ) )
    allocate( y_ExtCosLatN( jexmin:jexmax ) )


    call SLTTExtArrInit(                            &
      & x_LonS, y_LatS, x_LonN, y_LatN,             & ! (in )
      & x_ExtLonS, y_ExtLatS, x_ExtLonN, y_ExtLatN  & ! (out)
      & )

    y_ExtCosLatS = cos( y_ExtLatS )
    y_ExtCosLatN = cos( y_ExtLatN )

    allocate( p_LonS   (0-1:imax-1+1) )
    allocate( q_LatS   (1-1:jmax/2+1) )
    allocate( q_CosLatS(1-1:jmax/2+1) )
    allocate( p_LonN   (0-1:imax-1+1) )
    allocate( q_LatN   (1-1:jmax/2+1) )
    allocate( q_CosLatN(1-1:jmax/2+1) )


    do i = 0-1, imax-1+1
      p_LonS(i) = ( x_ExtLonS(i) + x_ExtLonS(i+1) ) / 2.0_DP
      p_LonN(i) = ( x_ExtLonN(i) + x_ExtLonN(i+1) ) / 2.0_DP
    end do
    do j = 1-1, jmax/2+1
      q_LatS(j) = ( y_ExtLatS(j) + y_ExtLatS(j+1) ) / 2.0_DP
    end do
    do j = 1-1, jmax/2+1
      q_LatN(j) = ( y_ExtLatN(j) + y_ExtLatN(j+1) ) / 2.0_DP
    end do
    if ( myrank == nprocs-1 ) then
      j = 0
      q_LatS(j) = - PI / 2.0_DP
      j = jmax/2+1
      q_LatN(j) =   PI / 2.0_DP
    end if

    q_CosLatS = cos( q_LatS )
    q_CosLatN = cos( q_LatN )



    ! ヒストリデータ出力のためのへの変数登録
    ! Register of variables for history data output
    !
!!$    do n = 1, ncmax
!!$      call HistoryAutoAddVariable( 'SLD'//trim(a_QMixName(n))//'DtHorMassFix', &
!!$        & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
!!$        & 'tendency of horizontal mass fix of '//trim(a_QMixName(n)), 's-1' )
!!$      call HistoryAutoAddVariable( 'SLD'//trim(a_QMixName(n))//'DtVerMassFix', &
!!$        & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
!!$        & 'tendency of vertical mass fix of '//trim(a_QMixName(n)), 's-1' )
!!$      call HistoryAutoAddVariable( 'SLD'//trim(a_QMixName(n))//'DtTotMassFix', &
!!$        & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
!!$        & 'tendency of mass fix of '//trim(a_QMixName(n)), 's-1' )
!!$    end do


    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, '  SOSeaIceDiffCoef              = %f', d = (/ SOSeaIceDiffCoef /) )
!!$    call MessageNotify( 'M', module_name, '  FlagSLTTArcsineVer       = %b', l = (/ FlagSLTTArcsineVer /) )
!!$    call MessageNotify( 'M', module_name, '  SLTTArcsineFactor        = %f', d = (/ SLTTArcsineFactor /) )
!!$    call MessageNotify( 'M', module_name, '  SLTTIntHor               = %c', c1 = trim( SLTTIntHor ) )
!!$    call MessageNotify( 'M', module_name, '  SLTTIntVer               = %c', c1 = trim( SLTTIntVer ) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    sosi_dynamics_inited = .true.

  end subroutine SOSIDynamicsInit

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

end module sosi_dynamics