!= 定常状態における \sigma 座標での TEM 計算 ! != Calculation of TEM variables in the \sigma coordinate system ! under steady-state conditions ! ! Authors:: Shin-ich Takehiro ! Version:: $Id: tem.F90,v 1.0 2026/10/13 00:00:00 takepiro Exp $ ! Tag Name:: $Name: $ ! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved. ! License:: See COPYRIGHT[link:../../../COPYRIGHT] ! module tem ! 格子点設定 ! 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 ! 座標データ設定 ! Axes data settings ! use axesset, only: & & y_Lat, & ! $ \varphi $ [rad.] . ! 緯度. Latitude & z_Sigma, & ! $ \sigma $ レベル (整数). ! Full $ \sigma $ level & r_Sigma, & ! $ \sigma $ レベル (半整数). ! Half $ \sigma $ level & z_DelSigma, & ! $ \Delta \sigma $ (整数). ! $ \Delta \sigma $ (Full) & r_DelSigma ! $ \Delta \sigma $ (半整数). ! $ \Delta \sigma $ (Half) ! 組成に関わる配列の設定 ! Settings of array for atmospheric composition ! use composition, only: ncmax, & & IndexH2OVap ! 種別型パラメタ ! Kind type parameter ! use dc_types, only: DP, & ! 倍精度実数型. Double precision. & STRING ! 文字列. Strings. ! メッセージ出力 ! Message output ! use dc_message, only: MessageNotify ! 宣言文 ; Declaration statements ! implicit none private ! 公開手続き ! Public procedure ! public TEMOutput, TEMOutputInit, TEMAnalysis contains ! !------------------------------------------------------------------ ! subroutine TEMOutput( & & xyz_U, xyz_V, xyr_SigDot, xyz_QVap, xyz_Temp, xy_Ps & ! (in) ) ! 時刻管理 ! Time control ! use timeset, only: TimeN ! ステップ $ t $ の時刻. ! Time of step $ t $. use gtool_historyauto, only: HistoryAutoPut ! ! 入力変数 ! real(DP), intent(in) :: xyz_U (0:imax-1, 1:jmax, 1:kmax) ! $ u $ . 東西風速. ! Eastward wind real(DP), intent(in) :: xyz_V (0:imax-1, 1:jmax, 1:kmax) ! $ v $ . 南北風速. ! Northward wind real(DP), intent(in) :: xyr_SigDot (0:imax-1, 1:jmax, 0:kmax) ! $ \dot{\sigma} $ . ! 鉛直流. Vertical flow real(DP), intent(in) :: xyz_QVap(0:imax-1, 1:jmax, 1:kmax) ! $ q (t-\Delta t) $ . ! 比湿. Specific humidity real(DP), intent(in) :: xyz_Temp (0:imax-1, 1:jmax, 1:kmax) ! $ T $ . 温度. ! Temperature real(DP), intent(in) :: xy_Ps (0:imax-1, 1:jmax) ! $ p_s $ . 地表面気圧. ! Surface pressure ! ! TEM 計算用変数 ! real(DP) :: xyz_EpLat (0:imax-1, 1:jmax, 1:kmax) ! $ Ep_lat $ . EP フラックス緯度成分. ! latitudinal EP flux real(DP) :: xyr_EpSig (0:imax-1, 1:jmax, 0:kmax) ! $ Ep_lat $ . EP フラックス鉛直成分. ! vertical EP flux real(DP) :: xyz_PsVRes (0:imax-1, 1:jmax, 1:kmax) ! $ vres $ . 残差南北風速. ! Residual northward wind real(DP) :: xyr_PsSigDotRes (0:imax-1, 1:jmax, 0:kmax) ! $ \dot\sigma_res $ . 残差鉛直風速. ! Residual vertical wind real(DP) :: xyr_StrRes (0:imax-1, 1:jmax, 0:kmax) ! $ \psi_res $ . 残差質量流線関数. ! Residual mass stream fkuction call TEMAnalysis( & & xyz_U, xyz_V, xyr_SigDot, xyz_QVap, xyz_Temp, xy_Ps, & ! (in) & xyz_EPLat, xyr_EPSig, xyz_PsVRes, xyr_PsSigDotRes, & ! (out) & xyr_StrRes & ! (out) & ) call HistoryAutoPut( TimeN, 'VRes', xyz_PsVRes ) call HistoryAutoPut( TimeN, 'SigDotRes', xyr_PsSigDotRes ) call HistoryAutoPut( TimeN, 'StrRes', xyr_StrRes ) call HistoryAutoPut( TimeN, 'EPLat', xyz_EPLat ) call HistoryAutoPut( TimeN, 'EPSig', xyr_EPSig ) end subroutine TEMOutput ! ! 初期化 ! subroutine TEMOutputInit use gtool_historyauto, only: HistoryAutoAddVariable call HistoryAutoAddVariable( 'VRes' , & ! (in) & (/ 'lon ', 'lat ', 'sig ', 'time' /), & ! (in) & 'residual northward wind', 'm s-1' ) ! (in) call HistoryAutoAddVariable( 'SigDotRes', & ! (in) & (/ 'lon ', 'lat ', 'sigm', 'time' /), & ! (in) & 'residual vertical wind', 'm s-1' ) ! (in) call HistoryAutoAddVariable( 'StrRes' , & ! (in) & (/ 'lon ', 'lat ', 'sigm', 'time' /), & ! (in) & 'residual mass stream function', 'kg s-1' ) ! (in) call HistoryAutoAddVariable( 'EPLat' , & ! (in) & (/ 'lon ', 'lat ', 'sig ', 'time' /), & ! (in) & 'northward EP flux', 'm s-1' ) ! (in) call HistoryAutoAddVariable( 'EPSig', & ! (in) & (/ 'lon ', 'lat ', 'sigm', 'time' /), & ! (in) & 'vertical EP flux', 'm s-1' ) ! (in) end subroutine TEMOutputInit ! !------------------------------------------------------------------ ! subroutine TEMAnalysis( & & xyz_U, xyz_V, xyr_SigDot, xyz_QVap, xyz_Temp, xy_Ps, & ! (in) & xyz_EPLat, xyr_EPSig, xyz_PsVRes, xyr_PsSigDotRes, & ! (out) & xyr_StrRes & ! (out) & ) ! 物理・数学定数 use constants0, only : PI ! $ \pi $ . ! 円周率. Circular constant ! 物理定数設定 ! Physical constants settings ! use constants, only: & & RPlanet,& ! $ a $ [m]. ! 惑星半径. ! Radius of planet & CpDry,& ! $ C_p $ [J kg-1 K-1]. ! 乾燥大気の定圧比熱. ! Specific heat of air at constant pressure & Omega,& ! $ \Omega $ [s-1]. ! 回転角速度. ! Angular velocity & Grav ! $ g $ [m s-2]. ! 重力加速度. ! Gravitational acceleration ! 温度の半整数σレベルの補間, 気圧と高度の算出 ! Interpolate temperature on half sigma level, ! and calculate pressure and height ! use auxiliary, only: AuxVars ! ! 入力用変数 ! real(DP), intent(in) :: xyz_U (0:imax-1, 1:jmax, 1:kmax) ! $ u $ . 東西風速. ! Eastward wind real(DP), intent(in) :: xyz_V (0:imax-1, 1:jmax, 1:kmax) ! $ v $ . 南北風速. ! Northward wind real(DP), intent(in) :: xyr_SigDot (0:imax-1, 1:jmax, 0:kmax) ! $ \dot{\sigma} $ . ! 鉛直流. Vertical flow real(DP), intent(in) :: xyz_QVap(0:imax-1, 1:jmax, 1:kmax) ! $ q (t-\Delta t) $ . ! 比湿. Specific humidity real(DP), intent(in) :: xyz_Temp (0:imax-1, 1:jmax, 1:kmax) ! $ T $ . 温度. ! Temperature real(DP), intent(in) :: xy_Ps (0:imax-1, 1:jmax) ! $ p_s $ . 地表面気圧. ! Surface pressure ! ! TEM 計算用変数 ! real(DP), intent(out) :: xyz_EpLat (0:imax-1, 1:jmax, 1:kmax) ! $ Ep_lat $ . EP フラックス緯度成分. ! latitudinal EP flux real(DP), intent(out) :: xyr_EpSig (0:imax-1, 1:jmax, 0:kmax) ! $ Ep_lat $ . EP フラックス鉛直成分. ! vertical EP flux real(DP), intent(out) :: xyz_PsVRes (0:imax-1, 1:jmax, 1:kmax) ! $ vres $ . 残差南北風速. ! Residual northward wind real(DP), intent(out) :: xyr_PsSigDotRes (0:imax-1, 1:jmax, 0:kmax) ! $ \dot\sigma_res $ . 残差鉛直風速. ! Residual vertical wind real(DP), intent(out) :: xyr_StrRes (0:imax-1, 1:jmax, 0:kmax) ! $ \psi_res $ . 残差質量流線関数. ! Residual mass stream fkuction integer :: j, k real(DP) :: xyz_Exner (0:imax-1, 1:jmax, 1:kmax) ! Exner function real(DP) :: xyr_Exner (0:imax-1, 1:jmax, 0:kmax) ! Exner function real(DP) :: xyr_Temp (0:imax-1, 1:jmax, 0:kmax) ! $ T $ . 温度. ! Temperature real(DP) :: xyz_PotTemp (0:imax-1, 1:jmax, 1:kmax) ! Potential temperature real(DP) :: xyz_PsV (0:imax-1, 1:jmax, 1:kmax) ! $ p_s v $ . 南北風速. ! Northward wind real(DP) :: xyr_PsV (0:imax-1, 1:jmax, 0:kmax) ! $ p_s v $ . 南北風速. ! Northward wind real(DP) :: xyz_PsVMean (0:imax-1, 1:jmax, 1:kmax) ! $ v $ . 帯状平均南北風速. ! Zonal mean northward wind real(DP) :: xyz_PsVDist (0:imax-1, 1:jmax, 1:kmax) ! $ v $ . 擾乱南北風速. ! Disturbance northward wind real(DP) :: xyr_PsVMean (0:imax-1, 1:jmax, 0:kmax) ! $ v $ . 帯状平均南北風速. ! Zonal mean northward wind real(DP) :: xyr_PsVDist (0:imax-1, 1:jmax, 0:kmax) ! $ v $ . 擾乱南北風速. ! Disturbance northward wind real(DP) :: xyr_PsSigDot (0:imax-1, 1:jmax, 0:kmax) ! $ \dot{\sigma} $ . ! 鉛直流. vertical flow real(DP) :: xyr_PsSigDotMean (0:imax-1, 1:jmax, 0:kmax) ! $ \dot{\sigma} $ . ! 帯状平均鉛直流. Zonal mean vertical flow real(DP) :: xyr_PsSigDotDist (0:imax-1, 1:jmax, 0:kmax) ! $ \dot{\sigma} $ . ! 擾乱鉛直流. Zonal mean vertical flow real(DP) :: xyz_Entro (0:imax-1, 1:jmax, 1:kmax) ! Entropy real(DP) :: xyr_Entro (0:imax-1, 1:jmax, 0:kmax) ! Entropy real(DP) :: xyz_EntroMean (0:imax-1, 1:jmax, 1:kmax) ! Zonal Mean Entropy real(DP) :: xyz_EntroDist (0:imax-1, 1:jmax, 1:kmax) ! Disturbance Entropy real(DP) :: xyr_EntroMean (0:imax-1, 1:jmax, 0:kmax) ! Zonal Mean Entropy real(DP) :: xyr_EntroDist (0:imax-1, 1:jmax, 0:kmax) ! Disturbance Entropy real(DP) :: xyz_DEntroMeanDSigma (0:imax-1, 1:jmax, 1:kmax) ! Vertical derivative of Zonal Mean Entropy real(DP) :: xyr_DEntroMeanDSigma (0:imax-1, 1:jmax, 0:kmax) ! Vertical derivative of Zonal Mean Entropy real(DP) :: xyz_EntroVDist (0:imax-1, 1:jmax, 1:kmax) ! $ s'v' $ . 擾乱南北エントロピーフラックス ! Disturbance northern entropy flux real(DP) :: xyr_EntroVDist (0:imax-1, 1:jmax, 0:kmax) ! $ s'v' $ . 擾乱南北エントロピーフラックス ! Disturbance northern entropy flux real(DP) :: xyz_Amp (0:imax-1, 1:jmax, 1:kmax) ! planetary angular momentum real(DP) :: xyz_Am (0:imax-1, 1:jmax, 1:kmax) ! local angular momentum real(DP) :: xyz_AmMean (0:imax-1, 1:jmax, 1:kmax) ! zonal mean local angular momentum real(DP) :: xyz_AmDist (0:imax-1, 1:jmax, 1:kmax) ! disuturbance local angular momentum ! 温度の半整数σレベルの補間, 気圧と高度の算出 ! Interpolate temperature on half sigma level, ! and calculate pressure and height ! call AuxVars( & & xy_Ps, xyz_Temp, xyz_QVap, & ! (in ) & xyr_Temp = xyr_Temp, & ! (out) optional & xyz_Exner = xyz_Exner, & ! (out) optional & xyr_Exner = xyr_Exner & ! (out) optional & ) !---- \overline{V}, \overline{W} の計算 xyz_PsV = spread(xy_Ps,3,kmax)*xyz_V xyz_PsVMean = xya_AvrLon_xya(xyz_PsV) xyz_PsVDist = xyz_PsV - xyz_PsVMean xyr_PsSigDot = spread(xy_Ps,3,kmax+1)*xyr_SigDot xyr_PsSigDotMean = xya_AvrLon_xya(xyr_PsSigDot) xyr_PsSigDotDist = xyr_PsSigDot-xyr_PsSigDotMean !---- \bar{s}_\sigma の計算 xyz_Entro = CpDry * log(xyz_Temp / xyz_Exner) xyz_EntroMean = xya_AvrLon_xya(xyz_Entro) xyz_EntroDist = xyz_Entro-xyz_EntroMean xyr_DEntroMeanDSigma = xyr_DSigma_xyz(xyz_EntroMean) xyz_DentroMeanDSigma = xyz_xyr(xyr_DentroMeanDSigma) !---- \overline{s'V'} の計算 xyz_EntroVDist = xya_AvrLon_xya(xyz_PsVDist*xyz_EntroDist) xyr_EntroVDist = xyr_xyz(xyz_EntroVDist) !---- 残差循環計算 xyz_PsVRes = xyz_PsVMean & - xyz_xyr(xyr_DSigma_xyz(xyz_EntroVDist/xyz_DentroMeanDSigma)) xyr_PsSigDotRes = xyr_PsSigDotMean & + xyr_DivLat_xyr(xyr_EntroVDist/xyr_DentroMeanDSigma) xyr_StrRes(:,:,kmax) = 0.0D0 do j=1,jmax do k=kmax-1,0,-1 xyr_StrRes(:,j,k) & = xyr_StrRes(:,j,k+1) & + 2*PI*xyz_PsVRes(:,j,k+1) * z_DelSigma(k+1) end do xyr_StrRes(:,j,:) & = xyr_StrRes(:,j,:) * Rplanet * cos(y_Lat(j))/Grav end do ! !---- 角運動量計算 ! do j=1,jmax xyz_Amp(:,j,:) = Omega * (Rplanet * cos(y_Lat(j)))**2 xyz_Am(:,j,:) = xyz_U(:,j,:) * Rplanet * cos(y_Lat(j)) end do xyz_AmMean = xya_AvrLon_xya(xyz_Am) xyz_AmDist = xyz_Am - xyz_AmMean !---- EP flux 計算 xyz_EPlat = - xya_AvrLon_xya(xyz_AmDist*xyz_PsVDist) & + xyz_Dsigma_xyr(xyr_xyz(xyz_AmMean)) & * (xyz_EntroVDist/xyz_DentroMeanDSigma) xyr_EPSig = - xya_AvrLon_xya(xyr_xyz(xyz_AmDist)*xyr_PsSigDotDist) & - xyr_GradLat_xyr(xyr_xyz(xyz_Am+xyz_Amp)) & * (xyr_EntroVDist/xyr_DentroMeanDSigma) end subroutine TEMAnalysis ! !====================================================================== ! !---- 経度平均 (経度座標をキープ) ! function xya_AvrLon_xya(xya_Data) result(xya_Avr) use intavr_operate, only : ya_AvrLon_xya real(DP), intent(in) :: xya_Data (:, :, :) real(DP) :: xya_Avr (0:imax-1, 1:jmax, size(xya_Data,3)) xya_Avr = spread(ya_AvrLon_xya(xya_Data),1,imax) end function xya_AvrLon_xya ! !---- 鉛直補間 ! function xyr_xyz(xyz_Data) result(xyr_Data) real(DP), intent(in) :: xyz_Data (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyr_Data (0:imax-1, 1:jmax, 0:kmax) integer :: i,j,k real(DP) :: z1, z2, z3 z1 = r_Sigma(0)-z_Sigma(1) z2 = r_Sigma(0)-z_Sigma(2) z3 = r_Sigma(0)-z_Sigma(3) xyr_Data(:,:,0) & = z2*z3*xyz_Data(:,:,1)/((z2-z1)*(z3-z1)) & + z3*z1*xyz_Data(:,:,2)/((z1-z2)*(z3-z2)) & + z1*z2*xyz_Data(:,:,3)/((z1-z3)*(z2-z3)) do k=1,kmax-1 xyr_Data(:,:,k) & = ( (r_Sigma(k)-z_Sigma(k+1))*xyz_Data(:,:,k) & + (z_Sigma(k)-r_Sigma(k))*xyz_Data(:,:,k+1) ) & /r_DelSigma(k) end do z1 = r_Sigma(kmax) - z_Sigma(kmax) z2 = r_Sigma(kmax) - z_Sigma(kmax-1) z3 = r_Sigma(kmax) - z_Sigma(kmax-2) xyr_Data(:,:,kmax) & = z2*z3*xyz_Data(:,:,kmax) /((z2-z1)*(z3-z1)) & + z3*z1*xyz_Data(:,:,kmax-1)/((z1-z2)*(z3-z2)) & + z1*z2*xyz_Data(:,:,kmax-2)/((z1-z3)*(z2-z3)) end function xyr_xyz ! !---- 鉛直補間 (r->z) ! function xyz_xyr(xyr_Data) result(xyz_Data) real(DP), intent(in) :: xyr_Data (0:imax-1, 1:jmax, 0:kmax) real(DP) :: xyz_Data (0:imax-1, 1:jmax, 1:kmax) integer :: i,j,k do k=1,kmax xyz_Data(:,:,k) & = ( (z_Sigma(k)-r_Sigma(k))*xyr_Data(:,:,k-1) & + (r_Sigma(k-1)-z_Sigma(k))*xyr_Data(:,:,k) ) & /z_DelSigma(k) end do end function xyz_xyr ! !---- 鉛直微分 (r->z) ! function xyz_DSigma_xyr(xyr_Data) result(xyz_Data) real(DP), intent(in) :: xyr_Data (0:imax-1, 1:jmax, 0:kmax) real(DP) :: xyz_Data (0:imax-1, 1:jmax, 1:kmax) integer :: i,j,k do k=1,kmax xyz_Data(:,:,k) & = (xyr_Data(:,:,k-1)-xyr_Data(:,:,k))/z_DelSigma(k) end do end function xyz_DSigma_xyr ! !---- 鉛直微分 (z->r) ! function xyr_DSigma_xyz(xyz_In) result(xyr_Out) real(DP), intent(in) :: xyz_In (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyr_Out (0:imax-1, 1:jmax, 0:kmax) integer :: i,j,k real(DP) :: z1, z2, z3 do k=1,kmax-1 xyr_Out(:,:,k) = (xyz_In(:,:,k)-xyz_In(:,:,k+1))/r_DelSigma(k) end do z1 = r_Sigma(0)-z_Sigma(1) z2 = r_Sigma(0)-z_Sigma(2) z3 = r_Sigma(0)-z_Sigma(3) xyr_Out(:,:,0) & = - (z2+z3) *xyz_In(:,:,1)/((z1-z2)*(z3-z1)) & - (z3+z1) *xyz_In(:,:,2)/((z1-z2)*(z2-z3)) & - (z1+z2) *xyz_In(:,:,3)/((z2-z3)*(z3-z1)) z1 = r_Sigma(kmax) - z_Sigma(kmax) z2 = r_Sigma(kmax) - z_Sigma(kmax-1) z3 = r_Sigma(kmax) - z_Sigma(kmax-2) xyr_Out(:,:,kmax) & = - (z2+z3) *xyz_In(:,:,kmax)/((z1-z2)*(z3-z1)) & - (z3+z1) *xyz_In(:,:,kmax-1)/((z1-z2)*(z2-z3)) & - (z1+z2) *xyz_In(:,:,kmax-2)/((z2-z3)*(z3-z1)) end function xyr_DSigma_xyz ! !---- 緯度微分 ! function xyr_GradLat_xyr(xyr_In) result(xyr_Out) ! 物理定数設定 ! Physical constants settings ! use constants, only: & & RPlanet ! $ a $ [m]. ! 惑星半径. ! Radius of planet #ifdef AXISYMMETRY use wa_zonal_module, only: & & wa_xya, & & xya_GradLat_wa #elif LIB_MPI use ua_mpi_module_base, only: & & wa_xya => ua_pva use ua_mpi_module_deriv, only: & & xya_Gradlat_wa => pva_GradLat_ua #else use wa_module, only: & & wa_xya, & & xya_Gradlat_wa #endif real(DP), intent(in) :: xyr_In (0:imax-1, 1:jmax, 0:kmax) real(DP) :: xyr_Out (0:imax-1, 1:jmax, 0:kmax) xyr_Out = xya_GradLat_wa(wa_xya(xyr_In))/Rplanet end function xyr_GradLat_xyr ! !---- 緯度微分 ! function xyr_DivLat_xyr(xyr_In) result(xyr_Out) ! 物理定数設定 ! Physical constants settings ! use constants, only: & & RPlanet ! $ a $ [m]. ! 惑星半径. ! Radius of planet #ifdef AXISYMMETRY use wa_zonal_module, only: & & xya_wa, & & wa_DivLat_xyaa #elif LIB_MPI use ua_mpi_module_base, only: & & xya_wa => pva_ua use ua_mpi_module_deriv, only: & & wa_Divlat_xya => ua_DivLat_pva #else use wa_module, only: & & xya_wa, & & wa_DivLat_xya #endif real(DP), intent(in) :: xyr_In (0:imax-1, 1:jmax, 0:kmax) real(DP) :: xyr_Out (0:imax-1, 1:jmax, 0:kmax) xyr_Out = xya_wa(wa_DivLat_xya(xyr_In))/Rplanet end function xyr_DivLat_xyr end module tem