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