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\'ecipitation stratiforme, sch\'ema de
    ! nuage, sch\'ema de condensation \`a grande \'echelle (pluie).

    USE comfisrtilp, ONLY: cld_lc_con, cld_lc_lsc, cld_tau_con, &
         cld_tau_lsc, coef_eva, ffallv_con, ffallv_lsc, iflag_pdf, reevap_ice
    USE dimphy, ONLY: klev, klon
    USE fcttre, ONLY: dqsatl, dqsats, foede, foeew, qsatl, qsats
    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\'evu, voir Z. X. Li", dtime
          PRINT *, "Je pr\'ef\`ere un sous-intervalle de 6 minutes."
       END IF
       appel1er = .FALSE.

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

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

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

    ztglace = rtt - 15.0
    nexpo = 6

    ! Initialiser les sorties:

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

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

    ! Initialiser le flux de precipitation a zero

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

    ! Pour plus de securite

    zalpha_tr = 0.
    zfrac_lessi = 0.

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

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

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

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

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

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

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

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

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

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

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

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

          ! Version avec les ratqs

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  contains

    ! vitesse de chute pour cristaux de glace

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

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

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

  END SUBROUTINE fisrtilp

end module fisrtilp_m
1	guez	68	module fisrtilp_m
2	guez	3
3	guez	68	IMPLICIT NONE
4	guez	3
5	guez	68	contains
6	guez	3
7	guez	68	SUBROUTINE fisrtilp(dtime, paprs, pplay, t, q, ptconv, ratqs, d_t, d_q, &
8			d_ql, rneb, radliq, rain, snow, pfrac_impa, pfrac_nucl, pfrac_1nucl, &
9			frac_impa, frac_nucl, prfl, psfl, rhcl)
10	guez	3
11	guez	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	207	! Objet : condensation et pr\'ecipitation stratiforme, sch\'ema de
16			! nuage, sch\'ema de condensation \`a grande \'echelle (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	207	USE fcttre, ONLY: dqsatl, dqsats, foede, foeew, qsatl, qsats
22	guez	78	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	78	!---------------------------------------------------------------
117
118	guez	68	zdelq = 0.0
119	guez	22
120	guez	68	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	guez	207	PRINT *, "fisrtilp : ce n'est pas pr\'evu, voir Z. X. Li", dtime
126			PRINT *, "Je pr\'ef\`ere un sous-intervalle de 6 minutes."
127	guez	68	END IF
128			appel1er = .FALSE.
129	guez	22
130	guez	78	! initialiation provisoire
131	guez	68	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	guez	22
136	guez	78	! Initialisation a 1 des coefs des fractions lessivees
137	guez	68	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	guez	78	END IF
145	guez	22
146	guez	78	! Initialisation a 0 de zoliq
147	guez	68	DO i = 1, klon
148			zoliq(i) = 0.
149			END DO
150			! Determiner les nuages froids par leur temperature
151	guez	78	! nexpo regle la raideur de la transition eau liquide / eau glace.
152	guez	22
153	guez	68	ztglace = rtt - 15.0
154			nexpo = 6
155	guez	22
156	guez	68	! Initialiser les sorties:
157	guez	22
158	guez	68	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	guez	22
165	guez	68	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	guez	22
181	guez	68	! Initialiser le flux de precipitation a zero
182	guez	22
183	guez	68	DO i = 1, klon
184			zrfl(i) = 0.0
185			zneb(i) = seuil_neb
186			END DO
187	guez	22
188	guez	78	! Pour plus de securite
189	guez	22
190	guez	68	zalpha_tr = 0.
191			zfrac_lessi = 0.
192	guez	22
193	guez	78	loop_vertical: DO k = klev, 1, -1
194	guez	68	DO i = 1, klon
195			zt(i) = t(i, k)
196			zq(i) = q(i, k)
197			END DO
198	guez	22
199	guez	68	! 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	guez	22
205	guez	68	DO i = 1, klon
206	guez	78	IF (k <= klev - 1) THEN
207	guez	68	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	guez	22
216	guez	68	IF (evap_prec) THEN
217	guez	78	! Calculer l'evaporation de la precipitation
218	guez	68	DO i = 1, klon
219			IF (zrfl(i)>0.) THEN
220	guez	207	zqs(i) = r2es*foeew(zt(i), rtt >= zt(i))/pplay(i, k)
221			zqs(i) = min(0.5, zqs(i))
222			zcor = 1./(1.-retv*zqs(i))
223			zqs(i) = zqs(i)*zcor
224	guez	68	zqev = max(0.0, (zqs(i)-zq(i))*zneb(i))
225			zqevt = coef_eva(1.0-zq(i)/zqs(i))sqrt(zrfl(i))* &
226			(paprs(i, k)-paprs(i, k+1))/pplay(i, k)zt(i)rd/rg
227			zqevt = max(0.0, min(zqevt, zrfl(i)))rgdtime/ &
228			(paprs(i, k)-paprs(i, k+1))
229			zqev = min(zqev, zqevt)
230			zrfln(i) = zrfl(i) - zqev*(paprs(i, k)-paprs(i, k+1))/rg/dtime
231	guez	22
232	guez	207	! pour la glace, on r\'e\'evapore toute la pr\'ecip dans la
233	guez	68	! couche du dessous la glace venant de la couche du
234			! dessus est simplement dans la couche du dessous.
235	guez	22
236	guez	68	IF (zt(i)<t_coup .AND. reevap_ice) zrfln(i) = 0.
237	guez	3
238	guez	68	zq(i) = zq(i) - (zrfln(i)-zrfl(i))*(rg/(paprs(i, k)-paprs(i, &
239			k+1)))*dtime
240			zt(i) = zt(i) + (zrfln(i)-zrfl(i))*(rg/(paprs(i, k)-paprs(i, &
241			k+1)))dtimerlvtt/rcpd/(1.0+rvtmp2*zq(i))
242			zrfl(i) = zrfln(i)
243			END IF
244			END DO
245			END IF
246	guez	3
247	guez	68	! Calculer Qs et L/Cp*dQs/dT:
248	guez	3
249	guez	207	DO i = 1, klon
250			zdelta = rtt >= zt(i)
251			zcvm5 = merge(r5iesrlstt, r5lesrlvtt, zdelta)
252			zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zq(i))
253			zqs(i) = r2es*foeew(zt(i), zdelta)/pplay(i, k)
254			zqs(i) = min(0.5, zqs(i))
255			zcor = 1./(1.-retv*zqs(i))
256			zqs(i) = zqs(i)*zcor
257			zdqs(i) = foede(zt(i), zdelta, zcvm5, zqs(i), zcor)
258			END DO
259	guez	3
260	guez	68	! Determiner la condensation partielle et calculer la quantite
261			! de l'eau condensee:
262	guez	3
263	guez	68	IF (cpartiel) THEN
264	guez	78	! Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
265			! nuageuse a partir des PDF de Sandrine Bony.
266			! rneb : fraction nuageuse
267			! zqn : eau totale dans le nuage
268			! zcond : eau condensee moyenne dans la maille.
269	guez	3
270	guez	207	! on prend en compte le r\'echauffement qui diminue
271			! la partie condens\'ee
272	guez	3
273	guez	78	! Version avec les ratqs
274	guez	3
275	guez	68	IF (iflag_pdf==0) THEN
276			DO i = 1, klon
277			zdelq = min(ratqs(i, k), 0.99)*zq(i)
278			rneb(i, k) = (zq(i)+zdelq-zqs(i))/(2.0*zdelq)
279			zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
280			END DO
281			ELSE
282	guez	78	! Version avec les nouvelles PDFs.
283	guez	68	DO i = 1, klon
284	guez	78	IF (zq(i) < 1E-15) THEN
285			zq(i) = 1E-15
286	guez	68	END IF
287			END DO
288			DO i = 1, klon
289			zpdf_sig(i) = ratqs(i, k)*zq(i)
290			zpdf_k(i) = -sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
291			zpdf_delta(i) = log(zq(i)/zqs(i))
292			zpdf_a(i) = zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
293			zpdf_b(i) = zpdf_k(i)/(2.*sqrt(2.))
294			zpdf_e1(i) = zpdf_a(i) - zpdf_b(i)
295			zpdf_e1(i) = sign(min(abs(zpdf_e1(i)), 5.), zpdf_e1(i))
296			zpdf_e1(i) = 1. - nr_erf(zpdf_e1(i))
297			zpdf_e2(i) = zpdf_a(i) + zpdf_b(i)
298			zpdf_e2(i) = sign(min(abs(zpdf_e2(i)), 5.), zpdf_e2(i))
299			zpdf_e2(i) = 1. - nr_erf(zpdf_e2(i))
300			IF (zpdf_e1(i)<1.E-10) THEN
301			rneb(i, k) = 0.
302			zqn(i) = zqs(i)
303			ELSE
304			rneb(i, k) = 0.5*zpdf_e1(i)
305			zqn(i) = zq(i)*zpdf_e2(i)/zpdf_e1(i)
306			END IF
307			END DO
308	guez	78	END IF
309	guez	22
310	guez	68	DO i = 1, klon
311			IF (rneb(i, k)<=0.0) zqn(i) = 0.0
312			IF (rneb(i, k)>=1.0) zqn(i) = zq(i)
313	guez	78	rneb(i, k) = max(0., min(1., rneb(i, k)))
314	guez	207	! On ne divise pas par 1 + zdqs pour forcer \`a avoir l'eau
315			! pr\'edite par la convection. Attention : il va falloir
316	guez	78	! verifier tout ca.
317			zcond(i) = max(0., zqn(i)-zqs(i))*rneb(i, k)
318	guez	68	rhcl(i, k) = (zqs(i)+zq(i)-zdelq)/2./zqs(i)
319	guez	78	IF (rneb(i, k) <= 0.) rhcl(i, k) = zq(i) / zqs(i)
320			IF (rneb(i, k) >= 1.) rhcl(i, k) = 1.
321	guez	68	END DO
322			ELSE
323			DO i = 1, klon
324			IF (zq(i)>zqs(i)) THEN
325			rneb(i, k) = 1.0
326			ELSE
327			rneb(i, k) = 0.0
328			END IF
329			zcond(i) = max(0.0, zq(i)-zqs(i))/(1.+zdqs(i))
330			END DO
331			END IF
332	guez	22
333	guez	68	DO i = 1, klon
334			zq(i) = zq(i) - zcond(i)
335			zt(i) = zt(i) + zcond(i)rlvtt/rcpd/(1.0+rvtmp2zq(i))
336			END DO
337	guez	22
338	guez	68	! Partager l'eau condensee en precipitation et eau liquide nuageuse
339	guez	22
340	guez	68	DO i = 1, klon
341			IF (rneb(i, k)>0.0) THEN
342			zoliq(i) = zcond(i)
343			zrho(i) = pplay(i, k)/zt(i)/rd
344			zdz(i) = (paprs(i, k)-paprs(i, k+1))/(zrho(i)*rg)
345			zfice(i) = 1.0 - (zt(i)-ztglace)/(273.13-ztglace)
346			zfice(i) = min(max(zfice(i), 0.0), 1.0)
347			zfice(i) = zfice(i)**nexpo
348			zneb(i) = max(rneb(i, k), seuil_neb)
349			radliq(i, k) = zoliq(i)/real(ninter+1)
350			END IF
351			END DO
352	guez	22
353	guez	68	DO n = 1, ninter
354			DO i = 1, klon
355			IF (rneb(i, k)>0.0) THEN
356			zrhol(i) = zrho(i)*zoliq(i)/zneb(i)
357	guez	22
358	guez	68	IF (ptconv(i, k)) THEN
359			zcl(i) = cld_lc_con
360			zct(i) = 1./cld_tau_con
361			ELSE
362			zcl(i) = cld_lc_lsc
363			zct(i) = 1./cld_tau_lsc
364			END IF
365	guez	208	! quantit\'e d'eau \`a \'eliminer
366	guez	68	zchau(i) = zct(i)dtime/real(ninter)zoliq(i)* &
367			(1.0-exp(-(zoliq(i)/zneb(i)/zcl(i))*2))(1.-zfice(i))
368	guez	208	! m\^eme chose pour la glace
369	guez	68	IF (ptconv(i, k)) THEN
370			zfroi(i) = dtime/real(ninter)/zdz(i)zoliq(i) &
371			fallvc(zrhol(i))*zfice(i)
372			ELSE
373			zfroi(i) = dtime/real(ninter)/zdz(i)zoliq(i) &
374			fallvs(zrhol(i))*zfice(i)
375			END IF
376			ztot(i) = zchau(i) + zfroi(i)
377			IF (zneb(i)==seuil_neb) ztot(i) = 0.0
378			ztot(i) = min(max(ztot(i), 0.0), zoliq(i))
379			zoliq(i) = max(zoliq(i)-ztot(i), 0.0)
380			radliq(i, k) = radliq(i, k) + zoliq(i)/real(ninter+1)
381			END IF
382			END DO
383			END DO
384	guez	22
385	guez	68	DO i = 1, klon
386			IF (rneb(i, k)>0.0) THEN
387			d_ql(i, k) = zoliq(i)
388			zrfl(i) = zrfl(i) + max(zcond(i) - zoliq(i), 0.) &
389			* (paprs(i, k) - paprs(i, k + 1)) / (rg * dtime)
390			END IF
391			IF (zt(i)<rtt) THEN
392			psfl(i, k) = zrfl(i)
393			ELSE
394			prfl(i, k) = zrfl(i)
395			END IF
396			END DO
397	guez	22
398	guez	208	! Calculer les tendances de q et de t :
399	guez	68	DO i = 1, klon
400			d_q(i, k) = zq(i) - q(i, k)
401			d_t(i, k) = zt(i) - t(i, k)
402			END DO
403	guez	22
404	guez	78	! Calcul du lessivage stratiforme
405	guez	68	DO i = 1, klon
406			zprec_cond(i) = max(zcond(i)-zoliq(i), 0.0)* &
407			(paprs(i, k)-paprs(i, k+1))/rg
408			IF (rneb(i, k)>0.0 .AND. zprec_cond(i)>0.) THEN
409	guez	78	! lessivage nucleation LMD5 dans la couche elle-meme
410	guez	68	IF (t(i, k)>=ztglace) THEN
411			zalpha_tr = a_tr_sca(3)
412			ELSE
413			zalpha_tr = a_tr_sca(4)
414			END IF
415			zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))
416			pfrac_nucl(i, k) = pfrac_nucl(i, k)(1.-zneb(i)zfrac_lessi)
417			frac_nucl(i, k) = 1. - zneb(i)*zfrac_lessi
418	guez	22
419	guez	68	! nucleation avec un facteur -1 au lieu de -0.5
420			zfrac_lessi = 1. - exp(-zprec_cond(i)/zneb(i))
421			pfrac_1nucl(i, k) = pfrac_1nucl(i, k)(1.-zneb(i)zfrac_lessi)
422			END IF
423	guez	78	END DO
424	guez	22
425	guez	78	! Lessivage par impaction dans les couches en-dessous
426			! boucle sur i
427	guez	68	DO kk = k - 1, 1, -1
428			DO i = 1, klon
429	guez	78	IF (rneb(i, k)>0. .AND. zprec_cond(i)>0.) THEN
430	guez	68	IF (t(i, kk)>=ztglace) THEN
431			zalpha_tr = a_tr_sca(1)
432			ELSE
433			zalpha_tr = a_tr_sca(2)
434			END IF
435			zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))
436			pfrac_impa(i, kk) = pfrac_impa(i, kk)(1.-zneb(i)zfrac_lessi)
437			frac_impa(i, kk) = 1. - zneb(i)*zfrac_lessi
438			END IF
439			END DO
440			END DO
441	guez	78	end DO loop_vertical
442	guez	22
443	guez	68	! Pluie ou neige au sol selon la temperature de la 1ere couche
444
445			DO i = 1, klon
446			IF ((t(i, 1)+d_t(i, 1))<rtt) THEN
447			snow(i) = zrfl(i)
448			zlh_solid(i) = rlstt - rlvtt
449			ELSE
450			rain(i) = zrfl(i)
451			zlh_solid(i) = 0.
452			END IF
453			END DO
454
455			! For energy conservation: when snow is present, the solification
456			! latent heat is considered.
457			DO k = 1, klev
458			DO i = 1, klon
459			zcpair = rcpd(1.0+rvtmp2(q(i, k)+d_q(i, k)))
460			zmair = (paprs(i, k)-paprs(i, k+1))/rg
461			zm_solid = (prfl(i, k)-prfl(i, k+1)+psfl(i, k)-psfl(i, k+1))*dtime
462			d_t(i, k) = d_t(i, k) + zlh_solid(i)zm_solid/(zcpairzmair)
463			END DO
464			END DO
465
466	guez	108	contains
467
468	guez	207	! vitesse de chute pour cristaux de glace
469	guez	108
470			REAL function fallvs(zzz)
471	guez	207	REAL, intent(in):: zzz
472	guez	108	fallvs = 3.29/2.0((zzz)0.16)ffallv_lsc
473			end function fallvs
474
475	guez	207	!********************************************************
476
477	guez	108	real function fallvc(zzz)
478	guez	207	REAL, intent(in):: zzz
479	guez	108	fallvc = 3.29/2.0((zzz)0.16)ffallv_con
480			end function fallvc
481
482	guez	68	END SUBROUTINE fisrtilp
483
484			end module fisrtilp_m