!--------------------------------------------------------------------- ! Copyright (C) GFD Dennou Club, 2005. All rights reserved. !--------------------------------------------------------------------- ! physics_implicit.f90 ! ! History ! 2005/09/21 Yamada Yukiko create ! module physics_implicit_mod ! != 物理過程 陰解法用モジュール ! !== 概要 ! ! 地表面フラックス, 放射フラックス, 表面温度を陰解法で ! 計算する. ! !== 履歴 (ここに書くので良いのかなあ?) ! (2006-7-30 石渡) subroutine physics_implicit_fluxcorrection 追加 ! implicit none private public :: physics_implicit_init public :: physics_implicit_integrate public :: physics_implicit_fluxcorrection contains subroutine physics_implicit_init( & & xyr_VelLonFlux , & ! (out) 速度経度成分フラックス & xyr_VelLatFlux , & ! (out) 速度緯度成分フラックス & xyr_TempFlux , & ! (out) 温度フラックス & xyr_QvapFlux , & ! (out) 比湿フラックス & xyzo_VelMatrix , & ! (out) 速度陰解行列 & xyzo_TempMatrix , & ! (out) 温度陰解行列 & xyzo_QvapMatrix , & ! (out) 比湿陰解行列 & xyr_Press , & ! (in) 気圧 (半整数) & DelTimePhy , & ! (in) ２Δt & xy_SurfHeatCapacity , & ! (in) 地表熱容量 & xy_SurfCondition ) ! (in) 地表状態 ! use type_mod, only: REKIND, DBKIND, INTKIND, TOKEN, STRING use grid_3d_mod, only: im, jm, km use constants_mod, only: Cp, Grav use dc_trace, only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump implicit none real(DBKIND), intent(out) :: & & xyr_VelLonFlux(im*jm,km+1) , & !速度経度成分フラックス & xyr_VelLatFlux(im*jm,km+1) , & !速度緯度成分フラックス & xyr_TempFlux(im*jm,km+1) , & !温度フラックス & xyr_QvapFlux(im*jm,km+1) , & !比湿フラックス & xyzo_VelMatrix(im*jm,km,-1:1) , & !速度陰解行列 & xyzo_TempMatrix(im*jm,0:km,-1:1) , & !温度陰解行列 & xyzo_QvapMatrix(im*jm,km,-1:1) !比湿陰解行列 real(DBKIND), intent(in) :: & & xyr_Press(im*jm,km+1) , & ! 気圧 (半整数) & DelTimePhy , & ! Δt & xy_SurfHeatCapacity(im*jm) ! 地表熱容量 integer(INTKIND), intent(in) :: & & xy_SurfCondition(im*jm) ! 地表状態 character(STRING), parameter:: subname = "physics_implicit_init" integer(INTKIND) :: ij, k ! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用) continue ! 開始処理 call BeginSub(subname) !---------------------------------------------------------------- ! 陰解行列初期化 !---------------------------------------------------------------- ! ---- 1. 質量, 熱容量の項 ---- xyzo_VelMatrix = 0.0 xyzo_TempMatrix = 0.0 xyzo_QvapMatrix = 0.0 do k = 1, km xyzo_VelMatrix(:,k,0) = ( xyr_Press(:,k) - xyr_Press(:,k+1) ) & & / Grav / DelTimePhy xyzo_TempMatrix(:,k,0) = xyzo_VelMatrix(:,k,0) * Cp xyzo_QvapMatrix(:,k,0) = xyzo_VelMatrix(:,k,0) * Cp end do do ij = 1, im*jm if ( xy_SurfCondition(ij) .GE. 1 ) then xyzo_TempMatrix(ij,0,0) = xy_SurfHeatCapacity(ij) / DelTimePhy else xyzo_TempMatrix(ij,0,0) = 1. end if end do ! ---- 2. フラックスをリセット ---- xyr_VelLonFlux = 0.0 xyr_VelLatFlux = 0.0 xyr_TempFlux = 0.0 xyr_QvapFlux = 0.0 ! 終了処理 call EndSub(subname) end subroutine physics_implicit_init ! (2006-7-28 石渡) ! Verdiff という名前は良くない. T の計算には放射も入っているから. ! subroutine physics_implicit_integrate( & & xyz_DVerdiffVelLonDt , & ! (out) 経度成分 鉛直拡散加速度 & xyz_DVerdiffVelLatDt , & ! (out) 緯度成分 鉛直拡散加速度 & xyz_DVerdiffTempDt , & ! (out) 鉛直拡散加熱率 & xyz_DVerdiffSurfTempDt , & ! (out) 地表面 鉛直拡散加熱率 & xyz_DVerdiffQvapDt , & ! (out) 鉛直拡散加湿率 & xyr_VelLonFlux , & ! (in) 速度経度成分フラックス & xyr_VelLatFlux , & ! (in) 速度緯度成分フラックス & xyr_TempFlux , & ! (in) 温度フラックス & xyr_SurfRadSFlux , & ! (in) 日射フラックス & xyr_SurfRadLFlux , & ! (in) 長波フラックス & xy_GroundTempFlux , & ! (in) 地中熱フラックス & xyr_QvapFlux , & ! (in) 比湿フラックス & xyzo_VelMatrix , & ! (in) 速度陰解行列 & xyzo_TempMatrix , & ! (in) 温度陰解行列 & xyzo_QvapMatrix , & ! (in) 比湿陰解行列 & xy_SurfVelMatrix , & ! (in) 速度陰解行列: 地表 & xyoo_SurfTempMatrix , & ! (in) 温度陰解行列: 地表 & xyoo_SurfQvapMatrix , & ! (in) 比湿陰解行列: 地表 & xyo_SurfRadLMatrix , & ! (in) Ｔ陰解行列：放射 & DelTimePhy , & ! (in) ２Δt & xy_SurfCondition ) ! (in) 地表状態 ! ! 時間変化率の計算 (implicit) use type_mod, only: REKIND, DBKIND, INTKIND, TOKEN, STRING use grid_3d_mod, only: im, jm, km use constants_mod, only: EL, Cp use dc_trace, only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump implicit none real(DBKIND), intent(out) :: & & xyz_DVerdiffVelLonDt(im*jm,km), & ! 経度成分 鉛直拡散加速度 & xyz_DVerdiffVelLatDt(im*jm,km), & ! 緯度成分 鉛直拡散加速度 & xyz_DVerdiffTempDt(im*jm,km), & ! 鉛直拡散加熱率 & xyz_DVerdiffSurfTempDt(im*jm), & ! 地表面 鉛直拡散加熱率 & xyz_DVerdiffQvapDt(im*jm,km) ! 鉛直拡散加湿率 real(DBKIND), intent(in) :: & & xyr_VelLonFlux(im*jm,km+1), & ! 速度経度成分フラックス & xyr_VelLatFlux(im*jm,km+1), & ! 速度緯度成分フラックス & xyr_TempFlux(im*jm,km+1), & ! 温度フラックス & xyr_SurfRadSFlux(im*jm), & ! 日射フラックス & xyr_SurfRadLFlux(im*jm), & ! 長波フラックス & xy_GroundTempFlux(im*jm) , & ! 地中熱フラックス & xyr_QvapFlux(im*jm,km+1) , & ! 比湿フラックス & xyzo_VelMatrix(im*jm,km,-1:1) , & ! 速度陰解行列 & xyzo_TempMatrix(im*jm,0:km,-1:1) , & ! 温度陰解行列 & xyzo_QvapMatrix(im*jm,km,-1:1) , & ! 比湿陰解行列 & xy_SurfVelMatrix(im*jm) , & ! 速度陰解行列: 地表 & xyoo_SurfTempMatrix(im*jm,0:1,-1:1) , & ! 温度陰解行列: 地表 & xyoo_SurfQvapMatrix(im*jm,0:1,-1:1) , & ! 比湿陰解行列: 地表 & xyo_SurfRadLMatrix(im*jm,-1:1) , & ! Ｔ陰解行列：放射 & DelTimePhy ! ２Δt integer(INTKIND) :: xy_SurfCondition(im*jm) ! 地表状態 !----- 作業用内部変数 ----- character(STRING), parameter:: subname = "physics_implicit_integrate" integer(INTKIND) :: ij, k, l ! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用) real(DBKIND) :: & & xyz_DelTempQvap(im*jm,-km:km) , & ! Ｔｑ時間変化 & xyzo_TempQvapLUMatrix(im*jm,-km:km,-1:1) , & ! ＬＵ行列 & xyzo_VelLUMatrix(im*jm,km,-1:1) ! ＬＵ行列 continue ! 開始処理 call BeginSub(subname) !---------------------------------------------------------------- ! 陰解行列計算 !---------------------------------------------------------------- ! ---- 1. 速度 (Vlon, Vlat) の解 ---- xyzo_VelLUMatrix = xyzo_VelMatrix xyzo_VelLUMatrix(:,1,0) = xyzo_VelLUMatrix(:,1,0) + xy_SurfVelMatrix(:) call lu_decomposition_tridiagonal( & & xyzo_VelLUMatrix, im*jm, km ) do k = 1, km xyz_DVerdiffVelLonDt(:,k) = xyr_VelLonFlux(:,k) - xyr_VelLonFlux(:,k+1) xyz_DVerdiffVelLatDt(:,k) = xyr_VelLatFlux(:,k) - xyr_VelLatFlux(:,k+1) end do call lu_solve_tridiagonal( & & xyz_DVerdiffVelLonDt , & & xyzo_VelLUMatrix , & & 1, im*jm, km ) call lu_solve_tridiagonal( & & xyz_DVerdiffVelLatDt , & & xyzo_VelLUMatrix , & & 1, im*jm, km ) ! ---- 2. 温度と比湿の解 ---- do l = -1, 1 do k = 1, km xyzo_TempQvapLUMatrix(:,k,l) = xyzo_TempMatrix(:,k,l) xyzo_TempQvapLUMatrix(:,-k,-l) = xyzo_QvapMatrix(:,k,l) end do xyzo_TempQvapLUMatrix(:,1,l) = xyzo_TempMatrix(:,1,l) & & + xyoo_SurfTempMatrix(:,1,l) xyzo_TempQvapLUMatrix(:,-1,-l) = xyzo_QvapMatrix(:,1,l) & & + xyoo_SurfQvapMatrix(:,1,l) end do xyzo_TempQvapLUMatrix(:,0,0) = xyzo_TempMatrix(:,0,0) & & + xyoo_SurfTempMatrix(:,0,0) + xyoo_SurfQvapMatrix(:,0,0) & & + xyo_SurfRadLMatrix(:,0) xyzo_TempQvapLUMatrix(:,0,1) = & & + xyoo_SurfTempMatrix(:,0,1) + xyo_SurfRadLMatrix(:,1) xyzo_TempQvapLUMatrix(:,0,-1) = xyoo_SurfQvapMatrix(:,0,1) call lu_decomposition_tridiagonal( & & xyzo_TempQvapLUMatrix, im*jm, 2*km+1 ) do k = 1, km xyz_DelTempQvap(:,k) = xyr_TempFlux(:,k) - xyr_TempFlux(:,k+1) xyz_DelTempQvap(:,-k) = xyr_QvapFlux(:,k) - xyr_QvapFlux(:,k+1) end do xyz_DelTempQvap(:,0) = - xyr_SurfRadSFlux - xyr_SurfRadLFlux & & - xyr_TempFlux(:,1) - xyr_QvapFlux(:,1) & & + xy_GroundTempFlux call lu_solve_tridiagonal( & & xyz_DelTempQvap , & & xyzo_TempQvapLUMatrix , & & 1, im*jm, 2*km+1 ) ! ---- 2. 時間変化率 ---- do k = 1, km xyz_DVerdiffVelLonDt(:,k) = xyz_DVerdiffVelLonDt(:,k) / DeltimePhy xyz_DVerdiffVelLatDt(:,k) = xyz_DVerdiffVelLatDt(:,k) / DeltimePhy xyz_DVerdiffTempDt(:,k) = xyz_DelTempQvap(:,k) / DeltimePhy xyz_DVerdiffQvapDt(:,k) = xyz_DelTempQvap(:,-k) / DeltimePhy / EL * Cp end do do ij = 1, im*jm if ( xy_SurfCondition(ij) .GE. 1 ) then xyz_DVerdiffSurfTempDt(ij) = xyz_DelTempQvap(ij,0) / DeltimePhy else xyz_DVerdiffSurfTempDt(ij) = 0. end if end do !---------------------------------------------------------------- ! 終了処理 !---------------------------------------------------------------- call EndSub(subname) end subroutine physics_implicit_integrate subroutine lu_decomposition_tridiagonal( & & jno_LUMatrix, JDimMax, NDimMax ) !==== Dependency use type_mod, only: REKIND, DBKIND, INTKIND, TOKEN, STRING use dc_trace, only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump implicit none !==== Parameter ! integer(INTKIND), intent(in) :: & JDimMax, NDimMax !==== In/Out ! real(DBKIND), intent(inout) :: & jno_LUMatrix(JDimMax, NDimMax, -1:1) ! 入力／ＬＵ行列 !----- 作業用内部変数 ----- character(STRING), parameter:: subname = "lu_decomposition_tridiagonal" integer(INTKIND) :: j, n ! do ループ用作業変数 continue !---------------------------------------------------------------- ! 開始処理 !---------------------------------------------------------------- call BeginSub(subname) !---------------------------------------------------------------- ! 行列のＬＵ分解 [ ３重対角行列 ] !---------------------------------------------------------------- do j = 1, JDimMax jno_LUMatrix(j,1,1) = jno_LUMatrix(j,1,1) / jno_LUMatrix(j,1,0) end do do n = 2, NDimMax-1 do j = 1, JDimMax jno_LUMatrix(j,n,0) = jno_LUMatrix(j,n,0) & & - jno_LUMatrix(j,n,-1) * jno_LUMatrix(j,n-1,1) jno_LUMatrix(j,n,1) = jno_LUMatrix(j,n,1) /jno_LUMatrix(j,n,0) end do end do do j = 1, JDimMax jno_LUMatrix(j,NDimMax,0) = jno_LUMatrix(j,NDimMax,0) & & - jno_LUMatrix(j,NDimMax,-1) * jno_LUMatrix(j,NDimMax-1,1) end do !---------------------------------------------------------------- ! 終了処理 !---------------------------------------------------------------- call EndSub(subname) end subroutine lu_decomposition_tridiagonal subroutine lu_solve_tridiagonal( & & ijn_Vector , & & jno_LUMatrix , & & IDimMax, JDimMax, NDimMax ) !==== Dependency use type_mod, only: REKIND, DBKIND, INTKIND, TOKEN, STRING use dc_trace, only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump implicit none !==== Parameter ! integer(INTKIND), intent(in) :: & IDimMax, JDimMax, NDimMax !==== In/Out ! real(DBKIND), intent(inout) :: & & ijn_Vector(IDimMax, JDimMax, NDimMax) ! 右辺ベクトル／解 !==== Input ! real(DBKIND), intent(in) :: & jno_LUMatrix(JDimMax, NDimMax, -1:1) ! ＬＵ行列 !----- 作業用内部変数 ----- character(STRING), parameter:: subname = "lu_solve_tridiagonal" integer(INTKIND) :: i, j, n ! do ループ用作業変数 continue !---------------------------------------------------------------- ! 開始処理 !---------------------------------------------------------------- call BeginSub(subname) !---------------------------------------------------------------- ! ＬＵ分解による解の計算 [ ３重対角行列 ] !---------------------------------------------------------------- ! ---- 1. 前進代入 ---- do i = 1, IDimMax do j = 1, JDimMax ijn_Vector(i,j,1) = ijn_Vector(i,j,1) / jno_LUMatrix(j,1,0) end do end do do n = 2, NDimMax do i = 1, IDimMax do j = 1, JDimMax ijn_Vector(i,j,n) = ( ijn_Vector(i,j,n) & & - ijn_Vector(i,j,n-1) * jno_LUMatrix(j,n,-1) ) & & / jno_LUMatrix(j,n,0) end do end do end do ! ---- 2. 後退代入 ---- do n = NDimMax-1, 1, -1 do i = 1, IDimMax do j = 1, JDimMax ijn_Vector(i,j,n) = ijn_Vector(i,j,n) & & - ijn_Vector(i,j,n+1) * jno_LUMatrix(j,n,1) end do end do end do !---------------------------------------------------------------- ! 終了処理 !---------------------------------------------------------------- call EndSub(subname) end subroutine lu_solve_tridiagonal subroutine physics_implicit_fluxcorrection ( & & xyr_VelLonFlux, & !(inout) & xyr_VelLatFlux, & !(inout) & xyr_TempFlux, & !(inout) & xyr_QvapFlux, & !(inout) & xyz_DVerdiffVelLonDt, & !(in) & xyz_DVerdiffVelLatDt, & !(in) & xyz_DVerdiffTempDt, & !(in) & xyz_DVerdiffSurfTempDt, & !(in) & xyz_DVerdiffQvapDt, & !(in) & xyzo_VelMatrix, & !(in) & xyzo_TempMatrix, & !(in) & xyzo_QvapMatrix, & !(in) & xy_SurfVelMatrix, & !(in) & xyoo_SurfTempMatrix, & !(in) & xyoo_SurfQvapMatrix, & !(in) & DelTimePhy & !(in) & ) ! ! フラックスの補正 ! ! use type_mod, only: DBKIND, INTKIND, STRING use grid_3d_mod, only: im, jm, km use constants_mod, only: EL, Cp use dc_trace, only: BeginSub, EndSub real(DBKIND), intent(inout) :: xyr_VelLonFlux(im,jm,km+1) ! Ｕのフラックス real(DBKIND), intent(inout) :: xyr_VelLatFlux(im,jm,km+1) ! Ｖのフラックス real(DBKIND), intent(inout) :: xyr_TempFlux(im,jm,km+1) ! Ｔのフラックス real(DBKIND), intent(inout) :: xyr_QvapFlux(im,jm,km+1) ! ｑのフラックス real(DBKIND), intent(in) :: xyz_DVerdiffVelLonDt(im,jm,km) ! 東西運動量変化項ＵＡ real(DBKIND), intent(in) :: xyz_DVerdiffVelLatDt(im,jm,km) ! 南北運動量変化項ＶＡ real(DBKIND), intent(in) :: xyz_DVerdiffTempDt(im,jm,km) ! 温度時間変化項Ｈ real(DBKIND), intent(in) :: xyz_DVerdiffSurfTempDt(im,jm) ! 地表温度変化率 real(DBKIND), intent(in) :: xyz_DVerdiffQvapDt(im,jm,km) ! 比湿時間変化項Ｒ real(DBKIND), intent(in) :: xyzo_VelMatrix(im,jm,km,-1:1) ! ｕ陰解行列 real(DBKIND), intent(in) :: xyzo_TempMatrix(im,jm,0:km,-1:1) ! Ｔ陰解行列 real(DBKIND), intent(in) :: xyzo_QvapMatrix(im,jm,km,-1:1) ! ｑ陰解行列 real(DBKIND), intent(in) :: xy_SurfVelMatrix(im,jm) ! ｕ陰解行列：地表 real(DBKIND), intent(in) :: xyoo_SurfTempMatrix(im,jm,0:1,-1:1) ! Ｔ陰解行列：地表 real(DBKIND), intent(in) :: xyoo_SurfQvapMatrix(im,jm,0:1,-1:1) ! ｑ陰解行列：地表 real(DBKIND), intent(in) :: DelTimePhy ! 時間刻みΔt real(DBKIND) :: ELF INTEGER(INTKIND) :: k character(STRING), parameter:: subname = "physics_implicit_fluxcorrection" ! 開始処理 call BeginSub(subname) ELF = EL/Cp DO k = 2, km xyr_VelLonFlux(:,:,k) = xyr_VelLonFlux(:,:,K) & & - ( xyzo_VelMatrix(:,:,k,-1)* xyz_DVerdiffVelLonDt(:,:,k-1) & & - xyzo_VelMatrix(:,:,k-1,1)* xyz_DVerdiffVelLonDt(:,:,k) & & ) * DelTimePhy xyr_VelLatFlux(:,:,k) = xyr_VelLatFlux(:,:,k) & & - ( xyzo_VelMatrix(:,:,k,-1) * xyz_DVerdiffVelLatDt(:,:,k-1) & & - xyzo_VelMatrix(:,:,k-1,1) * xyz_DVerdiffVelLatDt(:,:,k) & & ) * DelTimePhy xyr_TempFlux(:,:,k) = xyr_TempFlux(:,:,k) & & - ( xyzo_TempMatrix(:,:,k,-1) * xyz_DVerdiffTempDt(:,:,k-1) & & - xyzo_TempMatrix(:,:,k-1,1) * xyz_DVerdiffTempDt(:,:,k) & & ) * DelTimePhy xyr_QvapFlux(:,:,k) = xyr_QvapFlux(:,:,k) & & - ( xyzo_QvapMatrix(:,:,k,-1) * xyz_DVerdiffQvapDt(:,:,k-1) & & - xyzo_QvapMatrix(:,:,k-1,1) * xyz_DVerdiffQvapDt(:,:,k) & & ) * DelTimePhy * ELF end do xyr_VelLonFlux(:,:,1) = xyr_VelLonFlux(:,:,1) & & - xy_SurfVelMatrix(:,:) * xyz_DVerdiffVelLonDt(:,:,1) * DelTimePhy xyr_VelLatFlux(:,:,1) = xyr_VelLatFlux(:,:,1) & & - xy_SurfVelMatrix(:,:) * xyz_DVerdiffVelLatDt(:,:,1) * DelTimePhy xyr_TempFlux(:,:,1) = xyr_TempFlux(:,:,1) & & - ( xyoo_SurfTempMatrix(:,:,1,-1) * xyz_DVerdiffSurfTempDt(:,:) & & + xyoo_SurfTempMatrix(:,:,1,0) * xyz_DVerdiffTempDt(:,:,1) ) & & * DelTimePhy xyr_QvapFlux(:,:,1) = xyr_QvapFlux(:,:,1) & & - ( xyoo_SurfQvapMatrix(:,:,1,-1) * xyz_DVerdiffSurfTempDt(:,:) & & + xyoo_SurfQvapMatrix(:,:,1,0) * xyz_DVerdiffQvapDt(:,:,1) & & * ELF ) * DelTimePhy ! 終了処理 call EndSub(subname) end subroutine physics_implicit_fluxcorrection end module physics_implicit_mod