source: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/fisrtilp.F @ 1643

Last change on this file since 1643 was 1399, checked in by idelkadi, 14 years ago

Modifications relatives a l'inclusion des nouveaux nuages d'Arnaud Jam

  • calltherm.F90 : remontee de variables suplementaires.

Passage de ztv, zpspsk, zthla, zthl en argument

  • fisrtilp.F : le programme de nuages avec appel a cloudth

Ajout des arguments : zqta, fraca, ztv, zpspsk, ztla, ztlh,
iflag_cldcon

  • physiq.F : appel modifie a calltherm et fisrtilp
  • thermcell_main.F90 : changement dans les flag_thermals_ed
  • thermcell_plume.F90 : version modifiee entrainement et detrainement

(on retrouve la precendente pour iflag_thermals_ed=10)

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