!========================================================================= ! 二次元トレーサー移流モデル: 時間増分計算モジュール ! ! 履歴 1997/10/16 小高正嗣 ! 1997/10/22 小高正嗣 ! 1997/11/13 小高正嗣 ! 1997/11/18 小高正嗣 ! 1997/11/20 小高正嗣 ! 1998/01/20 小高正嗣 ! !========================================================================= MODULE cal_module USE param_module PUBLIC CONTAINS !========================================================================= SUBROUTINE get_dtrc_C2_2D( u, v, TRC, DTRC )! 二次中央差分 REAL,DIMENSION(:,:),INTENT(in) :: u REAL,DIMENSION(:,:),INTENT(in) :: v REAL,DIMENSION(0:,0:),INTENT(in):: TRC REAL,DIMENSION(0:,0:),INTENT(out):: DTRC INTEGER :: dim, xdim, ydim REAL,DIMENSION(:,:),ALLOCATABLE :: FLUX_X! フラックス1 REAL,DIMENSION(:,:),ALLOCATABLE :: FLUX_Y! フラックス2 REAL,DIMENSION(:,:),ALLOCATABLE :: ADVE_X! 移流項1 REAL,DIMENSION(:,:),ALLOCATABLE :: ADVE_Y! 移流項2 REAL,DIMENSION(:,:),ALLOCATABLE :: TRCM_X! 変数(u 格子上) REAL,DIMENSION(:,:),ALLOCATABLE :: TRCM_Y! 変数(v 格子上) xdim = SIZE( TRC, 1 ) - 1 ydim = SIZE( TRC, 2 ) - 1 dim = xdim - 1 ALLOCATE( TRCM_X(xdim,0:ydim), TRCM_Y(0:xdim,ydim) ) ALLOCATE( FLUX_X(xdim,ydim), FLUX_Y(xdim,ydim) ) ALLOCATE( ADVE_X(0:xdim,0:ydim), ADVE_Y(0:xdim,0:ydim) ) CALL get_mval_2D( TRC, TRCM_X, TRCM_Y ) FLUX_X = u * TRCM_X(:,1:) FLUX_Y = v * TRCM_Y(1:,:) ADVE_X(1:dim,1:dim) = - ( FLUX_X(2:,:dim) - FLUX_X(:dim,:dim) )/dx ADVE_Y(1:dim,1:dim) = - ( FLUX_Y(:dim,2:) - FLUX_Y(:dim,:dim) )/dy call bound2( ADVE_X ) call bound2( ADVE_Y ) DTRC = dt * ( ADVE_X + ADVE_Y ) !- 境界上での人工減衰 DTRC(1,:) = - 5.*dt*TRC(1,:) DTRC(dim,:) = - 5.*dt*TRC(dim,:) DTRC(:,1) = - 5.*dt*TRC(:,1) DTRC(:,dim) = - 5.*dt*TRC(:,dim) DEALLOCATE( TRCM_X, TRCM_Y, FLUX_X, FLUX_Y, ADVE_X, ADVE_Y ) END SUBROUTINE get_dtrc_C2_2D !========================================================================= SUBROUTINE get_dtrc_C4_2D( u, v, TRC, DTRC )! 四次中央差分 REAL,DIMENSION(:,:),INTENT(in) :: u REAL,DIMENSION(:,:),INTENT(in) :: v REAL,DIMENSION(0:,0:),INTENT(in):: TRC REAL,DIMENSION(0:,0:),INTENT(out):: DTRC INTEGER :: dim, xdim, ydim REAL,DIMENSION(:,:),ALLOCATABLE :: FLUX_X! フラックス1 REAL,DIMENSION(:,:),ALLOCATABLE :: FLUX_Y! フラックス2 REAL,DIMENSION(:,:),ALLOCATABLE :: ADVE_X! 移流項1 REAL,DIMENSION(:,:),ALLOCATABLE :: ADVE_Y! 移流項2 REAL,DIMENSION(:,:),ALLOCATABLE :: TRCM_X! 変数(u 格子上) REAL,DIMENSION(:,:),ALLOCATABLE :: TRCM_Y! 変数(v 格子上) xdim = SIZE( TRC, 1 ) - 1 ydim = SIZE( TRC, 2 ) - 1 dim = xdim - 1 ALLOCATE( TRCM_X(xdim,0:ydim), TRCM_Y(0:xdim,ydim) ) ALLOCATE( FLUX_X(0:xdim+1,ydim), FLUX_Y(xdim,0:ydim+1) ) ALLOCATE( ADVE_X(0:xdim,0:ydim), ADVE_Y(0:xdim,0:ydim) ) CALL get_mval_2D( TRC, TRCM_X, TRCM_Y ) FLUX_X(1:xdim,:) = u * TRCM_X(:,1:) FLUX_Y(:,1:ydim) = v * TRCM_Y(1:,:) FLUX_X(0,:) = FLUX_X(2,:) FLUX_X(xdim+1,:) = FLUX_X(dim,:) FLUX_Y(:,0) = FLUX_Y(:,2) FLUX_Y(:,ydim+1) = FLUX_Y(:,dim) ADVE_X(1:dim,1:dim) = & & ( - 1.125 * ( FLUX_X(2:xdim,1:dim) - FLUX_X(1:dim,1:dim) ) & & + 0.125/3.* ( FLUX_X(3:,1:dim) - FLUX_X(:dim-1,1:dim) ) )/dx ADVE_Y(1:dim,1:dim) = & & ( - 1.125 * ( FLUX_Y(1:dim,2:ydim) - FLUX_Y(1:dim,1:dim) ) & & + 0.125/3.* ( FLUX_Y(1:dim,3:) - FLUX_Y(1:dim,:dim-1) ) )/dy call bound2( ADVE_X ) call bound2( ADVE_Y ) DTRC = dt * ( ADVE_X + ADVE_Y ) !- 境界上での人工減衰 DTRC(1,:) = - 5.*dt*TRC(1,:) DTRC(dim,:) = - 5.*dt*TRC(dim,:) DTRC(:,1) = - 5.*dt*TRC(:,1) DTRC(:,dim) = - 5.*dt*TRC(:,dim) DEALLOCATE( TRCM_X, TRCM_Y, FLUX_X, FLUX_Y, ADVE_X, ADVE_Y ) END SUBROUTINE get_dtrc_C4_2D !========================================================================= SUBROUTINE get_dtrc_u1_2D( u, v, TRC, DTRC )! 上流差分 REAL,DIMENSION(:,:),INTENT(in) :: u REAL,DIMENSION(:,:),INTENT(in) :: v REAL,DIMENSION(0:,0:),INTENT(in):: TRC REAL,DIMENSION(0:,0:),INTENT(out):: DTRC INTEGER :: dim, xdim, ydim REAL,DIMENSION(:,:),ALLOCATABLE :: FLUX_X! フラックス1 REAL,DIMENSION(:,:),ALLOCATABLE :: FLUX_Y! フラックス2 REAL,DIMENSION(:,:),ALLOCATABLE :: ADVE_X! 移流項1 REAL,DIMENSION(:,:),ALLOCATABLE :: ADVE_Y! 移流項2 xdim = SIZE( TRC, 1 ) - 1 ydim = SIZE( TRC, 2 ) - 1 dim = xdim - 1 ALLOCATE( FLUX_X(xdim,ydim), FLUX_Y(xdim,ydim) ) ALLOCATE( ADVE_X(0:xdim,0:ydim), ADVE_Y(0:xdim,0:ydim) ) FLUX_X = MAX( u, 0. )*TRC(:dim,1:) + MIN( u, 0. )*TRC(1:,1:) FLUX_Y = MAX( v, 0. )*TRC(1:,:dim) + MIN( v, 0. )*TRC(1:,1:) ADVE_X(1:dim,1:dim) = - ( FLUX_X(2:,:dim) - FLUX_X(:dim,:dim) )/dx ADVE_Y(1:dim,1:dim) = - ( FLUX_Y(:dim,2:) - FLUX_Y(:dim,:dim) )/dy call bound2( ADVE_X ) call bound2( ADVE_Y ) DTRC = dt * ( ADVE_X + ADVE_Y ) DEALLOCATE( FLUX_X, FLUX_Y, ADVE_X, ADVE_Y ) END SUBROUTINE get_dtrc_u1_2D !========================================================================= SUBROUTINE get_dtrc_u_mpdata_2D( u, v, TRC, udif, vdif, DTRC )! MPDATA 上流差分 REAL,DIMENSION(:,:),INTENT(in) :: u REAL,DIMENSION(:,:),INTENT(in) :: v REAL,DIMENSION(0:,0:),INTENT(in):: TRC REAL,DIMENSION(:,:),INTENT(out) :: udif REAL,DIMENSION(:,:),INTENT(out) :: vdif REAL,DIMENSION(0:,0:),INTENT(out):: DTRC INTEGER :: dim, xdim, ydim REAL,DIMENSION(:,:),ALLOCATABLE :: TRCM_X! 変数(u 格子上) REAL,DIMENSION(:,:),ALLOCATABLE :: TRCM_Y! 変数(v 格子上) REAL,DIMENSION(:,:),ALLOCATABLE :: DTRC_X! 変数X方向差分(u 格子上) REAL,DIMENSION(:,:),ALLOCATABLE :: DTRC_Y! 変数Y方向差分(v 格子上) REAL,DIMENSION(:,:),ALLOCATABLE :: DTRC_Xvar! 変数X方向差分(v 格子上) REAL,DIMENSION(:,:),ALLOCATABLE :: DTRC_Yvar! 変数Y方向差分(u 格子上) REAL,DIMENSION(:,:),ALLOCATABLE :: u_var ! u(v格子上) REAL,DIMENSION(:,:),ALLOCATABLE :: v_var ! v(u格子上) REAL,PARAMETER :: epsilon = 1.e-8 xdim = SIZE( TRC, 1 ) - 1 ydim = SIZE( TRC, 2 ) - 1 dim = xdim - 1 ALLOCATE( TRCM_X(xdim,0:ydim), TRCM_Y(0:xdim,ydim) ) ALLOCATE( DTRC_X(xdim,ydim), DTRC_Y(xdim,ydim) ) ALLOCATE( DTRC_Xvar(xdim,ydim), DTRC_Yvar(xdim,ydim) ) ALLOCATE( u_var(xdim,ydim), v_var(xdim,ydim) ) CALL get_mval_2D( TRC, TRCM_X, TRCM_Y ) DTRC_X = TRC(1:,1:) - TRC(:dim,1:) DTRC_Y = TRC(1:,1:) - TRC(1:,:dim) CALL get_mval_utov( DTRC_X, DTRC_Xvar ) CALL get_mval_vtou( DTRC_Y, DTRC_Yvar ) CALL get_mval_utov( u, u_var ) CALL get_mval_vtou( v, v_var ) udif = 0.5*( ( ABS(u)*dx - dt*u**2 ) * DTRC_X / dx & & - dt*u*v_var*DTRC_Yvar / dy ) / ( TRCM_X(:,1:) + epsilon ) vdif = 0.5*( ( ABS(v)*dx - dt*v**2 ) * DTRC_Y / dy & & - dt*u_var*v*DTRC_Xvar / dx ) / ( TRCM_Y(1:,:) + epsilon ) CALL get_dtrc_u1_2D( udif, vdif, TRC, DTRC ) ! 上流差分 DEALLOCATE( TRCM_X, TRCM_Y, DTRC_X, DTRC_Y, & & DTRC_Xvar, DTRC_Yvar, u_var, v_var ) END SUBROUTINE get_dtrc_u_mpdata_2D !========================================================================= SUBROUTINE get_dtrc_fct_2D( u, v, TRC, DTRC )! FCT REAL,DIMENSION(:,:),INTENT(in) :: u REAL,DIMENSION(:,:),INTENT(in) :: v REAL,DIMENSION(0:,0:),INTENT(in):: TRC REAL,DIMENSION(0:,0:),INTENT(out):: DTRC INTEGER :: dim, xdim, ydim REAL,DIMENSION(:,:),ALLOCATABLE :: AFLUX_X! antidifuusiveフラックス1 REAL,DIMENSION(:,:),ALLOCATABLE :: AFLUX_Y! antidifuusiveフラックス2 REAL,DIMENSION(:,:),ALLOCATABLE :: ADVE_X! 移流項1 REAL,DIMENSION(:,:),ALLOCATABLE :: ADVE_Y! 移流項2 REAL,DIMENSION(:,:),ALLOCATABLE :: TRC_L! 低次スキームによる計算値 REAL,DIMENSION(:,:),ALLOCATABLE :: TRC_MAX REAL,DIMENSION(:,:),ALLOCATABLE :: TRC_MIN REAL,DIMENSION(:,:),ALLOCATABLE :: CF_X ! 補正係数1 REAL,DIMENSION(:,:),ALLOCATABLE :: CF_Y ! 補正係数2 xdim = SIZE( TRC, 1 ) - 1 ydim = SIZE( TRC, 2 ) - 1 dim = xdim - 1 ALLOCATE( ADVE_X(0:xdim,0:ydim), ADVE_Y(0:xdim,0:ydim) ) ALLOCATE( AFLUX_X(xdim,ydim), AFLUX_Y(xdim,ydim) ) ALLOCATE( TRC_L(0:xdim,0:xdim) ) ALLOCATE( TRC_MAX(0:xdim,0:xdim), TRC_MIN(0:xdim,0:xdim) ) ALLOCATE( CF_X(xdim,xdim), CF_Y(xdim,xdim) ) !- 低次スキームフラックスの計算: 上流差分 CALL get_dtrc_u1_2D( u, v, TRC, TRC_L ) TRC_L = TRC + TRC_L !- antidiffusive フラックスの計算 call get_aflux( TRC, AFLUX_X, AFLUX_Y ) !- リミッタの設定(1) call limit_aflux( TRC_L, AFLUX_X, AFLUX_Y ) !- リミット値の設定 call get_limitval( TRC, TRC_L, TRC_MAX, TRC_MIN ) !- 補正係数の計算 call get_cfact( TRC_L, TRC_MAX, TRC_MIN, AFLUX_X, AFLUX_Y, CF_X, CF_Y ) !- 補正フラックスの計算 AFLUX_X = CF_X * AFLUX_X AFLUX_Y = CF_Y * AFLUX_Y ADVE_X(1:dim,1:dim) = - ( AFLUX_X(2:,:dim) - AFLUX_X(:dim,:dim) )/dx ADVE_Y(1:dim,1:dim) = - ( AFLUX_Y(:dim,2:) - AFLUX_Y(:dim,:dim) )/dy call bound2( ADVE_X ) call bound2( ADVE_Y ) !- 増分計算 DTRC = TRC_L - TRC + dt * ( ADVE_X + ADVE_Y ) DEALLOCATE( ADVE_X, ADVE_Y, AFLUX_X, AFLUX_Y, TRC_L ) DEALLOCATE( TRC_MAX, TRC_MIN, CF_X, CF_Y ) CONTAINS !------------------------------------------------------------------------- SUBROUTINE get_aflux( TRC, AFLUX_X, AFLUX_Y ) REAL,DIMENSION(0:,0:),INTENT(in):: TRC REAL,DIMENSION(:,:),INTENT(out) :: AFLUX_X, AFLUX_Y REAL,DIMENSION(:,:),ALLOCATABLE :: FLUX_XH! フラックス1(高次スキーム) REAL,DIMENSION(:,:),ALLOCATABLE :: FLUX_YH! フラックス2(高次スキーム) REAL,DIMENSION(:,:),ALLOCATABLE :: FLUX_XL! フラックス1(低次スキーム) REAL,DIMENSION(:,:),ALLOCATABLE :: FLUX_YL! フラックス2(低次スキーム) REAL,DIMENSION(:,:),ALLOCATABLE :: TRCM_X! 変数(u 格子上) REAL,DIMENSION(:,:),ALLOCATABLE :: TRCM_Y! 変数(v 格子上) ALLOCATE( FLUX_XL(xdim,ydim), FLUX_YL(xdim,ydim) ) ALLOCATE( FLUX_XH(xdim,ydim), FLUX_YH(xdim,ydim) ) ALLOCATE( TRCM_X(xdim,0:ydim), TRCM_Y(0:xdim,ydim) ) !- 低次スキームフラックスの計算: 上流差分 FLUX_XL = MAX( u, 0. )*TRC(:dim,1:) + MIN( u, 0. )*TRC(1:,1:) FLUX_YL = MAX( v, 0. )*TRC(1:,:dim) + MIN( v, 0. )*TRC(1:,1:) !- 高次スキームフラックスの計算: 二次中央差分 CALL get_mval_2D( TRC, TRCM_X, TRCM_Y ) FLUX_XH = u * TRCM_X(:,1:) FLUX_YH = v * TRCM_Y(1:,:) !- antidiffusive フラックスの計算 AFLUX_X = FLUX_XH - FLUX_XL AFLUX_Y = FLUX_YH - FLUX_YL DEALLOCATE( FLUX_XL, FLUX_YL, FLUX_XH, FLUX_YH ) END SUBROUTINE get_aflux !------------------------------------------------------------------------- SUBROUTINE limit_aflux( TRC_L, AFLUX_X, AFLUX_Y ) REAL,DIMENSION(0:,0:),INTENT(in):: TRC_L REAL,DIMENSION(:,:),INTENT(inout) :: AFLUX_X, AFLUX_Y INTEGER :: i,j REAL,DIMENSION(:,:),ALLOCATABLE :: DIFF_X! 差分1 REAL,DIMENSION(:,:),ALLOCATABLE :: DIFF_Y! 差分2 ALLOCATE( DIFF_X(0:xdim+1,ydim), DIFF_Y(xdim,0:ydim+1) ) DIFF_X = 0. DIFF_Y = 0. DIFF_X(1:xdim,:) = TRC_L(1:,1:) - TRC_L(:dim,1:) DIFF_Y(:,1:ydim) = TRC_L(1:,1:) - TRC_L(1:,:dim) DO i = 1,xdim DO j = 1,ydim IF ( AFLUX_X(i,j)*DIFF_X(i,j) < 0. .AND. & & AFLUX_X(i,j)*DIFF_X(i-1,j) < 0. .OR. & & AFLUX_X(i,j)*DIFF_X(i+1,j) < 0. ) AFLUX_X(i,j) = 0. IF ( AFLUX_Y(i,j)*DIFF_Y(i,j) < 0. .AND. & & AFLUX_Y(i,j)*DIFF_Y(i,j-1) < 0. .OR. & & AFLUX_Y(i,j)*DIFF_Y(i,j+1) < 0. ) AFLUX_Y(i,j) = 0. END DO END DO DEALLOCATE( DIFF_X, DIFF_Y ) END SUBROUTINE limit_aflux !------------------------------------------------------------------------- SUBROUTINE get_limitval( TRC, TRC_L, TRC_MAX, TRC_MIN ) REAL,DIMENSION(0:,0:),INTENT(in):: TRC, TRC_L REAL,DIMENSION(0:,0:),INTENT(out) :: TRC_MAX, TRC_MIN REAL,DIMENSION(:,:),ALLOCATABLE :: TRC_A REAL,DIMENSION(:,:),ALLOCATABLE :: TRC_B ALLOCATE( TRC_A(0:xdim,0:xdim), TRC_B(0:xdim,0:xdim) ) TRC_A = MAX( TRC, TRC_L ) TRC_B = MIN( TRC, TRC_L ) TRC_MAX(1:dim,1:dim) = MAX( TRC_A(1:dim, 1:dim), & & TRC_A(0:dim-1,1:dim), & & TRC_A(2:xdim, 1:dim), & & TRC_A(1:dim, 0:dim-1), & & TRC_A(1:dim, 2:ydim ) ) TRC_MIN(1:dim,1:dim) = MIN( TRC_B(1:dim, 1:dim), & & TRC_B(0:dim-1,1:dim), & & TRC_B(2:xdim, 1:dim), & & TRC_B(1:dim, 0:dim-1), & & TRC_B(1:dim, 2:ydim ) ) call bound2( TRC_MAX ) call bound2( TRC_MIN ) DEALLOCATE( TRC_A, TRC_B ) END SUBROUTINE get_limitval !------------------------------------------------------------------------- SUBROUTINE get_cfact( TRC_L, TRC_MAX, TRC_MIN, AFLUX_X, AFLUX_Y, CF_X, CF_Y ) REAL,DIMENSION(0:,0:),INTENT(in):: TRC_L REAL,DIMENSION(0:,0:),INTENT(in) :: TRC_MAX, TRC_MIN REAL,DIMENSION(:,:),INTENT(in) :: AFLUX_X, AFLUX_Y REAL,DIMENSION(:,:),INTENT(out):: CF_X, CF_Y REAL,DIMENSION(:,:),ALLOCATABLE :: SN REAL,DIMENSION(:,:),ALLOCATABLE :: PP, QP, RP REAL,DIMENSION(:,:),ALLOCATABLE :: PN, QN, RN ALLOCATE( SN(xdim,ydim) ) ALLOCATE( PP(0:xdim,0:ydim), QP(0:xdim,0:ydim), RP(0:xdim,0:ydim) ) ALLOCATE( PN(0:xdim,0:ydim), QN(0:xdim,0:ydim), RN(0:xdim,0:ydim) ) PP(1:dim,1:dim) = MAX(0., AFLUX_X(:dim,:))*dy & & - MIN(0., AFLUX_X(2:,:))*dy & & + MAX(0., AFLUX_Y(:,:dim))*dx & & - MIN(0., AFLUX_Y(:,2:))*dx PP(0,:) = - MIN( 0., AFLUX_X(1,:) )*dy PP(xdim,:) = MAX( 0., AFLUX_X(xdim,:) )*dy PP(:,0) = - MIN( 0., AFLUX_Y(:,1) )*dx PP(:,ydim) = MAX( 0., AFLUX_X(:,ydim) )*dx QP = ( TRC_MAX - TRC_L )*dx*dy RP = 0. RP = MIN( 1., QP/(PP + 1.e-8)) PN(1:dim,1:dim) = MAX(0., AFLUX_X(2:,:))*dy & & - MIN(0., AFLUX_X(:dim,:))*dy & & + MAX(0., AFLUX_Y(:,2:))*dx & & - MIN(0., AFLUX_Y(:,:dim))*dx PN(0,:) = MAX( 0., AFLUX_X(1,:) )*dy PN(xdim,:) = - MIN( 0., AFLUX_X(xdim,:) )*dy PN(:,0) = MAX( 0., AFLUX_Y(:,1) )*dx PN(:,ydim) = - MIN( 0., AFLUX_X(:,ydim) )*dx QN = ( TRC_L - TRC_MIN )*dx*dy RN = 0. RN = MIN( 1., QP/(PN + 1.e-8)) SN = INT( SIGN( 1., AFLUX_X ) ) CF_X = MIN( MAX( 0.5*( 1 + SN )*RP(1:,1:), 0.5*( 1 - SN )*RN(1:,1:) ) & & , MAX( 0.5*( 1 + SN )*RN(:dim,1:), 0.5*( 1 - SN )*RP(:dim,1:) ) ) SN = INT( SIGN( 1., AFLUX_Y ) ) CF_Y = MIN( MAX( 0.5*( 1 + SN )*RP(1:,1:), 0.5*( 1 - SN )*RN(1:,1:) ) & & , MAX( 0.5*( 1 + SN )*RN(1:,:dim), 0.5*( 1 - SN )*RP(1:,:dim) ) ) DEALLOCATE( SN, PP, QP, RP, PN, QN, RN ) END SUBROUTINE get_cfact !------------------------------------------------------------------------- END SUBROUTINE get_dtrc_fct_2D !========================================================================= SUBROUTINE get_mval_2D( val, mval_x, mval_y ) REAL,DIMENSION(0:,0:),INTENT(in):: val REAL,DIMENSION(:,0:),INTENT(out):: mval_x REAL,DIMENSION(0:,:),INTENT(out):: mval_y INTEGER :: xdim,ydim,dim xdim = SIZE( mval_x, 1 ) ydim = SIZE( mval_y, 2 ) dim = xdim - 1 !- 2 点平均 ! mval_x = ( val(1:,:) + val(:xdim-1,:) )/2 ! mval_y = ( val(:,1:) + val(:,:ydim-1) )/2 !- 6 点平均 mval_x(:,1:dim) = ( val(1:,0:dim-1) + val(:dim,0:dim-1) + & & 2.0*val(1:,1:dim) + 2.0*val(:dim,1:dim) + & & val(1:,2:ydim) + val(:dim,2:ydim) )/8. mval_y(1:dim,:) = ( val(0:dim-1,1:) + val(0:dim-1,:dim) + & & 2.0*val(1:dim,1:) + 2.0*val(1:dim,:dim) + & & val(2:xdim,1:) + val(2:xdim,:dim) )/8. END SUBROUTINE get_mval_2D !========================================================================= SUBROUTINE get_mval_utov( val_u, val_v ) REAL,DIMENSION(:,:),INTENT(in) :: val_u ! u 格子上変数値 REAL,DIMENSION(:,:),INTENT(out) :: val_v ! v 格子上変数値 INTEGER :: xdim,ydim,dim xdim = SIZE( val_u, 1 ) ydim = SIZE( val_u, 2 ) dim = xdim - 1 val_v(:dim,2:) = ( val_u(:dim,2:) + val_u(2:xdim,2:) + & & val_u(:dim,1:dim)+ val_u(2:xdim,1:dim) )/4. val_v(:dim,1) = ( val_u(:dim,1) + val_u(2:xdim,1) )/2. val_v(xdim,2:) = ( val_u(xdim,:dim) + val_u(xdim,2:) )/2. val_v(xdim,1) = val_u(xdim,1) END SUBROUTINE get_mval_utov !========================================================================= SUBROUTINE get_mval_vtou( val_v, val_u ) REAL,DIMENSION(:,:),INTENT(in) :: val_v ! v 格子上変数値 REAL,DIMENSION(:,:),INTENT(out) :: val_u ! u 格子上変数値 INTEGER :: xdim,ydim,dim xdim = SIZE( val_v, 1 ) ydim = SIZE( val_v, 2 ) dim = xdim - 1 val_u(2:,:dim) = ( val_v(2:,:dim) + val_v(2:,2:ydim) + & & val_v(1:dim,:dim)+ val_v(1:dim,2:ydim) )/4. val_u(1,:dim) = ( val_v(1,:dim) + val_v(1,2:ydim) )/2. val_u(2:,ydim) = ( val_v(:dim,ydim) + val_v(1:,ydim) )/2. val_u(1,ydim) = val_v(1,ydim) END SUBROUTINE get_mval_vtou !========================================================================= SUBROUTINE get_plotval( val, pval ) REAL,DIMENSION(0:,0:),INTENT(in):: val REAL,DIMENSION(:,:),INTENT(out) :: pval INTEGER :: xdim,ydim xdim = SIZE( pval, 1 ) ydim = SIZE( pval, 2 ) !- 4 点平均 pval = ( val(1:,1:) + val(1:,:ydim-1) + & & val(:xdim-1,1:) + val(:xdim-1,:ydim-1) )/4. END SUBROUTINE get_plotval !========================================================================= SUBROUTINE bound2( val ) REAL,DIMENSION(0:,0:),INTENT(inout):: val INTEGER :: dim,xdim,ydim xdim = SIZE( val, 1 ) - 1 ydim = SIZE( val, 2 ) - 1 dim = xdim - 1 !- 周期境界条件 ! val(0,:) = val(dim,:) ! val(xdim,:) = val(1,:) ! val(:,0) = val(:,dim) ! val(:,ydim) = val(:,1) ! val(0,0) = val(dim,dim) ! val(0,ydim) = val(dim,1) ! val(xdim,0) = val(1,dim) ! val(xdim,ydim) = val(1,1) !- 壁で対称(1階微分が0) val(0,:) = val(1,:) val(xdim,:) = val(dim,:) val(:,0) = val(:,1) val(:,ydim) = val(:,dim) val(0,0) = val(1,1) val(0,ydim) = val(1,dim) val(xdim,0) = val(dim,1) val(xdim,ydim) = val(dim,dim) END SUBROUTINE bound2 !========================================================================= END MODULE cal_module