| Class | sosi_dynamics |
| In: |
sosi/sosi_dynamics.F90
|
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.
| !$ ! SLTTMain : | 移流計算 |
| !$ ! SLTTInit : | 初期化 |
| !$ ! SLTTTest : | 移流テスト用 |
| !$ ! ——————— : | ———— |
| !$ ! SLTTMain : | Main subroutine for SLTT |
| !$ ! SLTTInit : | Initialization for SLTT |
| !$ ! SLTTTest : | Generate velocity for SLTT Test |
NAMELIST#
モジュール引用 ; USE statements
種別型パラメタ Kind type parameter
| Subroutine : | |||
| xy_SurfType(0:imax-1, 1:jmax) : | integer , intent(in )
| ||
| xy_SurfTemp(0:imax-1, 1:jmax) : | real(DP), intent(inout) | ||
| xy_SOSeaIceMass(0:imax-1, 1:jmax) : | real(DP), intent(inout)
| ||
| 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 )
| ||
| xyz_DSOSeaIceTempDtPhy(0:imax-1, 1:jmax, 1:ksimax) : | real(DP), intent(in ) |
Calculates slab sea ice horizontal transports by diffusion
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 .
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
| 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 )
| ||
| xy_ExtSOSeaIceMassH(iexmin:iexmax, jexmin:jexmax) : | real(DP), intent(in )
| ||
| xy_DSOSeaIceMassDtH(0:imax-1, 1:jmax/2) : | real(DP), intent(out)
|
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
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 )
| ||
| xyz_ExtSOSeaIceMassH(iexmin:iexmax, jexmin:jexmax, 1:ksimax) : | real(DP), intent(in )
| ||
| xyz_DSOSeaIceMassDtH(0:imax-1, 1:jmax/2, 1:ksimax) : | real(DP), intent(out)
| ||
| xy_ExtSOSeaIceMassH(iexmin:iexmax, jexmin:jexmax) : | real(DP), intent(in ), optional |
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
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
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
| Constant : | |||
| module_name = ‘sosi_dynamics‘ : | character(*), parameter
|
| Variable : | |||
| sosi_dynamics_inited = .false. : | logical, save
|
| Variable : | |||
| y_ExtCosLatN(:) : | real(DP) , save, allocatable
|
| Variable : | |||
| y_ExtCosLatS(:) : | real(DP) , save, allocatable
|
| Variable : | |||
| y_ExtLatN(:) : | real(DP) , save, allocatable
|
| Variable : | |||
| y_ExtLatS(:) : | real(DP) , save, allocatable
|