source: LMDZ5/trunk/libf/phylmd/fisrtilp.F90 @ 4799

Last change on this file since 4799 was 2969, checked in by fhourdin, 7 years ago

Petit bug dans iflag_t_glace>=1 vu uniquement en debug en 1D quand
exposant < 1. C'est bien de l'avoir quand même dans les sources
au cas où.

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 57.2 KB
RevLine 
[524]1!
[1403]2! $Id: fisrtilp.F90 2969 2017-07-25 17:22:52Z nfevrier $
[524]3!
[1472]4!
5SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
[2086]6     d_t, d_q, d_ql, d_qi, rneb, radliq, rain, snow,          &
[1742]7     pfrac_impa, pfrac_nucl, pfrac_1nucl,               &
8     frac_impa, frac_nucl, beta,                        &
9     prfl, psfl, rhcl, zqta, fraca,                     &
[2236]10     ztv, zpspsk, ztla, zthl, iflag_cld_th,             &
[1849]11     iflag_ice_thermo)
[524]12
[1472]13  !
14  USE dimphy
[2109]15  USE icefrac_lsc_mod ! compute ice fraction (JBM 3/14)
[2311]16  USE print_control_mod, ONLY: prt_level, lunout
[2686]17  USE cloudth_mod
[2703]18  USE ioipsl_getin_p_mod, ONLY : getin_p
[2807]19  USE phys_local_var_mod, ONLY: ql_seri,qs_seri
[2945]20  USE phys_local_var_mod, ONLY: rneblsvol
[2807]21  ! flag to include modifications to ensure energy conservation (if flag >0)
22  USE add_phys_tend_mod, only : fl_cor_ebil
[1472]23  IMPLICIT none
24  !======================================================================
25  ! Auteur(s): Z.X. Li (LMD/CNRS)
26  ! Date: le 20 mars 1995
27  ! Objet: condensation et precipitation stratiforme.
28  !        schema de nuage
[2500]29  ! Fusion de fisrt (physique sursaturation, P. LeVan K. Laval)
30  !             et ilp (il pleut, L. Li)
31  ! Principales parties:
[2807]32  ! P0> Thermalisation des precipitations venant de la couche du dessus
[2500]33  ! P1> Evaporation de la precipitation (qui vient du niveau k+1)
34  ! P2> Formation du nuage (en k)
[2807]35  ! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
36  ! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour
37  ! les valeurs de T et Q initiales
38  ! P2.A.2> Prise en compte du couplage entre eau condensee et T.
39  ! P2.A.3> Calcul des valeures finales associees a la formation des nuages
40  ! P2.B> Nuage "tout ou rien"
41  ! P2.C> Prise en compte de la Chaleur latente apres formation nuage
[2500]42  ! P3> Formation de la precipitation (en k)
[1472]43  !======================================================================
[2807]44  ! JLD:
45  ! * Routine probablement fausse (au moins incoherente) si thermcep = .false.
46  ! * fl_cor_ebil doit etre > 0 ;
47  !   fl_cor_ebil= 0 pour reproduire anciens bugs
[1472]48  !======================================================================
49  include "YOMCST.h"
50  include "fisrtilp.h"
[2006]51  include "nuage.h" ! JBM (3/14)
[1506]52
[1472]53  !
[2500]54  ! Principaux inputs:
[1472]55  !
[2814]56  REAL, INTENT(IN)                              :: dtime  ! intervalle du temps (s)
57  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: paprs  ! pression a inter-couche
58  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay  ! pression au milieu de couche
59  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: t      ! temperature (K)
60  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: q      ! humidite specifique (kg/kg)
61  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv ! points ou le schema de conv. prof. est actif
62  INTEGER,                         INTENT(IN)   :: iflag_cld_th
63  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo
[2500]64  !
[2814]65  ! Inputs lies aux thermiques
66  !
67  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: ztv
68  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zqta, fraca
69  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zpspsk, ztla
70  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zthl
71  !
72  !  Input/output
73  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: ratqs  ! determine la largeur de distribution de vapeur
74  !
[2500]75  ! Principaux outputs:
76  !
[2814]77  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_t  ! incrementation de la temperature (K)
78  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_q  ! incrementation de la vapeur d'eau
79  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_ql ! incrementation de l'eau liquide
80  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_qi ! incrementation de l'eau glace
81  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneb ! fraction nuageuse
82  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radliq ! eau liquide utilisee dans rayonnements
83  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rhcl ! humidite relative en ciel clair
84  REAL, DIMENSION(klon),           INTENT(OUT)  :: rain
85  REAL, DIMENSION(klon),           INTENT(OUT)  :: snow
86  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: prfl
87  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: psfl
[1403]88
[1472]89  !AA
90  ! Coeffients de fraction lessivee : pour OFF-LINE
91  !
[2814]92  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfrac_nucl
93  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfrac_1nucl
94  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfrac_impa
[1472]95  !
96  ! Fraction d'aerosols lessivee par impaction et par nucleation
97  ! POur ON-LINE
98  !
[2814]99  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_impa
100  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl
[1472]101  !AA
[2814]102  ! --------------------------------------------------------------------------------
[1472]103  !
104  ! Options du programme:
105  !
[2923]106  REAL, SAVE :: seuil_neb=0.001 ! un nuage existe vraiment au-dela
107  !$OMP THREADPRIVATE(seuil_neb)
[524]108
[2923]109
[1472]110  INTEGER ninter ! sous-intervals pour la precipitation
111  PARAMETER (ninter=5)
[2923]112  INTEGER,SAVE :: iflag_evap_prec=1 ! evaporation de la pluie
113  !$OMP THREADPRIVATE(iflag_evap_prec)
[1472]114  !
115  LOGICAL cpartiel ! condensation partielle
116  PARAMETER (cpartiel=.TRUE.)
117  REAL t_coup
118  PARAMETER (t_coup=234.0)
[2814]119  REAL DDT0
120  PARAMETER (DDT0=.01)
121  REAL ztfondue
122  PARAMETER (ztfondue=278.15)
123  ! --------------------------------------------------------------------------------
[1472]124  !
125  ! Variables locales:
126  !
127  INTEGER i, k, n, kk
[2923]128  INTEGER,save::itap=0
129  !$OMP THREADPRIVATE(itap)
130
[2814]131  REAL qsl, qsi
132  real zct      ,zcl
133  INTEGER ncoreczq 
134  REAL ctot(klon,klev)
[2945]135  REAL ctot_vol(klon,klev)
[1901]136  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 
[2807]137  REAL zdqsdT_raw(klon)
[1901]138  REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
[2814]139
140  logical lognormale(klon)
141  logical ice_thermo
[1901]142  LOGICAL convergence(klon)
[2086]143  INTEGER n_i(klon), iter
144  REAL cste
[2814]145
146  real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
147  real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
148  real erf   
149  REAL qcloud(klon)
[1901]150 
[1849]151  REAL zrfl(klon), zrfln(klon), zqev, zqevt
152  REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
153  REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
154  REAL zoliqp(klon), zoliqi(klon)
[2006]155  REAL zt(klon)
156! JBM (3/14) nexpo is replaced by exposant_glace
157! REAL nexpo ! exponentiel pour glace/eau
158! INTEGER, PARAMETER :: nexpo=6
159  INTEGER exposant_glace_old
160  REAL t_glace_min_old
[1472]161  REAL zdz(klon),zrho(klon),ztot      , zrhol(klon)
[2923]162  REAL zchau      ,zfroi      ,zfice(klon),zneb(klon),znebprecip(klon)
[2814]163  REAL zmelt, zpluie, zice
[2086]164  REAL dzfice(klon)
[2415]165  REAL zsolid
[2466]166!!!!
167!  Variables pour Bergeron
[2807]168  REAL zcp, coef1, DeltaT, Deltaq, Deltaqprecl
[2466]169  REAL zqpreci(klon), zqprecl(klon)
[2807]170! Variable pour conservation enegie des precipitations
171  REAL zmqc(klon)
[1472]172  !
173  LOGICAL appel1er
174  SAVE appel1er
175  !$OMP THREADPRIVATE(appel1er)
176  !
[2703]177! iflag_oldbug_fisrtilp=0 enleve le BUG par JYG : tglace_min -> tglace_max
178! iflag_oldbug_fisrtilp=1 ajoute le BUG
179  INTEGER,SAVE :: iflag_oldbug_fisrtilp=0 !=0 sans bug
180  !$OMP THREADPRIVATE(iflag_oldbug_fisrtilp)
[1472]181  !---------------------------------------------------------------
182  !
183  !AA Variables traceurs:
184  !AA  Provisoire !!! Parametres alpha du lessivage
185  !AA  A priori on a 4 scavenging # possibles
186  !
187  REAL a_tr_sca(4)
188  save a_tr_sca
189  !$OMP THREADPRIVATE(a_tr_sca)
190  !
191  ! Variables intermediaires
192  !
193  REAL zalpha_tr
194  REAL zfrac_lessi
195  REAL zprec_cond(klon)
196  !AA
[1742]197! RomP >>> 15 nov 2012
198  REAL   beta(klon,klev) ! taux de conversion de l'eau cond
199! RomP <<<
[2807]200  REAL zmair(klon), zcpair, zcpeau
[1472]201  !     Pour la conversion eau-neige
202  REAL zlh_solid(klon), zm_solid
203  !---------------------------------------------------------------
204  !
205  ! Fonctions en ligne:
206  !
[2500]207  REAL fallvs,fallvc ! Vitesse de chute pour cristaux de glace
208                     ! (Heymsfield & Donner, 1990)
[1472]209  REAL zzz
[2807]210
[1472]211  include "YOETHF.h"
212  include "FCTTRE.h"
213  fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
214  fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc
215  !
216  DATA appel1er /.TRUE./
217  !ym
[2086]218!CR: pour iflag_ice_thermo=2, on active que la convection
219!  ice_thermo = iflag_ice_thermo .GE. 1
[2923]220
[2945]221  itap=itap+1
222  znebprecip(:)=0.
[2923]223
[2945]224  ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
[1472]225  zdelq=0.0
[2945]226  ctot_vol(1:klon,1:klev)=0.0
227  rneblsvol(1:klon,1:klev)=0.0
[524]228
[1506]229  if (prt_level>9)write(lunout,*)'NUAGES4 A. JAM'
[1472]230  IF (appel1er) THEN
[2703]231     CALL getin_p('iflag_oldbug_fisrtilp',iflag_oldbug_fisrtilp)
[2923]232     CALL getin_p('iflag_evap_prec',iflag_evap_prec)
233     CALL getin_p('seuil_neb',seuil_neb)
[2703]234     write(lunout,*)' iflag_oldbug_fisrtilp =',iflag_oldbug_fisrtilp
[1472]235     !
[1575]236     WRITE(lunout,*) 'fisrtilp, ninter:', ninter
[2923]237     WRITE(lunout,*) 'fisrtilp, iflag_evap_prec:', iflag_evap_prec
[1575]238     WRITE(lunout,*) 'fisrtilp, cpartiel:', cpartiel
[2923]239     
[1472]240     IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
[1575]241        WRITE(lunout,*) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
242        WRITE(lunout,*) 'Je prefere un sous-intervalle de 6 minutes'
[1472]243        !         CALL abort
244     ENDIF
245     appel1er = .FALSE.
246     !
247     !AA initialiation provisoire
248     a_tr_sca(1) = -0.5
249     a_tr_sca(2) = -0.5
250     a_tr_sca(3) = -0.5
251     a_tr_sca(4) = -0.5
252     !
253     !AA Initialisation a 1 des coefs des fractions lessivees
254     !
255     !cdir collapse
256     DO k = 1, klev
257        DO i = 1, klon
258           pfrac_nucl(i,k)=1.
259           pfrac_1nucl(i,k)=1.
260           pfrac_impa(i,k)=1.
[1742]261           beta(i,k)=0.  !RomP initialisation
[1472]262        ENDDO
263     ENDDO
[524]264
[1472]265  ENDIF          !  test sur appel1er
266  !
267  !MAf Initialisation a 0 de zoliq
268  !      DO i = 1, klon
269  !         zoliq(i)=0.
270  !      ENDDO
271  ! Determiner les nuages froids par leur temperature
272  !  nexpo regle la raideur de la transition eau liquide / eau glace.
273  !
[2086]274!CR: on est oblige de definir des valeurs fisrt car les valeurs de newmicro ne sont pas les memes par defaut
[2006]275  IF (iflag_t_glace.EQ.0) THEN
276!   ztglace = RTT - 15.0
277    t_glace_min_old = RTT - 15.0
278    !AJ<
279    IF (ice_thermo) THEN
280!     nexpo = 2
281      exposant_glace_old = 2
282    ELSE
283!     nexpo = 6
284      exposant_glace_old = 6
285    ENDIF
[2086]286   
[1849]287  ENDIF
[2006]288 
[1849]289!!  RLVTT = 2.501e6 ! pas de redefinition des constantes physiques (jyg)
290!!  RLSTT = 2.834e6 ! pas de redefinition des constantes physiques (jyg)
291!>AJ
[1472]292  !cc      nexpo = 1
293  !
294  ! Initialiser les sorties:
295  !
296  !cdir collapse
297  DO k = 1, klev+1
298     DO i = 1, klon
299        prfl(i,k) = 0.0
300        psfl(i,k) = 0.0
301     ENDDO
302  ENDDO
[524]303
[1472]304  !cdir collapse
305  DO k = 1, klev
306     DO i = 1, klon
307        d_t(i,k) = 0.0
308        d_q(i,k) = 0.0
309        d_ql(i,k) = 0.0
[2086]310        d_qi(i,k) = 0.0
[1472]311        rneb(i,k) = 0.0
312        radliq(i,k) = 0.0
313        frac_nucl(i,k) = 1.
314        frac_impa(i,k) = 1.
315     ENDDO
316  ENDDO
317  DO i = 1, klon
318     rain(i) = 0.0
319     snow(i) = 0.0
320     zoliq(i)=0.
321     !     ENDDO
322     !
323     ! Initialiser le flux de precipitation a zero
324     !
325     !     DO i = 1, klon
326     zrfl(i) = 0.0
[1849]327     zifl(i) = 0.0
[1472]328     zneb(i) = seuil_neb
329  ENDDO
330  !
331  !
332  !AA Pour plus de securite
[524]333
[1472]334  zalpha_tr   = 0.
335  zfrac_lessi = 0.
[524]336
[2500]337  !AA==================================================================
[1472]338  !
339  ncoreczq=0
[2500]340  ! BOUCLE VERTICALE (DU HAUT VERS LE BAS)
[1472]341  !
342  DO k = klev, 1, -1
343     !
[2500]344     !AA===============================================================
[1472]345     !
[2500]346     ! Initialisation temperature et vapeur
[1472]347     DO i = 1, klon
348        zt(i)=t(i,k)
349        zq(i)=q(i,k)
350     ENDDO
351     !
[2807]352     ! ----------------------------------------------------------------
353     ! P0> Thermalisation des precipitations venant de la couche du dessus
354     ! ----------------------------------------------------------------
[1472]355     ! Calculer la varition de temp. de l'air du a la chaleur sensible
[2807]356     ! transporter par la pluie. On thermalise la pluie avec l'air de la couche.
357     ! Cette quantite de pluie qui est thermalisee, et devra continue a l'etre lors
358     ! des differentes transformations thermodynamiques. Cette masse d'eau doit
359     ! donc etre ajoute a l'humidite de la couche lorsque l'on calcule la variation
360     ! de l'enthalpie  de la couche avec la temperature
361     ! Variables calculees ou modifiees:
362     !   -  zt: temperature de la cocuhe
363     !   - zmqc: masse de precip qui doit etre thermalisee
[1472]364     !
365     IF(k.LE.klevm1) THEN         
366        DO i = 1, klon
367           !IM
[2807]368           zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
369           ! il n'y a pas encore d'eau liquide ni glace dans la maiille, donc zq suffit
[1472]370           zcpair=RCPD*(1.0+RVTMP2*zq(i))
371           zcpeau=RCPD*RVTMP2
[2807]372         if (fl_cor_ebil .GT. 0) then
373           ! zmqc: masse de precip qui doit etre thermalisee avec l'air de la couche atm
374           ! pour s'assurer que la precip arrivant au sol aura bien la temperature de la
375           ! derniere couche
376           zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair(i)
377           ! t(i,k+1)+d_t(i,k+1): nouvelle temp de la couche au dessus
378           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) &
379                 / (zcpair + zmqc(i)*zcpeau)
380         else ! si on maintient les anciennes erreurs
[1472]381           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau &
[2807]382                + zmair(i)*zcpair*zt(i) ) &
383                / (zmair(i)*zcpair + zrfl(i)*dtime*zcpeau)
384         end if
[1472]385        ENDDO
[2885]386     ELSE  ! IF(k.LE.klevm1)
387        DO i = 1, klon
388           zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
389           zmqc(i) = 0.
390        ENDDO
[2807]391     ENDIF ! end IF(k.LE.klevm1)
392     !
[2500]393     ! ----------------------------------------------------------------
[2807]394     ! P1> Calcul de l'evaporation de la precipitation
[2500]395     ! ----------------------------------------------------------------
[2807]396     ! On evapore une partie des precipitations venant de la maille du dessus.
397     ! On calcule l'evaporation et la sublimation des precipitations, jusqu'a
398     ! ce que la fraction de cette couche qui est sous le nuage soit saturee.
399     ! Variables calculees ou modifiees:
400     !   - zrfl et zifl: flux de precip liquide et glace
401     !   - zq, zt: humidite et temperature de la cocuhe
402     !   - zmqc: masse de precip qui doit etre thermalisee
403     !
[2923]404     IF (iflag_evap_prec>=1) THEN
[1472]405        DO i = 1, klon
[2807]406!          S'il y a des precipitations
[1849]407           IF (zrfl(i)+zifl(i).GT.0.) THEN
[2500]408              ! Calcul du qsat
[1472]409              IF (thermcep) THEN
410                 zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
411                 zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
412                 zqs(i)=MIN(0.5,zqs(i))
413                 zcor=1./(1.-RETV*zqs(i))
414                 zqs(i)=zqs(i)*zcor
415              ELSE
416                 IF (zt(i) .LT. t_coup) THEN
417                    zqs(i) = qsats(zt(i)) / pplay(i,k)
418                 ELSE
419                    zqs(i) = qsatl(zt(i)) / pplay(i,k)
420                 ENDIF
421              ENDIF
[1849]422           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
423        ENDDO
424!AJ<
[2807]425
[1849]426       IF (.NOT. ice_thermo) THEN
427        DO i = 1, klon
[2807]428!          S'il y a des precipitations
[1849]429           IF (zrfl(i)+zifl(i).GT.0.) THEN
[2500]430                ! Evap max pour ne pas saturer la fraction sous le nuage
[2807]431                ! Evap max jusqu'à atteindre la saturation dans la partie
432                ! de la maille qui est sous le nuage de la couche du dessus
433                !!! On ne tient compte de cette fraction que sous une seule
434                !!! couche sous le nuage
[1849]435                zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
[2807]436             ! Ajout de la prise en compte des precip a thermiser
437             ! avec petite reecriture
438             if  (fl_cor_ebil .GT. 0) then ! nouveau
[2500]439                ! Calcul de l'evaporation du flux de precip herite
440                !   d'au-dessus
[1849]441                zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
[2807]442                     * zmair(i)/pplay(i,k)*zt(i)*RD
443                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) * dtime/zmair(i)
444
445                ! Seuil pour ne pas saturer la fraction sous le nuage
446                zqev = MIN (zqev, zqevt)
447                ! Nouveau flux de precip
448                zrfln(i) = zrfl(i) - zqev*zmair(i)/dtime
449                ! Aucun flux liquide pour T < t_coup, on reevapore tout.
450                IF (zt(i) .LT. t_coup.and.reevap_ice) THEN
451                  zrfln(i)=0.
452                  zqev = (zrfl(i)-zrfln(i))/zmair(i)*dtime
453                END IF
454                ! Nouvelle vapeur
455                zq(i) = zq(i) + zqev
456                zmqc(i) = zmqc(i)-zqev
457                ! Nouvelle temperature (chaleur latente)
458                zt(i) = zt(i) - zqev &
459                     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
[2814]460!!JLD debut de partie a supprimer a terme
[2807]461            else ! if  (fl_cor_ebil .GT. 0)
462                ! Calcul de l'evaporation du flux de precip herite
463                !   d'au-dessus
464                zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
[1849]465                     * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
466                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
467                     * RG*dtime/(paprs(i,k)-paprs(i,k+1))
[2500]468                ! Seuil pour ne pas saturer la fraction sous le nuage
[1849]469                zqev = MIN (zqev, zqevt)
[2500]470                ! Nouveau flux de precip
[1849]471                zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
472                     /RG/dtime
[2500]473                ! Aucun flux liquide pour T < t_coup
[1849]474                IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
[2500]475                ! Nouvelle vapeur
[1849]476                zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
477                     * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
[2500]478                ! Nouvelle temperature (chaleur latente)
[1849]479                zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
480                     * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
481                     * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
[2807]482              end if ! if  (fl_cor_ebil .GT. 0)
[2814]483!!JLD fin de partie a supprimer a terme
[1849]484                zrfl(i) = zrfln(i)
485                zifl(i) = 0.
486           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
487        ENDDO
488!
489       ELSE ! (.NOT. ice_thermo)
[2807]490!      ================================
491!      Avec thermodynamique de la glace
492!      ================================
[1849]493        DO i = 1, klon
494!AJ<
[2807]495!        S'il y a des precipitations
496         IF (zrfl(i)+zifl(i).GT.0.) THEN
[2923]497
498         IF (iflag_evap_prec==1) THEN
499            znebprecip(i)=zneb(i)
500         ELSE
501            znebprecip(i)=MAX(zneb(i),znebprecip(i))
502         ENDIF
[1849]503     
[2807]504        ! Evap max pour ne pas saturer la fraction sous le nuage
[2923]505         zqev0 = MAX (0.0, (zqs(i)-zq(i))*znebprecip(i) )
[524]506
[2807]507         !JAM
508         ! On differencie qsat pour l'eau et la glace
509         ! Si zdelta=1. --> glace
510         ! Si zdelta=0. --> eau liquide
[2500]511       
512         ! Calcul du qsat par rapport a l'eau liquide
[1849]513         qsl= R2ES*FOEEW(zt(i),0.)/pplay(i,k)
514         qsl= MIN(0.5,qsl)
515         zcor= 1./(1.-RETV*qsl)
516         qsl= qsl*zcor
517         
[2807]518         ! Calcul de l'evaporation du flux de precip venant du dessus
[2500]519         ! Formulation en racine du flux de precip
520         ! (Klemp & Wilhelmson, 1978; Sundqvist, 1988)
[2962]521         IF (iflag_evap_prec==3) THEN
522         zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl) &
523              *SQRT(zrfl(i)/max(1.e-4,znebprecip(i))) &
524              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
525         ELSE
[1849]526         zqevt = 1.*coef_eva*(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
527              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
[2962]528         ENDIF
529
530
[1849]531         zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
532              *RG*dtime/(paprs(i,k)-paprs(i,k+1))
[2500]533         
534         ! Calcul du qsat par rapport a la glace
[1849]535         qsi= R2ES*FOEEW(zt(i),1.)/pplay(i,k)
536         qsi= MIN(0.5,qsi)
537         zcor= 1./(1.-RETV*qsi)
538         qsi= qsi*zcor
[1472]539
[2500]540         ! Calcul de la sublimation du flux de precip solide herite
541         !   d'au-dessus
[2962]542         IF (iflag_evap_prec==3) THEN
543         zqevti = znebprecip(i)*coef_eva*(1.0-zq(i)/qsi) &
544              *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
545              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
546         ELSE
[1849]547         zqevti = 1.*coef_eva*(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
548              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
[2962]549         ENDIF
[1849]550         zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
551              *RG*dtime/(paprs(i,k)-paprs(i,k+1))   
552
[2807]553        !JAM
554        ! Limitation de l'evaporation. On s'assure qu'on ne sature pas
555        ! la fraction de la couche sous le nuage sinon on repartit zqev0
556        ! en conservant la proportion liquide / glace
[1849]557     
558         IF (zqevt+zqevti.GT.zqev0) THEN
[2807]559            zqev=zqev0*zqevt/(zqevt+zqevti)
560            zqevi=zqev0*zqevti/(zqevt+zqevti)
[1849]561         ELSE
[2807]562!JLD je ne comprends pas les lignes ci-dessous. On repartit les precips
563!       liquides et solides meme si on ne sature pas la couche.
564!       A mon avis, le test est inutile, et il faudrait tout remplacer par:
565!            zqev=zqevt
566!            zqevi=zqevti
[1849]567             IF (zqevt+zqevti.GT.0.) THEN
[2807]568                zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
569                zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
[1849]570             ELSE
571             zqev=0.
572             zqevi=0.
573             ENDIF
574         ENDIF
[2500]575         ! Nouveaux flux de precip liquide et solide
[1849]576         zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
577                                 /RG/dtime)
578         zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
579                                 /RG/dtime)
[2500]580
581         ! Mise a jour de la vapeur, temperature et flux de precip
[1849]582         zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
583                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
[2807]584       if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips
585         zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
586                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
[1849]587         zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
588                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
[2807]589                  * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
590                  + (zifln(i)-zifl(i)) &
591                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
592                  * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
593       else ! sans correction thermalisation des precips
594         zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
595                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
[1849]596                  * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
597                  + (zifln(i)-zifl(i)) &
598                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
599                  * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
[2807]600       end if
601        ! Nouvelles vaeleurs des precips liquides et solides
[1849]602         zrfl(i) = zrfln(i)
603         zifl(i) = zifln(i)
[2923]604!        print*,'REEVAP ',itap,k,znebprecip(1),zqev0,zqev,zqevi,zrfl(1)
[1849]605
[2086]606!CR ATTENTION: deplacement de la fonte de la glace
[2466]607!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
608!!!        zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2  !!!!!!!!! jyg
609!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
610           zmelt = ((zt(i)-273.15)/(ztfondue-273.15))             ! jyg
[2807]611           ! fraction de la precip solide qui est fondue
[2086]612           zmelt = MIN(MAX(zmelt,0.),1.)
[2500]613           ! Fusion de la glace
[2086]614           zrfl(i)=zrfl(i)+zmelt*zifl(i)
[2807]615           if (fl_cor_ebil .LE. 0) then
616             ! the following line should not be here. Indeed, if zifl is modified
617             ! now, zifl(i)*zmelt is no more the amount of ice that has melt
618             ! and therefore the change in temperature computed below is wrong
619             zifl(i)=zifl(i)*(1.-zmelt)
620           end if
[2500]621           ! Chaleur latente de fusion
[2807]622        if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips
[2086]623           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
[2807]624                      *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
625        else ! sans correction thermalisation des precips
626           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
[2086]627                      *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
[2807]628        end if
629           if (fl_cor_ebil .GT. 0) then ! correction bug, deplacement ligne precedente
630             zifl(i)=zifl(i)*(1.-zmelt)
631           end if
[2086]632
[2923]633           ELSE
634              ! Si on n'a plus de pluies, on reinitialise a 0 la farcion
635              ! sous nuageuse utilisee pour la pluie.
636              znebprecip(i)=0.
[1849]637           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
[1472]638        ENDDO
[1849]639
640      ENDIF ! (.NOT. ice_thermo)
641     
[2500]642     ! ----------------------------------------------------------------
643     ! Fin evaporation de la precipitation
644     ! ----------------------------------------------------------------
[2923]645     ENDIF ! (iflag_evap_prec>=1)
[1472]646     !
647     ! Calculer Qs et L/Cp*dQs/dT:
648     !
649     IF (thermcep) THEN
650        DO i = 1, klon
[524]651           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
652           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
[2807]653       if  (fl_cor_ebil .GT. 0) then ! nouveau
654           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
655       else   
[524]656           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
[2807]657       endif
[524]658           zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
659           zqs(i) = MIN(0.5,zqs(i))
660           zcor = 1./(1.-RETV*zqs(i))
661           zqs(i) = zqs(i)*zcor
662           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
[2807]663           zdqsdT_raw(i) = zdqs(i)*  &
664         &         RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
[1472]665        ENDDO
666     ELSE
667        DO i = 1, klon
668           IF (zt(i).LT.t_coup) THEN
669              zqs(i) = qsats(zt(i))/pplay(i,k)
670              zdqs(i) = dqsats(zt(i),zqs(i))
671           ELSE
672              zqs(i) = qsatl(zt(i))/pplay(i,k)
673              zdqs(i) = dqsatl(zt(i),zqs(i))
674           ENDIF
675        ENDDO
676     ENDIF
677     !
678     ! Determiner la condensation partielle et calculer la quantite
679     ! de l'eau condensee:
680     !
[1901]681!verification de la valeur de iflag_fisrtilp_qsat pour iflag_ice_thermo=1
[2086]682!       if ((iflag_ice_thermo.eq.1).and.(iflag_fisrtilp_qsat.ne.0)) then
683!         write(*,*) " iflag_ice_thermo==1 requires iflag_fisrtilp_qsat==0", &
684!        " but iflag_fisrtilp_qsat=",iflag_fisrtilp_qsat, ". Might as well stop here."
685!         stop
686!       endif
[1403]687
[2500]688     ! ----------------------------------------------------------------
689     ! P2> Formation du nuage
690     ! ----------------------------------------------------------------
[2807]691     ! Variables calculees:
692     !   rneb  : fraction nuageuse
693     !   zcond : eau condensee moyenne dans la maille.
694     !   rhcl: humidite relative ciel-clair
695     !   zt : temperature de la maille
696     ! ----------------------------------------------------------------
697     !
[1472]698     IF (cpartiel) THEN
[2807]699        ! -------------------------
700        ! P2.A> Nuage fractionnaire
701        ! -------------------------
[1472]702        !
703        !   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
704        !   nuageuse a partir des PDF de Sandrine Bony.
705        !   rneb  : fraction nuageuse
706        !   zqn   : eau totale dans le nuage
707        !   zcond : eau condensee moyenne dans la maille.
708        !  on prend en compte le réchauffement qui diminue la partie
709        ! condensee
710        !
711        !   Version avec les raqts
[524]712
[2807]713        ! ----------------------------------------------------------------
714        ! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
715        ! ----------------------------------------------------------------
[1472]716        if (iflag_pdf.eq.0) then
[524]717
[2500]718           ! version creneau de (Li, 1998)
[524]719           do i=1,klon
[1472]720              zdelq = min(ratqs(i,k),0.99) * zq(i)
721              rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
722              zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
[524]723           enddo
724
[2807]725        else !  if (iflag_pdf.eq.0)
726           ! ----------------------------------------------------------------
727           ! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour
728           ! les valeurs de T et Q initiales
729           ! ----------------------------------------------------------------
[524]730           do i=1,klon
731              if(zq(i).lt.1.e-15) then
[1472]732                 ncoreczq=ncoreczq+1
733                 zq(i)=1.e-15
[524]734              endif
[1472]735           enddo
[1403]736
[2236]737           if (iflag_cld_th>=5) then
[1403]738
[2696]739              if (iflag_cloudth_vert<=2) then
[2686]740               call cloudth(klon,klev,k,ztv, &
[1472]741                   zq,zqta,fraca, &
742                   qcloud,ctot,zpspsk,paprs,ztla,zthl, &
743                   ratqs,zqs,t)
[2911]744              elseif (iflag_cloudth_vert>=3) then
[2686]745               call cloudth_v3(klon,klev,k,ztv, &
746                   zq,zqta,fraca, &
[2945]747                   qcloud,ctot,ctot_vol,zpspsk,paprs,ztla,zthl, &
[2686]748                   ratqs,zqs,t)
749              endif
[1472]750              do i=1,klon
[1403]751                 rneb(i,k)=ctot(i,k)
[2945]752                 rneblsvol(i,k)=ctot_vol(i,k)
[1403]753                 zqn(i)=qcloud(i)
[1472]754              enddo
[1403]755
[1472]756           endif
757
[2236]758           if (iflag_cld_th <= 4) then
[1472]759              lognormale = .true.
[2236]760           elseif (iflag_cld_th >= 6) then
[1472]761              ! lognormale en l'absence des thermiques
762              lognormale = fraca(:,k) < 1e-10
763           else
[2236]764              ! Dans le cas iflag_cld_th=5, on prend systématiquement la
[1472]765              ! bi-gaussienne
766              lognormale = .false.
767           end if
768
[2500]769!CR: variation de qsat avec T en presence de glace ou non
[2086]770!initialisations
[1472]771           do i=1,klon
[2086]772              DT(i) = 0.
773              n_i(i)=0
[1901]774              Tbef(i)=zt(i)
[2086]775              qlbef(i)=0.
776           enddo
777
[2807]778        ! ----------------------------------------------------------------
779        ! P2.A.2> Prise en compte du couplage entre eau condensee et T.
780        ! Calcul des grandeurs nuageuses en tenant compte de l'effet de
781        ! la condensation sur T, et donc sur qsat et sur les grandeurs nuageuses
782        ! qui en dependent. Ce changement de temperature est provisoire, et
783        ! la valeur definitive sera calcule plus tard.
784        ! Variables calculees:
785        !   rneb : nebulosite
786        !   zcond: eau condensee en moyenne dans la maille
787        ! note JLD: si on n'a pas de pdf lognormale, ce qui se passe ne me semble
788        ! pas clair, il n'y a probablement pas de prise en compte de l'effet de
789        ! T sur qsat
790        ! ----------------------------------------------------------------
[2086]791
792!Boucle iterative: ATTENTION, l'option -1 n'est plus activable ici
793           if (iflag_fisrtilp_qsat.ge.0) then
[2500]794             ! Iteration pour condensation avec variation de qsat(T)
795             ! -----------------------------------------------------
[2086]796             do iter=1,iflag_fisrtilp_qsat+1
797               
798               do i=1,klon
799!                 do while ((abs(DT(i)).gt.DDT0).or.(n_i(i).eq.0))
[2807]800                 ! !! convergence = .true. tant que l'on n'a pas converge !!
801                 ! ------------------------------
[2086]802                 convergence(i)=abs(DT(i)).gt.DDT0
803                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
[2807]804                 ! si on n'a pas converge
805                 !
806                 ! P2.A.2.1> Calcul de la fraction nuageuse et de la quantite d'eau condensee
807                 ! ---------------------------------------------------------------
808                 ! Variables calculees:
809                 ! rneb : nebulosite
810                 ! zqn : eau condensee, dans le nuage (in cloud water content)
811                 ! zcond: eau condensee en moyenne dans la maille
812                 ! rhcl: humidite relative ciel-clair
813                 !
814                 Tbef(i)=Tbef(i)+DT(i) ! nouvelle temperature
[2086]815                 if (.not.ice_thermo) then   
816                    zdelta = MAX(0.,SIGN(1.,RTT-Tbef(i)))
817                 else
818                    if (iflag_t_glace.eq.0) then
819                    zdelta = MAX(0.,SIGN(1.,t_glace_min_old-Tbef(i)))
[2507]820                    else if (iflag_t_glace.ge.1) then
[2703]821                       if (iflag_oldbug_fisrtilp.EQ.0) then
822                          zdelta = MAX(0.,SIGN(1.,t_glace_max-Tbef(i)))
823                       else
[2807]824                          !avec bug : zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
[2703]825                          zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
826                       endif
[2086]827                    endif
828                 endif
[2807]829                 ! Calcul de rneb, qzn et zcond pour les PDF lognormales
[2086]830                 zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
[2807]831               if (fl_cor_ebil .GT. 0) then
832                 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
833               else
[2086]834                 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
[2807]835               end if
[2086]836                 zqs(i) = R2ES*FOEEW(Tbef(i),zdelta)/pplay(i,k)
837                 zqs(i) = MIN(0.5,zqs(i))
838                 zcor = 1./(1.-RETV*zqs(i))
839                 zqs(i) = zqs(i)*zcor
840                 zdqs(i) = FOEDE(Tbef(i),zdelta,zcvm5,zqs(i),zcor)
[1472]841                 zpdf_sig(i)=ratqs(i,k)*zq(i)
842                 zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
843                 zpdf_delta(i)=log(zq(i)/zqs(i))
844                 zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
845                 zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
846                 zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
847                 zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
848                 zpdf_e1(i)=1.-erf(zpdf_e1(i))
849                 zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
850                 zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
851                 zpdf_e2(i)=1.-erf(zpdf_e2(i))
[1901]852             
853                 if (zpdf_e1(i).lt.1.e-10) then
854                    rneb(i,k)=0.
855                    zqn(i)=zqs(i)
856                 else
857                    rneb(i,k)=0.5*zpdf_e1(i)
858                    zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
859                 endif
860
[2956]861                 ! If vertical heterogeneity, change fraction by volume as well
862                 if (iflag_cloudth_vert>=3) then
863                   ctot_vol(i,k)=rneb(i,k)
864                   rneblsvol(i,k)=ctot_vol(i,k)
865                 endif
866
[2086]867                 endif !convergence
[2956]868
[2086]869               enddo ! boucle en i
870
[2807]871                 ! P2.A.2.2> Calcul APPROCHE de la variation de temperature DT
872                 !         due a la condensation.
873                 ! ---------------------------------------------------------------
874                 ! Variables calculees:
875                 ! DT : variation de temperature due a la condensation
876
[2086]877                 if (.not. ice_thermo) then
[2807]878                 ! --------------------------
[2086]879                 do i=1,klon
880                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
881
[1901]882                 qlbef(i)=max(0.,zqn(i)-zqs(i))
[2807]883               if (fl_cor_ebil .GT. 0) then
884                 num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
885               else
[1901]886                 num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
[2807]887               end if
[1901]888                 denom=1.+rneb(i,k)*zdqs(i)
889                 DT(i)=num/denom
[2086]890                 n_i(i)=n_i(i)+1
891                 endif
892                 enddo
[1403]893
[2807]894                 else ! if (.not. ice_thermo)
895                 ! --------------------------
[2507]896                 if (iflag_t_glace.ge.1) then
[2109]897                 CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
[1472]898                 endif
[1411]899
[2086]900                 do i=1,klon
901                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
902                 
903                 if (iflag_t_glace.eq.0) then
904                    zfice(i) = 1.0 - (Tbef(i)-t_glace_min_old) / (RTT-t_glace_min_old)
905                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
906                    zfice(i) = zfice(i)**exposant_glace_old
[2807]907                    dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) &
908          &                     / (t_glace_min_old - RTT)
[1901]909                 endif
[2086]910                 
[2969]911                 if (iflag_t_glace.ge.1.and.zfice(i)>0.) then
[2807]912                 dzfice(i)= exposant_glace * zfice(i)**(exposant_glace-1) &
913          &                    / (t_glace_min - t_glace_max)
[2086]914                 endif
915               
916                 if ((zfice(i).eq.0).or.(zfice(i).eq.1)) then
917                    dzfice(i)=0.
918                 endif
[1411]919
[2086]920                 if (zfice(i).lt.1) then
921                    cste=RLVTT
922                 else
923                    cste=RLSTT
924                 endif
925
926                 qlbef(i)=max(0.,zqn(i)-zqs(i))
[2807]927               if (fl_cor_ebil .GT. 0) then
928                 num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
929           &          +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
930                 denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
931                         -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)    &
932           &               *qlbef(i)*dzfice(i)
933               else
934                 num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
935           &         +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
936                 denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
[2086]937                         -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))*rneb(i,k)*qlbef(i)*dzfice(i)
[2807]938               end if
[2086]939                 DT(i)=num/denom
940                 n_i(i)=n_i(i)+1
941
942                 endif ! fin convergence
943                 enddo ! fin boucle i
944
945                 endif !ice_thermo
946
[2500]947             enddo ! iter=1,iflag_fisrtilp_qsat+1
948             ! Fin d'iteration pour condensation avec variation de qsat(T)
949             ! -----------------------------------------------------------
[2807]950           endif !  if (iflag_fisrtilp_qsat.ge.0)
951     ! ----------------------------------------------------------------
952     ! Fin de P2.A.2> la prise en compte du couplage entre eau condensee et T
953     ! ----------------------------------------------------------------
[524]954
955        endif ! iflag_pdf
956
[2086]957!        if (iflag_fisrtilp_qsat.eq.-1) then
[2500]958!------------------------------------------
959!CR: ATTENTION option fausse mais a existe:
960! pour la re-activer, prendre iflag_fisrtilp_qsat=0 et
961! activer les lignes suivantes:
[2086]962       IF (1.eq.0) THEN
963       DO i=1,klon
[1146]964           IF (rneb(i,k) .LE. 0.0) THEN
965              zqn(i) = 0.0
966              rneb(i,k) = 0.0
967              zcond(i) = 0.0
968              rhcl(i,k)=zq(i)/zqs(i)
969           ELSE IF (rneb(i,k) .GE. 1.0) THEN
970              zqn(i) = zq(i)
[1901]971              rneb(i,k) = 1.0                 
972              zcond(i) = MAX(0.0,zqn(i)-zqs(i))/(1+zdqs(i))
[1146]973              rhcl(i,k)=1.0
974           ELSE
[1901]975              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1+zdqs(i))
[1146]976              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
977           ENDIF
[2500]978       ENDDO
979       ENDIF
980!------------------------------------------
[1901]981
[2086]982!        ELSE
[2807]983        ! ----------------------------------------------------------------
984        ! P2.A.3> Calcul des valeures finales associees a la formation des nuages
985        ! Variables calculees:
986        !   rneb : nebulosite
987        !   zcond: eau condensee en moyenne dans la maille
988        !   zq : eau vapeur dans la maille
989        !   zt : temperature de la maille
990        !   rhcl: humidite relative ciel-clair
991        ! ----------------------------------------------------------------
992        !
993        ! Bornage de l'eau in-cloud (zqn) et de la fraction nuageuse (rneb)
994        ! Calcule de l'eau condensee moyenne dans la maille (zcond),
995        ! et de l'humidite relative ciel-clair (rhcl)
[1901]996        DO i=1,klon
997           IF (rneb(i,k) .LE. 0.0) THEN
998              zqn(i) = 0.0
999              rneb(i,k) = 0.0
1000              zcond(i) = 0.0
1001              rhcl(i,k)=zq(i)/zqs(i)
1002           ELSE IF (rneb(i,k) .GE. 1.0) THEN
1003              zqn(i) = zq(i)
1004              rneb(i,k) = 1.0
1005              zcond(i) = MAX(0.0,zqn(i)-zqs(i))
1006              rhcl(i,k)=1.0
1007           ELSE
1008              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
1009              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
1010           ENDIF
1011        ENDDO
[2956]1012        ! If vertical heterogeneity, change fraction by volume as well
1013        if (iflag_cloudth_vert>=3) then
1014          ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
1015          rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
1016        endif
[1901]1017
[2086]1018!        ENDIF
[1901]1019
[2500]1020     ELSE ! de IF (cpartiel)
[2807]1021        ! -------------------------
1022        ! P2.B> Nuage "tout ou rien"
1023        ! -------------------------
1024        ! note JLD: attention, rhcl non calcule. Ca peut avoir des consequences?
[1472]1025        DO i = 1, klon
1026           IF (zq(i).GT.zqs(i)) THEN
1027              rneb(i,k) = 1.0
1028           ELSE
1029              rneb(i,k) = 0.0
1030           ENDIF
1031           zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
1032        ENDDO
[2807]1033     ENDIF ! de IF (cpartiel)
[1472]1034     !
[2500]1035     ! Mise a jour vapeur d'eau
[2807]1036     ! -------------------------
[1472]1037     DO i = 1, klon
1038        zq(i) = zq(i) - zcond(i)
1039        !         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
1040     ENDDO
[1849]1041!AJ<
[2807]1042     ! ------------------------------------
1043     ! P2.C> Prise en compte de la Chaleur latente apres formation nuage
[2500]1044     ! -------------------------------------
[2807]1045     ! Variable calcule:
1046     !   zt : temperature de la maille
1047     !
[1849]1048     IF (.NOT. ice_thermo) THEN
[1901]1049        if (iflag_fisrtilp_qsat.lt.1) then
1050           DO i = 1, klon
1051             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
1052           ENDDO
1053        else if (iflag_fisrtilp_qsat.gt.0) then
1054           DO i= 1, klon
[2807]1055    if (fl_cor_ebil .GT. 0) then
1056             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1057    else
[1901]1058             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
[2807]1059    end if
[1901]1060           ENDDO
1061        endif
[1849]1062     ELSE
[2507]1063         if (iflag_t_glace.ge.1) then
[2109]1064            CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
[1901]1065         endif
[2006]1066         if (iflag_fisrtilp_qsat.lt.1) then
1067           DO i = 1, klon
[2109]1068! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
[2086]1069!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
1070!                                     t_glace_max, exposant_glace)
1071              if (iflag_t_glace.eq.0) then
[2223]1072                    zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
[2086]1073                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1074                    zfice(i) = zfice(i)**exposant_glace_old
1075              endif
[2006]1076              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
1077                       +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
1078           ENDDO
1079         else
1080           DO i=1, klon
[2109]1081! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
[2086]1082!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
1083!                                     t_glace_max, exposant_glace)
1084              if (iflag_t_glace.eq.0) then
[2223]1085                    zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
[2086]1086                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1087                    zfice(i) = zfice(i)**exposant_glace_old
1088              endif
[2807]1089        if (fl_cor_ebil .GT. 0) then
1090              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
1091           &             * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
1092                      +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1093        else
[2006]1094              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
[2807]1095                      +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
1096        end if
[2006]1097           ENDDO
1098         endif
1099!         print*,zt(i),zrfl(i),zifl(i),'temp1'
1100       ENDIF
[1849]1101!>AJ
[2807]1102
[2500]1103     ! ----------------------------------------------------------------
1104     ! P3> Formation des precipitations
1105     ! ----------------------------------------------------------------
[1472]1106     !
1107     ! Partager l'eau condensee en precipitation et eau liquide nuageuse
1108     !
[2500]1109
1110     ! Initialisation de zoliq (eau condensee moyenne dans la maille)
[1472]1111     DO i = 1, klon
1112        IF (rneb(i,k).GT.0.0) THEN
1113           zoliq(i) = zcond(i)
1114           zrho(i) = pplay(i,k) / zt(i) / RD
1115           zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
[1849]1116        ENDIF
1117     ENDDO
1118!AJ<
1119     IF (.NOT. ice_thermo) THEN
[2006]1120       IF (iflag_t_glace.EQ.0) THEN
1121         DO i = 1, klon
1122            IF (rneb(i,k).GT.0.0) THEN
1123               zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.13-t_glace_min_old)
1124               zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1125               zfice(i) = zfice(i)**exposant_glace_old
1126!              zfice(i) = zfice(i)**nexpo
1127         !!      zfice(i)=0.
1128            ENDIF
1129         ENDDO
1130       ELSE ! of IF (iflag_t_glace.EQ.0)
[2109]1131         CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
[2086]1132!         DO i = 1, klon
1133!            IF (rneb(i,k).GT.0.0) THEN
[2109]1134! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
[2086]1135!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
1136!                                     t_glace_max, exposant_glace)
1137!            ENDIF
1138!         ENDDO
[2006]1139       ENDIF
[1849]1140     ENDIF
[2500]1141
1142     ! Calcul de radliq (eau condensee pour le rayonnement)
1143     ! Iteration pour realiser une moyenne de l'eau nuageuse lors de la precip
1144     ! Remarque: ce n'est donc pas l'eau restante en fin de precip mais une
1145     ! eau moyenne restante dans le nuage sur la duree du pas de temps qui est
1146     ! transmise au rayonnement;
1147     ! ----------------------------------------------------------------
[1849]1148     DO i = 1, klon
1149        IF (rneb(i,k).GT.0.0) THEN
[1472]1150           zneb(i) = MAX(rneb(i,k), seuil_neb)
[1849]1151     !      zt(i) = zt(i)+zcond(i)*zfice(i)*RLMLT/RCPD/(1.0+RVTMP2*zq(i)) 
1152!           print*,zt(i),'fractionglace'
1153!>AJ
[1472]1154           radliq(i,k) = zoliq(i)/REAL(ninter+1)
1155        ENDIF
1156     ENDDO
1157     !
1158     DO n = 1, ninter
1159        DO i = 1, klon
1160           IF (rneb(i,k).GT.0.0) THEN
1161              zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
[1855]1162              ! Initialization of zpluie and zice:
1163              zpluie=0
1164              zice=0
[1472]1165              IF (zneb(i).EQ.seuil_neb) THEN
1166                 ztot = 0.0
1167              ELSE
[2500]1168                 !  quantite d'eau a eliminer: zchau (Sundqvist, 1978)
1169                 !  meme chose pour la glace: zfroi (Zender & Kiehl, 1997)
[1472]1170                 if (ptconv(i,k)) then
1171                    zcl   =cld_lc_con
1172                    zct   =1./cld_tau_con
1173                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
1174                         *fallvc(zrhol(i)) * zfice(i)
1175                 else
1176                    zcl   =cld_lc_lsc
1177                    zct   =1./cld_tau_lsc
1178                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
1179                         *fallvs(zrhol(i)) * zfice(i)
1180                 endif
[2945]1181
1182                 ! si l'heterogeneite verticale est active, on utilise
1183                 ! la fraction volumique "vraie" plutot que la fraction
1184                 ! surfacique modifiee, qui est plus grande et reduit
1185                 ! sinon l'eau in-cloud de facon artificielle
1186                 if ((iflag_cloudth_vert>=3).AND.(iflag_rain_incloud_vol==1)) then
1187                    zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
1188                         *(1.0-EXP(-(zoliq(i)/ctot_vol(i,k)/zcl   )**2)) *(1.-zfice(i))
1189                 else
1190                    zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
1191                         *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl   )**2)) *(1.-zfice(i))
1192                 endif
[1849]1193!AJ<
1194                 IF (.NOT. ice_thermo) THEN
1195                   ztot    = zchau + zfroi
1196                 ELSE
1197                   zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
1198                   zice = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i))
1199                   ztot    = zpluie    + zice
1200                 ENDIF
1201!>AJ
[1472]1202                 ztot    = MAX(ztot   ,0.0)
1203              ENDIF
1204              ztot    = MIN(ztot,zoliq(i))
[1849]1205!AJ<
1206     !         zoliqp = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie   , 0.0)
1207     !         zoliqi = MAX(zoliq(i)*zfice(i)-1.*zice   , 0.0)
[2807]1208!JLD : les 2 variables zoliqp et zoliqi crorresponent a des pseudo precip
1209!      temporaires et ne doivent pas etre calcule (alors qu'elles le sont
1210!      si iflag_bergeron <> 2
1211!      A SUPPRIMER A TERME
[1849]1212              zoliqp(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie  , 0.0)
1213              zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zice  , 0.0)
[1472]1214              zoliq(i) = MAX(zoliq(i)-ztot   , 0.0)
[1849]1215!>AJ
[1472]1216              radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
1217           ENDIF
[2466]1218        ENDDO  ! i = 1,klon
1219     ENDDO     ! n = 1,ninter
[2807]1220
[2500]1221     ! ----------------------------------------------------------------
[1472]1222     !
[2466]1223     IF (.NOT. ice_thermo) THEN
[1849]1224       DO i = 1, klon
1225         IF (rneb(i,k).GT.0.0) THEN
[1472]1226           d_ql(i,k) = zoliq(i)
1227           zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
1228                * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
[1849]1229         ENDIF
1230       ENDDO
1231     ELSE
[2466]1232!
1233!CR&JYG<
1234! On prend en compte l'effet Bergeron dans les flux de precipitation :
1235! Si T < 0 C, alors les precipitations liquides sont converties en glace, ce qui
1236! provoque un accroissement de temperature DeltaT. L'effet de DeltaT sur le condensat
1237! et les precipitations est grossierement pris en compte en linearisant les equations
1238! et en approximant le processus de precipitation liquide par un processus a seuil.
1239! On fait l'hypothese que le condensat nuageux n'est pas modifié dans cette opération.
1240! Le condensat precipitant liquide est supprime (dans la limite DeltaT<273-T).
1241! Le condensat precipitant solide est augmente.
1242! La vapeur d'eau est augmentee.
1243!
1244       IF ((iflag_bergeron .EQ. 2)) THEN
1245         DO i = 1, klon
1246           IF (rneb(i,k) .GT. 0.0) THEN
1247             zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
1248             zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
[2807]1249           if (fl_cor_ebil .GT. 0) then
1250             zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1251             coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i)
1252!            Calcul de DT si toute les precips liquides congelent
1253             DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
1254!            On ne veut pas que T devienne superieur a la temp. de congelation.
1255!            donc que Delta > RTT-zt(i
1256             DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
1257             zt(i)      = zt(i)      + DeltaT
1258!            Eau vaporisee du fait de l'augmentation de T
1259             Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT
1260!            on reajoute cette eau vaporise a la vapeur et on l'enleve des precips
1261             zq(i) = zq(i) +  Deltaq
1262!            Les 3 max si dessous prtotegent uniquement des erreurs d'arrondies
1263             zcond(i)   = max( zcond(i)- Deltaq, 0. )
1264!            precip liquide qui congele ou qui s'evapore
1265             Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
1266             zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
1267!            bilan eau glacee
1268             zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
1269           else ! if (fl_cor_ebil .GT. 0)
1270!            ancien calcul
[2466]1271             zcp=RCPD*(1.0+RVTMP2*(zq(i)+zcond(i)))
1272             coef1 = RLMLT*zdqs(i)/RLVTT
1273             DeltaT = max( min( RTT-zt(i), RLMLT*zqprecl(i)/zcp/(1.+coef1) ) , 0.)
1274             zqpreci(i) = zqpreci(i) + zcp/RLMLT*DeltaT
1275             zqprecl(i) = max( zqprecl(i) - zcp/RLMLT*(1.+coef1)*DeltaT, 0. )
1276             zcond(i)   = max( zcond(i)   - zcp/RLVTT*zdqs(i)*DeltaT, 0. )
1277             zq(i)      = zq(i)      + zcp/RLVTT*zdqs(i)*DeltaT
1278             zt(i)      = zt(i)      + DeltaT
[2807]1279           end if ! if (fl_cor_ebil .GT. 0)
[2466]1280           ENDIF  ! rneb(i,k) .GT. 0.0
1281         ENDDO
1282         DO i = 1, klon
1283           IF (rneb(i,k).GT.0.0) THEN
1284             d_ql(i,k) = (1-zfice(i))*zoliq(i)
1285             d_qi(i,k) = zfice(i)*zoliq(i)
1286             zrfl(i) = zrfl(i)+ zqprecl(i) &
1287                 *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1288             zifl(i) = zifl(i)+ zqpreci(i) &
1289                      *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
1290           ENDIF                     
1291         ENDDO
1292!!
1293       ELSE  ! iflag_bergeron
1294!>CR&JYG
1295!!
[1849]1296       DO i = 1, klon
1297         IF (rneb(i,k).GT.0.0) THEN
[2086]1298!CR on prend en compte la phase glace
[2807]1299!JLD inutile car on ne passe jamais ici si .not.ice_thermo
1300!           if (.not.ice_thermo) then
1301!           d_ql(i,k) = zoliq(i)
1302!           d_qi(i,k) = 0.
1303!           else
[2086]1304           d_ql(i,k) = (1-zfice(i))*zoliq(i)
1305           d_qi(i,k) = zfice(i)*zoliq(i)
[2807]1306!           endif
[1849]1307!AJ<
1308           zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
1309               *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1310           zifl(i) = zifl(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
1311                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
1312     !      zrfl(i) = zrfl(i)+  zpluie                         &
1313     !          *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1314     !      zifl(i) = zifl(i)+  zice                    &
1315     !               *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)                                   
1316
[2415]1317!CR : on prend en compte l'effet Bergeron dans les flux de precipitation
[2466]1318           IF ((iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)) THEN
[2415]1319              zsolid = zrfl(i)
1320              zifl(i) = zifl(i)+zrfl(i)
1321              zrfl(i) = 0.
[2807]1322           if (fl_cor_ebil .GT. 0) then
[2415]1323              zt(i)=zt(i)+zsolid*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
[2807]1324                      *(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
1325           else
1326              zt(i)=zt(i)+zsolid*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
[2415]1327                      *(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))
[2807]1328           end if
[2466]1329           ENDIF  ! (iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)
[2415]1330!RC   
1331
[2466]1332         ENDIF ! rneb(i,k).GT.0.0
[1849]1333       ENDDO
1334
[2466]1335       ENDIF  ! iflag_bergeron .EQ. 2
1336     ENDIF  ! .NOT. ice_thermo
1337
[2086]1338!CR: la fonte est faite au debut
1339!      IF (ice_thermo) THEN
1340!       DO i = 1, klon
1341!           zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
1342!           zmelt = MIN(MAX(zmelt,0.),1.)
1343!           zrfl(i)=zrfl(i)+zmelt*zifl(i)
1344!           zifl(i)=zifl(i)*(1.-zmelt)
[1849]1345!           print*,zt(i),'octavio1'
[2086]1346!           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
1347!                      *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
[1849]1348!           print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
[2086]1349!       ENDDO
1350!     ENDIF
[1849]1351
1352       
1353     IF (.NOT. ice_thermo) THEN
1354       DO i = 1, klon
1355         IF (zt(i).LT.RTT) THEN
[1472]1356           psfl(i,k)=zrfl(i)
[1849]1357         ELSE
[1472]1358           prfl(i,k)=zrfl(i)
[1849]1359         ENDIF
1360       ENDDO
1361     ELSE
1362     ! JAM*************************************************
[2500]1363     ! Revoir partie ci-dessous: a quoi servent psfl et prfl?
[1849]1364     ! *****************************************************
1365
1366       DO i = 1, klon
1367     !   IF (zt(i).LT.RTT) THEN
1368           psfl(i,k)=zifl(i)
1369     !   ELSE
1370           prfl(i,k)=zrfl(i)
1371     !   ENDIF
1372!>AJ
1373       ENDDO
1374     ENDIF
[2500]1375     ! ----------------------------------------------------------------
1376     ! Fin de formation des precipitations
1377     ! ----------------------------------------------------------------
[1472]1378     !
1379     ! Calculer les tendances de q et de t:
1380     !
1381     DO i = 1, klon
1382        d_q(i,k) = zq(i) - q(i,k)
1383        d_t(i,k) = zt(i) - t(i,k)
1384     ENDDO
1385     !
1386     !AA--------------- Calcul du lessivage stratiforme  -------------
[524]1387
[1472]1388     DO i = 1,klon
1389        !
[1742]1390        if(zcond(i).gt.zoliq(i)+1.e-10) then
1391         beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1392        else
1393         beta(i,k) = 0.
1394        endif
[1472]1395        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
1396             * (paprs(i,k)-paprs(i,k+1))/RG
1397        IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
1398           !AA lessivage nucleation LMD5 dans la couche elle-meme
[2006]1399          IF (iflag_t_glace.EQ.0) THEN
1400           if (t(i,k) .GE. t_glace_min_old) THEN
[1472]1401              zalpha_tr = a_tr_sca(3)
1402           else
1403              zalpha_tr = a_tr_sca(4)
1404           endif
[2006]1405          ELSE ! of IF (iflag_t_glace.EQ.0)
1406           if (t(i,k) .GE. t_glace_min) THEN
1407              zalpha_tr = a_tr_sca(3)
1408           else
1409              zalpha_tr = a_tr_sca(4)
1410           endif
1411          ENDIF
[1472]1412           zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1413           pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
1414           frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
1415           !
1416           ! nucleation avec un facteur -1 au lieu de -0.5
1417           zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
1418           pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
1419        ENDIF
1420        !
1421     ENDDO      ! boucle sur i
1422     !
1423     !AA Lessivage par impaction dans les couches en-dessous
1424     DO kk = k-1, 1, -1
[524]1425        DO i = 1, klon
[1472]1426           IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
[2006]1427             IF (iflag_t_glace.EQ.0) THEN
1428              if (t(i,kk) .GE. t_glace_min_old) THEN
[1472]1429                 zalpha_tr = a_tr_sca(1)
1430              else
1431                 zalpha_tr = a_tr_sca(2)
1432              endif
[2006]1433             ELSE ! of IF (iflag_t_glace.EQ.0)
1434              if (t(i,kk) .GE. t_glace_min) THEN
1435                 zalpha_tr = a_tr_sca(1)
1436              else
1437                 zalpha_tr = a_tr_sca(2)
1438              endif
1439             ENDIF
[1472]1440              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1441              pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
1442              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
1443           ENDIF
[524]1444        ENDDO
[1472]1445     ENDDO
1446     !
[2500]1447     !AA===============================================================
1448     !                     FIN DE LA BOUCLE VERTICALE 
[1472]1449  end DO
1450  !
[2500]1451  !AA==================================================================
[1472]1452  !
1453  ! Pluie ou neige au sol selon la temperature de la 1ere couche
1454  !
[2086]1455!CR: si la thermo de la glace est active, on calcule zifl directement
1456  IF (.NOT.ice_thermo) THEN
[1472]1457  DO i = 1, klon
1458     IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
[1849]1459!AJ<
[2086]1460!        snow(i) = zrfl(i)
[1849]1461        snow(i) = zrfl(i)+zifl(i)
1462!>AJ
[1472]1463        zlh_solid(i) = RLSTT-RLVTT
1464     ELSE
1465        rain(i) = zrfl(i)
1466        zlh_solid(i) = 0.
1467     ENDIF
1468  ENDDO
[2086]1469
1470  ELSE
1471     DO i = 1, klon
1472        snow(i) = zifl(i)
1473        rain(i) = zrfl(i)
1474     ENDDO
1475   
1476   ENDIF
[1472]1477  !
1478  ! For energy conservation : when snow is present, the solification
1479  ! latent heat is considered.
[2086]1480!CR: si thermo de la glace, neige deja prise en compte
1481  IF (.not.ice_thermo) THEN
[1472]1482  DO k = 1, klev
1483     DO i = 1, klon
1484        zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
[2807]1485        zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
[1472]1486        zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
[2807]1487        d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair(i))
[1472]1488     END DO
1489  END DO
[2086]1490  ENDIF
[1472]1491  !
[883]1492
[1472]1493  if (ncoreczq>0) then
[1575]1494     WRITE(lunout,*)'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
[1472]1495  endif
1496
1497END SUBROUTINE fisrtilp
Note: See TracBrowser for help on using the repository browser.