* PACKAGE P2GRAV !" 物理過程(2) 重力波抵抗
* 
********************************************************************* 
      SUBROUTINE GRAVTY          !" 重力波抵抗
     O         ( UTG    , VTG    , TAUXG  , TAUYG  , 
     I           GDZSD  , IDSFC  , 
     I           GDU    , GDV    , GDT    , GDP    , GDTG   , 
     I           GDZ    , GDTM   , GDPM                       )
* 
*" MCFARLANE(1987) SCHEME IN NCAR CCM 
*
*" AROUND 150DAYS DIVERGENCE PROBLEM, ESTIMATE OF BN2 ATBOTTOM MUST BE 
*" CHANGED BY BOVILLE SUGGESTION, AND FC2=2(1:OLD)AS SAME AS NCAR MODEL
* 
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                     !" GRID AND WAVE NUMBER
      INCLUDE   (ZCCOM)                     !" BASIC CONSTANTS
#else
#include        "zcdim.F"                   !" GRID AND WAVE NUMBER
#include        "zccom.F"                   !" BASIC CONSTANTS
#endif
* 
*   [OUTPUT] 
      REAL      UTG   ( IDIM*JDIM, KMAX )   !" U TENDENCY BY GRAVITY 
      REAL      VTG   ( IDIM*JDIM, KMAX )   !" V TENDENCY BY GRAVITY 
      REAL      TAUXG ( IDIM*JDIM, KMAX+1 ) !" U GRAVITY MOMENTUM 
      REAL      TAUYG ( IDIM*JDIM, KMAX+1 ) !" V GRAVITY MOMENTUM 
* 
*   [INPUT] 
      REAL      GDZSD ( IDIM*JDIM )         !" STANDARD DEVIATION 
      INTEGER   IDSFC ( IDIM*JDIM )         !" SURFACE INDEX 
      REAL      GDU   ( IDIM*JDIM, KMAX )   !" WESTERLY 
      REAL      GDV   ( IDIM*JDIM, KMAX )   !" SOUTHERLY WIND 
      REAL      GDT   ( IDIM*JDIM, KMAX )   !" TEMPERATURE 
      REAL      GDP   ( IDIM*JDIM, KMAX )   !" PRESSURE 
      REAL      GDTG  ( IDIM*JDIM )         !" SURFACE TEMPERATURE!使ってない
      REAL      GDZ   ( IDIM*JDIM, KMAX )   !" HEIGHT 
* 
      REAL      GDTM  ( IDIM*JDIM, KMAX+1 ) !" TEMPERATURE(INTERFACE) 
      REAL      GDPM  ( IDIM*JDIM, KMAX+1 ) !" PRESSURE(INTERFACE) 
* 
*   [INTERNAL PARM] 
      REAL       VMIN
      REAL       V0
      REAL       ZSDMIN
      REAL       FC2
      REAL       BN2MIN
      REAL       ALP
      REAL       TNLIM
      REAL       P00
      REAL       THS
      DATA    VMIN   /  2.0     / !" MIN. REFERENCE WIND FOR WAVE DRAG 
      DATA    V0     /  1.0E-05 / !" NEARLY ZERO WIND 
      DATA    ZSDMIN /  10.0    / !" MIN. ZSD FOR GRAVITY WAVE DRAG 
      DATA    FC2    /  2.0     / !" CRITICAL FROUDE NUMBER SQUARED 
      DATA    BN2MIN /  2.5E-05 / !" MIN. BRUNT VAISALA FREQUENCY**2 
      DATA    ALP    /  8.0E-06 / !" EFFICIENCY*WAVENUMBER/2 
      DATA    TNLIM  /  2.9E-03 / !" MAX. V TENDENCY 
      DATA    P00    /  100000. / !" 1000MB 
      DATA    THS    /  300.    / !" BASIC POTENTIAL TEMPERATURE 
* 
      NAMELIST /NMGRAV/
     &         VMIN  , V0    , ZSDMIN  , FC2   , BN2MIN,
     &         ALP   , TNLIM , THS
*
      LOGICAL OFIRST 
      DATA    OFIRST /  .TRUE.  / 
      SAVE
