*%表題   ジオポテンシャル高度の格子点データの作成
*%
*%履歴   91/05/11 佐藤正樹       C:\FORT\OORT\GRDDGP.FOR
*%       91/05/16 佐藤正樹       C:\FORT\OORT\GRDGPZ.FOR
*%
*%
*-----------------------------------------------------------------------
      PROGRAM  GRDGPZ

      INTEGER  NLAT           ! 緯度配列の次元
      PARAMETER ( NLAT = 15 )
      INTEGER  NPRS           ! 圧力配列の次元
      PARAMETER ( NPRS = 8 )

      REAL  DGP(-NLAT:NLAT, 1:NPRS)
     &                         ! ジオポテンシャル高度の平均値からの差
      REAL  GPZ(-NLAT:NLAT, 1:NPRS) ! ジオポテンシャル高度
      REAL  GPZA(1:NPRS)       ! ジオポテンシャル高度の平均値


      INTEGER  NLAT2           ! 緯度配列の次元
      PARAMETER ( NLAT2 = 18 )
      INTEGER  NPRS2           ! 圧力配列の次元
      PARAMETER ( NPRS2 = 20 )

      REAL  GPZ2(-NLAT2:NLAT2, 1:NPRS)  ! ジオポテンシャル高度
      REAL  GPZ3(-NLAT2:NLAT2, 0:NPRS2) ! ジオポテンシャル高度

      INTEGER  ILAT           ! 緯度配列の次元, DOループ
      INTEGER  IPRS           ! 圧力配列の次元, DOループ

      INTEGER  ILAT2          ! 緯度配列の次元, DOループ



*     < 入出力パラメーター >

      INTEGER   IR   / 20 /          ! 入力ユニット番号
      CHARACTER HIFL*16              ! 入力ファイル名

      INTEGER   OR   / 30 /          ! 出力ユニット番号
      CHARACTER HOFL*16              ! 出力ファイル名

      CHARACTER  HNFORM*16           ! 文字型変数の入力形式
      CHARACTER  HXFORM*16           ! 実数型変数の入力形式
      CHARACTER  HXNAME*16           ! 変数名

*     < GTOOL >

      INTEGER  NCC
      INTEGER  NDC
      PARAMETER  ( NCC=16, NDC=64 )
*
      CHARACTER  HHEAD ( NDC )*(NCC)
*
      CHARACTER  HDSET *(NCC)
      CHARACTER  HITEM *(NCC)
      CHARACTER  HTITL *(NCC*2)
      CHARACTER  HUNIT *(NCC)
      CHARACTER  HUTIM *(NCC)
      CHARACTER  HAITM1*(NCC)
      CHARACTER  HAITM2*(NCC)
      CHARACTER  HAITM3*(NCC)
      CHARACTER  HEDIT*(NCC)
      CHARACTER  HETTL*(NCC)

      INTEGER    IFNUM
      INTEGER    ITDUR
      REAL       DMIN
      REAL       DMAX
      REAL       DIVS
      REAL       DIVL
      INTEGER    ISTYP

      INTEGER    IYEAR
      INTEGER    IMONTH
      INTEGER    IDAY
      INTEGER    IHOUR
      INTEGER    IMIN
      INTEGER    ISEC

*-----------------------------------------------------------------------

      DATA GPZA /   113.,   1457.,  3011.,  5572.,
     &             9159.,  11784., 16206.,  20632.  /

*-----------------------------------------------------------------------
*     < 出力形式 >

      HNFORM='(A16)'
      HXFORM='(5E16.6)'

*-----------------------------------------------------------------------
*     < 入力データファイル名の入力 >

      WRITE(6,*) '入力ファイル名を入力してください'
      READ(5,*) HIFL
      WRITE(6,*) '出力ファイル名を入力してください'
      READ(5,*) HOFL
      WRITE(6,*) 'HIFL=', HIFL
      WRITE(6,*) 'HOFL=', HOFL

*-----------------------------------------------------------------------
*     < データの入力 >

      CALL  READORG
     I           ( IR,       HIFL,
     O             HXNAME,   DGP,
     D             NLAT,     NPRS )

*-----------------------------------------------------------------------
*     < 平均値を加算 >

      DO 300 ILAT = -NLAT, NLAT
      DO 310 IPRS = 1, NPRS
        GPZ(ILAT, IPRS) = DGP(ILAT, IPRS) + GPZA(IPRS)
 310  CONTINUE
 300  CONTINUE

