!===================================================================== ! 温度分布計算モジュール ! pending ! - x = 1 and IMAX のエネルギー流入量はカウントしてない ! 2003/12/16 ! 2004/01/01 energy conservation を考慮した考え方に変更 ! 2004/01/15 wgt_dammy(KMAX+1)を付け加える ! 2004/01/19 拡散項の導入 !===================================================================== 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 energy2Temp 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) :: energy_tot(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) :: ratio_ab(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 energy2 = 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 の時のエネルギー(格子調節前) end do ! 消耗している場合 if(mass_Ms(i,j) < 0.0D0)then do k = KMAX, 1, -1 if(mass_org(i, j, k)+mass_fl(i, j, k)+mass_Ms(i,j) > 0.0D0)then ! それぞれの質量比で消耗量を出す ratio_ab(i, j, k) = mass_org(i, j, k) & & /(mass_org(i, j, k)+mass_fl(i, j, k)) mass_org(i, j, k) = mass_org(i,j, k)& & + ratio_ab(i, j, k)*mass_Ms(i, j) energy_org(i, j,k)= & & mass_org(i, j, k)*c_p*Temp1(i, j, k) energy_fl(i, j, k) = & & mass_fl(i, j, k) & & +(1.0D0-ratio_ab(i,j,k))* mass_Ms(i, j)& & /mass_fl(i, j, k) *energy_fl(i, j, k) mass_fl(i, j, k) = mass_fl(i,j, k) & & +(1.0D0-ratio_ab(i,j,k))* mass_Ms(i, j) ! 値をリセット mass_Ms(i, j) = 0.0D0 energy_Ms(i, j) =0.0D0 else ! 2 層以上にわたる場合 mass_Ms(i, j)=mass_org(i, j, k)+mass_fl(i, j, k)+mass_Ms(i,j) mass_org(i,j ,k) = 0.0D0 mass_fl(i, j, k) = 0.0D0 energy_org(i, j, k) = 0.0D0 energy_fl(i, j,k) = 0.0D0 end if end do end if ! 調節前のエネルギー energy1(i, j, 1:KMAX)=energy_org(i, j, 1:KMAX)+energy_fl(i, j, 1:KMAX) energy1(i, j, KMAX+1)=energy_Ms(i, j) energy_tot(i, j) = sum(energy1(i, j, :)) ! if((i == 2).and.(j==1))then ! write(6, *)'energy, before', energy1(i, j, :) ! end if ! 氷床の出来はじめ 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==5).and.(j==3))then ! write(6, *)i, j, 'energy', sum(energy1(i,j,:)), energy1(i, j, :) ! 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, 1, -1 ! 調節後の格子 do kk = kk_str, 1, -1 ! 調節前の格子 if(wgt(k) >= wgt_sum)then ! wgt_dammy(kk)の全部 energy2(i, j, k) = energy2(i, j, k)& & +energy1(i, j, kk) ! 入るなら下の層を足してまたチェック if(kk >= 2)then wgt_sum = wgt_sum + wgt_dammy(kk-1) else !kk=1 を足しても余る場合 if((wgt_sum/wgt(k) > 1.000000001).or.& (wgt_sum/wgt(k) > 0.9999999999))then write(6, *)' error!! mass is loss', wgt_sum stop end if end if else ! wgt_dammy(kk)の一部 energy2(i, j, k) = energy2(i, j, k)& & +(wgt(k)-(wgt_sum-wgt_dammy(kk)) )/wgt_dammy(kk) & & * energy1(i, j, kk) ! 残り energy1(i, j, kk) = & & (wgt_sum-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, :))/energy_tot(i,j) > 1.0000001).or. & &(sum(energy2(i, j, :))/energy_tot(i,j) < 0.9999999))then write(6, *)'energy is not conserved. postion is', i, j write(6, *)'energy before', energy_tot(i, j) write(6, *)'energy after', sum(energy2(i, j, :)) stop end if end do end do end subroutine adjust_energy subroutine energy2Temp(energy, Hsfc2, Temp) implicit none real(8), INTENT(IN) :: energy(1:IMAX, 1:JMAX, 1:KMAX) real(8), INTENT(IN) :: Hsfc2(1:IMAX, 1:JMAX) real(8), INTENT(OUT):: Temp(1:IMAX, 1:JMAX, 1:KMAX) real(8) :: mass(1:IMAX, 1:JMAX, 1:KMAX) integer :: i, j, k do i =1, IMAX do j = 1, JMAX if(Hsfc2(i, j) > 0.0D0)then do k =1, KMAX mass(i, j, k) = Hsfc2(i, j)*wgt(k)*dx*rho Temp(i, j, k) = energy(i, j, k)/mass(i, j,k)/c_p end do else Temp(i, j, :) = Ts(i, j) end if end do end do end subroutine energy2Temp end module T_pack