source: LMDZ4/branches/LMDZ4_par_0/libf/phylmd/fisrtilp.F @ 1420

Last change on this file since 1420 was 634, checked in by Laurent Fairhead, 20 years ago

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