/[lmdze]/trunk/phylmd/fisrtilp.f90
ViewVC logotype

Diff of /trunk/phylmd/fisrtilp.f90

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/phylmd/fisrtilp.f90 revision 73 by guez, Fri Nov 15 17:48:30 2013 UTC trunk/phylmd/fisrtilp.f90 revision 339 by guez, Thu Sep 26 17:08:42 2019 UTC
# Line 4  module fisrtilp_m Line 4  module fisrtilp_m
4    
5  contains  contains
6    
7    SUBROUTINE fisrtilp(dtime, paprs, pplay, t, q, ptconv, ratqs, d_t, d_q, &    SUBROUTINE fisrtilp(paprs, pplay, t, q, ptconv, ratqs, d_t, d_q, d_ql, rneb, &
8         d_ql, rneb, radliq, rain, snow, pfrac_impa, pfrac_nucl, pfrac_1nucl, &         cldliq, rain, snow, pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, &
9         frac_impa, frac_nucl, prfl, psfl, rhcl)         frac_nucl, prfl, psfl, rhcl)
10    
11      ! From phylmd/fisrtilp.F, version 1.2 2004/11/09 16:55:40      ! From phylmd/fisrtilp.F, version 1.2, 2004/11/09 16:55:40
12      ! Author: Z. X. Li (LMD/CNRS), 20 mars 1995      ! First author: Z. X. Li (LMD/CNRS), 20 mars 1995
13    
14      ! Objet : condensation et précipitation stratiforme, schéma de      ! Objet : condensation et pr\'ecipitation stratiforme, sch\'ema de
15      ! nuage, schéma de condensation à grande échelle (pluie).      ! nuage, sch\'ema de condensation \`a grande \'echelle (pluie).
16    
17        USE numer_rec_95, ONLY: nr_erf
18    
19        use comconst, only: dtphys
20        USE comfisrtilp, ONLY: cld_lc_con, cld_lc_lsc, cld_tau_con, &
21             cld_tau_lsc, coef_eva, ffallv_con, ffallv_lsc, iflag_pdf, reevap_ice
22      USE dimphy, ONLY: klev, klon      USE dimphy, ONLY: klev, klon
23        USE fcttre, ONLY: foede, foeew
24      USE suphec_m, ONLY: rcpd, rd, retv, rg, rlstt, rlvtt, rtt      USE suphec_m, ONLY: rcpd, rd, retv, rg, rlstt, rlvtt, rtt
25      USE yoethf_m, ONLY: r2es, r5ies, r5les, rvtmp2      USE yoethf_m, ONLY: r2es, r5ies, r5les, rvtmp2
     USE fcttre, ONLY: dqsatl, dqsats, foede, foeew, qsatl, qsats, thermcep  
     USE comfisrtilp, ONLY: cld_lc_con, cld_lc_lsc, cld_tau_con, &  
          cld_tau_lsc, coef_eva, ffallv_con, ffallv_lsc, iflag_pdf, reevap_ice  
     USE numer_rec_95, ONLY: nr_erf  
26    
27      ! Arguments:      REAL, INTENT (IN):: paprs(klon, klev+1) ! pression a inter-couche
   
     REAL, INTENT (IN):: dtime ! intervalle du temps (s)                  
     REAL, INTENT (IN):: paprs(klon, klev+1) ! pression a inter-couche    
