!= ! != Slab ocean sea ice horizontal transport ! ! Authors:: Yoshiyuki O. TAKAHASHI ! Version:: ! Tag Name:: $Name: $ ! Copyright:: Copyright (C) GFD Dennou Club, 2013. All rights reserved. ! License:: See COPYRIGHT[link:../../../COPYRIGHT] ! module sosi_dynamics ! != ! != 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 ! use dc_types, only: DP, & ! 倍精度実数型. Double precision. & TOKEN ! キーワード. Keywords. ! メッセージ出力 ! Message output ! use dc_message, only: MessageNotify ! ! MPI ! use mpi_wrapper, only : MPIWrapperFindMaxVal ! 格子点設定 ! Grid points settings ! use gridset, only: & & imax, & ! 経度格子点数. ! Number of grid points in longitude & jmax, & ! 緯度格子点数. ! Number of grid points in latitude & kmax, & ! 鉛直層数. ! Number of vertical level & ksimax ! 海氷の鉛直層数. ! Number of sea ice vertical level ! 組成に関わる配列の設定 ! Settings of array for atmospheric composition ! use composition, only: & & ncmax ! 成分の数 ! Number of composition ! 質量の補正 ! Mass fixer ! !!$ use mass_fixer, only: & !!$ & MassFixerBC02, MassFixerBC02Layer, MassFixerBC02Column, & !!$ & MassFixer, MassFixerR95, MassFixerWO94, MassFixerColumn!, MassFixerLayer ! 宣言文 ; Declaration statements ! implicit none private ! 公開手続き ! Public procedure ! public :: SOSIDynamics public :: SOSIDynamicsInit ! 公開変数 ! Public variables ! ! 非公開変数 ! Private variables ! logical, save :: FlagSlabOcean ! flag for use of slab ocean real(DP) , save :: SOSeaIceDiffCoef real(DP) , save, allocatable :: y_CosLat(:) ! $\cos\varphai$ real(DP) , save, allocatable :: x_LonS (:) ! $\lambda_S$ 南半球の経度。 ! longitude in SH. real(DP) , save, allocatable :: x_SinLonS(:) ! $\sin\lambda_S$ real(DP) , save, allocatable :: x_CosLonS(:) ! $\cos\lambda_S$ real(DP) , save, allocatable :: y_LatS (:) ! $\varphi_S$ 南半球の緯度。 ! latitude in SH. real(DP) , save, allocatable :: y_SinLatS(:) ! $\sin\varphai_S$ real(DP) , save, allocatable :: y_CosLatS(:) ! $\cos\varphai_S$ real(DP) , save, allocatable :: x_ExtLonS(:) ! $ x_LonSの拡張配列。 !Extended array of x_LonS. real(DP) , save, allocatable :: y_ExtLatS(:) ! $ x_LatSの拡張配列。 !Extended array of x_LatS. real(DP) , save, allocatable :: y_ExtCosLatS(:) ! $ y_CosLatS の拡張配列。 !Extended array of y_CosLatS. real(DP) , save, allocatable :: x_LonN (:) ! $\lambda_N$ 北半球の経度。 ! longitude in NH. real(DP) , save, allocatable :: x_SinLonN(:) ! $\sin\lambda_N$ real(DP) , save, allocatable :: x_CosLonN(:) ! $\cos\lambda_N$ real(DP) , save, allocatable :: y_LatN (:) ! $\varphi_N$ 北半球の緯度。 ! latitude in NH. real(DP) , save, allocatable :: y_SinLatN(:) ! $\sin\varphai_N$ real(DP) , save, allocatable :: y_CosLatN(:) ! $\cos\varphai_N$ real(DP) , save, allocatable :: x_ExtLonN(:) ! $ x_LonNの拡張配列。 !Extended array of x_LonN. real(DP) , save, allocatable :: y_ExtLatN(:) ! $ x_LatNの拡張配列。 !Extended array of x_LatN. real(DP) , save, allocatable :: y_ExtCosLatN(:) ! $ y_CosLatNの拡張配列。 !Extended array of y_CosLatN. real(DP) , save, allocatable :: p_LonS (:) real(DP) , save, allocatable :: q_LatS (:) real(DP) , save, allocatable :: q_CosLatS(:) real(DP) , save, allocatable :: p_LonN (:) real(DP) , save, allocatable :: q_LatN (:) real(DP) , save, allocatable :: q_CosLatN(:) logical, save :: sosi_dynamics_inited = .false. ! 初期設定フラグ. ! Initialization flag character(*), parameter:: module_name = 'sosi_dynamics' ! モジュールの名称. ! Module name character(*), parameter:: version = & & '$Name: $' // & & '$Id: sltt.F90,v 1.8 2014/06/29 07:21:28 yot Exp $' ! モジュールのバージョン ! Module version !-------------------------------------------------------------------------------------- contains !-------------------------------------------------------------------------------------- subroutine SOSIDynamics( & & xy_SurfType, & !(in ) & xy_SurfTemp, xy_SOSeaIceMass, xyz_SOSeaIceTemp, & !(inout) & xy_DSOSeaIceMassDtPhyTop, xy_DSOSeaIceMassDtPhyBot, & !(in ) & xyz_DSOSeaIceTempDtPhy & !(in ) & ) ! ! Calculates slab sea ice horizontal transports by diffusion ! ヒストリデータ出力 ! History data output ! use gtool_historyauto, only: HistoryAutoPut use timeset , only : & & DelTime ! $\Delta t$ ! 物理・数学定数設定 ! Physical and mathematical constants settings ! use constants0, only: & & WaterHeatCap ! 物理定数設定 ! Physical constants settings ! use constants, only: & & SOMass ! Slab ocean mass ! 雪と海氷の定数の設定 ! Setting constants of snow and sea ice ! use constants_snowseaice, only: & & TempCondWater, & & SeaIceDen, & & SeaIceHeatCap, & & TempBelowSeaIce, & & SOSeaIceThresholdMass, & & LatentHeatFusionBelowSeaIce ! ! Slab ocean sea ice utility module ! use sosi_utils, only : & & SOSIUtilsChkSOSeaIce, & & SOSIUtilsSetSOSeaIceLevels, & & SOSeaIceMassNegativeThreshold, & & SOSIUtilsAddPhysics ! 宣言文 ; Declaration statements ! integer , intent(in ) :: xy_SurfType (0:imax-1, 1:jmax) ! 土地利用. ! Surface index real(DP), intent(inout) :: xy_SurfTemp (0:imax-1, 1:jmax) real(DP), intent(inout) :: xy_SOSeaIceMass (0:imax-1, 1:jmax) ! $ M_si $ . 海氷質量 (kg m-2) ! Slab ocean sea ice mass (kg m-2) real(DP), intent(inout) :: xyz_SOSeaIceTemp (0:imax-1, 1:jmax, 1:ksimax) real(DP), intent(in ) :: xy_DSOSeaIceMassDtPhyTop(0:imax-1, 1:jmax) real(DP), intent(in ) :: xy_DSOSeaIceMassDtPhyBot(0:imax-1, 1:jmax) ! ! Slab sea ice mass at next time step real(DP), intent(in ) :: xyz_DSOSeaIceTempDtPhy(0:imax-1, 1:jmax, 1:ksimax) ! 作業変数 ! Work variables ! real(DP) :: xy_SurfTempSave (0:imax-1, 1:jmax) real(DP) :: xy_SOSeaIceMassSave (0:imax-1, 1:jmax) real(DP) :: xyz_SOSeaIceTempSave(0:imax-1, 1:jmax, ksimax) real(DP) :: xy_SOSeaIceThickness(0:imax-1, 1:jmax) real(DP) :: xyz_SOSeaIceThickness(0:imax-1, 1:jmax, 1:ksimax) ! ! Sea ice thickness integer :: xy_SOSILocalKMax (0:imax-1, 1:jmax) real(DP) :: xyr_SOSILocalDepth(0:imax-1, 1:jmax, 0:ksimax) real(DP) :: xyz_SOSILocalDepth(0:imax-1, 1:jmax, 1:ksimax) !!$ real(DP) :: xy_SeaIceThicknessA(0:imax-1, 1:jmax) ! ! Sea ice thickness !!$ integer :: xy_SOSILocalKMaxA (0:imax-1, 1:jmax) !!$ real(DP) :: xyr_SOSILocalDepthA(0:imax-1, 1:jmax, 0:ksimax) !!$ real(DP) :: xyz_SOSILocalDepthA(0:imax-1, 1:jmax, 1:ksimax) logical :: xy_FlagSlabOcean(0:imax-1, 1:jmax) real(DP) :: xyz_DSOSeaIceThicknessDt(0:imax-1, 1:jmax, 1:ksimax) real(DP) :: xyz_DSOSeaIceTempDt (0:imax-1, 1:jmax, 1:ksimax) real(DP) :: xy_DSOTempDt (0:imax-1, 1:jmax) !!$ real(DP) :: xy_DSOSeaIceMassDt(0:imax-1, 1:jmax) !!$ ! !!$ ! Slab sea ice mass tendency real(DP) :: xy_SOTemp (0:imax-1, 1:jmax) real(DP) :: xyz_SOSeaIceThicknessAdv(0:imax-1, 1:jmax, 1:ksimax) real(DP) :: xyz_SOSeaIceTempAdv (0:imax-1, 1:jmax, 1:ksimax) real(DP) :: xy_SOTempAdv (0:imax-1, 1:jmax) !!$ real(DP) :: xy_TempSlabOcean (0:imax-1, 1:jmax) !!$ ! !!$ ! Slab ocean temperature !!$ real(DP) :: xy_TempSlabSeaIce(0:imax-1, 1:jmax) !!$ ! !!$ ! Slab sea ice temperature real(DP) :: xyz_SOSIMassEachLayer(0:imax-1, 1:jmax, 1:ksimax) real(DP) :: SOSIMass real(DP) :: SOSIMass1L real(DP) :: DelSOSIMass !!$ real(DP) :: SurfTempTent real(DP) :: SOTempTent real(DP) :: SOTempTent1st logical, parameter :: FlagSOSIAdv = .false. logical :: FlagCalc integer:: i ! 東西方向に回る DO ループ用作業変数 ! Work variables for DO loop in zonal direction integer:: j ! 南北方向に回る DO ループ用作業変数 ! Work variables for DO loop in meridional direction integer:: k if ( .not. sosi_dynamics_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if ! add physics tendency call SOSIUtilsAddPhysics( & & xy_SOSeaIceMass, xyz_SOSeaIceTemp, & !(inout) & xy_DSOSeaIceMassDtPhyTop, xy_DSOSeaIceMassDtPhyBot, & !(in ) & xyz_DSOSeaIceTempDtPhy & !(in ) & ) !!$ xy_SOSeaIceMassA = xy_SOSeaIceMassB & !!$ & + xy_DSOSeaIceMassDtPhy * ( 2.0_DP * DelTime ) !!$ !!$ !!$ ! !!$ ! Calcuate sea ice thickness !!$ ! !!$ xy_SeaIceThickness = xy_SOSeaIceMassB / SeaIceDen !!$ ! !!$ ! Set slab ocean sea ice levels !!$ ! !!$ call SOSIUtilsSetSOSeaIceLevels( & !!$ & xy_SeaIceThickness, & ! (in ) !!$ & xy_SOSILocalKMax, xyr_SOSILocalDepth, xyz_SOSILocalDepth & ! (out) !!$ & ) !!$ !!$ ! !!$ ! Calcuate sea ice thickness !!$ ! !!$ xy_SeaIceThicknessA = xy_SOSeaIceMassA / SeaIceDen !!$ ! !!$ ! Set slab ocean sea ice levels !!$ ! !!$ call SOSIUtilsSetSOSeaIceLevels( & !!$ & xy_SeaIceThicknessA, & ! (in ) !!$ & xy_SOSILocalKMaxA, xyr_SOSILocalDepthA, xyz_SOSILocalDepthA & ! (out) !!$ & ) !!$ !!$ !!$ ! 海氷温度時間積分 !!$ ! Time integration of sea ice temperature !!$ ! !!$ xyz_SOSeaIceTemp = xyz_SOSeaIceTemp + xyz_DSOSeaIceTempDtPhy * DelTime !!$ !!$ !!$ ! Adjust levels !!$ ! !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( xy_SOSILocalKMaxA(i,j) > xy_SOSILocalKMax(i,j) ) then !!$ ! sea ice thickness increases !!$ do k = 1, xy_SOSILocalKMax(i,j) !!$ xyz_SOSeaIceTemp(i,j,k) = xyz_SOSeaIceTemp(i,j,k) & !!$ & + xyz_DSOSeaIceTempDtPhy(i,j,k) * DelTime !!$ end do !!$ do k = xy_SOSILocalKMax(i,j)+1, xy_SOSILocalKMaxA(i,j) !!$ kk = xy_SOSILocalKMax(i,j) !!$ xyz_SOSeaIceTemp(i,j,k) = xyz_SOSeaIceTemp(i,j,kk) !!$ end do !!$ else if ( xy_SOSILocalKMaxA(i,j) < xy_SOSILocalKMax(i,j) ) then !!$ ! sea ice thickness decreases !!$ ! Do nothing !!$ ! Melted sea ice had freezing temperature !!$ end if !!$ end do !!$ end do xy_FlagSlabOcean = .false. FlagCalc = .false. if ( FlagSlabOcean ) then do j = 1, jmax do i = 0, imax-1 if ( xy_SurfType(i,j) <= 0 ) then ! slab ocean xy_FlagSlabOcean(i,j) = .true. FlagCalc = .true. end if end do end do end if if ( SOSeaIceDiffCoef <= 0.0_DP ) then FlagCalc = .false. !!$ else !!$ call MessageNotify( 'E', module_name, & !!$ & ' Now, SOSeaIceDiffCoef has to be zero.' ) end if if ( .not. FlagCalc ) then return end if ! Save values ! xy_SurfTempSave = xy_SurfTemp xy_SOSeaIceMassSave = xy_SOSeaIceMass xyz_SOSeaIceTempSave = xyz_SOSeaIceTemp ! ! Calcuate sea ice thickness ! xy_SOSeaIceThickness = xy_SOSeaIceMass / SeaIceDen ! ! Set slab ocean sea ice levels ! call SOSIUtilsSetSOSeaIceLevels( & & xy_SOSeaIceThickness, & ! (in ) & xy_SOSILocalKMax, xyr_SOSILocalDepth, xyz_SOSILocalDepth & ! (out) & ) ! Slab ocean temperature ! do j = 1, jmax do i = 0, imax-1 !!$ if ( xy_SOSILocalKMax(i,j) <= 0 ) then if ( xy_SOSeaIceMass(i,j) < SOSeaIceThresholdMass ) then xy_SOTemp(i,j) = xy_SurfTemp(i,j) else xy_SOTemp(i,j) = TempBelowSeaIce end if end do end do do j = 1, jmax do i = 0, imax-1 do k = 1, xy_SOSILocalKMax(i,j) xyz_SOSeaIceThicknessAdv(i,j,k) = & & xyr_SOSILocalDepth(i,j,k-1) - xyr_SOSILocalDepth(i,j,k) xyz_SOSeaIceTempAdv (i,j,k) = xyz_SOSeaIceTemp(i,j,k) end do do k = xy_SOSILocalKMax(i,j)+1, ksimax xyz_SOSeaIceThicknessAdv(i,j,k) = 0.0_DP xyz_SOSeaIceTempAdv (i,j,k) = TempCondWater !!$ xyz_SOSeaIceTempAdv (i,j,k) = TempBelowSeaIce end do end do end do ! xy_SOTempAdv = xy_SOTemp !!$ call SOSIUtilsSetMissingValue( & !!$ & xy_SOSeaIceMass, & !(in ) !!$ & xyz_SOSeaIceTemp, & !(inout) !!$ & SOSeaIceValue & !(in ) optional !!$ & ) call SOSIHorTransportDiff( & & xy_FlagSlabOcean, & ! (in) & xy_SOSeaIceMass, & ! (in) & xy_SOTempAdv, & ! (in) & xyz_SOSeaIceThicknessAdv, xyz_SOSeaIceTempAdv, & ! (in) & xy_DSOTempDt, & ! (out) & xyz_DSOSeaIceThicknessDt, xyz_DSOSeaIceTempDt & ! (out) & ) xyz_SOSeaIceThickness = xyz_SOSeaIceThicknessAdv xyz_SOSeaIceTemp = xyz_SOSeaIceTempAdv xy_SOTemp = xy_SOTempAdv ! if ( FlagSOSIAdv ) then xyz_SOSeaIceThickness = xyz_SOSeaIceThickness & & + xyz_DSOSeaIceThicknessDt * DelTime xyz_SOSeaIceTemp = xyz_SOSeaIceTemp & & + xyz_DSOSeaIceTempDt * DelTime else xy_SOTemp = xy_SOTemp & & + xy_DSOTempDt * DelTime end if ! Adjustment ! !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( xy_FlagSlabOcean(i,j) ) then !!$ do k = xy_SOSILocalKMax(i,j), 1, -1 !!$ if ( xyz_SOSeaIceThickness(i,j,k) > 0.0_DP ) then !!$ if ( xy_SurfTemp(i,j) > TempCondWater ) then !!$ ! Check whether all sea ice melt !!$ SOSIMass = SeaIceDen * xyz_SOSeaIceThickness(i,j,k) !!$ SurfTempTent = & !!$ & ( WaterHeatCap * ( SOMass - SOSIMass ) * xy_SurfTemp(i,j) & !!$ & + SeaIceHeatCap * SOSIMass * xyz_SOSeaIceTemp(i,j,k) & !!$ & - LatentHeatFusion * SOSIMass ) & !!$ & / ( WaterHeatCap * SOMass ) !!$ !!$ if ( SurfTempTent >= TempCondWater ) then !!$ ! Case in which all sea ice melt !!$ xy_SurfTemp (i,j) = SurfTempTent !!$ xyz_SOSeaIceThickness(i,j,k) = 0.0_DP !!$ xyz_SOSeaIceTemp (i,j,k) = TempCondWater !!$ else !!$ ! Case in which a part of sea ice melt !!$ DelSOSIMass = & !!$ & - WaterHeatCap * ( SOMass - SOSIMass ) & !!$ & * ( xy_SurfTemp(i,j) - TempBelowSeaIce ) & !!$ & / ( & !!$ & LatentHeatFusion & !!$ & + SeaIceHeatCap & !!$ & * ( TempBelowSeaIce - xyz_SOSeaIceTemp(i,j,k) ) & !!$ & ) !!$ !!$ xy_SurfTemp(i,j) = TempCondWater !!$ xyz_SOSeaIceThickness(i,j,k) = & !!$ & xyz_SOSeaIceThickness(i,j,k) & !!$ & + DelSOSIMass / SeaIceDen !!$ end if !!$ end if !!$ end if !!$ end do !!$ end if !!$ end do !!$ end do !!$ ! !!$ xy_SOSeaIceMass = 0.0_DP !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ do k = 1, xy_SOSILocalKMax(i,j) !!$ xy_SOSeaIceMass(i,j) = xy_SOSeaIceMass(i,j) & !!$ & + SeaIceDen * xyz_SOSeaIceThickness(i,j,k) !!$ end do !!$ end do !!$ end do do j = 1, jmax do i = 0, imax-1 do k = 1, xy_SOSILocalKMax(i,j) xyz_SOSIMassEachLayer(i,j,k) = SeaIceDen * xyz_SOSeaIceThickness(i,j,k) end do do k = xy_SOSILocalKMax(i,j)+1, ksimax xyz_SOSIMassEachLayer(i,j,k) = 0.0_DP end do end do end do xy_SOSeaIceMass = 0.0_DP do j = 1, jmax do i = 0, imax-1 do k = 1, xy_SOSILocalKMax(i,j) xy_SOSeaIceMass(i,j) = xy_SOSeaIceMass(i,j) & & + xyz_SOSIMassEachLayer(i,j,k) end do end do end do do j = 1, jmax do i = 0, imax-1 if ( xy_FlagSlabOcean(i,j) ) then do k = xy_SOSILocalKMax(i,j), 1, -1 if ( xy_SOTemp(i,j) > TempBelowSeaIce ) then ! Check whether all sea ice melt SOSIMass = xy_SOSeaIceMass(i,j) SOSIMass1L = xyz_SOSIMassEachLayer(i,j,k) DelSOSIMass = - SOSIMass1L !!$ SOTempTent = & !!$ & ( WaterHeatCap * ( SOMass - SOSIMass ) * xy_SOTemp(i,j) & !!$ & + SeaIceHeatCap * SOSIMass1L * xyz_SOSeaIceTemp(i,j,k) & !!$ & + LatentHeatFusion * DelSOSIMass ) & !!$ & / ( WaterHeatCap * ( ( SOMass - SOSIMass ) - DelSOSIMass ) ) SOTempTent1st = TempBelowSeaIce SOTempTent = & & ( WaterHeatCap * ( SOMass - SOSIMass ) * xy_SOTemp(i,j) & & - WaterHeatCap * DelSOSIMass * SOTempTent1st & & + SeaIceHeatCap * DelSOSIMass & & * ( SOTempTent1st - xyz_SOSeaIceTemp(i,j,k) ) & & + LatentHeatFusionBelowSeaice * DelSOSIMass ) & & / ( WaterHeatCap * ( ( SOMass - SOSIMass ) - DelSOSIMass ) ) if ( SOTempTent >= TempBelowSeaIce ) then ! Case in which all sea ice melt xy_SOTemp (i,j) = SOTempTent xyz_SOSeaIceTemp (i,j,k) = TempBelowSeaIce else ! Case in which a part of sea ice melt SOTempTent = TempBelowSeaIce !!$ DelSOSIMass = & !!$ & ( & !!$ & - WaterHeatCap * ( SOMass - SOSIMass ) & !!$ & * ( SOTempTent - xy_SOTemp(i,j) ) & !!$ & - SeaIceHeatCap * SOSIMass1L & !!$ & * ( xyz_SOSeaIceTemp(i,j,k) - xyz_SOSeaIceTemp(i,j,k) ) & !!$ & ) & !!$ & / ( & !!$ & - WaterHeatCap * SOTempTent & !!$ & + SeaIceHeatCap * xyz_SOSeaIceTemp(i,j,k) & !!$ & - LatentHeatFusion & !!$ & ) SOTempTent1st = TempBelowSeaIce DelSOSIMass = & & ( & & WaterHeatCap * ( SOMass - SOSIMass ) & & * ( SOTempTent - xy_SOTemp(i,j) ) & & ) & & / ( & & WaterHeatCap * ( SOTempTent - SOTempTent1st ) & & + SeaIceHeatCap * ( SOTempTent1st - xyz_SOSeaIceTemp(i,j,k) ) & & + LatentHeatFusionBelowSeaIce & & ) xy_SOTemp(i,j) = SOTempTent end if xyz_SOSIMassEachLayer(i,j,k) = & & xyz_SOSIMassEachLayer(i,j,k) + DelSOSIMass xy_SOSeaIceMass(i,j) = xy_SOSeaIceMass(i,j) + DelSOSIMass end if end do end if end do end do do j = 1, jmax do i = 0, imax-1 if ( xy_SOSeaIceMass(i,j) > SOSeaIceThresholdMass ) then ! in case with sea ice remained xy_SurfTemp(i,j) = xy_SurfTemp(i,j) else if ( xy_SOSeaIceMassSave(i,j) > SOSeaIceThresholdMass ) then ! in case with sea ice was present but melts completely xy_SurfTemp(i,j) = xy_SOTemp(i,j) else ! in case with sea ice was/is not present !!$ xy_SurfTemp(i,j) = xy_SurfTemp(i,j) xy_SurfTemp(i,j) = xy_SOTemp(i,j) end if end do end do do j = 1, jmax do i = 0, imax-1 if ( xy_FlagSlabOcean(i,j) ) then if ( xy_SOSeaIceMass(i,j) < 0 ) then if ( xy_SOSeaIceMass(i,j) < SOSeaIceMassNegativeThreshold ) then call MessageNotify( 'M', module_name, & & ' Slab sea ice mass is negative after diffusion, %f, and this is set to zero.', & & d = (/ xy_SOSeaIceMass(i,j) /) ) end if xy_SOSeaIceMass(i,j) = 0.0_DP end if end if end do end do ! Check ! xy_SOSeaIceThickness = xy_SOSeaIceMass / SeaIceDen ! call SOSIUtilsChkSOSeaIce( & & xy_SOSeaIceThickness, xyz_SOSeaIceTemp, & ! (in) & "SOSIDynamics" & ! (in) & ) end subroutine SOSIDynamics !---------------------------------------------------------------------------- subroutine SOSIHorTransportFFSL( & & xy_FlagSlabOcean, & ! (in) & xy_SOSeaIceMass, & ! (in) & xy_SurfTemp, & ! (in) & xyz_SOSeaIceThickness, xyz_SOSeaIceTemp, & ! (in) & xy_DSurfTempDt, & ! (out) & xyz_DSOSeaIceThicknessDt, xyz_DSOSeaIceTempDt & ! (out) & ) ! ! Calculates slab sea ice transport by diffusion use axesset , only : x_Lon ! $\lambda, \varphai$ lon and lat use sltt_const , only : dtjw, iexmin, iexmax, jexmin, jexmax use sosi_dynamics_extarr, only : SLTTExtArrExt2 ! 配列拡張ルーチン ! Expansion of arrays ! 雪と海氷の定数の設定 ! Setting constants of snow and sea ice ! use constants_snowseaice, only: & & SeaIceDen logical , intent(in ) :: xy_FlagSlabOcean(0:imax-1, 1:jmax) real(DP), intent(in ) :: xy_SOSeaIceMass (0:imax-1, 1:jmax) real(DP), intent(in ) :: xy_SurfTemp (0:imax-1, 1:jmax) real(DP), intent(in ) :: xyz_SOSeaIceThickness(0:imax-1, 1:jmax, 1:ksimax) real(DP), intent(in ) :: xyz_SOSeaIceTemp (0:imax-1, 1:jmax, 1:ksimax) real(DP), intent(out) :: xy_DSurfTempDt (0:imax-1, 1:jmax) real(DP), intent(out) :: xyz_DSOSeaIceThicknessDt(0:imax-1, 1:jmax, 1:ksimax) real(DP), intent(out) :: xyz_DSOSeaIceTempDt (0:imax-1, 1:jmax, 1:ksimax) ! ! local variables ! real(DP) :: xy_SurfTempA (0:imax-1, 1:jmax) real(DP) :: xyz_SOSeaIceThicknessA(0:imax-1, 1:jmax, 1:ksimax) real(DP) :: xyz_SOSeaIceTempA (0:imax-1, 1:jmax, 1:ksimax) real(DP) :: xy_SOSeaIceMassA (0:imax-1, 1:jmax) real(DP) :: xyza_TMP4DArray(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) real(DP) :: xyza_ExtTMP4DArrayS(iexmin:iexmax, jexmin:jexmax, 1:kmax, 1:ncmax) ! ! Extended 4D array (SH) real(DP) :: xyza_ExtTMP4DArrayN(iexmin:iexmax, jexmin:jexmax, 1:kmax, 1:ncmax) ! ! Extended 4D array (NH) real(DP) :: xyz_ExtSOSeaIceThicknessS(iexmin:iexmax, jexmin:jexmax, 1:ksimax) real(DP) :: xyz_ExtSOSeaIceTempS (iexmin:iexmax, jexmin:jexmax, 1:ksimax) real(DP) :: xyz_ExtSurfTempS (iexmin:iexmax, jexmin:jexmax, 1:ksimax) real(DP) :: xy_ExtSOSeaIceMassS (iexmin:iexmax, jexmin:jexmax) ! ! Extended 4D array (SH) real(DP) :: xyz_ExtSOSeaIceThicknessN(iexmin:iexmax, jexmin:jexmax, 1:ksimax) real(DP) :: xyz_ExtSOSeaIceTempN (iexmin:iexmax, jexmin:jexmax, 1:ksimax) real(DP) :: xyz_ExtSurfTempN (iexmin:iexmax, jexmin:jexmax, 1:ksimax) real(DP) :: xy_ExtSOSeaIceMassN (iexmin:iexmax, jexmin:jexmax) ! ! Extended 4D array (NH) logical :: xy_ExtFlagSlabOceanS(iexmin:iexmax, jexmin:jexmax) ! ! Extended 4D array (SH) logical :: xy_ExtFlagSlabOceanN(iexmin:iexmax, jexmin:jexmax) ! ! Extended 4D array (NH) real(DP) :: xyz_DSOSeaIceThicknessDtS(0:imax-1, 1:jmax/2, 1:ksimax) real(DP) :: xyz_DSOSeaIceTempDtS (0:imax-1, 1:jmax/2, 1:ksimax) real(DP) :: xya_DSurfTempDtS (0:imax-1, 1:jmax/2, 1:ksimax) real(DP) :: xyz_DSOSeaIceThicknessDtN(0:imax-1, 1:jmax/2, 1:ksimax) real(DP) :: xyz_DSOSeaIceTempDtN (0:imax-1, 1:jmax/2, 1:ksimax) real(DP) :: xya_DSurfTempDtN (0:imax-1, 1:jmax/2, 1:ksimax) real(DP) :: PM ! 配列拡張する際、極ごえ後に符号が変わる場合は -1.0を与える。 ! そうでない場合は1.0を与える。 ! Sign change flag for array extension; -1.0 for sign change ! over the pole, 1.0 for no sign change integer:: i ! 東西方向に回る DO ループ用作業変数 ! Work variables for DO loop in zonal direction integer:: j ! 南北方向に回る DO ループ用作業変数 ! Work variables for DO loop in meridional direction integer:: k ! 鉛直方向に回る DO ループ用作業変数 ! Work variables for DO loop in vertical direction integer:: n ! 組成方向に回る DO ループ用作業変数 ! Work variables for DO loop in dimension of constituents integer:: kk ! 実行文 ; Executable statement ! ! 初期化確認 ! Initialization check ! if ( .not. sosi_dynamics_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if ! ! Longitudinal diffusion ! #if defined(AXISYMMETRY) || defined(AXISYMMETRY_SJPACK) xy_SurfTempA = xy_SurfTemp xyz_SOSeaIceThicknessA = xyz_SOSeaIceThickness xyz_SOSeaIceTempA = xyz_SOSeaIceTemp #else call SOSIDiffusionX( & & x_Lon(1)-x_Lon(0), y_CosLat, xy_FlagSlabOcean, & ! (in) & xy_SOSeaIceMass, & ! (in) & xy_SurfTemp , xyz_SOSeaIceThickness , xyz_SOSeaIceTemp, & ! (in) & xy_SurfTempA, xyz_SOSeaIceThicknessA, xyz_SOSeaIceTempA & ! (out) & ) #endif xy_SOSeaIceMassA = 0.0_DP do k = 1, ksimax xy_SOSeaIceMassA = xy_SOSeaIceMassA & & + SeaIceDen * xyz_SOSeaIceThicknessA(:,:,k) end do ! ! Latitudinal diffusion ! ! 配列の分割と拡張 ! Division and extension of arrays ! if ( ksimax > kmax ) then call MessageNotify( 'E', module_name, 'ksimax > kmax.' ) end if if ( kmax < 3 ) then call MessageNotify( 'E', module_name, 'kmax < 3.' ) end if if ( ncmax < 3 ) then call MessageNotify( 'E', module_name, 'ncmax < 3.' ) end if n = 1 do k = 1, ksimax xyza_TMP4DArray(:,:,k,n) = xyz_SOSeaIceThicknessA(:,:,k) end do do k = ksimax+1, kmax xyza_TMP4DArray(:,:,k,n) = 0.0_DP end do ! n = 2 do k = 1, ksimax xyza_TMP4DArray(:,:,k,n) = xyz_SOSeaIceTempA(:,:,k) end do do k = ksimax+1, kmax xyza_TMP4DArray(:,:,k,n) = 0.0_DP end do ! n = 3 k = 1 do j = 1, jmax do i = 0, imax-1 if ( xy_FlagSlabOcean(i,j) ) then xyza_TMP4DArray(i,j,k,n) = 1.0_DP else xyza_TMP4DArray(i,j,k,n) = 0.0_DP end if end do end do k = 2 xyza_TMP4DArray(:,:,k,n) = xy_SurfTempA k = 3 xyza_TMP4DArray(:,:,k,n) = xy_SOSeaIceMassA do k = 3+1, kmax xyza_TMP4DArray(:,:,k,n) = 0.0_DP end do ! do n = 3+1, ncmax xyza_TMP4DArray(:,:,:,n) = 0.0_DP end do PM = 1.0_DP call SLTTExtArrExt2( & & x_SinLonS, x_CosLonS, x_SinLonN, x_CosLonN, & ! (in) & xyza_TMP4DArray, PM, & ! (in) & xyza_ExtTMP4DArrayS, xyza_ExtTMP4DArrayN & ! (out) & ) n = 1 do k = 1, ksimax xyz_ExtSOSeaIceThicknessS(:,:,k) = xyza_ExtTMP4DArrayS(:,:,k,n) xyz_ExtSOSeaIceThicknessN(:,:,k) = xyza_ExtTMP4DArrayN(:,:,k,n) end do n = 2 do k = 1, ksimax xyz_ExtSOSeaIceTempS (:,:,k) = xyza_ExtTMP4DArrayS(:,:,k,n) xyz_ExtSOSeaIceTempN (:,:,k) = xyza_ExtTMP4DArrayN(:,:,k,n) end do n = 3 k = 1 do j = jexmin, jexmax do i = iexmin, iexmax if ( xyza_ExtTMP4DArrayS(i,j,k,n) == 1.0_DP ) then xy_ExtFlagSlabOceanS(i,j) = .true. else xy_ExtFlagSlabOceanS(i,j) = .false. end if if ( xyza_ExtTMP4DArrayN(i,j,k,n) == 1.0_DP ) then xy_ExtFlagSlabOceanN(i,j) = .true. else xy_ExtFlagSlabOceanN(i,j) = .false. end if end do end do k = 2 do kk = 1, ksimax xyz_ExtSurfTempS(:,:,kk) = xyza_ExtTMP4DArrayS(:,:,k,n) xyz_ExtSurfTempN(:,:,kk) = xyza_ExtTMP4DArrayN(:,:,k,n) end do k = 3 do kk = 1, ksimax xy_ExtSOSeaIceMassS = xyza_ExtTMP4DArrayS(:,:,k,n) xy_ExtSOSeaIceMassN = xyza_ExtTMP4DArrayN(:,:,k,n) end do !!$ call SSIDiffusion( & !!$ & x_ExtLonS, y_ExtLatS, y_ExtCosLatS, & ! (in) !!$ & p_LonS, q_LatS, q_CosLatS, & ! (in) !!$ & xy_ExtFlagSlabOceanS, xy_ExtSOSeaIceMassS, & ! (in) !!$ & xy_DSOSeaIceMassDtS & ! (out) !!$ & ) !!$ call SSIDiffusion( & !!$ & x_ExtLonN, y_ExtLatN, y_ExtCosLatN, & ! (in) !!$ & p_LonN, q_LatN, q_CosLatN, & ! (in) !!$ & xy_ExtFlagSlabOceanN, xy_ExtSOSeaIceMassN, & ! (in) !!$ & xy_DSOSeaIceMassDtN & ! (out) !!$ & ) call SOSIDiffusionY( & & y_ExtLatS, y_ExtCosLatS, & ! (in) & q_LatS, q_CosLatS, & ! (in) & xy_ExtFlagSlabOceanS, xyz_ExtSOSeaIceThicknessS, & ! (in) & xyz_DSOSeaIceThicknessDtS & ! (out) & ) call SOSIDiffusionY( & & y_ExtLatN, y_ExtCosLatN, & ! (in) & q_LatN, q_CosLatN, & ! (in) & xy_ExtFlagSlabOceanN, xyz_ExtSOSeaIceThicknessN, & ! (in) & xyz_DSOSeaIceThicknessDtN & ! (out) & ) ! call SOSIDiffusionY( & & y_ExtLatS, y_ExtCosLatS, & ! (in) & q_LatS, q_CosLatS, & ! (in) & xy_ExtFlagSlabOceanS, xyz_ExtSOSeaIceTempS, & ! (in) & xyz_DSOSeaIceTempDtS & ! (out) & ) call SOSIDiffusionY( & & y_ExtLatN, y_ExtCosLatN, & ! (in) & q_LatN, q_CosLatN, & ! (in) & xy_ExtFlagSlabOceanN, xyz_ExtSOSeaIceTempN, & ! (in) & xyz_DSOSeaIceTempDtN & ! (out) & ) ! call SOSIDiffusionY( & & y_ExtLatS, y_ExtCosLatS, & ! (in) & q_LatS, q_CosLatS, & ! (in) & xy_ExtFlagSlabOceanS, xyz_ExtSurfTempS, & ! (in) & xya_DSurfTempDtS, & ! (out) & xy_ExtSOSeaIceMassS & ! (in) optional & ) call SOSIDiffusionY( & & y_ExtLatN, y_ExtCosLatN, & ! (in) & q_LatN, q_CosLatN, & ! (in) & xy_ExtFlagSlabOceanN, xyz_ExtSurfTempN, & ! (in) & xya_DSurfTempDtN, & ! (out) & xy_ExtSOSeaIceMassN & ! (in) optional & ) xyz_DSOSeaIceThicknessDt(:,1:jmax/2 ,:) = xyz_DSOSeaIceThicknessDtS(:,1:jmax/2,:) xyz_DSOSeaIceThicknessDt(:,jmax/2+1:jmax,:) = xyz_DSOSeaIceThicknessDtN(:,1:jmax/2,:) ! xyz_DSOSeaIceTempDt(:,1:jmax/2 ,:) = xyz_DSOSeaIceTempDtS(:,1:jmax/2,:) xyz_DSOSeaIceTempDt(:,jmax/2+1:jmax,:) = xyz_DSOSeaIceTempDtN(:,1:jmax/2,:) ! xy_DSurfTempDt(:,1:jmax/2 ) = xya_DSurfTempDtS(:,1:jmax/2,1) xy_DSurfTempDt(:,jmax/2+1:jmax) = xya_DSurfTempDtN(:,1:jmax/2,1) !!$ ! sea ice mass at next time step is calculated temporarily !!$ xy_SOSeaIceMassA = xy_SOSeaIceMassA & !!$ & + xy_DSOSeaIceMassDt * ( 2.0_DP * DelTime ) !!$ !!$ ! tendency is calculated !!$ xy_DSOSeaIceMassDt = & !!$ & ( xy_SOSeaIceMassA - xy_SOSeaIceMassB ) / ( 2.0_DP * DelTime ) !!$ py_ExtSOSeaIceMassXFlux = 0.0_DP !!$ xq_ExtSOSeaIceMassYFlux = 0.0_DP !!$ do j = jexmin, jexmax-1 !!$ do i = iexmin, iexmax-1 !!$ py_ExtSOSeaIceMassXFlux(i,j) = & !!$ & - SOSeaIceDiffCoef & !!$ & * ( xy_ExtSOSeaIceMassS(i+1,j) - xy_ExtSOSeaIceMassS(i,j) ) & !!$ & / ( RPlanet * y_CosLatS(j) * ( x_ExtLonS(i+1) - x_ExtLonS(i) ) ) !!$ xq_ExtSOSeaIceMassYFlux(i,j) = & !!$ & - SOSeaIceDiffCoef & !!$ & * ( xy_ExtSOSeaIceMassS(i,j+1) - xy_ExtSOSeaIceMassS(i,j) ) & !!$ & / ( RPlanet * ( y_ExtLatS(j+1) - y_ExtLatS(j) ) ) !!$ end do !!$ end do !!$ if ( myrank == nprocs-1 ) then !!$ j = 0 !!$ do i = iexmin, iexmax-1 !!$ xq_ExtSOSeaIceMassYFlux(i,j) = 0.0_DP !!$ end do !!$ end if !!$ !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ xy_DSOSeaIceMassDt(i,j) = & !!$ & - ( py_ExtSOSeaIceMassXFlux(i,j) - py_ExtSOSeaIceMassXFlux(i-1,j) ) & !!$ & / ( RPlanet * y_CosLatS(j) * ( p_Lon(i) - p_Lon(i-1) ) ) & !!$ & - ( xq_ExtSOSeaIceMassYFlux(i,j ) * q_CosLatS(j ) & !!$ & - xq_ExtSOSeaIceMassYFlux(i,j-1) * q_CosLatS(j-1) ) & !!$ & / ( RPlanet * y_CosLatS(j) * ( q_LatS(j) - q_LatS(j-1) ) ) !!$ end do !!$ end do end subroutine SOSIHorTransportFFSL !---------------------------------------------------------------------------- !!$ subroutine SOSIFFSLX( & !!$ & DelLon, y_CosLat, xy_FlagSlabOcean, & ! (in) !!$ & xy_SOSeaIceMass, & ! (in) !!$ & xy_SurfTemp , xyz_SOSeaIceThickness , xyz_SOSeaIceTemp, & ! (in) !!$ & xy_SurfTempA, xyz_SOSeaIceThicknessA, xyz_SOSeaIceTempA & ! (out) !!$ & ) !!$ ! !!$ ! Calculates slab sea ice transport by longitudinal diffusion !!$ !!$ !!$ ! !!$ ! !!$ use ludecomp_module, only : & !!$ & ludecomp_prep_simple_many, & !!$ & ludecomp_solve_simple_many !!$ !!$ use constants, only: & !!$ & RPlanet, & !!$ ! $ a $ [m]. !!$ ! 惑星半径. !!$ ! Radius of planet !!$ & SOMass !!$ ! Slab ocean mass !!$ !!$ real(DP), intent(in ) :: DelLon !!$ real(DP), intent(in ) :: y_CosLat(1:jmax) !!$ logical , intent(in ) :: xy_FlagSlabOcean(0:imax-1, 1:jmax) !!$ real(DP), intent(in ) :: xy_SOSeaIceMass (0:imax-1, 1:jmax) !!$ real(DP), intent(in ) :: xy_SurfTemp (0:imax-1, 1:jmax) !!$ real(DP), intent(in ) :: xyz_SOSeaIceThickness (0:imax-1, 1:jmax, 1:ksimax) !!$ real(DP), intent(in ) :: xyz_SOSeaIceTemp (0:imax-1, 1:jmax, 1:ksimax) !!$ real(DP), intent(out) :: xy_SurfTempA (0:imax-1, 1:jmax) !!$ real(DP), intent(out) :: xyz_SOSeaIceThicknessA(0:imax-1, 1:jmax, 1:ksimax) !!$ real(DP), intent(out) :: xyz_SOSeaIceTempA (0:imax-1, 1:jmax, 1:ksimax) !!$ !!$ ! !!$ ! local variables !!$ ! !!$ real(DP) :: SOSeaIceU0 = 0.1_DP !!$ !!$ real(DP) :: y_Factor(1:jmax) !!$ !!$ real(DP) :: py_SOSeaIceU(0:imax-1, 1:jmax) !!$ ! !!$ ! Longitudional Flux of slab sea ice !!$ !!$ integer:: i ! 東西方向に回る DO ループ用作業変数 !!$ ! Work variables for DO loop in zonal direction !!$ integer:: j ! 南北方向に回る DO ループ用作業変数 !!$ ! Work variables for DO loop in meridional direction !!$ integer:: k !!$ !!$ integer:: l !!$ integer:: ii !!$ integer:: iprev !!$ integer:: inext !!$ !!$ !!$ ! 実行文 ; Executable statement !!$ ! !!$ !!$ ! 初期化確認 !!$ ! Initialization check !!$ ! !!$ if ( .not. sosi_dynamics_inited ) then !!$ call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) !!$ end if !!$ !!$ !!$ y_Factor = 1.0_DP / ( RPlanet * y_CosLat * DelLon )**2 !!$ !!$ ! grid arrangement !!$ ! !!$ ! | . | . | . | . | . | . !!$ ! !!$ ! x 0 1 2 3 4 !!$ ! p 0 1 2 3 4 !!$ ! !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( i == imax-1 ) then !!$ iprev = i !!$ inext = 0 !!$ else !!$ iprev = i !!$ inext = i+1 !!$ end if !!$ if ( xy_FlagSlabOcean(iprev,j) .and. xy_FlagSlabOcean(inext,j) ) then !!$ py_SOSeaIceU(i,j) = SOSeaIceU0 !!$ else !!$ py_SOSeaIceU(i,j) = 0.0_DP !!$ end if !!$ end do !!$ end do !!$ !!$! do i = 0, imax-1 !!$! p_Lon(i) = DelLon / 2.0_DP + DelLon * i !!$! end do !!$ p_Lon = x_Lon + ( x_Lon(1) - x_Lon(0) ) / 2.0_DP !!$ !!$ ! Departure point !!$ do j = 1, jmax !!$ py_DPLon(:,j) = & !!$ & p_Lon - py_SOSeaIceU(:,j) * DelTime / ( RPlanet * y_CosLat(j) ) !!$ end do !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( py_DPLon(i,j) < 0.0_DP ) then !!$ py_DPLon(i,j) = py_DPLon(i,j) + 2.0_DP * PI !!$ else if ( py_DPLon(i,j) > 2.0_DP * PI ) then !!$ py_DPLon(i,j) = py_DPLon(i,j) - 2.0_DP * PI !!$ end if !!$ end do !!$ end do !!$ !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ dp_search : do ii = 0, imax-1 !!$ if ( py_DPLon(i,j) <= p_Lon(ii) ) then !!$ exit dp_search !!$ end if !!$ end do dp_search !!$ py_DPIndex !!$ if ( !!$ py_SOSeaIceDPLon(i,j) = p_Lon !!$ end do !!$ end do !!$ !!$ !!$ !!$ !!$ end subroutine SOSIFFSLX !---------------------------------------------------------------------------- subroutine SOSIHorTransportDiff( & & xy_FlagSlabOcean, & ! (in) & xy_SOSeaIceMass, & ! (in) & xy_SurfTemp, & ! (in) & xyz_SOSeaIceThickness, xyz_SOSeaIceTemp, & ! (in) & xy_DSurfTempDt, & ! (out) & xyz_DSOSeaIceThicknessDt, xyz_DSOSeaIceTempDt & ! (out) & ) ! ! Calculates slab sea ice transport by diffusion use axesset , only : x_Lon ! $\lambda, \varphai$ lon and lat use sltt_const , only : dtjw, iexmin, iexmax, jexmin, jexmax use sosi_dynamics_extarr, only : SLTTExtArrExt2 ! 配列拡張ルーチン ! Expansion of arrays ! 雪と海氷の定数の設定 ! Setting constants of snow and sea ice ! use constants_snowseaice, only: & & SeaIceDen logical , intent(in ) :: xy_FlagSlabOcean(0:imax-1, 1:jmax) real(DP), intent(in ) :: xy_SOSeaIceMass (0:imax-1, 1:jmax) real(DP), intent(in ) :: xy_SurfTemp (0:imax-1, 1:jmax) real(DP), intent(in ) :: xyz_SOSeaIceThickness(0:imax-1, 1:jmax, 1:ksimax) real(DP), intent(in ) :: xyz_SOSeaIceTemp (0:imax-1, 1:jmax, 1:ksimax) real(DP), intent(out) :: xy_DSurfTempDt (0:imax-1, 1:jmax) real(DP), intent(out) :: xyz_DSOSeaIceThicknessDt(0:imax-1, 1:jmax, 1:ksimax) real(DP), intent(out) :: xyz_DSOSeaIceTempDt (0:imax-1, 1:jmax, 1:ksimax) ! ! local variables ! real(DP) :: xy_SurfTempA (0:imax-1, 1:jmax) real(DP) :: xyz_SOSeaIceThicknessA(0:imax-1, 1:jmax, 1:ksimax) real(DP) :: xyz_SOSeaIceTempA (0:imax-1, 1:jmax, 1:ksimax) real(DP) :: xy_SOSeaIceMassA (0:imax-1, 1:jmax) real(DP) :: xyza_TMP4DArray(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) real(DP) :: xyza_ExtTMP4DArrayS(iexmin:iexmax, jexmin:jexmax, 1:kmax, 1:ncmax) ! ! Extended 4D array (SH) real(DP) :: xyza_ExtTMP4DArrayN(iexmin:iexmax, jexmin:jexmax, 1:kmax, 1:ncmax) ! ! Extended 4D array (NH) real(DP) :: xyz_ExtSOSeaIceThicknessS(iexmin:iexmax, jexmin:jexmax, 1:ksimax) real(DP) :: xyz_ExtSOSeaIceTempS (iexmin:iexmax, jexmin:jexmax, 1:ksimax) real(DP) :: xyz_ExtSurfTempS (iexmin:iexmax, jexmin:jexmax, 1:ksimax) real(DP) :: xy_ExtSOSeaIceMassS (iexmin:iexmax, jexmin:jexmax) ! ! Extended 4D array (SH) real(DP) :: xyz_ExtSOSeaIceThicknessN(iexmin:iexmax, jexmin:jexmax, 1:ksimax) real(DP) :: xyz_ExtSOSeaIceTempN (iexmin:iexmax, jexmin:jexmax, 1:ksimax) real(DP) :: xyz_ExtSurfTempN (iexmin:iexmax, jexmin:jexmax, 1:ksimax) real(DP) :: xy_ExtSOSeaIceMassN (iexmin:iexmax, jexmin:jexmax) ! ! Extended 4D array (NH) logical :: xy_ExtFlagSlabOceanS(iexmin:iexmax, jexmin:jexmax) ! ! Extended 4D array (SH) logical :: xy_ExtFlagSlabOceanN(iexmin:iexmax, jexmin:jexmax) ! ! Extended 4D array (NH) real(DP) :: xyz_DSOSeaIceThicknessDtS(0:imax-1, 1:jmax/2, 1:ksimax) real(DP) :: xyz_DSOSeaIceTempDtS (0:imax-1, 1:jmax/2, 1:ksimax) real(DP) :: xya_DSurfTempDtS (0:imax-1, 1:jmax/2, 1:ksimax) real(DP) :: xyz_DSOSeaIceThicknessDtN(0:imax-1, 1:jmax/2, 1:ksimax) real(DP) :: xyz_DSOSeaIceTempDtN (0:imax-1, 1:jmax/2, 1:ksimax) real(DP) :: xya_DSurfTempDtN (0:imax-1, 1:jmax/2, 1:ksimax) real(DP) :: PM ! 配列拡張する際、極ごえ後に符号が変わる場合は -1.0を与える。 ! そうでない場合は1.0を与える。 ! Sign change flag for array extension; -1.0 for sign change ! over the pole, 1.0 for no sign change integer:: i ! 東西方向に回る DO ループ用作業変数 ! Work variables for DO loop in zonal direction integer:: j ! 南北方向に回る DO ループ用作業変数 ! Work variables for DO loop in meridional direction integer:: k ! 鉛直方向に回る DO ループ用作業変数 ! Work variables for DO loop in vertical direction integer:: n ! 組成方向に回る DO ループ用作業変数 ! Work variables for DO loop in dimension of constituents integer:: kk ! 実行文 ; Executable statement ! ! 初期化確認 ! Initialization check ! if ( .not. sosi_dynamics_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if ! ! Longitudinal diffusion ! #if defined(AXISYMMETRY) || defined(AXISYMMETRY_SJPACK) xy_SurfTempA = xy_SurfTemp xyz_SOSeaIceThicknessA = xyz_SOSeaIceThickness xyz_SOSeaIceTempA = xyz_SOSeaIceTemp #else call SOSIDiffusionX( & & x_Lon(1)-x_Lon(0), y_CosLat, xy_FlagSlabOcean, & ! (in) & xy_SOSeaIceMass, & ! (in) & xy_SurfTemp , xyz_SOSeaIceThickness , xyz_SOSeaIceTemp, & ! (in) & xy_SurfTempA, xyz_SOSeaIceThicknessA, xyz_SOSeaIceTempA & ! (out) & ) #endif xy_SOSeaIceMassA = 0.0_DP do k = 1, ksimax xy_SOSeaIceMassA = xy_SOSeaIceMassA & & + SeaIceDen * xyz_SOSeaIceThicknessA(:,:,k) end do ! ! Latitudinal diffusion ! ! 配列の分割と拡張 ! Division and extension of arrays ! if ( ksimax > kmax ) then call MessageNotify( 'E', module_name, 'ksimax > kmax.' ) end if if ( kmax < 3 ) then call MessageNotify( 'E', module_name, 'kmax < 3.' ) end if if ( ncmax < 3 ) then call MessageNotify( 'E', module_name, 'ncmax < 3.' ) end if n = 1 do k = 1, ksimax xyza_TMP4DArray(:,:,k,n) = xyz_SOSeaIceThicknessA(:,:,k) end do do k = ksimax+1, kmax xyza_TMP4DArray(:,:,k,n) = 0.0_DP end do ! n = 2 do k = 1, ksimax xyza_TMP4DArray(:,:,k,n) = xyz_SOSeaIceTempA(:,:,k) end do do k = ksimax+1, kmax xyza_TMP4DArray(:,:,k,n) = 0.0_DP end do ! n = 3 k = 1 do j = 1, jmax do i = 0, imax-1 if ( xy_FlagSlabOcean(i,j) ) then xyza_TMP4DArray(i,j,k,n) = 1.0_DP else xyza_TMP4DArray(i,j,k,n) = 0.0_DP end if end do end do k = 2 xyza_TMP4DArray(:,:,k,n) = xy_SurfTempA k = 3 xyza_TMP4DArray(:,:,k,n) = xy_SOSeaIceMassA do k = 3+1, kmax xyza_TMP4DArray(:,:,k,n) = 0.0_DP end do ! do n = 3+1, ncmax xyza_TMP4DArray(:,:,:,n) = 0.0_DP end do PM = 1.0_DP call SLTTExtArrExt2( & & x_SinLonS, x_CosLonS, x_SinLonN, x_CosLonN, & ! (in) & xyza_TMP4DArray, PM, & ! (in) & xyza_ExtTMP4DArrayS, xyza_ExtTMP4DArrayN & ! (out) & ) n = 1 do k = 1, ksimax xyz_ExtSOSeaIceThicknessS(:,:,k) = xyza_ExtTMP4DArrayS(:,:,k,n) xyz_ExtSOSeaIceThicknessN(:,:,k) = xyza_ExtTMP4DArrayN(:,:,k,n) end do n = 2 do k = 1, ksimax xyz_ExtSOSeaIceTempS (:,:,k) = xyza_ExtTMP4DArrayS(:,:,k,n) xyz_ExtSOSeaIceTempN (:,:,k) = xyza_ExtTMP4DArrayN(:,:,k,n) end do n = 3 k = 1 do j = jexmin, jexmax do i = iexmin, iexmax if ( xyza_ExtTMP4DArrayS(i,j,k,n) == 1.0_DP ) then xy_ExtFlagSlabOceanS(i,j) = .true. else xy_ExtFlagSlabOceanS(i,j) = .false. end if if ( xyza_ExtTMP4DArrayN(i,j,k,n) == 1.0_DP ) then xy_ExtFlagSlabOceanN(i,j) = .true. else xy_ExtFlagSlabOceanN(i,j) = .false. end if end do end do k = 2 do kk = 1, ksimax xyz_ExtSurfTempS(:,:,kk) = xyza_ExtTMP4DArrayS(:,:,k,n) xyz_ExtSurfTempN(:,:,kk) = xyza_ExtTMP4DArrayN(:,:,k,n) end do k = 3 do kk = 1, ksimax xy_ExtSOSeaIceMassS = xyza_ExtTMP4DArrayS(:,:,k,n) xy_ExtSOSeaIceMassN = xyza_ExtTMP4DArrayN(:,:,k,n) end do !!$ call SSIDiffusion( & !!$ & x_ExtLonS, y_ExtLatS, y_ExtCosLatS, & ! (in) !!$ & p_LonS, q_LatS, q_CosLatS, & ! (in) !!$ & xy_ExtFlagSlabOceanS, xy_ExtSOSeaIceMassS, & ! (in) !!$ & xy_DSOSeaIceMassDtS & ! (out) !!$ & ) !!$ call SSIDiffusion( & !!$ & x_ExtLonN, y_ExtLatN, y_ExtCosLatN, & ! (in) !!$ & p_LonN, q_LatN, q_CosLatN, & ! (in) !!$ & xy_ExtFlagSlabOceanN, xy_ExtSOSeaIceMassN, & ! (in) !!$ & xy_DSOSeaIceMassDtN & ! (out) !!$ & ) call SOSIDiffusionY( & & y_ExtLatS, y_ExtCosLatS, & ! (in) & q_LatS, q_CosLatS, & ! (in) & xy_ExtFlagSlabOceanS, xyz_ExtSOSeaIceThicknessS, & ! (in) & xyz_DSOSeaIceThicknessDtS & ! (out) & ) call SOSIDiffusionY( & & y_ExtLatN, y_ExtCosLatN, & ! (in) & q_LatN, q_CosLatN, & ! (in) & xy_ExtFlagSlabOceanN, xyz_ExtSOSeaIceThicknessN, & ! (in) & xyz_DSOSeaIceThicknessDtN & ! (out) & ) ! call SOSIDiffusionY( & & y_ExtLatS, y_ExtCosLatS, & ! (in) & q_LatS, q_CosLatS, & ! (in) & xy_ExtFlagSlabOceanS, xyz_ExtSOSeaIceTempS, & ! (in) & xyz_DSOSeaIceTempDtS & ! (out) & ) call SOSIDiffusionY( & & y_ExtLatN, y_ExtCosLatN, & ! (in) & q_LatN, q_CosLatN, & ! (in) & xy_ExtFlagSlabOceanN, xyz_ExtSOSeaIceTempN, & ! (in) & xyz_DSOSeaIceTempDtN & ! (out) & ) ! call SOSIDiffusionY( & & y_ExtLatS, y_ExtCosLatS, & ! (in) & q_LatS, q_CosLatS, & ! (in) & xy_ExtFlagSlabOceanS, xyz_ExtSurfTempS, & ! (in) & xya_DSurfTempDtS, & ! (out) & xy_ExtSOSeaIceMassS & ! (in) optional & ) call SOSIDiffusionY( & & y_ExtLatN, y_ExtCosLatN, & ! (in) & q_LatN, q_CosLatN, & ! (in) & xy_ExtFlagSlabOceanN, xyz_ExtSurfTempN, & ! (in) & xya_DSurfTempDtN, & ! (out) & xy_ExtSOSeaIceMassN & ! (in) optional & ) xyz_DSOSeaIceThicknessDt(:,1:jmax/2 ,:) = xyz_DSOSeaIceThicknessDtS(:,1:jmax/2,:) xyz_DSOSeaIceThicknessDt(:,jmax/2+1:jmax,:) = xyz_DSOSeaIceThicknessDtN(:,1:jmax/2,:) ! xyz_DSOSeaIceTempDt(:,1:jmax/2 ,:) = xyz_DSOSeaIceTempDtS(:,1:jmax/2,:) xyz_DSOSeaIceTempDt(:,jmax/2+1:jmax,:) = xyz_DSOSeaIceTempDtN(:,1:jmax/2,:) ! xy_DSurfTempDt(:,1:jmax/2 ) = xya_DSurfTempDtS(:,1:jmax/2,1) xy_DSurfTempDt(:,jmax/2+1:jmax) = xya_DSurfTempDtN(:,1:jmax/2,1) !!$ ! sea ice mass at next time step is calculated temporarily !!$ xy_SOSeaIceMassA = xy_SOSeaIceMassA & !!$ & + xy_DSOSeaIceMassDt * ( 2.0_DP * DelTime ) !!$ !!$ ! tendency is calculated !!$ xy_DSOSeaIceMassDt = & !!$ & ( xy_SOSeaIceMassA - xy_SOSeaIceMassB ) / ( 2.0_DP * DelTime ) !!$ py_ExtSOSeaIceMassXFlux = 0.0_DP !!$ xq_ExtSOSeaIceMassYFlux = 0.0_DP !!$ do j = jexmin, jexmax-1 !!$ do i = iexmin, iexmax-1 !!$ py_ExtSOSeaIceMassXFlux(i,j) = & !!$ & - SOSeaIceDiffCoef & !!$ & * ( xy_ExtSOSeaIceMassS(i+1,j) - xy_ExtSOSeaIceMassS(i,j) ) & !!$ & / ( RPlanet * y_CosLatS(j) * ( x_ExtLonS(i+1) - x_ExtLonS(i) ) ) !!$ xq_ExtSOSeaIceMassYFlux(i,j) = & !!$ & - SOSeaIceDiffCoef & !!$ & * ( xy_ExtSOSeaIceMassS(i,j+1) - xy_ExtSOSeaIceMassS(i,j) ) & !!$ & / ( RPlanet * ( y_ExtLatS(j+1) - y_ExtLatS(j) ) ) !!$ end do !!$ end do !!$ if ( myrank == nprocs-1 ) then !!$ j = 0 !!$ do i = iexmin, iexmax-1 !!$ xq_ExtSOSeaIceMassYFlux(i,j) = 0.0_DP !!$ end do !!$ end if !!$ !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ xy_DSOSeaIceMassDt(i,j) = & !!$ & - ( py_ExtSOSeaIceMassXFlux(i,j) - py_ExtSOSeaIceMassXFlux(i-1,j) ) & !!$ & / ( RPlanet * y_CosLatS(j) * ( p_Lon(i) - p_Lon(i-1) ) ) & !!$ & - ( xq_ExtSOSeaIceMassYFlux(i,j ) * q_CosLatS(j ) & !!$ & - xq_ExtSOSeaIceMassYFlux(i,j-1) * q_CosLatS(j-1) ) & !!$ & / ( RPlanet * y_CosLatS(j) * ( q_LatS(j) - q_LatS(j-1) ) ) !!$ end do !!$ end do end subroutine SOSIHorTransportDiff !---------------------------------------------------------------------------- subroutine SOSIDiffusion( & & x_ExtLonH, y_ExtLatH, y_ExtCosLatH, & ! (in) & p_LonH, q_LatH, q_CosLatH, & ! (in) & xy_ExtFlagSlabOceanH, xy_ExtSOSeaIceMassH, & ! (in) & xy_DSOSeaIceMassDtH & ! (out) & ) ! ! Calculates slab sea ice transport by diffusion use sltt_const , only : iexmin, iexmax, jexmin, jexmax ! 配列拡張ルーチン ! Expansion of arrays use constants, only: & & RPlanet ! $ a $ [m]. ! 惑星半径. ! Radius of planet real(DP), intent(in ) :: x_ExtLonH (iexmin:iexmax) real(DP), intent(in ) :: y_ExtLatH (jexmin:jexmax) real(DP), intent(in ) :: y_ExtCosLatH(jexmin:jexmax) real(DP), intent(in ) :: p_LonH (0-1:imax-1+1) real(DP), intent(in ) :: q_LatH (1-1:jmax/2+1) real(DP), intent(in ) :: q_CosLatH(1-1:jmax/2+1) real(DP), intent(in ) :: xy_ExtSOSeaIceMassH(iexmin:iexmax, jexmin:jexmax) ! ! Extended 4D array logical , intent(in ) :: xy_ExtFlagSlabOceanH(iexmin:iexmax, jexmin:jexmax) ! ! Extended 4D array real(DP), intent(out) :: xy_DSOSeaIceMassDtH(0:imax-1, 1:jmax/2) ! ! Slab sea ice mass tendency ! ! local variables ! real(DP) :: py_ExtSOSeaIceMassXFlux(iexmin-1:iexmax, jexmin:jexmax) ! ! Longitudional Flux of slab sea ice real(DP) :: xq_ExtSOSeaIceMassYFlux(iexmin:iexmax, jexmin-1:jexmax) ! ! Latitudinal Flux of slab sea ice integer:: i ! 東西方向に回る DO ループ用作業変数 ! Work variables for DO loop in zonal direction integer:: j ! 南北方向に回る DO ループ用作業変数 ! Work variables for DO loop in meridional direction ! 実行文 ; Executable statement ! ! 初期化確認 ! Initialization check ! if ( .not. sosi_dynamics_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if py_ExtSOSeaIceMassXFlux = 0.0_DP xq_ExtSOSeaIceMassYFlux = 0.0_DP do j = jexmin, jexmax-1 do i = iexmin, iexmax-1 if ( xy_ExtFlagSlabOceanH(i,j) .and. xy_ExtFlagSlabOceanH(i+1,j) ) then py_ExtSOSeaIceMassXFlux(i,j) = & & - SOSeaIceDiffCoef & & * ( xy_ExtSOSeaIceMassH(i+1,j) - xy_ExtSOSeaIceMassH(i,j) ) & & / ( RPlanet * y_ExtCosLatH(j) * ( x_ExtLonH(i+1) - x_ExtLonH(i) ) ) else py_ExtSOSeaIceMassXFlux(i,j) = 0.0_DP end if if ( xy_ExtFlagSlabOceanH(i,j) .and. xy_ExtFlagSlabOceanH(i,j+1) ) then xq_ExtSOSeaIceMassYFlux(i,j) = & & - SOSeaIceDiffCoef & & * ( xy_ExtSOSeaIceMassH(i,j+1) - xy_ExtSOSeaIceMassH(i,j) ) & & / ( RPlanet * ( y_ExtLatH(j+1) - y_ExtLatH(j) ) ) else xq_ExtSOSeaIceMassYFlux(i,j) = 0.0_DP end if end do end do !!$ if ( myrank == nprocs-1 ) then !!$ j = 0 !!$ do i = iexmin, iexmax-1 !!$ xq_ExtSOSeaIceMassYFlux(i,j) = 0.0_DP !!$ end do !!$ end if do j = 1, jmax/2 do i = 0, imax-1 xy_DSOSeaIceMassDtH(i,j) = & & - ( py_ExtSOSeaIceMassXFlux(i ,j) & & - py_ExtSOSeaIceMassXFlux(i-1,j) ) & & / ( RPlanet * y_ExtCosLatH(j) * ( p_LonH(i) - p_LonH(i-1) ) ) & & - ( xq_ExtSOSeaIceMassYFlux(i,j ) * q_CosLatH(j ) & & - xq_ExtSOSeaIceMassYFlux(i,j-1) * q_CosLatH(j-1) ) & & / ( RPlanet * y_ExtCosLatH(j) * ( q_LatH(j) - q_LatH(j-1) ) ) end do end do end subroutine SOSIDiffusion !---------------------------------------------------------------------------- subroutine SOSIDiffusionY( & & y_ExtLatH, y_ExtCosLatH, & ! (in) & q_LatH, q_CosLatH, & ! (in) & xy_ExtFlagSlabOceanH, xyz_ExtSOSeaIceMassH, & ! (in) & xyz_DSOSeaIceMassDtH, & ! (out) & xy_ExtSOSeaIceMassH & ! (in) optional & ) ! ! Calculates slab sea ice transport by diffusion use sltt_const , only : iexmin, iexmax, jexmin, jexmax ! 配列拡張ルーチン ! Expansion of arrays use constants, only: & & RPlanet, & ! $ a $ [m]. ! 惑星半径. ! Radius of planet & SOMass ! Slab ocean mass real(DP), intent(in ) :: y_ExtLatH (jexmin:jexmax) real(DP), intent(in ) :: y_ExtCosLatH(jexmin:jexmax) real(DP), intent(in ) :: q_LatH (1-1:jmax/2+1) real(DP), intent(in ) :: q_CosLatH(1-1:jmax/2+1) logical , intent(in ) :: xy_ExtFlagSlabOceanH(iexmin:iexmax, jexmin:jexmax) ! ! Extended 4D array real(DP), intent(in ) :: xyz_ExtSOSeaIceMassH(iexmin:iexmax, jexmin:jexmax, 1:ksimax) ! ! Extended 4D array real(DP), intent(out) :: xyz_DSOSeaIceMassDtH(0:imax-1, 1:jmax/2, 1:ksimax) ! ! Slab sea ice mass tendency real(DP), intent(in ), optional :: xy_ExtSOSeaIceMassH(iexmin:iexmax, jexmin:jexmax) ! ! local variables ! real(DP) :: xqz_ExtSODiffCoef (iexmin:iexmax, jexmin-1:jexmax, 1:ksimax) real(DP) :: xqz_ExtSOSeaIceMassYFlux(iexmin:iexmax, jexmin-1:jexmax, 1:ksimax) ! ! Latitudinal Flux of slab sea ice integer:: i ! 東西方向に回る DO ループ用作業変数 ! Work variables for DO loop in zonal direction integer:: j ! 南北方向に回る DO ループ用作業変数 ! Work variables for DO loop in meridional direction integer:: k ! 実行文 ; Executable statement ! ! 初期化確認 ! Initialization check ! if ( .not. sosi_dynamics_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if do k = 1, ksimax do j = jexmin, jexmax-1 do i = iexmin, iexmax-1 if ( xy_ExtFlagSlabOceanH(i,j) .and. xy_ExtFlagSlabOceanH(i,j+1) ) then xqz_ExtSODiffCoef(i,j,k) = SOSeaIceDiffCoef if ( present( xy_ExtSOSeaIceMassH ) ) then xqz_ExtSODiffCoef(i,j,k) = xqz_ExtSODiffCoef(i,j,k) & & * min( SOMass - xy_ExtSOSeaIceMassH(i,j ), & & SOMass - xy_ExtSOSeaIceMassH(i,j+1) ) end if else xqz_ExtSODiffCoef(i,j,k) = 0.0_DP end if end do end do end do do k = 1, ksimax do j = jexmin, jexmax-1 do i = iexmin, iexmax-1 xqz_ExtSOSeaIceMassYFlux(i,j,k) = & & - xqz_ExtSODiffCoef(i,j,k) & & * ( xyz_ExtSOSeaIceMassH(i,j+1,k) - xyz_ExtSOSeaIceMassH(i,j,k) ) & & / ( RPlanet * ( y_ExtLatH(j+1) - y_ExtLatH(j) ) ) end do end do end do !!$ if ( myrank == nprocs-1 ) then !!$ j = 0 !!$ do i = iexmin, iexmax-1 !!$ xq_ExtSOSeaIceMassYFlux(i,j) = 0.0_DP !!$ end do !!$ end if do k = 1, ksimax do j = 1, jmax/2 do i = 0, imax-1 xyz_DSOSeaIceMassDtH(i,j,k) = & & - ( xqz_ExtSOSeaIceMassYFlux(i,j ,k) * q_CosLatH(j ) & & - xqz_ExtSOSeaIceMassYFlux(i,j-1,k) * q_CosLatH(j-1) ) & & / ( RPlanet * y_ExtCosLatH(j) * ( q_LatH(j) - q_LatH(j-1) ) ) if ( present( xy_ExtSOSeaIceMassH ) ) then if ( SOMass - xy_ExtSOSeaIceMassH(i,j) > 0.0_DP ) then xyz_DSOSeaIceMassDtH(i,j,k) = xyz_DSOSeaIceMassDtH(i,j,k) & & / ( SOMass - xy_ExtSOSeaIceMassH(i,j) ) else xyz_DSOSeaIceMassDtH(i,j,k) = 0.0_DP end if end if end do end do end do end subroutine SOSIDiffusionY !---------------------------------------------------------------------------- subroutine SOSIDiffusionX( & & DelLon, y_CosLat, xy_FlagSlabOcean, & ! (in) & xy_SOSeaIceMass, & ! (in) & xy_SurfTemp , xyz_SOSeaIceThickness , xyz_SOSeaIceTemp, & ! (in) & xy_SurfTempA, xyz_SOSeaIceThicknessA, xyz_SOSeaIceTempA & ! (out) & ) ! ! Calculates slab sea ice transport by longitudinal diffusion ! ! use ludecomp_module, only : & & ludecomp_prep_simple_many, & & ludecomp_solve_simple_many use timeset , only : DelTime ! $\Delta t$ use constants, only: & & RPlanet, & ! $ a $ [m]. ! 惑星半径. ! Radius of planet & SOMass ! Slab ocean mass real(DP), intent(in ) :: DelLon real(DP), intent(in ) :: y_CosLat(1:jmax) logical , intent(in ) :: xy_FlagSlabOcean(0:imax-1, 1:jmax) real(DP), intent(in ) :: xy_SOSeaIceMass (0:imax-1, 1:jmax) real(DP), intent(in ) :: xy_SurfTemp (0:imax-1, 1:jmax) real(DP), intent(in ) :: xyz_SOSeaIceThickness (0:imax-1, 1:jmax, 1:ksimax) real(DP), intent(in ) :: xyz_SOSeaIceTemp (0:imax-1, 1:jmax, 1:ksimax) real(DP), intent(out) :: xy_SurfTempA (0:imax-1, 1:jmax) real(DP), intent(out) :: xyz_SOSeaIceThicknessA(0:imax-1, 1:jmax, 1:ksimax) real(DP), intent(out) :: xyz_SOSeaIceTempA (0:imax-1, 1:jmax, 1:ksimax) ! ! local variables ! real(DP) :: aax_LUMat(1:jmax*ksimax, 0:imax-1, 0:imax-1) real(DP) :: aa_LUVec (1:jmax*ksimax, 0:imax-1) real(DP) :: y_Factor(1:jmax) real(DP) :: pyz_SOSeaIceDiffCoef(0:imax-1, 1:jmax, 1:ksimax) ! ! Longitudional Flux of slab sea ice integer:: i ! 東西方向に回る DO ループ用作業変数 ! Work variables for DO loop in zonal direction integer:: j ! 南北方向に回る DO ループ用作業変数 ! Work variables for DO loop in meridional direction integer:: k integer:: l integer:: ii integer:: iprev integer:: inext ! 実行文 ; Executable statement ! ! 初期化確認 ! Initialization check ! if ( .not. sosi_dynamics_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if y_Factor = 1.0_DP / ( RPlanet * y_CosLat * DelLon )**2 do k = 1, ksimax do j = 1, jmax do i = 0, imax-1 if ( i == imax-1 ) then iprev = i inext = 0 else iprev = i inext = i+1 end if if ( xy_FlagSlabOcean(iprev,j) .and. xy_FlagSlabOcean(inext,j) ) then pyz_SOSeaIceDiffCoef(i,j,k) = SOSeaIceDiffCoef else pyz_SOSeaIceDiffCoef(i,j,k) = 0.0_DP end if end do end do end do aax_LUMat = 0.0_DP do k = 1, ksimax do j = 1, jmax l = (k-1)*jmax+j do ii = 0, 0 i = ii aax_LUMat(l,ii,imax-1) = & & - pyz_SOSeaIceDiffCoef(imax-1,j,k) & & * y_Factor(j) aax_LUMat(l,ii,i ) = & & 1.0_DP / ( 1.0_DP * DelTime ) & & + ( pyz_SOSeaIceDiffCoef(imax-1,j,k) + pyz_SOSeaIceDiffCoef(i,j,k) ) & & * y_Factor(j) aax_LUMat(l,ii,i+1) = & & - pyz_SOSeaIceDiffCoef(i ,j,k) & & * y_Factor(j) end do do ii = 0+1, (imax-1)-1 i = ii aax_LUMat(l,ii,i-1) = & & - pyz_SOSeaIceDiffCoef(i-1,j,k) & & * y_Factor(j) aax_LUMat(l,ii,i ) = & & 1.0_DP / ( 1.0_DP * DelTime ) & & + ( pyz_SOSeaIceDiffCoef(i-1,j,k) + pyz_SOSeaIceDiffCoef(i,j,k) ) & & * y_Factor(j) aax_LUMat(l,ii,i+1) = & & - pyz_SOSeaIceDiffCoef(i ,j,k) & & * y_Factor(j) end do do ii = imax-1, imax-1 i = ii aax_LUMat(l,ii,i-1) = & & - pyz_SOSeaIceDiffCoef(i-1,j,k) & & * y_Factor(j) aax_LUMat(l,ii,i ) = & & 1.0_DP / ( 1.0_DP * DelTime ) & & + ( pyz_SOSeaIceDiffCoef(i-1,j,k) + pyz_SOSeaIceDiffCoef(i,j,k) ) & & * y_Factor(j) aax_LUMat(l,ii,0 ) = & & - pyz_SOSeaIceDiffCoef(imax-1,j,k) & & * y_Factor(j) end do end do end do call ludecomp_prep_simple_many ( jmax*ksimax, imax, aax_LUMat, 1, jmax*ksimax ) do k = 1, ksimax do j = 1, jmax l = (k-1)*jmax+j do ii = 0, imax-1 i = ii aa_LUVec(l,ii) = xyz_SOSeaIceThickness(i,j,k) / ( 1.0_DP * DelTime ) end do end do end do call ludecomp_solve_simple_many( jmax*ksimax, imax, aax_LUMat, aa_LUVec, 1, jmax*ksimax ) do k = 1, ksimax do j = 1, jmax l = (k-1)*jmax+j do ii = 0, imax-1 i = ii xyz_SOSeaIceThicknessA(i,j,k) = aa_LUVec(l,ii) end do end do end do do k = 1, ksimax do j = 1, jmax l = (k-1)*jmax+j do ii = 0, imax-1 i = ii aa_LUVec(l,ii) = xyz_SOSeaIceTemp(i,j,k) / ( 1.0_DP * DelTime ) end do end do end do call ludecomp_solve_simple_many( jmax*ksimax, imax, aax_LUMat, aa_LUVec, 1, jmax*ksimax ) do k = 1, ksimax do j = 1, jmax l = (k-1)*jmax+j do ii = 0, imax-1 i = ii xyz_SOSeaIceTempA(i,j,k) = aa_LUVec(l,ii) end do end do end do ! Diffusion of slab ocean temperature ! do k = 1, ksimax do j = 1, jmax do i = 0, imax-1 if ( i == imax-1 ) then iprev = i inext = 0 else iprev = i inext = i+1 end if if ( xy_FlagSlabOcean(iprev,j) .and. xy_FlagSlabOcean(inext,j) ) then pyz_SOSeaIceDiffCoef(i,j,k) = SOSeaIceDiffCoef & & * min( SOMass - xy_SOSeaIceMass(iprev,j), & & SOMass - xy_SOSeaIceMass(inext,j) ) else pyz_SOSeaIceDiffCoef(i,j,k) = 0.0_DP end if end do end do end do aax_LUMat = 0.0_DP do k = 1, ksimax do j = 1, jmax l = (k-1)*jmax+j do ii = 0, 0 i = ii aax_LUMat(l,ii,imax-1) = & & - pyz_SOSeaIceDiffCoef(imax-1,j,k) & & * y_Factor(j) aax_LUMat(l,ii,i ) = & & 1.0_DP / ( 1.0_DP * DelTime ) & & * ( SOMass - xy_SOSeaIceMass(i,j) ) & & + ( pyz_SOSeaIceDiffCoef(imax-1,j,k) + pyz_SOSeaIceDiffCoef(i,j,k) ) & & * y_Factor(j) aax_LUMat(l,ii,i+1) = & & - pyz_SOSeaIceDiffCoef(i ,j,k) & & * y_Factor(j) end do do ii = 0+1, (imax-1)-1 i = ii aax_LUMat(l,ii,i-1) = & & - pyz_SOSeaIceDiffCoef(i-1,j,k) & & * y_Factor(j) aax_LUMat(l,ii,i ) = & & 1.0_DP / ( 1.0_DP * DelTime ) & & * ( SOMass - xy_SOSeaIceMass(i,j) ) & & + ( pyz_SOSeaIceDiffCoef(i-1,j,k) + pyz_SOSeaIceDiffCoef(i,j,k) ) & & * y_Factor(j) aax_LUMat(l,ii,i+1) = & & - pyz_SOSeaIceDiffCoef(i ,j,k) & & * y_Factor(j) end do do ii = imax-1, imax-1 i = ii aax_LUMat(l,ii,i-1) = & & - pyz_SOSeaIceDiffCoef(i-1,j,k) & & * y_Factor(j) aax_LUMat(l,ii,i ) = & & 1.0_DP / ( 1.0_DP * DelTime ) & & * ( SOMass - xy_SOSeaIceMass(i,j) ) & & + ( pyz_SOSeaIceDiffCoef(i-1,j,k) + pyz_SOSeaIceDiffCoef(i,j,k) ) & & * y_Factor(j) aax_LUMat(l,ii,0 ) = & & - pyz_SOSeaIceDiffCoef(imax-1,j,k) & & * y_Factor(j) end do end do end do ! do k = 1, ksimax do j = 1, jmax l = (k-1)*jmax+j do ii = 0, 0 i = ii if ( SOMass - xy_SOSeaIceMass(i,j) <= 0.0_DP ) then aax_LUMat(l,ii,imax-1) = 0.0_DP aax_LUMat(l,ii,i ) = 1.0_DP aax_LUMat(l,ii,i+1 ) = 0.0_DP end if end do do ii = 0+1, (imax-1)-1 i = ii if ( SOMass - xy_SOSeaIceMass(i,j) <= 0.0_DP ) then aax_LUMat(l,ii,i-1) = 0.0_DP aax_LUMat(l,ii,i ) = 1.0_DP aax_LUMat(l,ii,i+1) = 0.0_DP end if end do do ii = imax-1, imax-1 i = ii if ( SOMass - xy_SOSeaIceMass(i,j) <= 0.0_DP ) then aax_LUMat(l,ii,i-1) = 0.0_DP aax_LUMat(l,ii,i ) = 1.0_DP aax_LUMat(l,ii,0 ) = 0.0_DP end if end do end do end do call ludecomp_prep_simple_many ( jmax*ksimax, imax, aax_LUMat, 1, jmax*ksimax ) do k = 1, ksimax do j = 1, jmax l = (k-1)*jmax+j do ii = 0, imax-1 i = ii if ( SOMass - xy_SOSeaIceMass(i,j) <= 0.0_DP ) then aa_LUVec(l,ii) = xy_SurfTemp(i,j) else aa_LUVec(l,ii) = xy_SurfTemp(i,j) & & * ( SOMass - xy_SOSeaIceMass(i,j) ) & & / ( 1.0_DP * DelTime ) end if end do end do end do call ludecomp_solve_simple_many( jmax*ksimax, imax, aax_LUMat, aa_LUVec, 1, jmax*ksimax ) do k = 1, 1 do j = 1, jmax l = (k-1)*jmax+j do ii = 0, imax-1 i = ii xy_SurfTempA(i,j) = aa_LUVec(l,ii) end do end do end do end subroutine SOSIDiffusionX !------------------------------------------------- subroutine SOSIDynamicsInit( & & ArgFlagSlabOcean & & ) ! flag for use of slab ocean ! ! Initialization of module ! ! MPI ! ! 種別型パラメタ ! Kind type parameter ! use dc_types, only: & & STDOUT, & ! 標準出力の装置番号. Unit number of standard output & STRING ! 文字列. Strings. ! ファイル入出力補助 ! File I/O support ! use dc_iounit, only: FileOpen ! NAMELIST ファイル入力に関するユーティリティ ! Utilities for NAMELIST file input ! use namelist_util, only: namelist_filename, NmlutilMsg use mpi_wrapper , only : myrank, nprocs ! ヒストリデータ出力 ! History data output ! use gtool_historyauto, only: HistoryAutoAddVariable ! 組成に関わる配列の設定 ! Settings of array for atmospheric composition ! use composition, only: & & ncmax ! 成分の数 ! Number of composition ! 座標データ設定 ! Axes data settings ! use axesset, only: & & x_Lon, y_Lat, & & AxNameX, AxNameY, AxNameZ, AxNameT use constants0, only : PI ! 雪と海氷の定数の設定 ! Setting constants of snow and sea ice ! use constants_snowseaice, only: & & ConstantsSnowSeaIceInit ! ! Slab ocean sea ice utility module ! use sosi_utils, only : SOSIUtilsInit use sltt_const , only : SLTTConstInit use sosi_dynamics_extarr, only : SLTTExtArrInit use sltt_const , only : iexmin, iexmax, jexmin, jexmax logical, intent(in ) :: ArgFlagSlabOcean ! ! local variables ! integer:: i ! 東西方向に回る DO ループ用作業変数 ! Work variables for DO loop in zonal direction integer:: j ! 南北方向に回る DO ループ用作業変数 ! Work variables for DO loop in meridional direction integer:: unit_nml ! NAMELIST ファイルオープン用装置番号. ! Unit number for NAMELIST file open integer:: iostat_nml ! NAMELIST 読み込み時の IOSTAT. ! IOSTAT of NAMELIST read ! NAMELIST 変数群 ! NAMELIST group name ! namelist /sosi_dynamics_nml/ & & SOSeaIceDiffCoef ! 実行文 ; Executable statement ! if ( sosi_dynamics_inited ) return FlagSlabOcean = ArgFlagSlabOcean if ( mod( jmax, 2 ) /= 0 ) then stop 'jmax cannot be divided by 2.' end if ! Initialization of modules used in this module ! ! 雪と海氷の定数の設定 ! Setting constants of snow and sea ice ! call ConstantsSnowSeaIceInit ! ! Slab ocean sea ice utility module ! call SOSIUtilsInit call SLTTConstInit ! デフォルト値の設定 ! Default values settings ! SOSeaIceDiffCoef = 0.0e0_DP ! NAMELIST の読み込み ! NAMELIST is input ! if ( trim(namelist_filename) /= '' ) then call FileOpen( unit_nml, & ! (out) & namelist_filename, mode = 'r' ) ! (in) rewind( unit_nml ) read( unit_nml, & ! (in) & nml = sosi_dynamics_nml, & ! (out) & iostat = iostat_nml ) ! (out) close( unit_nml ) call NmlutilMsg( iostat_nml, module_name ) ! (in) if ( iostat_nml == 0 ) write( STDOUT, nml = sosi_dynamics_nml ) end if allocate( y_CosLat(1:jmax) ) y_CosLat = cos( y_Lat ) allocate( x_LonS (0:imax-1) ) allocate( x_SinLonS(0:imax-1) ) allocate( x_CosLonS(0:imax-1) ) allocate( y_latS (1:jmax/2) ) allocate( y_SinLatS(1:jmax/2) ) allocate( y_CosLatS(1:jmax/2) ) do i = 0, imax-1 x_LonS (i) = x_Lon(i) x_SinLonS(i) = sin( x_LonS(i) ) x_CosLonS(i) = cos( x_LonS(i) ) end do do j = 1, jmax/2 y_LatS (j) = y_Lat(j) y_SinLatS(j) = sin( y_LatS(j) ) y_CosLatS(j) = cos( y_LatS(j) ) end do allocate( x_LonN (0:imax-1) ) allocate( x_SinLonN(0:imax-1) ) allocate( x_CosLonN(0:imax-1) ) allocate( y_latN (1:jmax/2) ) allocate( y_SinLatN(1:jmax/2) ) allocate( y_CosLatN(1:jmax/2) ) do i = 0, imax-1 x_LonN (i) = x_Lon(i) x_SinLonN(i) = sin( x_LonN(i) ) x_CosLonN(i) = cos( x_LonN(i) ) end do do j = 1, jmax/2 y_LatN (j) = y_Lat(j+jmax/2) y_SinLatN(j) = sin( y_LatN(j) ) y_CosLatN(j) = cos( y_LatN(j) ) end do allocate( x_ExtLonS( iexmin:iexmax ) ) allocate( x_ExtLonN( iexmin:iexmax ) ) allocate( y_ExtLatS( jexmin:jexmax ) ) allocate( y_ExtLatN( jexmin:jexmax ) ) allocate( y_ExtCosLatS( jexmin:jexmax ) ) allocate( y_ExtCosLatN( jexmin:jexmax ) ) call SLTTExtArrInit( & & x_LonS, y_LatS, x_LonN, y_LatN, & ! (in ) & x_ExtLonS, y_ExtLatS, x_ExtLonN, y_ExtLatN & ! (out) & ) y_ExtCosLatS = cos( y_ExtLatS ) y_ExtCosLatN = cos( y_ExtLatN ) allocate( p_LonS (0-1:imax-1+1) ) allocate( q_LatS (1-1:jmax/2+1) ) allocate( q_CosLatS(1-1:jmax/2+1) ) allocate( p_LonN (0-1:imax-1+1) ) allocate( q_LatN (1-1:jmax/2+1) ) allocate( q_CosLatN(1-1:jmax/2+1) ) do i = 0-1, imax-1+1 p_LonS(i) = ( x_ExtLonS(i) + x_ExtLonS(i+1) ) / 2.0_DP p_LonN(i) = ( x_ExtLonN(i) + x_ExtLonN(i+1) ) / 2.0_DP end do do j = 1-1, jmax/2+1 q_LatS(j) = ( y_ExtLatS(j) + y_ExtLatS(j+1) ) / 2.0_DP end do do j = 1-1, jmax/2+1 q_LatN(j) = ( y_ExtLatN(j) + y_ExtLatN(j+1) ) / 2.0_DP end do if ( myrank == nprocs-1 ) then j = 0 q_LatS(j) = - PI / 2.0_DP j = jmax/2+1 q_LatN(j) = PI / 2.0_DP end if q_CosLatS = cos( q_LatS ) q_CosLatN = cos( q_LatN ) ! ヒストリデータ出力のためのへの変数登録 ! Register of variables for history data output ! !!$ do n = 1, ncmax !!$ call HistoryAutoAddVariable( 'SLD'//trim(a_QMixName(n))//'DtHorMassFix', & !!$ & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), & !!$ & 'tendency of horizontal mass fix of '//trim(a_QMixName(n)), 's-1' ) !!$ call HistoryAutoAddVariable( 'SLD'//trim(a_QMixName(n))//'DtVerMassFix', & !!$ & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), & !!$ & 'tendency of vertical mass fix of '//trim(a_QMixName(n)), 's-1' ) !!$ call HistoryAutoAddVariable( 'SLD'//trim(a_QMixName(n))//'DtTotMassFix', & !!$ & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), & !!$ & 'tendency of mass fix of '//trim(a_QMixName(n)), 's-1' ) !!$ end do ! 印字 ; Print ! call MessageNotify( 'M', module_name, '----- Initialization Messages -----' ) call MessageNotify( 'M', module_name, ' SOSeaIceDiffCoef = %f', d = (/ SOSeaIceDiffCoef /) ) !!$ call MessageNotify( 'M', module_name, ' FlagSLTTArcsineVer = %b', l = (/ FlagSLTTArcsineVer /) ) !!$ call MessageNotify( 'M', module_name, ' SLTTArcsineFactor = %f', d = (/ SLTTArcsineFactor /) ) !!$ call MessageNotify( 'M', module_name, ' SLTTIntHor = %c', c1 = trim( SLTTIntHor ) ) !!$ call MessageNotify( 'M', module_name, ' SLTTIntVer = %c', c1 = trim( SLTTIntVer ) ) call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) ) sosi_dynamics_inited = .true. end subroutine SOSIDynamicsInit !-------------------------------------------------------------------------------------- end module sosi_dynamics