!======================================= ! 2D cumulus model - kaminari ! - subroutine eddymx_sub ! ! Author : TAKAHASHI Koko ! Date : 2003/11/25 最終更新 ! 2003/11/05 新規作成 ! Note : 1ステップ目のサブグリッド ! 運動エネルギーを求めてから ! 渦粘性係数を求める ! 出力は渦粘性係数 ! !======================================= !--- 1,2,3,4,5,6行目入力変数, 7行目入出力変数 subroutine eddymx_sub(i,k,im,km,nuh,nuv,eps,lv,rd,cp,cm,ml, & & grv,dtb,dx,dz, & & g_sqrt,g_sqrt_z,g13x,g13z,g13xz, & & ptemp_bs,u_old,w_old,omg_old,pi_old, & & ptemp_old,qv_old,qc_old, & & km_sub_old,km_sub_old2, & & km_sub) implicit none integer(8), intent(in) :: i, k integer(8), intent(in) :: im, km real(8), intent(in) :: nuh,nuv,eps, lv, rd, cp, cm, ml real(8), intent(in) :: grv, dtb real(8), intent(in) :: dx, dz real(8), intent(in) :: g_sqrt(-2:im+1,-2:km+1) real(8), intent(in) :: g_sqrt_z(-2:im+1,-2:km+1) real(8), intent(in) :: g13x(-2:im+1,-2:km+1) real(8), intent(in) :: g13z(-2:im+1,-2:km+1) real(8), intent(in) :: g13xz(-2:im+1,-2:km+1) real(8), intent(in) :: ptemp_bs(-2:im+1,-2:km+1) real(8), intent(in) :: u_old(-2:im+1,-2:km+1) real(8), intent(in) :: w_old(-2:im+1,-2:km+1) real(8), intent(in) :: omg_old(-2:im+1,-2:km+1) real(8), intent(in) :: pi_old(-2:im+1,-2:km+1) real(8), intent(in) :: ptemp_old(-2:im+1,-2:km+1) real(8), intent(in) :: qv_old(-2:im+1,-2:km+1) real(8), intent(in) :: qc_old(-2:im+1,-2:km+1) real(8), intent(in) :: km_sub_old(-2:im+1,-2:km+1) real(8), intent(in) :: km_sub_old2(-2:im+1,-2:km+1) real(8), intent(inout) :: km_sub(-2:im+1,-2:km+1) real(8) :: temp_old(-2:im+1,-2:km+1) real(8) :: a(-2:im+1,-2:km+1) real(8) :: ptemp_e(-2:im+1,-2:km+1) real(8) :: ql(-2:im+1,-2:km+1) real(8) :: e_sub(-2:im+1,-2:km+1) real(8) :: e_sub_dif(-2:im+1,-2:km+1) real(8) :: e_sub_adv(-2:im+1,-2:km+1) real(8) :: e_sub_old(-2:im+1,-2:km+1) real(8) :: e_sub_old2(-2:im+1,-2:km+1) !--- サブグリッド運動エネルギー e_sub(i,k) = (km_sub(i,k)/(cm*ml))**2.0d0 e_sub_old(i,k) = (km_sub_old(i,k)/(cm*ml))**2.0d0 e_sub_old2(i,k) = (km_sub_old2(i,k)/(cm*ml))**2.0d0 !--- サブグリッド運動エネルギーの数値粘性項 e_sub_dif(i,k) = nuh & & *( & & e_sub_old(i+1,k) - 2.0d0*e_sub_old(i,k)& & + e_sub_old(i-1,k) & & )/dx**2.0d0 & & + nuv & & *( & & e_sub_old(i,k+1) - 2.0d0*e_sub_old(i,k)& & + e_sub_old(i,k-1) & & )/dz**2.0d0 !--- サブグリッド運動エネルギーの移流項 e_sub_adv(i,k) = (u_old(i+1,k) + u_old(i,k))/2.0d0 & & *( & & e_sub_old(i+1,k+1) + e_sub_old(i+1,k-1)& & - e_sub_old(i-1,k+1) & & - e_sub_old(i-1,k-1) & & )/(4.0d0*dx) & & + (omg_old(i,k+1) + omg_old(i,k))/2.0d0 & & *( & & e_sub_old(i+1,k+1) & & + e_sub_old(i-1,k+1) & & - e_sub_old(i+1,k-1) & & - e_sub_old(i-1,k-1) & & )/(4.0d0*dz) & & + e_sub_dif(i,k) !--- 温位 => 温度 temp_old(i,k) = ptemp_old(i,k)*pi_old(i,k) !--- 係数 A a(i,k) = 1.0d0/ptemp_bs(i,k) & & *( & & 1.0d0 & & + 1.61d0*eps*lv*qv_old(i,k)/(rd*temp_old(i,k)) & & )/(1.0d0 & & + ( & & eps*(lv**2.0d0)*qv_old(i,k) & & /(cp*rd*temp_old(i,k)**2.0d0) & & ) & & ) !--- 相当温位 ptemp_e(i,k) = ptemp_old(i,k) & & *(1.0d0 + lv*qv_old(i,k)/(cp*temp_old(i,k))) !--- 水蒸気と雲水混合比の和 ql(i,k) = qv_old(i,k) + qc_old(i,k) !-- サブグリッド運動エネルギー e_sub(i,k) = e_sub_old2(i,k) & & + 2.0d0*dtb & & *(- e_sub_adv(i,k) & & + 3.0d0*grv*cm*ml*sqrt(e_sub_old(i,k)) & & /g_sqrt(i,k) & & *(- a(i,k) & & *( & & ptemp_e(i+1,k+1) + ptemp_e(i-1,k+1) & & - ptemp_e(i+1,k-1) & & - ptemp_e(i-1,k-1) & & )/(4.0d0*dz) & & + ( & & ql(i+1,k+1) + ql(i-1,k+1) & & - ql(i+1,k-1) - ql(i-1,k-1) & & )/(4.0d0*dz) & & ) & & + 2.0d0*cm*ml*sqrt(e_sub_old(i,k)) & & *( & & ((u_old(i+1,k) - u_old(i,k))/dx & & + g13xz(i,k) & & *(u_old(i+1,k+1) + u_old(i-1,k+1) & & - u_old(i+1,k-1) - u_old(i-1,k-1)& & )/(4.0d0*dz) & & )**2.0d0 & & + ( & & 1.0d0/g_sqrt(i,k) & & *(w_old(i,k+1) - w_old(i,k))/dz & & )**2.0d0 & & ) & & + cm*ml*sqrt(e_sub_old(i,k)) & & *( & & ( & & w_old(i+1,k+1) + w_old(i+1,k-1) & & - w_old(i-1,k+1) - w_old(i-1,k-1) & & )/(4.0d0*dx) & & + g13xz(i,k) & & *(w_old(i,k+1) - w_old(i,k))/dz & & + 1.0d0/g_sqrt(i,k) & & *( & & u_old(i+1,k+1) + u_old(i-1,k+1) & & - u_old(i+1,k-1) - u_old(i-1,k-1)& & )/(4.0d0*dz) & & ) & & + cm*ml*sqrt(e_sub_old(i,k)) & & *( & & ( & & e_sub_old(i+1,k) & & - 2.0d0*e_sub_old(i,k) & & + e_sub_old(i-1,k) & & )/dx**2.0d0 & & + ( & & g13z(i+1,k) & & *( & & e_sub_old(i+1,k+1) & & + e_sub_old(i,k+1) & & - e_sub_old(i+1,k-1) & & - e_sub_old(i,k-1) & & )/(4.0d0*dz) & & - g13z(i,k) & & *( & & e_sub_old(i,k+1) & & + e_sub_old(i-1,k+1) & & - e_sub_old(i,k-1) & & - e_sub_old(i-1,k-1) & & )/(4.0d0*dz) & & )/dx & & + ( & & g13x(i,k+1) & & *( & & e_sub_old(i+1,k+1) & & + e_sub_old(i+1,k) & & - e_sub_old(i-1,k+1) & & - e_sub_old(i-1,k) & & )/(4.0d0*dx) & & - g13x(i,k) & & *( & & e_sub_old(i+1,k) & & + e_sub_old(i+1,k-1) & & - e_sub_old(i-1,k) & & - e_sub_old(i-1,k-1) & & )/(4.0d0*dx) & & )/dz & & + g13xz(i,k) & & *( & & g13x(i+1,k) & & *( & & e_sub_old(i,k+1) & & - e_sub_old(i,k) & & )/dz & & - g13x(i,k) & & *( & & e_sub_old(i,k) & & - e_sub_old(i,k-1) & & )/dz & & )/dz & & + 1.0d0/g_sqrt(i,k) & & *( & & 1.0d0/g_sqrt_z(i,k+1) & & *( & & e_sub_old(i,k+1) & & - e_sub_old(i,k) & & )/dz & & - 1.0d0/g_sqrt_z(i,k) & & *( & & e_sub_old(i,k) & & - e_sub_old(i,k-1) & & )/dz & & )/dz & & ) & & + 2.0d0*cm*ml*sqrt(e_sub_old(i,k)) & & *( & & ( & & ( & & sqrt(e_sub_old(i+1,k)) & & - sqrt(e_sub_old(i-1,k)) & & )/(2.0d0*dx) & & + g13xz(i,k) & & *( & & sqrt(e_sub_old(i,k+1)) & & - sqrt(e_sub_old(i,k-1)) & & )/(2.0d0*dz) & & )**2.0d0 & & + ( & & 1.0d0/g_sqrt_z(i,k) & & *( & & sqrt(e_sub_old(i,k+1)) & & - sqrt(e_sub_old(i,k-1)) & & )/(2.0d0*dz) & & )**2.0d0 & & ) & & - cm/ml*e_sub_old(i,k)**(3.0d0/2.0d0) & & ) !--- 渦混合係数 km_sub(i,k) = cm*ml*sqrt(e_sub(i,k)) write(*,*) i,k,e_sub(i,k),e_sub_old(i,k),e_sub_old2(i,k) write(*,*) e_sub_dif(i,k),e_sub_adv(i,k) end subroutine eddymx_sub