source: trunk/libf/phylmd/fisrtilp.F @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago

Import initial LMDZ5
File size: 18.2 KB
RevLine 
[1]1!
2! $Id: fisrtilp.F 1411 2010-07-09 11:06:15Z 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
52      logical lognormale(klon)
53
54cAA
55c Coeffients de fraction lessivee : pour OFF-LINE
56c
57      REAL pfrac_nucl(klon,klev)
58      REAL pfrac_1nucl(klon,klev)
59      REAL pfrac_impa(klon,klev)
60c
61c Fraction d'aerosols lessivee par impaction et par nucleation
62c POur ON-LINE
63c
64      REAL frac_impa(klon,klev)
65      REAL frac_nucl(klon,klev)
66      real zct      ,zcl
67cAA
68c
69c Options du programme:
70c
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)
87c
88      LOGICAL cpartiel ! condensation partielle
89      PARAMETER (cpartiel=.TRUE.)
90      REAL t_coup
91      PARAMETER (t_coup=234.0)
92c
93c Variables locales:
94c
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)
103c
104      LOGICAL appel1er
105      SAVE appel1er
106c$OMP THREADPRIVATE(appel1er)
107c
108c---------------------------------------------------------------
109c
110cAA Variables traceurs:
111cAA  Provisoire !!! Parametres alpha du lessivage
112cAA  A priori on a 4 scavenging # possibles
113c
114      REAL a_tr_sca(4)
115      save a_tr_sca
116c$OMP THREADPRIVATE(a_tr_sca)
117c
118c Variables intermediaires
119c
120      REAL zalpha_tr
121      REAL zfrac_lessi
122      REAL zprec_cond(klon)
123cAA
124      REAL zmair, zcpair, zcpeau
125C     Pour la conversion eau-neige
126      REAL zlh_solid(klon), zm_solid
127cIM
128cym      INTEGER klevm1
129c---------------------------------------------------------------
130c
131c Fonctions en ligne:
132c
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
139c
140      DATA appel1er /.TRUE./
141cym
142      zdelq=0.0
143     
144      print*,'NUAGES4 A. JAM'
145      IF (appel1er) THEN
146c
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'
153c         CALL abort
154         ENDIF
155         appel1er = .FALSE.
156c
157cAA 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
162c
163cAA Initialisation a 1 des coefs des fractions lessivees
164c
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
175c
176cMAf Initialisation a 0 de zoliq
177c      DO i = 1, klon
178c         zoliq(i)=0.
179c      ENDDO
180c Determiner les nuages froids par leur temperature
181c  nexpo regle la raideur de la transition eau liquide / eau glace.
182c
183      ztglace = RTT - 15.0
184      nexpo = 6
185ccc      nexpo = 1
186c
187c Initialiser les sorties:
188c
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.
213c     ENDDO
214c
215c Initialiser le flux de precipitation a zero
216c
217c     DO i = 1, klon
218         zrfl(i) = 0.0
219         zneb(i) = seuil_neb
220      ENDDO
221c
222c
223cAA Pour plus de securite
224
225      zalpha_tr   = 0.
226      zfrac_lessi = 0.
227
228cAA----------------------------------------------------------
229c
230      ncoreczq=0
231c Boucle verticale (du haut vers le bas)
232c
233cIM : klevm1
234cym      klevm1=klev-1
235      DO 9999 k = klev, 1, -1
236c
237cAA----------------------------------------------------------
238c
239      DO i = 1, klon
240         zt(i)=t(i,k)
241         zq(i)=q(i,k)
242      ENDDO
243c
244c Calculer la varition de temp. de l'air du a la chaleur sensible
245C transporter par la pluie.
246C Il resterait a rajouter cet effet de la chaleur sensible sur les
247C flux de surface, du a la diff. de temp. entre le 1er niveau et la
248C surface.
249C
250      IF(k.LE.klevm1) THEN         
251         DO i = 1, klon
252cIM
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)
259C     C        WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
260         ENDDO
261      ENDIF
262c
263c
264c Calculer l'evaporation de la precipitation
265c
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
293c pour la glace, on r�vapore toute la pr�ip dans la couche du dessous
294c la glace venant de la couche du dessus est simplement dans la couche
295c 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
308c
309c Calculer Qs et L/Cp*dQs/dT:
310c
311      IF (thermcep) THEN
312         DO i = 1, klon
313           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
314           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
315           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
316           zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
317           zqs(i) = MIN(0.5,zqs(i))
318           zcor = 1./(1.-RETV*zqs(i))
319           zqs(i) = zqs(i)*zcor
320           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
333c
334c Determiner la condensation partielle et calculer la quantite
335c de l'eau condensee:
336c
337
338      IF (cpartiel) THEN
339
340c        print*,'Dans partiel k=',k
341c
342c   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
343c   nuageuse a partir des PDF de Sandrine Bony.
344c   rneb  : fraction nuageuse
345c   zqn   : eau totale dans le nuage
346c   zcond : eau condensee moyenne dans la maille.
347c           on prend en compte le r�hauffement qui diminue la partie condensee
348c
349c   Version avec les raqts
350
351         if (iflag_pdf.eq.0) then
352
353           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
357           enddo
358
359         else
360c
361c   Version avec les nouvelles PDFs.
362           do i=1,klon
363              if(zq(i).lt.1.e-15) then
364                ncoreczq=ncoreczq+1
365                zq(i)=1.e-15
366              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
377                 rneb(i,k)=ctot(i,k)
378                 zqn(i)=qcloud(i)
379                 enddo
380
381              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
417
418
419        endif ! iflag_pdf
420
421        DO i=1,klon
422           IF (rneb(i,k) .LE. 0.0) THEN
423              zqn(i) = 0.0
424              rneb(i,k) = 0.0
425              zcond(i) = 0.0
426              rhcl(i,k)=zq(i)/zqs(i)
427           ELSE IF (rneb(i,k) .GE. 1.0) THEN
428              zqn(i) = zq(i)
429              rneb(i,k) = 1.0                 
430              zcond(i) = MAX(0.0,zqn(i)-zqs(i))
431              rhcl(i,k)=1.0
432           ELSE
433              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
434              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
435           ENDIF
436        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
463c
464      DO i = 1, klon
465         zq(i) = zq(i) - zcond(i)
466c         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
467         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
468      ENDDO
469c
470c Partager l'eau condensee en precipitation et eau liquide nuageuse
471c
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
484c
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
493c  quantite d'eau a eliminer: zchau
494c  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
517c
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
530c
531c Calculer les tendances de q et de t:
532c
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
537c
538cAA--------------- Calcul du lessivage stratiforme  -------------
539
540      DO i = 1,klon
541c
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
545cAA 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
554c
555c 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
559c
560      ENDDO      ! boucle sur i
561c
562cAA 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
577c
578cAA----------------------------------------------------------
579c                     FIN DE BOUCLE SUR K   
580 9999 CONTINUE
581c
582cAA-----------------------------------------------------------
583c
584c Pluie ou neige au sol selon la temperature de la 1ere couche
585c
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
595C
596C For energy conservation : when snow is present, the solification
597c 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
606c
607
608      if (ncoreczq>0) then
609         print*,'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
610      endif
611      RETURN
612      END
Note: See TracBrowser for help on using the repository browser.