!=========================================================================
! Bulk Coefficient Calculation Subroutine
!
! History 1998/07/13  Masatsugu Odaka
!         1999/01/20  Masatsugu Odaka
!=========================================================================

      SUBROUTINE CLDRAG
     I   ( TPOT0, 
     O     RICHD, CDRAGV, 
     G     ZGRID, 
     C     GRAV , TSFC , VSFC,
     D     NX   , NXB  , NZB   )

!*************************************************************************
      IMPLICIT REAL*8 ( A-H, O-Z )
!*************************************************************************

      DIMENSION
     \     CDRAGV( -2:NXB ),
     \     RICHD( -2:NXB ),
     \     TSFC ( -2:NXB )

      DIMENSION
     \     TPOT0( -2:NZB ),
     \     ZGRID( -2:NZB )

      REAL*8 KARMAN, HIGHT, UABS, URMIN

!*************************************************************************
! Set Parameter of Louis(1979)
!*************************************************************************

      HIGHT  = ZGRID(0)         ! lowest grid hight
      ROUGH  = 1.0D-2           ! surface roughness
      KARMAN = 0.35D0           ! von karman constant
      PARAM1 = 9.4D0            ! Louis parameter 1
      PARAM2 = PARAM1/2.0D0     ! Louis parameter 2
      URMIN  = 1.0D-2           ! minimum wind velocity for Ri


      CDRAGN = ( KARMAN / ( LOG ( HIGHT / ROUGH ) ) )**2
                                ! bulk coefficient at neutral stability

      PARAM3 = 7.4D0* CDRAGN * PARAM1 * SQRT( HIGHT / ROUGH ) 
                                ! Louis parameter 3

c      WRITE(*,*) 'HIGHT =',HIGHT
c      WRITE(*,*) 'ROUGH =',ROUGH
c      WRITE(*,*) 'KARMAN=',ARMAN      
c      WRITE(*,*) 'PARAM1=',PARAM1
c      WRITE(*,*) 'PARAM2=',PARAM2
c      WRITE(*,*) 'PARAM3=',PARAM3
c      WRITE(*,*) 'CDRAGN=',CDRAGN
c      WRITE(*,*) 'Z/Z0  =',HIGHT/ROUGH
      
!*************************************************************************
! Call Subroutine of Calculation of Bulk Richardson Number
!*************************************************************************

      CALL CLRICH
     I   ( TPOT0, VSFC, URMIN,
     C     GRAV , TSFC, HIGHT,
     O     RICHD,
     D     NX   , NXB , NZB    )

!*************************************************************************
! Calculation of Bulk Coeficiet
!*************************************************************************

      DO IX = 0, NX

         UABS = VSFC

         IF ( RICHD( IX ) .gt. 0.0 ) THEN

            CDRAGV( IX ) = CDRAGN 
     \           /( 1.0D0 + PARAM2 * RICHD( IX ) )**2
            
         ELSE

            CDRAGV( IX ) = CDRAGN *
     \           ( 1.0D0 - PARAM1 * RICHD( IX ) 
     \           /( 1.0D0 + PARAM3 * SQRT( ABS( RICHD( IX ) ) ) ) )

         END IF

      END DO

      RETURN
      END 
