Class sosi_dynamics
In: sosi/sosi_dynamics.F90

Slab sea ice horizontal transport

Note that Japanese and English are described in parallel.

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

Methods

Included Modules

dc_types dc_message mpi_wrapper timeset gridset composition gtool_historyauto constants0 constants constants_snowseaice sosi_utils axesset sltt_const sosi_dynamics_extarr ludecomp_module dc_iounit namelist_util

Public Instance methods

Subroutine :
xy_SurfType(0:imax-1, 1:jmax) :integer , intent(in )
: 土地利用. Surface index
xy_SurfTemp(0:imax-1, 1:jmax) :real(DP), intent(inout)
xy_SOSeaIceMass(0:imax-1, 1:jmax) :real(DP), intent(inout)
: $ M_si $ . 海氷質量 (kg m-2) Slab ocean sea ice mass (kg m-2)
xyz_SOSeaIceTemp(0:imax-1, 1:jmax, 1:ksimax) :real(DP), intent(inout)
xy_DSOSeaIceMassDtPhyTop(0:imax-1, 1:jmax) :real(DP), intent(in )
xy_DSOSeaIceMassDtPhyBot(0:imax-1, 1:jmax) :real(DP), intent(in )
: Slab sea ice mass at next time step
xyz_DSOSeaIceTempDtPhy(0:imax-1, 1:jmax, 1:ksimax) :real(DP), intent(in )

Calculates slab sea ice horizontal transports by diffusion

[Source]

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

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

    use timeset    , only : TimeN, 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
    integer:: kk


    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, xy_DSOSeaIceMassDtPhyTop, xy_DSOSeaIceMassDtPhyBot, xyz_DSOSeaIceTempDtPhy )


!!$    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, xy_SOSILocalKMax, xyr_SOSILocalDepth, xyz_SOSILocalDepth )

    ! 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, xy_SOSeaIceMass, xy_SOTempAdv, xyz_SOSeaIceThicknessAdv, xyz_SOSeaIceTempAdv, xy_DSOTempDt, xyz_DSOSeaIceThicknessDt, xyz_DSOSeaIceTempDt )


    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, "SOSIDynamics" )


  end subroutine SOSIDynamics
Subroutine :
ArgFlagSlabOcean :logical, intent(in )

flag for use of slab ocean

Initialization of module

This procedure input/output NAMELIST#sosi_dynamics_nml .

[Source]

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


    !
    ! MPI
    !
    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: STDOUT, 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: r_Sigma, z_Sigma, 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



    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:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
    integer:: n

    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, namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml, nml = sosi_dynamics_nml, iostat = iostat_nml )       ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
      if ( iostat_nml == 0 ) 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, x_ExtLonS, y_ExtLatS, x_ExtLonN, y_ExtLatN )

    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

Private Instance methods

FlagSlabOcean
Variable :
FlagSlabOcean :logical, save
: flag for use of slab ocean
Subroutine :
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_ExtFlagSlabOceanH(iexmin:iexmax, jexmin:jexmax) :logical , intent(in )
: Extended 4D array
xy_ExtSOSeaIceMassH(iexmin:iexmax, jexmin:jexmax) :real(DP), intent(in )
: Extended 4D array
xy_DSOSeaIceMassDtH(0:imax-1, 1:jmax/2) :real(DP), intent(out)
: Slab sea ice mass tendency

[Source]

  subroutine SOSIDiffusion( x_ExtLonH, y_ExtLatH, y_ExtCosLatH, p_LonH, q_LatH, q_CosLatH, xy_ExtFlagSlabOceanH, xy_ExtSOSeaIceMassH, xy_DSOSeaIceMassDtH )

    ! 
    ! 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 :
DelLon :real(DP), intent(in )
y_CosLat(1:jmax) :real(DP), intent(in )
xy_FlagSlabOcean(0:imax-1, 1:jmax) :logical , 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(in )
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) :real(DP), intent(out)

Calculates slab sea ice transport by longitudinal diffusion

[Source]

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


    !
    !
    use ludecomp_module, only : ludecomp_prep_simple_many, ludecomp_solve_simple_many

    use constants, only: RPlanet, 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 :
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) :real(DP), intent(in )
xy_ExtFlagSlabOceanH(iexmin:iexmax, jexmin:jexmax) :logical , intent(in )
: Extended 4D array
xyz_ExtSOSeaIceMassH(iexmin:iexmax, jexmin:jexmax, 1:ksimax) :real(DP), intent(in )
: Extended 4D array
xyz_DSOSeaIceMassDtH(0:imax-1, 1:jmax/2, 1:ksimax) :real(DP), intent(out)
: Slab sea ice mass tendency
xy_ExtSOSeaIceMassH(iexmin:iexmax, jexmin:jexmax) :real(DP), intent(in ), optional

