* PACKAGE XMKPHIS   !" 境界条件データ（地形)の作成(gtool3 形式)
*
*" [HIS]  96/11/15 (takepiro)
*******************************************************************
      PROGRAM MKPHIS
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)
      INCLUDE   (ZCCOM)
      INCLUDE   (GZSIZE)
#else
#include        "zcdim.F"
#include        "zccom.F"
#include        "gzsize.F"
#endif
*
*"     ========== 変数 =========================================
*
*"     ------- 格子点データ ------------------------------------
*
      REAL       GPHIS ( IDIM, JDIM )             !" 地表Φ
*
*"     ========= 定数 ==========================================
*
*"     ------ 座標格子定数 -------------------------------------
*
      REAL       ALON  ( IDIM )              !" 経度
      REAL       DLON  ( IDIM )              !" 経度荷重
      REAL       ALAT  ( JDIM )              !" 緯度
      REAL       DLAT  ( JDIM )              !" 緯度荷重
*
*"     ==========  内部変数 ====================================
*
*
*"     ------- 時刻変数 ----------------------------------------
*
      INTEGER    ISTEP
      INTEGER    IT
      INTEGER    IDATE( 6 )
*
*"     ------ 実験管理パラメーター -----------------------------
*
      REAL       BGEP                            !" 平均地表ジオポテンシャル
      REAL       DGEP                            !" 擾乱地表ジオポテンシャル
      REAL       ALAT0, ALON0, RAD               !" 擾乱の位置と大きさ
*
*"     ------ 入出力用変数 -------------------------------------
*
      DATA       ISTEP   / 0  /
      DATA       IT      / 0 /
      DATA       IDATE   / 1990, 1, 1, 0, 0, 0 /
*
      INTEGER    ITDUR, IOMODE, NOEND
      DATA       ITDUR   / 1 /
      DATA       IOMODE  / 1 /
      DATA       NOEND   / 0 /
*
      CHARACTER  HITEM*(NCC),HTITL*(NCC*2),HUNIT*(NCC)
*
      DATA       HITEM   /'GPHIS'/
      DATA       HTITL   /'surface geopotential'/
      DATA       HUNIT   /'m2/s2'/
*
      CHARACTER  FILE   *(NFILN)             !" 地表Φファイル
      CHARACTER  DFMT   *(NCC)               !" 
      DATA       FILE, DFMT    / 'shallow.phis', 'UR4' /
      NAMELIST  /NMPHIS/ FILE, DFMT
*
      INTEGER    IFPAR, JFPAR
      INTEGER    JFILE
*
*"          < 0. パラメター読み込み >
*
      CALL YPREP
*
      WRITE( 6, * ) ' STRENGTH OF MEAN SURFACE GEO-POT?'
      READ ( 5, * ) BGEP
      WRITE( 6, * ) ' STRENGTH OF DISTURBANCE SURFACE GEO-POT?'
      READ ( 5, * ) DGEP
      WRITE( 6, * ) ' CENTER OF DIST. IN DEGREE (lat,lon) ? ' 
      READ ( 5, * ) ALAT0, ALON0
      WRITE( 6, * ) ' RADIUS OF DIST. IN DEGREE '
      READ ( 5, * ) RAD
*
      CALL   REWNML ( IFPAR , JFPAR )
      READ   ( IFPAR, NMPHIS, END=1190 )
 1190 WRITE  ( JFPAR, NMPHIS )
*
*"          < 1. 定数のセット >
*
      CALL       PCONST
*
      CALL ACRGET
     O         ( ALON  , DLON  ,
     O           ALAT  , DLAT  )
*
*
*"       < 2. データ作成 >
*
      CALL CALPHS
     O         ( GPHIS,
     I           BGEP  , DGEP  , ALAT0 , ALON0 , RAD  ,
     C           ALAT  , ALON                           )
*
*"       < 3. データ書き込み >
*
      CALL IHSTRT
*
      CALL GFWOPN
     O         ( JFILE , FILE )
*
      CALL GDWRIT
     O         ( GPHIS ,
     I           HITEM , HTITL , HUNIT ,
     I           IT    , IDATE , ISTEP , ITDUR ,
     I           JFILE , IOMODE, NOEND , 'ASFC', DFMT  )
*
      STOP
      END
***********************************************************************
      SUBROUTINE CALPHS
     O         ( GPHIS ,
     I           BGEP  , DGEP  , ALAT0 , ALON0 , RAD  ,
     C           ALAT  , ALON                           )
*
*"           < 地表面ジオポテンシャル
*"                 擾乱 geo potential 強度 = DGEP 
*"                 中心 = (RIC,RJC) ; 半径 = RAD >
*
#if   SYS_IBMS
      INCLUDE   (ZCDIM)
#else
#include        "zcdim.F"
#endif
*
* [OUTPUT]
      REAL       GPHIS ( IDIM, JDIM )             !" 地表Φ
*
* [INPUT]
      REAL       BGEP                            !" 平均ジオポテンシャル
      REAL       DGEP                            !" 擾乱ジオポテンシャル
      REAL       ALAT0, ALON0, RAD               !" 擾乱の位置と大きさ
*
      REAL       ALON  ( IDIM )                  !" 経度
      REAL       ALAT  ( JDIM )                  !" 緯度
*
*  [INTERNAL WORK]
      INTEGER    I, J                            !" ループカウンター
      REAL       PI                              !" 円周率
*
*
      PI = ATAN( 1. )*4
      ALAT0 = ALAT0*PI/180.
      ALON0 = ALON0*PI/180.
      RAD   = RAD*PI/180.
*
      CALL RESET ( GPHIS , IDIM*JDIM )
*
      DO 100 J = 1, JDIM
        DO 100 I = 1, IDIM
              GPHIS ( I,J ) 
     &      =   BGEP 
     &        + DGEP * EXP( -(   ( ALAT(J) - ALAT0 )**2
     &                         + ( ALON(I) - ALON0 )**2  )/RAD**2 )
  100 CONTINUE
*
      RETURN
      END