!===================================================================== ! 温度分布計算モジュール ! pending ! - x = 1 and IMAX のエネルギー流入量はカウントしてない ! 2003/12/16 ! 2004/01/01 energy conservation を考慮した考え方に変更 !===================================================================== 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 public Ms_minus 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, energy) 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), intent(INOUT) :: energy(1:IMAX, 1:JMAX, 1:KMAX) ! J/m 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_tot(1:IMAX, 1:JMAX) ! 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) :: ratio_Ms(1:IMAX, 1:JMAX) real(8) :: wgt_dammy(1:KMAX) real(8) :: wgt_nonMs real(8), parameter :: non = -9.99E33 integer :: i, j, k integer :: kk_str, kk energy_tot = 0.0D0 ! x 方向境界 energy(1, :, :) = 0.0D0 energy(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-1 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 の時のエネルギー(格子調節前) energy_tot(i, j) = energy_tot(i, j)& & +energy_fl(i, j, k)+energy_org(i, j, k) end do if(energy_Ms(i,j) < 0.0D0)then mass_Ms(i, j) = 0.0D0 energy_Ms(i, j) =0.0D0 end if energy_tot(i, j) = energy_tot(i, j) + energy_Ms(i, j) if(Hsfc2(i, j) > 0.0D0)then wgt_dammy(1:KMAX-1)=& & (mass_fl(i,j,1:KMAX-1)+mass_org(i, j, 1:KMAX-1)) & & /sum(mass(i, j, :)) ! KMAX は上からかん養が加わる wgt_dammy(KMAX)=(mass_fl(i,j,KMAX)+mass_Ms(i, j) & & +mass_org(i, j, KMAX)) & & /sum(mass(i, j, :)) wgt_nonMs = mass_fl(i,j,KMAX)+mass_org(i, j, KMAX) & & /sum(mass(i, j, :)) ratio_Ms(i, j) = mass_Ms(i, j)/sum(mass(i, j, :)) else wgt_dammy(:) = non energy(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) ! write(6, *)i, j, 'energy_tot', energy_tot(i, j) ! end if do k=KMAX, 2, -1 ! k 層に入る質量が少ない場合(もらう) if(wgt_dammy(k) < wgt(k))then ! 下の層からもらう do kk = kk_str, 2, -1 ! kk 層を足しても少ない場合 if(wgt_dammy(k)+wgt_dammy(kk) < wgt(k))then wgt_dammy(k)= wgt_dammy(k)+wgt_dammy(kk) energy(i, j, k) = energy(i, j, k) & & + energy_fl(i, j, kk)+ energy_org(i, j, kk) ! kk 層を足したら多い場合 else wgt_dammy(k+1) = & & wgt_dammy(k+1)+wgt_dammy(k)+wgt_dammy(kk)-wgt(k) energy(i, j, k)= energy(i, j, k) & & + (wgt(k) - wgt_dammy(k)) & & /(wgt_dammy(kk)) & & *(energy_fl(i, j, kk) & & +energy_org(i, j, kk)) ! 残りはそのまま energy(i, j, k-1) = energy(i, j, k-1) & & + (wgt_dammy(k) + wgt_dammy(kk)-wgt(k)) & & /(wgt_dammy(kk)) & & *(energy_fl(i, j, kk) & & +energy_org(i, j, kk)) wgt_dammy(k)= wgt(k) kk_str = kk - 1 exit end if end do else !k 層に入る質量が多い場合(わたす) wgt_dammy(k-1) = wgt_dammy(k-1) + wgt_dammy(k) - wgt(k) ! そのうちかん養で補われる場合 if(wgt_dammy(k) <= ratio_Ms(i, j))then energy(i, j, k)=energy_Ms(i, j)*wgt(k)/(1.0D0 - wgt_nonMs) ! 入らなかったエネルギーは下へ受け渡す energy_fl(i,j,k-1) = energy_fl(i,j,k-1) + energy_fl(i,j,k) energy_org(i,j,k-1)= energy_org(i,j,k-1) + energy_org(i,j,k) ratio_Ms(i, j) = ratio_Ms(i, j) - wgt(k) else ! (かん養)+ 元のエネルギー + 移動によるエネルギー energy(i, j, k) = energy_Ms(i, j)& & * ratio_Ms(i, j)/(mass_Ms(i, j)/sum(mass(i, j, :)))& & + (energy_fl(i, j, k)+energy_org(i, j,k)) & & * (wgt(k)-ratio_Ms(i,j))/(wgt_dammy(k)-ratio_Ms(i,j)) energy_fl(i, j, k-1) = energy_fl(i, j, k-1) & & +energy_fl(i, j, k) & & * (1.0D0- (wgt(k)-ratio_Ms(i,j)) & & /(wgt_dammy(k)-ratio_Ms(i, j))) energy_org(i, j, k-1) = energy_org(i, j, k-1) & & +energy_org(i, j, k) & & * (1.0D0- (wgt(k)-ratio_Ms(i,j)) & & /(wgt_dammy(k)-ratio_Ms(i, j))) ratio_Ms(i, j) = 0.0D0 end if end if wgt_dammy(k) = wgt(k) end do ! k=1 のとき if(ratio_Ms(i, j) >0.0D0)then energy(i, j, 1) = energy_Ms(i, j)*ratio_Ms(i, j)/(1.0D0-wgt_nonMs) & & +energy_fl(i, j, 1)+energy_org(i, j, 1) else energy(i, j, 1)= energy_fl(i, j, 1)+energy_org(i, j, 1) end if ! 格子の切り分けの前後での保存のチェック if(energy_tot(i, j) ==0.0D0)cycle if((sum(energy(i, j, :))/energy_tot(i, j) > 1.0000001).or. & &(sum(energy(i, j, :))/energy_tot(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 subroutine Ms_minus(Hsfc1, Hsfc2, energy, energy2) implicit none real(8), intent(IN) :: Hsfc2(1:IMAX, 1:JMAX) ! m real(8), intent(IN) :: Hsfc1(1:IMAX, 1:JMAX) ! m real(8), intent(INOUT) :: energy(1:IMAX, 1:JMAX, 1:KMAX) ! J/m real(8), intent(INOUT) :: energy2(1:IMAX, 1:JMAX, 1:KMAX) ! J/m real(8) :: energy_tot(1:IMAX, 1:JMAX) ! J/m real(8) :: mass_Ms(1:IMAX, 1:JMAX) ! J/m real(8) :: rest real(8) :: mass(1:IMAX, 1:JMAX, 1:KMAX) real(8) :: mass1(1:IMAX, 1:JMAX, 1:KMAX) integer :: i, j, k, k1, k2 ! 失う質量 mass_Ms = - Ms*dt*dx*rho do i = 1, IMAX do j =1, JMAX mass1(i, j, :) = Hsfc1(i, j) * wgt(:)*dx*rho mass(i, j, :) = Hsfc2(i, j) * wgt(:)*dx*rho ! 総エネルギー energy_tot(i, j) = sum(energy(i, j, :)) if((Ms(i, j) < 0.0D0).and.(Hsfc2(i, j) > 0.0D0))then write(6, *)'Ms < 0', i, j do k = KMAX, 1, -1 ! 昇華が k 層に達する場合 if(mass_Ms(i, j) < mass(i, j, k))then energy(i, j, k) = mass_Ms(i, j)/mass1(i, j, k) & & * energy(i, j, k) energy_tot(i, j) = energy_tot(i, j) & & - energy(i, j, k)*(1.0D0 -mass_Ms(i,j)/mass(i, j, k)) rest = (1.0D0 - mass_Ms(i, j)/mass(i, j, k))& & *Hsfc1(i, j)*rho*dx exit else mass_Ms(i, j) = mass_Ms(i, j) -mass(i, j, k) energy_tot(i, j) = energy_tot(i, j) - energy(i, j, k) end if end do write(6, *)k ! 残りの質量が k2 層目には足りない場合 do k2 = KMAX, 2, -1 ! k2 層は切り直した層 do k1 = k, 2, -1 write(6, *)'rest=', rest, mass(i, j, k2) if( mass(i, j, k2) > rest)then energy2(i, j, k2) = energy2(i, j, k2)& & +energy(i, j, k1) write(6, *)'ene_1', i, j, k2, energy(i, j, k2) ! mass1(i, j, k)=0.0D0 rest = rest + mass1(i, j, k1-1) ! k1 が 2 まで行っても足りない場合 if(k1 == 2)then energy2(i, j, k2) = energy2(i, j, k2) & & + (mass(i, j, k2)+mass1(i, j, k1-1)-rest)& & /mass1(i, j, k1-1)*energy(i, j, k1-1) rest = rest - mass(i, j, k2) ! 残りを分配する energy2(i, j, 1:k2-1)=rest*wgt(1:k2-1)/sum(wgt(1:k2-1)) exit end if ! 残りの質量が k2 層より大きい場合 else ! rest のうち だけ加える energy2(i, j, k2) = energy2(i, j, k2) & & + (mass2(i, j, k2)-rest)/mass1(i,j, k1) * energy(i, j, k1) ! 余りは下の層へ渡す energy2(i, j, k2-1) = energy2(i, j, k2-1) & & + (1.0D0-(mass2(i, j, k2)-rest)/mass1(i, j, k1))*energy(i,j,k1) write(6, *)'ene_2', i, j, k2, energy2(i, j, k2) end if end do end do end if end do end do end subroutine Ms_minus end module T_pack