*-----------------------------------------------------------------------
*     < グリッドの値の作成 >

      GPZ(-15,1) = 3.   * GPZ(-14, 1) - 2.  * GPZ(-13, 1)
      GPZ(-15,2) = 3.   * GPZ(-14, 2) - 2.  * GPZ(-13, 2)
      DO 400 IPRS = 1, NPRS
      DO 410 ILAT = -14, 14
        GPZ2(ILAT, IPRS) = GPZ(ILAT, IPRS)
 410  CONTINUE
        GPZ2(-15, IPRS) = 0.5 * GPZ(-14, IPRS) + 0.5 * GPZ(-15, IPRS)
        GPZ2(15, IPRS)  = 0.5 * GPZ(14, IPRS)  + 0.5 * GPZ(15, IPRS)
        GPZ2(-16, IPRS) =       GPZ(-15, IPRS)
        GPZ2(16, IPRS)  =       GPZ(15, IPRS)
        GPZ2(-17, IPRS) = 1.5 * GPZ(-15, IPRS) - 0.5 * GPZ(-14, IPRS)
        GPZ2(17, IPRS)  = 1.5 * GPZ(15, IPRS)  - 0.5 * GPZ(14, IPRS)
        GPZ2(-18, IPRS) = 2.  * GPZ(-15, IPRS) -       GPZ(-14, IPRS)
        GPZ2(18, IPRS)  = 2.  * GPZ(15, IPRS)  -       GPZ(14, IPRS)
 400  CONTINUE

      DO 500 ILAT2 = -18, 18
        GPZ3(ILAT2, 0) = -999.
        GPZ3(ILAT2, 1) =        GPZ2(ILAT2,8)
        GPZ3(ILAT2, 2) =        GPZ2(ILAT2,7)
        GPZ3(ILAT2, 3) = 0.5  * GPZ2(ILAT2,7) + 0.5  * GPZ2(ILAT2,6)
        GPZ3(ILAT2, 4) =        GPZ2(ILAT2,6)
        GPZ3(ILAT2, 5) = 0.5  * GPZ2(ILAT2,6) + 0.5  * GPZ2(ILAT2,5)
        GPZ3(ILAT2, 6) =        GPZ2(ILAT2,5)
        GPZ3(ILAT2, 7) = 0.75 * GPZ2(ILAT2,5) + 0.25 * GPZ2(ILAT2,4)
        GPZ3(ILAT2, 8) = 0.5  * GPZ2(ILAT2,5) + 0.5  * GPZ2(ILAT2,4)
        GPZ3(ILAT2, 9) = 0.25 * GPZ2(ILAT2,5) + 0.75 * GPZ2(ILAT2,4)
        GPZ3(ILAT2,10) =        GPZ2(ILAT2,4)
        GPZ3(ILAT2,11) = 0.75 * GPZ2(ILAT2,4) + 0.25 * GPZ2(ILAT2,3)
        GPZ3(ILAT2,12) = 0.5  * GPZ2(ILAT2,4) + 0.5  * GPZ2(ILAT2,3)
        GPZ3(ILAT2,13) = 0.25 * GPZ2(ILAT2,4) + 0.75 * GPZ2(ILAT2,3)
        GPZ3(ILAT2,14) =        GPZ2(ILAT2,3)
        GPZ3(ILAT2,15) = 2./3.* GPZ2(ILAT2,3) + 1./3.* GPZ2(ILAT2,2)
        GPZ3(ILAT2,16) = 1./3.* GPZ2(ILAT2,3) + 2./3.* GPZ2(ILAT2,2)
        GPZ3(ILAT2,17) =        GPZ2(ILAT2,2)
        GPZ3(ILAT2,18) = 2./3.* GPZ2(ILAT2,2) + 1./3.* GPZ2(ILAT2,1)
        GPZ3(ILAT2,19) = 1./3.* GPZ2(ILAT2,2) + 2./3.* GPZ2(ILAT2,1)
        GPZ3(ILAT2,20) =        GPZ2(ILAT2,1)
 500  CONTINUE

*-----------------------------------------------------------------------
*     < パラメーターの設定 >

      HDSET  = 'Oort and Peixoto'
      HITEM  = 'GPZ'
      HTITL  = 'Geopotential Height'
      HUNIT  = 'gpm'
      IFNUM  = 1
      HUTIM  = 'SEC'
      ITDUR  = 31536000.
      HAITM1 = 'GLON1'
      HAITM2 = 'GLAT37'
      HAITM3 = 'GPRS21'
      DMIN   = 0.
      DMAX   = 20000.
      DIVS   = 2000.
      DIVL   = 4000.
      ISTYP  =   1

      HEDIT  = 'YR'
      HETTL  = 'YEAR 63-73'

      IYEAR  = 1963
      IMONTH = 0
      IDAY   = 0
      IHOUR  = 0
      IMIN   = 0
      ISEC   = 0

*-----------------------------------------------------------------------
*     < GTOOL 形式で出力 >

      CALL GTDATA
     I           ( OR,     HOFL,     HXFORM, 
     I             GPZ3,
     I             HDSET , HITEM , IFNUM ,
     I             HTITL , HUNIT , HUTIM , ITDUR ,
     I             HEDIT,  HETTL, 
     I             HAITM1, HAITM2, HAITM3,
     I             DMIN  , DMAX  , DIVS  , DIVL  , ISTYP,
     I             IYEAR , IMONTH, IDAY  , IHOUR , IMIN  , ISEC ,
     D             1,      2*NLAT2+1, NPRS2+1 )

*-----------------------------------------------------------------------
      STOP
      END
