*" σ面上で PV (ちゃんと 3 次元)を求める     96/09/02 akahori
*********************************************************************
      PROGRAM GTPV2
      IMPLICIT NONE
*
#ifdef SYS_IBMS
      INCLUDE    (GTSINC)
      INCLUDE    (GZSIZE)
#else
#include         "gtsinc.F"
#include         "gzsize.F"
#endif
      COMMON     /GMWORK/ MWORK
      REAL       MWORK  ( IJKDIM )
*
*"[入出力データの領域確保]
      CHARACTER  HTHETA ( NDC )*(NCC)  !" θ
      REAL       GTHETA ( IJKDIM )
      CHARACTER  HU ( NDC )*(NCC)      !" ｕ
      REAL       GU ( IJKDIM )
      CHARACTER  HV ( NDC )*(NCC)      !" ｖ
      REAL       GV ( IJKDIM )
      CHARACTER  HPS ( NDC )*(NCC)     !" Ps
      REAL       GPS ( IJKDIM )
      CHARACTER  HPV ( NDC )*(NCC)     !" Ps
      REAL       GPV ( IJKDIM )
*
*"[途中の計算に必要なワーク領域]
      CHARACTER  HUS ( NDC )*(NCC)     !" ∂ｕ∂σ
      REAL       GUS ( IJKDIM )
      CHARACTER  HVS ( NDC )*(NCC)     !" ∂ｖ/∂σ
      REAL       GVS ( IJKDIM )
      CHARACTER  HVOR ( NDC )*(NCC)    !" ∂ζ/∂σ
      REAL       GVOR ( IJKDIM )
      CHARACTER  HTHETAX ( NDC )*(NCC)   !" ∂θ/∂ｘ
      REAL       GTHETAX ( IJKDIM )
      CHARACTER  HTHETAY ( NDC )*(NCC)   !" ∂θ/∂ｙ
      REAL       GTHETAY ( IJKDIM )
      CHARACTER  HTHETAS ( NDC )*(NCC)   !" ∂θ/∂σ
      REAL       GTHETAS ( IJKDIM )
*
*"[データに初期値を与える]
      INTEGER ITHETA, IU, IV, IPS, JPV
      DATA       ITHETA / 50 /
      DATA       IU     / 52 /
      DATA       IV     / 54 /
      DATA       IPS    / 56 /
      DATA       JPV    / 64 /
*
      CHARACTER  THT    *(NFILN)
      DATA       THT    / '$GTTMPDIR/THETA' /
      CHARACTER  U      *(NFILN)
      DATA       U      / '$GTTMPDIR/U' /
      CHARACTER  V      *(NFILN)
      DATA       V      / '$GTTMPDIR/V' /
      CHARACTER  PS     *(NFILN)
      DATA       PS     / '$GTTMPDIR/PS' /
      CHARACTER  PV     *(NFILN)
      DATA       PV     / '$GTTMPDIR/PV' /
*
*"[球面調和関数による微分のための定義]
      LOGICAL    DX  !" T: x 微分 (DY=Tの方が優先される)
      LOGICAL    DY  !" T: y 微分
      LOGICAL    A   !" T: 地球半径で割る
      INTEGER    WB !"(cosφ)**WB の重みを微分前に掛ける   (x微分の場合WBと
      INTEGER    WA !"(cosφ)**WA の重みを微分後に掛ける    WAは同じ働き)
      LOGICAL    ZONAL !" T: 一番目の軸を y 軸とみなす(波数０成分を仮定)
*
      CHARACTER  ITEM   *(NCC)
      DATA       ITEM   / 'PV' /
      CHARACTER  UNIT   *(NCC)
      DATA       UNIT   / 'PVU' /
      CHARACTER  TITLE  *(NCC*2)
      DATA       TITLE  / 'PV' /
      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. /
*
      EXTERNAL   GSPLIN2, CORIOL, DSPHE, U2VOR, PTLVOR
