!= dcpam 主プログラム
!
!= dcpam main program
!
! Authors::   Yasuhiro Morikawa, Satoshi Noda, Yoshiyuki O. Takahashi, Shin-ichi Takehiro
! Version::   $Id: dcpam_main.f90,v 1.67 2026/05/21 22:00:00 takepiro Exp $ 
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008-2025. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

program dcpam_main_spml_test
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! モデルの使い方については {チュートリアル}[link:../../../doc/tutorial/rakuraku/] を
  ! 参照してください.
  !
  ! See {Tutorial}[link:../../../doc/tutorial/rakuraku/index.htm.en] for usage of the 
  ! model. 
  !

  use dc_message, only: MessageNotify

  ! メッセージ制御
  ! Message control
  !
  use mpi_messagecntl, only : DoesOutputMPIMessage, MPIMessageCntlInit

  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! 倍精度実数型. Double precision. 
    &                 STDOUT, &  ! 標準出力の装置番号. Unit number of standard output
    &                 STRING, &  ! 文字列.       Strings. 
    &                 TOKEN      ! キーワード.   Keywords. 

  ! ファイル入出力補助
  ! File I/O support
  !
  use dc_iounit, only: FileOpen

  use dc_test, only : AssertEqual

  ! モジュール引用 ; USE statements
  !
  ! MPI
  !
  use mpi_wrapper, only : MPIWrapperInit, MPIWrapperFinalize, myrank

  ! コマンドライン引数処理
  ! Command line option parser
  !
  use option_parser, only: OptParseInit

  ! NAMELIST ファイル入力に関するユーティリティ
  ! Utilities for NAMELIST file input
  !
  use namelist_util, only: NmlutilInit, NmlutilMsg

  ! 時刻管理
  ! Time control
  !
  use timeset, only: TimesetInit, TimesetProgress, &
    & TimeB, &                ! ステップ $ t - \Delta t $ の時刻.
                              ! Time of step $ t - \Delta t $.
    & TimeN, &                ! ステップ $ t $ の時刻.
                              ! Time of step $ t $.
    & TimeA, &                ! ステップ $ t + \Delta t $ の時刻.
                              ! Time of step $ t + \Delta t $.
    & EndTime, &              ! 計算終了時刻.
                              ! End time of calculation
    & DelTime                 ! $ \Delta t $ [s]

  ! 物理定数設定
  ! Physical constants settings
  !
  use constants0, only: PI
  use constants, only: ConstantsInit, Rplanet

  ! 格子点設定
  ! Grid points settings
  !
  use gridset, only: GridSetInit, &
                     imax, &  ! 経度格子点数. 
                              ! Number of grid points in longitude
    &                jmax, &  ! 緯度格子点数. 
                              ! Number of grid points in latitude
    &                kmax     ! 鉛直層数. 
                              ! Number of vertical level

  use axesset    , only : AxesSetInit, x_Lon, y_Lat, z_Sigma, r_Sigma

  use composition, only:  CompositionInit, ncmax
  !
  !

  ! 出力ファイルの基本情報管理
  ! Management basic information for output files
  !
  use fileset, only: FilesetInit
  
  ! ヒストリデータファイルの初期化
  ! Initialization of history data files
  !
  use history_file_io, only: HistoryFileOpen, HistoryFileClose

  
  use sltt       , only: SLTTinit, SLTTHorAdv, SLTTVerAdv, SLTTCheckInit

  use gtool_historyauto, only: &
       & HistoryAutoAddVariable, HistoryAutoPut

  use gtool_history, only: HistoryClose

  ! 宣言文 ; Declaration statements
  !
  implicit none

  character(*), parameter:: prog_name = 'dcpam_main'
                            ! 主プログラム名. 
                            ! Main program name

  character(STRING)      :: namelist_filename
                            ! NAMELIST ファイルの名称. 
                            ! NAMELIST file name

  ! 予報変数 (ステップ $ t-\Delta t $ , $ t $ , $ t+\Delta t $ )
  ! Prediction variables  (Step $ t-\Delta t $ , $ t $ , $ t+\Delta t $ )
  !
  real(DP), allocatable:: xyz_U (:,:,:)
                              ! $ u (t-\Delta t) $ .   東西風速. Eastward wind (m s-1)
  real(DP), allocatable:: xyz_V (:,:,:)
                              ! $ v (t-\Delta t) $ .   南北風速. Northward wind (m s-1)
  real(DP), allocatable:: xyr_SigDot (:,:,:)
                              ! $ u (t-\Delta t) $ .   鉛直風速. vertial wind (m s-1)
  real(DP), allocatable:: xyzf_QMixA(:,:,:,:)
                              ! $ q (t+\Delta t) $ .   混合比. Mass mixing ratio of constituents (1)
  real(DP), allocatable:: xyzf_QMixLinA(:,:,:,:)
                              ! $ q (t+\Delta t) $ .   混合比. Mass mixing ratio
  real(DP), allocatable:: xyzf_QMixN(:,:,:,:)
                              ! $ q (t) $ .            混合比. Mass mixing ratio of constituents (1)
   real(DP), allocatable:: xyzf_QMixB(:,:,:,:)
                              ! $ q (t-\Delta t) $ .   混合比. Mass mixing ratio of constituents (1)
