#undef DIM_
#undef SIZ_
#undef SUB_

#if   (DIMS == 1)

#define DIM_ (:)
#define SIZ_ (size(TVA,1))
#define SUB_ GETCDM_1D

#elif (DIMS == 2)

#define DIM_ (:,:)
#define SIZ_ (size(TVA,1),size(TVA,2))
#define SUB_ GETCDM_2D

#endif                      


! !INTERFACE:

  subroutine SUB_ (TVA,UA,DZ,TVS,Z0, CM,CN,RI)


! !INPUT PARAMETERS: 

    real,    intent(IN )           :: TVA  DIM_ ! SURFACE AIR TEMPERATURE  (K)
    real,    intent(IN )           :: UA   DIM_ ! SURFACE WIND SPEED       (M/SEC)
    real,    intent(IN )           :: DZ   DIM_ ! HEIGHT OF AIR VALUES     (M)
    real,    intent(IN )           :: TVS  DIM_ ! CANOPY TEMPERATURE       (K)
    real,    intent(IN )           :: Z0   DIM_ ! THE SURFACE ROUGHNESS    (M)
                                                                       
! !OUTPUT PARAMETERS:                                                  
                                                                             
    real,    intent(OUT)           :: CM   DIM_ ! MOMENTUM DRAG COEFF.     (N.D.)
    real,    intent(OUT)           :: CN   DIM_ ! NEUTRAL DRAG COEFF       (N.D.)
    real,    intent(OUT)           :: RI   DIM_ ! RICHARDSON NUMBER        (N.D.)


!  Locals

    real    :: FM     SIZ_
    real    :: PSI    SIZ_
    real    :: R      SIZ_
    real    :: RM     SIZ_
    real    :: DRIDTV SIZ_

    real    :: LOUIS_B
    real    :: LOUIS_C
    real    :: LOUIS_D

    LOUIS_B = LOUIS
    LOUIS_C = LOUIS
    LOUIS_D = LOUIS

!  NEUTRAL COEFFICIENTS

    RM = DZ/Z0 + 1.0
    CN = (MAPL_KARMAN/ALOG(RM))

    CN = CN*CN

!   RICHARDSON NUMBER

    DRIDTV = (MAPL_GRAV*DZ)/(TVA*(UA*UA))
    RI = min(DRIDTV * ( TVA - TVS + (MAPL_GRAV/MAPL_CP)*DZ),MAXRI)

!   FH AND FM FACTORS - - FUNCTIONS OF RI AND DZ/ZO.

    where ( RI < 0.0 ) !   UNSTABLE CASE
       R    = (3.0*LOUIS_B*LOUIS_C)*(CN*sqrt(-RI*RM))
       PSI  = RI / (1.0 + R)
       FM   = (1.0 - (2.0*LOUIS_B)*PSI)
    end where

    where ( RI >= 0.0 )  !   STABLE CASE
       PSI  = sqrt(1.0+LOUIS_D*RI)
       R    = RI/PSI
       FM   = 1.0 /(1.0 + (2.0*LOUIS_B)*R     )
    end where

    CM   = CN*FM

    return

    end subroutine SUB_