28      REAL, INTENT (IN):: pplay(klon, klev) ! pression au milieu de couche      REAL, INTENT (IN):: pplay(klon, klev) ! pression au milieu de couche
29      REAL, INTENT (IN):: t(klon, klev) ! temperature (K)      REAL, INTENT (IN):: t(klon, klev) ! temperature (K)
30      REAL, INTENT (IN):: q(klon, klev) ! humidite specifique (kg/kg)      REAL, INTENT (IN):: q(klon, klev) ! humidite specifique (kg/kg)
31      LOGICAL ptconv(klon, klev) ! determine la largeur de distribution de vapeur      LOGICAL, INTENT (IN):: ptconv(klon, klev)
32      REAL ratqs(klon, klev) ! determine la largeur de distribution de vapeur  
33      REAL d_t(klon, klev) ! incrementation de la temperature (K)              REAL, INTENT (IN):: ratqs(klon, klev)
34      REAL d_q(klon, klev) ! incrementation de la vapeur d'eau                ! determine la largeur de distribution de vapeur
35      REAL d_ql(klon, klev) ! incrementation de l'eau liquide              
36      REAL rneb(klon, klev) ! fraction nuageuse                                REAL, INTENT (out):: d_t(klon, klev) ! incrementation de la temperature (K)
37      REAL radliq(klon, klev) ! eau liquide utilisee dans rayonnements        REAL, INTENT (out):: d_q(klon, klev) ! incrementation de la vapeur d'eau
38      REAL rain(klon) ! pluies (mm/s)                                        REAL, INTENT (out):: d_ql(klon, klev) ! incrementation de l'eau liquide
39      REAL snow(klon) ! neige (mm/s)                                          REAL, INTENT (out):: rneb(klon, klev) ! fraction nuageuse
40    
41      ! Coeffients de fraction lessivee : pour OFF-LINE      REAL, INTENT (out):: cldliq(klon, klev)
42      REAL pfrac_impa(klon, klev)      ! eau liquide utilisee dans rayonnement
43      REAL pfrac_nucl(klon, klev)  
44      REAL pfrac_1nucl(klon, klev)      REAL, INTENT (out):: rain(klon) ! pluies (mm/s)
45        REAL, INTENT (out):: snow(klon) ! neige (mm/s)
46      ! Fraction d'aerosols lessivee par impaction et par nucleation  
47      ! POur ON-LINE      ! Coeffients de fraction lessivee :
48      REAL frac_nucl(klon, klev)      REAL, INTENT (inout):: pfrac_impa(klon, klev)
49        REAL, INTENT (inout):: pfrac_nucl(klon, klev)
50      REAL prfl(klon, klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)      REAL, INTENT (inout):: pfrac_1nucl(klon, klev)
51      REAL psfl(klon, klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)  
52      REAL rhcl(klon, klev) ! humidite relative en ciel clair                  ! Fraction d'aerosols lessivee par impaction
53        REAL, INTENT (out):: frac_impa(klon, klev)
54    
55        ! Fraction d'aerosols lessivee par nucleation
56        REAL, INTENT (out):: frac_nucl(klon, klev)
57    
58        REAL, INTENT (out):: prfl(klon, klev+1)
59        ! flux d'eau precipitante aux interfaces (kg/m2/s)
60    
61        REAL, INTENT (out):: psfl(klon, klev+1)
62        ! flux d'eau precipitante aux interfaces (kg/m2/s)
63    
64        REAL, INTENT (out):: rhcl(klon, klev) ! humidite relative en ciel clair
65    
66      ! Local:      ! Local:
67    
     ! Fraction d'aerosols lessivee par impaction et par nucleation  
     ! POur ON-LINE  
     REAL frac_impa(klon, klev)  
68      REAL zct(klon), zcl(klon)      REAL zct(klon), zcl(klon)
     !AA  
69    
70      ! Options du programme:      ! Options du programme:
71    
# Line 67  contains Line 74  contains
74    
75      INTEGER ninter ! sous-intervals pour la precipitation      INTEGER ninter ! sous-intervals pour la precipitation
76      PARAMETER (ninter=5)      PARAMETER (ninter=5)
77      LOGICAL evap_prec ! evaporation de la pluie                            LOGICAL evap_prec ! evaporation de la pluie
78      PARAMETER (evap_prec=.TRUE.)      PARAMETER (evap_prec=.TRUE.)
79      REAL zpdf_sig(klon), zpdf_k(klon), zpdf_delta(klon)      REAL zpdf_sig(klon), zpdf_k(klon), zpdf_delta(klon)
80      REAL zpdf_a(klon), zpdf_b(klon), zpdf_e1(klon), zpdf_e2(klon)      REAL zpdf_a(klon), zpdf_b(klon), zpdf_e1(klon), zpdf_e2(klon)
81    
82      LOGICAL cpartiel ! condensation partielle                              LOGICAL cpartiel ! condensation partielle
83      PARAMETER (cpartiel=.TRUE.)      PARAMETER (cpartiel=.TRUE.)
84      REAL t_coup      REAL t_coup
85      PARAMETER (t_coup=234.0)      PARAMETER (t_coup=234.0)
86    
     ! Variables locales:  
   
87      INTEGER i, k, n, kk      INTEGER i, k, n, kk
88      REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5      REAL zqs(klon), zdqs(klon), zcor, zcvm5
89        logical zdelta
90      REAL zrfl(klon), zrfln(klon), zqev, zqevt      REAL zrfl(klon), zrfln(klon), zqev, zqevt
91      REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq      REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
92      REAL ztglace, zt(klon)      REAL ztglace, zt(klon)
93      INTEGER nexpo ! exponentiel pour glace/eau                              INTEGER nexpo ! exponentiel pour glace/eau
94      REAL zdz(klon), zrho(klon), ztot(klon), zrhol(klon)      REAL zdz(klon), zrho(klon), ztot(klon), zrhol(klon)
95      REAL zchau(klon), zfroi(klon), zfice(klon), zneb(klon)      REAL zchau(klon), zfroi(klon), zfice(klon), zneb(klon)
96    
97      LOGICAL appel1er      LOGICAL:: appel1er = .TRUE.
     SAVE appel1er  
   
     !---------------------------------------------------------------  
