!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Program init
!
!   * Developers: Morikawa Yasuhiro
!   * Version: $Id: init.f90,v 1.7 2005/01/03 09:34:35 morikawa Exp $
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview
!
!Generate Various Initial Data
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!
!
                                                                 !=end

program init
                                                                 !=begin
  !== Dependency
  !----- 構造体参照モジュール -----
  use type_mod, only: STRING, INTKIND, REKIND, DBKIND
  use nmlfile_mod, only: nmlfile_init
  use constants_mod, only: constants_init, pi, R0

  !----- 格子点取得モジュール -----
  use grid_3d_mod,         only: grid_3d_init, im, jm, km, grid_3d_end
  use grid_wavenumber_mod, only: grid_wavenumber_init, nm, grid_wavenumber_end

  !----- 座標データ生成モジュール -----
  use axis_type_mod, only: AXISINFO
  use axis_x_mod, only: axis_x_init, axis_x_weight, axis_x_spectral, axis_x_end
  use axis_y_mod, only: axis_y_init, axis_y_weight, axis_y_spectral, axis_y_end
  use axis_z_mod, only: axis_z_init, axis_z_sigmahalf_manual, axis_z_end

  !----- データI/Oモジュール -----
  use io_gt4_out_mod, only : io_gt4_out_init   , io_gt4_out_SetDims, &
       &                     io_gt4_out_SetVars, io_gt4_out_Put    , &
       &                     io_gt4_out_end

  !----- 時刻管理モジュール -----
  use time_mod, only: time_init, tvar, ttype, tname, tunit, time_end

  !----- SPMODEL モジュール (本当はこんな上から呼ぶべきではない) -----
  use spml_mod,  only: spml_init, xy_Lon, xy_Lat, xy_w, w_Div_xy_xy &
       &              ,wa_LaplaInv_wa, wa_xya, xya_GradLat_wa, xya_GradLon_wa

  !----- デバッグ・汎用ツール -----
  use dc_trace,  only: SetDebug, DbgMessage, BeginSub, EndSub
  use dc_string, only: toChar
  use dc_error,  only: StoreError, USR_ECHAR
                                                                 !=end
  implicit none

  !-------------------------------------------------------------------
  !   変数定義
  !-------------------------------------------------------------------
  type(AXISINFO)::          &
       & x_Lon            , & ! 格子点座標(経度)
       & x_Lon_Weight     , & ! 格子点座標(経度重み)
       & y_Lat            , & ! 格子点座標(緯度)
       & y_Lat_Weight     , & ! 格子点座標(緯度重み)
       & z_Sigma          , & ! 格子点座標(σレベル(整数))
       & z_SigmaHalf          ! 格子点座標(σレベル(半整数))


  real(DBKIND), pointer     ::  &
       & xyz_VLon(:,:,:)      , & ! 格子点データ(速度経度成分)
       & xyz_VLat(:,:,:)      , & ! 格子点データ(速度緯度成分)
       & xyz_Vor(:,:,:)       , & ! 格子点データ(渦度)
       & xyz_Div(:,:,:)       , & ! 格子点データ(発散)
       & xyz_Temp(:,:,:)      , & ! 格子点データ(温度)
       & xyz_QVap(:,:,:)      , & ! 格子点データ(比湿)
       & 
       & xy_Ps(:,:)               ! 格子点データ(地表面気圧)

  !----- 作業用内部変数 -----
  integer(INTKIND)   :: k

  !----- 「お山を1つ」で利用 -----
  real(DBKIND) :: PsTop    ! Ps の最大値
  real(DBKIND) :: Lat_Ps   ! 緯度の裾野
  real(DBKIND) :: Lon_Ps   ! 経度の裾野
  real(DBKIND) :: R_Ps     ! 半径

  !----- 「温度の山を1つ」で利用 -----
  real(DBKIND) :: TempTop    ! Temp の最大値
  real(DBKIND) :: Lat_Temp   ! 緯度の裾野
  real(DBKIND) :: Lon_Temp   ! 経度の裾野
  real(DBKIND) :: R_Temp     ! 半径


  !----- 「剛体回転流」で利用 -----
  real(DBKIND), allocatable :: xy_UVFact(:,:) ! u→U, v→V のファクター
  real(DBKIND), allocatable :: wz_Psi(:,:)    ! スペクトル(流線関数)
  real(DBKIND), allocatable :: wz_Chi(:,:)    ! スペクトル(ポテンシャル)


  character(STRING),  parameter:: version = &
       & '$Id: init.f90,v 1.7 2005/01/03 09:34:35 morikawa Exp $'
  character(STRING),  parameter:: tagname = '$Name:  $'
  character(STRING),  parameter:: subname = "init"

  !-------------------------------------------------------------------
  !   デバッグモード開始
  !-------------------------------------------------------------------
  call SetDebug

  call BeginSub(subname)

  !----------------------------------------------------------------
  !   Version identifier
  !----------------------------------------------------------------
  call DbgMessage('%c :: %c', c1=trim(version), c2=trim(tagname))

  !-------------------------------------------------------------------
  !   NAMELIST file Setting
  !-------------------------------------------------------------------
  call nmlfile_init('init.nml')

  !-------------------------------------------------------------------
  !   初期化
  !-------------------------------------------------------------------
  call constants_init        ! 定数の設定
  call grid_3d_init          ! 空間格子点数の取得
  call grid_wavenumber_init  ! 波数格子点数の取得
  call axis_x_init           ! 座標軸データの生成
  call axis_y_init           ! 座標軸データの生成
  call axis_z_init           ! 座標軸データの生成
  call io_gt4_out_init       ! データ出力の初期設定

  ! 初期値生成プログラムを別モジュールにした場合はここから外せる。
  call constants_init          ! 物理定数の取得
  call spml_init             ! SPMODEL の初期設定

  !-------------------------------------------------------------------
  !  軸データ生成
  !-------------------------------------------------------------------