!!$  real(DP), allocatable:: xyzf_QMixS(:,:,:,:)
!!$                              ! $ q (t-\Delta t) $ .   混合比. Mass mixing ratio of constituents (1)

  real(DP) :: V0               ! 子午面流の速度振幅
  
  integer :: i, j, k, n
  
  integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号. 
                            ! Unit number for NAMELIST file open
  integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT. 
                            ! IOSTAT of NAMELIST read

  ! NAMELIST 変数群
  ! NAMELIST group name
  !
  namelist /sltt_meridadv_test_nml/  V0

  ! 判定誤差設定
  integer, parameter :: check_digits = 3
  integer, parameter :: ignore = -4

  ! SLTT ルーチンチェック出力スイッチ
  ! logical :: sltt_check = .false.
  logical :: sltt_check = .true.
  
  !======================================================================
  ! 実行文 ; Executable statement
  !
  ! Initialize MPI
  !
  call MPIWrapperInit

  ! コマンドライン引数処理
  ! Command line option parser
  !
  call OptParseInit(       &
       & namelist_filename,    & ! (out)
       & prog_name            & ! (in )
       & )

  ! NAMELIST ファイル名入力
  ! Input NAMELIST file name
  !
  call NmlutilInit( &
       & namelist_filename  & ! (in)
       & )

  ! 時刻管理
  ! Time control
  !
  call TimesetInit

  ! 出力ファイルの基本情報管理
  ! Management basic information for output files
  !
  call FilesetInit

  ! 格子点設定
  ! Grid points settings
  !
  call GridsetInit
  
  ! 組成に関わる配列の設定
  ! Settings of array for atmospheric composition
  !
  call CompositionInit

  ! 物理定数設定
  ! Physical constants settings
  !
  call ConstantsInit

  ! 座標データ設定
  ! Axes data settings
  !
  call AxessetInit

  ! ヒストリデータファイルの初期化
  ! Initialization of history data files
  !
  call HistoryFileOpen
  
  call SLTTInit
  if ( sltt_check ) call SLTTCheckInit
 
  ! NAMELIST の読み込み
  ! NAMELIST is input
  !
  if ( trim(namelist_filename) /= '' ) then
     call FileOpen( unit_nml, &          ! (out)
          & namelist_filename, mode = 'r' ) ! (in)

     rewind( unit_nml )
     read( unit_nml, &                ! (in)
          & nml = sltt_meridadv_test_nml, &  ! (out)
          & iostat = iostat_nml )        ! (out)
     close( unit_nml )

     call NmlutilMsg( iostat_nml, 'main' ) ! (in)
     if ( iostat_nml == 0 .AND. DoesOutputMPIMessage() ) write( STDOUT, nml = sltt_meridadv_test_nml )
  end if

  call HistoryAutoAddVariable( 'U' ,            &  ! (in)
        & (/ 'lon ', 'lat ', 'sig ','time'/),   &  ! (in)
        & 'U', 'm/s' )                             ! (in)
  call HistoryAutoAddVariable( 'V' ,            &  ! (in)
        & (/ 'lon ', 'lat ', 'sig ','time'/),   &  ! (in)
        & 'V', 'm/s' )                             ! (in)
  call HistoryAutoAddVariable( 'SigDot' ,       &  ! (in)
        & (/ 'lon ', 'lat ', 'sigm','time'/),   &  ! (in)
        & 'SigDot', '1/s' )                        ! (in)
