trunk/phylmd/fisrtilp.f90

module fisrtilp_m

  IMPLICIT NONE

contains

  SUBROUTINE fisrtilp(dtime, paprs, pplay, t, q, ptconv, ratqs, d_t, d_q, &
       d_ql, rneb, radliq, rain, snow, pfrac_impa, pfrac_nucl, pfrac_1nucl, &
       frac_impa, frac_nucl, prfl, psfl, rhcl)

    ! From phylmd/fisrtilp.F, version 1.2, 2004/11/09 16:55:40
    ! First author: Z. X. Li (LMD/CNRS), 20 mars 1995

    ! Objet : condensation et pr\'ecipitation stratiforme, sch\'ema de
    ! nuage, sch\'ema de condensation \`a grande \'echelle (pluie).

    USE comfisrtilp, ONLY: cld_lc_con, cld_lc_lsc, cld_tau_con, &
         cld_tau_lsc, coef_eva, ffallv_con, ffallv_lsc, iflag_pdf, reevap_ice
    USE dimphy, ONLY: klev, klon
    USE fcttre, ONLY: foede, foeew
    USE numer_rec_95, ONLY: nr_erf
    USE suphec_m, ONLY: rcpd, rd, retv, rg, rlstt, rlvtt, rtt
    USE yoethf_m, ONLY: r2es, r5ies, r5les, rvtmp2

    REAL, INTENT (IN):: dtime ! intervalle du temps (s) 
    REAL, INTENT (IN):: paprs(klon, klev+1) ! pression a inter-couche 
    REAL, INTENT (IN):: pplay(klon, klev) ! pression au milieu de couche
    REAL, INTENT (IN):: t(klon, klev) ! temperature (K)
    REAL, INTENT (IN):: q(klon, klev) ! humidite specifique (kg/kg)
    LOGICAL, INTENT (IN):: ptconv(klon, klev)

    REAL, INTENT (IN):: ratqs(klon, klev)
    ! determine la largeur de distribution de vapeur

    REAL, INTENT (out):: d_t(klon, klev) ! incrementation de la temperature (K)
    REAL, INTENT (out):: d_q(klon, klev) ! incrementation de la vapeur d'eau
    REAL, INTENT (out):: d_ql(klon, klev) ! incrementation de l'eau liquide
    REAL, INTENT (out):: rneb(klon, klev) ! fraction nuageuse

    REAL, INTENT (out):: radliq(klon, klev)
    ! eau liquide utilisee dans rayonnement

    REAL, INTENT (out):: rain(klon) ! pluies (mm/s)
    REAL, INTENT (out):: snow(klon) ! neige (mm/s)

    ! Coeffients de fraction lessivee :
    REAL, INTENT (inout):: pfrac_impa(klon, klev)
    REAL, INTENT (inout):: pfrac_nucl(klon, klev)
    REAL, INTENT (inout):: pfrac_1nucl(klon, klev)

    ! Fraction d'aerosols lessivee par impaction
    REAL, INTENT (out):: frac_impa(klon, klev)

    ! Fraction d'aerosols lessivee par nucleation
    REAL, INTENT (out):: frac_nucl(klon, klev)

    REAL, INTENT (out):: prfl(klon, klev+1)
    ! flux d'eau precipitante aux interfaces (kg/m2/s)

    REAL, INTENT (out):: psfl(klon, klev+1)
    ! flux d'eau precipitante aux interfaces (kg/m2/s)

    REAL, INTENT (out):: rhcl(klon, klev) ! humidite relative en ciel clair 

    ! Local:

    REAL zct(klon), zcl(klon)

    ! Options du programme:

    REAL seuil_neb ! un nuage existe vraiment au-dela
    PARAMETER (seuil_neb=0.001)

    INTEGER ninter ! sous-intervals pour la precipitation
    PARAMETER (ninter=5)
    LOGICAL evap_prec ! evaporation de la pluie 
    PARAMETER (evap_prec=.TRUE.)
    REAL zpdf_sig(klon), zpdf_k(klon), zpdf_delta(klon)
    REAL zpdf_a(klon), zpdf_b(klon), zpdf_e1(klon), zpdf_e2(klon)

    LOGICAL cpartiel ! condensation partielle 
    PARAMETER (cpartiel=.TRUE.)
    REAL t_coup
    PARAMETER (t_coup=234.0)

    INTEGER i, k, n, kk
    REAL zqs(klon), zdqs(klon), zcor, zcvm5
    logical zdelta
    REAL zrfl(klon), zrfln(klon), zqev, zqevt
    REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
    REAL ztglace, zt(klon)
    INTEGER nexpo ! exponentiel pour glace/eau 
    REAL zdz(klon), zrho(klon), ztot(klon), zrhol(klon)
    REAL zchau(klon), zfroi(klon), zfice(klon), zneb(klon)

    LOGICAL:: appel1er = .TRUE.

    ! Variables traceurs:
    ! Provisoire !!! Parametres alpha du lessivage
    ! A priori on a 4 scavenging numbers possibles

    REAL, save:: a_tr_sca(4)

    ! Variables intermediaires

    REAL zalpha_tr
    REAL zfrac_lessi
    REAL zprec_cond(klon)
    REAL zmair, zcpair, zcpeau
    ! Pour la conversion eau-neige
    REAL zlh_solid(klon), zm_solid

    !---------------------------------------------------------------

    zdelq = 0.0

    IF (appel1er) THEN
       PRINT *, 'fisrtilp, ninter:', ninter
       PRINT *, 'fisrtilp, evap_prec:', evap_prec
       PRINT *, 'fisrtilp, cpartiel:', cpartiel
       IF (abs(dtime / real(ninter) - 360.) > 0.001) THEN
          PRINT *, "fisrtilp : ce n'est pas pr\'evu, voir Z. X. Li", dtime
          PRINT *, "Je pr\'ef\`ere un sous-intervalle de 6 minutes."
       END IF
       appel1er = .FALSE.

       ! initialiation provisoire
       a_tr_sca(1) = -0.5
       a_tr_sca(2) = -0.5
       a_tr_sca(3) = -0.5
       a_tr_sca(4) = -0.5

       ! Initialisation a 1 des coefs des fractions lessivees
       DO k = 1, klev
          DO i = 1, klon
             pfrac_nucl(i, k) = 1.
             pfrac_1nucl(i, k) = 1.
             pfrac_impa(i, k) = 1.
          END DO
       END DO
    END IF

    ! Initialisation a 0 de zoliq
    DO i = 1, klon
       zoliq(i) = 0.
    END DO
    ! Determiner les nuages froids par leur temperature
    ! nexpo regle la raideur de la transition eau liquide / eau glace.

    ztglace = rtt - 15.0
    nexpo = 6

    ! Initialiser les sorties:

    DO k = 1, klev + 1
       DO i = 1, klon
          prfl(i, k) = 0.0
          psfl(i, k) = 0.0
       END DO
    END DO

    DO k = 1, klev
       DO i = 1, klon
          d_t(i, k) = 0.0
          d_q(i, k) = 0.0
          d_ql(i, k) = 0.0
          rneb(i, k) = 0.0
          radliq(i, k) = 0.0
          frac_nucl(i, k) = 1.
          frac_impa(i, k) = 1.
       END DO
    END DO
    DO i = 1, klon
       rain(i) = 0.0
       snow(i) = 0.0
    END DO

    ! Initialiser le flux de precipitation a zero

    DO i = 1, klon
       zrfl(i) = 0.0
       zneb(i) = seuil_neb
    END DO

    ! Pour plus de securite

    zalpha_tr = 0.
    zfrac_lessi = 0.

    loop_vertical: DO k = klev, 1, -1
       DO i = 1, klon
          zt(i) = t(i, k)
          zq(i) = q(i, k)
       END DO

       ! Calculer la varition de temp. de l'air du a la chaleur sensible
       ! transporter par la pluie.
       ! Il resterait a rajouter cet effet de la chaleur sensible sur les
       ! flux de surface, du a la diff. de temp. entre le 1er niveau et la
       ! surface.

       DO i = 1, klon
          IF (k <= klev - 1) THEN
             zmair = (paprs(i, k)-paprs(i, k+1))/rg
             zcpair = rcpd*(1.0+rvtmp2*zq(i))
             zcpeau = rcpd*rvtmp2
             zt(i) = ((t(i, k + 1) + d_t(i, k + 1)) * zrfl(i) * dtime &
                  * zcpeau + zmair * zcpair* zt(i)) &
                  / (zmair * zcpair + zrfl(i) * dtime * zcpeau)
          END IF
       END DO

       IF (evap_prec) THEN
          ! Calculer l'evaporation de la precipitation
          DO i = 1, klon
             IF (zrfl(i)>0.) THEN
                zqs(i) = r2es*foeew(zt(i), rtt >= zt(i))/pplay(i, k)
                zqs(i) = min(0.5, zqs(i))
                zcor = 1./(1.-retv*zqs(i))
                zqs(i) = zqs(i)*zcor
                zqev = max(0.0, (zqs(i)-zq(i))*zneb(i))
                zqevt = coef_eva*(1.0-zq(i)/zqs(i))*sqrt(zrfl(i))* &
                     (paprs(i, k)-paprs(i, k+1))/pplay(i, k)*zt(i)*rd/rg
                zqevt = max(0.0, min(zqevt, zrfl(i)))*rg*dtime/ &
                     (paprs(i, k)-paprs(i, k+1))
                zqev = min(zqev, zqevt)
                zrfln(i) = zrfl(i) - zqev*(paprs(i, k)-paprs(i, k+1))/rg/dtime

                ! pour la glace, on r\'e\'evapore toute la pr\'ecip dans la
                ! couche du dessous la glace venant de la couche du
                ! dessus est simplement dans la couche du dessous.

                IF (zt(i)<t_coup .AND. reevap_ice) zrfln(i) = 0.

                zq(i) = zq(i) - (zrfln(i)-zrfl(i))*(rg/(paprs(i, k)-paprs(i, &
                     k+1)))*dtime
                zt(i) = zt(i) + (zrfln(i)-zrfl(i))*(rg/(paprs(i, k)-paprs(i, &
                     k+1)))*dtime*rlvtt/rcpd/(1.0+rvtmp2*zq(i))
                zrfl(i) = zrfln(i)
             END IF
          END DO
       END IF

       ! Calculer Qs et L/Cp*dQs/dT:

       DO i = 1, klon
          zdelta = rtt >= zt(i)
          zcvm5 = merge(r5ies*rlstt, r5les*rlvtt, zdelta)
          zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zq(i))
          zqs(i) = r2es*foeew(zt(i), zdelta)/pplay(i, k)
          zqs(i) = min(0.5, zqs(i))
          zcor = 1./(1.-retv*zqs(i))
          zqs(i) = zqs(i)*zcor
          zdqs(i) = foede(zt(i), zdelta, zcvm5, zqs(i), zcor)
       END DO

       ! Determiner la condensation partielle et calculer la quantite
       ! de l'eau condensee:

       IF (cpartiel) THEN
          ! Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
          ! nuageuse a partir des PDF de Sandrine Bony.
          ! rneb : fraction nuageuse
          ! zqn : eau totale dans le nuage
          ! zcond : eau condensee moyenne dans la maille.

          ! on prend en compte le r\'echauffement qui diminue
          ! la partie condens\'ee

          ! Version avec les ratqs

          IF (iflag_pdf==0) THEN
             DO i = 1, klon
                zdelq = min(ratqs(i, k), 0.99)*zq(i)
                rneb(i, k) = (zq(i)+zdelq-zqs(i))/(2.0*zdelq)
                zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
             END DO
          ELSE
             ! Version avec les nouvelles PDFs.
             DO i = 1, klon
                IF (zq(i) < 1E-15) THEN
                   zq(i) = 1E-15
                END IF
             END DO
             DO i = 1, klon
                zpdf_sig(i) = ratqs(i, k)*zq(i)
                zpdf_k(i) = -sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
                zpdf_delta(i) = log(zq(i)/zqs(i))
                zpdf_a(i) = zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
                zpdf_b(i) = zpdf_k(i)/(2.*sqrt(2.))
                zpdf_e1(i) = zpdf_a(i) - zpdf_b(i)
                zpdf_e1(i) = sign(min(abs(zpdf_e1(i)), 5.), zpdf_e1(i))
                zpdf_e1(i) = 1. - nr_erf(zpdf_e1(i))
                zpdf_e2(i) = zpdf_a(i) + zpdf_b(i)
                zpdf_e2(i) = sign(min(abs(zpdf_e2(i)), 5.), zpdf_e2(i))
                zpdf_e2(i) = 1. - nr_erf(zpdf_e2(i))
                IF (zpdf_e1(i)<1.E-10) THEN
                   rneb(i, k) = 0.
                   zqn(i) = zqs(i)
                ELSE
                   rneb(i, k) = 0.5*zpdf_e1(i)
                   zqn(i) = zq(i)*zpdf_e2(i)/zpdf_e1(i)
                END IF
             END DO
          END IF

          DO i = 1, klon
             IF (rneb(i, k)<=0.0) zqn(i) = 0.0
             IF (rneb(i, k)>=1.0) zqn(i) = zq(i)
             rneb(i, k) = max(0., min(1., rneb(i, k)))
             ! On ne divise pas par 1 + zdqs pour forcer \`a avoir l'eau
             ! pr\'edite par la convection. Attention : il va falloir
             ! verifier tout ca.
             zcond(i) = max(0., zqn(i)-zqs(i))*rneb(i, k)
             rhcl(i, k) = (zqs(i)+zq(i)-zdelq)/2./zqs(i)
             IF (rneb(i, k) <= 0.) rhcl(i, k) = zq(i) / zqs(i)
             IF (rneb(i, k) >= 1.) rhcl(i, k) = 1.
          END DO
       ELSE
          DO i = 1, klon
             IF (zq(i)>zqs(i)) THEN
                rneb(i, k) = 1.0
             ELSE
                rneb(i, k) = 0.0
             END IF
             zcond(i) = max(0.0, zq(i)-zqs(i))/(1.+zdqs(i))
          END DO
       END IF

       DO i = 1, klon
          zq(i) = zq(i) - zcond(i)
          zt(i) = zt(i) + zcond(i)*rlvtt/rcpd/(1.0+rvtmp2*zq(i))
       END DO

       ! Partager l'eau condensee en precipitation et eau liquide nuageuse

       DO i = 1, klon
          IF (rneb(i, k)>0.0) THEN
             zoliq(i) = zcond(i)
             zrho(i) = pplay(i, k)/zt(i)/rd
             zdz(i) = (paprs(i, k)-paprs(i, k+1))/(zrho(i)*rg)
             zfice(i) = 1.0 - (zt(i)-ztglace)/(273.13-ztglace)
             zfice(i) = min(max(zfice(i), 0.0), 1.0)
             zfice(i) = zfice(i)**nexpo
             zneb(i) = max(rneb(i, k), seuil_neb)
             radliq(i, k) = zoliq(i)/real(ninter+1)
          END IF
       END DO

       DO n = 1, ninter
          DO i = 1, klon
             IF (rneb(i, k)>0.0) THEN
                zrhol(i) = zrho(i)*zoliq(i)/zneb(i)

                IF (ptconv(i, k)) THEN
                   zcl(i) = cld_lc_con
                   zct(i) = 1./cld_tau_con
                ELSE
                   zcl(i) = cld_lc_lsc
                   zct(i) = 1./cld_tau_lsc
                END IF
                ! quantit\'e d'eau \`a \'eliminer
                zchau(i) = zct(i)*dtime/real(ninter)*zoliq(i)* &
                     (1.0-exp(-(zoliq(i)/zneb(i)/zcl(i))**2))*(1.-zfice(i))
                ! m\^eme chose pour la glace
                IF (ptconv(i, k)) THEN
                   zfroi(i) = dtime/real(ninter)/zdz(i)*zoliq(i)* &
                        fallvc(zrhol(i))*zfice(i)
                ELSE
                   zfroi(i) = dtime/real(ninter)/zdz(i)*zoliq(i)* &
                        fallvs(zrhol(i))*zfice(i)
                END IF
                ztot(i) = zchau(i) + zfroi(i)
                IF (zneb(i)==seuil_neb) ztot(i) = 0.0
                ztot(i) = min(max(ztot(i), 0.0), zoliq(i))
                zoliq(i) = max(zoliq(i)-ztot(i), 0.0)
                radliq(i, k) = radliq(i, k) + zoliq(i)/real(ninter+1)
             END IF
          END DO
       END DO

       DO i = 1, klon
          IF (rneb(i, k)>0.0) THEN
             d_ql(i, k) = zoliq(i)
             zrfl(i) = zrfl(i) + max(zcond(i) - zoliq(i), 0.) &
                  * (paprs(i, k) - paprs(i, k + 1)) / (rg * dtime)
          END IF
          IF (zt(i)<rtt) THEN
             psfl(i, k) = zrfl(i)
          ELSE
             prfl(i, k) = zrfl(i)
          END IF
       END DO

       ! Calculer les tendances de q et de t :
       DO i = 1, klon
          d_q(i, k) = zq(i) - q(i, k)
          d_t(i, k) = zt(i) - t(i, k)
       END DO

       ! Calcul du lessivage stratiforme
       DO i = 1, klon
          zprec_cond(i) = max(zcond(i) - zoliq(i), 0.0) &
               * (paprs(i, k)-paprs(i, k+1))/rg
          IF (rneb(i, k)>0.0 .AND. zprec_cond(i)>0.) THEN
             ! lessivage nucleation LMD5 dans la couche elle-meme
             IF (t(i, k)>=ztglace) THEN
                zalpha_tr = a_tr_sca(3)
             ELSE
                zalpha_tr = a_tr_sca(4)
             END IF
             zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))
             pfrac_nucl(i, k) = pfrac_nucl(i, k)*(1.-zneb(i)*zfrac_lessi)
             frac_nucl(i, k) = 1. - zneb(i)*zfrac_lessi

             ! nucleation avec un facteur -1 au lieu de -0.5
             zfrac_lessi = 1. - exp(-zprec_cond(i)/zneb(i))
             pfrac_1nucl(i, k) = pfrac_1nucl(i, k)*(1.-zneb(i)*zfrac_lessi)
          END IF
       END DO

       ! Lessivage par impaction dans les couches en-dessous
       ! boucle sur i 
       DO kk = k - 1, 1, -1
          DO i = 1, klon
             IF (rneb(i, k)>0. .AND. zprec_cond(i)>0.) THEN
                IF (t(i, kk)>=ztglace) THEN
                   zalpha_tr = a_tr_sca(1)
                ELSE
                   zalpha_tr = a_tr_sca(2)
                END IF
                zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))
                pfrac_impa(i, kk) = pfrac_impa(i, kk)*(1.-zneb(i)*zfrac_lessi)
                frac_impa(i, kk) = 1. - zneb(i)*zfrac_lessi
             END IF
          END DO
       END DO
    end DO loop_vertical

    ! Pluie ou neige au sol selon la temperature de la 1ere couche

    DO i = 1, klon
       IF ((t(i, 1)+d_t(i, 1))<rtt) THEN
          snow(i) = zrfl(i)
          zlh_solid(i) = rlstt - rlvtt
       ELSE
          rain(i) = zrfl(i)
          zlh_solid(i) = 0.
       END IF
    END DO

    ! For energy conservation: when snow is present, the solification
    ! latent heat is considered.
    DO k = 1, klev
       DO i = 1, klon
          zcpair = rcpd*(1.0+rvtmp2*(q(i, k)+d_q(i, k)))
          zmair = (paprs(i, k)-paprs(i, k+1))/rg
          zm_solid = (prfl(i, k)-prfl(i, k+1)+psfl(i, k)-psfl(i, k+1))*dtime
          d_t(i, k) = d_t(i, k) + zlh_solid(i)*zm_solid/(zcpair*zmair)
       END DO
    END DO

  contains

    ! vitesse de chute pour cristaux de glace

    REAL function fallvs(zzz)
      REAL, intent(in):: zzz
      fallvs = 3.29/2.0*((zzz)**0.16)*ffallv_lsc
    end function fallvs

    !********************************************************

    real function fallvc(zzz)
      REAL, intent(in):: zzz
      fallvc = 3.29/2.0*((zzz)**0.16)*ffallv_con
    end function fallvc

  END SUBROUTINE fisrtilp

end module fisrtilp_m
1	guez	68	module fisrtilp_m
2	guez	3
3	guez	68	IMPLICIT NONE
4	guez	3
5	guez	68	contains
6	guez	3
7	guez	68	SUBROUTINE fisrtilp(dtime, paprs, pplay, t, q, ptconv, ratqs, d_t, d_q, &
8			d_ql, rneb, radliq, rain, snow, pfrac_impa, pfrac_nucl, pfrac_1nucl, &
9			frac_impa, frac_nucl, prfl, psfl, rhcl)
10	guez	3
11	guez	213	! From phylmd/fisrtilp.F, version 1.2, 2004/11/09 16:55:40
12	guez	78	! First author: Z. X. Li (LMD/CNRS), 20 mars 1995
13	guez	3
14	guez	207	! Objet : condensation et pr\'ecipitation stratiforme, sch\'ema de
15			! nuage, sch\'ema de condensation \`a grande \'echelle (pluie).
16	guez	3
17	guez	78	USE comfisrtilp, ONLY: cld_lc_con, cld_lc_lsc, cld_tau_con, &
18			cld_tau_lsc, coef_eva, ffallv_con, ffallv_lsc, iflag_pdf, reevap_ice
19	guez	68	USE dimphy, ONLY: klev, klon
20	guez	221	USE fcttre, ONLY: foede, foeew
21	guez	78	USE numer_rec_95, ONLY: nr_erf
22	guez	68	USE suphec_m, ONLY: rcpd, rd, retv, rg, rlstt, rlvtt, rtt
23			USE yoethf_m, ONLY: r2es, r5ies, r5les, rvtmp2
24	guez	3
25	guez	78	REAL, INTENT (IN):: dtime ! intervalle du temps (s)
26			REAL, INTENT (IN):: paprs(klon, klev+1) ! pression a inter-couche
27	guez	68	REAL, INTENT (IN):: pplay(klon, klev) ! pression au milieu de couche
28			REAL, INTENT (IN):: t(klon, klev) ! temperature (K)
29	guez	73	REAL, INTENT (IN):: q(klon, klev) ! humidite specifique (kg/kg)
30	guez	78	LOGICAL, INTENT (IN):: ptconv(klon, klev)
31	guez	73
32	guez	78	REAL, INTENT (IN):: ratqs(klon, klev)
33			! determine la largeur de distribution de vapeur
34	guez	3
35	guez	78	REAL, INTENT (out):: d_t(klon, klev) ! incrementation de la temperature (K)
36			REAL, INTENT (out):: d_q(klon, klev) ! incrementation de la vapeur d'eau
37			REAL, INTENT (out):: d_ql(klon, klev) ! incrementation de l'eau liquide
38			REAL, INTENT (out):: rneb(klon, klev) ! fraction nuageuse
39	guez	3
40	guez	78	REAL, INTENT (out):: radliq(klon, klev)
41			! eau liquide utilisee dans rayonnement
42	guez	73
43	guez	78	REAL, INTENT (out):: rain(klon) ! pluies (mm/s)
44			REAL, INTENT (out):: snow(klon) ! neige (mm/s)
45
46			! Coeffients de fraction lessivee :
47			REAL, INTENT (inout):: pfrac_impa(klon, klev)
48			REAL, INTENT (inout):: pfrac_nucl(klon, klev)
49			REAL, INTENT (inout):: pfrac_1nucl(klon, klev)
50
51			! Fraction d'aerosols lessivee par impaction
52			REAL, INTENT (out):: frac_impa(klon, klev)
53
54			! Fraction d'aerosols lessivee par nucleation
55			REAL, INTENT (out):: frac_nucl(klon, klev)
56
57			REAL, INTENT (out):: prfl(klon, klev+1)
58			! flux d'eau precipitante aux interfaces (kg/m2/s)
59
60			REAL, INTENT (out):: psfl(klon, klev+1)
61			! flux d'eau precipitante aux interfaces (kg/m2/s)
62
63			REAL, INTENT (out):: rhcl(klon, klev) ! humidite relative en ciel clair
64
65	guez	73	! Local:
66
67	guez	68	REAL zct(klon), zcl(klon)
68	guez	3
69	guez	68	! Options du programme:
70	guez	22
71	guez	68	REAL seuil_neb ! un nuage existe vraiment au-dela
72			PARAMETER (seuil_neb=0.001)
73	guez	22
74	guez	68	INTEGER ninter ! sous-intervals pour la precipitation
75			PARAMETER (ninter=5)
76	guez	78	LOGICAL evap_prec ! evaporation de la pluie
77	guez	68	PARAMETER (evap_prec=.TRUE.)
78			REAL zpdf_sig(klon), zpdf_k(klon), zpdf_delta(klon)
79			REAL zpdf_a(klon), zpdf_b(klon), zpdf_e1(klon), zpdf_e2(klon)
80	guez	22
81	guez	78	LOGICAL cpartiel ! condensation partielle
82	guez	68	PARAMETER (cpartiel=.TRUE.)
83			REAL t_coup
84			PARAMETER (t_coup=234.0)
85	guez	22
86	guez	68	INTEGER i, k, n, kk
87	guez	103	REAL zqs(klon), zdqs(klon), zcor, zcvm5
88			logical zdelta
89	guez	68	REAL zrfl(klon), zrfln(klon), zqev, zqevt
90			REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
91			REAL ztglace, zt(klon)
92	guez	78	INTEGER nexpo ! exponentiel pour glace/eau
93	guez	68	REAL zdz(klon), zrho(klon), ztot(klon), zrhol(klon)
94			REAL zchau(klon), zfroi(klon), zfice(klon), zneb(klon)
95	guez	22
96	guez	78	LOGICAL:: appel1er = .TRUE.
97	guez	22
98	guez	78	! Variables traceurs:
99			! Provisoire !!! Parametres alpha du lessivage
100			! A priori on a 4 scavenging numbers possibles
101	guez	22
102	guez	78	REAL, save:: a_tr_sca(4)
103	guez	22
104	guez	68	! Variables intermediaires
105	guez	22
106	guez	68	REAL zalpha_tr
107			REAL zfrac_lessi
108			REAL zprec_cond(klon)
109			REAL zmair, zcpair, zcpeau
110	guez	78	! Pour la conversion eau-neige
111	guez	68	REAL zlh_solid(klon), zm_solid
112	guez	22
113	guez	78	!---------------------------------------------------------------
114
115	guez	68	zdelq = 0.0
116	guez	22
117	guez	68	IF (appel1er) THEN
118			PRINT *, 'fisrtilp, ninter:', ninter
119			PRINT *, 'fisrtilp, evap_prec:', evap_prec
120			PRINT *, 'fisrtilp, cpartiel:', cpartiel
121			IF (abs(dtime / real(ninter) - 360.) > 0.001) THEN
122	guez	207	PRINT *, "fisrtilp : ce n'est pas pr\'evu, voir Z. X. Li", dtime
123			PRINT *, "Je pr\'ef\`ere un sous-intervalle de 6 minutes."
124	guez	68	END IF
125			appel1er = .FALSE.
126	guez	22
127	guez	78	! initialiation provisoire
128	guez	68	a_tr_sca(1) = -0.5
129			a_tr_sca(2) = -0.5
130			a_tr_sca(3) = -0.5
131			a_tr_sca(4) = -0.5
132	guez	22
133	guez	78	! Initialisation a 1 des coefs des fractions lessivees
134	guez	68	DO k = 1, klev
135			DO i = 1, klon
136			pfrac_nucl(i, k) = 1.
137			pfrac_1nucl(i, k) = 1.
138			pfrac_impa(i, k) = 1.
139			END DO
140			END DO
141	guez	78	END IF
142	guez	22
143	guez	78	! Initialisation a 0 de zoliq
144	guez	68	DO i = 1, klon
145			zoliq(i) = 0.
146			END DO
147			! Determiner les nuages froids par leur temperature
148	guez	78	! nexpo regle la raideur de la transition eau liquide / eau glace.
149	guez	22
150	guez	68	ztglace = rtt - 15.0
151			nexpo = 6
152	guez	22
153	guez	68	! Initialiser les sorties:
154	guez	22
155	guez	68	DO k = 1, klev + 1
156			DO i = 1, klon
157			prfl(i, k) = 0.0
158			psfl(i, k) = 0.0
159			END DO
160			END DO
161	guez	22
162	guez	68	DO k = 1, klev
163			DO i = 1, klon
164			d_t(i, k) = 0.0
165			d_q(i, k) = 0.0
166			d_ql(i, k) = 0.0
167			rneb(i, k) = 0.0
168			radliq(i, k) = 0.0
169			frac_nucl(i, k) = 1.
170			frac_impa(i, k) = 1.
171			END DO
172			END DO
173			DO i = 1, klon
174			rain(i) = 0.0
175			snow(i) = 0.0
176			END DO
177	guez	22
178	guez	68	! Initialiser le flux de precipitation a zero
179	guez	22
180	guez	68	DO i = 1, klon
181			zrfl(i) = 0.0
182			zneb(i) = seuil_neb
183			END DO
184	guez	22
185	guez	78	! Pour plus de securite
186	guez	22
187	guez	68	zalpha_tr = 0.
188			zfrac_lessi = 0.
189	guez	22
190	guez	78	loop_vertical: DO k = klev, 1, -1
191	guez	68	DO i = 1, klon
192			zt(i) = t(i, k)
193			zq(i) = q(i, k)
194			END DO
195	guez	22
196	guez	68	! Calculer la varition de temp. de l'air du a la chaleur sensible
197			! transporter par la pluie.
198			! Il resterait a rajouter cet effet de la chaleur sensible sur les
199			! flux de surface, du a la diff. de temp. entre le 1er niveau et la
200			! surface.
201	guez	22
202	guez	68	DO i = 1, klon
203	guez	78	IF (k <= klev - 1) THEN
204	guez	68	zmair = (paprs(i, k)-paprs(i, k+1))/rg
205			zcpair = rcpd(1.0+rvtmp2zq(i))
206			zcpeau = rcpd*rvtmp2
207			zt(i) = ((t(i, k + 1) + d_t(i, k + 1)) * zrfl(i) * dtime &
208			* zcpeau + zmair * zcpair* zt(i)) &
209			/ (zmair * zcpair + zrfl(i) * dtime * zcpeau)
210			END IF
211			END DO
212	guez	22
213	guez	68	IF (evap_prec) THEN
214	guez	78	! Calculer l'evaporation de la precipitation
215	guez	68	DO i = 1, klon
216			IF (zrfl(i)>0.) THEN
217	guez	207	zqs(i) = r2es*foeew(zt(i), rtt >= zt(i))/pplay(i, k)
218			zqs(i) = min(0.5, zqs(i))
219			zcor = 1./(1.-retv*zqs(i))
220			zqs(i) = zqs(i)*zcor
221	guez	68	zqev = max(0.0, (zqs(i)-zq(i))*zneb(i))
222			zqevt = coef_eva(1.0-zq(i)/zqs(i))sqrt(zrfl(i))* &
223			(paprs(i, k)-paprs(i, k+1))/pplay(i, k)zt(i)rd/rg
224			zqevt = max(0.0, min(zqevt, zrfl(i)))rgdtime/ &
225			(paprs(i, k)-paprs(i, k+1))
226			zqev = min(zqev, zqevt)
227			zrfln(i) = zrfl(i) - zqev*(paprs(i, k)-paprs(i, k+1))/rg/dtime
228	guez	22
229	guez	207	! pour la glace, on r\'e\'evapore toute la pr\'ecip dans la
230	guez	68	! couche du dessous la glace venant de la couche du
231			! dessus est simplement dans la couche du dessous.
232	guez	22
233	guez	68	IF (zt(i)<t_coup .AND. reevap_ice) zrfln(i) = 0.
234	guez	3
235	guez	68	zq(i) = zq(i) - (zrfln(i)-zrfl(i))*(rg/(paprs(i, k)-paprs(i, &
236			k+1)))*dtime
237			zt(i) = zt(i) + (zrfln(i)-zrfl(i))*(rg/(paprs(i, k)-paprs(i, &
238			k+1)))dtimerlvtt/rcpd/(1.0+rvtmp2*zq(i))
239			zrfl(i) = zrfln(i)
240			END IF
241			END DO
242			END IF
243	guez	3
244	guez	68	! Calculer Qs et L/Cp*dQs/dT:
245	guez	3
246	guez	207	DO i = 1, klon
247			zdelta = rtt >= zt(i)
248			zcvm5 = merge(r5iesrlstt, r5lesrlvtt, zdelta)
249			zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zq(i))
250			zqs(i) = r2es*foeew(zt(i), zdelta)/pplay(i, k)
251			zqs(i) = min(0.5, zqs(i))
252			zcor = 1./(1.-retv*zqs(i))
253			zqs(i) = zqs(i)*zcor
254			zdqs(i) = foede(zt(i), zdelta, zcvm5, zqs(i), zcor)
255			END DO
256	guez	3
257	guez	68	! Determiner la condensation partielle et calculer la quantite
258			! de l'eau condensee:
259	guez	3
260	guez	68	IF (cpartiel) THEN
261	guez	78	! Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
262			! nuageuse a partir des PDF de Sandrine Bony.
263			! rneb : fraction nuageuse
264			! zqn : eau totale dans le nuage
265			! zcond : eau condensee moyenne dans la maille.
266	guez	3
267	guez	207	! on prend en compte le r\'echauffement qui diminue
268			! la partie condens\'ee
269	guez	3
270	guez	78	! Version avec les ratqs
271	guez	3
272	guez	68	IF (iflag_pdf==0) THEN
273			DO i = 1, klon
274			zdelq = min(ratqs(i, k), 0.99)*zq(i)
275			rneb(i, k) = (zq(i)+zdelq-zqs(i))/(2.0*zdelq)
276			zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
277			END DO
278			ELSE
279	guez	78	! Version avec les nouvelles PDFs.
280	guez	68	DO i = 1, klon
281	guez	78	IF (zq(i) < 1E-15) THEN
282			zq(i) = 1E-15
283	guez	68	END IF
284			END DO
285			DO i = 1, klon
286			zpdf_sig(i) = ratqs(i, k)*zq(i)
287			zpdf_k(i) = -sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
288			zpdf_delta(i) = log(zq(i)/zqs(i))
289			zpdf_a(i) = zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
290			zpdf_b(i) = zpdf_k(i)/(2.*sqrt(2.))
291			zpdf_e1(i) = zpdf_a(i) - zpdf_b(i)
292			zpdf_e1(i) = sign(min(abs(zpdf_e1(i)), 5.), zpdf_e1(i))
293			zpdf_e1(i) = 1. - nr_erf(zpdf_e1(i))
294			zpdf_e2(i) = zpdf_a(i) + zpdf_b(i)
295			zpdf_e2(i) = sign(min(abs(zpdf_e2(i)), 5.), zpdf_e2(i))
296			zpdf_e2(i) = 1. - nr_erf(zpdf_e2(i))
297			IF (zpdf_e1(i)<1.E-10) THEN
298			rneb(i, k) = 0.
299			zqn(i) = zqs(i)
300			ELSE
301			rneb(i, k) = 0.5*zpdf_e1(i)
302			zqn(i) = zq(i)*zpdf_e2(i)/zpdf_e1(i)
303			END IF
304			END DO
305	guez	78	END IF
306	guez	22
307	guez	68	DO i = 1, klon
308			IF (rneb(i, k)<=0.0) zqn(i) = 0.0
309			IF (rneb(i, k)>=1.0) zqn(i) = zq(i)
310	guez	78	rneb(i, k) = max(0., min(1., rneb(i, k)))
311	guez	207	! On ne divise pas par 1 + zdqs pour forcer \`a avoir l'eau
312			! pr\'edite par la convection. Attention : il va falloir
313	guez	78	! verifier tout ca.
314			zcond(i) = max(0., zqn(i)-zqs(i))*rneb(i, k)
315	guez	68	rhcl(i, k) = (zqs(i)+zq(i)-zdelq)/2./zqs(i)
316	guez	78	IF (rneb(i, k) <= 0.) rhcl(i, k) = zq(i) / zqs(i)
317			IF (rneb(i, k) >= 1.) rhcl(i, k) = 1.
318	guez	68	END DO
319			ELSE
320			DO i = 1, klon
321			IF (zq(i)>zqs(i)) THEN
322			rneb(i, k) = 1.0
323			ELSE
324			rneb(i, k) = 0.0
325			END IF
326			zcond(i) = max(0.0, zq(i)-zqs(i))/(1.+zdqs(i))
327			END DO
328			END IF
329	guez	22
330	guez	68	DO i = 1, klon
331			zq(i) = zq(i) - zcond(i)
332			zt(i) = zt(i) + zcond(i)rlvtt/rcpd/(1.0+rvtmp2zq(i))
333			END DO
334	guez	22
335	guez	68	! Partager l'eau condensee en precipitation et eau liquide nuageuse
336	guez	22
337	guez	68	DO i = 1, klon
338			IF (rneb(i, k)>0.0) THEN
339			zoliq(i) = zcond(i)
340			zrho(i) = pplay(i, k)/zt(i)/rd
341			zdz(i) = (paprs(i, k)-paprs(i, k+1))/(zrho(i)*rg)
342			zfice(i) = 1.0 - (zt(i)-ztglace)/(273.13-ztglace)
343			zfice(i) = min(max(zfice(i), 0.0), 1.0)
344			zfice(i) = zfice(i)**nexpo
345			zneb(i) = max(rneb(i, k), seuil_neb)
346			radliq(i, k) = zoliq(i)/real(ninter+1)
347			END IF
348			END DO
349	guez	22
350	guez	68	DO n = 1, ninter
351			DO i = 1, klon
352			IF (rneb(i, k)>0.0) THEN
353			zrhol(i) = zrho(i)*zoliq(i)/zneb(i)
354	guez	22
355	guez	68	IF (ptconv(i, k)) THEN
356			zcl(i) = cld_lc_con
357			zct(i) = 1./cld_tau_con
358			ELSE
359			zcl(i) = cld_lc_lsc
360			zct(i) = 1./cld_tau_lsc
361			END IF
362	guez	208	! quantit\'e d'eau \`a \'eliminer
363	guez	68	zchau(i) = zct(i)dtime/real(ninter)zoliq(i)* &
364			(1.0-exp(-(zoliq(i)/zneb(i)/zcl(i))*2))(1.-zfice(i))
365	guez	208	! m\^eme chose pour la glace
366	guez	68	IF (ptconv(i, k)) THEN
367			zfroi(i) = dtime/real(ninter)/zdz(i)zoliq(i) &
368			fallvc(zrhol(i))*zfice(i)
369			ELSE
370			zfroi(i) = dtime/real(ninter)/zdz(i)zoliq(i) &
371			fallvs(zrhol(i))*zfice(i)
372			END IF
373			ztot(i) = zchau(i) + zfroi(i)
374			IF (zneb(i)==seuil_neb) ztot(i) = 0.0
375			ztot(i) = min(max(ztot(i), 0.0), zoliq(i))
376			zoliq(i) = max(zoliq(i)-ztot(i), 0.0)
377			radliq(i, k) = radliq(i, k) + zoliq(i)/real(ninter+1)
378			END IF
379			END DO
380			END DO
381	guez	22
382	guez	68	DO i = 1, klon
383			IF (rneb(i, k)>0.0) THEN
384			d_ql(i, k) = zoliq(i)
385			zrfl(i) = zrfl(i) + max(zcond(i) - zoliq(i), 0.) &
386			* (paprs(i, k) - paprs(i, k + 1)) / (rg * dtime)
387			END IF
388			IF (zt(i)<rtt) THEN
389			psfl(i, k) = zrfl(i)
390			ELSE
391			prfl(i, k) = zrfl(i)
392			END IF
393			END DO
394	guez	22
395	guez	208	! Calculer les tendances de q et de t :
396	guez	68	DO i = 1, klon
397			d_q(i, k) = zq(i) - q(i, k)
398			d_t(i, k) = zt(i) - t(i, k)
399			END DO
400	guez	22
401	guez	78	! Calcul du lessivage stratiforme
402	guez	68	DO i = 1, klon
403	guez	213	zprec_cond(i) = max(zcond(i) - zoliq(i), 0.0) &
404			* (paprs(i, k)-paprs(i, k+1))/rg
405	guez	68	IF (rneb(i, k)>0.0 .AND. zprec_cond(i)>0.) THEN
406	guez	78	! lessivage nucleation LMD5 dans la couche elle-meme
407	guez	68	IF (t(i, k)>=ztglace) THEN
408			zalpha_tr = a_tr_sca(3)
409			ELSE
410			zalpha_tr = a_tr_sca(4)
411			END IF
412			zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))
413			pfrac_nucl(i, k) = pfrac_nucl(i, k)(1.-zneb(i)zfrac_lessi)
414			frac_nucl(i, k) = 1. - zneb(i)*zfrac_lessi
415	guez	22
416	guez	68	! nucleation avec un facteur -1 au lieu de -0.5
417			zfrac_lessi = 1. - exp(-zprec_cond(i)/zneb(i))
418			pfrac_1nucl(i, k) = pfrac_1nucl(i, k)(1.-zneb(i)zfrac_lessi)
419			END IF
420	guez	78	END DO
421	guez	22
422	guez	78	! Lessivage par impaction dans les couches en-dessous
423			! boucle sur i
424	guez	68	DO kk = k - 1, 1, -1
425			DO i = 1, klon
426	guez	78	IF (rneb(i, k)>0. .AND. zprec_cond(i)>0.) THEN
427	guez	68	IF (t(i, kk)>=ztglace) THEN
428			zalpha_tr = a_tr_sca(1)
429			ELSE
430			zalpha_tr = a_tr_sca(2)
431			END IF
432			zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))
433			pfrac_impa(i, kk) = pfrac_impa(i, kk)(1.-zneb(i)zfrac_lessi)
434			frac_impa(i, kk) = 1. - zneb(i)*zfrac_lessi
435			END IF
436			END DO
437			END DO
438	guez	78	end DO loop_vertical
439	guez	22
440	guez	68	! Pluie ou neige au sol selon la temperature de la 1ere couche
441
442			DO i = 1, klon
443			IF ((t(i, 1)+d_t(i, 1))<rtt) THEN
444			snow(i) = zrfl(i)
445			zlh_solid(i) = rlstt - rlvtt
446			ELSE
447			rain(i) = zrfl(i)
448			zlh_solid(i) = 0.
449			END IF
450			END DO
451
452			! For energy conservation: when snow is present, the solification
453			! latent heat is considered.
454			DO k = 1, klev
455			DO i = 1, klon
456			zcpair = rcpd(1.0+rvtmp2(q(i, k)+d_q(i, k)))
457			zmair = (paprs(i, k)-paprs(i, k+1))/rg
458			zm_solid = (prfl(i, k)-prfl(i, k+1)+psfl(i, k)-psfl(i, k+1))*dtime
459			d_t(i, k) = d_t(i, k) + zlh_solid(i)zm_solid/(zcpairzmair)
460			END DO
461			END DO
462
463	guez	108	contains
464
465	guez	207	! vitesse de chute pour cristaux de glace
466	guez	108
467			REAL function fallvs(zzz)
468	guez	207	REAL, intent(in):: zzz
469	guez	108	fallvs = 3.29/2.0((zzz)0.16)ffallv_lsc
470			end function fallvs
471
472	guez	207	!********************************************************
473
474	guez	108	real function fallvc(zzz)
475	guez	207	REAL, intent(in):: zzz
476	guez	108	fallvc = 3.29/2.0((zzz)0.16)ffallv_con
477			end function fallvc
478
479	guez	68	END SUBROUTINE fisrtilp
480
481			end module fisrtilp_m