*
      INTEGER NFILE, NOPT, IOS, IAX, IEOTHT, IEOU, IEOPS, NDIF,
     &        IXDIM, IXSTR, IXEND,
     &        IYDIM, IYSTR, IYEND,
     &        IZDIM, IZSTR, IZEND, 
     &        IEOV

      NAMELIST  /OPTION/
     &     THT, U, V, PS,                    !" 入力ファイル名
     &     PV,                               !" 出力ファイル名
     &     ITEM, UNIT, TITLE, DSET, EDIT, ETTL, 
     &     APND ,GRESET, HELP
*"[オプション]
      CALL OPTARG ( 91, 'OPTION', 'HFILE', NOPT, NFILE )
      READ (91,OPTION,IOSTAT=IOS)
      CLOSE(91)
      IF ( IOS.NE.0 .OR. HELP ) THEN
         WRITE(6,OPTION)
         STOP
      ENDIF
*"[gtool をはじめるよ]
      CALL GTOPEN
*"[データ入力のための前処理]
      CALL GTSIZE ( HTHETA , IJKDIM )
      CALL GTSIZE ( HU     , IJKDIM )
      CALL GTSIZE ( HV     , IJKDIM )
      CALL GTSIZE ( HPS    , IJKDIM )
      CALL GTSIZE ( HPV    , IJKDIM )
      CALL GTSIZE ( HUS    , IJKDIM )
      CALL GTSIZE ( HVS    , IJKDIM )
      CALL GTSIZE ( HVOR   , IJKDIM )
      CALL GTSIZE ( HTHETAX, IJKDIM )
      CALL GTSIZE ( HTHETAY, IJKDIM )
      CALL GTSIZE ( HTHETAS, IJKDIM )
      CALL GMSIZE ( IJKDIM  )
*"[データ入力]
      CALL GFROPN ( ITHETA , THT )
      CALL GFROPN ( IU     , U   )
      CALL GFROPN ( IV     , V   )
      CALL GFROPN ( IPS    , PS  )
      CALL GFOOPN ( JPV    , PV  , APND )
*"[鉛直補間のためのパラメータ設定]
      IAX=3                    !" 3 番目の軸を補間
      CALL SETPAR2 ( IAX )
*
 1100 CONTINUE
*
*"[データ入力]
         CALL   GFREAD                      !" θ 入力
     O        ( HTHETA , GTHETA , IEOTHT ,
     I          ITHETA , 1                )
*
         CALL   GFREAD                      !" Ｕ 入力
     O        ( HU     , GU     , IEOU   ,
     I          IU     , 1                )
*
         CALL   GFREAD                      !" Ｖ 入力
     O        ( HV     , GV     , IEOV   ,
     I          IV     , 1                )
*
         CALL   GFREAD                      !" Ps 入力
     O        ( HPS    , GPS    , IEOPS  ,
     I          IPS    , 1               )
*
         IF (      ( IEOTHT .EQ. 0 ) .AND. ( IEOU  .EQ. 0 )
     +       .AND. ( IEOV   .EQ. 0 ) .AND. ( IEOPS .EQ. 0 ) ) THEN
*
*"[θの鉛直微分]
           CALL GPFSET
     +          ( HTHETA , GTHETA ,
     +             '   ' , '   '  ,
     +            HTHETAS, GTHETAS )
            NDIF=1                       !" 一階微分
            CALL SETDIF ( NDIF )
            CALL GMCAL1                  !" dθ/dσ を求める
     I         ( GSPLIN2,
     I           HTHETAS, GTHETAS,
     I           '   ' , '   ' )
*"[uの鉛直微分]
           CALL GPFSET
     +          ( HU     , GU     ,
     +             '   ' , '   '  ,
     +            HUS    , GUS     )
            NDIF=1                       !" 一階微分
            CALL SETDIF ( NDIF )
            CALL GMCAL1                  !" dｕ/dσ を求める
     I         ( GSPLIN2,
     I           HUS   , GUS     ,
     I           '   ' , '   ' )
*"[vの鉛直微分]
           CALL GPFSET
     +          ( HV     , GV     ,
     +             '   ' , '   '  ,
     +            HVS    , GVS     )
            NDIF=1                       !" 一階微分
            CALL SETDIF ( NDIF )
            CALL GMCAL1                  !" dｖ/dσ を求める
     I         ( GSPLIN2,
     I           HVS   , GVS     ,
     I           '   ' , '   ' )
