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

    ! Fonctions en ligne:

    REAL fallvs, fallvc ! vitesse de chute pour crystaux de glace 
    REAL zzz

    fallvc(zzz) = 3.29/2.0*((zzz)**0.16)*ffallv_con
    fallvs(zzz) = 3.29/2.0*((zzz)**0.16)*ffallv_lsc

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

    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

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