*%表題   角運動量, 角速度を計算
*%
*%履歴   91/05/14 佐藤正樹       C:\FORT\OORT\MKAMOM.FOR
*%       91/05/25 佐藤正樹       C:\FORT\OORT\MKAMOM.FOR
*%
*%
*-----------------------------------------------------------------------
      PROGRAM  MKAMOM

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

      REAL*8  LAT(-NLAT:NLAT)         ! 緯度
      REAL  UUU(-NLAT:NLAT, 0:NPRS) ! 東西風

      REAL  UUA(-NLAT:NLAT, 0:NPRS) ! 絶対速度
      REAL  AMA(-NLAT:NLAT, 0:NPRS) ! 絶対角運動量
      REAL  AMR(-NLAT:NLAT, 0:NPRS) ! 相対角運動量
      REAL  OMA(-NLAT:NLAT, 0:NPRS) ! 絶対角速度
      REAL  OMR(-NLAT:NLAT, 0:NPRS) ! 相対角速度

      REAL*8  R                     ! 惑星半径
      REAL*8  OMEGA0                ! 自転角速度

      REAL*8  X                     ! 自転軸からの距離


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



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

      INTEGER   IPR  / 10 /          ! 入力ユニット番号
      CHARACTER HIPFL*13 / 'MKAMOM.PAR' / 
     &                               ! 入力パラメーターファイル名
      LOGICAL   LIPEX                ! 入力パラメーターファイルの存在


      INTEGER   IR   / 20 /          ! 入力ユニット番号
      INTEGER  IIFL                  ! 入力ファイル数
      PARAMETER( IIFL=1 )
      CHARACTER HIFL(IIFL)*16        ! 入力ファイル名

      INTEGER   OR   / 30 /          ! 出力ユニット番号
      INTEGER  IOFL                  ! 出力ファイル数
      PARAMETER( IOFL=5 )
      CHARACTER HOFL(IOFL)*16        ! 出力ファイル名

      CHARACTER  HXFORM*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

      INTEGER  IDIM
      INTEGER  JDIM
      INTEGER  KDIM
      INTEGER  IJKDIM



*     < 定数 >

      DOUBLE PRECISION  PI
      PARAMETER( PI = 3.141592653589793 D 0 ) ! 円周率



*     < NAMELIST >

      NAMELIST / FLPARA /  HIFL,
     &                     HOFL

*-----------------------------------------------------------------------
*     < アンダーフローのスイッチの設定 >

      CALL UNDER0(.TRUE.)

*-----------------------------------------------------------------------
*     < 入力形式 >

      HXFORM='(5E16.6)'

*-----------------------------------------------------------------------
*     < 定数 >

      R      = 6378.D3
      OMEGA0 = 2.D0 * PI / 24.D0 / 60.D0 / 60.D0

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

      INQUIRE(FILE=HIPFL, EXIST=LIPEX)
      IF( .NOT. LIPEX ) THEN
        WRITE(6,*) '入力パラメーターファイルが存在しません'
        WRITE(6,*) ' ファイル名: ', HIPFL
        STOP
      END IF


      OPEN(IPR,FILE=HIPFL)
        READ(IPR,FLPARA)
      CLOSE(IPR)

      WRITE(6,*) 'ファイルを読み込みます', HIFL

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

      CALL  RDGTDT
     I           ( IR,       HIFL(1),
     O             HHEAD,    UUU,
     D             IDIM,     JDIM,     KDIM,     IJKDIM  )

      IF( ( IDIM .NE. 1 ) .OR. ( JDIM .NE. (2*NLAT+1) ) 
     &  .OR.  ( KDIM .NE. NPRS+1 ) 
     &  .OR. ( IJKDIM .NE. (2*NLAT+1) * ( NPRS+1) ) ) THEN

         WRITE(6,*) '配列の数がおかしい'
         STOP
       ENDIF

*-----------------------------------------------------------------------
*     < 緯度の作成 >
*     精度が必要なのでここで作る

      DO 10 ILAT = -NLAT, NLAT
        LAT(ILAT) = DBLE(ILAT) / DBLE(NLAT) * PI / 2.
 10   CONTINUE