!!$  call HistoryAutoAddVariable( 'QMixS' ,        &  ! (in)
!!$        & (/ 'lon ', 'lat ', 'sig ', 'time'/),  &  ! (in)
!!$        & 'QMixS', 'kg kg-1' )                     ! (in)
  call HistoryAutoAddVariable( 'QVap',          &  ! (in)
        & (/ 'lon ', 'lat ', 'sig ','time'/),   &  ! (in)
        & 'QVap', 'kg kg-1' )                      ! (in)
  call HistoryAutoAddVariable( 'QLiq',          &  ! (in)
        & (/ 'lon ', 'lat ', 'sig ','time'/),   &  ! (in)
        & 'QLiq', 'kg kg-1' )                      ! (in)
  call HistoryAutoAddVariable( 'QSol',          &  ! (in)
        & (/ 'lon ', 'lat ', 'sig ','time'/),   &  ! (in)
        & 'QSol', 'kg kg-1' )                      ! (in)

  
  allocate( xyz_U (0:imax-1, 1:jmax, 1:kmax) )
  allocate( xyz_V (0:imax-1, 1:jmax, 1:kmax) )
  allocate( xyr_SigDot (0:imax-1, 1:jmax, 0:kmax) )
  allocate( xyzf_QMixB(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )
  allocate( xyzf_QMixA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )
  allocate( xyzf_QMixLinA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )
  allocate( xyzf_QMixN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )
!!$  allocate( xyzf_QMixS(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )

  do k = 1, kmax
     do j = 1, jmax
        do i = 0, imax-1
           xyz_U(i,j,k) = 0.0D0
!           xyz_V(i,j,k) = -V0*PI*cos(PI*z_Sigma(k))*cos(y_Lat(j))
           xyz_V(i,j,k) = -V0*PI*cos(PI*z_Sigma(k))*3*cos(y_Lat(j))*sin(y_Lat(j))
        end do
     end do
  end do
  do k = 0, kmax
     do j = 1, jmax
        do i = 0, imax-1
!           xyr_SigDot(i,j,k) = - 2.0*V0/RPlanet * sin(PI*r_Sigma(k))*sin(y_Lat(j))
           xyr_SigDot(i,j,k) = - 3.0*V0/RPlanet * sin(PI*r_Sigma(k)) &
                * ( 3*sin(y_Lat(j))**2 -1.0)
        end do
     end do
  end do
  
  do k = 1, kmax
     do j = 1, jmax
        do i = 0, imax-1
           xyzf_QMixB(i,j,k,:) = f_distribution(y_Lat(j),z_Sigma(k))
           xyzf_QMixN(i,j,k,:) = f_distribution(y_Lat(j),z_Sigma(k))
        end do
     end do
  end do

  call HistoryAutoPut( TimeN, 'U', xyz_U )
  call HistoryAutoPut( TimeN, 'V', xyz_V )
  call HistoryAutoPut( TimeN, 'SigDot', xyr_SigDot )
  call HistoryAutoPut( TimeN, 'QVap',  xyzf_QMixN(:,:,:,1) )
  call HistoryAutoPut( TimeN, 'QLiq', xyzf_QMixN(:,:,:,2) )
  call HistoryAutoPut( TimeN, 'QSol', xyzf_QMixN(:,:,:,3) )
  
  ! 時間積分
  ! Time integration
  !
  loop_time : do while ( TimeB < EndTime )

     xyzf_QMixA = QmixAdv( xyz_U, xyz_V, xyr_SigDot, xyzf_QMixN, xyzf_QMixB )

     call HistoryAutoPut( TimeA, 'U', xyz_U )
     call HistoryAutoPut( TimeA, 'V', xyz_V )
     call HistoryAutoPut( TimeA, 'SigDot', xyr_SigDot )
     call HistoryAutoPut( TimeN, 'QVap',  xyzf_QMixN(:,:,:,1) )
     call HistoryAutoPut( TimeN, 'QLiq', xyzf_QMixN(:,:,:,2) )
     call HistoryAutoPut( TimeN, 'QSol', xyzf_QMixN(:,:,:,3) )

     ! 次の時間ステップに向けて予報変数の入れ替え
     ! Exchange prediction variables for the next time step
     !
     xyzf_QMixB = xyzf_QMixN ; xyzf_QMixN = xyzf_QMixA

     ! 時刻の進行
     ! Progress time
     !
     call TimesetProgress

  end do loop_time

  call HistoryFileClose
  if ( sltt_check ) call HistoryClose
  
  !----------------------------------------------------------------------
  ! Finalize MPI
  !
  call MPIWrapperFinalize

contains
  !
  ! 移流計算 (スペクトル)
  !
  function QmixAdv( xyz_UN, xyz_VN, xyr_SigDotN, xyzf_QMixN, xyzf_QMixB ) result(xyzf_QMixA)

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: lmax, &  ! スペクトルデータの配列サイズ
                                ! Size of array for spectral data
         &             imax, &  ! 経度格子点数. 
                                ! Number of grid points in longitude
         &             jmax, &  ! 緯度格子点数. 
                                ! Number of grid points in latitude
         &             kmax     ! 鉛直層数. 
                                ! Number of vertical level

    use axesset    , only : y_Lat

    ! 組成に関わる配列の設定
    ! Settings of array for atmospheric composition
    !
    use composition, only:  ncmax  ! 成分の数
                                   ! Number of composition

#ifdef AXISYMMETRY
    use wa_zonal_module, only:   &
      & wa_Lapla_wa,       &
      & w_Lapla_w,         &
      & wa_LaplaInv_wa,    &
      & wa_DivLambda_xya,  &
      & wa_DivMu_xya,      &
      & xy_GradLambda_w,   &
      & xya_GradLambda_wa, &
      & xy_GradMu_w,       &
      & xya_GradMu_wa,     &
      & w_xy, xy_w,        &
      & wa_xya, xya_wa
#elif LIB_MPI
    use ua_mpi_module_base, only:                &
      & w_xy              => u_pv,               &
      & xy_w              => pv_u,               &
      & wa_xya            => ua_pva,             &
      & xya_wa            => pva_ua
    use ua_mpi_module_deriv, only:               &
      & wa_Lapla_wa       => ua_Lapla_ua,        &
      & w_Lapla_w         => u_Lapla_u,          &
      & wa_LaplaInv_wa    => ua_LaplaInv_ua,     &
      & wa_DivLambda_xya  => ua_DivLambda_pva,   &
      & wa_DivMu_xya      => ua_DivMu_pva,       &
      & xy_GradLambda_w   => pv_GradLambda_u,    &
      & xya_GradLambda_wa => pva_GradLambda_ua,  &
      & xy_GradMu_w       => pv_GradMu_u,        &
      & xya_GradMu_wa     => pva_GradMu_ua
#else
    use wa_module, only:   &
      & wa_Lapla_wa,       &
      & w_Lapla_w,         &
      & wa_LaplaInv_wa,    &
      & wa_DivLambda_xya,  &
      & wa_DivMu_xya,      &
      & xy_GradLambda_w,   &
      & xya_GradLambda_wa, &
      & xy_GradMu_w,       &
      & xya_GradMu_wa,     &
      & w_xy, xy_w,        &
      & wa_xya, xya_wa
#endif

    real(DP), intent(IN) :: xyz_UN   (0:imax-1, 1:jmax, 1:kmax)
                          ! $ U (t-\Delta t) = u (t-\Delta t) $ .
    real(DP), intent(IN) :: xyz_VN   (0:imax-1, 1:jmax, 1:kmax)
                          ! $ V (t-\Delta t) = v (t-\Delta t) $
    real(DP), intent(IN) :: xyr_SigDotN (0:imax-1, 1:jmax, 0:kmax)
                          ! $ \dot{\sigma} $ 鉛直流. Vertical flow
    real(DP), intent(IN):: xyzf_QMixN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                          ! $ q (t) $ .   比湿. Specific humidity
    real(DP), intent(IN):: xyzf_QMixB(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                          ! $ q (t-\Delta t) $ .   比湿. Specific humidity
    real(DP) ::  xyzf_QMixA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                          ! $ q (t+\Delta t) $ .   比湿. Specific humidity

    ! 作業変数
    ! Work variables
    !
    real(DP):: xyz_UCosLatN   (0:imax-1, 1:jmax, 1:kmax)
                          ! $ U (t) = u (t) \cos \varphi $ .
    real(DP):: xyz_VCosLatN   (0:imax-1, 1:jmax, 1:kmax)
                          ! $ V (t) = v (t) \cos \varphi $ .
    real(DP):: xyz_DivN       (0:imax-1, 1:jmax, 1:kmax)
                          ! $ D (t) $ .     発散. Divergence
    real(DP):: wzf_QMixA      (lmax, 1:kmax, 1:ncmax)
                          ! $ q (t+\Delta t) $ . 比湿 (スペクトル). 
                          ! Specific humidity (spectral)
    real(DP):: wzf_QMixB      (lmax, 1:kmax, 1:ncmax)
                          ! $ q (t-\Delta t) $ . 比湿 (スペクトル). 
                          ! Specific humidity (spectral)
    real(DP):: xyzf_QMixUAdvN (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                          ! $ Uq (t) $ . 比湿東西移流項. 
                          ! Eastward advection of specific humidity
    real(DP):: xyzf_QMixVAdvN (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                          ! $ Vq (t) $ . 比湿南北移流項. 
                          ! Northward advection of specific humidity
    real(DP):: wzf_DQMixDtN(lmax, 1:kmax, 1:ncmax)
                          ! $ \DD{q}{t} (t) $ . 比湿変化 (スペクトル). 
                          ! Specific humidity tendency (spectral)
    real(DP):: xyzf_QMixNonLinearN (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                          ! $ R (t) $ . 比湿時間変化項. 
                          ! Specific humidity tendency

    integer :: j, n

    do j=1,jmax
       xyz_UCosLatN(:,j,:) = xyz_UN(:,j,:)*cos(y_Lat(j))
       xyz_VCosLatN(:,j,:) = xyz_VN(:,j,:)*cos(y_Lat(j))
    enddo

    xyz_DivN  =  xya_wa (   wa_DivLambda_xya( xyz_UCosLatN )  &
                          + wa_DivMu_xya    ( xyz_VCosLatN ) ) / RPlanet
    
    do n = 1, ncmax
       wzf_QMixB(:,:,n) = wa_xya( xyzf_QMixB(:,:,:,n) )
    end do
    !
    call NonLinearOnGridQMix( &
         & xyz_UCosLatN, xyz_VCosLatN, xyz_DivN, xyr_SigDotN, & ! (in)
         & xyzf_QMixN,                                        & ! (in)
                                !
         & xyzf_QMixNonLinearN, &             ! (out)
         & xyzf_QMixUAdvN, xyzf_QMixVAdvN &   ! (out)
         & )

    do n = 1, ncmax
       wzf_DQMixDtN(:,:,n) &
            = - (   wa_DivLambda_xya( xyzf_QMixUAdvN(:,:,:,n) )   &
            + wa_DivMu_xya    ( xyzf_QMixVAdvN(:,:,:,n) ) ) &
            / RPlanet &
            + wa_xya(  xyzf_QMixNonLinearN(:,:,:,n) )
    end do

    do n = 1, ncmax
       do k = 1, kmax
          wzf_QMixA(:,k,n) = &
               wzf_QMixB(:,k,n) + 2.0_DP * DelTime * wzf_DQMixDtN(:,k,n) 
       end do
    end do

    do n = 1, ncmax
       xyzf_QMixA(:,:,:,n) = xya_wa( wzf_QMixA(:,:,n) )
    end do
    
  end function QmixAdv
  

  !--------------------------------------------------------------------------------------

  subroutine NonLinearOnGridQMix( &
    & xyz_UCosLatN, xyz_VCosLatN, xyz_DivN, xyr_SigDotN, & ! (in)
    & xyzf_QMixN,                                        & ! (in)
!
    & xyzf_QMixNonLinearN, &             ! (out)
    & xyzf_QMixUAdvN, xyzf_QMixVAdvN &   ! (out)
    & )
    !
    ! 非線形項 (非重力波項) を格子点上で計算します.
    !
    ! Non-linear terms (non gravitational terms) are calculated on
    ! grid points
    !

    ! モジュール引用 ; USE statements
    !

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: &
      & z_DelSigma            ! $ \Delta \sigma $ (整数). 
                              ! $ \Delta \sigma $ (Full)

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, & ! 経度格子点数. 
                               ! Number of grid points in longitude
      &                jmax, & ! 緯度格子点数. 
                               ! Number of grid points in latitude
      &                kmax    ! 鉛直層数. 
                               ! Number of vertical level

    ! 組成に関わる配列の設定
    ! Settings of array for atmospheric composition
    !
    use composition, only: &
      &                    ncmax
                               ! 成分の数
                               ! Number of composition

    implicit none
    real(DP), intent(in):: xyz_UCosLatN   (0:imax-1, 1:jmax, 1:kmax)
                              ! $ U (t) = u (t) \cos \varphi $ .
    real(DP), intent(in):: xyz_VCosLatN   (0:imax-1, 1:jmax, 1:kmax)
                              ! $ V (t) = v (t) \cos \varphi $ .
    real(DP), intent(in):: xyz_DivN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ D (t) $ . 発散. Divergence
    real(DP), intent(in ):: xyr_SigDotN (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \dot{\sigma} (t) $ .
                              ! 鉛直流. Vertical flow
    real(DP), intent(in):: xyzf_QMixN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ q (t) $ . 比湿. Specific humidity
    real(DP), intent(out):: xyzf_QMixNonLinearN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ R (t) $ . 比湿時間変化項. 
                              ! Specific humidity tendency
    real(DP), intent(out):: xyzf_QMixUAdvN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ Uq (t) $ . 比湿東西移流項. 
                              ! Eastward advection of specific humidity
    real(DP), intent(out):: xyzf_QMixVAdvN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ Vq (t) $ . 比湿南北移流項. 
                              ! Northward advection of specific humidity

    !-----------------------------------
    !  作業変数
    !  Work variables
    !
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
    integer:: n               ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituents

    ! 実行文 ; Executable statement
    !


    ! 比湿東西移流項, 比湿南北移流項の計算
    ! Calculate eastward and northward advection of specific humidity
    !
    do n = 1, ncmax
      do k = 1, kmax
        xyzf_QMixUAdvN(:,:,k,n) =  xyz_UCosLatN(:,:,k) * xyzf_QMixN(:,:,k,n)
        xyzf_QMixVAdvN(:,:,k,n) =  xyz_VCosLatN(:,:,k) * xyzf_QMixN(:,:,k,n)
      end do
    end do

    ! 比湿時間変化項 $ R $ の計算
    ! Calculate specific humidity tendency $ R $
    !
    do n = 1, ncmax
      xyzf_QMixNonLinearN(:,:,1,n) =   &
        &    xyzf_QMixN(:,:,1,n) * xyz_DivN(:,:,1) &
        &  - 1.0_DP / ( 2.0_DP * z_DelSigma(1) ) &
        &       * xyr_SigDotN(:,:,1) &
        &       * ( xyzf_QMixN(:,:,1,n) - xyzf_QMixN(:,:,2,n) )

      do k = 2, kmax - 1
        xyzf_QMixNonLinearN(:,:,k,n) =   &
          &    xyzf_QMixN(:,:,k,n) * xyz_DivN(:,:,k) &
          &  - 1.0_DP / ( 2.0_DP * z_DelSigma(k) ) &
          &       * ( &
          &             xyr_SigDotN(:,:,k-1) &
          &               * ( xyzf_QMixN(:,:,k-1,n) - xyzf_QMixN(:,:,k  ,n)   ) &
          &           + xyr_SigDotN(:,:,k) &
          &               * ( xyzf_QMixN(:,:,k  ,n) - xyzf_QMixN(:,:,k+1,n) ) )
      end do

      xyzf_QMixNonLinearN(:,:,kmax,n) =   &
        &    xyzf_QMixN(:,:,kmax,n) * xyz_DivN(:,:,kmax) &
        &  - 1.0_DP / ( 2.0_DP * z_DelSigma(kmax) ) &
        &       * xyr_SigDotN(:,:,kmax-1) &
        &       * ( xyzf_QMixN(:,:,kmax-1,n) - xyzf_QMixN(:,:,kmax,n) )
    end do


  end subroutine NonLinearOnGridQMix


  !----------------------------------------------------------------------
  !
  ! トレーサー分布設定
  !
  real(dp) function f_distribution(phi, sig)
    real(dp), intent(in) :: phi, sig 
    character(2), save :: disttype='CB'  ! 分布タイプ
                                         ! CB : cosine Bell
                                         ! 11,21 : (l,n)
    real(dp), save  :: radius            ! cosine Bell 型の半径

    namelist /fdist_nml/ disttype, radius

    logical, save :: first = .true.
    
    
    if ( first ) then
       first = .false.
       radius = 1.0D0/3.0D0         ! デフォルト値
       
       ! NAMELIST の読み込み
       ! NAMELIST is input
       !
       if ( trim(namelist_filename) /= '' ) then
          call FileOpen( unit_nml, &          ! (out)
               & namelist_filename, mode = 'r' ) ! (in)

          rewind( unit_nml )
          read( unit_nml, &                ! (in)
               & nml = fdist_nml, &  ! (out)
               & iostat = iostat_nml )        ! (out)
          close( unit_nml )

          call NmlutilMsg( iostat_nml, 'main' ) ! (in)
          if ( iostat_nml == 0 .AND. DoesOutputMPIMessage() ) write( STDOUT, nml = fdist_nml )
       end if
    endif

    if ( disttype == 'CB' ) then
       f_distribution = CosBell(phi, sig, radius)
    else
       f_distribution = Distln(phi, sig, disttype)
    end if

  end function f_distribution

  ! --- 元の分布関数 (cos Bell) ---
  real(dp) function CosBell(phi, sig, radius)
    real(dp), intent(in) :: phi, sig, radius
    real(dp) :: dist, center_phi, center_sig

    center_sig = 0.5_dp
    center_phi = 0.0_dp

    ! 角度距離の計算 (単純化のためユークリッド距離、厳密には大円距離も可)
    dist = sqrt((sig - center_sig)**2 + (phi - center_phi)**2)

    if (dist < radius) then
       CosBell = cos(0.5_dp * pi * dist / radius)**2
    else
       CosBell = 0.0_dp
    end if
    
  end function CosBell

  ! --- 元の分布関数 (Y_l^m) ---
  real(dp) function Distln(phi, sig, ln)
    real(dp), intent(in)     :: phi, sig
    character(2), intent(in) :: ln

    if ( ln == '11') then
       distln = sin(PI*sig) * cos(phi)**2
!       distln = sin(PI*sig) * cos(phi)
!       distln = 1.0D0
    else if ( ln == '21') then
       distln = 3*sin(PI*sig) * sin(phi) * cos(phi)**2
    else
       call MessageNotify('E','Distln','l,m values not supported')
    end if
    
  end function Distln

end program dcpam_main_spml_test