!!$  call axis_x_weight(x_Lon_Weight)    ! 経度座標重みデータ取得
  call axis_x_spectral(x_Lon)         ! 経度座標データ取得
!!$  call axis_y_weight(y_Lat_Weight)    ! 緯度座標重みデータ取得
  call axis_y_spectral(y_Lat)         ! 緯度座標データ取得
  call axis_z_sigmahalf_manual &
       & ( z_Sigma           , &      ! 整数σレベル座標データ取得
       &   z_SigmaHalf  )             ! 半整数σレベル座標データ取得

  !-------------------------------------------------------------------
  !   出力用の軸データ設定
  !-------------------------------------------------------------------
  call io_gt4_out_SetDims(x_Lon)        ! 経度座標重みデータ取得
!!$  call io_gt4_out_SetDims(x_Lon_Weight) ! 経度座標データ取得
  call io_gt4_out_SetDims(y_Lat)        ! 緯度座標重みデータ取得
!!$  call io_gt4_out_SetDims(y_Lat_Weight) ! 緯度座標データ取得
  call io_gt4_out_SetDims(z_Sigma)      ! 整数σレベル座標データ取得
  call io_gt4_out_SetDims(z_SigmaHalf)  ! 半整数σレベル座標データ取得

  !-------------------------------------------------------------------
  !   出力用の変数データ設定
  !-------------------------------------------------------------------
  call io_gt4_out_SetVars('VLon')
  call io_gt4_out_SetVars('VLat')
  call io_gt4_out_SetVars('Vor')
  call io_gt4_out_SetVars('Div')
  call io_gt4_out_SetVars('Temp')
  call io_gt4_out_SetVars('QVap')
  call io_gt4_out_SetVars('Ps')

  !-------------------------------------------------------------------
  !   格子点データ作成
  !-------------------------------------------------------------------
  ! 平均値を代入
  allocate( xyz_VLon(im,jm,km) )
  allocate( xyz_VLat(im,jm,km) )
  allocate( xyz_Vor(im,jm,km)  )
  allocate( xyz_Div(im,jm,km)  )
  allocate( xyz_Temp(im,jm,km) )
  allocate( xyz_QVap(im,jm,km) )
  allocate( xy_Ps(im,jm)       )

  ! 等温無風
  xyz_VLon(:,:,:) = 0.0d0
  xyz_VLat(:,:,:) = 0.0d0
  xyz_Vor(:,:,:)  = 0.0d0
  xyz_Div(:,:,:)  = 0.0d0
  xyz_Temp(:,:,:) = 273.0d0
  xyz_QVap(:,:,:) = 0.0d0
  xy_Ps(:,:)      = 100000.0d0

!!$  ! お山を1つ
!!$  PsTop  = -200.0d0  ! Ps の最大値
!!$  Lat_Ps = pi/6.0   ! 緯度の裾野
!!$  Lon_Ps = 0.5*pi   ! 経度の裾野
!!$  R_Ps   = pi/9.0   ! 半径
!!$
!!$  xy_Ps(:,:) = xy_Ps(:,:) &
!!$       & + PsTop * (1.0D0 - sqrt(  &
!!$       &     min( R_Ps**2 ,                       &
!!$!         &       (x_Lon%a_Dim(:)-Lon_Ps)**2       &
!!$!         &         + (y_Lat%a_Dim(:)-Lat_Ps)**2 ) &
!!$       &         (xy_Lon(:,:)-Lon_Ps)**2       &
!!$       &           + (xy_Lat(:,:)-Lat_Ps)**2 ) &
!!$       &                         ) / R_Ps  )

  ! 東西風 10 m
!!$  do k = 1, km
!!$     xyz_VLon(:,:,k) = 15.0d0 * cos(xy_Lat(:,:))
!!$  enddo