*"[θの X 微分]
           CALL GPFSET
     +          ( HTHETA , GTHETA ,
     +           '   '   ,'   '  ,
     +            HTHETAX, GTHETAX )
           DX=.TRUE.
           DY=.FALSE.
           ZONAL=.FALSE.
           WB=0
           WA=-1
           A=.TRUE.
           CALL EDSPHE(DX,DY,WA,WB,A,ZONAL)
           CALL GMCAL1
     I         ( DSPHE  ,
     M           HTHETAX, GTHETAX,
     I           '  '   , '  '   )
*"[θの Y 微分]
           CALL GPFSET
     +          ( HTHETA , GTHETA ,
     +           '   '   ,'   '  ,
     +            HTHETAY, GTHETAY )
           DX=.FALSE.
           DY=.TRUE.
           ZONAL=.FALSE.
           WB=0
           WA=0
           A=.TRUE.
           CALL EDSPHE(DX,DY,WA,WB,A,ZONAL)
           CALL GMCAL1
     I         ( DSPHE  ,
     M           HTHETAY, GTHETAY,
     I           '  '   , '  '   )
*"[相対渦度の鉛直成分]
           CALL GPFSET
     +          ( HU     , GU    ,
     +           '   '   ,'   '  ,
     +            HVOR   , GVOR   )
           CALL GMCAL2
     I          ( U2VOR,
     M            HVOR , GVOR,    !" input U ; output VOR
     I            HV   , GV  ,
     I            '  ' , '  ' )
*"[相対渦度-->絶対渦度]
           CALL GMCAL1            !" 絶対渦度を求める
     I          ( CORIOL,
     M            HVOR  , GVOR,
     I            '   ' , '   ' )
*"[合体]
           CALL GUQSIZ
     I          ( HVOR  ,
     O            IXSTR , IXEND, IXDIM,
     O            IYSTR , IYEND, IYDIM,
     O            IZSTR , IZEND, IZDIM  )
           CALL GHCOPY ( HVOR  , HPV )
           CALL PTLVOR
     I          ( GPS    ,
     I            GUS    ,
     I            GVS    ,
     I            GVOR   ,
     I            GTHETAS,
     I            GTHETAX,
     I            GTHETAY,
     O            GPV    ,
     I            IXDIM, IYDIM, IZDIM )
*"[ PV データ出力]
                  IF ( ITEM .NE. ' ' ) THEN
                     CALL GHCSET( HPV , 'ITEM', ITEM )
                  ENDIF
                  IF ( UNIT .NE. ' ' ) THEN
                     CALL GHCSET( HPV , 'UNIT', UNIT )
                  ENDIF
                  IF ( TITLE .NE. ' ' ) THEN
                     CALL GHCSTS( HPV , 'TITL', TITLE )
                  ENDIF
                  IF ( DSET .NE. ' ' ) THEN
                     CALL GHCSET( HPV , 'DSET', DSET )
                  ENDIF
                  IF ( GRESET ) THEN
                     CALL GHRSGP( HPV  )
                  ENDIF

              IF ( GRESET ) CALL GHRSGP( HPV )
              CALL GFWRIT
     I             ( HPV , GPV   ,
     I               JPV , 1     , 0       )
*
      GOTO 1100
         ENDIF
* commented out by toyoda
*     CALL GTCLSE
      STOP
      END



********************************************************************
*" (旧 95/01/16 horinout  CUBIC SPLINE INTERPOLATIO)
*" 改良している
*" 補間か微分かを決める NDIF は ENTRY 文で代入
*===================================================================
      SUBROUTINE GSPLIN2
     I         ( HHEAD1, GDATA1,
     O           HHEADO, GDATAO,
     D           IDIM1 , JDIM1 , KDIM1 )
      IMPLICIT NONE
*
#ifdef SYS_IBMS
      INCLUDE    (GZSIZE)
#else
#include         "gzsize.F"
#endif
*
      INTEGER    IDIM1, JDIM1, KDIM1
      CHARACTER  HHEAD1 ( * )*(*)
      REAL       GDATA1 ( IDIM1,JDIM1,KDIM1 )
      CHARACTER  HHEADO ( * )*(*)
      REAL       GDATAO ( * )
