* PACKAGE DBNDR    !" 境界条件の適用, 流線関数の計算
*
*"  [HIS] 92/01/09 (momoko)
*"        95/08/10 (takepiro) 境界値をどう記憶するか?
*"        95/08/29 (takepiro) psi=const に対応
*"        96/02/05 (takepiro) bug fix
*"        96/03/08 (takepiro) 熱境界条件を配列記憶
*"        96/03/13 (takepiro) NAMELIST による熱境界条件設定
***********************************************************************
      SUBROUTINE DBNSET         !" 境界条件の設定
     I           ( PSI  , 
     I             ZETA , T   ,
     I             DX   , DZ   )
*
* [PARAM]
#include        "zcdim.F"            !" 格子点数, 波数
*
* [INPUT]
      REAL       PSI  ( 0:NX+1 , 0:NZ+1 )  !" 流線関数
      REAL       ZETA ( 0:NX+1 , 0:NZ+1 )  !" 渦度
      REAL       T    ( 0:NX+1 , 0:NZ+1 )  !" 温度
*
      REAL       DX             !" X方向の格子点間隔
      REAL       DZ             !" Z方向の格子点間隔
*
* [INTERNAL WORK]
      INTEGER    IFPAR, JFPAR   !" NAMELIST 読み書き装置番号
      INTEGER    IUNIT          !" メッセージ出力装置番号
      REAL       PSIBMX, PSIBMN !" 配列最大最小ルーチン用
      REAL       TFBMX, TFBMN   !" 熱境界条件
      INTEGER    IPMIN, JPMIN, KPMIN !" 配列最大最小ルーチン用
      INTEGER    IPMAX, JPMAX, KPMAX !" 配列最大最小ルーチン用
      INTEGER    IX             !" ループ変数
      REAL       VMISS          !" 欠損値
*
* [INTERNAL ONCE]
      LOGICAL    RIGIDB         !" 下端力学境界条件スイッチ
      LOGICAL    RIGIDT         !" 上端力学境界条件スイッチ
      LOGICAL    TFIXB          !" 下端熱境界条件スイッチ
      LOGICAL    TFIXT          !" 上端熱境界条件スイッチ
      REAL       TFBTM          !" 下端境界の温度/温度傾度
      REAL       TFTOP          !" 上端境界の温度/温度傾度
*
      NAMELIST / NMBNDR / 
     &           RIGIDB, RIGIDT, TFIXB, TFIXT, 
     &           TFBTM , TFTOP
*
*  [INTERNAL SAVE]
      LOGICAL    ORIGDB         !" 下端力学境界条件(記憶)スイッチ
      LOGICAL    ORIGDT         !" 上端力学境界条件(記憶)スイッチ
      LOGICAL    OTFIXB         !" 下端熱境界条件(記憶)スイッチ
      LOGICAL    OTFIXT         !" 上端熱境界条件(記憶)スイッチ
      REAL       PSIBTM         !" 下端境界の流線関数値
      REAL       PSITOP         !" 下端境界の流線関数値
      REAL       TFBBTM(NX)     !" 下端境界の温度/温度傾度
      REAL       TFBTOP(NX)     !" 上端境界の温度/温度傾度
*
      DATA       RIGIDB, RIGIDT      !" デフォルトは stress free
     &        / .FALSE., .FALSE./
      DATA       TFIXB , TFIXT       !" デフォルトは温度固定
     &        / .TRUE. , .TRUE. /
*
      SAVE   ORIGDB, ORIGDT, OTFIXB, OTFIXT, 
     &       PSIBTM, PSITOP, TFBBTM, TFBTOP
*
      CALL ZTPGET 
     I         ( 'MISS',
     O            VMISS   )
*
      TFBTM = VMISS
      TFTOP = VMISS
*
*"        < 0.1 境界条件スイッチの読みこみ >
*
      CALL GLPGET( 'MSGUNIT', IUNIT )
      WRITE ( IUNIT, * ) ' BOUNDARY CONDITIONS : DATE=96/03/13' 
