!----------------------------------------------------------------------
!  Copyright (C) 2001--2018 SPMODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
! Sample program for gt4_history/gt4f90io and ISPACK   2002/08/21 S.Takehiro
!                                                      2004/01/26 M.Odaka
!                                                      2010/03/09 S.Takehiro
!                                                      2018/07/18 S.Takehiro
!                                                      2020/10/24 R.Suzuki spml2 mpi用に修正
!
! Solving a linear 2-D shallow water system on a sphere 
!     with an isolated mountain
!     du/dt + u0/a\cos\phi du/d\lambda + v/a du0/d\phi - u0v\tan\phi/a 
!           - 2\Omega\sin\phi v = -g/a\cos\phi dh/d\lambda - \nu\lapla^4 u,
!     dv/dt + u0/a\cos\phi dv/d\lambda + 2 u0 u \tan\phi/a 
!           + 2\Omega\sin\phi u  =  -g/a  dh/d\phi - \nu\lapla^4 v,
!     dh/dt + ( 1/\cos\phi d( (H+h0)u+(h-ht)u0 )/d\lambda 
!               + 1/\cos\phi d( (H+h0)v\cos\phi)/d\phi ) = - \nu\lapla^4 h.
!
! A setup is similar to the experiment of Grose and Hoskins (1979) 
! with a superrotating rigid-rotation zonal wind profile. 
!
  program spshallow_topo_linear
    use mpi           ! add rsuzuki 
    use w_mpi_module  ! add rsuzuki
!    use dc_string, only : StoA, Printf ! add rsuzuki
!    use w_module
    use dc_types      ! add rsuzuki ファイルを分割出力する際，
                      ! 必要となる文字列cip を保持するために必要
    use dc_string, only : StoA, Printf ! add rsuzuki
    use gt4_history
    implicit none

  !---- 空間解像度設定 ----
!!$    integer, parameter :: im=512, jm=256             ! 格子点の設定(X,Y)
!!$    integer, parameter :: nm=170
    integer, parameter :: im=128, jm=64             ! 格子点の設定(X,Y)
    integer, parameter :: nm=42
!!$    integer, parameter :: im=96, jm=48             ! 格子点の設定(X,Y)
!!$    integer, parameter :: nm=31
!    integer, parameter :: im=64, jm=32             ! 格子点の設定(X,Y)
!    integer, parameter :: nm=21

  !---- 変数 ----         
    real(8), allocatable :: xv_U(:,:)          ! 格子点データ(速度経度成分)
    real(8), allocatable :: xv_V(:,:)          ! 格子点データ(速度緯度成分)
    real(8), allocatable :: xv_H(:,:)          ! 格子点データ(変位)
!    real(8) :: xv_U(0:im-1,jc)          ! 格子点データ(速度経度成分)
!    real(8) :: xv_V(0:im-1,jc)          ! 格子点データ(速度緯度成分)
!    real(8) :: xv_H(0:im-1,jc)          ! 格子点データ(変位)

    real(8), allocatable :: w_U(:)       ! スペクトル(速度経度成分)
    real(8), allocatable :: w_V(:)       ! スペクトル(速度緯度成分)
    real(8), allocatable :: w_H(:)       ! スペクトル(変位)
!    real(8)            :: w_U((nm+1)*(nm+1))       ! スペクトル(速度経度成分)
!    real(8)            :: w_V((nm+1)*(nm+1))       ! スペクトル(速度緯度成分)
!    real(8)            :: w_H((nm+1)*(nm+1))       ! スペクトル(変位)

    real(8), allocatable :: xv_U0(:,:)         ! 格子点データ(速度経度成分)
    real(8), allocatable :: xv_H0(:,:)         ! 格子点データ(変位)
!    real(8) :: xv_U0(0:im-1,jc)         ! 格子点データ(速度経度成分)
!    real(8) :: xv_H0(0:im-1,jc)         ! 格子点データ(変位)

    real(8), allocatable :: w_U0(:)      ! スペクトル(速度経度成分)
    real(8), allocatable :: w_H0(:)      ! スペクトル(変位)
!    real(8)            :: w_U0((nm+1)*(nm+1))      ! スペクトル(速度経度成分)
 !   real(8)            :: w_H0((nm+1)*(nm+1))      ! スペクトル(変位)