* internal work
      CHARACTER  HAX1 ( NDC )*(NCC)
      INTEGER NAXD
      PARAMETER  (NAXD=2048)
      REAL       AX1(NAXD)
      REAL       X(NAXD),Y(NAXD),Y2(NAXD)
*
      INTEGER   IAX, NDIF,  IAXE, NDIFE
      SAVE  IAX, NDIF
      INTEGER I, J, IJKDO, N, K, IJK, IEOD1
      REAL VMISS, AX, YB
*
* / read axis /
*
      CALL GTSIZE(HAX1,NAXD)

      CALL GUQAXV   !" original axis
     I       ( HHEAD1, IAX, 'LOC' ,
     O         HAX1, AX1 , IEOD1 )      
*
*  / header /
*
      CALL GHPGET( HHEAD1,'MISS',VMISS )
*
      IJKDO=IDIM1*JDIM1*KDIM1
      CALL GTSIZE( HHEADO,IJKDO )
*
* / cubic spline interpolation /
*
        DO 1300 J=1,JDIM1
        DO 1300 I=1,IDIM1
          N=0
          DO 1310 K=1,KDIM1
            IF(GDATA1(I,J,K).NE.VMISS) THEN
              N=N+1
              X(N)=AX1(K)
              Y(N)=GDATA1(I,J,K)
            ENDIF
 1310     CONTINUE
          CALL SPLINE(X,Y,N,1E31,1E31,Y2)
          DO 1320 K=1,KDIM1
            IJK=I+IDIM1*(J-1)+IDIM1*JDIM1*(K-1)
            AX=AX1(K)
            CALL SPLINT(X,Y,Y2,N,AX,NDIF,YB)
            GDATAO(IJK)=YB
 1320     CONTINUE
 1300   CONTINUE
*
      RETURN
*===========
      ENTRY SETPAR2(IAXE)
      IAX  = IAXE
      RETURN
*===========
      ENTRY SETDIF (NDIFE)
      NDIF=NDIFE
      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



**************************************************************
* 相対渦度 ---> 絶対渦度
**************************************************************
      SUBROUTINE CORIOL
     I         ( HHEAD , GDATA ,
     O           HHEADS, GDATAS,
     D           IXDIM , IYDIM , IZDIM )
*
      IMPLICIT NONE
#ifdef SYS_IBMS
      INCLUDE    (GZSIZE)           !" NCC, NDC, NED, NFILN
#else
#include         "gzsize.F"
#endif
*
      INTEGER    IXDIM, IYDIM, IZDIM
      CHARACTER  HHEAD  ( * ) *(*)
      CHARACTER  HHEADS ( * ) *(*)
      REAL       GDATA  ( IXDIM, IYDIM, IZDIM )
      REAL       GDATAS ( IXDIM, IYDIM, IZDIM )
      REAL       VMISS
*
      CHARACTER  HHEADA ( NDC ) * ( NCC )
      INTEGER    JMAX
      PARAMETER  ( JMAX=200 )
      REAL       AXIS ( JMAX )
      REAL       F    ( JMAX )

      INTEGER IEOD, J, K, I
      REAL    PI, OMEGA
*
      IF(IYDIM.GT.JMAX) THEN
        CALL MSGDMP('E','CORIOL','ARRAY IS TOO SMALL')
      ENDIF
*
      CALL GTSIZE(HHEADA, JMAX)
      CALL GUQAXV
     I     ( HHEAD , 2    , 'LOC' ,
     O       HHEADA, AXIS , IEOD   )
*
      PI     = ATAN ( 1. )*4.
      OMEGA  = 2. *PI / 86400.
      DO 2000 J=1, IYDIM
        F(J) = 2.*OMEGA*SIN(AXIS(J)*PI/180.) !" コリオリパラメータ
 2000 CONTINUE
*
      CALL GHPGET( HHEAD, 'MISS', VMISS )
