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), zdelta, zcor, zcvm5
    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
                   zdelta = max(0., sign(1., rtt-zt(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
                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 = max(0., sign(1., rtt-zt(i)))
             zcvm5 = r5les*rlvtt*(1.-zdelta) + r5ies*rlstt*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			REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
91			REAL zrfl(klon), zrfln(klon), zqev, zqevt
92			REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
93			REAL ztglace, zt(klon)
94	guez	78	INTEGER nexpo ! exponentiel pour glace/eau
95	guez	68	REAL zdz(klon), zrho(klon), ztot(klon), zrhol(klon)
96			REAL zchau(klon), zfroi(klon), zfice(klon), zneb(klon)
97	guez	22
98	guez	78	LOGICAL:: appel1er = .TRUE.
99	guez	22
100	guez	78	! Variables traceurs:
101			! Provisoire !!! Parametres alpha du lessivage
102			! A priori on a 4 scavenging numbers possibles
103	guez	22
104	guez	78	REAL, save:: a_tr_sca(4)
105	guez	22
106	guez	68	! Variables intermediaires
107	guez	22
108	guez	68	REAL zalpha_tr
109			REAL zfrac_lessi
110			REAL zprec_cond(klon)
111			REAL zmair, zcpair, zcpeau
112	guez	78	! Pour la conversion eau-neige
113	guez	68	REAL zlh_solid(klon), zm_solid
114	guez	22
115	guez	68	! Fonctions en ligne:
116	guez	22
117	guez	78	REAL fallvs, fallvc ! vitesse de chute pour crystaux de glace
118	guez	68	REAL zzz
119	guez	22
120	guez	68	fallvc(zzz) = 3.29/2.0((zzz)0.16)ffallv_con
121			fallvs(zzz) = 3.29/2.0((zzz)0.16)ffallv_lsc
122	guez	22
123	guez	78	!---------------------------------------------------------------
124
125	guez	68	zdelq = 0.0
126	guez	22
127	guez	68	IF (appel1er) THEN
128			PRINT *, 'fisrtilp, ninter:', ninter
129			PRINT *, 'fisrtilp, evap_prec:', evap_prec
130			PRINT *, 'fisrtilp, cpartiel:', cpartiel
131			IF (abs(dtime / real(ninter) - 360.) > 0.001) THEN
132			PRINT *, "fisrtilp : ce n'est pas prévu, voir Z. X. Li", dtime
133			PRINT *, 'Je préfère un sous-intervalle de 6 minutes.'
134			END IF
135			appel1er = .FALSE.
136	guez	22
137	guez	78	! initialiation provisoire
138	guez	68	a_tr_sca(1) = -0.5
139			a_tr_sca(2) = -0.5
140			a_tr_sca(3) = -0.5
141			a_tr_sca(4) = -0.5
142	guez	22
143	guez	78	! Initialisation a 1 des coefs des fractions lessivees
144	guez	68	DO k = 1, klev
145			DO i = 1, klon
146			pfrac_nucl(i, k) = 1.
147			pfrac_1nucl(i, k) = 1.
148			pfrac_impa(i, k) = 1.
149			END DO
150			END DO
151	guez	78	END IF
152	guez	22
153	guez	78	! Initialisation a 0 de zoliq
154	guez	68	DO i = 1, klon
155			zoliq(i) = 0.
156			END DO
157			! Determiner les nuages froids par leur temperature
158	guez	78	! nexpo regle la raideur de la transition eau liquide / eau glace.
159	guez	22
160	guez	68	ztglace = rtt - 15.0
161			nexpo = 6
162	guez	22
163	guez	68	! Initialiser les sorties:
164	guez	22
165	guez	68	DO k = 1, klev + 1
166			DO i = 1, klon
167			prfl(i, k) = 0.0
168			psfl(i, k) = 0.0
169			END DO
170			END DO
171	guez	22
172	guez	68	DO k = 1, klev
173			DO i = 1, klon
174			d_t(i, k) = 0.0
175			d_q(i, k) = 0.0
176			d_ql(i, k) = 0.0
177			rneb(i, k) = 0.0
178			radliq(i, k) = 0.0
179			frac_nucl(i, k) = 1.
180			frac_impa(i, k) = 1.
181			END DO
182			END DO
183			DO i = 1, klon
184			rain(i) = 0.0
185			snow(i) = 0.0
186			END DO
187	guez	22
188	guez	68	! Initialiser le flux de precipitation a zero
189	guez	22
190	guez	68	DO i = 1, klon
191			zrfl(i) = 0.0
192			zneb(i) = seuil_neb
193			END DO
194	guez	22
195	guez	78	! Pour plus de securite
196	guez	22
197	guez	68	zalpha_tr = 0.
198			zfrac_lessi = 0.
199	guez	22
200	guez	78	loop_vertical: DO k = klev, 1, -1
201	guez	68	DO i = 1, klon
202			zt(i) = t(i, k)
203			zq(i) = q(i, k)
204			END DO
205	guez	22
206	guez	68	! Calculer la varition de temp. de l'air du a la chaleur sensible
207			! transporter par la pluie.
208			! Il resterait a rajouter cet effet de la chaleur sensible sur les
209			! flux de surface, du a la diff. de temp. entre le 1er niveau et la
210			! surface.
211	guez	22
212	guez	68	DO i = 1, klon
213	guez	78	IF (k <= klev - 1) THEN
214	guez	68	zmair = (paprs(i, k)-paprs(i, k+1))/rg
215			zcpair = rcpd(1.0+rvtmp2zq(i))
216			zcpeau = rcpd*rvtmp2
217			zt(i) = ((t(i, k + 1) + d_t(i, k + 1)) * zrfl(i) * dtime &
218			* zcpeau + zmair * zcpair* zt(i)) &
219			/ (zmair * zcpair + zrfl(i) * dtime * zcpeau)
220			END IF
221			END DO
222	guez	22
223	guez	68	IF (evap_prec) THEN
224	guez	78	! Calculer l'evaporation de la precipitation
225	guez	68	DO i = 1, klon
226			IF (zrfl(i)>0.) THEN
227			IF (thermcep) THEN
228			zdelta = max(0., sign(1., rtt-zt(i)))
229			zqs(i) = r2es*foeew(zt(i), zdelta)/pplay(i, k)
230			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			zdelta = max(0., sign(1., rtt-zt(i)))
268			zcvm5 = r5lesrlvtt(1.-zdelta) + r5iesrlsttzdelta
269			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