*"  96/09/24
*"  gtonpv を改良
*********************************************************************
      PROGRAM GTONTH2
*
#ifdef SYS_IBMS
      INCLUDE    (GTSINC)
      INCLUDE    (GZSIZE)
#else
#include         "gtsinc.F"
#include         "gzsize.F"
#endif
      PARAMETER  ( NFILE=9, NIDX=24 )
      COMMON     /GMWORK/ MWORK
      REAL       MWORK  ( IJKDIM )
*"[入出力データの領域確保]
      CHARACTER  HTHETA   ( NDC )*(NCC)
      CHARACTER  HHEAD ( NDC )*(NCC)
      REAL       GTHETA   ( IJKDIM )
      REAL       GDATA ( IJKDIM )
*"[途中の計算に必要なワーク領域]
      CHARACTER  HWORK ( NDC )*(NCC)
      REAL       GWORK ( IJKDIM )
*"[データに初期値を与える]
      DATA       ITHETA / 50 /
      INTEGER    IFILE( NFILE )
      DATA       IFILE / 51,52,53,54,55,56,57,58,59 /
      INTEGER    JFILE( NFILE )
      DATA       JFILE / 61,62,63,64,65,66,67,68,69 /
      DATA       JSIGTH  / 70 /
*
      CHARACTER  THT    *(NFILN)
      DATA       THT    / '$GTTMPDIR/THETA' /
      PARAMETER  (NFM1=NFILE-1)
      CHARACTER  HFILE(NFILE)*(NFILN)
      DATA       HFILE / '$GTTMPDIR/gtool.out',NFM1*' ' /
      CHARACTER  HFILEW(NFILE)*(NFILN)
      CHARACTER  SIGTH  *(NFILN)
      DATA       SIGTH  / '$GTTMPDIR/SIGTH' /
*
      CHARACTER  ITEM   *(NCC)
      DATA       ITEM   / ' ' /
      CHARACTER  UNIT   *(NCC)
      DATA       UNIT   / ' ' /
      CHARACTER  TITLE  *(NCC*2)
      DATA       TITLE  / ' ' /
      CHARACTER  DSET   *(NCC)
      DATA       DSET   / ' ' /
      CHARACTER  EDIT   *(NCC)
      DATA       EDIT   / ' ' /
      CHARACTER  ETTL   *(NCC)
      DATA       ETTL   / ' ' /
      LOGICAL    APND
      DATA       APND   / .FALSE. /
      LOGICAL    GRESET
      DATA       GRESET / .TRUE. /
      LOGICAL    HELP
      DATA       HELP   / .FALSE. /
*
      CHARACTER  HAXTH*(NCC)
      INTEGER    KTHETA
      DATA       KTHETA / 4 /
      REAL       VALTH(20)
      DATA      (VALTH(I),I=1,4) / 350., 330., 310., 290. /
      LOGICAL    OSIG
      DATA       OSIG  / .FALSE. /
*
      EXTERNAL   GSPLIN1, GSPLIN2
*
      NAMELIST  /OPTION/
     &     THT, HFILE,              !" 入力ファイル名
     &     SIGTH,                  !" 出力ファイル名
     &     ITEM, UNIT, TITLE, DSET, EDIT, ETTL, 
     &     APND ,GRESET, HELP,
     &     KTHETA, VALTH,             !" THETA 軸
     &     OSIG                    !" データ出力の有無
*
      CALL OPTARG ( 91, 'OPTION', 'HFILE', NOPT, NFILES )
      READ (91,OPTION,IOSTAT=IOS)
      CLOSE(91)
      IF ( IOS.NE.0 .OR. HELP ) THEN
         WRITE(6,OPTION)
         STOP
      ENDIF
*
      CALL GTOPEN
      CALL GTSIZE ( HTHETA   , IJKDIM )
      CALL GTSIZE ( HHEAD , IJKDIM )
      CALL GTSIZE ( HWORK , IJKDIM )
      CALL GMSIZE ( IJKDIM  )