*
      DO 6100 K=1,IZDIM
        DO 6200 J=1,IYDIM
          DO 6300 I=1,IXDIM
            IF ( GDATA( I, J, K ).NE.VMISS ) THEN
              GDATAS( I, J, K ) = GDATA( I, J, K ) + F ( J )
            ELSE 
              GDATAS( I, J, K ) = VMISS
            ENDIF
 6300     CONTINUE
 6200   CONTINUE
 6100 CONTINUE
*
      RETURN
      END


*********************************************************************
*" cp from gtdsph
*" ただし地球半径で割る部分を改良した。
*********************************************************************
      SUBROUTINE DSPHE
     I         ( HH    , GDATA , 
     O           HHF   , GDATAF,
     I           IDIMI , JDIMI , KDIMI )
*
      IMPLICIT NONE
      INTEGER    IDIMI, JDIMI, KDIMI
      CHARACTER  HH    ( * )*(*)             !" ヘッダー(入力)
      REAL       GDATA ( IDIMI*JDIMI*KDIMI )  !" データ(入力)
      CHARACTER  HHF   ( * )*(*)             !" ヘッダー(出力)
      REAL       GDATAF( IDIMI*JDIMI*KDIMI )  !" データ(出力)
*
      LOGICAL    LDX,DX  !" T: x 微分 (DY=Tの方が優先される)
      LOGICAL    LDY,DY  !" T: y 微分
      INTEGER    IWB,WB !"(cosφ)**WB の重みを微分前に掛ける
      INTEGER    IWA,WA !"(cosφ)**WA の重みを微分後に掛ける
      LOGICAL    LA,A  !" T: 地球半径で割る
      LOGICAL    LZONAL,ZONAL  !" T: 一番目の軸を y 軸とみなす(波数０成分用)
      SAVE       LDX,LDY,IWB,IWA,LA,LZONAL
*
      REAL       VMISS
      CHARACTER  HGRAD*4
*
#ifdef SYS_IBMS
      INCLUDE    (GTSINC)
      INCLUDE    (gpainc)              !" /GPAWRK/
#else
#include         "gtsinc.F"
#include         "gpainc.F"
#endif
*"    COMMON /GPAWK?/ WDATA , ZDATA , ZWORK ,
*"                    PNM   , DPNM  , TRIGS , IFAX  ,
*"                    NMO   , GW    , COSLAT,
*"                    QPNM  , QDPNM , QGW   , QSINLA
*"    REAL ER / 6370.E3 /
*
*"         < 配列の次元 (incl.どこが y 軸 ?) >
*
      INTEGER IDIMS, JDIMS, KDIMS, IMAXS, JMAXS, KMAXS, JMXHF,
     &        LMAX, MINT, MMAX, NMAX, NMDIM,
     &        I, J, K, IJK

      IF(.NOT.LZONAL) THEN
        IDIMS=IDIMI
        JDIMS=JDIMI
        KDIMS=KDIMI
        IMAXS=2*(IDIMI/2) !" IDIMI が奇数なら IMAX=IDIMI-1 と仮定する
        JMAXS=JDIMI
        KMAXS=KDIMI
      ELSE
        IDIMS=1
        JDIMS=IDIMI
        KDIMS=JDIMI
        IMAXS=1
        JMAXS=IDIMI
        KMAXS=JDIMI
        IF(KDIMI.GT.1) CALL MSGDMP('E','DSPHE',
     +       '3D DATA CANNOT BE PROCESSED IF -ZONAL')
      ENDIF
*
*"         < 球関数の準備 >
*
      CALL GPASPW !" 作業配列の大きさを GPASPS に知らせる
     I         ( IDIM,JDIM,KDIM )
*
      CALL GPASPS !" 三角形切断しか扱えない！
     I         ( IMAXS , JMAXS , KMAXS , ER    ,
     O           PNM   , DPNM  , TRIGS , IFAX  ,
     O           NMO   , GW    , COSLAT,
     O           LMAX  , MMAX  , NMAX  , MINT  ,
     O           NMDIM , JMXHF ,
     W           QPNM  , QDPNM , QGW   , QSINLA         )