98    
99      !AA Variables traceurs:      ! Variables traceurs:
100      !AA  Provisoire !!! Parametres alpha du lessivage      ! Provisoire !!! Parametres alpha du lessivage
101      !AA  A priori on a 4 scavenging numbers possibles      ! A priori on a 4 scavenging numbers possibles
102    
103      REAL a_tr_sca(4)      REAL, save:: a_tr_sca(4)
     SAVE a_tr_sca  
104    
105      ! Variables intermediaires      ! Variables intermediaires
106    
107      REAL zalpha_tr      REAL zalpha_tr
108      REAL zfrac_lessi      REAL zfrac_lessi
109      REAL zprec_cond(klon)      REAL zprec_cond(klon)
     !AA  
110      REAL zmair, zcpair, zcpeau      REAL zmair, zcpair, zcpeau
111      !     Pour la conversion eau-neige      ! Pour la conversion eau-neige
112      REAL zlh_solid(klon), zm_solid      REAL zlh_solid(klon), zm_solid
     !IM  
     INTEGER klevm1  
     !---------------------------------------------------------------  
   
     ! Fonctions en ligne:  
113    
114      REAL fallvs, fallvc ! vitesse de chute pour crystaux de glace            !---------------------------------------------------------------
     REAL zzz  
   
     fallvc(zzz) = 3.29/2.0*((zzz)**0.16)*ffallv_con  
     fallvs(zzz) = 3.29/2.0*((zzz)**0.16)*ffallv_lsc  
115    
     DATA appel1er/ .TRUE./  
     !ym  
116      zdelq = 0.0      zdelq = 0.0
117    
118      IF (appel1er) THEN      IF (appel1er) THEN
   
119         PRINT *, 'fisrtilp, ninter:', ninter         PRINT *, 'fisrtilp, ninter:', ninter
120         PRINT *, 'fisrtilp, evap_prec:', evap_prec         PRINT *, 'fisrtilp, evap_prec:', evap_prec
121         PRINT *, 'fisrtilp, cpartiel:', cpartiel         PRINT *, 'fisrtilp, cpartiel:', cpartiel
122         IF (abs(dtime / real(ninter) - 360.) > 0.001) THEN         IF (abs(dtphys / real(ninter) - 360.) > 0.001) THEN
123            PRINT *, "fisrtilp : ce n'est pas prévu, voir Z. X. Li", dtime            PRINT *, "fisrtilp : ce n'est pas pr\'evu, voir Z. X. Li", dtphys
124            PRINT *, 'Je préfère un sous-intervalle de 6 minutes.'            PRINT *, "Je pr\'ef\`ere un sous-intervalle de 6 minutes."
125         END IF         END IF
126         appel1er = .FALSE.         appel1er = .FALSE.
127    
128         !AA initialiation provisoire         ! initialiation provisoire
129         a_tr_sca(1) = -0.5         a_tr_sca(1) = -0.5
130         a_tr_sca(2) = -0.5         a_tr_sca(2) = -0.5
131         a_tr_sca(3) = -0.5         a_tr_sca(3) = -0.5
132         a_tr_sca(4) = -0.5         a_tr_sca(4) = -0.5
133    
134         !AA Initialisation a 1 des coefs des fractions lessivees         ! Initialisation a 1 des coefs des fractions lessivees
   
135         DO k = 1, klev         DO k = 1, klev
136            DO i = 1, klon            DO i = 1, klon
137               pfrac_nucl(i, k) = 1.               pfrac_nucl(i, k) = 1.
# Line 151  contains Line 139  contains
139               pfrac_impa(i, k) = 1.               pfrac_impa(i, k) = 1.
140            END DO            END DO
141         END DO         END DO
142        END IF
143    
144        ! Initialisation a 0 de zoliq
     END IF !  test sur appel1er  
     !MAf Initialisation a 0 de zoliq  
