trunk/phylmd/fisrtilp.f

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
    ! Other authors: Olivier, AA, IM, YM, MAF

    ! Objet : condensation et précipitation stratiforme, schéma de
    ! nuage, schéma de condensation à grande échelle (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: dqsatl, dqsats, foede, foeew, qsatl, qsats, thermcep
    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

    ! Arguments:

    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évu, voir Z. X. Li", dtime
          PRINT *, 'Je préfère 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
                IF (thermcep) 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
                ELSE
                   IF (zt(i)<t_coup) THEN
                      zqs(i) = qsats(zt(i))/pplay(i, k)
                   ELSE
                      zqs(i) = qsatl(zt(i))/pplay(i, k)
                   END IF
                END IF
                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éévapore toute la précip 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:

       IF (thermcep) THEN
          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
       ELSE
          DO i = 1, klon
             IF (zt(i)<t_coup) THEN
                zqs(i) = qsats(zt(i))/pplay(i, k)
                zdqs(i) = dqsats(zt(i), zqs(i))
             ELSE
                zqs(i) = qsatl(zt(i))/pplay(i, k)
                zdqs(i) = dqsatl(zt(i), zqs(i))
             END IF
          END DO
       END IF

       ! 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échauffement qui diminue
          ! la partie condensée

          ! 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 à avoir l'eau
             ! prédite 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é d'eau à élminier.
                zchau(i) = zct(i)*dtime/real(ninter)*zoliq(i)* &
                     (1.0-exp(-(zoliq(i)/zneb(i)/zcl(i))**2))*(1.-zfice(i))
                ! meme 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 crystaux de glace

    REAL function fallvs(zzz)
      REAL zzz
      fallvs = 3.29/2.0*((zzz)**0.16)*ffallv_lsc
    end function fallvs

    real function fallvc(zzz)
      REAL zzz
      fallvc = 3.29/2.0*((zzz)**0.16)*ffallv_con
    end function fallvc

  END SUBROUTINE fisrtilp

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