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
    ! 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, INTENT (IN):: q(klon, klev) ! humidite specifique (kg/kg)
    LOGICAL ptconv(klon, klev) ! determine la largeur de distribution de vapeur
    REAL ratqs(klon, klev) ! determine la largeur de distribution de vapeur
    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 rain(klon) ! pluies (mm/s)                                   
    REAL snow(klon) ! neige (mm/s)                                    

    ! Coeffients de fraction lessivee : pour OFF-LINE
    REAL pfrac_impa(klon, klev)
    REAL pfrac_nucl(klon, klev)
    REAL pfrac_1nucl(klon, klev)

    ! Fraction d'aerosols lessivee par impaction et par nucleation
    ! POur ON-LINE
    REAL frac_nucl(klon, klev)

    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)
    REAL rhcl(klon, klev) ! humidite relative en ciel clair            

    ! Local:

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