145      DO i = 1, klon      DO i = 1, klon
146         zoliq(i) = 0.         zoliq(i) = 0.
147      END DO      END DO
148      ! Determiner les nuages froids par leur temperature      ! Determiner les nuages froids par leur temperature
149      !  nexpo regle la raideur de la transition eau liquide / eau glace.      ! nexpo regle la raideur de la transition eau liquide / eau glace.
150    
151      ztglace = rtt - 15.0      ztglace = rtt - 15.0
152      nexpo = 6      nexpo = 6
     !cc      nexpo = 1  
153    
154      ! Initialiser les sorties:      ! Initialiser les sorties:
155    
# Line 180  contains Line 166  contains
166            d_q(i, k) = 0.0            d_q(i, k) = 0.0
167            d_ql(i, k) = 0.0            d_ql(i, k) = 0.0
168            rneb(i, k) = 0.0            rneb(i, k) = 0.0
169            radliq(i, k) = 0.0            cldliq(i, k) = 0.0
170            frac_nucl(i, k) = 1.            frac_nucl(i, k) = 1.
171            frac_impa(i, k) = 1.            frac_impa(i, k) = 1.
172         END DO         END DO
# Line 197  contains Line 183  contains
183         zneb(i) = seuil_neb         zneb(i) = seuil_neb
184      END DO      END DO
185    
186        ! Pour plus de securite
     !AA Pour plus de securite  
187    
188      zalpha_tr = 0.      zalpha_tr = 0.
189      zfrac_lessi = 0.      zfrac_lessi = 0.
190    
191      !AA----------------------------------------------------------      loop_vertical: DO k = klev, 1, -1
   
     ! Boucle verticale (du haut vers le bas)  
   
     !IM : klevm1  
     klevm1 = klev - 1  
     DO  k = klev, 1, -1  
   
        !AA----------------------------------------------------------  
   
192         DO i = 1, klon         DO i = 1, klon
193            zt(i) = t(i, k)            zt(i) = t(i, k)
194            zq(i) = q(i, k)            zq(i) = q(i, k)
# Line 225  contains Line 201  contains
201         ! surface.         ! surface.
202    
203         DO i = 1, klon         DO i = 1, klon
204            IF (k<=klevm1) THEN            IF (k <= klev - 1) THEN
205               zmair = (paprs(i, k)-paprs(i, k+1))/rg               zmair = (paprs(i, k)-paprs(i, k+1))/rg
206               zcpair = rcpd*(1.0+rvtmp2*zq(i))               zcpair = rcpd*(1.0+rvtmp2*zq(i))
207               zcpeau = rcpd*rvtmp2               zcpeau = rcpd*rvtmp2
208               zt(i) = ((t(i, k + 1) + d_t(i, k + 1)) * zrfl(i) * dtime &               zt(i) = ((t(i, k + 1) + d_t(i, k + 1)) * zrfl(i) * dtphys &
209                    * zcpeau + zmair * zcpair* zt(i)) &                    * zcpeau + zmair * zcpair* zt(i)) &
210                    / (zmair * zcpair + zrfl(i) * dtime * zcpeau)                    / (zmair * zcpair + zrfl(i) * dtphys * zcpeau)
211            END IF            END IF
212         END DO         END DO
213    
        ! Calculer l'evaporation de la precipitation  
   
   
   
214         IF (evap_prec) THEN         IF (evap_prec) THEN
215              ! Calculer l'evaporation de la precipitation
216            DO i = 1, klon            DO i = 1, klon
217               IF (zrfl(i)>0.) THEN               IF (zrfl(i)>0.) THEN
218                  IF (thermcep) THEN                  zqs(i) = r2es*foeew(zt(i), rtt >= zt(i))/pplay(i, k)
219                     zdelta = max(0., sign(1., rtt-zt(i)))                  zqs(i) = min(0.5, zqs(i))
220                     zqs(i) = r2es*foeew(zt(i), zdelta)/pplay(i, k)                  zcor = 1./(1.-retv*zqs(i))
221                     zqs(i) = min(0.5, zqs(i))                  zqs(i) = zqs(i)*zcor
                    zcor = 1./(1.-retv*zqs(i))  
                    zqs(i) = zqs(i)*zcor  
                 ELSE  
                    IF (zt(i)<t_coup) THEN  
                       zqs(i) = qsats(zt(i))/pplay(i, k)  
                    ELSE  
                       zqs(i) = qsatl(zt(i))/pplay(i, k)  
                    END IF  
                 END IF  