*
      CALL   REWNML ( IFPAR , JFPAR )
      READ   ( IFPAR, NMBNDR, END=1190 )
 1190 WRITE  ( JFPAR, NMBNDR  )
*
      ORIGDB = RIGIDB
      ORIGDT = RIGIDT
      OTFIXB = TFIXB
      OTFIXT = TFIXT
*
*"        < 0.2 上端力学境界条件の記憶 >
*
      CALL RVMXMN               !" 最大・最小を求める
     O     ( PSIBMX, PSIBMN, 
     O       IPMAX , JPMAX , KPMAX , 
     O       IPMIN , JPMIN , KPMIN , 
     I       PSI   , NX+2  , 1  , NZ+2 , 
     I         2   , NX+1  , 1  , 1 , NZ+2 , NZ+2  )
*
      IF ( PSIBMX .EQ. PSIBMN )THEN
         PSITOP = PSIBMX
      ELSE
         CALL MSGDMP( 'W','DBNSET','TOP BOUNDARY PSI IS NOT CONSTANT' )
         WRITE ( IUNIT, * ) 'PSI(TOP) MAX, MIN = ', PSIBMN, PSIBMX
         PSITOP = PSIBMX
      ENDIF
      CALL MSGDMP( 'M','DBNSET','TOP BOUNDARY PSI WAS SET INTERNALLY' )
      WRITE ( IUNIT, * ) 'PSI(TOP) = ', PSITOP
*
*"        < 0.3 下端力学境界条件の記憶 >
*
      CALL RVMXMN               !" 最大・最小を求める
     O     ( PSIBMX, PSIBMN, 
     O       IPMAX , JPMAX , KPMAX , 
     O       IPMIN , JPMIN , KPMIN , 
     I       PSI   , NX+2  , 1  , NZ+2 , 
     I       2   , NX+1  , 1  , 1 , 1 , 1  )
*
      IF ( PSIBMX .EQ. PSIBMN )THEN
         PSIBTM = PSIBMX
      ELSE
         CALL MSGDMP
     &    ( 'W','DBNSET','BOTTOM BOUNDARY PSI IS NOT CONSTANT' )
         WRITE ( IUNIT, * ) 'PSI(BOTTOM) MAX, MIN = ', PSIBMN, PSIBMX
         PSIBTM = PSIBMX
      ENDIF
      CALL MSGDMP
     &     ( 'M','DBNSET','BOTTOM BOUNDARY PSI WAS SET INTERNALLY' )
      WRITE ( IUNIT, * ) 'PSI(BOTTOM) = ', PSIBTM
*
*"        < 0.4 上端熱境界条件の記憶 >
*
      IF ( TFTOP .NE. VMISS ) THEN
         CALL SETV
     O        ( TFBTOP,
     I          TFTOP ,
     D          NX     )
         CALL MSGDMP
     &     ( 'M','DBNSET','TOP BOUNDARY COND. WAS SET BY NAMELIST' )
      ELSE
         IF ( OTFIXT ) THEN
            DO 2300 IX = 1, NX
               TFBTOP(IX) = T( IX, NZ+1 )
 2300       CONTINUE 
         ELSE
            DO 2400 IX = 1, NX
               TFBTOP(IX) = ( T(IX,NZ+1) - T(IX,NZ) ) / DZ
 2400       CONTINUE 
         ENDIF
         CALL MSGDMP
     &     ( 'M','DBNSET','TOP BOUNDARY COND. WAS SET INTERNALLY' )
      ENDIF
*
      CALL RVMXMN               !" 最大・最小を求める
     O     ( TFBMX , TFBMN , 
     O       IPMAX , JPMAX , KPMAX , 
     O       IPMIN , JPMIN , KPMIN , 
     I       TFBTOP ,   NX  , 1  ,  1 , 
     I         1   ,   NX  , 1  ,  1 , 1 , 1  )
