#if defined(L08_1A)
! *****************************COPYRIGHT*******************************
! (C) Crown copyright Met Office. All rights reserved.
! For further details please refer to the file COPYRIGHT.txt
! which you should have received as part of this distribution.
! *****************************COPYRIGHT*******************************
!   SUBROUTINE FCDCH-------------------------------------------------

!  Purpose: Calculate surface transfer coefficients at one or more
!           gridpoints.


!  Documentation: UM Documentation Paper No 24, section P243.

!--------------------------------------------------------------------
SUBROUTINE fcdch (                                                &
 row_length,rows,cor_mo_iter,points,tile_pts,                     &
 tile_index,pts_index,                                            &
 db,vshr,z0m,z0h,zh,z1_uv,z1_uv_top,z1_tq,z1_tq_top,              &
 wind_profile_factor,                                             &
 ddmfx,                                                           &
 cdv,chv,cdv_std,v_s,v_s_std,recip_l_mo,u_s_std,                  &
 ltimer                                                           &
)

USE c_vkman
USE surf_param, ONLY :                                            &
                       beta                                       &
                      ,third                                      &
                      ,beta_cndd                                  &
                      ,cnst_cndd_0                                &
                      ,cnst_cndd_1                                &
                      ,cnst_cndd_2                                &
                      ,min_wind                                   &
                      ,min_ustar


#if defined(UM_RUN)
USE blopt8a, ONLY :                                               &
                    off                                           &
                   ,on                                            &
                   ,ip_srfexwithcnv
USE bl_option_mod, ONLY : isrfexcnvgust
#else
USE blopt8a, ONLY :                                               &
                    off                                           &
                   ,on                                            &
                   ,ip_srfexwithcnv                               &
                   ,isrfexcnvgust
#endif

USE switches, ONLY :                                            &
                    i_modiscopt

USE parkind1, ONLY: jprb, jpim
USE yomhook, ONLY: lhook, dr_hook
IMPLICIT NONE

INTEGER                                                           &
 cor_mo_iter          ! IN Switch for MO iteration correction

INTEGER                                                           &
 row_length                                                       &
                      ! IN Number of X points?
,rows                                                             &
                      ! IN Number of Y points?
,points                                                           &
                      ! IN Number of points.
,tile_pts                                                         &
                      ! IN Number of tile points.
,tile_index(points)                                               &
                      ! IN Index of tile points.
,pts_index(points) ! IN Index of land points.

LOGICAL                                                           &
 ltimer


REAL                                                              &
 db(points)                                                       &
               ! IN Buoyancy difference between surface and lowest
!                    !    temperature and humidity level in the
!                    !    atmosphere (m/s^2).
,vshr(row_length,rows)                                            &
                       ! IN Wind speed difference between the surf
!                    !    the lowest wind level in the atmosphere (m/s).
,z0m(points)                                                      &
               ! IN Roughness length for momentum transport (m).
,z0h(points)                                                      &
               ! IN Roughness length for heat and moisture (m).
,zh(row_length,rows)                                              &
                       ! IN Depth of boundary layer (m).
,z1_uv(row_length,rows)                                           &
                       ! IN Height of lowest wind level (m).
,z1_tq(row_length,rows)                                           &
                       ! IN Height of lowest temperature and
!                    !    humidity level (m).
,wind_profile_factor(points)
!                    ! IN for adjusting the surface transfer
!                    !    coefficients to remove form drag effects.
REAL, INTENT(IN) ::                                               &
  z1_uv_top(row_length, rows)
                     ! Height of top of lowest uv-layer
REAL, INTENT(IN) ::                                               &
  z1_tq_top(row_length, rows)
                     !  Height of top of lowest Tq-layer
REAL, INTENT(IN) :: ddmfx(row_length, rows)
!                    !    Convective downdraught mass-flux
!                    !    at cloud-base.
REAL                                                              &
 cdv(points)                                                      &
               ! OUT Surface transfer coefficient for momentum
!                    !     including orographic form drag (m/s).
,chv(points)                                                      &
               ! OUT Surface transfer coefficient for
!                    !     heat, moisture & other scalars (m/s).
,cdv_std(points)                                                  &
!                    ! OUT Surface transfer coefficient for momentum
!                    !     excluding orographic form drag (m/s).
,v_s(points)                                                      &
               ! OUT Surface layer scaling velocity
