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

Last change on this file since 1080 was 883, checked in by Laurent Fairhead, 17 years ago

modifications pour faire de l'aquaplanète FH
LF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.3 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
233      DO i = 1, klon
234cIM
235       IF(k.LE.klevm1) THEN         
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)
242CC        WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
243       ENDIF
244      ENDDO
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
374         do i=1,klon
375            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
376            IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
377            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
378c           zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
379c  On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
380c  la convection.
381c  ATTENTION !!! Il va falloir verifier tout ca.
382            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
383c           print*,'ZDQS ',zdqs(i)
384c--Olivier
385            rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
386            IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
387            IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
388c--fin
389           ENDDO
390      ELSE
391         DO i = 1, klon
392            IF (zq(i).GT.zqs(i)) THEN
393               rneb(i,k) = 1.0
394            ELSE
395               rneb(i,k) = 0.0
396            ENDIF
397            zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
398         ENDDO
399      ENDIF
400c
401      DO i = 1, klon
402         zq(i) = zq(i) - zcond(i)
403c         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
404         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
405      ENDDO
406c
407c Partager l'eau condensee en precipitation et eau liquide nuageuse
408c
409      DO i = 1, klon
410      IF (rneb(i,k).GT.0.0) THEN
411         zoliq(i) = zcond(i)
412         zrho(i) = pplay(i,k) / zt(i) / RD
413         zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
414         zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
415         zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
416         zfice(i) = zfice(i)**nexpo
417         zneb(i) = MAX(rneb(i,k), seuil_neb)
418         radliq(i,k) = zoliq(i)/FLOAT(ninter+1)
419      ENDIF
420      ENDDO
421c
422      DO n = 1, ninter
423      DO i = 1, klon
424      IF (rneb(i,k).GT.0.0) THEN
425         zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
426
427         if (ptconv(i,k)) then
428            zcl(i)=cld_lc_con
429            zct(i)=1./cld_tau_con
430         else
431            zcl(i)=cld_lc_lsc
432            zct(i)=1./cld_tau_lsc
433         endif
[766]434c  quantit�d'eau ��minier.
[524]435         zchau(i) = zct(i)*dtime/FLOAT(ninter) * zoliq(i)
436     .         *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl(i))**2)) *(1.-zfice(i))
437c  meme chose pour la glace.
438         if (ptconv(i,k)) then
439            zfroi(i) = dtime/FLOAT(ninter)/zdz(i)*zoliq(i)
440     .              *fallvc(zrhol(i)) * zfice(i)
441         else
442            zfroi(i) = dtime/FLOAT(ninter)/zdz(i)*zoliq(i)
443     .              *fallvs(zrhol(i)) * zfice(i)
444         endif
445         ztot(i) = zchau(i) + zfroi(i)
446         IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
447         ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
448         zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
449         radliq(i,k) = radliq(i,k) + zoliq(i)/FLOAT(ninter+1)
450      ENDIF
451      ENDDO
452      ENDDO
453c
454      DO i = 1, klon
455      IF (rneb(i,k).GT.0.0) THEN
456         d_ql(i,k) = zoliq(i)
457         zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0)
458     .                    * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
459      ENDIF
460      IF (zt(i).LT.RTT) THEN
461        psfl(i,k)=zrfl(i)
462      ELSE
463        prfl(i,k)=zrfl(i)
464      ENDIF
465      ENDDO
466c
467c Calculer les tendances de q et de t:
468c
469      DO i = 1, klon
470         d_q(i,k) = zq(i) - q(i,k)
471         d_t(i,k) = zt(i) - t(i,k)
472      ENDDO
473c
474cAA--------------- Calcul du lessivage stratiforme  -------------
475
476      DO i = 1,klon
477c
478         zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)
479     .                * (paprs(i,k)-paprs(i,k+1))/RG
480         IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
481cAA lessivage nucleation LMD5 dans la couche elle-meme
482            if (t(i,k) .GE. ztglace) THEN
483               zalpha_tr = a_tr_sca(3)
484            else
485               zalpha_tr = a_tr_sca(4)
486            endif
487            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
488            pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
489            frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
490c
491c nucleation avec un facteur -1 au lieu de -0.5
492            zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
493            pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
494         ENDIF
495c
496      ENDDO      ! boucle sur i
497c
498cAA Lessivage par impaction dans les couches en-dessous
499      DO kk = k-1, 1, -1
500        DO i = 1, klon
501          IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
502            if (t(i,kk) .GE. ztglace) THEN
503              zalpha_tr = a_tr_sca(1)
504            else
505              zalpha_tr = a_tr_sca(2)
506            endif
507            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
508            pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
509            frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
510          ENDIF
511        ENDDO
512      ENDDO
513c
514cAA----------------------------------------------------------
515c                     FIN DE BOUCLE SUR K   
516 9999 CONTINUE
517c
518cAA-----------------------------------------------------------
519c
520c Pluie ou neige au sol selon la temperature de la 1ere couche
521c
522      DO i = 1, klon
523      IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
524         snow(i) = zrfl(i)
525         zlh_solid(i) = RLSTT-RLVTT
526      ELSE
527         rain(i) = zrfl(i)
528         zlh_solid(i) = 0.
529      ENDIF
530      ENDDO
531C
532C For energy conservation : when snow is present, the solification
533c latent heat is considered.
534      DO k = 1, klev
535        DO i = 1, klon
536          zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
537          zmair=(paprs(i,k)-paprs(i,k+1))/RG
538          zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
539          d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair)
540        END DO
541      END DO
542c
[883]543
544      if (ncoreczq>0) then
545         print*,'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
546      endif
[524]547      RETURN
548      END
Note: See TracBrowser for help on using the repository browser.