trunk/phylmd/fisrtilp.f

module fisrtilp_m

  IMPLICIT NONE

contains

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

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

    ! Objet : condensation et précipitation stratiforme, schéma de
    ! nuage, schéma de condensation à grande échelle (pluie).

    USE comfisrtilp, ONLY: cld_lc_con, cld_lc_lsc, cld_tau_con, &
         cld_tau_lsc, coef_eva, ffallv_con, ffallv_lsc, iflag_pdf, reevap_ice
    USE dimphy, ONLY: klev, klon
    USE fcttre, ONLY: dqsatl, dqsats, foede, foeew, qsatl, qsats, thermcep
    USE numer_rec_95, ONLY: nr_erf
    USE suphec_m, ONLY: rcpd, rd, retv, rg, rlstt, rlvtt, rtt
    USE yoethf_m, ONLY: r2es, r5ies, r5les, rvtmp2

    ! Arguments:

    REAL, INTENT (IN):: dtime ! intervalle du temps (s) 
    REAL, INTENT (IN):: paprs(klon, klev+1) ! pression a inter-couche 
    REAL, INTENT (IN):: pplay(klon, klev) ! pression au milieu de couche
    REAL, INTENT (IN):: t(klon, klev) ! temperature (K)
    REAL, INTENT (IN):: q(klon, klev) ! humidite specifique (kg/kg)
    LOGICAL, INTENT (IN):: ptconv(klon, klev)

    REAL, INTENT (IN):: ratqs(klon, klev)
    ! determine la largeur de distribution de vapeur

    REAL, INTENT (out):: d_t(klon, klev) ! incrementation de la temperature (K)
    REAL, INTENT (out):: d_q(klon, klev) ! incrementation de la vapeur d'eau
    REAL, INTENT (out):: d_ql(klon, klev) ! incrementation de l'eau liquide
    REAL, INTENT (out):: rneb(klon, klev) ! fraction nuageuse

    REAL, INTENT (out):: radliq(klon, klev)
    ! eau liquide utilisee dans rayonnement

    REAL, INTENT (out):: rain(klon) ! pluies (mm/s)
    REAL, INTENT (out):: snow(klon) ! neige (mm/s)

    ! Coeffients de fraction lessivee :
    REAL, INTENT (inout):: pfrac_impa(klon, klev)
    REAL, INTENT (inout):: pfrac_nucl(klon, klev)
    REAL, INTENT (inout):: pfrac_1nucl(klon, klev)

    ! Fraction d'aerosols lessivee par impaction
    REAL, INTENT (out):: frac_impa(klon, klev)

    ! Fraction d'aerosols lessivee par nucleation
    REAL, INTENT (out):: frac_nucl(klon, klev)

    REAL, INTENT (out):: prfl(klon, klev+1)
    ! flux d'eau precipitante aux interfaces (kg/m2/s)

    REAL, INTENT (out):: psfl(klon, klev+1)
    ! flux d'eau precipitante aux interfaces (kg/m2/s)

    REAL, INTENT (out):: rhcl(klon, klev) ! humidite relative en ciel clair 

    ! Local:

    REAL zct(klon), zcl(klon)

    ! Options du programme:

    REAL seuil_neb ! un nuage existe vraiment au-dela
    PARAMETER (seuil_neb=0.001)

    INTEGER ninter ! sous-intervals pour la precipitation
    PARAMETER (ninter=5)
    LOGICAL evap_prec ! evaporation de la pluie 
    PARAMETER (evap_prec=.TRUE.)
    REAL zpdf_sig(klon), zpdf_k(klon), zpdf_delta(klon)
    REAL zpdf_a(klon), zpdf_b(klon), zpdf_e1(klon), zpdf_e2(klon)

    LOGICAL cpartiel ! condensation partielle 
    PARAMETER (cpartiel=.TRUE.)
    REAL t_coup
    PARAMETER (t_coup=234.0)

    INTEGER i, k, n, kk
    REAL zqs(klon), zdqs(klon), zcor, zcvm5
    logical zdelta
    REAL zrfl(klon), zrfln(klon), zqev, zqevt
    REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
    REAL ztglace, zt(klon)
    INTEGER nexpo ! exponentiel pour glace/eau 
    REAL zdz(klon), zrho(klon), ztot(klon), zrhol(klon)
    REAL zchau(klon), zfroi(klon), zfice(klon), zneb(klon)

    LOGICAL:: appel1er = .TRUE.

    ! Variables traceurs:
    ! Provisoire !!! Parametres alpha du lessivage
    ! A priori on a 4 scavenging numbers possibles

    REAL, save:: a_tr_sca(4)

    ! Variables intermediaires

    REAL zalpha_tr
    REAL zfrac_lessi
    REAL zprec_cond(klon)
    REAL zmair, zcpair, zcpeau
    ! Pour la conversion eau-neige
    REAL zlh_solid(klon), zm_solid

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

    zdelq = 0.0

    IF (appel1er) THEN
       PRINT *, 'fisrtilp, ninter:', ninter
       PRINT *, 'fisrtilp, evap_prec:', evap_prec
       PRINT *, 'fisrtilp, cpartiel:', cpartiel
       IF (abs(dtime / real(ninter) - 360.) > 0.001) THEN
          PRINT *, "fisrtilp : ce n'est pas prévu, voir Z. X. Li", dtime
          PRINT *, 'Je préfère un sous-intervalle de 6 minutes.'
       END IF
       appel1er = .FALSE.

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

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

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

    ztglace = rtt - 15.0
    nexpo = 6

    ! Initialiser les sorties:

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

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

    ! Initialiser le flux de precipitation a zero

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

    ! Pour plus de securite

    zalpha_tr = 0.
    zfrac_lessi = 0.

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

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

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

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

                ! pour la glace, on réévapore toute la précip dans la
                ! couche du dessous la glace venant de la couche du
                ! dessus est simplement dans la couche du dessous.

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

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

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

       IF (thermcep) THEN
          DO i = 1, klon
             zdelta = rtt >= zt(i)
             zcvm5 = merge(r5ies*rlstt, r5les*rlvtt, zdelta)
             zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zq(i))
             zqs(i) = r2es*foeew(zt(i), zdelta)/pplay(i, k)
             zqs(i) = min(0.5, zqs(i))
             zcor = 1./(1.-retv*zqs(i))
             zqs(i) = zqs(i)*zcor
             zdqs(i) = foede(zt(i), zdelta, zcvm5, zqs(i), zcor)
          END DO
       ELSE
          DO i = 1, klon
             IF (zt(i)<t_coup) THEN
                zqs(i) = qsats(zt(i))/pplay(i, k)
                zdqs(i) = dqsats(zt(i), zqs(i))
             ELSE
                zqs(i) = qsatl(zt(i))/pplay(i, k)
                zdqs(i) = dqsatl(zt(i), zqs(i))
             END IF
          END DO
       END IF

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

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

          ! on prend en compte le réchauffement qui diminue
          ! la partie condensée

          ! Version avec les ratqs

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  contains

    ! vitesse de chute pour crystaux de glace

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

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

  END SUBROUTINE fisrtilp

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