*
      CALL GFROPN ( ITHETA   , THT )
      IF ( OSIG ) CALL GFWOPN ( JSIGTH  , SIGTH  )
      NFL = 1
      DO 300 IFL = 1, NFILES
         IF ( HFILE(IFL) .NE. ' ' ) THEN
            CALL GFROPN ( IFILE(IFL), HFILE(IFL) )
            LH=LENC(HFILE(IFL))
            HFILEW(IFL)=HFILE(IFL)(1:LENC(HFILE(IFL)))//'TH'
            WRITE(6,*) 'HFILEW=',HFILEW(IFL)
            CALL GFWOPN ( JFILE(IFL), HFILEW(IFL) )
            NFL = IFL
         ENDIF
 300  CONTINUE
*
      CALL GHPGET ( HTHETA , 'MISS' , VMISS )
*
      CALL MKAX                       !" θ 軸をつくる
     &     ( VALTH , KTHETA , VMISS , HAXTH )
*
      IAX=3                    !" 3 番目の軸を補間
      CALL SETPAR1 ( IAX , VALTH , KTHETA , HAXTH )
      CALL SETPAR2 ( IAX )
*
 1100 CONTINUE
        CALL   GFREAD                   !" THETA データ入力
     O        ( HTHETA    , GTHETA   , IEOTHETA  ,
     I          ITHETA    , 1               )
*
        IF (IEOTHETA.EQ.0) THEN
*
          CALL GMCAL1           !" σ on THETA を求める
     I         ( GSPLIN1,
     M           HTHETA    , GTHETA   ,
     I           '   '  , '   '  )
          CALL GHCSET ( HTHETA   , 'ITEM' , 'SIG' )
          CALL GHCSET ( HTHETA   , 'UNIT' , ' ' )
          CALL GHCSTS ( HTHETA   , 'TITL', 'sigma on THETA' )
*
          IF ( OSIG ) THEN
            CALL  GFWRIT        !" σ on THETA を出力
     I             ( HTHETA  , GTHETA   ,
     I             JSIGTH , 1     , 0       )
          ENDIF
*
          IJ  = 0
 1200       CONTINUE
            IJ = IJ + 1
            IF ( IJ.LE.NFL ) THEN
              IFILEX = IFILE( IJ )
              JFILEX = JFILE( IJ )
              CALL GFREAD       !" DATA on σ の入力
     O            ( HHEAD , GDATA , IEOD  ,
     I              IFILEX, 1               )
*
              CALL GPFSET       !" HWORK に σ on THETA を代入
     I           ( HTHETA  , GTHETA   ,
     I           '   '  , '  '  ,
     O            HWORK, GWORK  )
*
              NDIF=0            !" 補間
              CALL SETDIF ( NDIF )
              CALL GMCAL2       !" DATA on THETA を求める
     I             ( GSPLIN2,
     M             HWORK, GWORK,
     I             HHEAD, GDATA,
     I             '   ', '   ' )