*
*"         < 微分前後の重み掛けなど >
*
**      IF (LA) THEN  !" 地球半径で割る   !" commented out by akahori
**        DO 100 J=1,JMAXS
**          GW(J)=GW(J)/ER
**  100   CONTINUE
**      ENDIF
**
      CALL GHPGET( HH,'MISS',VMISS )
      IF(IWB.NE.0) THEN
        DO 200 K=1,KDIMS
        DO 200 J=1,JDIMS
        DO 200 I=1,IDIMS
          IJK = I + (J-1)*IDIMS + (K-1)*IDIMS*JDIMS
          IF( GDATA(IJK).NE.VMISS ) THEN
            GDATA(IJK) = COSLAT(J)**(IWB) * GDATA(IJK)
          ENDIF
  200   CONTINUE
      ENDIF
*
*"         < 微分フラッグ >
*
      IF(LDY) THEN
        HGRAD='YGRA'
      ELSE IF(LDX) THEN
        HGRAD='XGRA'
      ELSE
        HGRAD='    '
      ENDIF
*
*"         < スペクトルに (微分も) >
*
      IF(.NOT.LZONAL) THEN
        CALL SPG2W
     M         ( WDATA ,
     I           GDATA ,
     C           PNM   , NMO   , TRIGS , IFAX  , GW    ,
     F           '    ','POSO' ,
     D           IMAXS , JMAXS , KMAXS , IDIMS , JDIMS ,
     D           LMAX  , MMAX  , NMAX  , MINT  , NMDIM , JMXHF ,
     W           ZDATA , ZWORK                                   )
      ELSE
        CALL SPG2WM
     M         ( WDATA ,
     I           GDATA ,
     C           PNM   , NMO   , TRIGS , IFAX  , GW    ,
     F           '    ','POSO' ,
     D           IMAXS , JMAXS , KMAXS , IDIMS , JDIMS ,
     D           LMAX  , MMAX  , NMAX  , MINT  , NMDIM , JMXHF ,
     W           ZDATA , ZWORK  ) !" in this subroutine ZDATA is dummy
      ENDIF
*
*"         < 格子点値に戻す >
*
      IF(.NOT.LZONAL) THEN
        IF(.NOT.LDY) THEN
          CALL SPW2G
     M         ( GDATAF,
     I           WDATA ,
     C           PNM   , NMO   , TRIGS , IFAX  ,
     F           HGRAD , 'POS' ,
     D           IMAXS , JMAXS , KMAXS , IDIMS , JDIMS ,
     D           LMAX  , MMAX  , NMAX  , MINT  , NMDIM , JMXHF , 
     W           ZDATA , ZWORK                                   )
        ELSE
          CALL SPW2G
     M         ( GDATAF,
     I           WDATA ,
     C           DPNM  , NMO   , TRIGS , IFAX  ,
     F           HGRAD , 'POSO' ,
     D           IMAXS , JMAXS , KMAXS , IDIMS , JDIMS ,
     D           LMAX  , MMAX  , NMAX  , MINT  , NMDIM , JMXHF , 
     W           ZDATA , ZWORK                                   )
        ENDIF
      ELSE
        IF(.NOT.LDY) THEN
          CALL SPW2GM
     M         ( GDATAF,
     I           WDATA ,
     C           PNM   , NMO   , TRIGS , IFAX  ,
     F           HGRAD , 'POS' ,
     D           IMAXS , JMAXS , KMAXS , IDIMS , JDIMS ,
     D           LMAX  , MMAX  , NMAX  , MINT  , NMDIM , JMXHF , 
     W           ZDATA , ZWORK                                   )
        ELSE
          CALL SPW2GM
     M         ( GDATAF,
     I           WDATA ,
     C           DPNM  , NMO   , TRIGS , IFAX  ,
     F           HGRAD , 'POSO' ,
     D           IMAXS , JMAXS , KMAXS , IDIMS , JDIMS ,
     D           LMAX  , MMAX  , NMAX  , MINT  , NMDIM , JMXHF , 
     W           ZDATA , ZWORK                                   )
        ENDIF
      ENDIF