!                    !     including orographic form drag (m/s).
,v_s_std(points)                                                  &
!                    ! OUT Surface layer scaling velocity
!                    !     excluding orographic form drag (m/s).
,u_s_std(points)                                                  &
               ! OUT Scaling velocity from middle of MO iteration
!                    !     - picked up in error by dust code!
,recip_l_mo(points)
!                    ! OUT Reciprocal of the Monin-Obukhov length
!                    !     (m^-1).

REAL :: wstrcnvgust(row_length, rows)
!                   ! Turbulent velocity scale for convective gustiness
REAL :: wgst_tmp

!    Workspace usage----------------------------------------------------

!     Local work arrays.

REAL                                                              &
 phi_m(points)                                                    &
                 ! Monin-Obukhov stability function for momentum
!                      ! integrated to the model's lowest wind level.
,phi_h(points) ! Monin-Obukhov stability function for scalars
!                      ! integrated to the model's lowest temperature
!                      ! and humidity level.

!*----------------------------------------------------------------------

EXTERNAL phi_m_h, phi_m_h_vol
EXTERNAL timer

!  Define local variables

INTEGER i,j,k,l ! Loop counter; horizontal field index.
INTEGER it    ! Iteration loop counter.
INTEGER n_its ! Number of iterations for Monin-Obukhov length
!                   ! and stability functions.

REAL                                                              &
 b_flux                                                           &
              ! Surface buoyancy flux over air density.
,u_s                                                              &
              ! Iteration surface friction velocity
              ! (effective value)
,u_s2                                                             &
              ! Iteration surface friction velocity squared
              ! (effective value)
,u_s_std2                                                         &
              ! Non-effective version of U_S2
,w_s          ! Surface turbulent convective scaling velocity.

INTEGER(KIND=jpim), PARAMETER :: zhook_in  = 0
INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
REAL(KIND=jprb)               :: zhook_handle

IF (lhook) CALL dr_hook('FCDCH',zhook_in,zhook_handle)

IF (ltimer) THEN
! DEPENDS ON: timer
  CALL timer('FCDCH   ',3)
END IF

!-----------------------------------------------------------------------
!! 0. Initialise values
!-----------------------------------------------------------------------
cdv(:)        = 0.0
chv(:)        = 0.0
cdv_std(:)    = 0.0
v_s(:)        = 0.0
v_s_std(:)    = 0.0
u_s_std(:)    = 0.0
recip_l_mo(:) = 0.0
phi_m(:)      = 0.0
phi_h(:)      = 0.0

!-----------------------------------------------------------------------
!! 1. Set initial values for the iteration.
!-----------------------------------------------------------------------
IF (cor_mo_iter == off) THEN
  n_its=5   ! original iteration count
!                 ! Found typically 2% from converged value
ELSE
  n_its=8   ! Found typically 0.2% from converged value
END IF

IF (isrfexcnvgust == ip_srfexwithcnv) THEN
!CDIR NODEP
  DO k=1,tile_pts
    l = tile_index(k)
    j=(pts_index(l)-1)/row_length + 1
    i = pts_index(l) - (j-1)*row_length

!         Compute convective gustiness, which is unchanged during
!         the iteration. Redelsperger's adjustment is to U_10; here it
!         is scaled to apply to the friction velocity.
    wgst_tmp = MIN( 0.6, MAX(0.0, ddmfx(i, j) ) )
    wstrcnvgust(i,j) = 0.067* LOG ( cnst_cndd_0 +                 &
                                     cnst_cndd_1 * wgst_tmp +      &
                                     cnst_cndd_2 * wgst_tmp ** 2 )
  END DO
!print*,"1maxmin",maxval(wstrcnvgust),minval(wstrcnvgust)
!print*,"2maxmin",maxval(ddmfx),minval(ddmfx)
!call flush(6)

END IF
!CDIR NODEP
DO k=1,tile_pts
  l = tile_index(k)
  j=(pts_index(l)-1)/row_length + 1
  i = pts_index(l) - (j-1)*row_length
  IF (db(l)  <   0.0 .AND. vshr(i,j)  <   2.0) THEN
!-----------------------------------------------------------------------
!         Start the iteration from the convective limit.
!-----------------------------------------------------------------------
    recip_l_mo(l) = -vkman/(beta*beta*beta*zh(i,j))
  ELSE
!-----------------------------------------------------------------------
!         Start the iteration from neutral values.
!-----------------------------------------------------------------------
    recip_l_mo(l) = 0.0
  END IF
