!===================================================================== ! 温度分布計算モジュール ! pending ! - x = 1 and IMAX のエネルギー流入量はカウントしてない ! 2003/12/16 ! 2004/01/01 energy conservation を考慮した考え方に変更 ! 2004/01/15 wgt_dammy(KMAX+1)を付け加える !===================================================================== module T_pack use prm, only : IMAX, JMAX, KMAX,dt, dx use prmphys, only: k_ice, gamma, rho, c_p use com_var, only: zk, wgt, zk_half use Tbound, only: Ts, Ms implicit none private public get_T public mk_energy public mk_temp1 public adjust_energy contains !===================================================================== ! 各格子のエネルギーの式 !===================================================================== subroutine mk_energy(Temp1, energy) implicit none real(8), intent(IN) :: Temp1(1:IMAX, 1:JMAX, 1:KMAX) real(8), intent(OUT) :: energy(1:IMAX, 1:JMAX, 1:KMAX) energy = Temp1 * rho * c_p end subroutine mk_energy !===================================================================== subroutine get_T(Hsfc1,Temp1, Temp2) implicit none real(8), INTENT(IN) :: Hsfc1(1:IMAX, 1:JMAX) real(8), INTENT(IN) :: Temp1(1:IMAX, 1:JMAX, 1:KMAX) real(8), INTENT(OUT) :: Temp2(1:IMAX, 1:JMAX, 1:KMAX) real(8) :: kappa real(8) :: Tbnd(1:IMAX, 1:JMAX) ! 最下層の下の温度 real(8) :: WORK_aft(1:IMAX, 1:JMAX, 1:KMAX)! 係数 k+1 real(8) :: WORK_cnt(1:IMAX, 1:JMAX, 1:KMAX)! 係数 k real(8) :: WORK_bef(1:IMAX, 1:JMAX, 1:KMAX)! 係数 k-1 integer :: i, j, k kappa = k_ice/(rho*c_p) ! 上部境界値は表面温度 Temp2(:, :, KMAX) = Ts ! 下部境界値は地殻熱流量熱フラックス Tbnd(:, :) = Temp1(:, :, 1)+wgt(1)/k_ice * gamma do j = 1, JMAX do i = 1, IMAX ! 氷床がない場合は氷床温度 if(Hsfc1(i, j) == 0.0)then Temp2(i, j, :) = Ts(i, j) cycle end if ! レベル 1 〜 KMAX-1 までの係数 do k = 1, KMAX-1 WORK_aft(i, j, k) = kappa *dt & & /Hsfc1(i, j)**2/wgt(k)/(zk(k+1) - zk(k)) end do WORK_cnt(i, j, 1) = 1 - kappa *dt & & /Hsfc1(i, j)**2/wgt(1) & & *(1/(zk(2) - zk(1)) + 1/wgt(1)) do k = 2, KMAX-1 WORK_cnt(i, j, k) = 1 - kappa *dt & & /Hsfc1(i, j)**2/wgt(k) & & *(1/(zk(k+1) - zk(k)) + 1/(zk(k) - zk(k-1))) end do WORK_bef(i, j, 1) = kappa *dt & & /Hsfc1(i, j)**2/wgt(1)/wgt(1) do k=2, KMAX WORK_bef(i, j, k) = kappa *dt & & /Hsfc1(i, j)**2/wgt(k)/(zk(k) - zk(k-1)) end do ! 拡散項の計算 Temp2(i, j, 1) = Temp1(i, j, 1) *WORK_cnt(i, j, 1) & & + Tbnd(i, j) *WORK_bef(i, j, 1) & & + Temp1(i, j, 2)*WORK_aft(i, j, 1) do k=2, KMAX-1 Temp2(i, j, k) = Temp1(i, j, k)*WORK_cnt(i, j, k) & & + Temp1(i, j, k-1)*WORK_bef(i, j, k) & & + Temp1(i, j, k+1)*WORK_aft(i, j, k) end do end do end do end subroutine get_T ! subroutine mk_temp1(Hsfc1, Hsfc2, temp1) implicit none real(8), intent(IN) :: Hsfc1(1:IMAX, 1:JMAX) real(8), intent(IN) :: Hsfc2(1:IMAX, 1:JMAX) real(8), intent(INOUT) :: temp1(1:IMAX, 1:JMAX, 1:KMAX) integer :: i, j do i = 1, IMAX do j =1, JMAX if((Hsfc1(i, j) == 0.0D0 ).and.(Hsfc2(i, j) > 0.0D0))then temp1(i, j, :) = Ts(i, j) end if end do end do end subroutine mk_temp1 subroutine adjust_energy(Temp1, Uflgrd, Hsfc1, Hsfc2, energy2) implicit none real(8), intent(IN) :: Uflgrd( 0:IMAX, 1:JMAX, 1:KMAX) ! m^2/s real(8), intent(IN) :: Temp1(1:IMAX, 1:JMAX, 1:KMAX) ! K real(8), intent(IN) :: Hsfc1(1:IMAX, 1:JMAX) ! m real(8), intent(IN) :: Hsfc2(1:IMAX, 1:JMAX) ! m real(8) :: energy1(1:IMAX, 1:JMAX, 1:KMAX+1)! 調節前のエネルギー real(8),intent(OUT) :: energy2(1:IMAX, 1:JMAX, 1:KMAX) ! 調節後 real(8) :: energy_fl(1:IMAX, 1:JMAX, 1:KMAX) ! J/m real(8) :: energy_org(1:IMAX, 1:JMAX, 1:KMAX) ! J/m real(8) :: energy_Ms(1:IMAX, 1:JMAX) ! J/m real(8) :: Uflgrdl_pls(0:IMAX, 1:JMAX, 1:KMAX) !m^2/s real(8) :: Uflgrdl_mns(0:IMAX, 1:JMAX, 1:KMAX) !m^2/s real(8) :: Uflgrdr_pls(0:IMAX, 1:JMAX, 1:KMAX) !m^2/s real(8) :: Uflgrdr_mns(0:IMAX, 1:JMAX, 1:KMAX) !m^2/s real(8) :: mass_fl(1:IMAX, 1:JMAX, 1:KMAX) real(8) :: mass_org(1:IMAX, 1:JMAX, 1:KMAX) real(8) :: mass(1:IMAX, 1:JMAX, 1:KMAX) real(8) :: mass_Ms(1:IMAX, 1:JMAX) real(8) :: wgt_dammy(1:KMAX+1) real(8) :: wgt_sum real(8), parameter :: non = -9.99E33 integer :: i, j, k integer :: kk_str, kk ! x 方向境界 energy1(1, :, :) = 0.0D0 energy1(IMAX, :, :) = 0.0D0 !かん養による質量変化 mass_Ms = Ms*dt*dx*rho energy_Ms = Ms*dt*dx*rho*c_p*Ts do j =1, JMAX do i = 2, IMAX-1 kk_str = KMAX do k = 1, KMAX ! 質量フラックスの符号を分ける ! 格子 i の左側 if(Uflgrd(i-1, j, k) >= 0.0D0)then ! 格子 iへ流入 Uflgrdl_pls(i, j, k) = Uflgrd(i-1, j, k) !* Temp1(i-1, j, k)*c_p Uflgrdl_mns(i, j, k) = 0.0D0 else ! 流出 Uflgrdl_pls(i, j, k) = 0.0D0 Uflgrdl_mns(i, j, k) = Uflgrd(i-1, j, k) !* Temp1(i, j, k)*c_p end if ! 格子 i の右側 if(Uflgrd(i, j, k) >= 0.0D0)then ! 流出 Uflgrdr_pls(i, j, k) = 0.0D0 Uflgrdr_mns(i, j, k) = - Uflgrd(i, j, k) ! * Temp1(i, j, k)*c_p else ! 流入 Uflgrdr_pls(i, j, k) = - Uflgrd(i, j, k) ! * Temp1(i+1, j, k)c_p Uflgrdr_mns(i, j, k) = 0.0D0 end if energy_fl(i, j, k) = & & + (Uflgrdl_pls(i, j, k) * Temp1(i-1, j, k) & & + (Uflgrdl_mns(i, j, k)+Uflgrdr_mns(i, j, k))*Temp1(i,j,k)& & + Uflgrdr_pls(i, j, k) * Temp1(i+1, j, k)) & & * dt * rho * c_p ! 流動による質量変化 mass_fl(i, j, k) = & & + (Uflgrdl_pls(i, j, k)+Uflgrdl_mns(i, j, k) & & + Uflgrdr_pls(i, j, k)+Uflgrdr_mns(i, j, k)) & & * dt * rho ! 時刻 t の時の質量 mass_org(i, j, k) = Hsfc1(i, j)*wgt(k)*dx*rho ! 時刻 t の時のエネルギー energy_org(i, j, k)= mass_org(i, j, k)*c_p*Temp1(i, j, k) ! 時刻 t+1の時の質量 mass(i, j, k) = Hsfc2(i, j)*wgt(k)*dx*rho ! 時刻 t+1 の時のエネルギー(格子調節前) end do if(energy_Ms(i,j) < 0.0D0)then mass_Ms(i, j) = 0.0D0 energy_Ms(i, j) =0.0D0 end if ! 調節前のエネルギー energy1(i, j, 1:KMAX)=energy_org(i, j, 1:KAMX)+energy_fl(i, j, 1:KMAX) energy1(i, j, KMAX+1)=energy_Ms(i, j) if(Hsfc2(i, j) > 0.0D0)then wgt_dammy(1:KMAX)=& & (mass_fl(i,j,1:KMAX)+mass_org(i, j, 1:KMAX)) & & /sum(mass(i, j, :)) ! KMAX+1 層目は上からかん養が加わる wgt_dammy(KMAX+1) = mass_Ms(i, j)/sum(mass(i, j, :)) else wgt_dammy(:) = non energy1(i, j, :) = 0.0D0 cycle end if ! チェック ! if((i==2).and.(j==2))then ! write(6, *)i, j, 'wgt_dammy', sum(wgt_dammy), wgt_dammy ! write(6, *)i, j, 'Hsfc1', Hsfc1(i, j) ! write(6, *)i, j, 'Hsfc2', Hsfc2(i, j) ! end if wgt_sum = wgt_dammy(kk_str) do k=KMAX, 2, -1 ! 調節後の格子 do kk = kk_str, 2, -1 ! 調節前の格子 if(wgt(k) > wgt_sum)then ! wgt_dammy(kk)の全部 energy2(i, j, k) = energy2(i, j, k)& & +energy1(i, j, kk) ! 入るなら下の層を足してまたチェック wgt_sum = wgt_sum + wgt_dammy(kk-1) if(kk==2)then if(wgt(k) > wgt_sum)then energy2(i, j, k) = energy2(i, j, k)+energy1(i, j, kk-1) else energy2(i, j, k) = energy2(i, j, k)& & + (wgt_sum - wgt_dammy(kk-1)-wgt(k))& & /wgt_dammy(kk-1)* energy1(i, j, kk-1) ! 最後の余りは決めた割合に分割 energy2(i, j, 1:k-1) = & & (1.0D0 -(wgt_sum - wgt_dammy(kk-1)-wgt(k))& & /wgt_dammy(kk-1))*energy1(i, j, kk-1) end if else ! wgt_dammy(kk)の一部 energy2(i, j, k) = energy2(i, j, k)& & +(wgt_sum-wgt_dammy(kk)-wgt(k))/wgt_dammy(kk) & & * energy1(i, j, kk) ! 残り energy1(i, j, kk) = & & (1.0D0-(wgt_sum-wgt_dammy(kk)-wgt(k))/wgt_dammy(kk)) & & * energy1(i, j, kk) wgt_sum = wgt_sum-wgt(k) ! 足しすぎた余り kk_str = kk exit end if end do end do ! 格子の切り分けの前後での保存のチェック if(sum(energy1(i, j, :)) ==0.0D0)cycle if((sum(energy2(i, j, :))/sum(energy1(i, j, :)) > 1.0000001).or. & &(sum(energy2(i, j, :))/sum(energy1(i, j, :)) < 0.9999999))then write(6, *)'energy is not conserved. postion is', i, j write(6, *)'energy total', energy_tot(i, j) write(6, *)'energy sum', sum(energy(i, j, :)) stop end if end do end do end subroutine adjust_energy end module T_pack