source: LMDZ5/branches/LF-private/libf/phylmd/fisrtilp.F90 @ 5456

Last change on this file since 5456 was 1855, checked in by Ehouarn Millour, 11 years ago

Bug fix, zpluie and zice are sometimes used without being set.
EM

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.5 KB
Line 
1!
2! $Id: fisrtilp.F90 1855 2013-08-30 13:16:02Z aborella $
3!
4!
5SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
6     d_t, d_q, d_ql, rneb, radliq, rain, snow,          &
7     pfrac_impa, pfrac_nucl, pfrac_1nucl,               &
8     frac_impa, frac_nucl, beta,                        &
9     prfl, psfl, rhcl, zqta, fraca,                     &
10     ztv, zpspsk, ztla, zthl, iflag_cldcon,             &
11     iflag_ice_thermo)
12
13  !
14  USE dimphy
15  IMPLICIT none
16  !======================================================================
17  ! Auteur(s): Z.X. Li (LMD/CNRS)
18  ! Date: le 20 mars 1995
19  ! Objet: condensation et precipitation stratiforme.
20  !        schema de nuage
21  !======================================================================
22  !======================================================================
23  !ym include "dimensions.h"
24  !ym include "dimphy.h"
25  include "YOMCST.h"
26  include "tracstoke.h"
27  include "fisrtilp.h"
28  include "iniprint.h"
29
30  !
31  ! Arguments:
32  !
33  REAL dtime ! intervalle du temps (s)
34  REAL paprs(klon,klev+1) ! pression a inter-couche
35  REAL pplay(klon,klev) ! pression au milieu de couche
36  REAL t(klon,klev) ! temperature (K)
37  REAL q(klon,klev) ! humidite specifique (kg/kg)
38  REAL d_t(klon,klev) ! incrementation de la temperature (K)
39  REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
40  REAL d_ql(klon,klev) ! incrementation de l'eau liquide
41  REAL rneb(klon,klev) ! fraction nuageuse
42  REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
43  REAL rhcl(klon,klev) ! humidite relative en ciel clair
44  REAL rain(klon) ! pluies (mm/s)
45  REAL snow(klon) ! neige (mm/s)
46  REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
47  REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
48  REAL ztv(klon,klev)
49  REAL zqta(klon,klev),fraca(klon,klev)
50  REAL sigma1(klon,klev),sigma2(klon,klev)
51  REAL qltot(klon,klev),ctot(klon,klev)
52  REAL zpspsk(klon,klev),ztla(klon,klev)
53  REAL zthl(klon,klev)
54  REAL ztfondue, qsl, qsi
55
56  logical lognormale(klon)
57  logical ice_thermo
58
59  !AA
60  ! Coeffients de fraction lessivee : pour OFF-LINE
61  !
62  REAL pfrac_nucl(klon,klev)
63  REAL pfrac_1nucl(klon,klev)
64  REAL pfrac_impa(klon,klev)
65  !
66  ! Fraction d'aerosols lessivee par impaction et par nucleation
67  ! POur ON-LINE
68  !
69  REAL frac_impa(klon,klev)
70  REAL frac_nucl(klon,klev)
71  real zct      ,zcl
72  !AA
73  !
74  ! Options du programme:
75  !
76  REAL seuil_neb ! un nuage existe vraiment au-dela
77  PARAMETER (seuil_neb=0.001)
78
79  INTEGER ninter ! sous-intervals pour la precipitation
80  INTEGER ncoreczq 
81  INTEGER iflag_cldcon
82  INTEGER iflag_ice_thermo
83  PARAMETER (ninter=5)
84  LOGICAL evap_prec ! evaporation de la pluie
85  PARAMETER (evap_prec=.TRUE.)
86  REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
87  logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur
88
89  real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
90  real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
91  real erf   
92  REAL qcloud(klon)
93  !
94  LOGICAL cpartiel ! condensation partielle
95  PARAMETER (cpartiel=.TRUE.)
96  REAL t_coup
97  PARAMETER (t_coup=234.0)
98  !
99  ! Variables locales:
100  !
101  INTEGER i, k, n, kk
102  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5   
103  REAL zrfl(klon), zrfln(klon), zqev, zqevt
104  REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
105  REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
106  REAL zoliqp(klon), zoliqi(klon)
107  REAL ztglace, zt(klon)
108  INTEGER nexpo ! exponentiel pour glace/eau
109  REAL zdz(klon),zrho(klon),ztot      , zrhol(klon)
110  REAL zchau      ,zfroi      ,zfice(klon),zneb(klon)
111  REAL zmelt, zpluie, zice, zcondold
112  PARAMETER (ztfondue=278.15)
113  !
114  LOGICAL appel1er
115  SAVE appel1er
116  !$OMP THREADPRIVATE(appel1er)
117  !
118  !---------------------------------------------------------------
119  !
120  !AA Variables traceurs:
121  !AA  Provisoire !!! Parametres alpha du lessivage
122  !AA  A priori on a 4 scavenging # possibles
123  !
124  REAL a_tr_sca(4)
125  save a_tr_sca
126  !$OMP THREADPRIVATE(a_tr_sca)
127  !
128  ! Variables intermediaires
129  !
130  REAL zalpha_tr
131  REAL zfrac_lessi
132  REAL zprec_cond(klon)
133  !AA
134! RomP >>> 15 nov 2012
135  REAL   beta(klon,klev) ! taux de conversion de l'eau cond
136! RomP <<<
137  REAL zmair, zcpair, zcpeau
138  !     Pour la conversion eau-neige
139  REAL zlh_solid(klon), zm_solid
140  !IM
141  !ym      INTEGER klevm1
142  !---------------------------------------------------------------
143  !
144  ! Fonctions en ligne:
145  !
146  REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace
147  REAL zzz
148  include "YOETHF.h"
149  include "FCTTRE.h"
150  fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
151  fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc
152  !
153  DATA appel1er /.TRUE./
154  !ym
155  ice_thermo = iflag_ice_thermo .GE. 1
156  zdelq=0.0
157
158  if (prt_level>9)write(lunout,*)'NUAGES4 A. JAM'
159  IF (appel1er) THEN
160     !
161     WRITE(lunout,*) 'fisrtilp, ninter:', ninter
162     WRITE(lunout,*) 'fisrtilp, evap_prec:', evap_prec
163     WRITE(lunout,*) 'fisrtilp, cpartiel:', cpartiel
164     IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
165        WRITE(lunout,*) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
166        WRITE(lunout,*) 'Je prefere un sous-intervalle de 6 minutes'
167        !         CALL abort
168     ENDIF
169     appel1er = .FALSE.
170     !
171     !AA initialiation provisoire
172     a_tr_sca(1) = -0.5
173     a_tr_sca(2) = -0.5
174     a_tr_sca(3) = -0.5
175     a_tr_sca(4) = -0.5
176     !
177     !AA Initialisation a 1 des coefs des fractions lessivees
178     !
179     !cdir collapse
180     DO k = 1, klev
181        DO i = 1, klon
182           pfrac_nucl(i,k)=1.
183           pfrac_1nucl(i,k)=1.
184           pfrac_impa(i,k)=1.
185           beta(i,k)=0.  !RomP initialisation
186        ENDDO
187     ENDDO
188
189  ENDIF          !  test sur appel1er
190  !
191  !MAf Initialisation a 0 de zoliq
192  !      DO i = 1, klon
193  !         zoliq(i)=0.
194  !      ENDDO
195  ! Determiner les nuages froids par leur temperature
196  !  nexpo regle la raideur de la transition eau liquide / eau glace.
197  !
198  ztglace = RTT - 15.0
199!AJ<
200  IF (ice_thermo) THEN
201    nexpo = 2
202  ELSE
203    nexpo = 6
204  ENDIF
205!!  RLVTT = 2.501e6 ! pas de redefinition des constantes physiques (jyg)
206!!  RLSTT = 2.834e6 ! pas de redefinition des constantes physiques (jyg)
207!>AJ
208  !cc      nexpo = 1
209  !
210  ! Initialiser les sorties:
211  !
212  !cdir collapse
213  DO k = 1, klev+1
214     DO i = 1, klon
215        prfl(i,k) = 0.0
216        psfl(i,k) = 0.0
217     ENDDO
218  ENDDO
219
220  !cdir collapse
221  DO k = 1, klev
222     DO i = 1, klon
223        d_t(i,k) = 0.0
224        d_q(i,k) = 0.0
225        d_ql(i,k) = 0.0
226        rneb(i,k) = 0.0
227        radliq(i,k) = 0.0
228        frac_nucl(i,k) = 1.
229        frac_impa(i,k) = 1.
230     ENDDO
231  ENDDO
232  DO i = 1, klon
233     rain(i) = 0.0
234     snow(i) = 0.0
235     zoliq(i)=0.
236     !     ENDDO
237     !
238     ! Initialiser le flux de precipitation a zero
239     !
240     !     DO i = 1, klon
241     zrfl(i) = 0.0
242     zifl(i) = 0.0
243     zneb(i) = seuil_neb
244  ENDDO
245  !
246  !
247  !AA Pour plus de securite
248
249  zalpha_tr   = 0.
250  zfrac_lessi = 0.
251
252  !AA----------------------------------------------------------
253  !
254  ncoreczq=0
255  ! Boucle verticale (du haut vers le bas)
256  !
257  !IM : klevm1
258  !ym      klevm1=klev-1
259  DO k = klev, 1, -1
260     !
261     !AA----------------------------------------------------------
262     !
263     DO i = 1, klon
264        zt(i)=t(i,k)
265        zq(i)=q(i,k)
266     ENDDO
267     !
268     ! Calculer la varition de temp. de l'air du a la chaleur sensible
269     ! transporter par la pluie.
270     ! Il resterait a rajouter cet effet de la chaleur sensible sur les
271     ! flux de surface, du a la diff. de temp. entre le 1er niveau et la
272     ! surface.
273     !
274     IF(k.LE.klevm1) THEN         
275        DO i = 1, klon
276           !IM
277           zmair=(paprs(i,k)-paprs(i,k+1))/RG
278           zcpair=RCPD*(1.0+RVTMP2*zq(i))
279           zcpeau=RCPD*RVTMP2
280           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau &
281                + zmair*zcpair*zt(i) ) &
282                / (zmair*zcpair + zrfl(i)*dtime*zcpeau)
283           !     C        WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
284        ENDDO
285     ENDIF
286     !
287     !
288     ! Calculer l'evaporation de la precipitation
289     !
290
291
292     ! Calculer l'evaporation de la precipitation
293     !
294
295
296     IF (evap_prec) THEN
297        DO i = 1, klon
298!AJ<
299!!           IF (zrfl(i) .GT.0.) THEN
300           IF (zrfl(i)+zifl(i).GT.0.) THEN
301!>AJ
302              IF (thermcep) THEN
303                 zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
304                 zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
305                 zqs(i)=MIN(0.5,zqs(i))
306                 zcor=1./(1.-RETV*zqs(i))
307                 zqs(i)=zqs(i)*zcor
308              ELSE
309                 IF (zt(i) .LT. t_coup) THEN
310                    zqs(i) = qsats(zt(i)) / pplay(i,k)
311                 ELSE
312                    zqs(i) = qsatl(zt(i)) / pplay(i,k)
313                 ENDIF
314              ENDIF
315           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
316        ENDDO
317!AJ<
318       IF (.NOT. ice_thermo) THEN
319        DO i = 1, klon
320!AJ<
321!!           IF (zrfl(i) .GT.0.) THEN
322           IF (zrfl(i)+zifl(i).GT.0.) THEN
323!>AJ
324                zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
325                zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
326                     * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
327                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
328                     * RG*dtime/(paprs(i,k)-paprs(i,k+1))
329                zqev = MIN (zqev, zqevt)
330                zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
331                     /RG/dtime
332       
333                ! pour la glace, on ré-évapore toute la précip dans la
334                ! couche du dessous
335                ! la glace venant de la couche du dessus est simplement
336                ! dans la couche du dessous.
337       
338                IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
339       
340                zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
341                     * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
342                zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
343                     * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
344                     * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
345                zrfl(i) = zrfln(i)
346                zifl(i) = 0.
347           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
348        ENDDO
349!
350       ELSE ! (.NOT. ice_thermo)
351!
352        DO i = 1, klon
353!AJ<
354!!           IF (zrfl(i) .GT.0.) THEN
355           IF (zrfl(i)+zifl(i).GT.0.) THEN
356!>AJ
357     !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
358     ! Modification de l'évaporation avec la glace
359     ! Différentiation entre précipitation liquide et solide
360     ! On suppose que coef_evai=2*coef_eva
361     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
362     
363         zqev0 = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
364     !    zqev0 = MAX (0.0, zqs(i)-zq(i) )
365
366     !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
367     ! On différencie qsat pour l'eau et la glace
368     ! Si zdelta=1. --> glace
369     ! Si zdelta=0. --> eau liquide
370     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
371         
372         qsl= R2ES*FOEEW(zt(i),0.)/pplay(i,k)
373         qsl= MIN(0.5,qsl)
374         zcor= 1./(1.-RETV*qsl)
375         qsl= qsl*zcor
376         
377         zqevt = 1.*coef_eva*(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
378              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
379         zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
380              *RG*dtime/(paprs(i,k)-paprs(i,k+1))
381
382         qsi= R2ES*FOEEW(zt(i),1.)/pplay(i,k)
383         qsi= MIN(0.5,qsi)
384         zcor= 1./(1.-RETV*qsi)
385         qsi= qsi*zcor
386
387         zqevti = 1.*coef_eva*(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
388              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
389         zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
390              *RG*dtime/(paprs(i,k)-paprs(i,k+1))   
391
392             
393     !JAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
394     ! Vérification sur l'évaporation
395     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
396     
397         IF (zqevt+zqevti.GT.zqev0) THEN
398                  zqev=zqev0*zqevt/(zqevt+zqevti)
399                  zqevi=zqev0*zqevti/(zqevt+zqevti)
400             
401         ELSE
402             IF (zqevt+zqevti.GT.0.) THEN
403                  zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
404                  zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
405             ELSE
406             zqev=0.
407             zqevi=0.
408             ENDIF
409         ENDIF
410     
411         zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
412                                 /RG/dtime)
413         zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
414                                 /RG/dtime)
415         
416     ! Pour la glace, on révapore toute la précip dans la couche du dessous
417     ! la glace venant de la couche du dessus est simplement dans la couche
418     ! du dessous.
419     
420     !    IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
421!         print*,zrfl(i),zrfln(i),zqevt,zqevti,RLMLT,'fluxdeprecip'
422         zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
423                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
424         zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
425                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
426                  * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
427                  + (zifln(i)-zifl(i)) &
428                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
429                  * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
430     
431         zrfl(i) = zrfln(i)
432         zifl(i) = zifln(i)
433
434           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
435        ENDDO
436
437      ENDIF ! (.NOT. ice_thermo)
438     
439     ENDIF ! (evap_prec)
440     !
441     ! Calculer Qs et L/Cp*dQs/dT:
442     !
443     IF (thermcep) THEN
444        DO i = 1, klon
445           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
446           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
447           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
448           zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
449           zqs(i) = MIN(0.5,zqs(i))
450           zcor = 1./(1.-RETV*zqs(i))
451           zqs(i) = zqs(i)*zcor
452           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
453        ENDDO
454     ELSE
455        DO i = 1, klon
456           IF (zt(i).LT.t_coup) THEN
457              zqs(i) = qsats(zt(i))/pplay(i,k)
458              zdqs(i) = dqsats(zt(i),zqs(i))
459           ELSE
460              zqs(i) = qsatl(zt(i))/pplay(i,k)
461              zdqs(i) = dqsatl(zt(i),zqs(i))
462           ENDIF
463        ENDDO
464     ENDIF
465     !
466     ! Determiner la condensation partielle et calculer la quantite
467     ! de l'eau condensee:
468     !
469
470     IF (cpartiel) THEN
471
472        !        print*,'Dans partiel k=',k
473        !
474        !   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
475        !   nuageuse a partir des PDF de Sandrine Bony.
476        !   rneb  : fraction nuageuse
477        !   zqn   : eau totale dans le nuage
478        !   zcond : eau condensee moyenne dans la maille.
479        !  on prend en compte le réchauffement qui diminue la partie
480        ! condensee
481        !
482        !   Version avec les raqts
483
484        if (iflag_pdf.eq.0) then
485
486           do i=1,klon
487              zdelq = min(ratqs(i,k),0.99) * zq(i)
488              rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
489              zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
490           enddo
491
492        else
493           !
494           !   Version avec les nouvelles PDFs.
495           do i=1,klon
496              if(zq(i).lt.1.e-15) then
497                 ncoreczq=ncoreczq+1
498                 zq(i)=1.e-15
499              endif
500           enddo
501
502           if (iflag_cldcon>=5) then
503
504              call cloudth(klon,klev,k,ztv, &
505                   zq,zqta,fraca, &
506                   qcloud,ctot,zpspsk,paprs,ztla,zthl, &
507                   ratqs,zqs,t)
508
509              do i=1,klon
510                 rneb(i,k)=ctot(i,k)
511                 zqn(i)=qcloud(i)
512              enddo
513
514           endif
515
516           if (iflag_cldcon <= 4) then
517              lognormale = .true.
518           elseif (iflag_cldcon >= 6) then
519              ! lognormale en l'absence des thermiques
520              lognormale = fraca(:,k) < 1e-10
521           else
522              ! Dans le cas iflag_cldcon=5, on prend systématiquement la
523              ! bi-gaussienne
524              lognormale = .false.
525           end if
526
527           do i=1,klon
528              if (lognormale(i)) then
529                 zpdf_sig(i)=ratqs(i,k)*zq(i)
530                 zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
531                 zpdf_delta(i)=log(zq(i)/zqs(i))
532                 zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
533                 zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
534                 zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
535                 zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
536                 zpdf_e1(i)=1.-erf(zpdf_e1(i))
537                 zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
538                 zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
539                 zpdf_e2(i)=1.-erf(zpdf_e2(i))
540              endif
541           enddo
542
543           do i=1,klon
544              if (lognormale(i)) then
545                 if (zpdf_e1(i).lt.1.e-10) then
546                    rneb(i,k)=0.
547                    zqn(i)=zqs(i)
548                 else
549                    rneb(i,k)=0.5*zpdf_e1(i)
550                    zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
551                 endif
552              endif
553
554           enddo
555
556
557        endif ! iflag_pdf
558
559        DO i=1,klon
560           IF (rneb(i,k) .LE. 0.0) THEN
561              zqn(i) = 0.0
562              rneb(i,k) = 0.0
563              zcond(i) = 0.0
564              rhcl(i,k)=zq(i)/zqs(i)
565           ELSE IF (rneb(i,k) .GE. 1.0) THEN
566              zqn(i) = zq(i)
567              rneb(i,k) = 1.0                 
568              zcond(i) = MAX(0.0,zqn(i)-zqs(i))/(1+zdqs(i))
569              rhcl(i,k)=1.0
570           ELSE
571              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1+zdqs(i))
572              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
573           ENDIF
574        ENDDO
575        !         do i=1,klon
576        !            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
577        !            IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
578        !            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
579        !c           zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
580        !c  On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
581        !c  la convection.
582        !c  ATTENTION !!! Il va falloir verifier tout ca.
583        !            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
584        !c           print*,'ZDQS ',zdqs(i)
585        !c--Olivier
586        !            rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
587        !            IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
588        !            IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
589        !c--fin
590        !           ENDDO
591     ELSE
592        DO i = 1, klon
593           IF (zq(i).GT.zqs(i)) THEN
594              rneb(i,k) = 1.0
595           ELSE
596              rneb(i,k) = 0.0
597           ENDIF
598           zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
599        ENDDO
600     ENDIF
601     !
602     DO i = 1, klon
603        zq(i) = zq(i) - zcond(i)
604        !         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
605     ENDDO
606!AJ<
607     IF (.NOT. ice_thermo) THEN
608       DO i = 1, klon
609         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
610       ENDDO
611     ELSE
612       DO i = 1, klon
613         zfice(i) = 1.0 - (zt(i)-ztglace) / (273.15-ztglace)
614         zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
615         zfice(i) = zfice(i)**nexpo
616         zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
617                       +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
618!         print*,zt(i),zrfl(i),zifl(i),'temp1'
619       ENDDO
620     ENDIF
621!>AJ
622     !
623     ! Partager l'eau condensee en precipitation et eau liquide nuageuse
624     !
625     DO i = 1, klon
626        IF (rneb(i,k).GT.0.0) THEN
627           zoliq(i) = zcond(i)
628           zrho(i) = pplay(i,k) / zt(i) / RD
629           zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
630        ENDIF
631     ENDDO
632!AJ<
633     IF (.NOT. ice_thermo) THEN
634     DO i = 1, klon
635        IF (rneb(i,k).GT.0.0) THEN
636           zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
637           zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
638           zfice(i) = zfice(i)**nexpo
639     !!      zfice(i)=0.
640        ENDIF
641     ENDDO
642     ENDIF
643     DO i = 1, klon
644        IF (rneb(i,k).GT.0.0) THEN
645           zneb(i) = MAX(rneb(i,k), seuil_neb)
646     !      zt(i) = zt(i)+zcond(i)*zfice(i)*RLMLT/RCPD/(1.0+RVTMP2*zq(i)) 
647!           print*,zt(i),'fractionglace'
648!>AJ
649           radliq(i,k) = zoliq(i)/REAL(ninter+1)
650        ENDIF
651     ENDDO
652     !
653     DO n = 1, ninter
654        DO i = 1, klon
655           IF (rneb(i,k).GT.0.0) THEN
656              zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
657              ! Initialization of zpluie and zice:
658              zpluie=0
659              zice=0
660              IF (zneb(i).EQ.seuil_neb) THEN
661                 ztot = 0.0
662              ELSE
663                 !  quantite d'eau a eliminer: zchau
664                 !  meme chose pour la glace: zfroi
665                 if (ptconv(i,k)) then
666                    zcl   =cld_lc_con
667                    zct   =1./cld_tau_con
668                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
669                         *fallvc(zrhol(i)) * zfice(i)
670                 else
671                    zcl   =cld_lc_lsc
672                    zct   =1./cld_tau_lsc
673                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
674                         *fallvs(zrhol(i)) * zfice(i)
675                 endif
676                 zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
677                      *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl   )**2)) *(1.-zfice(i))
678!AJ<
679                 IF (.NOT. ice_thermo) THEN
680                   ztot    = zchau + zfroi
681                 ELSE
682                   zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
683                   zice = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i))
684                   ztot    = zpluie    + zice
685                 ENDIF
686!>AJ
687                 ztot    = MAX(ztot   ,0.0)
688              ENDIF
689              ztot    = MIN(ztot,zoliq(i))
690!AJ<
691     !         zoliqp = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie   , 0.0)
692     !         zoliqi = MAX(zoliq(i)*zfice(i)-1.*zice   , 0.0)
693              zoliqp(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie  , 0.0)
694              zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zice  , 0.0)
695              zoliq(i) = MAX(zoliq(i)-ztot   , 0.0)
696!>AJ
697              radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
698           ENDIF
699        ENDDO
700     ENDDO
701     !
702         IF (.NOT. ice_thermo) THEN
703       DO i = 1, klon
704         IF (rneb(i,k).GT.0.0) THEN
705           d_ql(i,k) = zoliq(i)
706           zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
707                * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
708         ENDIF
709       ENDDO
710     ELSE
711       DO i = 1, klon
712         IF (rneb(i,k).GT.0.0) THEN
713           d_ql(i,k) = zoliq(i)
714!AJ<
715           zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
716               *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
717           zifl(i) = zifl(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
718                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
719     !      zrfl(i) = zrfl(i)+  zpluie                         &
720     !          *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
721     !      zifl(i) = zifl(i)+  zice                    &
722     !               *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)                                   
723
724         ENDIF                     
725       ENDDO
726     ENDIF
727
728     IF (ice_thermo) THEN
729       DO i = 1, klon
730           zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
731           zmelt = MIN(MAX(zmelt,0.),1.)
732           zrfl(i)=zrfl(i)+zmelt*zifl(i)
733           zifl(i)=zifl(i)*(1.-zmelt)
734!           print*,zt(i),'octavio1'
735           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
736                      *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
737!           print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
738       ENDDO
739     ENDIF
740
741       
742     IF (.NOT. ice_thermo) THEN
743       DO i = 1, klon
744         IF (zt(i).LT.RTT) THEN
745           psfl(i,k)=zrfl(i)
746         ELSE
747           prfl(i,k)=zrfl(i)
748         ENDIF
749       ENDDO
750     ELSE
751     ! JAM*************************************************
752     ! Revoir partie ci-dessous: à quoi servent psfl et prfl?
753     ! *****************************************************
754
755       DO i = 1, klon
756     !   IF (zt(i).LT.RTT) THEN
757           psfl(i,k)=zifl(i)
758     !   ELSE
759           prfl(i,k)=zrfl(i)
760     !   ENDIF
761!>AJ
762       ENDDO
763     ENDIF
764     !
765     !
766     ! Calculer les tendances de q et de t:
767     !
768     DO i = 1, klon
769        d_q(i,k) = zq(i) - q(i,k)
770        d_t(i,k) = zt(i) - t(i,k)
771     ENDDO
772     !
773     !AA--------------- Calcul du lessivage stratiforme  -------------
774
775     DO i = 1,klon
776        !
777        if(zcond(i).gt.zoliq(i)+1.e-10) then
778         beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
779        else
780         beta(i,k) = 0.
781        endif
782        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
783             * (paprs(i,k)-paprs(i,k+1))/RG
784        IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
785           !AA lessivage nucleation LMD5 dans la couche elle-meme
786           if (t(i,k) .GE. ztglace) THEN
787              zalpha_tr = a_tr_sca(3)
788           else
789              zalpha_tr = a_tr_sca(4)
790           endif
791           zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
792           pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
793           frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
794           !
795           ! nucleation avec un facteur -1 au lieu de -0.5
796           zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
797           pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
798        ENDIF
799        !
800     ENDDO      ! boucle sur i
801     !
802     !AA Lessivage par impaction dans les couches en-dessous
803     DO kk = k-1, 1, -1
804        DO i = 1, klon
805           IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
806              if (t(i,kk) .GE. ztglace) THEN
807                 zalpha_tr = a_tr_sca(1)
808              else
809                 zalpha_tr = a_tr_sca(2)
810              endif
811              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
812              pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
813              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
814           ENDIF
815        ENDDO
816     ENDDO
817     !
818     !AA----------------------------------------------------------
819     !                     FIN DE BOUCLE SUR K   
820  end DO
821  !
822  !AA-----------------------------------------------------------
823  !
824  ! Pluie ou neige au sol selon la temperature de la 1ere couche
825  !
826  DO i = 1, klon
827     IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
828!AJ<
829!!        snow(i) = zrfl(i)
830        snow(i) = zrfl(i)+zifl(i)
831!>AJ
832        zlh_solid(i) = RLSTT-RLVTT
833     ELSE
834        rain(i) = zrfl(i)
835        zlh_solid(i) = 0.
836     ENDIF
837  ENDDO
838  !
839  ! For energy conservation : when snow is present, the solification
840  ! latent heat is considered.
841  DO k = 1, klev
842     DO i = 1, klon
843        zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
844        zmair=(paprs(i,k)-paprs(i,k+1))/RG
845        zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
846        d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair)
847     END DO
848  END DO
849  !
850
851  if (ncoreczq>0) then
852     WRITE(lunout,*)'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
853  endif
854
855END SUBROUTINE fisrtilp
Note: See TracBrowser for help on using the repository browser.