*
      IF ( TFBMN .EQ. TFBMX ) THEN
         CALL MSGDMP
     &  ( 'M','DBNSET','TOP BOUNDARY T or dT/dz IS CONSTANT' )
         WRITE ( IUNIT, * ) 'TF(TOP) = ', TFBMN
      ELSE
         CALL MSGDMP
     &  ( 'M','DBNSET','TOP BOUNDARY T or dT/dz IS NOT CONSTANT' )
         WRITE ( IUNIT, * ) 'TF(TOP) MAX, MIN = ', TFBMN, TFBMX
      ENDIF
*
*"        < 0.5 下端熱境界条件の記憶 >
*
      IF ( TFBTM .NE. VMISS ) THEN
         CALL SETV
     O          ( TFBBTM,
     I            TFBTM ,
     D            NX     )
         CALL MSGDMP
     &     ( 'M','DBNSET','BOTTOM BOUNDARY COND. WAS SET BY NAMELIST' )
      ELSE
         IF ( OTFIXB ) THEN
            DO 2100 IX = 1, NX
               TFBBTM(IX) = T( IX, 0 )
 2100       CONTINUE 
         ELSE
            DO 2200 IX = 1, NX
               TFBBTM(IX) = ( T(IX,1) - T(IX,0) ) / DZ
 2200       CONTINUE 
         ENDIF
         CALL MSGDMP
     &     ( 'M','DBNSET','BOTTOM BOUNDARY COND. WAS SET INTERNALLY' )
      ENDIF
*
      CALL RVMXMN               !" 最大・最小を求める
     O     ( TFBMX , TFBMN, 
     O       IPMAX , JPMAX , KPMAX , 
     O       IPMIN , JPMIN , KPMIN , 
     I       TFBBTM ,   NX  , 1  ,  1 , 
     I         1   ,   NX  , 1  ,  1 , 1 , 1  )
*
      IF ( TFBMN .EQ. TFBMX ) THEN
         CALL MSGDMP
     &  ( 'M','DBNSET','BOTTOM BOUNDARY T or dT/dz IS CONSTANT' )
         WRITE ( IUNIT, * ) 'TF(BOTTOM) = ', TFBMN
      ELSE
         CALL MSGDMP
     &  ( 'M','DBNSET','BOTTOM BOUNDARY T or dT/dz IS NOT CONSTANT' )
         WRITE ( IUNIT, * ) 'TF(BOTTOM) MAX, MIN = ', TFBMN, TFBMX
      ENDIF
*
      RETURN
* --------------------------------------------------
*"     << 1. 境界条件の適用と流線関数の計算 >>
*
      ENTRY DBNDR
     O       ( PSI  , 
     M         ZETA , T   ,
     I         DX   , DZ   )
*
*"        < 1.1  温度の境界条件 >
*
      CALL BTEMP
     M     ( T      ,
     I       TFBBTM , TFBTOP ,
     I       OTFIXB , OTFIXT , 
     I       DX     , DZ      )
*
*     ------------------------
      ENTRY DBNDRZ      !" 運動学/力学境界条件の適用と流線関数の計算
     O       ( PSI  , 
     M         ZETA , 
     I         DX   , DZ  )
*     ------------------------
*
*"        < 1.2 渦度の境界条件 >
*
      CALL BZETA
     M      ( ZETA   , 
     I        ORIGDB , ORIGDT )
*
*"        < 1.3 流線関数の計算 >
*
      CALL DZ2PSI
     O     ( PSI   , 
     I       ZETA  ,
     I       PSIBTM, PSITOP, 
     I       DX    , DZ      )
*
      RETURN
* --------------------------------------------------------
*"        < 2. 温度だけの境界条件 >
*
      ENTRY DBNDRT      !" 熱的境界条件の適用
     M       ( T    ,
     I         DX   , DZ  )
*
      CALL BTEMP
     M     ( T      ,
     I       TFBBTM , TFBTOP ,
     I       OTFIXB , OTFIXT , 
     I       DX     , DZ       )
*
      RETURN
      END
************************************************************************
      SUBROUTINE BZETA                     !" 渦度について境界条件を適用
     M           ( ZETA   , 
     I             ORIGDB , ORIGDT )
*
* [PARAM]
#include        "zcdim.F"            !" 格子点数, 波数
*
* [MODIFY]
      REAL       ZETA( 0:NX+1 , 0:NZ+1 )   !" 渦度