!    real(8) :: xv_Zeta(0:im-1,jc)       ! 格子点データ(渦度)
!    real(8) :: xv_Htopo(0:im-1,jc)      ! 格子点データ(地形)
    real(8), allocatable :: xv_Zeta(:,:)       ! 格子点データ(渦度)
    real(8), allocatable :: xv_Htopo(:,:)      ! 格子点データ(地形)

  !---- 時間積分パラメター ----
    real(8), parameter :: dt = 100                 ! 時間ステップ間隔
    integer, parameter :: nt = 15000, ndisp = 500  ! 時間積分数, 表示ステップ
 
  !---- 物理パラメター ----
    real(8), parameter :: pi = 3.141592653589793   ! 円周率
    real(8), parameter :: R0 = 6.371D6             ! 球の半径
    real(8), parameter :: H0 = 10D3                ! 平均変位
    real(8), parameter :: Omega = 7.292D-5         ! 回転角速度
    real(8), parameter :: Grav = 9.8               ! 重力加速度

    real(8), parameter :: U0   = 15                ! 平均流
    real(8), parameter :: Lon0 = pi, Lat0=pi/6     ! 山の中心位置
    real(8), parameter :: W0 = 8                   ! 山の広がりの逆数
    real(8), parameter :: Htopo0= 2500             ! 山の高さ

    integer, parameter :: ndiff=4                  ! 超粘性の次数
    real(8), parameter :: Nu=2.338D16              ! 超粘性係数

  !------ 作業変数 ------
    integer            :: i, j, it=0
    real(8)            :: Delta
  
  !------ MPI 用変数 ------ add rsuzuki
    integer :: np, ip, ierr  ! MPI 用変数 np:プロセス数 ip:プロセスID
    character(len=STRING) :: cip ! ip 文字変数

  !------ MPI ヒストリー出力 ----- add rsuzuki 
    character(len=100) :: hst_file

  !------ 時間計測 ----- add rsuzuki
    real(8) :: TIM0, TIM1 ! start, end time

  !------ OpenMP 並列化のスレッド数の取得 ---- add rsuzuki
    integer(8) :: NTH

  !------ MPI スタート----- add rsuzuki
    call MPI_INIT(IERR)
    call MPI_COMM_RANK(MPI_COMM_WORLD,IP,IERR)
    call MPI_COMM_SIZE(MPI_COMM_WORLD,NP,IERR)

  !------ OpenMP 並列化のスレッド数の出力 ---- add rsuzuki
    call MXGOMP(NTH) 
    print *, 'NUMBER OF THREADS :', NTH 

  !----- MPI 用出力ファイル----- add rsuzuki
    call Printf(unit=cip,fmt="%04d",i=(/ip/))  
    hst_file='spshallow_topo_linear-ip'//trim(adjustl(cip))//'.nc'

  !---------------- 座標値の設定 ---------------------
    call w_mpi_Initial(nm,im,jm)                ! ISPACK初期化 w_Initial -> w_mpi_Initial 

  !------ 変数の割り付け----- add rsuzuki
    allocate(xv_U(0:im-1,jc), xv_V(0:im-1,jc), xv_H(0:im-1,jc))
    allocate(xv_U0(0:im-1,jc), xv_H0(0:im-1,jc))
    allocate(xv_Zeta(0:im-1,jc), xv_Htopo(0:im-1,jc))
 
    allocate(w_U(nc), w_V(nc), w_H(nc))
    allocate(w_U0(nc), w_H0(nc))

  !------ 計算開始時間の取得 ------ add rsuzuki
    call MXTIME(TIM0)

  !------------------- 初期値設定 ----------------------
    xv_U0  = U0*cos(xv_Lat)
    xv_H0  = ( Omega*R0*U0/(2*Grav) + U0**2/(4*Grav) )*cos(2*xv_Lat)

    xv_U  = 0 ; xv_V  = 0 ; xv_H  = 0

    w_U0 = w_xv(xv_U0) ; w_H0 = w_xv(xv_H0)
    w_U = w_xv(xv_U) ; w_V = w_xv(xv_V) ; w_H = w_xv(xv_H)

  !------------------- 地形設定 ----------------------

    do j=1,jc
       do i=0,im-1
          delta = acos( cos(xv_Lat(i,j))*cos(Lat0)*cos(xv_Lon(i,j)-Lon0) &
                       + sin(xv_Lat(i,j))*sin(Lat0) )
          if ( abs(delta) <= pi/W0 ) then 
             xv_Htopo(i,j)= htopo0 * (1+cos(W0*Delta))/2
          else
             xv_Htopo(i,j)= 0
          endif
       enddo
    enddo

    call output_gtool4_init                        ! ヒストリー初期化
    call output_gtool4                             ! 初期値出力

  !------------------- 時間積分 ----------------------
    do it=1,nt
       w_U = ( w_U &
               + dt * w_xv( - xv_U0 * xv_GradLon_w(w_U) / R0   &
                            - xv_V  * xv_GradLat_w(w_U0) / R0  &
                            + xv_U0 * xv_V * tan(xv_Lat) / R0  &
                            + 2 * Omega * sin(xv_Lat) * xv_V   &
                           - Grav * xv_GradLon_w(w_H)/ R0   ) &
             )/(1+Nu*(-rn(:,1)/R0**2)**(ndiff/2)*dt)

       xv_U = xv_w(w_U)

       w_V = ( w_V &
              + dt * w_xv( - xv_U0 * xv_GradLon_w(w_V) / R0   &
                           - xv_U * xv_U0 * tan(xv_Lat) / R0  &
                           - xv_U0 * xv_U * tan(xv_Lat) / R0  &
                           - 2 * Omega * sin(xv_Lat) * xv_U   &
                           - Grav * xv_GradLat_w(w_H) / R0 )  &
             )/(1+Nu*(-rn(:,1)/r0**2)**(ndiff/2)*dt)

       xv_V = xv_w(w_V)

       w_H = ( w_H &
              + dt * ( - w_Div_xv_xv( xv_U*(H0+xv_H0)+xv_U0*(xv_H-xv_Htopo),   &
                                      xv_V*(H0+xv_H0) ) / R0      ) &
              )/(1+Nu*(-rn(:,1)/R0**2)**(ndiff/2)*dt)

       xv_H = xv_w(w_H)

       if(mod(it,ndisp) .eq. 0)then                        ! 出力
          call output_gtool4
       endif
    enddo

  !------ 計算終了時間の取得 ------ add rsuzuki
    call MXTIME(TIM1)

    print *, '計算時間[sec] : ', TIM1-TIM0

    call output_gtool4_close
    call w_mpi_Finalize()   ! add rsuzuki
    call MPI_FINALIZE(IERR) ! add rsuzuki
    stop

  contains
    subroutine output_gtool4_init
      call HistoryCreate( &                               ! ヒストリー作成
           file=trim(hst_file), & ! add rsuzuki
!           file='spshallow_topo_linear.nc',                         &
           title='Shallow water equation on a sphere',            &
           source='Sample program of gtool_History/gtool4',       &
           institution='GFD_Dennou Club davis/spmodel project',   &
           dims=(/'lon','lat','t  '/), dimsizes=(/im,jc,0/),      &
           longnames=(/'longitude','latitude ','time     '/),     &
           units=(/'degree_west ','degree_north','s           '/),&
           origin=0.0, interval=real(ndisp*dt) )

      call HistoryPut('lon',x_Lon*180/pi)              ! 変数出力
      call HistoryAddattr('lon','topology','circular')   ! 周期属性
      call HistoryAddattr('lon','modulo',360.0)          ! 周期属性
      call HistoryPut('lat',v_Lat*180/pi)              ! 変数出力 
                        ! y_Lat -> v_Lat  add rsuzuki
      call HistoryAddVariable( &                          ! 変数定義
           varname='h', dims=(/'lon','lat','t  '/), & 
           longname='surface displacement ', units='m', xtype='double')
      call HistoryAddVariable( &                          ! 変数定義
           varname='u', dims=(/'lon','lat','t  '/), & 
           longname='velocity(longitude) ', units='m/s', xtype='double')
      call HistoryAddVariable( &                          ! 変数定義
           varname='v', dims=(/'lon','lat','t  '/), & 
           longname='velocity(latitude) ', units='m/s', xtype='double')

      call HistoryAddVariable( &                          ! 変数定義
           varname='zeta', dims=(/'lon','lat','t  '/), & 
           longname='vorticity', units='1/s', xtype='double')

      call HistoryAddVariable( &                          ! 変数定義
           varname='htopo', dims=(/'lon','lat'/), & 
           longname='orography', units='m', xtype='double')
      call HistoryPut('htopo',xv_Htopo)
    end subroutine output_gtool4_init

    subroutine output_gtool4
      
      write(6,*) 'it = ',it, 'ip = ', ip  ! add rsuzuki ip
      call HistoryPut('t',real(it*dt))
      call HistoryPut('u',xv_U)
      call HistoryPut('v',xv_V)
      call HistoryPut('h',xv_H)
      xv_Zeta = xv_w(w_Divlon_xv(xv_V) - w_Divlat_xv(xv_U))/r0
      call HistoryPut('zeta',xv_Zeta)
    end subroutine output_gtool4

    subroutine output_gtool4_close
      call HistoryClose
    end subroutine output_gtool4_close

  end program spshallow_topo_linear
