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

Diff of /trunk/phylmd/fisrtilp.f

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

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

Legend:
Removed from v.10  
changed lines
  Added in v.134

  ViewVC Help
Powered by ViewVC 1.1.21