!----------------------------------------------------------------------
!     Copyright (c) 2002 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!表題  SPCPACK/ISPACK-F90
!
!履歴  2002/02/02  竹広真一 
!
module spcpack

  !use lumatrix
  use sppack
  use ctpack, g_rad => g_x, t_dr_t => t_dx_t, at_dr_at => at_dx_at

  implicit none
  private

  public spcinitial

  public g_lon, g_lat, g_lon_weight, g_lat_weight
  public l_nm, nm_l
  public gg_lon, gg_lat, g_rad
  public ggg_lon, ggg_lat, ggg_rad
  public yg_rad

  public y_gg, gg_y
  public at_dr_at, t_dr_t, ag_at, at_ag

  public ggg_yt, yt_ggg, ggg_yg, yg_ggg, yg_yt, yt_yg
  public yt_drad_yt, yt_divrad_yt, yt_lapla_yt
  public ggg_gradlon_yt, ggg_gradlat_yt
  public yt_divlon_ggg, yt_divlat_ggg, yt_div_ggg_ggg_ggg
  public yt_kxrgrad_yt, ggg_kgrad_yt, yt_l2_yt, yt_qoperator_yt
  public yt_radrot_ggg_ggg, yt_radrotrot_ggg_ggg_ggg
  public ytboundaries


  integer            :: im=64, jm=32, km=16  ! 格子点の設定(経度, 緯度, 動径)
  integer            :: nm=21, lm=16         ! 切断波数の設定(水平, 動径)
  real(8)            :: ri=0.0, ro=1.0       ! 球殻内外半径
  real(8), parameter :: pi=3.1415926535897932385D0

  real(8), dimension(:,:,:), allocatable :: ggg_lon, ggg_lat, ggg_rad
  real(8), dimension(:,:), allocatable   :: yg_rad

  save im, jm, km, nm, lm, ri, ro

  contains
  !--------------- 初期化 -----------------
    subroutine spcinitial(i,j,k,n,l,r_in,r_out)

     integer,intent(in) :: i, j, k        ! 格子点の設定(経度, 緯度, 動径)
     integer,intent(in) :: n, l           ! 切断波数の設定(水平, 動径)

     real(8),intent(in) :: r_in, r_out    ! 球殻内外半径

     im = i  ; jm = j ; km = k
     nm = n  ; lm = l
     ri = r_in ; ro = r_out

     call spinitial(nm,im,jm,km+1)
     call ctinitial(km,lm,r_in,r_out)

     allocate(ggg_lon(im,jm,0:km))
     allocate(ggg_lat(im,jm,0:km))
     allocate(ggg_rad(im,jm,0:km))
     allocate(yg_rad((nm+1)*(nm+1),0:km))

     ggg_lon = spread(gg_lon,3,km+1)
     ggg_lat = spread(gg_lat,3,km+1)
     ggg_rad = spread(spread(g_rad,1,jm),1,im)

     yg_rad = spread(g_rad,1,(nm+1)*(nm+1))

   end subroutine spcinitial

  !--------------- 基本変換 -----------------

    function ggg_yt(yt)  ! スペクトル -> 格子
      real(8), dimension(im,jm,0:km)                     :: ggg_yt
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: yt

      ggg_yt = gga_ya(ag_at(yt))

    end function ggg_yt

    function yt_ggg(ggg)  ! 格子 -> スペクトル
      real(8), dimension((nm+1)*(nm+1),0:lm)             :: yt_ggg
      real(8), dimension(im,jm,0:km), intent(in)         :: ggg

      yt_ggg = at_ag(ya_gga(ggg))

    end function yt_ggg

    function ggg_yg(yg)  ! 水平スペクトル -> 格子
      real(8), dimension(im,jm,0:km)                     :: ggg_yg
      real(8), dimension((nm+1)*(nm+1),0:km), intent(in) :: yg

      ggg_yg = gga_ya(yg)

    end function ggg_yg

    function yg_ggg(ggg)  ! 格子 -> 水平スペクトル
      real(8), dimension((nm+1)*(nm+1),0:km)             :: yg_ggg
      real(8), dimension(im,jm,0:km), intent(in)         :: ggg

      yg_ggg = ya_gga(ggg)

    end function yg_ggg

    function yg_yt(yt)  ! スペクトル -> 水平スペクトル
      real(8), dimension((nm+1)*(nm+1),0:km)             :: yg_yt
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: yt

      yg_yt = ag_at(yt)

    end function yg_yt

    function yt_yg(yg)  ! 水平スペクトル -> スペクトル
      real(8), dimension((nm+1)*(nm+1),0:km)             :: yt_yg
      real(8), dimension((nm+1)*(nm+1),0:km), intent(in) :: yg

      yt_yg = at_ag(yg)

    end function yt_yg

  !--------------- 微分計算 -----------------
    function yt_drad_yt(yt)     ! 動径微分
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: yt
      real(8), dimension((nm+1)*(nm+1),0:lm)             :: yt_drad_yt

      yt_drad_yt = at_dr_at(yt)

    end function yt_drad_yt

    function yt_divrad_yt(yt)    ! 動径発散微分 
                                 ! 1/r^2 ∂/∂r r^2 = ∂/∂r + 2/r 
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: yt
      real(8), dimension((nm+1)*(nm+1),0:lm)             :: yt_divrad_yt

      yt_divrad_yt = yt_drad_yt(yt) + yt_yg(2/yg_rad*yg_yt(yt))

    end function yt_divrad_yt

    function yt_lapla_yt(yt)     ! ラプラシアン
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: yt
      real(8), dimension((nm+1)*(nm+1),0:lm)             :: yt_lapla_yt

      yt_lapla_yt = yt_divrad_yt(yt_drad_yt(yt)) &
                   + yt_yg(yg_yt(ya_lapla_ya(yt))/yg_rad**2)

    end function yt_lapla_yt

    function ggg_gradlon_yt(yt) ! スペクトルに作用する勾配型経度微分
                                ! 1/rcosφ・∂/∂λ
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: yt
      real(8), dimension(im,jm,0:km)                     :: ggg_gradlon_yt

      ggg_gradlon_yt = gga_gradlon_ya(yg_yt(yt))/ggg_rad

    end function ggg_gradlon_yt

    function ggg_gradlat_yt(yt) ! スペクトルに作用する勾配型緯度微分
                                ! 1/r ∂/∂φ
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: yt
      real(8), dimension(im,jm,0:km)                     :: ggg_gradlat_yt

      ggg_gradlat_yt = gga_gradlat_ya(yg_yt(yt))/ggg_rad
    end function ggg_gradlat_yt

    function yt_divlon_ggg(ggg)   ! 格子に作用する発散型経度微分 
                                  ! 1/rcosφ・∂/∂λ

      real(8), dimension(im,jm,0:km), intent(in)   :: ggg
      real(8), dimension((nm+1)*(nm+1),0:lm)       :: yt_divlon_ggg

      yt_divlon_ggg = yt_yg(ya_divlon_gga(ggg/ggg_rad))
    end function yt_divlon_ggg

    function yt_divlat_ggg(ggg)   ! 格子に作用する発散型緯度微分
                                  ! 1/rcosφ・∂(f cosφ)/∂φ

      real(8), dimension(im,jm,0:km), intent(in)   :: ggg
      real(8), dimension((nm+1)*(nm+1),0:lm)       :: yt_divlat_ggg

      yt_divlat_ggg = yt_yg(ya_divlat_gga(ggg/ggg_rad))
    end function yt_divlat_ggg

    function yt_div_ggg_ggg_ggg(ggg_vlon,ggg_vlat,ggg_vrad) ! 発散

      real(8), dimension(im,jm,0:km), intent(in) :: ggg_vlon(:,:,:) ! 経度成分
      real(8), dimension(im,jm,0:km), intent(in) :: ggg_vlat(:,:,:) ! 緯度成分
      real(8), dimension(im,jm,0:km), intent(in) :: ggg_vrad(:,:,:) ! 動径成分
      real(8), dimension((nm+1)*(nm+1),0:lm)     :: yt_div_ggg_ggg_ggg

      yt_div_ggg_ggg_ggg =   yt_divlon_ggg(ggg_vlon) &
                           + yt_divlat_ggg(ggg_vlat) &
                           + yt_divrad_yt(yt_ggg(ggg_vrad))

    end function yt_div_ggg_ggg_ggg


  !--------------- ポロイダル/トロイダルモデル用微分 -----------------

    function yt_kxrgrad_yt(yt)    ! k×r・▽ = ∂/∂λ
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: yt
      real(8), dimension((nm+1)*(nm+1),0:lm)             :: yt_kxrgrad_yt

      yt_kxrgrad_yt =  ya_dlon_ya(yt)
    end function yt_kxrgrad_yt

    function ggg_kgrad_yt(yt)    ! k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: yt
      real(8), dimension(im,jm,0:km)                     :: ggg_kgrad_yt

      ggg_kgrad_yt =  cos(ggg_lat)*ggg_gradlat_yt(yt) &
                    + sin(ggg_lat)*ggg_yt(yt_drad_yt(yt))
    end function ggg_kgrad_yt

    function yt_l2_yt(yt)    ! L2 作用素 = -水平ラプラシアン 
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: yt
      real(8), dimension((nm+1)*(nm+1),0:lm)             :: yt_l2_yt

      yt_l2_yt = -ya_lapla_ya(yt)
    end function yt_l2_yt

    function yt_qoperator_yt(yt)    ! Q=(k・▽-1/2(L2 k・▽+ k・▽L2))
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: yt
      real(8), dimension((nm+1)*(nm+1),0:lm)             :: yt_qoperator_yt

      yt_qoperator_yt = &
             yt_ggg(ggg_kgrad_yt(yt) - ggg_kgrad_yt(yt_l2_yt(yt))/2) &
           - yt_l2_yt(yt_ggg(ggg_kgrad_yt(yt)))/2

    end function yt_qoperator_yt

    function yt_radrot_ggg_ggg(ggg_vlon,ggg_vlat)  ! r・(▽×v)
      real(8), dimension(im,jm,0:km), intent(in) :: ggg_vlon
      real(8), dimension(im,jm,0:km), intent(in) :: ggg_vlat
      real(8), dimension((nm+1)*(nm+1),0:lm)     :: yt_radrot_ggg_ggg

      yt_radrot_ggg_ggg = yt_yg(ya_divlon_gga(ggg_vlat) &
                                - ya_divlat_gga(ggg_vlon))
      
    end function yt_radrot_ggg_ggg

    function yt_radrotrot_ggg_ggg_ggg(ggg_vlon,ggg_vlat,ggg_vrad) 
                                                  ! r・(▽×▽×v)
      real(8), dimension(im,jm,0:km), intent(in) :: ggg_vlon
      real(8), dimension(im,jm,0:km), intent(in) :: ggg_vlat
      real(8), dimension(im,jm,0:km), intent(in) :: ggg_vrad
      real(8), dimension((nm+1)*(nm+1),0:lm)     :: yt_radrotrot_ggg_ggg_ggg

      yt_radrotrot_ggg_ggg_ggg = &
           yt_divrad_yt(yt_yg( &
              (ya_divlon_gga(ggg_vlon)+ ya_divlat_gga(ggg_vlat))/yg_rad)) &
         - ya_lapla_ya(yt_ggg(ggg_vrad/ggg_rad**2))

    end function yt_radrotrot_ggg_ggg_ggg

  !--------------- 境界値問題 -----------------

    subroutine ytboundaries(yt,values,cond)   ! ディリクレ, ノイマン条件
      ! Chebyshev 空間での境界条件適用

      real(8), dimension((nm+1)*(nm+1),0:lm),intent(inout)      :: yt
              ! 境界条件を適用するデータ. 修正された値を返す. 

      real(8), dimension((nm+1)*(nm+1),2), intent(in), optional :: values
              ! 境界での 値/勾配 分布を水平スペクトル変換したものを与える. 
              ! 省略時は値/勾配 0 となる. 

      character(len=2), intent(in), optional             :: cond
              ! 境界条件. 省略時は 'DD'
              !   DD : 両端ディリクレ
              !   DN,ND : ディリクレ/ノイマン条件
              !   NN : 両端ノイマン

      if (.not. present(cond)) then
         if (present(values)) then
            call ctboundaries_DD(yt,values)
         else
            call ctboundaries_DD(yt)
         endif
         return
      endif

      select case(cond)
      case ('NN')
         if (present(values)) then
            call ctboundaries_NN(yt,values)
         else
            call ctboundaries_NN(yt)
         endif
      case ('DN')
         if (present(values)) then
            call ctboundaries_DN(yt,values)
         else
            call ctboundaries_DN(yt)
         endif
      case ('ND')
         if (present(values)) then
            call ctboundaries_ND(yt,values)
         else
            call ctboundaries_ND(yt)
         endif
      case ('DD')
         if (present(values)) then
            call ctboundaries_DD(yt,values)
         else
            call ctboundaries_DD(yt)
         endif
      case default
         call msgdmp('E','ytboundaries','B.C. not supported')
      end select

    end subroutine ytboundaries

  end module spcpack