*
      IF(LDY) THEN
        DO 400 K=1,KDIMS !"微分後の重み掛け
        DO 400 J=1,JDIMS
        DO 400 I=1,IDIMS
          IJK = I + (J-1)*IDIMS + (K-1)*IDIMS*JDIMS
          IF( GDATA(IJK).NE.VMISS ) THEN
            GDATAF(IJK) = COSLAT(J)**(IWA-1) * GDATAF(IJK)
          ENDIF
  400   CONTINUE
      ELSE
        DO 500 K=1,KDIMS !"微分後の重み掛け
        DO 500 J=1,JDIMS
        DO 500 I=1,IDIMS
          IJK = I + (J-1)*IDIMS + (K-1)*IDIMS*JDIMS
          IF( GDATA(IJK).NE.VMISS ) THEN
            GDATAF(IJK) = COSLAT(J)**(IWA) * GDATAF(IJK)
          ENDIF
  500   CONTINUE
      ENDIF

      IF (LA) THEN   !" 地球半径で割る   !" moved by akahori
        DO 1000 K=1,KDIMS
        DO 1000 J=1,JDIMS
        DO 1000 I=1,IDIMS
          IJK = I + (J-1)*IDIMS + (K-1)*IDIMS*JDIMS
          IF( GDATA(IJK).NE.VMISS ) THEN
            GDATAF(IJK) = GDATAF(IJK)/ER
          ENDIF
 1000   CONTINUE
      ENDIF

      RETURN
*===========
      ENTRY EDSPHE(DX,DY,WA,WB,A,ZONAL)
      LDX=DX
      LDY=DY
      IWA=WA
      IWB=WB
      LA=A
      LZONAL=ZONAL
*
      RETURN
      END



*********************************************************************
*" cp from gtu2vor
*********************************************************************
      SUBROUTINE U2VOR
     I         ( HU    , GDU,
     I           HV    , GDV,
     O           HVOR  , GDVOR,
     I           IMAX  , JMAX  , KMAX,
     I           IMAX2  , JMAX2  , KMAX2  )
*
      IMPLICIT NONE
      INTEGER    IMAX, JMAX, KMAX
      INTEGER    IMAX2, JMAX2, KMAX2
      CHARACTER  HU     ( * )*(*)        !" ヘッダー(入力)
      REAL       GDU    ( * )            !" データ(入力)
      CHARACTER  HV     ( * )*(*)        !" ヘッダー(入力)
      REAL       GDV    ( * )            !" データ(入力)
      CHARACTER  HVOR   ( * )*(*)        !" ヘッダー(出力)
      REAL       GDVOR  ( * )            !" データ(出力)
*
#include         "gzsize.F"
#include         "gtsinc.F"
#include         "gpainc.F"
*
      REAL  UVFACT ( IDIM*JDIM )
      REAL  WDVOR  ( NMDIMC*KDIM    )       !" スペクトルデータ
      REAL  GWDEL  ( JDIM )              !" 微分用のガウス荷重

      INTEGER JMXHF, LMAX, MINT, MMAX, NMAX, NMDIM,
     &        I, J, IJ, K, IJK 

*
*"         < 球関数の準備 >
*
      CALL GPASPW !" 作業配列の大きさを GPASPS に知らせる
     I         ( IDIM,JDIM,KDIM )
      CALL GPASPS
     I         ( IMAX  , JMAX  , KMAX  , ER    ,
     O           PNM   , DPNM  , TRIGS , IFAX  ,
     O           NMO   , GW    , COSLAT,
     O           LMAX  , MMAX  , NMAX  , MINT  ,
     O           NMDIM , JMXHF ,
     W           QPNM  , QDPNM , QGW   , QSINLA         )
*
*"         < Uファクター >
*
      DO 3210 J = 1, JMAX
        DO 3200 I = 1, IMAX
          IJ=IMAX*(J-1)+I
          UVFACT ( IJ ) = COSLAT(J)
 3200   CONTINUE
 3210 CONTINUE
*
*"         < 微分用のガウス荷重 >
*
      DO 3100 J = 1, JMAX
         GWDEL( J ) = QGW( J ) / ER / ( 1. - QSINLA(J)**2 )
 3100 CONTINUE
*
*"         < 1. ζ のスペクトル >
*
      DO 1100 K = 1, KMAX
         DO 1100 IJ = 1, IMAX*JMAX
           IJK=IMAX*JMAX*(K-1)+IJ
            GDU ( IJK ) = GDU ( IJK ) * UVFACT( IJ )
            GDV ( IJK ) = GDV ( IJK ) * UVFACT( IJ )
 1100 CONTINUE
