!= 
!
!= 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_diffusion
  !
  != 
  !
  != 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_diffusion'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: sltt.F90,v 1.8 2025/09/17 22:00:00 takepiro Exp $'
                              ! モジュールのバージョン
                              ! Module version


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

contains

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

  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 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)

  end subroutine SOSIHorTransportFFSL

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

  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 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)



  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   (0:jmax+1)
    real(DP), intent(in ) :: y_ExtCosLatH(0:jmax+1)
    real(DP), intent(in ) :: q_LatH   (0:jmax)
    real(DP), intent(in ) :: q_CosLatH(0:jmax)
    logical , intent(in ) :: xy_ExtFlagSlabOceanH(0:imax-1, 0:jmax+1)
                              ! 
                              ! Extended 4D array
    real(DP), intent(in ) :: xyz_ExtSOSeaIceMassH(0:imax-1, 0:jmax+1, 1:ksimax)
                              ! 
                              ! Extended 4D array
    real(DP), intent(out) :: xyz_DSOSeaIceMassDtH(0:imax-1, 1:jmax, 1:ksimax)
                              ! 
                              ! Slab sea ice mass tendency
    real(DP), intent(in ), optional :: xy_ExtSOSeaIceMassH(0:imax-1, 0:jmax+1)


    !
    ! local variables
    !

    real(DP) :: xqz_ExtSODiffCoef       (0:imax-1, 0:jmax, 1:ksimax)
    real(DP) :: xqz_ExtSOSeaIceMassYFlux(0:imax-1, 0:jmax, 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 = 0, jmax
          do i = 0, imax-1
             if ( xy_ExtFlagSlabOceanH(i,j) .and. xy_ExtFlagSlabOceanH(i,j+1) ) then
                xqz_ExtSODiffCoef(i,j,) = 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
    enddo
    
    do k = 1, ksimax
      do j = 0, jmax
        do i = 0, imax-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

    do k = 1, ksimax
      do j = 1, jmax
        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

end module sosi_diffusion