trunk/phylmd/fisrtilp.f90

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
    ! Author: Z. X. Li (LMD/CNRS), 20 mars 1995

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

    USE dimphy, ONLY: klev, klon
    USE suphec_m, ONLY: rcpd, rd, retv, rg, rlstt, rlvtt, rtt
    USE yoethf_m, ONLY: r2es, r5ies, r5les, rvtmp2
    USE fcttre, ONLY: dqsatl, dqsats, foede, foeew, qsatl, qsats, thermcep
    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 numer_rec_95, ONLY: nr_erf

    ! 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 q(klon, klev) ! humidite specifique (kg/kg)                   
    REAL d_t(klon, klev) ! incrementation de la temperature (K)        
    REAL d_q(klon, klev) ! incrementation de la vapeur d'eau           
    REAL d_ql(klon, klev) ! incrementation de l'eau liquide            
    REAL rneb(klon, klev) ! fraction nuageuse                          
    REAL radliq(klon, klev) ! eau liquide utilisee dans rayonnements   
    REAL rhcl(klon, klev) ! humidite relative en ciel clair            
    REAL rain(klon) ! pluies (mm/s)                                   
    REAL snow(klon) ! neige (mm/s)                                    
    REAL prfl(klon, klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
    REAL psfl(klon, klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
    ! Coeffients de fraction lessivee : pour OFF-LINE

    REAL pfrac_nucl(klon, klev)
    REAL pfrac_1nucl(klon, klev)
    REAL pfrac_impa(klon, klev)

    ! Fraction d'aerosols lessivee par impaction et par nucleation
    ! POur ON-LINE

    REAL frac_impa(klon, klev)
    REAL frac_nucl(klon, klev)
    REAL zct(klon), zcl(klon)
    !AA

    ! 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 ratqs(klon, klev) ! determine la largeur de distribution de vapeur
    LOGICAL ptconv(klon, klev) ! determine la largeur de distribution de vapeur
    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)

    ! Variables locales:

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

    LOGICAL appel1er
    SAVE appel1er

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

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

    REAL a_tr_sca(4)
    SAVE a_tr_sca

    ! Variables intermediaires

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

    ! Fonctions en ligne:

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

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

    DATA appel1er/ .TRUE./
    !ym
    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.

       !AA 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

       !AA 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 !  test sur appel1er
    !MAf 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
    !cc      nexpo = 1

    ! 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


    !AA Pour plus de securite

    zalpha_tr = 0.
    zfrac_lessi = 0.

    !AA----------------------------------------------------------

    ! Boucle verticale (du haut vers le bas)

    !IM : klevm1
    klevm1 = klev - 1
    DO  k = klev, 1, -1

       !AA----------------------------------------------------------

       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<=klevm1) 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

       ! Calculer l'evaporation de la precipitation


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

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

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

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

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

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

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

       IF (cpartiel) THEN

          !        print*, 'Dans partiel k=', k

          !   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 condensee

          !   Version avec les raqts

          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)<1.E-15) THEN
                   zq(i) = 1.E-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
          ! iflag_pdf                                               
          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.0, min(1.0, rneb(i, k)))
             !  On ne divise pas par 1+zdqs pour forcer a avoir l'eau
             !  predite par la convection.  ATTENTION !!! Il va
             !  falloir verifier tout ca.
             zcond(i) = max(0.0, zqn(i)-zqs(i))*rneb(i, k)
             !           print*, 'ZDQS ', zdqs(i)
             !--Olivier
             rhcl(i, k) = (zqs(i)+zq(i)-zdelq)/2./zqs(i)
             IF (rneb(i, k)<=0.0) rhcl(i, k) = zq(i)/zqs(i)
             IF (rneb(i, k)>=1.0) rhcl(i, k) = 1.0
             !--fin
          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
          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

       !AA--------------- 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
             !AA 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
       !AA 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.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

       !AA----------------------------------------------------------
       !                     FIN DE BOUCLE SUR K
    end DO

    !AA-----------------------------------------------------------

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

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

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

  END SUBROUTINE fisrtilp

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