*%表題   比湿の格子点データの作成
*%
*%履歴   91/04/15 佐藤正樹       C:\FORT\OORT\GRDQHU.FOR
*%       91/05/16 佐藤正樹       C:\FORT\OORT\GRDQHU.FOR
*%
*%
*-----------------------------------------------------------------------
      PROGRAM  GRDQHU

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

      REAL  QHU(-NLAT:NLAT, 1:NPRS) ! 比湿


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

      REAL  QHU2(-NLAT2:NLAT2, 1:NPRS)  ! 温度
      REAL  QHU3(-NLAT2:NLAT2, 0:NPRS2) ! 温度

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

      INTEGER  ILAT2          ! 緯度配列の次元, DOループ
      INTEGER  IPRS2          ! 圧力配列の次元, 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

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

      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,   QHU,
     D             NLAT,     NPRS )

*-----------------------------------------------------------------------
*     < 比湿の単位の変更 : g/kg => kg/kg >

      DO 300 ILAT = -NLAT, NLAT
      DO 310 IPRS = 1, NPRS
        QHU(ILAT, IPRS) = QHU(ILAT, IPRS) / 1000.
 310  CONTINUE
 300  CONTINUE

*-----------------------------------------------------------------------
*     < グリッドの比湿の作成 >

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

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

      DO 600 IPRS2 = 0, NPRS2
      DO 610 ILAT2 = -NLAT2, NLAT2
        IF( QHU3(ILAT2, IPRS2) .LT. 0. ) THEN
          QHU3(ILAT2, IPRS2) = 0.
        ENDIF
 610  CONTINUE
 600  CONTINUE

*-----------------------------------------------------------------------
*     < テスト出力 >

C      WRITE(6,*) 'QHU3(ILAT2,20)',
C     &           (QHU3(ILAT2,20), ILAT2 = -NLAT2, NLAT2)
C      WRITE(6,*) 'QHU3', QHU3

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

      HDSET  = 'Oort and Peixoto'
      HITEM  = 'QHU'
      HTITL  = 'Specific Humidity'
      HUNIT  = 'kg/kg'
      IFNUM  = 1
      HUTIM  = 'SEC'
      ITDUR  = 31536000.
      HAITM1 = 'GLON1'
      HAITM2 = 'GLAT37'
      HAITM3 = 'GPRS21'
      DMIN   = 0.
      DMAX   = 0.02
      DIVS   =  0.001
      DIVL   =  0.005
      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             QHU3,
     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
