* PACKAGE XMKINI   !" 初期値ファイルの作成(球面非発散順圧モデル用)
*
*" [HIS]  92/02/22 
*         96/11/15 (takepiro)
*******************************************************************
      PROGRAM XMKINI
*
#if   SYS_IBMS
      INCLUDE   (ZCDIM)
      INCLUDE   (ZCCOM)
#else
#include        "zcdim.F"
#include        "zccom.F"
#endif
*
*"     ========== 変数 =========================================
*
*"     ------- 格子点データ ------------------------------------
*
      REAL       GAU   ( IDIM, JDIM, KMAX )  !" 西風  ｕ(t)
      REAL       GAV   ( IDIM, JDIM, KMAX )  !" 南風  ｖ(t)
      REAL       GAVOR ( IDIM, JDIM, KMAX )  !" 渦度  ζ(t)

      REAL       GBU   ( IDIM, JDIM, KMAX )  !" 西風  ｕ(t-Δt)
      REAL       GBV   ( IDIM, JDIM, KMAX )  !" 南風  ｖ(t-Δt)
      REAL       GBVOR ( IDIM, JDIM, KMAX )  !" 渦度  ζ(t-Δt)
*
*"     ========= 定数 ==========================================
*
*"     ------ 座標格子定数 -------------------------------------
*
      REAL       ALON  ( IDIM )              !" 経度
      REAL       DLON  ( IDIM )              !" 経度荷重
      REAL       ALAT  ( JDIM )              !" 緯度
      REAL       DLAT  ( JDIM )              !" 緯度荷重
*
*"     ==========  内部変数 ====================================
*
*"     ------- スペクトルデータ --------------------------------
*
      REAL       WDVOR ( NMDIM , KMAX     )  !" 渦度  ζ
*
*"     ======= 内部定数(DSETCで確定) ============================
*
*"     ------ 球関数展開定数 -----------------------------------
*
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX ) !" スペクトルの添字順番
      REAL       FLAPLA( NMDIM  )            !" ラプラシアンの係数
      REAL       EDEL  ( NMDIM  )            !" ζ，Ｄ→Ｕ，Ｖ
*
*"     ------ 力学定数 -----------------------------------------
*
      REAL       UVFACT ( IDIM, JDIM )       !" u→U のファクター
*
*"     ------- 時刻変数 ----------------------------------------
*
      INTEGER    ISTEP
      INTEGER    ITA, ITB
      INTEGER    IDATEA( 6 ) , IDATEB( 6 )
*
*"     ------ 実験管理パラメーター -----------------------------
*
      REAL       VORC                            !" 渦の強度
      REAL       ALAT0, ALON0, RAD               !" 擾乱の位置と大きさ
*
      DATA       ISTEP        / 0  /
      DATA       ITA, ITB     / 0 , 0 /
      DATA       IDATEA       / 1990, 1, 1, 0, 0, 0 /
      DATA       IDATEB       / 1990, 1, 1, 0, 0, 0 /
*
      INTEGER    K
*
      CALL YPREP
*
      WRITE( 6, * ) ' STRENGTH OF VORTEX ?'
      READ ( 5, * ) VORC
      WRITE( 6, * ) ' CENTER OF VORTEX IN DEGREE (lat,lon) ? ' 
      READ ( 5, * ) ALAT0, ALON0
      WRITE( 6, * ) ' RADIUS OF VORTEX IN DEGREE '
      READ ( 5, * ) RAD
*
*"          < 0. 定数のセット >
*
      CALL PCONST
*
      CALL ACRGET
     O         ( ALON  , DLON  ,
     O           ALAT  , DLAT    )
*
      CALL DCSPHE     !"  力学パートのスペクトル定数セット
     O     ( NMO   , FLAPLA , EDEL  , UVFACT, 
     C       ALAT  , DLAT                     )
*
      CALL INIVOR
     O         ( GAVOR,
     I           VORC , ALAT0 , ALON0 , RAD  ,
     C           ALAT  , ALON                )
*
      CALL GD2WD
     O         ( WDVOR ,
     I           GAVOR   )
*
      DO 5100 K = 1, KMAX
         WDVOR( NMO(1,0,0),K ) = 0.
 5100 CONTINUE
*
      CALL GENGD
     O         ( GAU   , GAV   ,
     O           GAVOR ,
     I           WDVOR ,
     C           EDEL  , UVFACT          )
*
      CALL COPY     ( GBU  , GAU  , IDIM*JDIM*KMAX )
      CALL COPY     ( GBV  , GAV  , IDIM*JDIM*KMAX )
      CALL COPY     ( GBVOR, GAVOR, IDIM*JDIM*KMAX )
*
      CALL       IHSTRT
*
      CALL IWRSTR
     I     ( GAU   , GAV   , GAVOR ,
     I       GBU   , GBV   , GBVOR ,
     I       ITA   , ITB   , IDATEA, IDATEB, ISTEP )
*
      CALL       IHEND
*
      STOP
      END
**********************************************************************
      SUBROUTINE INIVOR
     O          ( GAVOR ,
     I            VORC  , ALAT0 , ALON0 , RAD  ,
     C            ALAT  , ALON                   )
*
*"           < 初期渦 強度 = VORC ; 中心 = (RIC,RJC) ; 半径 = RAD
*
#if   SYS_IBMS
      INCLUDE   (ZCDIM)
#else
#include        "zcdim.F"
#endif
*
* [OUTPUT]
      REAL GAVOR ( IDIM , JDIM, KMAX )
*
* [INPUT]
      REAL       ALON  ( IDIM )              !" 経度
      REAL       ALAT  ( JDIM )              !" 緯度
*
      REAL       VORC                            !" 渦の強度
      REAL       ALAT0, ALON0, RAD               !" 擾乱の位置と大きさ
*
*  [INTERNAL WORK]
      INTEGER    I, J, K                         !" ループカウンター
      REAL       PI                              !" 円周率
*
      PI = ATAN( 1. )*4
      ALAT0 = ALAT0*PI/180.
      ALON0 = ALON0*PI/180.
      RAD   = RAD*PI/180.
*
      CALL RESET ( GAVOR , IDIM*JDIM*KMAX )
*
      DO 100 J = 1 , JDIM
        DO 100 I = 1, IDIM
           DO 100 K = 1, KMAX
              GAVOR ( I,J,K ) = VORC *
     &              ( ALON(I) - ALON0 )/RAD * 
     &              EXP( -(   ( ALAT(J) - ALAT0 )**2
     &                    + ( ALON(I) - ALON0 )**2  )/RAD**2 )
  100 CONTINUE
*
      RETURN
      END