[Source]

  subroutine SOSIDiffusionY( y_ExtLatH, y_ExtCosLatH, q_LatH, q_CosLatH, xy_ExtFlagSlabOceanH, xyz_ExtSOSeaIceMassH, xyz_DSOSeaIceMassDtH, xy_ExtSOSeaIceMassH )

    ! 
    ! Calculates slab sea ice transport by diffusion

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

    use constants, only: RPlanet, 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 :
xy_FlagSlabOcean(0:imax-1, 1:jmax) :logical , 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(in )
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) :real(DP), intent(out)

Calculates slab sea ice transport by diffusion

[Source]

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

    use timeset    , only : DelTime
                              ! $\Delta t$
    use axesset    , only : x_Lon, y_Lat
                              ! $\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, xy_SOSeaIceMass, xy_SurfTemp , xyz_SOSeaIceThickness , xyz_SOSeaIceTemp, xy_SurfTempA, xyz_SOSeaIceThicknessA, xyz_SOSeaIceTempA )
#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, xyza_TMP4DArray, PM, xyza_ExtTMP4DArrayS, xyza_ExtTMP4DArrayN )

    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, q_LatS, q_CosLatS, xy_ExtFlagSlabOceanS, xyz_ExtSOSeaIceThicknessS, xyz_DSOSeaIceThicknessDtS )
    call SOSIDiffusionY( y_ExtLatN, y_ExtCosLatN, q_LatN, q_CosLatN, xy_ExtFlagSlabOceanN, xyz_ExtSOSeaIceThicknessN, xyz_DSOSeaIceThicknessDtN )
    !
    call SOSIDiffusionY( y_ExtLatS, y_ExtCosLatS, q_LatS, q_CosLatS, xy_ExtFlagSlabOceanS, xyz_ExtSOSeaIceTempS, xyz_DSOSeaIceTempDtS )
    call SOSIDiffusionY( y_ExtLatN, y_ExtCosLatN, q_LatN, q_CosLatN, xy_ExtFlagSlabOceanN, xyz_ExtSOSeaIceTempN, xyz_DSOSeaIceTempDtN )
    !
    call SOSIDiffusionY( y_ExtLatS, y_ExtCosLatS, q_LatS, q_CosLatS, xy_ExtFlagSlabOceanS, xyz_ExtSurfTempS, xya_DSurfTempDtS, xy_ExtSOSeaIceMassS )
    call SOSIDiffusionY( y_ExtLatN, y_ExtCosLatN, q_LatN, q_CosLatN, xy_ExtFlagSlabOceanN, xyz_ExtSurfTempN, xya_DSurfTempDtN, xy_ExtSOSeaIceMassN )

    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 :
xy_FlagSlabOcean(0:imax-1, 1:jmax) :logical , 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(in )
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) :real(DP), intent(out)

Calculates slab sea ice transport by diffusion

[Source]

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

    use timeset    , only : DelTime
                              ! $\Delta t$
    use axesset    , only : x_Lon, y_Lat
                              ! $\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, xy_SOSeaIceMass, xy_SurfTemp , xyz_SOSeaIceThickness , xyz_SOSeaIceTemp, xy_SurfTempA, xyz_SOSeaIceThicknessA, xyz_SOSeaIceTempA )
#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, xyza_TMP4DArray, PM, xyza_ExtTMP4DArrayS, xyza_ExtTMP4DArrayN )

    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, q_LatS, q_CosLatS, xy_ExtFlagSlabOceanS, xyz_ExtSOSeaIceThicknessS, xyz_DSOSeaIceThicknessDtS )
    call SOSIDiffusionY( y_ExtLatN, y_ExtCosLatN, q_LatN, q_CosLatN, xy_ExtFlagSlabOceanN, xyz_ExtSOSeaIceThicknessN, xyz_DSOSeaIceThicknessDtN )
    !
    call SOSIDiffusionY( y_ExtLatS, y_ExtCosLatS, q_LatS, q_CosLatS, xy_ExtFlagSlabOceanS, xyz_ExtSOSeaIceTempS, xyz_DSOSeaIceTempDtS )
    call SOSIDiffusionY( y_ExtLatN, y_ExtCosLatN, q_LatN, q_CosLatN, xy_ExtFlagSlabOceanN, xyz_ExtSOSeaIceTempN, xyz_DSOSeaIceTempDtN )
    !
    call SOSIDiffusionY( y_ExtLatS, y_ExtCosLatS, q_LatS, q_CosLatS, xy_ExtFlagSlabOceanS, xyz_ExtSurfTempS, xya_DSurfTempDtS, xy_ExtSOSeaIceMassS )
    call SOSIDiffusionY( y_ExtLatN, y_ExtCosLatN, q_LatN, q_CosLatN, xy_ExtFlagSlabOceanN, xyz_ExtSurfTempN, xya_DSurfTempDtN, xy_ExtSOSeaIceMassN )

    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
