!======================================= ! 2D cumulus model - kaminari ! - subroutine non_dim_pi ! ! Author : TAKAHASHI Koko ! Date : 2003/12/09 最終更新 ! 2003/11/19 新規作成 ! Note : SSL2 使用 ! 記号の意味/表記 : ノート => ここ ! 対角成分 : A => d(k) ! 上副対角成分 : B => spd(k) ! 下副対角成分 : C => sbd(k) ! 定数ベクトル(行列) : D => b(k) ! !======================================= !--- 1,2,3行目入力変数, 「k」と4行目入出力変数 subroutine non_dim_pi(its,im,km,cp,cv,rd,alp,beta,dx,dz, & & dts,pi_bs,vptemp_bs,dens_bs,u,w,w_adv,& & w_buoy,w_trb,g_sqrt,g_sqrt_z,g13xz, & & div_z_w, & & pi) implicit none integer(8), intent(in) :: its integer(8), intent(in) :: im, km real(8), intent(in) :: cp,cv,rd,alp,beta real(8), intent(in) :: dx, dz, dts real(8), intent(in) :: pi_bs(-2:im+2,-2:km+2) real(8), intent(in) :: vptemp_bs(-2:im+2,-2:km+2) real(8), intent(in) :: dens_bs(-2:im+2,-2:km+2) real(8), intent(in) :: u(-2:im+2,-2:km+2) real(8), intent(in) :: w(-2:im+2,-2:km+3) real(8), intent(in) :: w_adv(-2:im+2,-2:km+3) real(8), intent(in) :: w_buoy(-2:im+2,-2:km+3) real(8), intent(in) :: w_trb(-2:im+2,-2:km+3) real(8), intent(in) :: g_sqrt(-2:im+2,-2:km+2) real(8), intent(in) :: g_sqrt_z(-2:im+2,-2:km+2) real(8), intent(in) :: g13xz(-2:im+2,-2:km+2) real(8), intent(in) :: div_z_w(0:im+1,0:km+1) real(8), intent(inout) :: pi(-2:im+2,-2:km+2) integer(8) :: i, k real(8) :: c_bs2(0:km) ! 基本場での音速 real(8) :: ssl_fctr(0:km) ! 行列成分の係数 real(8) :: ssl_fctr2(-1:km) ! 行列成分の係数 !--- SSL2 行列計算用の変数 real(8) :: d(0:km) ! 対角成分 real(8) :: sbd(0:km-1) ! 下副対角成分 real(8) :: spd(0:km-1) ! 上副対角成分 real(8) :: b(0:km) ! 定数ベクトル(行列) ! SSL2 ライブラリ用パラメータ real(8) :: vw(0:km) integer(8) :: isw, is, ip(0:km), icon !--- 連立1次方程式を解くためのパラメータ isw = 1 do i = 0,im !--------------------- do k = 0,km !--- 音速の2乗 c_bs2(k) = cp*rd/cv*pi_bs(i,k)*vptemp_bs(i,k) !--- 行列成分の係数 その1 ssl_fctr(k) = c_bs2(k)/g_sqrt(i,k) & & *(beta*dts/(vptemp_bs(i,k)*dz))**2.0d0 & & /dens_bs(i,k) end do !--------------------- !--- 行列成分の係数 その2 do k = -1,km ssl_fctr2(k) = dens_bs(i,k) & & *vptemp_bs(i,k)**2.0d0/g_sqrt(i,k) end do !--------------------- !--- 圧力方程式を解くための行列: 対角部分 (A) ! k = 0 d(0) = 1.0d0 + ssl_fctr(0)*ssl_fctr2(0) do k = 1,km-1 d(k) = 1.0d0 & & + ssl_fctr(k)*(ssl_fctr2(k) + ssl_fctr2(k-1)) end do ! k = km d(km) = 1.0d0 + ssl_fctr(km)*ssl_fctr2(km-1) !--------------------- do k = 0,km-1 !--- 圧力方程式を解くための行列: 上副対角部分 (B) spd(k) = - ssl_fctr(k)*ssl_fctr2(k) !--- 圧力方程式を解くための行列: 下副対角部分 (C) sbd(k) = - ssl_fctr(k)*ssl_fctr2(k-1) end do !--------------------- !--- 定数ベクトル (行列) (D) b(0) = pi(i,0) & & - (1.0d0 - beta)/g_sqrt(i,0)*c_bs2(0)*dts & & /(cp*dens_bs(i,0)*vptemp_bs(i,0)**2.0d0) & & *( & & (dens_bs(i,1) + dens_bs(i,0))/2.0d0 & & *(vptemp_bs(i,1) + vptemp_bs(i,0))/2.0d0 & & *w(i,1) & & - (dens_bs(i,0) + dens_bs(i,-1))/2.0d0 & & *(vptemp_bs(i,0) + vptemp_bs(i,-1))/2.0d0 & & *w(i,0) & & )/dz & & - c_bs2(0)*dts/(cp*vptemp_bs(i,0)) & & *( & & (u(i+1,0) - u(i,0))/dx & & + g13xz(i,0) & & *(u(i+1,1) + u(i,1) - u(i+1,-1) - u(i,-1)) & & /(4.0d0*dz) & & ) & & - beta/g_sqrt(i,0)*c_bs2(0)*dts & & /(cp*dens_bs(i,0)*vptemp_bs(i,0)**2.0d0) & & *( & & (dens_bs(i,1) + dens_bs(i,0))/2.0d0 & & *(vptemp_bs(i,1) + vptemp_bs(i,0))/2.0d0 & & *( & & w(i,1) & & - cp*dts/g_sqrt_z(i,1) & & *(vptemp_bs(i,1) + vptemp_bs(i,0))/2.0d0 & & *( & & - alp*div_z_w(i,1) & & + (1.0d0 - beta)*(pi(i,1) - pi(i,0))/dz & & ) & & - (- w_adv(i,1) + w_buoy(i,1) + w_trb(i,1)) & & *dts & & ) & & - (dens_bs(i,0) + dens_bs(i,-1))/2.0d0 & & *(vptemp_bs(i,0) + vptemp_bs(i,-1))/2.0d0 & & *( & & w(i,0) & & - cp*dts/g_sqrt_z(i,0) & & *(vptemp_bs(i,0) + vptemp_bs(i,-1))/2.0d0& & *( & & - alp*div_z_w(i,0) & & + (1.0d0 - beta) & & *(pi(i,0) - pi(i,-1))/dz & & ) & & - ( & & - w_adv(i,0) + w_buoy(i,0) + w_trb(i,0) & & )*dts & & ) & & )/dz & & - beta/g_sqrt(i,0)*c_bs2(0)*dts**2.0d0 & & /(dens_bs(i,0)*vptemp_bs(i,0)**2.0d0*dz) & & *cp & & *(vptemp_bs(i,0) + vptemp_bs(i,-1)/2.0d0) & & **2.0d0/g_sqrt_z(i,0) & & *( & & alp*div_z_w(i,0) & & - (1.0d0 - beta)*(pi(i,0) - pi(i,-1))/dz & & + g_sqrt_z(i,0) & & /cp & & /(vptemp_bs(i,0) + vptemp_bs(i,-1)/2.0d0) & & *(- w_adv(i,0) + w_buoy(i,0)+ w_trb(i,0)) & & ) do k = 1,km-1 b(k) = pi(i,k) & & - (1.0d0 - beta) & & /g_sqrt(i,k)*c_bs2(k)*dts & & /(cp*dens_bs(i,k)*vptemp_bs(i,k)**2.0d0) & & *( & & (dens_bs(i,k+1) + dens_bs(i,k))/2.0d0 & & *(vptemp_bs(i,k+1) + vptemp_bs(i,k))/2.0d0 & & *w(i,k+1) & & - (dens_bs(i,k) + dens_bs(i,k-1))/2.0d0 & & *(vptemp_bs(i,k) + vptemp_bs(i,k-1))/2.0d0 & & *w(i,k) & & )/dz & & - c_bs2(k)*dts/(cp*vptemp_bs(i,k)) & & *( & & (u(i+1,k) - u(i,k))/dx & & + g13xz(i,k) & & *( & & u(i+1,k+1) + u(i,k+1) & & - u(i+1,k-1) - u(i,k-1) & & )/(4.0d0*dz) & & ) & & - beta/g_sqrt(i,k)*c_bs2(k)*dts & & /(cp*dens_bs(i,k)*vptemp_bs(i,k)**2.0d0) & & *( & & (dens_bs(i,k+1) + dens_bs(i,k))/2.0d0 & & *(vptemp_bs(i,k+1) + vptemp_bs(i,k))/2.0d0 & & *( & & w(i,k+1) & & - cp*dts/g_sqrt_z(i,k+1) & & *(vptemp_bs(i,k+1) + vptemp_bs(i,k)) & & /2.0d0 & & *( & & - alp*div_z_w(i,k+1) & & + (1.0d0 - beta) & & *(pi(i,k+1) - pi(i,k))/dz & & ) & & - ( & & - w_adv(i,k+1) + w_buoy(i,k+1) & & + w_trb(i,k+1) & & )*dts & & ) & & - (dens_bs(i,k) + dens_bs(i,k-1))/2.0d0 & & *(vptemp_bs(i,k) + vptemp_bs(i,k-1))/2.0d0 & & *( & & w(i,k) & & - cp*dts/g_sqrt_z(i,k) & & *(vptemp_bs(i,k) + vptemp_bs(i,k-1)) & & /2.0d0 & & *( & & - alp*div_z_w(i,k-1) & & + (1.0d0 - beta) & & *(pi(i,k) - pi(i,k-1))/dz & & ) & & - ( & & - w_adv(i,k) + w_buoy(i,k) & & + w_trb(i,k) & & )*dts & & ) & & )/dz end do b(km) = pi(i,km) & & - (1.0d0 - beta)/g_sqrt(i,km)*c_bs2(km)*dts & & /(cp*dens_bs(i,km)*vptemp_bs(i,km)**2.0d0) & & *( & & (dens_bs(i,km+1) + dens_bs(i,km))/2.0d0 & & *(vptemp_bs(i,km+1) + vptemp_bs(i,km)) & & /2.0d0 & & *w(i,km+1) & & - (dens_bs(i,km) + dens_bs(i,km-1))/2.0d0 & & *(vptemp_bs(i,km) + vptemp_bs(i,km-1)) & & /2.0d0 & & *w(i,km) & & )/dz & & - c_bs2(km)*dts/(cp*vptemp_bs(i,km)) & & *( & & (u(i+1,km) - u(i,km))/dx & & + g13xz(i,km) & & *( & & u(i+1,km+1) + u(i,km+1) & & - u(i+1,km-1) - u(i,km-1) & & )/(4.0d0*dz) & & ) & & - beta/g_sqrt(i,km)*c_bs2(km)*dts & & /(cp*dens_bs(i,km)*vptemp_bs(i,km)**2.0d0) & & *( & & (dens_bs(i,km+1) + dens_bs(i,km))/2.0d0 & & *(vptemp_bs(i,km+1) + vptemp_bs(i,km)) & & /2.0d0 & & *( & & w(i,km+1) & & - cp*dts/g_sqrt_z(i,km+1) & & *(vptemp_bs(i,km+1) + vptemp_bs(i,km)) & & /2.0d0 & & *( & & - alp*div_z_w(i,km+1) & & + (1.0d0 - beta) & & *(pi(i,km+1) - pi(i,km))/dz & & ) & & - ( & & - w_adv(i,km+1) + w_buoy(i,km+1) & & + w_trb(i,km+1) & & )*dts & & ) & & - (dens_bs(i,km) + dens_bs(i,km-1))/2.0d0 & & *(vptemp_bs(i,km) + vptemp_bs(i,km-1)) & & /2.0d0 & & *( & & w(i,km) & & - cp*dts/g_sqrt_z(i,km) & & *(vptemp_bs(i,km) + vptemp_bs(i,km-1))& & /2.0d0 & & *( & & - alp*div_z_w(i,km) & & + (1.0d0 - beta) & & *(pi(i,km) - pi(i,km-1))/dz & & ) & & - ( & & - w_adv(i,km) + w_buoy(i,km) & & + w_trb(i,km) & & )*dts & & ) & & )/dz & & + beta/g_sqrt(i,km)*c_bs2(km)*dts**2.0d0 & & /(dens_bs(i,km)*vptemp_bs(i,km)**2.0d0*dz) & & *cp & & *(vptemp_bs(i,km+1)+vptemp_bs(i,km)/2.0d0) & & **2.0d0/g_sqrt_z(i,km+1) & & *( & & alp*div_z_w(i,km+1) & & - (1.0d0 - beta)*(pi(i,km+1) - pi(i,km))/dz& & + g_sqrt_z(i,km+1) & & /cp & & /(vptemp_bs(i,km+1)+vptemp_bs(i,km)/2.0d0)& & *( & & - w_adv(i,km+1) + w_buoy(i,km+1) & & + w_trb(i,km+1) & & ) & & ) !--------------------- !--- 無次元圧力 call dltx(sbd,d,spd,km+1,b,0.0d0,isw,is,ip,vw,icon) do k = 0,km pi(i,k) = b(k) !--- エラーチェック if (icon /= 0) then write(*,*) 'its =', its,', i =',i,', k =',k, & & ', Condition code:', icon stop ! else ! write(*,*) 'its =', its,', i =',i,', k =',k, & ! & ', Condition code:', icon end if end do end do end subroutine non_dim_pi