222                  zqev = max(0.0, (zqs(i)-zq(i))*zneb(i))                  zqev = max(0.0, (zqs(i)-zq(i))*zneb(i))
223                  zqevt = coef_eva*(1.0-zq(i)/zqs(i))*sqrt(zrfl(i))* &                  zqevt = coef_eva*(1.0-zq(i)/zqs(i))*sqrt(zrfl(i))* &
224                       (paprs(i, k)-paprs(i, k+1))/pplay(i, k)*zt(i)*rd/rg                       (paprs(i, k)-paprs(i, k+1))/pplay(i, k)*zt(i)*rd/rg
225                  zqevt = max(0.0, min(zqevt, zrfl(i)))*rg*dtime/ &                  zqevt = max(0.0, min(zqevt, zrfl(i)))*rg*dtphys/ &
226                       (paprs(i, k)-paprs(i, k+1))                       (paprs(i, k)-paprs(i, k+1))
227                  zqev = min(zqev, zqevt)                  zqev = min(zqev, zqevt)
228                  zrfln(i) = zrfl(i) - zqev*(paprs(i, k)-paprs(i, k+1))/rg/dtime                  zrfln(i) = zrfl(i) - zqev*(paprs(i, k)-paprs(i, k+1))/rg/dtphys
229    
230                  ! pour la glace, on réévapore toute la précip dans la                  ! pour la glace, on r\'e\'evapore toute la pr\'ecip dans la
231                  ! couche du dessous la glace venant de la couche du                  ! couche du dessous la glace venant de la couche du
232                  ! dessus est simplement dans la couche du dessous.                  ! dessus est simplement dans la couche du dessous.
233    
234                  IF (zt(i)<t_coup .AND. reevap_ice) zrfln(i) = 0.                  IF (zt(i)<t_coup .AND. reevap_ice) zrfln(i) = 0.
235    
236                  zq(i) = zq(i) - (zrfln(i)-zrfl(i))*(rg/(paprs(i, k)-paprs(i, &                  zq(i) = zq(i) - (zrfln(i)-zrfl(i))*(rg/(paprs(i, k)-paprs(i, &
237                       k+1)))*dtime                       k+1)))*dtphys
238                  zt(i) = zt(i) + (zrfln(i)-zrfl(i))*(rg/(paprs(i, k)-paprs(i, &                  zt(i) = zt(i) + (zrfln(i)-zrfl(i))*(rg/(paprs(i, k)-paprs(i, &
239                       k+1)))*dtime*rlvtt/rcpd/(1.0+rvtmp2*zq(i))                       k+1)))*dtphys*rlvtt/rcpd/(1.0+rvtmp2*zq(i))
240                  zrfl(i) = zrfln(i)                  zrfl(i) = zrfln(i)
241               END IF               END IF
242            END DO            END DO
# Line 280  contains Line 244  contains
244    
245         ! Calculer Qs et L/Cp*dQs/dT:         ! Calculer Qs et L/Cp*dQs/dT:
246    
247         IF (thermcep) THEN         DO i = 1, klon
248            DO i = 1, klon            zdelta = rtt >= zt(i)
249               zdelta = max(0., sign(1., rtt-zt(i)))            zcvm5 = merge(r5ies*rlstt, r5les*rlvtt, zdelta)
250               zcvm5 = r5les*rlvtt*(1.-zdelta) + r5ies*rlstt*zdelta            zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zq(i))
251               zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zq(i))            zqs(i) = r2es*foeew(zt(i), zdelta)/pplay(i, k)
252               zqs(i) = r2es*foeew(zt(i), zdelta)/pplay(i, k)            zqs(i) = min(0.5, zqs(i))
253               zqs(i) = min(0.5, zqs(i))            zcor = 1./(1.-retv*zqs(i))
254               zcor = 1./(1.-retv*zqs(i))            zqs(i) = zqs(i)*zcor
255               zqs(i) = zqs(i)*zcor            zdqs(i) = foede(zt(i), zdelta, zcvm5, zqs(i), zcor)
256               zdqs(i) = foede(zt(i), zdelta, zcvm5, zqs(i), zcor)         END DO
           END DO  
        ELSE  
           DO i = 1, klon  
              IF (zt(i)<t_coup) THEN  
                 zqs(i) = qsats(zt(i))/pplay(i, k)  
                 zdqs(i) = dqsats(zt(i), zqs(i))  
              ELSE  
                 zqs(i) = qsatl(zt(i))/pplay(i, k)  
                 zdqs(i) = dqsatl(zt(i), zqs(i))  
              END IF  
           END DO  
        END IF  