SOSeaIceDiffCoef
Variable :
SOSeaIceDiffCoef :real(DP) , save
module_name
Constant :
module_name = ‘sosi_dynamics :character(*), parameter
: モジュールの名称. Module name
p_LonN
Variable :
p_LonN(:) :real(DP) , save, allocatable
p_LonS
Variable :
p_LonS(:) :real(DP) , save, allocatable
q_CosLatN
Variable :
q_CosLatN(:) :real(DP) , save, allocatable
q_CosLatS
Variable :
q_CosLatS(:) :real(DP) , save, allocatable
q_LatN
Variable :
q_LatN(:) :real(DP) , save, allocatable
q_LatS
Variable :
q_LatS(:) :real(DP) , save, allocatable
sosi_dynamics_inited
Variable :
sosi_dynamics_inited = .false. :logical, save
: 初期設定フラグ. Initialization flag
version
Constant :
version = ’$Name: $’ // ’$Id: sltt.F90,v 1.8 2014/06/29 07:21:28 yot Exp $’ :character(*), parameter
: モジュールのバージョン Module version
x_CosLonN
Variable :
x_CosLonN(:) :real(DP) , save, allocatable
: $cos\lambda_N$
x_CosLonS
Variable :
x_CosLonS(:) :real(DP) , save, allocatable
: $cos\lambda_S$
x_ExtLonN
Variable :
x_ExtLonN(:) :real(DP) , save, allocatable
: $ x_LonNの拡張配列。 Extended array of x_LonN.
x_ExtLonS
Variable :
x_ExtLonS(:) :real(DP) , save, allocatable
: $ x_LonSの拡張配列。 Extended array of x_LonS.
x_LonN
Variable :
x_LonN(:) :real(DP) , save, allocatable
: $lambda_N$ 北半球の経度。 longitude in NH.
x_LonS
Variable :
x_LonS(:) :real(DP) , save, allocatable
: $lambda_S$ 南半球の経度。 longitude in SH.
x_SinLonN
Variable :
x_SinLonN(:) :real(DP) , save, allocatable
: $sin\lambda_N$
x_SinLonS
Variable :
x_SinLonS(:) :real(DP) , save, allocatable
: $sin\lambda_S$
y_CosLat
Variable :
y_CosLat(:) :real(DP) , save, allocatable
: $cos\varphai$
y_CosLatN
Variable :
y_CosLatN(:) :real(DP) , save, allocatable
: $cos\varphai_N$
y_CosLatS
Variable :
y_CosLatS(:) :real(DP) , save, allocatable
: $cos\varphai_S$
y_ExtCosLatN
Variable :
y_ExtCosLatN(:) :real(DP) , save, allocatable
: $ y_CosLatNの拡張配列。 Extended array of y_CosLatN.
y_ExtCosLatS
Variable :
y_ExtCosLatS(:) :real(DP) , save, allocatable
: $ y_CosLatS の拡張配列。 Extended array of y_CosLatS.
y_ExtLatN
Variable :
y_ExtLatN(:) :real(DP) , save, allocatable
: $ x_LatNの拡張配列。 Extended array of x_LatN.
y_ExtLatS
Variable :
y_ExtLatS(:) :real(DP) , save, allocatable
: $ x_LatSの拡張配列。 Extended array of x_LatS.
y_LatN
Variable :
y_LatN(:) :real(DP) , save, allocatable
: $varphi_N$ 北半球の緯度。 latitude in NH.
y_LatS
Variable :
y_LatS(:) :real(DP) , save, allocatable
: $varphi_S$ 南半球の緯度。 latitude in SH.
y_SinLatN
Variable :
y_SinLatN(:) :real(DP) , save, allocatable
: $sin\varphai_N$
y_SinLatS
Variable :
y_SinLatS(:) :real(DP) , save, allocatable
: $sin\varphai_S$