trunk/phylmd/fisrtilp.f90

module fisrtilp_m

  IMPLICIT NONE

contains

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

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

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

    USE numer_rec_95, ONLY: nr_erf

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

    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(dtphys / real(ninter) - 360.) > 0.001) THEN
          PRINT *, "fisrtilp : ce n'est pas pr\'evu, voir Z. X. Li", dtphys
          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) * dtphys &
                  * zcpeau + zmair * zcpair* zt(i)) &
                  / (zmair * zcpair + zrfl(i) * dtphys * 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*dtphys/ &
                     (paprs(i, k)-paprs(i, k+1))
                zqev = min(zqev, zqevt)
                zrfln(i) = zrfl(i) - zqev*(paprs(i, k)-paprs(i, k+1))/rg/dtphys

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