*
              IF ( ITEM .NE. ' ' ) THEN
                CALL GHCSET ( HWORK, 'ITEM', ITEM ) 
              ELSE
                CALL GHCSET ( HWORK, 'ITEM' ,
     +                 HFILEW(IJ)(1:LENC(HFILEW(IJ)))  )
              ENDIF
              IF ( UNIT .NE. ' ' ) THEN
                CALL GHCSET ( HWORK, 'UNIT' , UNIT )
              ELSE
                CALL GHCGET ( HHEAD, 'UNIT' , UNIT )
                CALL GHCSET ( HWORK, 'UNIT' , UNIT )
              ENDIF
              IF ( TITLE .NE. ' ' ) THEN 
                CALL GHCSTS ( HWORK, 'TITL' , TITLE )
              ELSE
                CALL GHCSTS ( HWORK, 'TITL', 
     +                 HFILE(IJ)(1:LENC(HFILE(IJ)))//' on THETA' )
              ENDIF
              IF ( DSET .NE. ' ' ) THEN
                CALL GHCSET ( HWORK, 'DSET', DSET )
              ELSE
                CALL GHCGET ( HHEAD, 'DSET' , DSET )
                CALL GHCSET ( HWORK, 'DSET' , DSET )
              ENDIF
              IF ( GRESET ) THEN
                CALL GHRSGP ( HWORK )
              ENDIF
*
              CALL  GFWRIT      !" DATA on THETA を出力
     I             ( HWORK , GWORK ,
     I               JFILEX, 1     , 0       )
*
              GOTO 1200
            ENDIF
*
          GOTO 1100
        ENDIF
* commented out by toyoda
*     CALL GTCLSE
*
      STOP
      END


*****************************************************
*" θ 軸を作成する
******************************************************
      SUBROUTINE MKAX
     I      ( VAL, LEVEL , VMISS, HAITM )
*
#ifdef SYS_IBMS
      INCLUDE    (GZSIZE)                !" NCC, NDC, NED
      INCLUDE    (GZIWK)                 !" NW, NAXIS
#else
#include         "gzsize.F"              !" NCC, NDC, NED
#include         "gziwrk.F"              !" NW, NAXIS
#endif
*
      INTEGER    LEVEL                  !" 与える θ 軸の格子数
      REAL       VAL(20)               !" 与える θ の値
*
      CHARACTER  HAITM   *(NCC)          !" 軸の名称 (例:CTHETA10 )
      REAL       AX  ( NAXIS )           !" 格子位置
      REAL       DAX ( NAXIS )           !" 軸の重み
      CHARACTER  HAX (NDC)*(NCC)         !" 軸ファイルのヘッダ
*
      HAITM = 'CTHETA'
      LH = LENC(HAITM)
      IF(LEVEL.LT.10) THEN
        WRITE( HAITM(LH+1:LH+1),'(I1)') LEVEL
      ELSE IF(LEVEL.LT.100) THEN
        WRITE( HAITM(LH+1:LH+2),'(I2)') LEVEL
      ELSE
        CALL MSGDMP ( 'E', 'MKAX', 'GRID NUMBER IS TOO MUCH!!')
      ENDIF
*
      DO 10 I = 1, LEVEL
        AX(I) = VAL(LEVEL+1-I)      !" 逆の順番で代入
        DAX(I)=1.                      !" 注意！この重みはいいかげん
   10 CONTINUE
*
      CALL    GHPACA                   !" 格子情報を与える
     O         ( HAX,
     I           'IAXLOC', HAITM,
     I           'theta', ' ',
     I           LEVEL,
     I           'UR4' , VMISS ,
     I           VMISS , VMISS , VMISS , VMISS , 1 )
*
      CALL    GUSIAX                   !" 内部設定格子の設定
     I         ( HAX , AX , DAX ,
     I           1                )
*
      CALL  GUWIAX                     !" 内部設定格子の出力
     I         ( HAX )
*
      RETURN
      END





********************************************************************
*" (旧 95/01/16 horinout  CUBIC SPLINE INTERPOLATIO)
*" 改良している
*" 補間か微分かを決める NDIF は ENTRY 文で代入
*===================================================================
      SUBROUTINE GSPLIN2
     I         ( HHEAD1, GDATA1,
     I           HHEAD2, GDATA2,
     O           HHEADO, GDATAO,
     D           IDIM1 , JDIM1 , KDIM1,
     D           IDIM2 , JDIM2 , KDIM2 )
*
#ifdef SYS_IBMS
      INCLUDE    (GZSIZE)
#else
#include         "gzsize.F"
#endif
*
      CHARACTER  HHEAD1 ( * )*(*)
      REAL       GDATA1 ( IDIM1,JDIM1,KDIM1 )
      CHARACTER  HHEAD2 ( * )*(*)
      REAL       GDATA2 ( IDIM2,JDIM2,KDIM2 )
      CHARACTER  HHEADO ( * )*(*)
      REAL       GDATAO ( * )
* internal work
      CHARACTER  HAX1 ( NDC )*(NCC)
      CHARACTER  HAX2 ( NDC )*(NCC)
      PARAMETER  (NAXD=2048)
      REAL       AX2(NAXD)
      REAL       X(NAXD),Y(NAXD),Y2(NAXD)
*
      INTEGER   IAX
      SAVE  IAX, NDIF
*
* / read axis /
*
      CALL GTSIZE(HAX1,NAXD)
      CALL GTSIZE(HAX2,NAXD)

      CALL GUQAXV   !" original axis
     I       ( HHEAD2, IAX, 'LOC' ,
     O         HAX2, AX2 , IEOD2 )      
*
*  / header /
*
      CALL GHPGET( HHEAD1,'MISS',VMISS )
*
      IJKDO=IDIM1*JDIM1*KDIM1
      CALL GTSIZE( HHEADO,IJKDO )
*
* / cubic spline interpolation /
*
        DO 1300 J=1,JDIM2
        DO 1300 I=1,IDIM2
          N=0
          DO 1310 K=1,KDIM2
            IF(GDATA2(I,J,K).NE.VMISS) THEN
              N=N+1
              X(N)=AX2(K)
              Y(N)=GDATA2(I,J,K)
            ENDIF
 1310     CONTINUE
          CALL SPLINE(X,Y,N,1E31,1E31,Y2)
          DO 1320 K=1,KDIM1
            IJK=I+IDIM2*(J-1)+IDIM2*JDIM2*(K-1)
            AX1=GDATA1(I,J,K)
            IF (AX1.NE.VMISS) THEN
              CALL SPLINT(X,Y,Y2,N,AX1,NDIF,YB)
              GDATAO(IJK)=YB
            ELSE
              GDATAO(IJK)=VMISS
            ENDIF
 1320     CONTINUE
 1300   CONTINUE
*
      RETURN
*===========
      ENTRY SETPAR2(IAXE)
      IAX  = IAXE
      RETURN
*===========
      ENTRY SETDIF (NDIFE)
      NDIF=NDIFE
      RETURN
      END


********************************************************************
*" (旧: 95/01/16 horinout  CUBIC SPLINE INTERPOLATION)
*" 適当な θ 面における σ の値を求める
*" θ 軸の情報は ENTRY 文から入手
*===================================================================
      SUBROUTINE GSPLIN1
     I         ( HDT   , DT    ,
     O           HPW   , PW    ,
     D           IDIM  , JDIM  , KDIM  )
*
#ifdef SYS_IBMS
      INCLUDE    (GZSIZE)
#else
#include         "gzsize.F"
#endif
*
      CHARACTER  HDT ( * )*(*)
      REAL       DT  ( IDIM, JDIM, KDIM )
      CHARACTER  HPW ( * )*(*)
      REAL       PW  ( * )
* internal work
      CHARACTER  HAITM*(NCC)
      CHARACTER  HAITM2*(NCC)
      CHARACTER  HAX ( NDC )*(NCC)
      PARAMETER  (NAXD=2048)
      REAL       AXA(NAXD),AXB(NAXD)
      LOGICAL    OCYCL
      REAL       X(NAXD),Y(NAXD),Y2(NAXD)
*
      SAVE  IAX
*
      REAL  SIGLVL(557056)   !" include/gtsinc.F の IJKDIM の大きさ
      CHARACTER  HAITMX*(NCC), HAXTH*(NCC)
      REAL  VALPV(20), VALPVX(20)   !" 求める PV のレベル
      SAVE  KPV,VALPV,HAXTH

      LOGICAL OFIRST
      DATA    OFIRST / .TRUE. /
      SAVE    OFIRST
*     
      LMAX=30                !" 鉛直にアバウトに攻めていく回数
      EPS=0.0001             !" PV の収束判定値(PVU)
      NDIF=0                 !" 補間
*
* / read axis /
*
      CALL GTSIZE(HAX,NAXD)
      CALL GUQAXV   !" original axis
     I       ( HDT, IAX, 'LOC' ,
     O         HAX, AXA , IEOD )
*
      CALL GHCGET(HDT,'AITM3',HAITM) !" original の軸の情報を得る
      CALL GHCGET(HDT,'AITM2',HAITM2) !" original の軸の情報を得る
      CALL GUQAXS
     I       ( HAITM, 'LOC',
     O         HAX  , AXB  , IEOD )
      CALL GUQAXD
     I       ( HAITM,
     O         NADIM, OCYCL )
*
*  / header /
*
      CALL GHPGET( HDT,'MISS',VMISS )
*
      CALL GHCOPY(HDT,HPW)           !" PV 軸の情報を与える
        CALL GHCSET(HPW,'AITM3',HAXTH)
        CALL GHPSET(HPW,'AEND3',KPV)
        IJKDO=IDIM*JDIM*KPV
      CALL GTSIZE( HPW,IJKDO )
*
* / cubic spline interpolation /
*
      SIGMAX=AXB(1)
      SIGMIN=AXB(NADIM)

      IF(SIGMIN.GT.SIGMAX) THEN   !" SIGMA の与え方の妥当性チェック
        CALL MSGDMP( 'E','GSPLIN','SIGMIN IS GREATER THAN SIGMAX')
      ENDIF

      DO 1300 J=1,JDIM
      DO 1300 I=1,IDIM
        N=0
        DO 1310 K=1,KDIM
          IF(DT(I,J,K).NE.VMISS) THEN
            N=N+1
            X(N)=AXA(K)
            Y(N)=DT(I,J,K)
          ENDIF
 1310   CONTINUE
        CALL SPLINE(X,Y,N,1E31,1E31,Y2)

        AXB2=0.
        LMIN=1
        DO 1320 K=1,KPV
          DO 4000 L=LMIN,LMAX
            IJK=I+IDIM*(J-1)+IDIM*JDIM*(K-1)
            AXB1=(SIGMAX-SIGMIN)*(L-1)/(LMAX-1)+SIGMIN
            CALL SPLINT(X,Y,Y2,N,AXB1,NDIF,YB)

            IF((YB-VALPV(K)).LT.0) THEN   !" 目的のレベルを越えたら
              IF (L.EQ.1) THEN    !" VALPV の与え方チェック
                CALL MSGDMP('W','GSPLIN','VALPV IS TOO HIGH')
                WRITE(6,*) '(I,J)=',I,J
                SIGLVL(IJK)=VMISS
                GOTO 1320
              ENDIF
 3000         CONTINUE                         !" 2 分法で正確に求める
              AXBMID = (AXB1+AXB2)/2
                CALL SPLINT(X,Y,Y2,N,AXBMID,NDIF,YB)
                IF((YB-VALPV(K)).LT.0.) THEN
                  AXB1=AXBMID
                ELSE
                  AXB2=AXBMID
                ENDIF
                IF(ABS(YB-VALPV(K)).GE.EPS) THEN !" 収束判定
                  GOTO 3000
                ENDIF
              SIGLVL(IJK)=AXBMID               !" SIGMA LEVEL を代入
              LMIN=L-2       !" 次のレベルはそこより下から探す。
              GOTO 1320
            ENDIF

            AXB2=AXB1
 4000     CONTINUE
          SIGLVL(IJK)=VMISS !" SIGMA LEVEL が見つからないときは欠損値を代入
          LMIN=LMAX         !" 次は一番下のレベルから探す
 1320   CONTINUE
 1300 CONTINUE

      DO 6100 I=1,IDIM
        DO 6100 J=1,JDIM
          DO 6100 K=1,KPV
            IJK=I+IDIM*(J-1)+IDIM*JDIM*(K-1)
            IJKR=I+IDIM*(J-1)+IDIM*JDIM*(KPV-K)
            PW(IJK)=SIGLVL(IJKR)  !" 鉛直の順番並べ替え
 6100 CONTINUE
*
      RETURN
*===========
      ENTRY SETPAR1(IAXE,VALPVX,KPVX,HAITMX)
      IAX  = IAXE
      KPV=KPVX
      HAXTH=HAITMX
      DO 5000 I=1,KPV
        VALPV(I)=VALPVX(I)
 5000 CONTINUE
      RETURN
      END


************************************************************************
*(Original Source Location: /home/horinout/src/numrecipes/splin.f)
*
* NUMERICAL  RECIPES
*
*" 3.3  3次スプライン補間
*" subroutines  SPLINE  準備 ( Y2 = y" )
*"              SPLINT  補間  仮定: X(1)<X(2)<..<X(N) or X(1)>X(2)>..>X(N)
*" 95/02/01  Ver.3 一、二階微分値の計算ができる(NDIF で指定)。
*" [注意]  SPLINT の引数は Ver.1, Ver.2 と違う. (入力 NDIF を追加)
************************************************************************
      SUBROUTINE SPLINE
     I            ( X,Y,N,YP1,YPN,
     O              Y2 )
      IMPLICIT  REAL  (A-H,O-Z)
      PARAMETER (NMAX=200)
      DIMENSION X(N),Y(N), Y2(N), U(NMAX)
