! $Id$
MODULE nuage_mod
  PRIVATE

  PUBLIC diagcld1, diagcld2

  CONTAINS

SUBROUTINE diagcld1(paprs, pplay, rain, snow, kbot, ktop, diafra, dialiq)
  USE dimphy
  USE yomcst_mod_h
IMPLICIT NONE

  ! Laurent Li (LMD/CNRS), le 12 octobre 1998
  ! (adaptation du code ECMWF)

  ! Dans certains cas, le schema pronostique des nuages n'est
  ! pas suffisament performant. On a donc besoin de diagnostiquer
  ! ces nuages. Je dois avouer que c'est une frustration.



  ! Arguments d'entree:
  REAL paprs(klon, klev+1) ! pression (Pa) a inter-couche
  REAL pplay(klon, klev) ! pression (Pa) au milieu de couche
  REAL t(klon, klev) ! temperature (K)
  REAL q(klon, klev) ! humidite specifique (Kg/Kg)
  REAL rain(klon) ! pluie convective (kg/m2/s)
  REAL snow(klon) ! neige convective (kg/m2/s)
  INTEGER ktop(klon) ! sommet de la convection
  INTEGER kbot(klon) ! bas de la convection

  ! Arguments de sortie:
  REAL diafra(klon, klev) ! fraction nuageuse diagnostiquee
  REAL dialiq(klon, klev) ! eau liquide nuageuse

  ! Constantes ajustables:
  REAL canva, canvb, canvh
  PARAMETER (canva=2.0, canvb=0.3, canvh=0.4)
  REAL cca, ccb, ccc
  PARAMETER (cca=0.125, ccb=1.5, ccc=0.8)
  REAL ccfct, ccscal
  PARAMETER (ccfct=0.400)
  PARAMETER (ccscal=1.0E+11)
  REAL cetahb, cetamb
  PARAMETER (cetahb=0.45, cetamb=0.80)
  REAL cclwmr
  PARAMETER (cclwmr=1.E-04)
  REAL zepscr
  PARAMETER (zepscr=1.0E-10)

  ! Variables locales:
  INTEGER i, k
  REAL zcc(klon)

  ! Initialisation:

  DO k = 1, klev
    DO i = 1, klon
      diafra(i, k) = 0.0
      dialiq(i, k) = 0.0
    END DO
  END DO

  DO i = 1, klon ! Calculer la fraction nuageuse
    zcc(i) = 0.0
    IF ((rain(i)+snow(i))>0.) THEN
      zcc(i) = cca*log(max(zepscr,(rain(i)+snow(i))*ccscal)) - ccb
      zcc(i) = min(ccc, max(0.0,zcc(i)))
    END IF
  END DO

  DO i = 1, klon ! pour traiter les enclumes
    diafra(i, ktop(i)) = max(diafra(i,ktop(i)), zcc(i)*ccfct)
    IF ((zcc(i)>=canvh) .AND. (pplay(i,ktop(i))<=cetahb*paprs(i, &
      1))) diafra(i, ktop(i)) = max(diafra(i,ktop(i)), max(zcc( &
      i)*ccfct,canva*(zcc(i)-canvb)))
    dialiq(i, ktop(i)) = cclwmr*diafra(i, ktop(i))
  END DO

  DO k = 1, klev ! nuages convectifs (sauf enclumes)
    DO i = 1, klon
      IF (k<ktop(i) .AND. k>=kbot(i)) THEN
        diafra(i, k) = max(diafra(i,k), zcc(i)*ccfct)
        dialiq(i, k) = cclwmr*diafra(i, k)
      END IF
    END DO
  END DO

  RETURN
END SUBROUTINE diagcld1

SUBROUTINE diagcld2(paprs, pplay, t, q, diafra, dialiq)
  USE dimphy
  USE yomcst_mod_h
  USE yoethf_mod_h
IMPLICIT NONE



  ! Arguments d'entree:
  REAL paprs(klon, klev+1) ! pression (Pa) a inter-couche
  REAL pplay(klon, klev) ! pression (Pa) au milieu de couche
  REAL t(klon, klev) ! temperature (K)
  REAL q(klon, klev) ! humidite specifique (Kg/Kg)

  ! Arguments de sortie:
  REAL diafra(klon, klev) ! fraction nuageuse diagnostiquee
  REAL dialiq(klon, klev) ! eau liquide nuageuse

  REAL cetamb
  PARAMETER (cetamb=0.80)
  REAL cloia, cloib, cloic, cloid
  PARAMETER (cloia=1.0E+02, cloib=-10.00, cloic=-0.6, cloid=5.0)
  ! cc      PARAMETER (CLOIA=1.0E+02, CLOIB=-10.00, CLOIC=-0.9, CLOID=5.0)
  REAL rgammas
  PARAMETER (rgammas=0.05)
  REAL crhl
  PARAMETER (crhl=0.15)
  ! cc      PARAMETER (CRHL=0.70)
  REAL t_coup
  PARAMETER (t_coup=234.0)

  ! Variables locales:
  INTEGER i, k, kb, invb(klon)
  REAL zqs, zrhb, zcll, zdthmin(klon), zdthdp
  REAL zdelta, zcor

  ! Fonctions thermodynamiques:
  include "FCTTRE.h"

  ! Initialisation:

  DO k = 1, klev
    DO i = 1, klon
      diafra(i, k) = 0.0
      dialiq(i, k) = 0.0
    END DO
  END DO

  DO i = 1, klon
    invb(i) = klev
    zdthmin(i) = 0.0
  END DO

  DO k = 2, klev/2 - 1
    DO i = 1, klon
      zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) - &
        rd*0.5*(t(i,k)+t(i,k+1))/rcpd/paprs(i, k+1)
      zdthdp = zdthdp*cloia
      IF (pplay(i,k)>cetamb*paprs(i,1) .AND. zdthdp<zdthmin(i)) THEN
        zdthmin(i) = zdthdp
        invb(i) = k
      END IF
    END DO
  END DO

  DO i = 1, klon
    kb = invb(i)
    IF (thermcep) THEN
      zdelta = max(0., sign(1.,rtt-t(i,kb)))
      zqs = r2es*foeew(t(i,kb), zdelta)/pplay(i, kb)
      zqs = min(0.5, zqs)
      zcor = 1./(1.-retv*zqs)
      zqs = zqs*zcor
    ELSE
      IF (t(i,kb)<t_coup) THEN
        zqs = qsats(t(i,kb))/pplay(i, kb)
      ELSE
        zqs = qsatl(t(i,kb))/pplay(i, kb)
      END IF
    END IF
    zcll = cloib*zdthmin(i) + cloic
    zcll = min(1.0, max(0.0,zcll))
    zrhb = q(i, kb)/zqs
    IF (zcll>0.0 .AND. zrhb<crhl) zcll = zcll*(1.-(crhl-zrhb)*cloid)
    zcll = min(1.0, max(0.0,zcll))
    diafra(i, kb) = max(diafra(i,kb), zcll)
    dialiq(i, kb) = diafra(i, kb)*rgammas*zqs
  END DO

  RETURN
END SUBROUTINE diagcld2

END MODULE nuage_mod