END DO
SELECT CASE (i_modiscopt)
  CASE(off)
! DEPENDS ON: phi_m_h
    CALL phi_m_h ( row_length,rows,points,tile_pts,               &
                   tile_index,pts_index,                          &
                   recip_l_mo,z1_uv,z1_tq,z0m,z0h,                &
                   phi_m,phi_h,                                   &
                   ltimer )
  CASE(on)
! DEPENDS ON: phi_m_h_vol
    CALL phi_m_h_vol ( row_length,rows,points,tile_pts,           &
                   tile_index,pts_index,                          &
                   recip_l_mo,z1_uv_top,z1_tq_top,z0m,z0h,        &
                   phi_m,phi_h,                                   &
                   ltimer )
END SELECT

!CDIR NODEP
DO k=1,tile_pts
  l = tile_index(k)
  j=(pts_index(l)-1)/row_length + 1
  i = pts_index(l) - (j-1)*row_length
  IF (db(l)  <   0.0 .AND. vshr(i,j)  <   2.0) THEN
!-----------------------------------------------------------------------
!         Start the iteration from the convective limit.
!-----------------------------------------------------------------------
    v_s_std(l) = beta *                                           &
        SQRT( beta * ( vkman / phi_h(l) ) * zh(i,j) * (-db(l)) )
    v_s(l) = v_s_std(l)
  ELSE
!-----------------------------------------------------------------------
!         Start the iteration from neutral values.
!-----------------------------------------------------------------------
    v_s(l) = ( vkman / phi_m(l) ) * vshr(i,j)
    v_s_std(l) = v_s(l) * wind_profile_factor(l)
  END IF
  chv(l) = ( vkman / phi_h(l) ) * v_s_std(l)
  cdv(l) = ( vkman / phi_m(l) ) * v_s(l)
  cdv_std(l) = cdv(l) * ( v_s_std(l) / v_s(l) ) *                 &
                        wind_profile_factor(l)
END DO
!-----------------------------------------------------------------------
!! 2. Iterate to obtain sucessively better approximations for CD & CH.
!-----------------------------------------------------------------------
DO it = 1,n_its

  IF (cor_mo_iter == off) THEN
!-----------------------------------------------------------------------
!  Original version with incorrect iteration of gustiness
!-----------------------------------------------------------------------
!CDIR NODEP
    DO k=1,tile_pts
      l = tile_index(k)
      j=(pts_index(l)-1)/row_length + 1
      i = pts_index(l) - (j-1)*row_length
      b_flux = -chv(l) * db(l)

      IF( vshr(i,j) > min_wind)THEN
        u_s = SQRT( cdv(l) * vshr(i,j) )
      ELSE
        u_s = min_ustar
      END IF

      u_s_std(l) = SQRT( cdv_std(l) * vshr(i,j) )

      IF (db(l)  <   -EPSILON(0.0)) THEN
        w_s = MAX( (zh(i,j) * b_flux)**third, EPSILON(0.0) )

        SELECT CASE(isrfexcnvgust)
          CASE(off)
            v_s(l) = SQRT(u_s*u_s + beta*beta*w_s*w_s)
            v_s_std(l) =                                            &
            SQRT( u_s_std(l)*u_s_std(l) + beta*beta*w_s*w_s )
          CASE(ip_srfexwithcnv)
            v_s(l) = SQRT(u_s*u_s + beta_cndd*beta_cndd*w_s*w_s +   &
                     wstrcnvgust(i,j)**2 )
            v_s_std(l) = SQRT( u_s_std(l)*u_s_std(l) +              &
            beta_cndd*beta_cndd*w_s*w_s + wstrcnvgust(i,j)**2 )
        END SELECT

      ELSE
        v_s(l) = u_s
        v_s_std(l) = u_s_std(l)
      END IF

      IF( v_s(l) > min_wind)THEN
        recip_l_mo(l) = -vkman * b_flux                           &
                      / (v_s(l)*v_s(l)*v_s(l))
      ELSE
        recip_l_mo(l) = SQRT(EPSILON(0.0))
      END IF
    END DO
  ELSE       ! cor_mo_iter