*
      IF(N.GT.NMAX) THEN
        WRITE(6,*) '***(SPLINE)*** NOT ENOUGH WORKING AREA SIZE NMAX'
        STOP
      ENDIF
*
      IF (YP1 .GT. .99E30) THEN !" the lower b.c. is set either to be "natural"
        Y2(1)=0.
        U(1)=0.
      ELSE
        Y2(1)=-0.5
        U(1)=(3./(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1)
      ENDIF
      DO 11 I=2,N-1
        SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1))
        P=SIG*Y2(I-1)+2.
        Y2(I)=(SIG-1.)/P
        U(I)=(6.*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1))
     &       /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P
   11 CONTINUE
      IF (YPN .GT. .99E30) THEN !" the upper b.c. is set either to be "natural"
        QN=0.
        UN=0.
      ELSE
        QN=0.5
        UN=(3./(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1)))
      ENDIF
      Y2(N)=(UN-QN*U(N-1))/(QN*Y2(N-1)+1.)
      DO 12 K=N-1,1,-1
        Y2(K)=Y2(K)*Y2(K+1)+U(K)
   12 CONTINUE
      RETURN
      END
************************************************************************
* NDIF=0:simple interpolation, NDIF=1,2:d(YA)/d(XA),d^2(YA)/d(XA)^2
*
      SUBROUTINE SPLINT
     I            ( XA,YA,Y2A,N,X, NDIF,
     O              Y )
      IMPLICIT  REAL  (A-H,O-Z)
      DIMENSION XA(N),YA(N),Y2A(N)
*
      KLO=1
      KHI=N
    1  IF (KHI-KLO.GT.1) THEN
         K=(KHI+KLO)/2
         IF((XA(N).GT.XA(1)).EQV.(XA(K).GT.X)) THEN
           KHI=K
         ELSE
           KLO=K
         ENDIF
      GOTO 1
       ENDIF
      H=XA(KHI)-XA(KLO)
      IF (H.EQ.0.) THEN
        WRITE(6,*) '***(SPLINT)*** Bad XA input.'
        STOP
      ENDIF
      A=(XA(KHI)-X)/H
      B=(X-XA(KLO))/H
      IF (NDIF.NE.1 .AND. NDIF.NE.2) THEN
        Y=A*YA(KLO)+B*YA(KHI)+
     &     ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6.
      ELSE IF (NDIF.EQ.1) THEN  !" 1st derivative
        Y=(YA(KHI)-YA(KLO))/H 
     &      - (3.*A*A-1.)/6.*H*Y2A(KLO) + (3.*B*B-1.)/6.*H*Y2A(KHI)
      ELSE IF (NDIF.EQ.2) THEN  !" 2nd derivative
        Y=A*Y2A(KLO)+B*Y2A(KHI)
      ENDIF
      RETURN
      END