* 
*   [INTERNAL WORK] 
      COMMON /COMWRK/ 
     &          UB    , VB    , VMB   , PTOP  ,
     &          ULI   , EXNER , VP    , VPM   ,
     &          GTAU  , BN2   , BN    , FTINT
      REAL      UB    ( IDIM*JDIM )         !" U AT REFERENCE LEVEL 
      REAL      VB    ( IDIM*JDIM )         !" V AT REFERENCE LEVEL 
      REAL      VMB   ( IDIM*JDIM )         !" MAG.AT REFERENCE LEVEL 
      REAL      PTOP  ( IDIM*JDIM )         !" TOP PRESSURE FOR DRAG 
      REAL      ULI   ( IDIM*JDIM )         !" U TENDENCY LIMIT 
      REAL      EXNER ( IDIM*JDIM, KMAX )   !" EXNER FUNCTION(MIDPOINT)
      REAL      VP    ( IDIM*JDIM, KMAX )   !" PROJECTED WIND 
      REAL      VPM   ( IDIM*JDIM, KMAX+1 ) !" PROJECTED WIND(MIDPOINT)
      REAL      GTAU  ( IDIM*JDIM, KMAX+1 ) !" WAVE MOMENTUM FLUX 
      REAL      BN2   ( IDIM*JDIM, KMAX+1 ) 
      REAL      BN    ( IDIM*JDIM, KMAX+1 ) 
      REAL      FTINT ( IDIM*JDIM, KMAX )   !" FORINTERPOLATION 
*
      INTEGER    IJ, K
      INTEGER    IFPAR, JFPAR
      REAL       AKAPPA, VABST
* 
*          < 0.0 GRAVITY WAVE DRAG PARAMETERS >
* 
      IF ( OFIRST ) THEN 
         WRITE ( 6,* ) ' GRAVITY WAVE DRAG - MCFARLENE(1987) 92/04/16' 
         OFIRST = .FALSE. 
*
         CALL   REWNML ( IFPAR , JFPAR )
         READ   ( IFPAR, NMGRAV, END=390 )
  390    WRITE  ( JFPAR, NMGRAV )
*
         AKAPPA = RAIR / CP 
      ENDIF 
* 
* 
*          < 1.1  SET LOWEST WIND MAGNITUDE AND UNIT VECTOR > 
* 
      DO 1100 IJ = 1, IDIM*JDIM 
         VMB( IJ ) = MAX( SQRT( GDU( IJ,1 )**2+GDV( IJ,1 )**2 ), V0 )
         UB ( IJ ) = GDU( IJ,1 ) / VMB ( IJ ) 
         VB ( IJ ) = GDV( IJ,1 ) / VMB ( IJ ) 
 1100 CONTINUE 
* 
*          < 1.2 PROJECTION OF WIND IN THE REFERENCE LEVEL 
*                NO NEGATIVE AND NO ZERO PERMITTED  >
* 
      DO 1200 K = 1, KMAX 
         DO 1200 IJ = 1, IDIM*JDIM
            VP   ( IJ, K ) = MAX(   GDU( IJ,K )*UB ( IJ ) 
     &                            + GDV( IJ,K )*VB ( IJ ),  V0 )
            EXNER( IJ, K ) = ( GDP( IJ,K )/P00 ) ** AKAPPA 
 1200 CONTINUE 
* 
*          < 2.0 INTEPORATION FUNCTION: THE SAME AS TEMPERATURE >
* 
      DO 2000 K = 2, KMAX 
         DO 2000 IJ = 1, IDIM*JDIM
            FTINT( IJ,K ) = 
     &               ( LOG( GDP ( IJ,K-1 ) ) - LOG( GDPM( IJ,K ) ) ) 
     &             / ( LOG( GDP ( IJ,K-1 ) ) - LOG( GDP ( IJ,K ) ) ) 
 2000 CONTINUE 
* 
*          < 2.1 INTERFACE WIND PROJECTIONS IS THE SAME AS TEMPERATURE* 
* 
      DO 2100 K = 2, KMAX
         DO 2100 IJ = 1, IDIM*JDIM
            VPM ( IJ,K ) = ( 1.  -FTINT( IJ,K ) ) * VP ( IJ,K-1 ) 
     &                   +        FTINT( IJ,K )   * VP ( IJ,K   ) 
 2100 CONTINUE 
* 
      DO 2110 IJ = 1, IDIM*JDIM
            VPM ( IJ,KMAX+1 ) = VP ( IJ,KMAX ) 
            VPM ( IJ,1      ) = VP ( IJ,1    ) 
 2110 CONTINUE
* 
*          < 2.2 INTERFACE BRUNT VAISALA FREQUENCY SQUARED >
* 
      DO 2200 K = 2, KMAX 
         DO 2200 IJ = 1, IDIM*JDIM
            BN2 ( IJ,K ) = GRAV / THS 
     &                     * ( GDT( IJ,K   )/EXNER( IJ,K   ) 
     &                        -GDT( IJ,K-1 )/EXNER( IJ,K-1 )  ) 
     &                     / ( GDZ( IJ,K ) - GDZ( IJ,K-1 )    ) 
            BN2 ( IJ,K ) = MAX( BN2( IJ,K ) , BN2MIN ) 
 2200 CONTINUE 
