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

Last change on this file since 1168 was 1146, checked in by Laurent Fairhead, 16 years ago

Réintegration dans le tronc des modifications issues de la branche LMDZ-dev
comprises entre la révision 1074 et 1145
Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes
résultats entre la trunk et la dev
LF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.8 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)
56      real zct(klon),zcl(klon)
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
89      REAL zdz(klon),zrho(klon),ztot(klon), zrhol(klon)
90      REAL zchau(klon),zfroi(klon),zfice(klon),zneb(klon)
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
152      DO k = 1, klev
153       DO i = 1, klon
154          pfrac_nucl(i,k)=1.
155          pfrac_1nucl(i,k)=1.
156          pfrac_impa(i,k)=1.
157       ENDDO
158      ENDDO
159
160      ENDIF          !  test sur appel1er
161c
162cMAf Initialisation a 0 de zoliq
163       DO i = 1, klon
164          zoliq(i)=0.
165       ENDDO
166c Determiner les nuages froids par leur temperature
167c  nexpo regle la raideur de la transition eau liquide / eau glace.
168c
169      ztglace = RTT - 15.0
170      nexpo = 6
171ccc      nexpo = 1
172c
173c Initialiser les sorties:
174c
175      DO k = 1, klev+1
176      DO i = 1, klon
177         prfl(i,k) = 0.0
178         psfl(i,k) = 0.0
179      ENDDO
180      ENDDO
181
182      DO k = 1, klev
183      DO i = 1, klon
184         d_t(i,k) = 0.0
185         d_q(i,k) = 0.0
186         d_ql(i,k) = 0.0
187         rneb(i,k) = 0.0
188         radliq(i,k) = 0.0
189         frac_nucl(i,k) = 1.
190         frac_impa(i,k) = 1.
191      ENDDO
192      ENDDO
193      DO i = 1, klon
194         rain(i) = 0.0
195         snow(i) = 0.0
196      ENDDO
197c
198c Initialiser le flux de precipitation a zero
199c
200      DO i = 1, klon
201         zrfl(i) = 0.0
202         zneb(i) = seuil_neb
203      ENDDO
204c
205c
206cAA Pour plus de securite
207
208      zalpha_tr   = 0.
209      zfrac_lessi = 0.
210
211cAA----------------------------------------------------------
212c
[883]213      ncoreczq=0
[524]214c Boucle verticale (du haut vers le bas)
215c
216cIM : klevm1
[766]217cym      klevm1=klev-1
[524]218      DO 9999 k = klev, 1, -1
219c
220cAA----------------------------------------------------------
221c
222      DO i = 1, klon
223         zt(i)=t(i,k)
224         zq(i)=q(i,k)
225      ENDDO
226c
227c Calculer la varition de temp. de l'air du a la chaleur sensible
228C transporter par la pluie.
229C Il resterait a rajouter cet effet de la chaleur sensible sur les
230C flux de surface, du a la diff. de temp. entre le 1er niveau et la
231C surface.
232C
[1146]233      IF(k.LE.klevm1) THEN         
234         DO i = 1, klon
[524]235cIM
[1146]236            zmair=(paprs(i,k)-paprs(i,k+1))/RG
237            zcpair=RCPD*(1.0+RVTMP2*zq(i))
238            zcpeau=RCPD*RVTMP2
239            zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau
240     $           + zmair*zcpair*zt(i) )
241     $           / (zmair*zcpair + zrfl(i)*dtime*zcpeau)
242C     C        WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
243         ENDDO
244      ENDIF
[524]245c
246c
247c Calculer l'evaporation de la precipitation
248c
249
250
251      IF (evap_prec) THEN
252      DO i = 1, klon
253      IF (zrfl(i) .GT.0.) THEN
254         IF (thermcep) THEN
255           zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
256           zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
257           zqs(i)=MIN(0.5,zqs(i))
258           zcor=1./(1.-RETV*zqs(i))
259           zqs(i)=zqs(i)*zcor
260         ELSE
261           IF (zt(i) .LT. t_coup) THEN
262              zqs(i) = qsats(zt(i)) / pplay(i,k)
263           ELSE
264              zqs(i) = qsatl(zt(i)) / pplay(i,k)
265           ENDIF
266         ENDIF
267         zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
268         zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i))
269     .         * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
270         zqevt = MAX(0.0,MIN(zqevt,zrfl(i)))
271     .         * RG*dtime/(paprs(i,k)-paprs(i,k+1))
272         zqev = MIN (zqev, zqevt)
273         zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))
274     .                            /RG/dtime
275
[766]276c pour la glace, on r�vapore toute la pr�ip dans la couche du dessous
[524]277c la glace venant de la couche du dessus est simplement dans la couche
278c du dessous.
279
280         IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
281
282         zq(i) = zq(i) - (zrfln(i)-zrfl(i))
283     .             * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
284         zt(i) = zt(i) + (zrfln(i)-zrfl(i))
285     .             * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
286     .             * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
287         zrfl(i) = zrfln(i)
288      ENDIF
289      ENDDO
290      ENDIF
291c
292c Calculer Qs et L/Cp*dQs/dT:
293c
294      IF (thermcep) THEN
295         DO i = 1, klon
296           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
297           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
298           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
299           zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
300           zqs(i) = MIN(0.5,zqs(i))
301           zcor = 1./(1.-RETV*zqs(i))
302           zqs(i) = zqs(i)*zcor
303           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
304         ENDDO
305      ELSE
306         DO i = 1, klon
307            IF (zt(i).LT.t_coup) THEN
308               zqs(i) = qsats(zt(i))/pplay(i,k)
309               zdqs(i) = dqsats(zt(i),zqs(i))
310            ELSE
311               zqs(i) = qsatl(zt(i))/pplay(i,k)
312               zdqs(i) = dqsatl(zt(i),zqs(i))
313            ENDIF
314         ENDDO
315      ENDIF
316c
317c Determiner la condensation partielle et calculer la quantite
318c de l'eau condensee:
319c
320      IF (cpartiel) THEN
321
322c        print*,'Dans partiel k=',k
323c
324c   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
325c   nuageuse a partir des PDF de Sandrine Bony.
326c   rneb  : fraction nuageuse
327c   zqn   : eau totale dans le nuage
328c   zcond : eau condensee moyenne dans la maille.
[766]329c           on prend en compte le r�hauffement qui diminue la partie condensee
[524]330c
331c   Version avec les raqts
332
333         if (iflag_pdf.eq.0) then
334
335           do i=1,klon
336            zdelq = min(ratqs(i,k),0.99) * zq(i)
337            rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
338            zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
339           enddo
340
341         else
342c
343c   Version avec les nouvelles PDFs.
344           do i=1,klon
345              if(zq(i).lt.1.e-15) then
[883]346                ncoreczq=ncoreczq+1
[524]347                zq(i)=1.e-15
348              endif
349           enddo
350           do i=1,klon
351            zpdf_sig(i)=ratqs(i,k)*zq(i)
352            zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
353            zpdf_delta(i)=log(zq(i)/zqs(i))
354            zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
355            zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
356            zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
357            zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
358            zpdf_e1(i)=1.-erf(zpdf_e1(i))
359            zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
360            zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
361            zpdf_e2(i)=1.-erf(zpdf_e2(i))
362            if (zpdf_e1(i).lt.1.e-10) then
363               rneb(i,k)=0.
364               zqn(i)=zqs(i)
365            else
366               rneb(i,k)=0.5*zpdf_e1(i)
367               zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
368            endif
369           
370           enddo
371
372        endif ! iflag_pdf
373
[1146]374        DO i=1,klon
375           IF (rneb(i,k) .LE. 0.0) THEN
376              zqn(i) = 0.0
377              rneb(i,k) = 0.0
378              zcond(i) = 0.0
379              rhcl(i,k)=zq(i)/zqs(i)
380           ELSE IF (rneb(i,k) .GE. 1.0) THEN
381              zqn(i) = zq(i)
382              rneb(i,k) = 1.0                 
383              zcond(i) = MAX(0.0,zqn(i)-zqs(i))
384              rhcl(i,k)=1.0
385           ELSE
386              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
387              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
388           ENDIF
389        ENDDO
390!         do i=1,klon
391!            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
392!            IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
393!            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
394!c           zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
395!c  On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
396!c  la convection.
397!c  ATTENTION !!! Il va falloir verifier tout ca.
398!            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
399!c           print*,'ZDQS ',zdqs(i)
400!c--Olivier
401!            rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
402!            IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
403!            IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
404!c--fin
405!           ENDDO
[524]406      ELSE
407         DO i = 1, klon
408            IF (zq(i).GT.zqs(i)) THEN
409               rneb(i,k) = 1.0
410            ELSE
411               rneb(i,k) = 0.0
412            ENDIF
413            zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
414         ENDDO
415      ENDIF
416c
417      DO i = 1, klon
418         zq(i) = zq(i) - zcond(i)
419c         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
420         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
421      ENDDO
422c
423c Partager l'eau condensee en precipitation et eau liquide nuageuse
424c
425      DO i = 1, klon
426      IF (rneb(i,k).GT.0.0) THEN
427         zoliq(i) = zcond(i)
428         zrho(i) = pplay(i,k) / zt(i) / RD
429         zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
430         zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
431         zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
432         zfice(i) = zfice(i)**nexpo
433         zneb(i) = MAX(rneb(i,k), seuil_neb)
434         radliq(i,k) = zoliq(i)/FLOAT(ninter+1)
435      ENDIF
436      ENDDO
437c
438      DO n = 1, ninter
439      DO i = 1, klon
440      IF (rneb(i,k).GT.0.0) THEN
441         zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
442
443         if (ptconv(i,k)) then
444            zcl(i)=cld_lc_con
445            zct(i)=1./cld_tau_con
446         else
447            zcl(i)=cld_lc_lsc
448            zct(i)=1./cld_tau_lsc
449         endif
[766]450c  quantit�d'eau ��minier.
[524]451         zchau(i) = zct(i)*dtime/FLOAT(ninter) * zoliq(i)
452     .         *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl(i))**2)) *(1.-zfice(i))
453c  meme chose pour la glace.
454         if (ptconv(i,k)) then
455            zfroi(i) = dtime/FLOAT(ninter)/zdz(i)*zoliq(i)
456     .              *fallvc(zrhol(i)) * zfice(i)
457         else
458            zfroi(i) = dtime/FLOAT(ninter)/zdz(i)*zoliq(i)
459     .              *fallvs(zrhol(i)) * zfice(i)
460         endif
461         ztot(i) = zchau(i) + zfroi(i)
462         IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
463         ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
464         zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
465         radliq(i,k) = radliq(i,k) + zoliq(i)/FLOAT(ninter+1)
466      ENDIF
467      ENDDO
468      ENDDO
469c
470      DO i = 1, klon
471      IF (rneb(i,k).GT.0.0) THEN
472         d_ql(i,k) = zoliq(i)
473         zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0)
474     .                    * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
475      ENDIF
476      IF (zt(i).LT.RTT) THEN
477        psfl(i,k)=zrfl(i)
478      ELSE
479        prfl(i,k)=zrfl(i)
480      ENDIF
481      ENDDO
482c
483c Calculer les tendances de q et de t:
484c
485      DO i = 1, klon
486         d_q(i,k) = zq(i) - q(i,k)
487         d_t(i,k) = zt(i) - t(i,k)
488      ENDDO
489c
490cAA--------------- Calcul du lessivage stratiforme  -------------
491
492      DO i = 1,klon
493c
494         zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)
495     .                * (paprs(i,k)-paprs(i,k+1))/RG
496         IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
497cAA lessivage nucleation LMD5 dans la couche elle-meme
498            if (t(i,k) .GE. ztglace) THEN
499               zalpha_tr = a_tr_sca(3)
500            else
501               zalpha_tr = a_tr_sca(4)
502            endif
503            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
504            pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
505            frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
506c
507c nucleation avec un facteur -1 au lieu de -0.5
508            zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
509            pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
510         ENDIF
511c
512      ENDDO      ! boucle sur i
513c
514cAA Lessivage par impaction dans les couches en-dessous
515      DO kk = k-1, 1, -1
516        DO i = 1, klon
517          IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
518            if (t(i,kk) .GE. ztglace) THEN
519              zalpha_tr = a_tr_sca(1)
520            else
521              zalpha_tr = a_tr_sca(2)
522            endif
523            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
524            pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
525            frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
526          ENDIF
527        ENDDO
528      ENDDO
529c
530cAA----------------------------------------------------------
531c                     FIN DE BOUCLE SUR K   
532 9999 CONTINUE
533c
534cAA-----------------------------------------------------------
535c
536c Pluie ou neige au sol selon la temperature de la 1ere couche
537c
538      DO i = 1, klon
539      IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
540         snow(i) = zrfl(i)
541         zlh_solid(i) = RLSTT-RLVTT
542      ELSE
543         rain(i) = zrfl(i)
544         zlh_solid(i) = 0.
545      ENDIF
546      ENDDO
547C
548C For energy conservation : when snow is present, the solification
549c latent heat is considered.
550      DO k = 1, klev
551        DO i = 1, klon
552          zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
553          zmair=(paprs(i,k)-paprs(i,k+1))/RG
554          zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
555          d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair)
556        END DO
557      END DO
558c
[883]559
560      if (ncoreczq>0) then
561         print*,'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
562      endif
[524]563      RETURN
564      END
Note: See TracBrowser for help on using the repository browser.