257    
258         ! Determiner la condensation partielle et calculer la quantite         ! Determiner la condensation partielle et calculer la quantite
259         ! de l'eau condensee:         ! de l'eau condensee:
260    
261         IF (cpartiel) THEN         IF (cpartiel) THEN
262              ! Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
263              ! nuageuse a partir des PDF de Sandrine Bony.
264              ! rneb : fraction nuageuse
265              ! zqn : eau totale dans le nuage
266              ! zcond : eau condensee moyenne dans la maille.
267    
268            !        print*, 'Dans partiel k=', k            ! on prend en compte le r\'echauffement qui diminue
269              ! la partie condens\'ee
270    
271            !   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau            ! Version avec les ratqs
           !   nuageuse a partir des PDF de Sandrine Bony.  
           !   rneb  : fraction nuageuse  
           !   zqn   : eau totale dans le nuage  
           !   zcond : eau condensee moyenne dans la maille.  
   
           !           on prend en compte le réchauffement qui diminue  
           !           la partie condensee  
   
           !   Version avec les raqts  
272    
273            IF (iflag_pdf==0) THEN            IF (iflag_pdf==0) THEN
   
274               DO i = 1, klon               DO i = 1, klon
275                  zdelq = min(ratqs(i, k), 0.99)*zq(i)                  zdelq = min(ratqs(i, k), 0.99)*zq(i)
276                  rneb(i, k) = (zq(i)+zdelq-zqs(i))/(2.0*zdelq)                  rneb(i, k) = (zq(i)+zdelq-zqs(i))/(2.0*zdelq)
277                  zqn(i) = (zq(i)+zdelq+zqs(i))/2.0                  zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
278               END DO               END DO
   
279            ELSE            ELSE
280                 ! Version avec les nouvelles PDFs.
              !   Version avec les nouvelles PDFs.  
281               DO i = 1, klon               DO i = 1, klon
282                  IF (zq(i)<1.E-15) THEN                  IF (zq(i) < 1E-15) THEN
283                     zq(i) = 1.E-15                     zq(i) = 1E-15
284                  END IF                  END IF
285               END DO               END DO
286               DO i = 1, klon               DO i = 1, klon
# Line 356  contains Line 302  contains
302                     rneb(i, k) = 0.5*zpdf_e1(i)                     rneb(i, k) = 0.5*zpdf_e1(i)
303                     zqn(i) = zq(i)*zpdf_e2(i)/zpdf_e1(i)                     zqn(i) = zq(i)*zpdf_e2(i)/zpdf_e1(i)
304                  END IF                  END IF
   
305               END DO               END DO
   
   
306            END IF            END IF
307            ! iflag_pdf                                                
308            DO i = 1, klon            DO i = 1, klon
309               IF (rneb(i, k)<=0.0) zqn(i) = 0.0               IF (rneb(i, k)<=0.0) zqn(i) = 0.0
310               IF (rneb(i, k)>=1.0) zqn(i) = zq(i)               IF (rneb(i, k)>=1.0) zqn(i) = zq(i)
311               rneb(i, k) = max(0.0, min(1.0, rneb(i, k)))               rneb(i, k) = max(0., min(1., rneb(i, k)))
312               !  On ne divise pas par 1+zdqs pour forcer a avoir l'eau               ! On ne divise pas par 1 + zdqs pour forcer \`a avoir l'eau
313               !  predite par la convection.  ATTENTION !!! Il va               ! pr\'edite par la convection. Attention : il va falloir
314               !  falloir verifier tout ca.               ! verifier tout ca.
315               zcond(i) = max(0.0, zqn(i)-zqs(i))*rneb(i, k)               zcond(i) = max(0., zqn(i)-zqs(i))*rneb(i, k)
              !           print*, 'ZDQS ', zdqs(i)  
              !--Olivier  
316               rhcl(i, k) = (zqs(i)+zq(i)-zdelq)/2./zqs(i)               rhcl(i, k) = (zqs(i)+zq(i)-zdelq)/2./zqs(i)
317               IF (rneb(i, k)<=0.0) rhcl(i, k) = zq(i)/zqs(i)               IF (rneb(i, k) <= 0.) rhcl(i, k) = zq(i) / zqs(i)
318               IF (rneb(i, k)>=1.0) rhcl(i, k) = 1.0               IF (rneb(i, k) >= 1.) rhcl(i, k) = 1.
              !--fin  
319            END DO            END DO
320         ELSE         ELSE
321            DO i = 1, klon            DO i = 1, klon
# Line 390  contains Line 330  contains
330    
331         DO i = 1, klon         DO i = 1, klon
332            zq(i) = zq(i) - zcond(i)            zq(i) = zq(i) - zcond(i)
           !         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD  
