source: LMDZ4/trunk/libf/phylmd/fisrtilp.F @ 1283

Last change on this file since 1283 was 1279, checked in by Laurent Fairhead, 15 years ago

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