Ignore:
Timestamp:
Dec 23, 2010, 5:38:42 PM (14 years ago)
Author:
lguez
Message:

Conversion to free source form for "disvert" and "iniacademic", no
other change. Bug fix in "fisrtilp": "fraca" could appear in an
expression while not defined.

File:
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/fisrtilp.F90

    r1463 r1472  
    22! $Id$
    33!
    4 c
    5       SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs,
    6      s                   d_t, d_q, d_ql, rneb, radliq, rain, snow,
    7      s                   pfrac_impa, pfrac_nucl, pfrac_1nucl,
    8      s                   frac_impa, frac_nucl,
    9      s                   prfl, psfl, rhcl, zqta, fraca,
    10      s                   ztv, zpspsk, ztla, zthl, iflag_cldcon)
    11 
    12 c
    13       USE dimphy
    14       IMPLICIT none
    15 c======================================================================
    16 c Auteur(s): Z.X. Li (LMD/CNRS)
    17 c Date: le 20 mars 1995
    18 c Objet: condensation et precipitation stratiforme.
    19 c        schema de nuage
    20 c======================================================================
    21 c======================================================================
    22 cym#include "dimensions.h"
    23 cym#include "dimphy.h"
    24 #include "YOMCST.h"
    25 #include "tracstoke.h"
    26 #include "fisrtilp.h"
    27 c
    28 c Arguments:
    29 c
    30       REAL dtime ! intervalle du temps (s)
    31       REAL paprs(klon,klev+1) ! pression a inter-couche
    32       REAL pplay(klon,klev) ! pression au milieu de couche
    33       REAL t(klon,klev) ! temperature (K)
    34       REAL q(klon,klev) ! humidite specifique (kg/kg)
    35       REAL d_t(klon,klev) ! incrementation de la temperature (K)
    36       REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
    37       REAL d_ql(klon,klev) ! incrementation de l'eau liquide
    38       REAL rneb(klon,klev) ! fraction nuageuse
    39       REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
    40       REAL rhcl(klon,klev) ! humidite relative en ciel clair
    41       REAL rain(klon) ! pluies (mm/s)
    42       REAL snow(klon) ! neige (mm/s)
    43       REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
    44       REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
    45       REAL ztv(klon,klev)
    46       REAL zqta(klon,klev),fraca(klon,klev)
    47       REAL sigma1(klon,klev),sigma2(klon,klev)
    48       REAL qltot(klon,klev),ctot(klon,klev)
    49       REAL zpspsk(klon,klev),ztla(klon,klev)
    50       REAL zthl(klon,klev)
    51 
    52       logical lognormale(klon)
    53 
    54 cAA
    55 c Coeffients de fraction lessivee : pour OFF-LINE
    56 c
    57       REAL pfrac_nucl(klon,klev)
    58       REAL pfrac_1nucl(klon,klev)
    59       REAL pfrac_impa(klon,klev)
    60 c
    61 c Fraction d'aerosols lessivee par impaction et par nucleation
    62 c POur ON-LINE
    63 c
    64       REAL frac_impa(klon,klev)
    65       REAL frac_nucl(klon,klev)
    66       real zct      ,zcl
    67 cAA
    68 c
    69 c Options du programme:
    70 c
    71       REAL seuil_neb ! un nuage existe vraiment au-dela
    72       PARAMETER (seuil_neb=0.001)
    73 
    74       INTEGER ninter ! sous-intervals pour la precipitation
    75       INTEGER ncoreczq 
    76       INTEGER iflag_cldcon
    77       PARAMETER (ninter=5)
    78       LOGICAL evap_prec ! evaporation de la pluie
    79       PARAMETER (evap_prec=.TRUE.)
    80       REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
    81       logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur
    82 
    83       real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
    84       real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
    85       real erf   
    86       REAL qcloud(klon)
    87 c
    88       LOGICAL cpartiel ! condensation partielle
    89       PARAMETER (cpartiel=.TRUE.)
    90       REAL t_coup
    91       PARAMETER (t_coup=234.0)
    92 c
    93 c Variables locales:
    94 c
    95       INTEGER i, k, n, kk
    96       REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5   
    97       REAL zrfl(klon), zrfln(klon), zqev, zqevt
    98       REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
    99       REAL ztglace, zt(klon)
    100       INTEGER nexpo ! exponentiel pour glace/eau
    101       REAL zdz(klon),zrho(klon),ztot      , zrhol(klon)
    102       REAL zchau      ,zfroi      ,zfice(klon),zneb(klon)
    103 c
    104       LOGICAL appel1er
    105       SAVE appel1er
    106 c$OMP THREADPRIVATE(appel1er)
    107 c
    108 c---------------------------------------------------------------
    109 c
    110 cAA Variables traceurs:
    111 cAA  Provisoire !!! Parametres alpha du lessivage
    112 cAA  A priori on a 4 scavenging # possibles
    113 c
    114       REAL a_tr_sca(4)
    115       save a_tr_sca
    116 c$OMP THREADPRIVATE(a_tr_sca)
    117 c
    118 c Variables intermediaires
    119 c
    120       REAL zalpha_tr
    121       REAL zfrac_lessi
    122       REAL zprec_cond(klon)
    123 cAA
    124       REAL zmair, zcpair, zcpeau
    125 C     Pour la conversion eau-neige
    126       REAL zlh_solid(klon), zm_solid
    127 cIM
    128 cym      INTEGER klevm1
    129 c---------------------------------------------------------------
    130 c
    131 c Fonctions en ligne:
    132 c
    133       REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace
    134       REAL zzz
    135 #include "YOETHF.h"
    136 #include "FCTTRE.h"
    137       fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
    138       fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc
    139 c
    140       DATA appel1er /.TRUE./
    141 cym
    142       zdelq=0.0
    143      
    144       print*,'NUAGES4 A. JAM'
    145       IF (appel1er) THEN
    146 c
    147          PRINT*, 'fisrtilp, ninter:', ninter
    148          PRINT*, 'fisrtilp, evap_prec:', evap_prec
    149          PRINT*, 'fisrtilp, cpartiel:', cpartiel
    150          IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
    151           PRINT*, 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
    152           PRINT*, 'Je prefere un sous-intervalle de 6 minutes'
    153 c         CALL abort
    154          ENDIF
    155          appel1er = .FALSE.
    156 c
    157 cAA initialiation provisoire
    158        a_tr_sca(1) = -0.5
    159        a_tr_sca(2) = -0.5
    160        a_tr_sca(3) = -0.5
    161        a_tr_sca(4) = -0.5
    162 c
    163 cAA Initialisation a 1 des coefs des fractions lessivees
    164 c
    165 !cdir collapse
    166       DO k = 1, klev
    167        DO i = 1, klon
    168           pfrac_nucl(i,k)=1.
    169           pfrac_1nucl(i,k)=1.
    170           pfrac_impa(i,k)=1.
    171        ENDDO
    172       ENDDO
    173 
    174       ENDIF          !  test sur appel1er
    175 c
    176 cMAf Initialisation a 0 de zoliq
    177 c      DO i = 1, klon
    178 c         zoliq(i)=0.
    179 c      ENDDO
    180 c Determiner les nuages froids par leur temperature
    181 c  nexpo regle la raideur de la transition eau liquide / eau glace.
    182 c
    183       ztglace = RTT - 15.0
    184       nexpo = 6
    185 ccc      nexpo = 1
    186 c
    187 c Initialiser les sorties:
    188 c
    189 !cdir collapse
    190       DO k = 1, klev+1
    191       DO i = 1, klon
    192          prfl(i,k) = 0.0
    193          psfl(i,k) = 0.0
    194       ENDDO
    195       ENDDO
    196 
    197 !cdir collapse
    198       DO k = 1, klev
    199       DO i = 1, klon
    200          d_t(i,k) = 0.0
    201          d_q(i,k) = 0.0
    202          d_ql(i,k) = 0.0
    203          rneb(i,k) = 0.0
    204          radliq(i,k) = 0.0
    205          frac_nucl(i,k) = 1.
    206          frac_impa(i,k) = 1.
    207       ENDDO
    208       ENDDO
    209       DO i = 1, klon
    210          rain(i) = 0.0
    211          snow(i) = 0.0
    212          zoliq(i)=0.
    213 c     ENDDO
    214 c
    215 c Initialiser le flux de precipitation a zero
    216 c
    217 c     DO i = 1, klon
    218          zrfl(i) = 0.0
    219          zneb(i) = seuil_neb
    220       ENDDO
    221 c
    222 c
    223 cAA Pour plus de securite
    224 
    225       zalpha_tr   = 0.
    226       zfrac_lessi = 0.
    227 
    228 cAA----------------------------------------------------------
    229 c
    230       ncoreczq=0
    231 c Boucle verticale (du haut vers le bas)
    232 c
    233 cIM : klevm1
    234 cym      klevm1=klev-1
    235       DO 9999 k = klev, 1, -1
    236 c
    237 cAA----------------------------------------------------------
    238 c
    239       DO i = 1, klon
    240          zt(i)=t(i,k)
    241          zq(i)=q(i,k)
    242       ENDDO
    243 c
    244 c Calculer la varition de temp. de l'air du a la chaleur sensible
    245 C transporter par la pluie.
    246 C Il resterait a rajouter cet effet de la chaleur sensible sur les
    247 C flux de surface, du a la diff. de temp. entre le 1er niveau et la
    248 C surface.
    249 C
    250       IF(k.LE.klevm1) THEN         
    251          DO i = 1, klon
    252 cIM
    253             zmair=(paprs(i,k)-paprs(i,k+1))/RG
    254             zcpair=RCPD*(1.0+RVTMP2*zq(i))
    255             zcpeau=RCPD*RVTMP2
    256             zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau
    257      $           + zmair*zcpair*zt(i) )
    258      $           / (zmair*zcpair + zrfl(i)*dtime*zcpeau)
    259 C     C        WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
    260          ENDDO
    261       ENDIF
    262 c
    263 c
    264 c Calculer l'evaporation de la precipitation
    265 c
    266 
    267 
    268       IF (evap_prec) THEN
    269       DO i = 1, klon
    270       IF (zrfl(i) .GT.0.) THEN
    271          IF (thermcep) THEN
    272            zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
    273            zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
    274            zqs(i)=MIN(0.5,zqs(i))
    275            zcor=1./(1.-RETV*zqs(i))
    276            zqs(i)=zqs(i)*zcor
    277          ELSE
    278            IF (zt(i) .LT. t_coup) THEN
    279               zqs(i) = qsats(zt(i)) / pplay(i,k)
    280            ELSE
    281               zqs(i) = qsatl(zt(i)) / pplay(i,k)
     4!
     5SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
     6     d_t, d_q, d_ql, rneb, radliq, rain, snow, &
     7     pfrac_impa, pfrac_nucl, pfrac_1nucl, &
     8     frac_impa, frac_nucl, &
     9     prfl, psfl, rhcl, zqta, fraca, &
     10     ztv, zpspsk, ztla, zthl, iflag_cldcon)
     11
     12  !
     13  USE dimphy
     14  IMPLICIT none
     15  !======================================================================
     16  ! Auteur(s): Z.X. Li (LMD/CNRS)
     17  ! Date: le 20 mars 1995
     18  ! Objet: condensation et precipitation stratiforme.
     19  !        schema de nuage
     20  !======================================================================
     21  !======================================================================
     22  !ym include "dimensions.h"
     23  !ym include "dimphy.h"
     24  include "YOMCST.h"
     25  include "tracstoke.h"
     26  include "fisrtilp.h"
     27  !
     28  ! Arguments:
     29  !
     30  REAL dtime ! intervalle du temps (s)
     31  REAL paprs(klon,klev+1) ! pression a inter-couche
     32  REAL pplay(klon,klev) ! pression au milieu de couche
     33  REAL t(klon,klev) ! temperature (K)
     34  REAL q(klon,klev) ! humidite specifique (kg/kg)
     35  REAL d_t(klon,klev) ! incrementation de la temperature (K)
     36  REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
     37  REAL d_ql(klon,klev) ! incrementation de l'eau liquide
     38  REAL rneb(klon,klev) ! fraction nuageuse
     39  REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
     40  REAL rhcl(klon,klev) ! humidite relative en ciel clair
     41  REAL rain(klon) ! pluies (mm/s)
     42  REAL snow(klon) ! neige (mm/s)
     43  REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
     44  REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
     45  REAL ztv(klon,klev)
     46  REAL zqta(klon,klev),fraca(klon,klev)
     47  REAL sigma1(klon,klev),sigma2(klon,klev)
     48  REAL qltot(klon,klev),ctot(klon,klev)
     49  REAL zpspsk(klon,klev),ztla(klon,klev)
     50  REAL zthl(klon,klev)
     51
     52  logical lognormale(klon)
     53
     54  !AA
     55  ! Coeffients de fraction lessivee : pour OFF-LINE
     56  !
     57  REAL pfrac_nucl(klon,klev)
     58  REAL pfrac_1nucl(klon,klev)
     59  REAL pfrac_impa(klon,klev)
     60  !
     61  ! Fraction d'aerosols lessivee par impaction et par nucleation
     62  ! POur ON-LINE
     63  !
     64  REAL frac_impa(klon,klev)
     65  REAL frac_nucl(klon,klev)
     66  real zct      ,zcl
     67  !AA
     68  !
     69  ! Options du programme:
     70  !
     71  REAL seuil_neb ! un nuage existe vraiment au-dela
     72  PARAMETER (seuil_neb=0.001)
     73
     74  INTEGER ninter ! sous-intervals pour la precipitation
     75  INTEGER ncoreczq 
     76  INTEGER iflag_cldcon
     77  PARAMETER (ninter=5)
     78  LOGICAL evap_prec ! evaporation de la pluie
     79  PARAMETER (evap_prec=.TRUE.)
     80  REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
     81  logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur
     82
     83  real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
     84  real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
     85  real erf   
     86  REAL qcloud(klon)
     87  !
     88  LOGICAL cpartiel ! condensation partielle
     89  PARAMETER (cpartiel=.TRUE.)
     90  REAL t_coup
     91  PARAMETER (t_coup=234.0)
     92  !
     93  ! Variables locales:
     94  !
     95  INTEGER i, k, n, kk
     96  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5   
     97  REAL zrfl(klon), zrfln(klon), zqev, zqevt
     98  REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
     99  REAL ztglace, zt(klon)
     100  INTEGER nexpo ! exponentiel pour glace/eau
     101  REAL zdz(klon),zrho(klon),ztot      , zrhol(klon)
     102  REAL zchau      ,zfroi      ,zfice(klon),zneb(klon)
     103  !
     104  LOGICAL appel1er
     105  SAVE appel1er
     106  !$OMP THREADPRIVATE(appel1er)
     107  !
     108  !---------------------------------------------------------------
     109  !
     110  !AA Variables traceurs:
     111  !AA  Provisoire !!! Parametres alpha du lessivage
     112  !AA  A priori on a 4 scavenging # possibles
     113  !
     114  REAL a_tr_sca(4)
     115  save a_tr_sca
     116  !$OMP THREADPRIVATE(a_tr_sca)
     117  !
     118  ! Variables intermediaires
     119  !
     120  REAL zalpha_tr
     121  REAL zfrac_lessi
     122  REAL zprec_cond(klon)
     123  !AA
     124  REAL zmair, zcpair, zcpeau
     125  !     Pour la conversion eau-neige
     126  REAL zlh_solid(klon), zm_solid
     127  !IM
     128  !ym      INTEGER klevm1
     129  !---------------------------------------------------------------
     130  !
     131  ! Fonctions en ligne:
     132  !
     133  REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace
     134  REAL zzz
     135  include "YOETHF.h"
     136  include "FCTTRE.h"
     137  fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
     138  fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc
     139  !
     140  DATA appel1er /.TRUE./
     141  !ym
     142  zdelq=0.0
     143
     144  print*,'NUAGES4 A. JAM'
     145  IF (appel1er) THEN
     146     !
     147     PRINT*, 'fisrtilp, ninter:', ninter
     148     PRINT*, 'fisrtilp, evap_prec:', evap_prec
     149     PRINT*, 'fisrtilp, cpartiel:', cpartiel
     150     IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
     151        PRINT*, 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
     152        PRINT*, 'Je prefere un sous-intervalle de 6 minutes'
     153        !         CALL abort
     154     ENDIF
     155     appel1er = .FALSE.
     156     !
     157     !AA initialiation provisoire
     158     a_tr_sca(1) = -0.5
     159     a_tr_sca(2) = -0.5
     160     a_tr_sca(3) = -0.5
     161     a_tr_sca(4) = -0.5
     162     !
     163     !AA Initialisation a 1 des coefs des fractions lessivees
     164     !
     165     !cdir collapse
     166     DO k = 1, klev
     167        DO i = 1, klon
     168           pfrac_nucl(i,k)=1.
     169           pfrac_1nucl(i,k)=1.
     170           pfrac_impa(i,k)=1.
     171        ENDDO
     172     ENDDO
     173
     174  ENDIF          !  test sur appel1er
     175  !
     176  !MAf Initialisation a 0 de zoliq
     177  !      DO i = 1, klon
     178  !         zoliq(i)=0.
     179  !      ENDDO
     180  ! Determiner les nuages froids par leur temperature
     181  !  nexpo regle la raideur de la transition eau liquide / eau glace.
     182  !
     183  ztglace = RTT - 15.0
     184  nexpo = 6
     185  !cc      nexpo = 1
     186  !
     187  ! Initialiser les sorties:
     188  !
     189  !cdir collapse
     190  DO k = 1, klev+1
     191     DO i = 1, klon
     192        prfl(i,k) = 0.0
     193        psfl(i,k) = 0.0
     194     ENDDO
     195  ENDDO
     196
     197  !cdir collapse
     198  DO k = 1, klev
     199     DO i = 1, klon
     200        d_t(i,k) = 0.0
     201        d_q(i,k) = 0.0
     202        d_ql(i,k) = 0.0
     203        rneb(i,k) = 0.0
     204        radliq(i,k) = 0.0
     205        frac_nucl(i,k) = 1.
     206        frac_impa(i,k) = 1.
     207     ENDDO
     208  ENDDO
     209  DO i = 1, klon
     210     rain(i) = 0.0
     211     snow(i) = 0.0
     212     zoliq(i)=0.
     213     !     ENDDO
     214     !
     215     ! Initialiser le flux de precipitation a zero
     216     !
     217     !     DO i = 1, klon
     218     zrfl(i) = 0.0
     219     zneb(i) = seuil_neb
     220  ENDDO
     221  !
     222  !
     223  !AA Pour plus de securite
     224
     225  zalpha_tr   = 0.
     226  zfrac_lessi = 0.
     227
     228  !AA----------------------------------------------------------
     229  !
     230  ncoreczq=0
     231  ! Boucle verticale (du haut vers le bas)
     232  !
     233  !IM : klevm1
     234  !ym      klevm1=klev-1
     235  DO k = klev, 1, -1
     236     !
     237     !AA----------------------------------------------------------
     238     !
     239     DO i = 1, klon
     240        zt(i)=t(i,k)
     241        zq(i)=q(i,k)
     242     ENDDO
     243     !
     244     ! Calculer la varition de temp. de l'air du a la chaleur sensible
     245     ! transporter par la pluie.
     246     ! Il resterait a rajouter cet effet de la chaleur sensible sur les
     247     ! flux de surface, du a la diff. de temp. entre le 1er niveau et la
     248     ! surface.
     249     !
     250     IF(k.LE.klevm1) THEN         
     251        DO i = 1, klon
     252           !IM
     253           zmair=(paprs(i,k)-paprs(i,k+1))/RG
     254           zcpair=RCPD*(1.0+RVTMP2*zq(i))
     255           zcpeau=RCPD*RVTMP2
     256           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau &
     257                + zmair*zcpair*zt(i) ) &
     258                / (zmair*zcpair + zrfl(i)*dtime*zcpeau)
     259           !     C        WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
     260        ENDDO
     261     ENDIF
     262     !
     263     !
     264     ! Calculer l'evaporation de la precipitation
     265     !
     266
     267
     268     IF (evap_prec) THEN
     269        DO i = 1, klon
     270           IF (zrfl(i) .GT.0.) THEN
     271              IF (thermcep) THEN
     272                 zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
     273                 zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
     274                 zqs(i)=MIN(0.5,zqs(i))
     275                 zcor=1./(1.-RETV*zqs(i))
     276                 zqs(i)=zqs(i)*zcor
     277              ELSE
     278                 IF (zt(i) .LT. t_coup) THEN
     279                    zqs(i) = qsats(zt(i)) / pplay(i,k)
     280                 ELSE
     281                    zqs(i) = qsatl(zt(i)) / pplay(i,k)
     282                 ENDIF
     283              ENDIF
     284              zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
     285              zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
     286                   * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
     287              zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
     288                   * RG*dtime/(paprs(i,k)-paprs(i,k+1))
     289              zqev = MIN (zqev, zqevt)
     290              zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
     291                   /RG/dtime
     292
     293              ! pour la glace, on ré-évapore toute la précip dans la
     294              ! couche du dessous
     295              ! la glace venant de la couche du dessus est simplement
     296              ! dans la couche du dessous.
     297
     298              IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
     299
     300              zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
     301                   * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
     302              zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
     303                   * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
     304                   * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
     305              zrfl(i) = zrfln(i)
    282306           ENDIF
    283          ENDIF
    284          zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
    285          zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i))
    286      .         * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
    287          zqevt = MAX(0.0,MIN(zqevt,zrfl(i)))
    288      .         * RG*dtime/(paprs(i,k)-paprs(i,k+1))
    289          zqev = MIN (zqev, zqevt)
    290          zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))
    291      .                            /RG/dtime
    292 
    293 c pour la glace, on r�vapore toute la pr�ip dans la couche du dessous
    294 c la glace venant de la couche du dessus est simplement dans la couche
    295 c du dessous.
    296 
    297          IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
    298 
    299          zq(i) = zq(i) - (zrfln(i)-zrfl(i))
    300      .             * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
    301          zt(i) = zt(i) + (zrfln(i)-zrfl(i))
    302      .             * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
    303      .             * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
    304          zrfl(i) = zrfln(i)
    305       ENDIF
    306       ENDDO
    307       ENDIF
    308 c
    309 c Calculer Qs et L/Cp*dQs/dT:
    310 c
    311       IF (thermcep) THEN
    312          DO i = 1, klon
     307        ENDDO
     308     ENDIF
     309     !
     310     ! Calculer Qs et L/Cp*dQs/dT:
     311     !
     312     IF (thermcep) THEN
     313        DO i = 1, klon
    313314           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
    314315           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
    … …  
    319320           zqs(i) = zqs(i)*zcor
    320321           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
    321          ENDDO
    322       ELSE
    323          DO i = 1, klon
    324             IF (zt(i).LT.t_coup) THEN
    325                zqs(i) = qsats(zt(i))/pplay(i,k)
    326                zdqs(i) = dqsats(zt(i),zqs(i))
    327             ELSE
    328                zqs(i) = qsatl(zt(i))/pplay(i,k)
    329                zdqs(i) = dqsatl(zt(i),zqs(i))
    330             ENDIF
    331          ENDDO
    332       ENDIF
    333 c
    334 c Determiner la condensation partielle et calculer la quantite
    335 c de l'eau condensee:
    336 c
    337 
    338       IF (cpartiel) THEN
    339 
    340 c        print*,'Dans partiel k=',k
    341 c
    342 c   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
    343 c   nuageuse a partir des PDF de Sandrine Bony.
    344 c   rneb  : fraction nuageuse
    345 c   zqn   : eau totale dans le nuage
    346 c   zcond : eau condensee moyenne dans la maille.
    347 c           on prend en compte le r�hauffement qui diminue la partie condensee
    348 c
    349 c   Version avec les raqts
    350 
    351          if (iflag_pdf.eq.0) then
     322        ENDDO
     323     ELSE
     324        DO i = 1, klon
     325           IF (zt(i).LT.t_coup) THEN
     326              zqs(i) = qsats(zt(i))/pplay(i,k)
     327              zdqs(i) = dqsats(zt(i),zqs(i))
     328           ELSE
     329              zqs(i) = qsatl(zt(i))/pplay(i,k)
     330              zdqs(i) = dqsatl(zt(i),zqs(i))
     331           ENDIF
     332        ENDDO
     333     ENDIF
     334     !
     335     ! Determiner la condensation partielle et calculer la quantite
     336     ! de l'eau condensee:
     337     !
     338
     339     IF (cpartiel) THEN
     340
     341        !        print*,'Dans partiel k=',k
     342        !
     343        !   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
     344        !   nuageuse a partir des PDF de Sandrine Bony.
     345        !   rneb  : fraction nuageuse
     346        !   zqn   : eau totale dans le nuage
     347        !   zcond : eau condensee moyenne dans la maille.
     348        !  on prend en compte le réchauffement qui diminue la partie
     349        ! condensee
     350        !
     351        !   Version avec les raqts
     352
     353        if (iflag_pdf.eq.0) then
    352354
    353355           do i=1,klon
    354             zdelq = min(ratqs(i,k),0.99) * zq(i)
    355             rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
    356             zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
     356              zdelq = min(ratqs(i,k),0.99) * zq(i)
     357              rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
     358              zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
    357359           enddo
    358360
    359          else
    360 c
    361 c   Version avec les nouvelles PDFs.
     361        else
     362           !
     363           !   Version avec les nouvelles PDFs.
    362364           do i=1,klon
    363365              if(zq(i).lt.1.e-15) then
    364                 ncoreczq=ncoreczq+1
    365                 zq(i)=1.e-15
     366                 ncoreczq=ncoreczq+1
     367                 zq(i)=1.e-15
    366368              endif
    367            enddo 
    368 
    369               if (iflag_cldcon>=5) then
    370 
    371                  call cloudth(klon,klev,k,ztv,
    372      .           zq,zqta,fraca,
    373      .           qcloud,ctot,zpspsk,paprs,ztla,zthl,
    374      .           ratqs,zqs,t)
    375 
    376                  do i=1,klon
     369           enddo
     370
     371           if (iflag_cldcon>=5) then
     372
     373              call cloudth(klon,klev,k,ztv, &
     374                   zq,zqta,fraca, &
     375                   qcloud,ctot,zpspsk,paprs,ztla,zthl, &
     376                   ratqs,zqs,t)
     377
     378              do i=1,klon
    377379                 rneb(i,k)=ctot(i,k)
    378380                 zqn(i)=qcloud(i)
    379                  enddo
    380 
     381              enddo
     382
     383           endif
     384
     385           if (iflag_cldcon <= 4) then
     386              lognormale = .true.
     387           elseif (iflag_cldcon == 6) then
     388              ! lognormale en l'absence des thermiques
     389              lognormale = fraca(:,k) < 1e-10
     390           else
     391              ! Dans le cas iflag_cldcon=5, on prend systématiquement la
     392              ! bi-gaussienne
     393              lognormale = .false.
     394           end if
     395
     396           do i=1,klon
     397              if (lognormale(i)) then
     398                 zpdf_sig(i)=ratqs(i,k)*zq(i)
     399                 zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
     400                 zpdf_delta(i)=log(zq(i)/zqs(i))
     401                 zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
     402                 zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
     403                 zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
     404                 zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
     405                 zpdf_e1(i)=1.-erf(zpdf_e1(i))
     406                 zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
     407                 zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
     408                 zpdf_e2(i)=1.-erf(zpdf_e2(i))
    381409              endif
    382 
    383 ! Pour iflag_cldcon<=4, on prend toujours la lognormale
    384 ! Dans le cas iflag_cldcon=5, on prend systématiquement la bi-gaussienne
    385 ! Dans le cas iflagÃ_cldcon=6, on prend la lognormale en absence des thermiques
    386 
    387             lognormale(:)=
    388      .      iflag_cldcon<=4.or.(iflag_cldcon==6.and.fraca(:,k)<1.e-10)
    389             do i=1,klon
    390             if (lognormale(i)) then
    391             zpdf_sig(i)=ratqs(i,k)*zq(i)
    392             zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
    393             zpdf_delta(i)=log(zq(i)/zqs(i))
    394             zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
    395             zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
    396             zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
    397             zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
    398             zpdf_e1(i)=1.-erf(zpdf_e1(i))
    399             zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
    400             zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
    401             zpdf_e2(i)=1.-erf(zpdf_e2(i))
    402             endif
    403             enddo
    404 
    405             do i=1,klon
    406             if (lognormale(i)) then
    407             if (zpdf_e1(i).lt.1.e-10) then
    408                rneb(i,k)=0.
    409                zqn(i)=zqs(i)
    410             else
    411                rneb(i,k)=0.5*zpdf_e1(i)
    412                zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
    413             endif
    414             endif
    415            
    416            enddo
     410           enddo
     411
     412           do i=1,klon
     413              if (lognormale(i)) then
     414                 if (zpdf_e1(i).lt.1.e-10) then
     415                    rneb(i,k)=0.
     416                    zqn(i)=zqs(i)
     417                 else
     418                    rneb(i,k)=0.5*zpdf_e1(i)
     419                    zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
     420                 endif
     421              endif
     422
     423           enddo
    417424
    418425
    … …  
    435442           ENDIF
    436443        ENDDO
    437 !         do i=1,klon
    438 !            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
    439 !            IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
    440 !            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
    441 !c           zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
    442 !c  On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
    443 !c  la convection.
    444 !c  ATTENTION !!! Il va falloir verifier tout ca.
    445 !            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
    446 !c           print*,'ZDQS ',zdqs(i)
    447 !c--Olivier
    448 !            rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
    449 !            IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
    450 !            IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
    451 !c--fin
    452 !           ENDDO
    453       ELSE
    454          DO i = 1, klon
    455             IF (zq(i).GT.zqs(i)) THEN
    456                rneb(i,k) = 1.0
    457             ELSE
    458                rneb(i,k) = 0.0
    459             ENDIF
    460             zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
    461          ENDDO
    462       ENDIF
    463 c
    464       DO i = 1, klon
    465          zq(i) = zq(i) - zcond(i)
    466 c         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
    467          zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
    468       ENDDO
    469 c
    470 c Partager l'eau condensee en precipitation et eau liquide nuageuse
    471 c
    472       DO i = 1, klon
    473       IF (rneb(i,k).GT.0.0) THEN
    474          zoliq(i) = zcond(i)
    475          zrho(i) = pplay(i,k) / zt(i) / RD
    476          zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
    477          zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
    478          zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
    479          zfice(i) = zfice(i)**nexpo
    480          zneb(i) = MAX(rneb(i,k), seuil_neb)
    481          radliq(i,k) = zoliq(i)/REAL(ninter+1)
    482       ENDIF
    483       ENDDO
    484 c
    485       DO n = 1, ninter
    486       DO i = 1, klon
    487       IF (rneb(i,k).GT.0.0) THEN
    488          zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
    489 
    490          IF (zneb(i).EQ.seuil_neb) THEN
    491              ztot = 0.0
    492          ELSE
    493 c  quantite d'eau a eliminer: zchau
    494 c  meme chose pour la glace: zfroi
    495              if (ptconv(i,k)) then
    496                 zcl   =cld_lc_con
    497                 zct   =1./cld_tau_con
    498                 zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i)
    499      .              *fallvc(zrhol(i)) * zfice(i)
    500              else
    501                 zcl   =cld_lc_lsc
    502                 zct   =1./cld_tau_lsc
    503                 zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i)
    504      .              *fallvs(zrhol(i)) * zfice(i)
    505              endif
    506              zchau    = zct   *dtime/REAL(ninter) * zoliq(i)
    507      .         *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl   )**2)) *(1.-zfice(i))
    508              ztot    = zchau    + zfroi
    509              ztot    = MAX(ztot   ,0.0)
    510          ENDIF
    511          ztot    = MIN(ztot,zoliq(i))
    512          zoliq(i) = MAX(zoliq(i)-ztot   , 0.0)
    513          radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
    514       ENDIF
    515       ENDDO
    516       ENDDO
    517 c
    518       DO i = 1, klon
    519       IF (rneb(i,k).GT.0.0) THEN
    520          d_ql(i,k) = zoliq(i)
    521          zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0)
    522      .                    * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
    523       ENDIF
    524       IF (zt(i).LT.RTT) THEN
    525         psfl(i,k)=zrfl(i)
    526       ELSE
    527         prfl(i,k)=zrfl(i)
    528       ENDIF
    529       ENDDO
    530 c
    531 c Calculer les tendances de q et de t:
    532 c
    533       DO i = 1, klon
    534          d_q(i,k) = zq(i) - q(i,k)
    535          d_t(i,k) = zt(i) - t(i,k)
    536       ENDDO
    537 c
    538 cAA--------------- Calcul du lessivage stratiforme  -------------
    539 
    540       DO i = 1,klon
    541 c
    542          zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)
    543      .                * (paprs(i,k)-paprs(i,k+1))/RG
    544          IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
    545 cAA lessivage nucleation LMD5 dans la couche elle-meme
    546             if (t(i,k) .GE. ztglace) THEN
    547                zalpha_tr = a_tr_sca(3)
    548             else
    549                zalpha_tr = a_tr_sca(4)
    550             endif
    551             zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
    552             pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
    553             frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
    554 c
    555 c nucleation avec un facteur -1 au lieu de -0.5
    556             zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
    557             pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
    558          ENDIF
    559 c
    560       ENDDO      ! boucle sur i
    561 c
    562 cAA Lessivage par impaction dans les couches en-dessous
    563       DO kk = k-1, 1, -1
    564         DO i = 1, klon
    565           IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
    566             if (t(i,kk) .GE. ztglace) THEN
    567               zalpha_tr = a_tr_sca(1)
    568             else
    569               zalpha_tr = a_tr_sca(2)
    570             endif
    571             zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
    572             pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
    573             frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
    574           ENDIF
    575         ENDDO
    576       ENDDO
    577 c
    578 cAA----------------------------------------------------------
    579 c                     FIN DE BOUCLE SUR K   
    580  9999 CONTINUE
    581 c
    582 cAA-----------------------------------------------------------
    583 c
    584 c Pluie ou neige au sol selon la temperature de la 1ere couche
    585 c
    586       DO i = 1, klon
    587       IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
    588          snow(i) = zrfl(i)
    589          zlh_solid(i) = RLSTT-RLVTT
    590       ELSE
    591          rain(i) = zrfl(i)
    592          zlh_solid(i) = 0.
    593       ENDIF
    594       ENDDO
    595 C
    596 C For energy conservation : when snow is present, the solification
    597 c latent heat is considered.
    598       DO k = 1, klev
    599         DO i = 1, klon
    600           zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
    601           zmair=(paprs(i,k)-paprs(i,k+1))/RG
    602           zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
    603           d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair)
    604         END DO
    605       END DO
    606 c
    607 
    608       if (ncoreczq>0) then
    609          print*,'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
    610       endif
    611       RETURN
    612       END
     444        !         do i=1,klon
     445        !            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
     446        !            IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
     447        !            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
     448        !c           zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
     449        !c  On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
     450        !c  la convection.
     451        !c  ATTENTION !!! Il va falloir verifier tout ca.
     452        !            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
     453        !c           print*,'ZDQS ',zdqs(i)
     454        !c--Olivier
     455        !            rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
     456        !            IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
     457        !            IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
     458        !c--fin
     459        !           ENDDO
     460     ELSE
     461        DO i = 1, klon
     462           IF (zq(i).GT.zqs(i)) THEN
     463              rneb(i,k) = 1.0
     464           ELSE
     465              rneb(i,k) = 0.0
     466           ENDIF
     467           zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
     468        ENDDO
     469     ENDIF
     470     !
     471     DO i = 1, klon
     472        zq(i) = zq(i) - zcond(i)
     473        !         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
     474        zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
     475     ENDDO
     476     !
     477     ! Partager l'eau condensee en precipitation et eau liquide nuageuse
     478     !
     479     DO i = 1, klon
     480        IF (rneb(i,k).GT.0.0) THEN
     481           zoliq(i) = zcond(i)
     482           zrho(i) = pplay(i,k) / zt(i) / RD
     483           zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
     484           zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
     485           zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
     486           zfice(i) = zfice(i)**nexpo
     487           zneb(i) = MAX(rneb(i,k), seuil_neb)
     488           radliq(i,k) = zoliq(i)/REAL(ninter+1)
     489        ENDIF
     490     ENDDO
     491     !
     492     DO n = 1, ninter
     493        DO i = 1, klon
     494           IF (rneb(i,k).GT.0.0) THEN
     495              zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
     496
     497              IF (zneb(i).EQ.seuil_neb) THEN
     498                 ztot = 0.0
     499              ELSE
     500                 !  quantite d'eau a eliminer: zchau
     501                 !  meme chose pour la glace: zfroi
     502                 if (ptconv(i,k)) then
     503                    zcl   =cld_lc_con
     504                    zct   =1./cld_tau_con
     505                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
     506                         *fallvc(zrhol(i)) * zfice(i)
     507                 else
     508                    zcl   =cld_lc_lsc
     509                    zct   =1./cld_tau_lsc
     510                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
     511                         *fallvs(zrhol(i)) * zfice(i)
     512                 endif
     513                 zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
     514                      *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl   )**2)) *(1.-zfice(i))
     515                 ztot    = zchau    + zfroi
     516                 ztot    = MAX(ztot   ,0.0)
     517              ENDIF
     518              ztot    = MIN(ztot,zoliq(i))
     519              zoliq(i) = MAX(zoliq(i)-ztot   , 0.0)
     520              radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
     521           ENDIF
     522        ENDDO
     523     ENDDO
     524     !
     525     DO i = 1, klon
     526        IF (rneb(i,k).GT.0.0) THEN
     527           d_ql(i,k) = zoliq(i)
     528           zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
     529                * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
     530        ENDIF
     531        IF (zt(i).LT.RTT) THEN
     532           psfl(i,k)=zrfl(i)
     533        ELSE
     534           prfl(i,k)=zrfl(i)
     535        ENDIF
     536     ENDDO
     537     !
     538     ! Calculer les tendances de q et de t:
     539     !
     540     DO i = 1, klon
     541        d_q(i,k) = zq(i) - q(i,k)
     542        d_t(i,k) = zt(i) - t(i,k)
     543     ENDDO
     544     !
     545     !AA--------------- Calcul du lessivage stratiforme  -------------
     546
     547     DO i = 1,klon
     548        !
     549        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
     550             * (paprs(i,k)-paprs(i,k+1))/RG
     551        IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
     552           !AA lessivage nucleation LMD5 dans la couche elle-meme
     553           if (t(i,k) .GE. ztglace) THEN
     554              zalpha_tr = a_tr_sca(3)
     555           else
     556              zalpha_tr = a_tr_sca(4)
     557           endif
     558           zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
     559           pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
     560           frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
     561           !
     562           ! nucleation avec un facteur -1 au lieu de -0.5
     563           zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
     564           pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
     565        ENDIF
     566        !
     567     ENDDO      ! boucle sur i
     568     !
     569     !AA Lessivage par impaction dans les couches en-dessous
     570     DO kk = k-1, 1, -1
     571        DO i = 1, klon
     572           IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
     573              if (t(i,kk) .GE. ztglace) THEN
     574                 zalpha_tr = a_tr_sca(1)
     575              else
     576                 zalpha_tr = a_tr_sca(2)
     577              endif
     578              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
     579              pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
     580              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
     581           ENDIF
     582        ENDDO
     583     ENDDO
     584     !
     585     !AA----------------------------------------------------------
     586     !                     FIN DE BOUCLE SUR K   
     587  end DO
     588  !
     589  !AA-----------------------------------------------------------
     590  !
     591  ! Pluie ou neige au sol selon la temperature de la 1ere couche
     592  !
     593  DO i = 1, klon
     594     IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
     595        snow(i) = zrfl(i)
     596        zlh_solid(i) = RLSTT-RLVTT
     597     ELSE
     598        rain(i) = zrfl(i)
     599        zlh_solid(i) = 0.
     600     ENDIF
     601  ENDDO
     602  !
     603  ! For energy conservation : when snow is present, the solification
     604  ! latent heat is considered.
     605  DO k = 1, klev
     606     DO i = 1, klon
     607        zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
     608        zmair=(paprs(i,k)-paprs(i,k+1))/RG
     609        zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
     610        d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair)
     611     END DO
     612  END DO
     613  !
     614
     615  if (ncoreczq>0) then
     616     print*,'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
     617  endif
     618
     619END SUBROUTINE fisrtilp
Note: See TracChangeset for help on using the changeset viewer.