333            zt(i) = zt(i) + zcond(i)*rlvtt/rcpd/(1.0+rvtmp2*zq(i))            zt(i) = zt(i) + zcond(i)*rlvtt/rcpd/(1.0+rvtmp2*zq(i))
334         END DO         END DO
335    
# Line 405  contains Line 344  contains
344               zfice(i) = min(max(zfice(i), 0.0), 1.0)               zfice(i) = min(max(zfice(i), 0.0), 1.0)
345               zfice(i) = zfice(i)**nexpo               zfice(i) = zfice(i)**nexpo
346               zneb(i) = max(rneb(i, k), seuil_neb)               zneb(i) = max(rneb(i, k), seuil_neb)
347               radliq(i, k) = zoliq(i)/real(ninter+1)               cldliq(i, k) = zoliq(i)/real(ninter+1)
348            END IF            END IF
349         END DO         END DO
350    
# Line 421  contains Line 360  contains
360                     zcl(i) = cld_lc_lsc                     zcl(i) = cld_lc_lsc
361                     zct(i) = 1./cld_tau_lsc                     zct(i) = 1./cld_tau_lsc
362                  END IF                  END IF
363                  !  quantité d'eau à élminier.                  ! quantit\'e d'eau \`a \'eliminer
364                  zchau(i) = zct(i)*dtime/real(ninter)*zoliq(i)* &                  zchau(i) = zct(i)*dtphys/real(ninter)*zoliq(i)* &
365                       (1.0-exp(-(zoliq(i)/zneb(i)/zcl(i))**2))*(1.-zfice(i))                       (1.0-exp(-(zoliq(i)/zneb(i)/zcl(i))**2))*(1.-zfice(i))
366                  !  meme chose pour la glace.                  ! m\^eme chose pour la glace
367                  IF (ptconv(i, k)) THEN                  IF (ptconv(i, k)) THEN
368                     zfroi(i) = dtime/real(ninter)/zdz(i)*zoliq(i)* &                     zfroi(i) = dtphys/real(ninter)/zdz(i)*zoliq(i)* &
369                          fallvc(zrhol(i))*zfice(i)                          fallvc(zrhol(i))*zfice(i)
370                  ELSE                  ELSE
371                     zfroi(i) = dtime/real(ninter)/zdz(i)*zoliq(i)* &                     zfroi(i) = dtphys/real(ninter)/zdz(i)*zoliq(i)* &
372                          fallvs(zrhol(i))*zfice(i)                          fallvs(zrhol(i))*zfice(i)
373                  END IF                  END IF
374                  ztot(i) = zchau(i) + zfroi(i)                  ztot(i) = zchau(i) + zfroi(i)
375                  IF (zneb(i)==seuil_neb) ztot(i) = 0.0                  IF (zneb(i)==seuil_neb) ztot(i) = 0.0
376                  ztot(i) = min(max(ztot(i), 0.0), zoliq(i))                  ztot(i) = min(max(ztot(i), 0.0), zoliq(i))
377                  zoliq(i) = max(zoliq(i)-ztot(i), 0.0)                  zoliq(i) = max(zoliq(i)-ztot(i), 0.0)
378                  radliq(i, k) = radliq(i, k) + zoliq(i)/real(ninter+1)                  cldliq(i, k) = cldliq(i, k) + zoliq(i)/real(ninter+1)
379               END IF               END IF
380            END DO            END DO
381         END DO         END DO
# Line 445  contains Line 384  contains
384            IF (rneb(i, k)>0.0) THEN            IF (rneb(i, k)>0.0) THEN
385               d_ql(i, k) = zoliq(i)               d_ql(i, k) = zoliq(i)
386               zrfl(i) = zrfl(i) + max(zcond(i) - zoliq(i), 0.) &               zrfl(i) = zrfl(i) + max(zcond(i) - zoliq(i), 0.) &
387                    * (paprs(i, k) - paprs(i, k + 1)) / (rg * dtime)                    * (paprs(i, k) - paprs(i, k + 1)) / (rg * dtphys)
388            END IF            END IF
389            IF (zt(i)<rtt) THEN            IF (zt(i)<rtt) THEN
390               psfl(i, k) = zrfl(i)               psfl(i, k) = zrfl(i)
# Line 454  contains Line 393  contains
393            END IF            END IF
394         END DO         END DO
395    
396         ! Calculer les tendances de q et de t:         ! Calculer les tendances de q et de t :
   
397         DO i = 1, klon         DO i = 1, klon
398            d_q(i, k) = zq(i) - q(i, k)            d_q(i, k) = zq(i) - q(i, k)
399            d_t(i, k) = zt(i) - t(i, k)            d_t(i, k) = zt(i) - t(i, k)
400         END DO         END DO
401    
402         !AA--------------- Calcul du lessivage stratiforme  -------------         ! Calcul du lessivage stratiforme
   
403         DO i = 1, klon         DO i = 1, klon
404            zprec_cond(i) = max(zcond(i)-zoliq(i), 0.0)* &            zprec_cond(i) = max(zcond(i) - zoliq(i), 0.0) &
405                 (paprs(i, k)-paprs(i, k+1))/rg                 * (paprs(i, k)-paprs(i, k+1))/rg
406            IF (rneb(i, k)>0.0 .AND. zprec_cond(i)>0.) THEN            IF (rneb(i, k)>0.0 .AND. zprec_cond(i)>0.) THEN
407               !AA lessivage nucleation LMD5 dans la couche elle-meme               ! lessivage nucleation LMD5 dans la couche elle-meme
408               IF (t(i, k)>=ztglace) THEN               IF (t(i, k)>=ztglace) THEN
409                  zalpha_tr = a_tr_sca(3)                  zalpha_tr = a_tr_sca(3)
410               ELSE               ELSE
# Line 481  contains Line 418  contains
418               zfrac_lessi = 1. - exp(-zprec_cond(i)/zneb(i))               zfrac_lessi = 1. - exp(-zprec_cond(i)/zneb(i))
419               pfrac_1nucl(i, k) = pfrac_1nucl(i, k)*(1.-zneb(i)*zfrac_lessi)               pfrac_1nucl(i, k) = pfrac_1nucl(i, k)*(1.-zneb(i)*zfrac_lessi)
420            END IF            END IF
   
   
421         END DO         END DO
422         !AA Lessivage par impaction dans les couches en-dessous  
423         ! boucle sur i                                                 ! Lessivage par impaction dans les couches en-dessous
424           ! boucle sur i
425         DO kk = k - 1, 1, -1         DO kk = k - 1, 1, -1
426            DO i = 1, klon            DO i = 1, klon
427               IF (rneb(i, k)>0.0 .AND. zprec_cond(i)>0.) THEN               IF (rneb(i, k)>0. .AND. zprec_cond(i)>0.) THEN
428                  IF (t(i, kk)>=ztglace) THEN                  IF (t(i, kk)>=ztglace) THEN
429                     zalpha_tr = a_tr_sca(1)                     zalpha_tr = a_tr_sca(1)
430                  ELSE                  ELSE
# Line 500  contains Line 436  contains
436               END IF               END IF
437            END DO            END DO
438         END DO         END DO
439        end DO loop_vertical
        !AA----------------------------------------------------------  
        !                     FIN DE BOUCLE SUR K  
     end DO  
   
     !AA-----------------------------------------------------------  
440    
441      ! Pluie ou neige au sol selon la temperature de la 1ere couche      ! Pluie ou neige au sol selon la temperature de la 1ere couche
442    
# Line 525  contains Line 456  contains
456         DO i = 1, klon         DO i = 1, klon
457            zcpair = rcpd*(1.0+rvtmp2*(q(i, k)+d_q(i, k)))            zcpair = rcpd*(1.0+rvtmp2*(q(i, k)+d_q(i, k)))
458            zmair = (paprs(i, k)-paprs(i, k+1))/rg            zmair = (paprs(i, k)-paprs(i, k+1))/rg
459            zm_solid = (prfl(i, k)-prfl(i, k+1)+psfl(i, k)-psfl(i, k+1))*dtime            zm_solid = (prfl(i, k)-prfl(i, k+1)+psfl(i, k)-psfl(i, k+1))*dtphys
460            d_t(i, k) = d_t(i, k) + zlh_solid(i)*zm_solid/(zcpair*zmair)            d_t(i, k) = d_t(i, k) + zlh_solid(i)*zm_solid/(zcpair*zmair)
461         END DO         END DO
462      END DO      END DO
463    
464      contains
465    
466        ! vitesse de chute pour cristaux de glace
467    
468        REAL function fallvs(zzz)
469          REAL, intent(in):: zzz
470          fallvs = 3.29/2.0*((zzz)**0.16)*ffallv_lsc
471        end function fallvs
472    
473        !********************************************************
474    
475        real function fallvc(zzz)
476          REAL, intent(in):: zzz
477          fallvc = 3.29/2.0*((zzz)**0.16)*ffallv_con
478        end function fallvc
479    
480    END SUBROUTINE fisrtilp    END SUBROUTINE fisrtilp
481    
482  end module fisrtilp_m  end module fisrtilp_m
Legend:
Removed from v.73  
changed lines
  Added in v.339
  ViewVC Help
Powered by ViewVC 1.1.21