*
      CALL SPG2W
     M         ( WDVOR ,
     I           GDU   ,
     C           DPNM  , NMO   , TRIGS , IFAX  , GWDEL ,
     F          'YGRA', 'POS ' ,
     D           IMAX  , JMAX  , KMAX  , IMAX  , JMAX  ,
     D           LMAX  , MMAX  , NMAX  , MINT  , NMDIM , JMXHF ,
     W           ZDATA , ZWORK                                   )
*
      IF ( MMAX .GE. 1 ) THEN
      CALL SPG2W
     M         ( WDVOR ,
     I           GDV   ,
     C           PNM   , NMO   , TRIGS , IFAX  , GWDEL ,
     F          'XGRA', 'ADD ' ,
     D           IMAX  , JMAX  , KMAX  , IMAX  , JMAX  ,
     D           LMAX  , MMAX  , NMAX  , MINT  , NMDIM , JMXHF ,
     W           ZDATA , ZWORK                                   )
      ENDIF
*
*"         < 格子点値に戻す >
*
      CALL SPW2G
     M         ( GDVOR ,
     I           WDVOR ,
     C           PNM   , NMO   , TRIGS , IFAX  ,
     F          '    ' , 'POS ',
     D           IMAX  , JMAX  , KMAX  , IMAX  , JMAX  ,
     D           LMAX  , MMAX  , NMAX  , MINT  , NMDIM , JMXHF , 
     W           ZDATA , ZWORK                                   )
*
      RETURN
      END






*******************************************************************
*"    Potential Vorticity 
*"  = - (g/Ps) * { - ∂ｖ/∂σ     * 1/(a cosφ) ∂θ/∂λ
*"                 + ∂ｕ/∂σ     * 1/a ∂θ/∂φ
*"                 + ( f + rotｖ ) * ∂θ/∂σ             }
*"      rot ｖ ≡ 1/(a cosφ){∂ｖ/∂λ - ∂(cosφ ｕ)/∂λ}
*"      単位は PVU = 10^{-6} m^2 K / ( kg s )
*******************************************************************
           SUBROUTINE PTLVOR
     I          ( GPS    ,
     I            GUS    ,
     I            GVS    ,
     I            GVOR   ,
     I            GTHETAS,
     I            GTHETAX,
     I            GTHETAY,
     O            GPV    ,
     I            IMAX, JMAX, KMAX )
      IMPLICIT NONE
*
      INTEGER    IMAX, JMAX, KMAX
#ifdef SYS_IBMS
      INCLUDE    (GTSINC)
      INCLUDE    (GZSIZE)
#else
#include         "gtsinc.F"
#include         "gzsize.F"
#endif
      COMMON     /GMWORK/ MWORK
      REAL       MWORK  ( IJKDIM )
*
*"[入出力データの領域確保]
      REAL  GPS     (IMAX*JMAX)
      REAL  GUS     (IMAX*JMAX, KMAX)
      REAL  GVS     (IMAX*JMAX, KMAX)
      REAL  GVOR    (IMAX*JMAX, KMAX)
      REAL  GTHETAS (IMAX*JMAX, KMAX)
      REAL  GTHETAX (IMAX*JMAX, KMAX)
      REAL  GTHETAY (IMAX*JMAX, KMAX)
      REAL  GPV     (IMAX*JMAX, KMAX)
*
      REAL PVX, PVY, PVZ
      REAL  G
      DATA  G / 9.8 /
*
      INTEGER IJ, K

      DO 10 K  = 1, KMAX
      DO 10 IJ = 1, IMAX*JMAX
        PVX = - GVS(IJ,K) * GTHETAX(IJ,K)
        PVY =   GUS(IJ,K) * GTHETAY(IJ,K)
        PVZ =   GVOR(IJ,K)* GTHETAS(IJ,K)
        GPV(IJ,K) = - G / ( GPS(IJ) * 1E2 ) * ( PVX+PVY+PVZ ) * 1E6
   10 CONTINUE
*
      RETURN
      END
