!===================================================================== ! 温度分布計算モジュール ! pending ! - x = 1 and IMAX のエネルギー流入量はカウントしてない ! 2003/12/16 ! 2004/01/01 energy conservation を考慮した考え方に変更 ! 2004/01/13 Ms < 0 となる場合が考慮されていないことがわかるので考える !===================================================================== 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 add_energy_Ms public add_energy_F public adjust_energy_Ms_ver3 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 add_energy_Ms(energy) implicit none real(8), intent(INOUT) :: energy(1:IMAX, 1:JMAX, 1:KMAX) integer :: i, j do j =1, JMAX do i = 1, IMAX energy(i, j, KMAX) = energy(i, j, KMAX) & & + Ts(i, j) * rho *c_p * Ms(i, j) * dx ! write(6, *) i, j, 'energy(KMAX)=', energy(i, j, KMAX) end do end do end subroutine add_energy_Ms 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 add_energy_F(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_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(1:KMAX) real(8) :: wgt_dammy(1:KMAX) real(8), parameter :: non = -9.99E33 integer :: i, j, k integer :: kk_str, kk ! col_mass = Hsfc2 * dx * rho ! 初期化 energy = 0.0D0 do j =1, JMAX do i = 2, IMAX-1 kk_str = 2 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 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 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 energy_tot(i, j) = energy_tot(i, j) + energy_fl(i, j, k) & & + Hsfc1(i, j)*wgt(k)*rho*dx*c_p*Temp1(i,j,k) end do if(Hsfc2(i, j) > 0.0D0)then wgt_dammy(:) = (mass_fl(i, j, :) + Hsfc1(i, j)*wgt(:)*dx*rho)& & /(Hsfc2(i, j)* rho*dx) ratio(:) = mass_fl(i, j, :)/(Hsfc2(i, j)*rho*dx) else wgt_dammy(:) = non end if if((i==2).and.(j==2))then write(6, *)i, j, 'Uflgrdl_pls', Uflgrdl_pls(i, j, :) write(6, *)i, j, 'Uflgrdr_pls', Uflgrdr_pls(i, j, :) write(6, *)i, j, 'Uflgrdl_mns', Uflgrdl_mns(i, j, :) write(6, *)i, j, 'Uflgrdr_mns', Uflgrdr_mns(i, j, :) write(6, *)i, j, 'energy_tot', energy_tot(i, j) write(6, *)i, j, 'energy_fl', sum(energy_fl(i, j, :)) write(6, *)i, j, 'mass', (Hsfc1(i, j)*rho*dx+sum(mass_fl(i, j, :)))/rho/dx 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, 'ratio', sum(ratio), ratio end if do k=1,KMAX-1 if(wgt_dammy(k) == non)cycle ! k 層に入る質量が少ない場合(もらう) if(wgt_dammy(k) < wgt(k))then ! 上の層からもらう do kk = kk_str, KMAX ! 本当は KMAX -1 までにしないとダメ ! kk 層を足しても少ない場合 if(wgt_dammy(k)+wgt_dammy(kk) < wgt(k))then wgt_dammy(k)= wgt_dammy(k)+wgt_dammy(kk) ! wgt_dammy(kk)のうちわけは ratio : 1-ratio ! ただし mass_fl(raito) < 0 のときはすべて元の氷から if(mass_fl(i, j, k)<0.0D0)then energy(i, j, k) = energy(i, j, k) & & +wgt_dammy(kk)*Hsfc1(i,j)*dx*rho & & *c_p*Temp1(i,j,kk) if((i==2).and.(j==2)) & & write(6, *)i, j, k, 'ene_1', energy(i, j, k) else energy(i, j, k) = energy(i, j, k) & & +wgt_dammy(kk)*(ratio(k)*energy_fl(i, j, k)& & +(1-ratio(k))*Hsfc1(i,j)*dx*rho & & *c_p*Temp1(i,j,kk)) if((i==2).and.(j==2)) & & write(6, *)i, j, k, 'ene_2',energy(i, j, k) end if ! kk 層を足したら多い場合 else wgt_dammy(k+1) = & & wgt_dammy(k+1)+wgt_dammy(k)+wgt_dammy(kk)-wgt(k) ! 元の氷のみ if(mass_fl(i, j, k) < 0.0D0)then energy(i, j, k)= energy(i, j, k) & & + (wgt(k) - wgt_dammy(kk)) & & *Hsfc2(i, j)*dx*rho *c_p*Temp1(i, j, kk) energy(i, j, k+1) = energy(i, j, k+1) & & + (wgt_dammy(k) + wgt_dammy(kk)-wgt(k)) & & *Hsfc1(i, j)*dx*rho *c_p*Temp1(i, j, kk) if((i==2).and.(j==2))& & write(6, *)i, j, k, 'ene_3',energy(i, j, k) ! 元の氷と移動した氷両方 else energy(i, j, k)= energy(i, j, k) & & + (wgt(k) - wgt_dammy(k)) & & /(wgt_dammy(kk)) & & *(energy_fl(i, j, kk) & & +Hsfc1(i, j)*wgt(kk)*dx*rho *c_p*Temp1(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) & & +Hsfc1(i, j)*wgt(kk)*dx*rho *c_p*Temp1(i, j, kk)) ! if((i==2).and.(j==2))& ! & write(6, *)i, j, k, 'ene_4',energy(i, j, k) end if wgt_dammy(k)= wgt(k) kk_str = kk + 1 exit end if end do else !k 層に入る質量が多い場合(わたす) if((i==2).and.(j==2))then write(6, *)"amari(k)", k, (wgt_dammy(k) -wgt(k))*Temp1(i, j, k)*c_p*dx*rho*Hsfc2(i, j) end if wgt_dammy(k+1) = wgt_dammy(k+1) + wgt_dammy(k) - wgt(k) ! 元の氷のみ if(mass_fl(i, j, k)< 0.0D0)then if((i==2).and.(j==2))then write(6, *)"before(k)", k+1, energy(i, j, k+1) end if energy(i, j, k+1) = energy(i, j, k+1) & & + (wgt_dammy(k) - wgt(k))*Hsfc2(i, j)*dx*rho & & * c_p*Temp1(i, j, k) ! write(6, *)'mass_fl <0', energy(i, j, k+1) energy(i, j, k) = wgt(k)*Hsfc2(i, j)*dx*rho*c_p*Temp1(i, j, k) if((i==2).and.(j==2)) & & write(6, *)i, j, k, 'ene_5',energy(i, j, k) if((i==2).and.(j==2))then write(6, *)"after(k)", k+1, energy(i, j, k+1) end if ! 元の氷と移動した氷両方 else energy(i, j,k+1) = energy(i, j, k+1) & & + (1.0D0-wgt(k)/wgt_dammy(k)) *(energy_fl(i, j, k) & & + Hsfc1(i, j)*wgt(k)*dx*rho*c_p*Temp1(i, j, k)) energy(i,j ,k) = energy(i, j, k)& & +wgt(k)/wgt_dammy(k)& & * (energy_fl(i, j, k)+& & Hsfc1(i, j)*wgt(k)*dx*rho*c_p*Temp1(i, j, k)) if((i==2).and.(j==2))& & write(6, *)i, j, k, 'ene_6',energy(i, j, k) end if wgt_dammy(k) = wgt(k) end if end do ! KMAX のみ別で計算 if(mass_fl(i, j, KMAX)<0.0D0)then if((i==2).and.(j==2))& write(6, *) 'ene_bef', energy(i, j, KMAX) !(energy(i, j, KMAX)/Temp1(i, j, KMAX)/c_p+mass_fl(i, j, KMAX)+wgt(KMAX)*Hsfc1(i, j)*dx*rho)/(Hsfc2(i, j)*wgt(KMAX)*dx*rho) energy(i, j, KMAX) = energy(i ,j, KMAX) & + (wgt(KMAX)*Hsfc1(i, j)*dx*rho+mass_fl(i, j, KMAX))& & * c_p*Temp1(i, j, KMAX) if((i==2).and.(j==2))& & write(6, *)i, j, k, 'ene_7',energy(i, j, KMAX), (mass_fl(i, j, KMAX)+wgt(KMAX)*Hsfc1(i, j)*dx*rho)/(Hsfc2(i, j)*wgt(KMAX)*dx*rho) else energy(i, j, KMAX) = energy(i, j, KMAX) + & & wgt(KMAX)*Hsfc1(i, j)*dx*rho & & * c_p*Temp1(i, j, KMAX) + energy_fl(i, j, KMAX) end if if((i==2).and.(j==2))then write(6, *)i, j, 'after energy', sum(energy(i, j,:)), energy(i, j, :) energy_tot(i, j) = 0.0D0 do k = 1, KMAX energy_tot(i, j) = Hsfc2(i, j)*dx*wgt(k)*rho*c_p*temp1(i, j, k) write(6, *)i, j, k, 'after energy2', Hsfc2(i, j)*dx*wgt(k)*rho*c_p*temp1(i, j, k) end do write(6, *)energy_tot(i, j) write(6, *)i, j, 'after wgt_dammy', sum(wgt_dammy), wgt_dammy end if end do end do end subroutine add_energy_F 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 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) if((i==2).and.(j==2))& & write(6, *)i, j, k, 'ene_5',energy(i, j, 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 if((i==2).and.(j==2))& & write(6, *)i, j, k, 'ene_6',energy(i, j, k), Hsfc2(i, j)*wgt(k)*rho*dx*c_p*Temp1(i, j, k), energy_fl(i, j, k-1), energy_org(i, j, k-1) 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((i==2).and.(j==2))then write(6, *)i, j, 'after energy', sum(energy(i, j,:)), energy(i, j, :) write(6, *)i, j, 'after wgt_dammy', sum(wgt_dammy), wgt_dammy end if end do end do end subroutine adjust_energy subroutine adjust_energy_Ms_ver3(Temp1, Hsfc1, Hsfc2, energy2) implicit none real(8), intent(IN) :: Temp1(1:IMAX, 1:JMAX, 1:KMAX) real(8) :: Temp(1:IMAX, 1:JMAX, 1:KMAX+1) real(8), intent(IN) :: Hsfc1(1:IMAX, 1:JMAX) real(8), intent(IN) :: Hsfc2(1:IMAX, 1:JMAX) real(8), intent(INOUT) :: energy2(1:IMAX, 1:JMAX, 1:KMAX) real(8) :: energy(1:IMAX, 1:JMAX, 1:KMAX) real(8) :: ttene(1:IMAX, 1:JMAX) real(8) :: comp_Ms(1:IMAX, 1:JMAX) ! かん養による氷の質量 real(8) :: comp_ice(1:IMAX, 1:JMAX, 1:KMAX) ! 元の氷の質量 real(8) :: rest real(8) :: col_mass(1:IMAX, 1:JMAX) ! カラム中の質量 real(8) :: mass_Ms(1:IMAX, 1:JMAX) ! かん養の質量 ! real(8) :: Mass(1:IMAX, 1:JMAX) ! かん養の質量 integer :: i, j, k , kk integer :: k_str, kk_str ! カラム中の質量 ! col_mass(:, :) = Hsfc1(:, :) * dx + Ms(:, :)*dt*dx ! col_mass(:, :) = Hsfc2(:, :) * dx ! かん養の質量 col_mass = 0.0D0 energy2 = 0.0D0 Temp(:, :, 1:KMAX) = Temp1(:, :, :) Temp(:, :, KMAX+1) = Ts(:, :) ttene(:, :) = ttene(:, :) + Ms(:, :) * dt * rho * dx * c_p * Ts(:, :) do j = 1,JMAX do i = 1, IMAX ! 初期設定 k_str = KMAX ! 時刻 t+1 の時の層番号 kk_str = KMAX ! 時刻 t の時の層番号 kk= kk_str ! かん養量設定 mass_Ms = Ms * dt * dx * rho comp_ice = 0.0D0 ! 時刻 t で氷がない場合 if((Hsfc1(i, j) == 0.0D0).and.(Hsfc2(i, j) > 0.0D0))then write(6, *) 'growing ice in position', i, j comp_ice(i, j, :) = Hsfc2(i, j)*wgt(:) *dx*rho comp_Ms(i, j) = 0.0D0 energy2(i, j, :) = energy2(i, j, :)& & +comp_ice(i, j, :)*c_p*Temp(i, j, KMAX+1) energy(i, j, :) = energy(i, j, :)+comp_ice(i, j, :) cycle end if ! 変化しない場合 if((Hsfc1(i, j) == 0.0D0).and.(Hsfc2(i, j) == 0.0D0))then comp_ice(i, j, :) = 0.0D0 comp_Ms(i, j) = 0.0D0 energy(i, j, :) = 0.0D0 energy2(i, j, :) = 0.0D0 cycle end if do k = k_str, 1, -1 ! かん養による氷と元からある氷がある格子 if(Hsfc2(i, j) * wgt(k)* rho * dx >= mass_Ms(i, j))then k_str = k - 1 ! 次の氷から始める ! かん養, 残りの氷 comp_Ms(i, j) = mass_Ms(i, j) if((i==2).and.(j ==2))then write(6, *)'k=', k end if energy2(i, j, k) = energy2(i, j, k)& + mass_Ms(i, j) * c_p * temp(i, j, k+1) ! k 層目の残り rest = Hsfc2(i, j)*wgt(k)*rho*dx - mass_Ms(i, j) !元からある氷の総エネルギー do kk = kk_str, 1, -1 ! 時刻 t+1 の時の k 層中に含まれる時刻 t の時の層番号kk if(Hsfc1(i, j)*wgt(kk)*rho*dx > rest)then kk_str = kk-1 ! kk まで調べた comp_ice(i, j, k)= comp_ice(i, j, k) & & + rest energy2(i, j, k) = energy2(i, j, k) & & + rest * c_p * Temp(i, j, kk) mass_Ms = Hsfc1(i, j)*wgt(kk)*rho*dx - rest exit else ! 質量が区別つかなくなるので注意 comp_ice(i, j, k)= comp_ice(i, j, k) & & + Hsfc1(i, j)*wgt(kk)*rho*dx energy2(i, j, k) = energy2(i, j, k) & & + Hsfc1(i, j)*wgt(kk)*rho*dx * c_p * Temp(i, j, kk) rest = rest - Hsfc1(i, j) * wgt(kk)*rho*dx end if end do energy(i, j, k) = comp_Ms(i, j) + comp_ice(i, j, k) else ! かん養によるエネルギーの増加 comp_ice(i, j, k) = Hsfc2(i, j)*wgt(k)*rho*dx mass_Ms(i, j) = mass_Ms(i, j) - comp_ice(i, j, k) energy(i, j, k) = energy(i, j, k)+comp_ice(i, j, k) energy2(i, j, k) = energy2(i, j, k)+comp_ice(i, j, k) * c_p * temp(i, j, kk+1) if((i==2).and.(j ==2))then write(6, *)'k=', k , kk, temp(i, j, kk+1) end if end if ! if((i==2).and.(j==2))then ! write(6, *)k, 'energy', energy2(i, j, k)/rho/dx ! write(6, *)k, 'energy', energy2(i, j, k)/rho/dx/c_p/Ts(i, j) ! end if end do do k = 1, KMAX col_mass(i, j) = col_mass(i,j)+energy(i, j, k) end do ! if((i==2).and.(j==2))then ! write(6, *)'total energy', col_mass(i,j)/rho/dx ! write(6 ,*)'total energy', sum(energy2(i, j, :))/c_p/Ts(i, j)/rho/dx ! write(6, *)'total energy', ttene(i, j)/c_p/Ts(i, j)/rho/dx ! end if end do end do end subroutine adjust_energy_Ms_ver3 end module T_pack