!!$  ! 剛体回転流
!!$  allocate( xy_UVFact(im, jm) )
!!$  allocate( wz_Psi((nm+1)*(nm+1), km) )
!!$  allocate( wz_Chi((nm+1)*(nm+1), km) )
!!$
!!$  xy_UVFact(:,:) = cos( xy_Lat(:,:) )
!!$
!!$  do k = 1, km
!!$     xyz_VLon(:,:,k) = 100.0d0 * cos( xy_Lat(:,:) )
!!$
!!$     ! 風速から渦度発散を生成
!!$     xyz_Vor(:,:,k) = &
!!$          & xy_w(                      &
!!$          &      w_Div_xy_xy(          &
!!$          &         xyz_VLat(:,:,k)  , &
!!$          &         - xyz_VLon(:,:,k)  &
!!$          &        ) / R0              &
!!$          &     )
!!$
!!$     xyz_Div(:,:,k) = &
!!$          & xy_w(                     &
!!$          &      w_Div_xy_xy(         &
!!$          &         xyz_VLon(:,:,k) , &
!!$          &         xyz_VLat(:,:,k)   &
!!$          &        ) / R0             &
!!$          &     )
!!$
!!$  enddo
!!$
!!$  ! 試しに戻してみただけ。
!!$  !wz_Psi(:,:) = wa_LaplaInv_wa(  wa_xya( xyz_Vor(:,:,:) )  ) * R0**2
!!$  !wz_Chi(:,:) = wa_LaplaInv_wa(  wa_xya( xyz_Div(:,:,:) )  ) * R0**2
!!$  !
!!$  !xyz_VLon(:,:,:) = (  xya_GradLon_wa( wz_Chi(:,:) ) &
!!$  !     &                - xya_GradLat_wa( wz_Psi(:,:) )  ) / R0
!!$  !
!!$  !xyz_VLat(:,:,:) = (  xya_GradLon_wa( wz_Psi(:,:) ) &
!!$  !     &                + xya_GradLat_wa( wz_Chi(:,:) )  ) / R0


  ! 温度のお山
  TempTop  = 5.0d0  ! Temp の最大値
  Lat_Temp = pi/6.0   ! 緯度の裾野
  Lon_Temp = 0.5*pi   ! 経度の裾野
  R_Temp   = pi/9.0   ! 半径

  do k = 1, km
     xyz_Temp(:,:,k) = xyz_Temp(:,:,k) &
          & + TempTop * (1.0D0 - sqrt(  &
          &     min( R_Temp**2 ,                       &
!          &       (x_Lon%a_Dim(:)-Lon_Temp)**2       &
!          &         + (y_Lat%a_Dim(:)-Lat_Temp)**2 ) &
          &         (xy_Lon(:,:)-Lon_Temp)**2       &
          &           + (xy_Lat(:,:)-Lat_Temp)**2 ) &
          &                         ) / R_Temp  )
  enddo

  !-------------------------------------------------------------------
  !   データの出力 (t-Δt)
  !-------------------------------------------------------------------
  call io_gt4_out_Put('VLon' , real(xyz_VLon(:,:,:), REKIND) )
  call io_gt4_out_Put('VLat' , real(xyz_VLat(:,:,:), REKIND) )
  call io_gt4_out_Put('Vor'  , real(xyz_Vor(:,:,:) , REKIND) )
  call io_gt4_out_Put('Div'  , real(xyz_Div(:,:,:) , REKIND) )
  call io_gt4_out_Put('Temp' , real(xyz_Temp(:,:,:), REKIND) )
  call io_gt4_out_Put('QVap' , real(xyz_QVap(:,:,:), REKIND) )
  call io_gt4_out_Put('Ps'   , real(xy_Ps(:,:)     , REKIND) )


  !-------------------------------------------------------------------
  !   データの出力 (t)
  !-------------------------------------------------------------------
  call io_gt4_out_Put('VLon' , real(xyz_VLon(:,:,:), REKIND) )
  call io_gt4_out_Put('VLat' , real(xyz_VLat(:,:,:), REKIND) )
  call io_gt4_out_Put('Vor'  , real(xyz_Vor(:,:,:) , REKIND) )
  call io_gt4_out_Put('Div'  , real(xyz_Div(:,:,:) , REKIND) )
  call io_gt4_out_Put('Temp' , real(xyz_Temp(:,:,:), REKIND) )
  call io_gt4_out_Put('QVap' , real(xyz_QVap(:,:,:), REKIND) )
  call io_gt4_out_Put('Ps'   , real(xy_Ps(:,:)     , REKIND) )

  !-------------------------------------------------------------------
  !   終了処理
  !-------------------------------------------------------------------
  call grid_3d_end          ! 空間格子点数の取得
  call grid_wavenumber_end  ! 波数格子点数の取得
  call axis_x_end           ! 座標軸データの生成
  call axis_y_end           ! 座標軸データの生成
  call axis_z_end           ! 座標軸データの生成
  call io_gt4_out_end       ! データ出力の初期設定

  call EndSub(subname)
end program init