!-----------------------------------------------------------------------
!  Corrected version
!-----------------------------------------------------------------------
!CDIR NODEP
    DO k=1,tile_pts
      l = tile_index(k)
      j=(pts_index(l)-1)/row_length + 1
      i = pts_index(l) - (j-1)*row_length
      b_flux = -chv(l) * db(l)

      IF( vshr(i,j) > min_wind)THEN
        u_s2 = cdv(l) * vshr(i,j)
      ELSE
        u_s2 = min_ustar
      END IF
      u_s_std2 = cdv_std(l) * vshr(i,j)
      u_s_std(l) = SQRT( u_s_std2 )

      IF (db(l)  <   -EPSILON(0.0)) THEN
        w_s = MAX( (zh(i,j) * b_flux)**third, EPSILON(0.0) )
!           ! Note that, during this iteration, CDV already includes
!           ! this gust enhancement and so U_S2 is not simply the
!           ! friction velocity arising from the mean wind, hence:
        SELECT CASE(isrfexcnvgust)
          CASE(off)
            v_s(l) = SQRT( 0.5*( beta*beta*w_s*w_s                  &
                   + SQRT( (beta*w_s)**4 + 4.0*u_s2*u_s2 )) )
            v_s_std(l) = SQRT( 0.5*( beta*beta*w_s*w_s              &
                       + SQRT( (beta*w_s)**4                        &
                       + 4.0*u_s_std2*u_s_std2 )) )
          CASE(ip_srfexwithcnv)
            v_s(l) = SQRT( 0.5*( beta_cndd*beta_cndd*w_s*w_s        &
                   + wstrcnvgust(i,j)**2                            &
                   + SQRT( ((beta_cndd*w_s)**2                      &
                   + wstrcnvgust(i,j)**2)**2                        &
                       + 4.0*u_s2*u_s2 )) )
            v_s_std(l) = SQRT( 0.5*( beta_cndd*beta_cndd*w_s*w_s    &
                   + wstrcnvgust(i,j)**2                            &
                   + SQRT( ((beta_cndd*w_s)**2                      &
                   + wstrcnvgust(i,j)**2)**2                        &
                   + 4.0*u_s_std2*u_s_std2 )) )
        END SELECT
      ELSE
        v_s(l) = SQRT( u_s2 )
        v_s_std(l) = SQRT( u_s_std2 )
      END IF

      IF( v_s(l) > min_wind)THEN
        recip_l_mo(l) = -vkman * b_flux                           &
                      / (v_s(l)*v_s(l)*v_s(l))
      ELSE
        recip_l_mo(l) = SQRT(EPSILON(0.0))
      END IF
    END DO

  END IF  ! test on COR_MO_ITER

  SELECT CASE (i_modiscopt)
    CASE (off)
! DEPENDS ON: phi_m_h
      CALL phi_m_h ( row_length,rows,points,tile_pts,             &
                     tile_index,pts_index,                        &
                     recip_l_mo,z1_uv,z1_tq,z0m,z0h,              &
                     phi_m,phi_h,                                 &
                     ltimer )
    CASE (on)
! DEPENDS ON: phi_m_h_vol
      CALL phi_m_h_vol ( row_length,rows,points,tile_pts,         &
                     tile_index,pts_index,                        &
                     recip_l_mo,z1_uv_top,z1_tq_top,z0m,z0h,      &
                     phi_m,phi_h,                                 &
                     ltimer )
  END SELECT


!CDIR NODEP
  DO k=1,tile_pts
    l = tile_index(k)
    j=(pts_index(l)-1)/row_length + 1
    i = pts_index(l) - (j-1)*row_length
    chv(l) = ( vkman / phi_h(l) ) * v_s_std(l)
    cdv(l) = ( vkman / phi_m(l) ) * v_s(l)
    cdv_std(l) = cdv(l) * ( v_s_std(l) / v_s(l) ) *               &
                          wind_profile_factor(l)
  END DO
END DO ! Iteration loop

!-----------------------------------------------------------------------
!! Set CD's and CH's to be dimensionless paremters
!-----------------------------------------------------------------------
!CDIR NODEP
DO k=1,tile_pts
  l = tile_index(k)
  j=(pts_index(l)-1)/row_length + 1
  i = pts_index(l) - (j-1)*row_length
  cdv(l) = cdv(l) / vshr(i,j)
  cdv_std(l) = cdv_std(l) / vshr(i,j)
  chv(l) = chv(l) / vshr(i,j)
END DO


IF (ltimer) THEN
! DEPENDS ON: timer
  CALL timer('FCDCH   ',4)
END IF

IF (lhook) CALL dr_hook('FCDCH',zhook_out,zhook_handle)
RETURN
END SUBROUTINE fcdch
#endif