*-----------------------------------------------------------------------
*     < 角運動量, 角速度の計算 >

      DO 100 ILAT = -NLAT, NLAT
        X = R * COS ( LAT(ILAT) )
        DO 200 IPRS = 0, NPRS
          UUA(ILAT, IPRS) = UUU(ILAT, IPRS) + OMEGA0 * X
          AMA(ILAT, IPRS) = ( UUU(ILAT, IPRS) + OMEGA0 * X ) *X 
          AMR(ILAT, IPRS) = UUU(ILAT, IPRS) * X

          IF( ( ILAT .EQ. -NLAT ) .OR. ( ILAT .EQ. NLAT ) ) THEN
            OMA(ILAT, IPRS) = -999.
            OMR(ILAT, IPRS) = -999.
          ELSE
            OMA(ILAT, IPRS) = UUU(ILAT, IPRS) / X + OMEGA0
            OMR(ILAT, IPRS) = UUU(ILAT, IPRS) / X
          ENDIF

 200  CONTINUE
 100  CONTINUE

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

      HDSET  = 'Oort and Peixoto'

      IFNUM  = 1
      HUTIM  = 'SEC'
      ITDUR  = 31536000.
      HAITM1 = 'GLON1'
      HAITM2 = 'GLAT37'
      HAITM3 = 'GPRS21'

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

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

*-----------------------------------------------------------------------
*     < UUA の出力 >

      HITEM  = 'UUA'
      HTITL  = 'Absolute Zonal  Wind'
      HUNIT  = 'm/s'

      DMIN   = 0.
      DMAX   = 500.
      DIVS   =  50.
      DIVL   =  100.
      ISTYP  =   1

      CALL GTDATA
     I           ( OR,     HOFL(1),    HXFORM, 
     I             UUA,
     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*NLAT+1, NPRS+1 )

*-----------------------------------------------------------------------
*     < AMA の出力 >

      HITEM  = 'AMA'
      HTITL  = 'Absolute Angular Momentum'
      HUNIT  = 'm**2/s'

      DMIN   = 0.
      DMAX   = 3.E9
      DIVS   = 2.E8
      DIVL   = 1.E9
      ISTYP  =   1

      CALL GTDATA
     I           ( OR,     HOFL(2),    HXFORM, 
     I             AMA,
     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*NLAT+1, NPRS+1 )

*-----------------------------------------------------------------------
*     < AMR の出力 >

      HITEM  = 'AMR'
      HTITL  = 'Relative Angular Momentum'
      HUNIT  = 'm**2/s'

      DMIN   = -0.5E8
      DMAX   =  1.5E8
      DIVS   =  0.2E8
      DIVL   =  1.E8
      ISTYP  =   1

      CALL GTDATA
     I           ( OR,     HOFL(3),    HXFORM, 
     I             AMR,
     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*NLAT+1, NPRS+1 )

*-----------------------------------------------------------------------
*     < OMA の出力 >

      HITEM  = 'OMA'
      HTITL  = 'Absolute Angular Velocity'
      HUNIT  = '1/s'

      DMIN   =  6.E-5
      DMAX   =  9.E-5
      DIVS   =  0.1E-6
      DIVL   =  0.5E-5
      ISTYP  =   1

      CALL GTDATA
     I           ( OR,     HOFL(4),    HXFORM, 
     I             OMA,
     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*NLAT+1, NPRS+1 )

*-----------------------------------------------------------------------
*     < OMR の出力 >

      HITEM  = 'OMR'
      HTITL  = 'Relative Angular Velocity'
      HUNIT  = '1/s'

      DMIN   =  -1.E-5
      DMAX   =  2.E-5
      DIVS   =  0.1E-6
      DIVL   =  0.5E-5
      ISTYP  =   1

      CALL GTDATA
     I           ( OR,     HOFL(5),    HXFORM, 
     I             OMR,
     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*NLAT+1, NPRS+1 )

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