* 
      DO 2210 IJ = 1, IDIM*JDIM
            BN2 ( IJ,1 )      = BN2 ( IJ,2 ) 
* 
            BN2 ( IJ,KMAX+1 ) = GRAV * GRAV / CP / GDT( IJ,KMAX ) 
            BN2 ( IJ,KMAX+1 ) = MAX( BN2( IJ,KMAX+1 ) , BN2MIN  ) 
 2210 CONTINUE   
* 
*          < 2.3 BRUNT VAISALA FREQUENCY >
* 
      DO 2300 K = 1, KMAX+1 
         DO 2300 IJ = 1, IDIM*JDIM
            BN ( IJ,K ) = SQRT( BN2( IJ,K ) ) 
 2300 CONTINUE 
* 
*          < 2.4 TOP PRESSURE FOR GRAVITY WAVE DRAG >
* 
      DO 2400 IJ = 1, IDIM*JDIM
         PTOP ( IJ ) = EXP( 2.  *LOG( GDPM( IJ,KMAX   ) ) 
     &                         - LOG( GDPM( IJ,KMAX-1 ) )   ) 
 2400 CONTINUE
* 
*          < 3.0 GRAVITY WAVE MOMENTUM FLUX >
* 
      K = 1 
      DO 3000 IJ = 1, IDIM*JDIM
            GTAU( IJ,K ) = ALP/RAIR * VMB( IJ ) * GDZSD( IJ )**2
     &                     * GDPM( IJ,K ) * BN( IJ,K ) / GDTM( IJ,K )
*
            GTAU( IJ,K ) = MIN( GTAU( IJ,K ) , 
     &                          ALP/RAIR/FC2 * VMB( IJ )**3
     &                          * GDPM( IJ,K )
     &                          / ( BN( IJ,K ) * GDTM( IJ,K ) )  ) 
 3000 CONTINUE
* 
      DO 3010 K = 2, KMAX 
         DO 3010 IJ = 1, IDIM*JDIM
            GTAU( IJ,K ) = MIN( GTAU( IJ,K-1 ) , 
     &                          ALP/RAIR/FC2 * VPM( IJ,K )**3
     &                          * GDPM( IJ,K )
     &                          / ( BN( IJ,K ) * GDTM( IJ,K ) )  ) 
 3010 CONTINUE 
* 
* 
      K = KMAX+1 
      DO 3020 IJ = 1, IDIM*JDIM
            GTAU( IJ,K ) = MIN( GTAU( IJ,K-1 ) , 
     &                          ALP/RAIR/FC2 * VPM( IJ,K )**3
     &                          * PTOP( IJ )
     &                          / ( BN( IJ,K ) * GDTM( IJ,K ) )  ) 
 3020 CONTINUE 
* 
*          < 3.1 MOMENTUM TENDENCY AND FLUX >
* 
*               A. OVER NOT SEA IDSFC >= 1 ( SEA=0 ) 
*               B. REF. LEVEL WIND > VMIN 
*               C. STANDARD DEVIATION > ZSDMIN 
* 
      DO 3100 K = 1, KMAX 
         DO 3100 IJ = 1, IDIM*JDIM
            IF ( IDSFC( IJ ) .GE. 1       .AND.
     &           VMB  ( IJ ) .GT. VMIN    .AND.
     &           GDZSD( IJ ) .GT. ZSDMIN        ) THEN
*
               VABST =   ( GTAU ( IJ,K ) - GTAU ( IJ,K+1 ) ) 
     &                 / ( GDPM ( IJ,K ) - GDPM ( IJ,K+1 ) ) * GRAV
               VABST = MIN( VABST, TNLIM )
*
               UTG( IJ,K )   = - UB( IJ ) * VABST   
               VTG( IJ,K )   = - VB( IJ ) * VABST
*
               TAUXG( IJ,K ) = - GTAU( IJ,K ) * UB ( IJ ) 
               TAUYG( IJ,K ) = - GTAU( IJ,K ) * VB ( IJ ) 
            ELSE
*
               UTG  ( IJ,K ) = 0.  
               VTG  ( IJ,K ) = 0.  
               TAUXG( IJ,K ) = 0.  
               TAUYG( IJ,K ) = 0.  
            ENDIF
 3100 CONTINUE 
* 
      RETURN 
      END 
