!----------------------------------------------------------------------
!     Copyright (c) 2002 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  SPCPACK/ISPACK-F90 
!      ポロイダル磁場拡散問題
!
!履歴  2002/02/07  竹広真一
!
program spcdiff3

  use lumatrix
  use spcpack
  use gtool_history
  implicit none

 !---- 空間解像度設定 ----
  integer,parameter  :: im=32, jm=16, km=16  ! 格子点の設定(経度, 緯度, 動径)
  integer,parameter  :: nm=10, lm=16         ! 切断波数の設定(水平, 動径)

 !---- 座標変数など ----
  real(8),parameter  :: ri=0.5, ro=1.5       ! 内外半径

 !---- 変数 ----
  real(8), dimension(im,jm,0:km)     :: ggg_pmag   ! ポロイダル磁場
  real(8), dimension((nm+1)**2,0:lm) :: yt_pmag    ! ポロイダル磁場

  real(8)             :: gg_pmagbndry(im,jm,2)     ! 境界値
  real(8)             :: y_pmagbndry((nm+1)**2,2)  ! 境界値

 !---- 時間積分パラメター ----
  real(8), parameter :: dt=1e-5                 ! 時間ステップ間隔
  integer, parameter :: nt=10000, ndisp=1000    ! 時間積分数, 表示ステップ

 !---- 物理パラメター ----
  character(len=2), parameter :: tempbc='DN'    ! 温度境界条件(DD/DN/ND/NN)
  real(8), parameter :: nu=1.0                  ! 粘性係数

  real(8), parameter :: pi=3.1415926535897932385D0
  real(8), parameter :: lon1=pi                 ! 初期分布緯度座標
  real(8), parameter :: lat1=0.0                ! 初期分布経度座標
  real(8), parameter :: rad1=(ri+ro)/2.0        ! 初期分布動径座標
  real(8), parameter :: sigma=0.2               ! 初期分布の大きさ

 !---- その他 ----
  integer :: it


 !---------------- 座標値の設定 ---------------------
  call spcinitial(im,jm,km,nm,lm,ri,ro)

 !------------------- 初期値設定 ----------------------
  ggg_pmag =  exp(-((ggg_lon-lon1)**2+(ggg_lat-lat1)**2+(ggg_rad-rad1)**2)&
                 /(2*sigma**2))
  yt_pmag = yt_ggg(ggg_pmag)


 !------------------- 境界値設定 ----------------------
  gg_pmagbndry(:,:,1) = 0
  gg_pmagbndry(:,:,2) = 0

  y_pmagbndry(:,1) = y_gg(gg_pmagbndry(:,:,1))
  y_pmagbndry(:,2) = y_gg(gg_pmagbndry(:,:,2))

  call output_gtool4_init
  call output_gtool4

 !------------------- 時間積分 ----------------------
  do it=1,nt                                         ! Euler 法による時間積分
     yt_pmag = yt_pmag + dt *( nu * yt_lapla_yt(yt_pmag) )
     call ytpmagboundaries(yt_pmag)

     if(mod(it,ndisp) .eq. 0)then                    ! 出力
        ggg_pmag = ggg_yt(yt_pmag)
        call output_gtool4
     endif
  enddo
  call output_gtool4_close

contains
  subroutine ytpmagboundaries(yt) ! スペクトル空間でのポロイダル磁場境界条件
    ! 周囲が非電気伝導体の場合
    !     dBP_l^m/dr + (l+1)BP_l^m/r = 0 at r=r_o
    !     dBP_l^m/dr - l BP_l^m/r = 0    at r=r_i

    real(8), dimension((nm+1)*(nm+1),0:lm),intent(inout)  :: yt ! データ

    real(8), dimension(:,:,:), allocatable  :: alu
    integer, dimension(:,:), allocatable    :: kp
    real(8), dimension(0:lm,0:lm)           :: tt_data
    real(8), dimension(0:lm,0:km)           :: tg_f
    real(8), dimension(0:lm,0:km)           :: tg_dfdr

    logical :: first = .true.
    integer :: l, nwave(2), m, nn
    save    :: alu, kp, first

    if ( first ) then
       first = .false.

       allocate(alu((nm+1)*(nm+1),0:lm,0:lm),kp((nm+1)*(nm+1),0:lm))

       tt_data = 0
       do l=0,lm
          tt_data(l,l)=1
       enddo
       alu = spread(tt_data,1,(nm+1)*(nm+1))

       tg_f    = ag_at(tt_data)
       tg_dfdr = ag_at(at_dr_at(tt_data))

       do nn=1,(nm+1)*(nm+1)
          nwave = nm_l(nn)
          l = nwave(1)
          alu(nn,lm-1,:) = tg_dfdr(:,0)  + (l+1) * tg_f(:,0)/g_rad(0)
          alu(nn,lm,:)   = tg_dfdr(:,km) - l * tg_f(:,km)/g_rad(km)
       end do

       call ludecomp(alu,kp)
    end if

    yt(:,km-1) = 0
    yt(:,km)   = 0
    yt = lusolve(alu,kp,yt)

  end subroutine ytpmagboundaries

  subroutine output_gtool4_init
    call HistoryCreate( &                                  ! ヒストリー作成
           file='spcdiff3.nc', title='Diffusion model in a spherical shell', &
           source='Sample program of gtool_history/gtool4', &
           institution='GFD_Dennou Club davis project',     &
           dims=(/'lon','lat','rad','t  '/), dimsizes=(/im,jm,km+1,0/),&
           longnames=(/'Latitude ','Longitude','Radius   ','time     '/),&
           units=(/'1','1','1','1'/),   &
           origin=0.0, interval=real(ndisp*dt) )

    call HistoryPut('lon',g_lon/pi*180)                       ! 変数出力
    call HistoryPut('lat',g_lat/pi*180)                       ! 変数出力
    call HistoryPut('rad',g_rad)                              ! 変数出力

    call HistoryAddVariable( &                                ! 変数定義
           varname='pmag', dims=(/'lon','lat','rad','t  '/), & 
           longname='poloidal magnetic potential', units='1', xtype='double')
  end subroutine output_gtool4_init

  subroutine output_gtool4
    write(6,*) 'it = ',it
    call HistoryPut('t',real(it*dt))
    call HistoryPut('pmag',ggg_pmag)
  end subroutine output_gtool4

  subroutine output_gtool4_close
    call HistoryClose
  end subroutine output_gtool4_close

end program spcdiff3
