source: LMDZ6/branches/Amaury_dev/libf/phylmd/coefcdrag.F90 @ 5220

SUBROUTINE coefcdrag(klon, knon, nsrf, zxli, &
        speed, t, q, zgeop, psol, &
        ts, qsurf, rugos, okri, ri1, &
        cdram, cdrah, cdran, zri1, pref)

  USE indice_sol_mod
  USE lmdz_abort_physic, ONLY: abort_physic
  USE lmdz_clesphys
  USE lmdz_yoethf
  USE lmdz_yomcst

  IMPLICIT NONE
  !-------------------------------------------------------------------------
  ! Objet : calcul des cdrags pour le moment (cdram) et les flux de chaleur
  !         sensible et latente (cdrah), du cdrag neutre (cdran),
  !         du nombre de Richardson entre la surface et le niveau de reference
  !         (zri1) et de la pression au niveau de reference (pref).

  ! I. Musat, 01.07.2002
  !-------------------------------------------------------------------------

  ! klon----input-I- dimension de la grille physique (= nb_pts_latitude X nb_pts_longitude)
  ! knon----input-I- nombre de points pour un type de surface
  ! nsrf----input-I- indice pour le type de surface; voir indice_sol_mod.F90
  ! zxli----input-L- TRUE si calcul des cdrags selon Laurent Li
  ! speed---input-R- module du vent au 1er niveau du modele
  ! t-------input-R- temperature de l'air au 1er niveau du modele
  ! q-------input-R- humidite de l'air au 1er niveau du modele
  ! zgeop---input-R- geopotentiel au 1er niveau du modele
  ! psol----input-R- pression au sol
  ! ts------input-R- temperature de l'air a la surface
  ! qsurf---input-R- humidite de l'air a la surface
  ! rugos---input-R- rugosite
  ! okri----input-L- TRUE si on veut tester le nb. Richardson entre la sfce
  !                  et zref par rapport au Ri entre la sfce et la 1ere couche
  ! ri1-----input-R- nb. Richardson entre la surface et la 1ere couche

  ! cdram--output-R- cdrag pour le moment
  ! cdrah--output-R- cdrag pour les flux de chaleur latente et sensible
  ! cdran--output-R- cdrag neutre
  ! zri1---output-R- nb. Richardson entre la surface et la couche zgeop/RG
  ! pref---output-R- pression au niveau zgeop/RG

  INTEGER, INTENT(IN) :: klon, knon, nsrf
  LOGICAL, INTENT(IN) :: zxli
  REAL, DIMENSION(klon), INTENT(IN) :: speed, t, q, zgeop, psol
  REAL, DIMENSION(klon), INTENT(IN) :: ts, qsurf, rugos, ri1
  LOGICAL, INTENT(IN) :: okri

  REAL, DIMENSION(klon), INTENT(OUT) :: cdram, cdrah, cdran, zri1, pref
  !-------------------------------------------------------------------------

  ! Quelques constantes :
  REAL, parameter :: RKAR = 0.40, CB = 5.0, CC = 5.0, CD = 5.0, cepdu2 = (0.1)**2

  ! Variables locales :
  INTEGER :: i
  REAL, DIMENSION(klon) :: zdu2, zdphi, ztsolv, ztvd
  REAL, DIMENSION(klon) :: zscf, friv, frih, zucf, zcr
  REAL, DIMENSION(klon) :: zcfm1, zcfh1
  REAL, DIMENSION(klon) :: zcfm2, zcfh2
  REAL, DIMENSION(klon) :: trm0, trm1

  CHARACTER (LEN = 80) :: abort_message
  CHARACTER (LEN = 20) :: modname = 'coefcdra'

  !-------------------------------------------------------------------------
  REAL :: fsta, fins, x
  fsta(x) = 1.0 / (1.0 + 10.0 * x * (1 + 8.0 * x))
  fins(x) = SQRT(1.0 - 18.0 * x)
  !-------------------------------------------------------------------------

  abort_message = 'obsolete, remplace par cdrag, use at you own risk'
  CALL abort_physic(modname, abort_message, 1)

  DO i = 1, knon

    zdphi(i) = zgeop(i)
    zdu2(i) = max(cepdu2, speed(i)**2)
    pref(i) = exp(log(psol(i)) - zdphi(i) / (RD * t(i) * &
            (1. + RETV * max(q(i), 0.0))))
    ztsolv(i) = ts(i)
    !       ztvd(i) = t(i) * (psol(i)/pref(i))**RKAPPA
    !       ztvd(i) = (t(i)+zdphi(i)/RCPD/(1.+RVTMP2*q(i))) &
    !          *(1.+RETV*q(i))
    ztvd(i) = (t(i) + zdphi(i) / RCPD / (1. + RVTMP2 * q(i)))
    trm0(i) = 1. + RETV * max(qsurf(i), 0.0)
    trm1(i) = 1. + RETV * max(q(i), 0.0)
    ztsolv(i) = ztsolv(i) * trm0(i)
    ztvd(i) = ztvd(i) * trm1(i)
    zri1(i) = zdphi(i) * (ztvd(i) - ztsolv(i)) / (zdu2(i) * ztvd(i))

    ! on teste zri1 par rapport au Richardson de la 1ere couche ri1

    !IM +++
    IF(1==0) THEN
      IF (okri) THEN
        IF (ri1(i)>=0.0.AND.zri1(i)<0.0) THEN
          zri1(i) = ri1(i)
        ELSE IF(ri1(i)<0.0.AND.zri1(i)>=0.0) THEN
          zri1(i) = ri1(i)
        ENDIF
      ENDIF
    ENDIF
    !IM ---

    cdran(i) = (RKAR / log(1. + zdphi(i) / (RG * rugos(i))))**2

    IF (zri1(i) >= 0.) THEN

      ! situation stable : pour eviter les inconsistances dans les cas
      ! tres stables on limite zri1 a 20. cf Hess et al. (1995)

      zri1(i) = min(20., zri1(i))

      IF (.NOT.zxli) THEN
        zscf(i) = SQRT(1. + CD * ABS(zri1(i)))
        friv(i) = max(1. / (1. + 2. * CB * zri1(i) / zscf(i)), f_ri_cd_min)
        zcfm1(i) = cdran(i) * friv(i)
        frih(i) = max(1. / (1. + 3. * CB * zri1(i) * zscf(i)), f_ri_cd_min)
        !           zcfh1(i) = cdran(i) * frih(i)
        zcfh1(i) = f_cdrag_ter * cdran(i) * frih(i)
        IF(nsrf==is_oce) zcfh1(i) = f_cdrag_oce * cdran(i) * frih(i)
        cdram(i) = zcfm1(i)
        cdrah(i) = zcfh1(i)
      ELSE
        cdram(i) = cdran(i) * fsta(zri1(i))
        cdrah(i) = cdran(i) * fsta(zri1(i))
      ENDIF

    ELSE

      ! situation instable

      IF (.NOT.zxli) THEN
        zucf(i) = 1. / (1. + 3.0 * CB * CC * cdran(i) * SQRT(ABS(zri1(i)) &
                * (1.0 + zdphi(i) / (RG * rugos(i)))))
        zcfm2(i) = cdran(i) * max((1. - 2.0 * CB * zri1(i) * zucf(i)), f_ri_cd_min)
        !           zcfh2(i) = cdran(i)*max((1.-3.0*CB*zri1(i)*zucf(i)),f_ri_cd_min)
        zcfh2(i) = f_cdrag_ter * cdran(i) * max((1. - 3.0 * CB * zri1(i) * zucf(i)), f_ri_cd_min)
        cdram(i) = zcfm2(i)
        cdrah(i) = zcfh2(i)
      ELSE
        cdram(i) = cdran(i) * fins(zri1(i))
        cdrah(i) = cdran(i) * fins(zri1(i))
      ENDIF

      ! cdrah sur l'ocean cf. Miller et al. (1992)

      zcr(i) = (0.0016 / (cdran(i) * SQRT(zdu2(i)))) * ABS(ztvd(i) - ztsolv(i)) &
              **(1. / 3.)
      !         IF (nsrf.EQ.is_oce) cdrah(i) = cdran(i)*(1.0+zcr(i)**1.25) &
      !                  **(1./1.25)
      IF (nsrf==is_oce) cdrah(i) = f_cdrag_oce * cdran(i) * (1.0 + zcr(i)**1.25) &
              **(1. / 1.25)
    ENDIF

  END DO

END SUBROUTINE coefcdrag