trunk/phylmd/fisrtilp.f

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,v 1.2 2004/11/09 16:55:40
  ! Auteur(s): Z.X. Li (LMD/CNRS)
  ! Date: le 20 mars 1995
  ! Objet: condensation et precipitation stratiforme.
  !        schema de nuage

  USE dimens_m
  USE dimphy
  USE tracstoke
  USE suphec_m
  USE yoethf_m
  USE fcttre
  USE comfisrtilp
  use numer_rec_95, only: nr_erf

  IMPLICIT NONE

  ! Arguments:

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