*
* [INPUT]
      LOGICAL    ORIGDB                    !" 下端力学境界条件スイッチ
      LOGICAL    ORIGDT                    !" 上端力学境界条件スイッチ
*
* [INTERNAL WORK]
      INTEGER    IX                        !" ループカウンタ
      INTEGER    IZ                        !" ループカウンタ
*
*"        < 1. 下端境界条件の適用 >
*
      IF ( ORIGDB ) THEN
         DO 1010 IX = 1, NX
           CALL MSGDMP( 'E', 'BZETA', 'RIGID CONDITION NOT SUPPORTED' )
CCC            ZETA( IX, 0 ) =????  !" rigid condition 
 1010    CONTINUE
      ELSE
         DO 1020 IX = 1, NX
            ZETA( IX, 0 ) = 0.0D0          !" sterss-free condition 
 1020    CONTINUE
      ENDIF
*
*"        < 2. 上端境界条件の適用 >
*
      IF ( ORIGDT ) THEN
         DO 2010 IX = 1, NX
           CALL MSGDMP( 'E', 'BZETA', 'RIGID CONDITION NOT SUPPORTED' )
CCC            ZETA( IX, NZ+1 ) =??? !" rigid condition 
 2010    CONTINUE
      ELSE
         DO 2020 IX = 1, NX
            ZETA( IX, NZ+1 ) = 0.0D0           !" sterss-free condition 
 2020    CONTINUE
      ENDIF
*
*"       < 3. 水平境界条件 の適用 >
*
      DO 3010 IZ = 0, NZ+1
        ZETA(0   ,IZ) = ZETA(NX,IZ)            !" periodic B.C.
        ZETA(NX+1,IZ) = ZETA(1 ,IZ)            !" periodic B.C.
 3010 CONTINUE
*
      RETURN
      END
************************************************************************
      SUBROUTINE BTEMP                     !" 温度について境界条件を適用
     M           ( T      ,
     I             TFB    , TFT    ,
     I             OTFIXB , OTFIXT , 
     I             DX     , DZ      )
*
* [PARAM]
#include        "zcdim.F"            !" 格子点数, 波数
*
* [MODIFY]
      REAL       T( 0:NX+1 , 0:NZ+1 )
*
* [INPUT]
      REAL       TFB(NX)             !" 下端境界の温度/温度傾度
      REAL       TFT(NX)             !" 上端境界の温度/温度傾度
*
      LOGICAL    OTFIXB              !" 下端熱境界条件スイッチ
      LOGICAL    OTFIXT              !" 上端熱境界条件スイッチ
*
      REAL       DX                  !" X方向の格子点間隔
      REAL       DZ                  !" Z方向の格子点間隔
*
* [INTERNAL WORK]
      INTEGER    IX
      INTEGER    IZ
*
*"        < 1. 下端境界条件の適用 >
*
      IF ( OTFIXB ) THEN
         DO 1010 IX = 1, NX
            T( IX ,0 ) = TFB(IX)
 1010    CONTINUE
      ELSE
         DO 1020 IX = 1, NX
            T( IX, 0 ) = T( IX, 1 ) - TFB(IX)*DZ
 1020    CONTINUE
      ENDIF
*
*"        < 2. 上端境界条件の適用 >
*
      IF ( OTFIXT ) THEN
         DO 2010 IX = 1, NX
            T( IX, NZ+1 ) = TFT(IX)
 2010    CONTINUE
      ELSE
         DO 2020 IX = 1, NX
            T( IX, NZ+1 ) =  T( IX, NZ ) + TFT(IX)*DZ
 2020    CONTINUE
      ENDIF
*
*"        < 3. 水平境界条件の適用 >
*
      DO 3010 IZ = 0, NZ+1
        T( 0   , IZ ) = T( NX ,IZ )          !" periodic B.C. 
        T( NX+1, IZ ) = T( 1  ,IZ )          !" periodic B.C. 
 3010 CONTINUE
*
      RETURN
      END

