*%表題   子午面流線関数を計算
*%       南北風から積分して求める
*%
*%履歴   91/05/14 佐藤正樹       C:\FORT\OORT\MKSFM.FOR
*%       91/05/25 佐藤正樹       C:\FORT\OORT\MKSFM.FOR
*%
*%
*-----------------------------------------------------------------------
      PROGRAM  MKSFM

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

      REAL*8  LAT(-NLAT:NLAT)         ! 緯度
      REAL  PRS(0:NPRS)               ! 圧力
      REAL  VVV(-NLAT:NLAT, 0:NPRS)   ! 南北風
      REAL*8  VSB(-NLAT:NLAT, 1:NPRS) ! 副格子点上での南北風

      REAL  SFM(-NLAT:NLAT, 0:NPRS)   ! 子午面流線関数



      REAL*8  R                       ! 惑星半径
      REAL*8  GRAV                    ! 重力加速度
      DOUBLE PRECISION  PI            ! 円周率
      PARAMETER( PI = 3.141592653589793 D 0 )



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



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

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

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

      INTEGER   OR   / 30 /          ! 出力ユニット番号
      INTEGER  IOFL                  ! 出力ファイル数
      PARAMETER( IOFL=1 )
      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



*     < NAMELIST >

      NAMELIST / FLPARA /  HIFL,
     &                     HOFL

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

      CALL UNDER0(.TRUE.)

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

      HXFORM='(5E16.6)'

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

      R      = 6378.D3
      GRAV   = 9.8

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

      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,    PRS,
     D             IDIM,     JDIM,     KDIM,     IJKDIM  )

      IF( IDIM .NE. NPRS+1 ) THEN

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



      CALL  RDGTDT
     I           ( IR,       HIFL(2),
     O             HHEAD,    VVV,
     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,*) '配列の数がおかしい', HIFL(2)
         STOP
       ENDIF

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

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

*-----------------------------------------------------------------------
*     < 流線関数の計算 >

      DO 100 ILAT = -NLAT, NLAT

*       << 副格子上での V >>
        DO 200 IPRS = 1, NPRS
          VSB(ILAT, IPRS) = ( VVV(ILAT, IPRS) + VVV(ILAT, IPRS-1) )
     &                      / 2.D0
 200  CONTINUE


*       << 鉛直積分 >>
        SFM(ILAT, 0) = 0.
        DO 300 IPRS = 1, NPRS

          SFM(ILAT, IPRS) = 2.D0 * PI * R * COS(LAT(ILAT)) / GRAV
     &                      * VSB(ILAT, IPRS)
     &                      * ( PRS(IPRS) - PRS(IPRS-1) ) 
     &                      + SFM(ILAT, IPRS-1)

 300  CONTINUE


*       << 補正 >>
        DO 400 IPRS = 1, NPRS

          SFM(ILAT, IPRS) = SFM(ILAT, IPRS)
     &                      - SFM(ILAT, NPRS)
     &                        * ( PRS(IPRS) - PRS(0) )
     &                        / ( PRS(NPRS) - PRS(0) )

 400  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

*-----------------------------------------------------------------------
*     < 出力 >

      HITEM  = 'SFM'
      HTITL  = 'Mass Stream Function'
      HUNIT  = 'kg/s'

      DMIN   = -5.E10
      DMAX   = 10.E10
      DIVS   = 1.E10
      DIVL   = 5.E10
      ISTYP  =   1

      CALL GTDATA
     I           ( OR,     HOFL(1),    HXFORM, 
     I             SFM,
     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
