!---------------------------------------------------------------------- ! Copyright (c) 2004--2005 YAMADA YUKIKO. All rights reserved. !---------------------------------------------------------------------- ! Sample program for SPML(based on gt4f90io and ISPACK) ! ! ! Solving equatorial beta plane wave-CISK equations ! ! history: 2004/06/23 YAMADA, Y. 作成 ! 2004/06/27 YAMADA, Y. ! 2005/03/31 YAMADA, Y. ! program equatorial_beta_wavecisk use esc_module use gt4_history use gt4f90io implicit none !---- 空間解像度設定 ---- integer, parameter :: km=16, lm=16, mm = 3 ! 切断波数の設定(X,Y,Z) integer, parameter :: im=64, jm=32, vm = mm+1 ! 格子点の設定(X,Y,Z) ! integer, parameter :: km=64, lm=64, mm = 3 ! 切断波数の設定(X,Y,Z) ! integer, parameter :: im=200, jm=100, vm = mm+1 ! 格子点の設定(X,Y,Z) !---- 加熱のパラメタリゼーション設定 ---- ! wave-CISK パラメタリゼーション: TRUE=上昇加熱型, FALSE=線形加熱型 logical :: positive_only_heat_parameta = .FALSE. !---- 変数 ---- real(8) :: yxz_U(0:jm,0:im-1,vm) ! 格子データ(U(t+Δt)) real(8) :: yxz_V(0:jm,0:im-1,vm) ! 格子データ(V(t+Δt)) real(8) :: yxz_H(0:jm,0:im-1,vm) ! 格子データ(H(t+Δt)) real(8) :: yxm_U(0:jm,0:im-1,mm) ! 格子データ(U(t+Δt)) real(8) :: yxm_V(0:jm,0:im-1,mm) ! 格子データ(V(t+Δt)) real(8) :: yxm_H(0:jm,0:im-1,mm) ! 格子データ(H(t+Δt)) real(8) :: ecm_U(-km:km,0:lm,mm) ! スペクトル(U(t)) real(8) :: esm_V(-km:km,lm,mm) ! スペクトル(V(t)) real(8) :: ecm_H(-km:km,0:lm,mm) ! スペクトル(H(t)) real(8) :: ecm_U_B(-km:km,0:lm,mm) ! スペクトル(U(t-Δt)) real(8) :: esm_V_B(-km:km,lm,mm) ! スペクトル(V(t-Δt)) real(8) :: ecm_H_B(-km:km,0:lm,mm) ! スペクトル(H(t-Δt)) real(8) :: ec_QH(-km:km,0:lm) ! スペクトル(熱源分布) real(8) :: ec_QH_A(-km:km,0:lm) ! スペクトル(熱源分布) real(8) :: gamma(mm) ! 鉛直モード固有値^0.5 real(8) :: eig(mm,0:vm) ! 鉛直モード固有関数 real(8) :: eta(mm) ! 加熱の鉛直分布 !---- 座標変数など ---- real(8), parameter :: xmin=-10.0, xmax=10.0 ! 領域範囲(X 方向) ! real(8), parameter :: xmin=-25.0, xmax=25.0 ! 領域範囲(X 方向) real(8), parameter :: ymin=-5.0, ymax=5.0 ! 領域範囲(Y 方向) !---- 物理パラメター ---- real(8), parameter :: Grav=1.0 ! 重力加速度 real(8), parameter :: Beta=1.0 ! 赤道βパラメタ real(8), parameter :: H0 = 1.0 ! 平均の深さ real(8), parameter :: Rm = 0.0 ! レイリー摩擦係数 real(8), parameter :: Rt = 0.0 ! 平均ニュートン冷却係数 real(8), parameter :: Num = 0.02 ! 粘性拡散係数 real(8), parameter :: Nut = 0.02 ! 熱拡散係数 real(8), parameter :: Q0=1.0 ! 初期熱源分布の振幅 real(8), parameter :: AX0=0.5 ! 初期熱源分布の半値幅 ! real(8), parameter :: AX0=0.5*4*2.5 ! 初期熱源分布の半値幅 real(8), parameter :: AY0=0.5 ! 初期熱源分布の半値幅 !---- 時間積分パラメタ ---- real(8), parameter :: dt=2e-2 ! 時間ステップ間隔 real(8), parameter :: tfil=0.05 ! 時間フィルター係数 integer, parameter :: nt=1500, ndisp=20 ! 時間積分数, 表示ステップ ! integer, parameter :: nt=3000, ndisp=20 ! 時間積分数, 表示ステップ !---- 作業配列 ---- integer :: i, j, it, m !---- スペクトル設定 ---- call esc_initial(im,jm,km,lm,xmin,xmax,ymin,ymax) ! ISPACK 初期化 call vmode_initial( eig, gamma ) ! 鉛直モード固有値固有関数の設定 !-------- 初期値設定 ----------- ! 初期値は静止流体 ! yxz_H = H0 ! 平均の深さ yxz_H = 0.0 ! 平均の深さ yxz_U = 0.0 yxz_V = 0.0 ! yxm_H = H0 ! 平均の深さ yxm_H = 0.0 ! 平均の深さ yxm_U = 0.0 yxm_V = 0.0 ! 加熱の鉛直モードの設定 eta(1) = 1.5 ; eta(2) = -1.5 do m=1,mm eta(m) = eta(m) / eig(m,1) * gamma(m) / vm * (2**0.5) enddo ! 初期熱源分布 (第一モードのみ加熱を与える) do i=0,im-1 do j=0,jm yxm_H(j,i,1) = - Q0 & * exp( -yx_Y(j,i)**2/AY0 ) & * exp( -yx_X(j,i)**2/AX0 ) enddo enddo ! 水平スペクトル変換 do m = 1, mm ecm_U(:,:,m) = ec_yx(yxm_U(:,:,m) ) esm_V(:,:,m) = es_yx(yxm_V(:,:,m) ) ecm_H(:,:,m) = ec_yx(yxm_H(:,:,m) ) enddo !------------------------------------------------------------ call output_gtool4_init ! ヒストリー初期化 call output_gtool4 ! 初期値出力 !----------------------- 積分 main ---------------------------- ! 1 ステップ目はオイラー法で時間積分 ecm_H_B = ecm_H ecm_U_B = ecm_U esm_V_B = esm_V do it=1,nt ! 鉛直流 (加熱項) 初期化 ec_QH_A = 0.0 do m = 1, mm ! 時間積分 call time_integral( gamma(m), eta(m), ec_QH, & ecm_U(:,:,m), esm_V(:,:,m), ecm_H(:,:,m), & ecm_U_B(:,:,m), esm_V_B(:,:,m), ecm_H_B(:,:,m) ) ! 鉛直流 (加熱項) の計算 ec_QH_A = ec_QH_A & - (ec_Dx_ec(ecm_U(:,:,m)) + ec_Dy_es(esm_V(:,:,m))) enddo ec_QH = ec_QH_A ! 上昇流域のみ加熱の場合 if (positive_only_heat_parameta .AND. .TRUE. ) then call positive_only_heat(ec_QH) endif ! 出力 if(mod(it,ndisp) .eq. 0) then call output_gtool4 endif enddo call output_gtool4_close stop contains !----------- 時間積分[leapflog and timefilter] ------------ subroutine time_integral( gamma_m, eta_m, ec_QH, & ec_U,es_V,ec_H,ec_U_B,es_V_B,ec_H_B) real(8) :: ec_U_A(-km:km,0:lm) ! スペクトル(U(t+Δt)) real(8) :: es_V_A(-km:km,lm) ! スペクトル(V(t+Δt)) real(8) :: ec_H_A(-km:km,0:lm) ! スペクトル(H(t+Δt)) real(8),intent(inout):: ec_U(-km:km,0:lm) ! スペクトル(U(t)) real(8),intent(inout):: es_V(-km:km,lm) ! スペクトル(V(t)) real(8),intent(inout):: ec_H(-km:km,0:lm) ! スペクトル(H(t)) real(8),intent(inout):: ec_U_B(-km:km,0:lm) ! スペクトル(U(t-Δt)) real(8),intent(inout):: es_V_B(-km:km,lm) ! スペクトル(V(t-Δt)) real(8),intent(inout):: ec_H_B(-km:km,0:lm) ! スペクトル(H(t-Δt)) real(8),intent(in) :: gamma_m ! 鉛直モード固有値^0.5 real(8),intent(in) :: eta_m ! 鉛直加熱分布 real(8),intent(in) :: ec_QH(-km:km,0:lm) ! スペクトル(熱源分布) ec_U_A = ec_U_B & & + dt * 2.0 * ( & !& - ec_yx(yx_ec(ec_U)*yx_ec(ec_Dx_ec(ec_U))) & !& - ec_yx(yx_es(es_V)*yx_es(es_Dy_ec(ec_U))) & & + ec_yx( Beta * yx_Y * (yx_es(es_V))) & & - Grav * ec_dx_ec(ec_H) & & + Num * ( ec_Dx_ec(ec_Dx_ec(ec_U_B)) & & + ec_Dy_es(es_Dy_ec(ec_U_B)) ) & & - Rm * ec_U & & ) es_V_A = es_V_B & & + dt * 2.0 * ( & !& - es_yx(yx_ec(ec_U)*yx_es(es_Dx_es(es_V))) & !& - es_yx(yx_es(es_V)*yx_ec(ec_Dy_es(es_V))) & & - es_yx( Beta * yx_Y * (yx_ec(ec_U))) & & - Grav * es_dy_ec(ec_H) & & + Num * ( es_Dx_es(es_Dx_es(es_V_B)) & & + es_Dy_ec(ec_Dy_es(es_V_B)) ) & & -Rm * es_V & & ) ec_H_A = ec_H_B & & + dt * 2.0 * ( & & + 1.0 / gamma_m **2 * ( & & - H0*(ec_dx_ec(ec_U) + ec_dy_es(es_V)) & !& - ec_dx_ec(ec_yx(yx_ec(ec_H)*yx_ec(ec_U))) & !& - ec_dy_es(es_yx(yx_ec(ec_H)*yx_es(es_V))) & & - eta_m * ec_QH & & ) & & + Nut * ( ec_Dx_ec(ec_Dx_ec(ec_H_B)) & & + ec_Dy_es(es_Dy_ec(ec_H_B)) ) & & - Rt * ec_H & & ) ! 時間フィルター ec_U_B = (1.0 - 2.0 * tfil) * ec_U + tfil * (ec_U_B + ec_U_A) ec_U = ec_U_A ec_H_B = (1.0 - 2.0 * tfil) * ec_H + tfil * (ec_H_B + ec_H_A) ec_H = ec_H_A es_V_B = (1.0 - 2.0 * tfil) * es_V + tfil * (es_V_B + es_V_A) es_V = es_V_A end subroutine time_integral ! 上昇流域のみ加熱の場合 subroutine positive_only_heat(ec_QH) real(8),intent(inout):: ec_QH(-km:km,0:lm) ! スペクトル(熱源分布) real(8) :: yx_QH(0:jm,0:im-1) ! 格子データ(熱源分布) yx_QH = yx_ec(ec_QH) do i=0,im-1 do j=0,jm yx_QH(j,i) = max(0,yx_QH(j,i)) enddo enddo ec_QH = ec_yx(yx_QH) end subroutine positive_only_heat ! 鉛直モードの固有値, 固有関数の設定 subroutine vmode_initial( eig, gamma ) real(8),intent(out):: gamma(mm) ! 鉛直モード固有値^0.5 real(8),intent(out):: eig(mm,0:vm) ! 鉛直モード固有関数 integer :: j, m real(8) :: pi ! 円周率 pi = atan(1.0D0)*4.0D0 ! 固有関数 上下端で 0 eig(:,0) = 0.0 ; eig(:,vm) = 0.0 do m=1,mm ! 固有値^0.5 の計算 gamma(m) = vm * ( 2 - 2 * cos( m*pi/(mm+1) ) )**0.5 ! 固有関数 (規格化してない) 計算 do j=1,vm-1 eig(m,j) = sin ( (m*j*pi)/(mm+1) ) enddo enddo end subroutine vmode_initial ! 鉛直モード逆変換 (u,v,phi) subroutine inv_vmode_transpote(yxm_data, yxz_data) real(8),intent(in) :: yxm_data(0:jm,0:im-1,mm) ! 鉛直モードデータ real(8),intent(out):: yxz_data(0:jm,0:im-1,vm) ! 鉛直格子データ yxz_data = 0.0 do j=1,vm do m=1,mm yxz_data(:,:,j) = yxz_data(:,:,j) & + ( eig(m,j) - eig(m,j-1) ) & * yxm_data(:,:,m) * vm / gamma(m) / (2**0.5) enddo enddo end subroutine inv_vmode_transpote !----------------------- gtool4 設定 ---------------------------- subroutine output_gtool4_init real(8) :: vm_ary(vm) ! 鉛直層の番号付 do j = 1,vm vm_ary(j) = j * 1.0 enddo ! ヒストリー作成 call HistoryCreate( & file='wavecisk.nc', & title='Equatorial beta-plane wave-CISK model ', & source='A sample program of ISPACK and gtool_history/gtool4', & institution='GFD_Dennou Club spmodel project', & dims=(/'x','y','z','t'/), dimsizes=(/im,jm+1,vm,0/), & longnames=(/'X-coordinate','Y-coordinate', & 'Z-level ','time '/), & units=(/'1','1','1','1'/), & origin=0.0, interval=real( (ndisp)*dt ) & ) call HistoryPut('x',x_X) ! 変数出力 call HistoryAddattr('x','topology','circular') ! 周期属性 call HistoryAddattr('x','modulo',xmax-xmin) ! 周期属性 call HistoryPut('y',y_Y) ! 変数出力 call HistoryPut('z',vm_ary) ! 変数出力 call HistoryAddVariable( & ! 変数定義 varname='uvel', dims=(/'x','y','z','t'/), & longname='X-velocity', units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='vvel', dims=(/'x','y','z','t'/), & longname='Y-velocity', units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='height', dims=(/'x','y','z','t'/), & longname='height', units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='wvel', dims=(/'x','y','t'/), & longname='Z-velocity_(3/2)', units='1', xtype='double') end subroutine output_gtool4_init subroutine output_gtool4 real(8) :: yxz_tr_U(0:im-1,0:jm,vm) ! 格子データ (X-Y 軸入れ換え) real(8) :: yxz_tr_V(0:im-1,0:jm,vm) ! 格子データ (X-Y 軸入れ換え) real(8) :: yxz_tr_H(0:im-1,0:jm,vm) ! 格子データ (X-Y 軸入れ換え) write(6,*) 'it = ',it ! 水平スペクトル逆変換 do m = 1, mm yxm_U(:,:,m) = yx_ec(ecm_U(:,:,m)) yxm_V(:,:,m) = yx_es(esm_V(:,:,m)) yxm_H(:,:,m) = yx_ec(ecm_H(:,:,m)) enddo ! 鉛直モード逆変換 call inv_vmode_transpote(yxm_U, yxz_U) call inv_vmode_transpote(yxm_V, yxz_V) call inv_vmode_transpote(yxm_H, yxz_H) ! X-Y 軸入れ換え do m = 1, vm yxz_tr_U(:,:,m) = transpose(yxz_U(:,:,m)) yxz_tr_V(:,:,m) = transpose(yxz_V(:,:,m)) yxz_tr_H(:,:,m) = transpose(yxz_H(:,:,m)) enddo call HistoryPut('t',real( it*dt ) ) call HistoryPut('uvel', yxz_tr_U ) call HistoryPut('vvel', yxz_tr_V ) call HistoryPut('height', yxz_tr_H ) call HistoryPut('wvel', transpose( yx_ec(ec_QH_A)) ) end subroutine output_gtool4 subroutine output_gtool4_close call HistoryClose end subroutine output_gtool4_close end program equatorial_beta_wavecisk