LMDZ
fisrtilp.F90
Go to the documentation of this file.
1 !
2 ! $Id: fisrtilp.F90 2346 2015-08-21 15:13:46Z emillour $
3 !
4 !
5 SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
6  d_t, d_q, d_ql, d_qi, 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_cld_th, &
11  iflag_ice_thermo)
12 
13  !
14  USE dimphy
15  USE icefrac_lsc_mod ! compute ice fraction (JBM 3/14)
17  IMPLICIT none
18  !======================================================================
19  ! Auteur(s): Z.X. Li (LMD/CNRS)
20  ! Date: le 20 mars 1995
21  ! Objet: condensation et precipitation stratiforme.
22  ! schema de nuage
23  !======================================================================
24  !======================================================================
25  include "YOMCST.h"
26  include "fisrtilp.h"
27  include "nuage.h" ! JBM (3/14)
28 
29  !
30  ! Arguments:
31  !
32  REAL dtime ! intervalle du temps (s)
33  REAL paprs(klon,klev+1) ! pression a inter-couche
34  REAL pplay(klon,klev) ! pression au milieu de couche
35  REAL t(klon,klev) ! temperature (K)
36  REAL q(klon,klev) ! humidite specifique (kg/kg)
37  REAL d_t(klon,klev) ! incrementation de la temperature (K)
38  REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
39  REAL d_ql(klon,klev) ! incrementation de l'eau liquide
40  REAL d_qi(klon,klev) ! incrementation de l'eau glace
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_cld_th
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 Tbef(klon),qlbef(klon),DT(klon),num,denom
104  LOGICAL convergence(klon)
105  REAL DDT0
106  parameter(ddt0=.01)
107  INTEGER n_i(klon), iter
108  REAL cste
109 
110  REAL zrfl(klon), zrfln(klon), zqev, zqevt
111  REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
112  REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
113  REAL zoliqp(klon), zoliqi(klon)
114  REAL zt(klon)
115 ! JBM (3/14) nexpo is replaced by exposant_glace
116 ! REAL nexpo ! exponentiel pour glace/eau
117 ! INTEGER, PARAMETER :: nexpo=6
118  INTEGER exposant_glace_old
119  REAL t_glace_min_old
120  REAL zdz(klon),zrho(klon),ztot , zrhol(klon)
121  REAL zchau ,zfroi ,zfice(klon),zneb(klon)
122  REAL zmelt, zpluie, zice, zcondold
123  parameter(ztfondue=278.15)
124  REAL dzfice(klon)
125  !
126  LOGICAL appel1er
127  SAVE appel1er
128  !$OMP THREADPRIVATE(appel1er)
129  !
130  !---------------------------------------------------------------
131  !
132  !AA Variables traceurs:
133  !AA Provisoire !!! Parametres alpha du lessivage
134  !AA A priori on a 4 scavenging # possibles
135  !
136  REAL a_tr_sca(4)
137  save a_tr_sca
138  !$OMP THREADPRIVATE(a_tr_sca)
139  !
140  ! Variables intermediaires
141  !
142  REAL zalpha_tr
143  REAL zfrac_lessi
144  REAL zprec_cond(klon)
145  !AA
146 ! RomP >>> 15 nov 2012
147  REAL beta(klon,klev) ! taux de conversion de l'eau cond
148 ! RomP <<<
149  REAL zmair, zcpair, zcpeau
150  ! Pour la conversion eau-neige
151  REAL zlh_solid(klon), zm_solid
152  !IM
153  !ym INTEGER klevm1
154  !---------------------------------------------------------------
155  !
156  ! Fonctions en ligne:
157  !
158  REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace
159  REAL zzz
160  include "YOETHF.h"
161  include "FCTTRE.h"
162  fallvc(zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
163  fallvs(zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc
164  !
165  DATA appel1er /.true./
166  !ym
167 !CR: pour iflag_ice_thermo=2, on active que la convection
168 ! ice_thermo = iflag_ice_thermo .GE. 1
169  ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
170  zdelq=0.0
171 
172  if (prt_level>9)write(lunout,*)'NUAGES4 A. JAM'
173  IF (appel1er) THEN
174  !
175  WRITE(lunout,*) 'fisrtilp, ninter:', ninter
176  WRITE(lunout,*) 'fisrtilp, evap_prec:', evap_prec
177  WRITE(lunout,*) 'fisrtilp, cpartiel:', cpartiel
178  IF (abs(dtime/REAL(ninter)-360.0).GT.0.001) then
179  WRITE(lunout,*) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
180  WRITE(lunout,*) 'Je prefere un sous-intervalle de 6 minutes'
181  ! CALL abort
182  ENDIF
183  appel1er = .false.
184  !
185  !AA initialiation provisoire
186  a_tr_sca(1) = -0.5
187  a_tr_sca(2) = -0.5
188  a_tr_sca(3) = -0.5
189  a_tr_sca(4) = -0.5
190  !
191  !AA Initialisation a 1 des coefs des fractions lessivees
192  !
193  !cdir collapse
194  DO k = 1, klev
195  DO i = 1, klon
196  pfrac_nucl(i,k)=1.
197  pfrac_1nucl(i,k)=1.
198  pfrac_impa(i,k)=1.
199  beta(i,k)=0. !RomP initialisation
200  ENDDO
201  ENDDO
202 
203  ENDIF ! test sur appel1er
204  !
205  !MAf Initialisation a 0 de zoliq
206  ! DO i = 1, klon
207  ! zoliq(i)=0.
208  ! ENDDO
209  ! Determiner les nuages froids par leur temperature
210  ! nexpo regle la raideur de la transition eau liquide / eau glace.
211  !
212 !CR: on est oblige de definir des valeurs fisrt car les valeurs de newmicro ne sont pas les memes par defaut
213  IF (iflag_t_glace.EQ.0) THEN
214 ! ztglace = RTT - 15.0
215  t_glace_min_old = rtt - 15.0
216  !AJ<
217  IF (ice_thermo) THEN
218 ! nexpo = 2
219  exposant_glace_old = 2
220  ELSE
221 ! nexpo = 6
222  exposant_glace_old = 6
223  ENDIF
224 
225  ENDIF
226 
227 !! RLVTT = 2.501e6 ! pas de redefinition des constantes physiques (jyg)
228 !! RLSTT = 2.834e6 ! pas de redefinition des constantes physiques (jyg)
230  !cc nexpo = 1
231  !
232  ! Initialiser les sorties:
233  !
234  !cdir collapse
235  DO k = 1, klev+1
236  DO i = 1, klon
237  prfl(i,k) = 0.0
238  psfl(i,k) = 0.0
239  ENDDO
240  ENDDO
241 
242  !cdir collapse
243  DO k = 1, klev
244  DO i = 1, klon
245  d_t(i,k) = 0.0
246  d_q(i,k) = 0.0
247  d_ql(i,k) = 0.0
248  d_qi(i,k) = 0.0
249  rneb(i,k) = 0.0
250  radliq(i,k) = 0.0
251  frac_nucl(i,k) = 1.
252  frac_impa(i,k) = 1.
253  ENDDO
254  ENDDO
255  DO i = 1, klon
256  rain(i) = 0.0
257  snow(i) = 0.0
258  zoliq(i)=0.
259  ! ENDDO
260  !
261  ! Initialiser le flux de precipitation a zero
262  !
263  ! DO i = 1, klon
264  zrfl(i) = 0.0
265  zifl(i) = 0.0
266  zneb(i) = seuil_neb
267  ENDDO
268  !
269  !
270  !AA Pour plus de securite
271 
272  zalpha_tr = 0.
273  zfrac_lessi = 0.
274 
275  !AA----------------------------------------------------------
276  !
277  ncoreczq=0
278  ! Boucle verticale (du haut vers le bas)
279  !
280  !IM : klevm1
281  !ym klevm1=klev-1
282  DO k = klev, 1, -1
283  !
284  !AA----------------------------------------------------------
285  !
286  DO i = 1, klon
287  zt(i)=t(i,k)
288  zq(i)=q(i,k)
289  ENDDO
290  !
291  ! Calculer la varition de temp. de l'air du a la chaleur sensible
292  ! transporter par la pluie.
293  ! Il resterait a rajouter cet effet de la chaleur sensible sur les
294  ! flux de surface, du a la diff. de temp. entre le 1er niveau et la
295  ! surface.
296  !
297  IF(k.LE.klevm1) THEN
298  DO i = 1, klon
299  !IM
300  zmair=(paprs(i,k)-paprs(i,k+1))/rg
301  zcpair=rcpd*(1.0+rvtmp2*zq(i))
302  zcpeau=rcpd*rvtmp2
303  zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau &
304  + zmair*zcpair*zt(i) ) &
305  / (zmair*zcpair + zrfl(i)*dtime*zcpeau)
306  ! C WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
307  ENDDO
308  ENDIF
309  !
310  !
311  ! Calculer l'evaporation de la precipitation
312  !
313 
314 
315  ! Calculer l'evaporation de la precipitation
316  !
317 
318 
319  IF (evap_prec) THEN
320  DO i = 1, klon
321 !AJ<
322 !! IF (zrfl(i) .GT.0.) THEN
323  IF (zrfl(i)+zifl(i).GT.0.) THEN
325  IF (thermcep) THEN
326  zdelta=max(0.,sign(1.,rtt-zt(i)))
327  zqs(i)= r2es*foeew(zt(i),zdelta)/pplay(i,k)
328  zqs(i)=min(0.5,zqs(i))
329  zcor=1./(1.-retv*zqs(i))
330  zqs(i)=zqs(i)*zcor
331  ELSE
332  IF (zt(i) .LT. t_coup) THEN
333  zqs(i) = qsats(zt(i)) / pplay(i,k)
334  ELSE
335  zqs(i) = qsatl(zt(i)) / pplay(i,k)
336  ENDIF
337  ENDIF
338  ENDIF ! (zrfl(i)+zifl(i).GT.0.)
339  ENDDO
340 !AJ<
341  IF (.NOT. ice_thermo) THEN
342  DO i = 1, klon
343 !AJ<
344 !! IF (zrfl(i) .GT.0.) THEN
345  IF (zrfl(i)+zifl(i).GT.0.) THEN
347  zqev = max(0.0, (zqs(i)-zq(i))*zneb(i) )
348  zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * sqrt(zrfl(i)) &
349  * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*rd/rg
350  zqevt = max(0.0,min(zqevt,zrfl(i))) &
351  * rg*dtime/(paprs(i,k)-paprs(i,k+1))
352  zqev = min(zqev, zqevt)
353  zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
354  /rg/dtime
355 
356  ! pour la glace, on ré-évapore toute la précip dans la
357  ! couche du dessous
358  ! la glace venant de la couche du dessus est simplement
359  ! dans la couche du dessous.
360 
361  IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
362 
363  zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
364  * (rg/(paprs(i,k)-paprs(i,k+1)))*dtime
365  zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
366  * (rg/(paprs(i,k)-paprs(i,k+1)))*dtime &
367  * rlvtt/rcpd/(1.0+rvtmp2*zq(i))
368  zrfl(i) = zrfln(i)
369  zifl(i) = 0.
370  ENDIF ! (zrfl(i)+zifl(i).GT.0.)
371  ENDDO
372 !
373  ELSE ! (.NOT. ice_thermo)
374 !
375  DO i = 1, klon
376 !AJ<
377 !! IF (zrfl(i) .GT.0.) THEN
378  IF (zrfl(i)+zifl(i).GT.0.) THEN
380  !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
381  ! Modification de l'évaporation avec la glace
382  ! Différentiation entre précipitation liquide et solide
383  ! On suppose que coef_evai=2*coef_eva
384  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
385 
386  zqev0 = max(0.0, (zqs(i)-zq(i))*zneb(i) )
387  ! zqev0 = MAX (0.0, zqs(i)-zq(i) )
388 
389  !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
390  ! On différencie qsat pour l'eau et la glace
391  ! Si zdelta=1. --> glace
392  ! Si zdelta=0. --> eau liquide
393  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
394 
395  qsl= r2es*foeew(zt(i),0.)/pplay(i,k)
396  qsl= min(0.5,qsl)
397  zcor= 1./(1.-retv*qsl)
398  qsl= qsl*zcor
399 
400  zqevt = 1.*coef_eva*(1.0-zq(i)/qsl)*sqrt(zrfl(i)) &
401  *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*rd/rg
402  zqevt = max(0.0,min(zqevt,zrfl(i))) &
403  *rg*dtime/(paprs(i,k)-paprs(i,k+1))
404 
405  qsi= r2es*foeew(zt(i),1.)/pplay(i,k)
406  qsi= min(0.5,qsi)
407  zcor= 1./(1.-retv*qsi)
408  qsi= qsi*zcor
409 
410  zqevti = 1.*coef_eva*(1.0-zq(i)/qsi)*sqrt(zifl(i)) &
411  *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*rd/rg
412  zqevti = max(0.0,min(zqevti,zifl(i))) &
413  *rg*dtime/(paprs(i,k)-paprs(i,k+1))
414 
415 
416  !JAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
417  ! Vérification sur l'évaporation
418  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
419 
420  IF (zqevt+zqevti.GT.zqev0) THEN
421  zqev=zqev0*zqevt/(zqevt+zqevti)
422  zqevi=zqev0*zqevti/(zqevt+zqevti)
423 
424  ELSE
425  IF (zqevt+zqevti.GT.0.) THEN
426  zqev=min(zqev0*zqevt/(zqevt+zqevti),zqevt)
427  zqevi=min(zqev0*zqevti/(zqevt+zqevti),zqevti)
428  ELSE
429  zqev=0.
430  zqevi=0.
431  ENDIF
432  ENDIF
433 
434  zrfln(i) = max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
435  /rg/dtime)
436  zifln(i) = max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
437  /rg/dtime)
438 
439  ! Pour la glace, on révapore toute la précip dans la couche du dessous
440  ! la glace venant de la couche du dessus est simplement dans la couche
441  ! du dessous.
442 
443  ! IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
444 ! print*,zrfl(i),zrfln(i),zqevt,zqevti,RLMLT,'fluxdeprecip'
445  zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
446  * (rg/(paprs(i,k)-paprs(i,k+1)))*dtime
447  zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
448  * (rg/(paprs(i,k)-paprs(i,k+1)))*dtime &
449  * rlvtt/rcpd/(1.0+rvtmp2*zq(i)) &
450  + (zifln(i)-zifl(i)) &
451  * (rg/(paprs(i,k)-paprs(i,k+1)))*dtime &
452  * rlstt/rcpd/(1.0+rvtmp2*zq(i))
453 
454  zrfl(i) = zrfln(i)
455  zifl(i) = zifln(i)
456 
457 !CR ATTENTION: deplacement de la fonte de la glace
458  zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
459  zmelt = min(max(zmelt,0.),1.)
460  zrfl(i)=zrfl(i)+zmelt*zifl(i)
461  zifl(i)=zifl(i)*(1.-zmelt)
462 ! print*,zt(i),'octavio1'
463  zt(i)=zt(i)-zifl(i)*zmelt*(rg*dtime)/(paprs(i,k)-paprs(i,k+1)) &
464  *rlmlt/rcpd/(1.0+rvtmp2*zq(i))
465 ! print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
466 !fin CR
467 
468 
469 
470  ENDIF ! (zrfl(i)+zifl(i).GT.0.)
471  ENDDO
472 
473  ENDIF ! (.NOT. ice_thermo)
474 
475  ENDIF ! (evap_prec)
476  !
477  ! Calculer Qs et L/Cp*dQs/dT:
478  !
479  IF (thermcep) THEN
480  DO i = 1, klon
481  zdelta = max(0.,sign(1.,rtt-zt(i)))
482  zcvm5 = r5les*rlvtt*(1.-zdelta) + r5ies*rlstt*zdelta
483  zcvm5 = zcvm5 /rcpd/(1.0+rvtmp2*zq(i))
484  zqs(i) = r2es*foeew(zt(i),zdelta)/pplay(i,k)
485  zqs(i) = min(0.5,zqs(i))
486  zcor = 1./(1.-retv*zqs(i))
487  zqs(i) = zqs(i)*zcor
488  zdqs(i) = foede(zt(i),zdelta,zcvm5,zqs(i),zcor)
489  ENDDO
490  ELSE
491  DO i = 1, klon
492  IF (zt(i).LT.t_coup) THEN
493  zqs(i) = qsats(zt(i))/pplay(i,k)
494  zdqs(i) = dqsats(zt(i),zqs(i))
495  ELSE
496  zqs(i) = qsatl(zt(i))/pplay(i,k)
497  zdqs(i) = dqsatl(zt(i),zqs(i))
498  ENDIF
499  ENDDO
500  ENDIF
501  !
502  ! Determiner la condensation partielle et calculer la quantite
503  ! de l'eau condensee:
504  !
505 !verification de la valeur de iflag_fisrtilp_qsat pour iflag_ice_thermo=1
506 ! if ((iflag_ice_thermo.eq.1).and.(iflag_fisrtilp_qsat.ne.0)) then
507 ! write(*,*) " iflag_ice_thermo==1 requires iflag_fisrtilp_qsat==0", &
508 ! " but iflag_fisrtilp_qsat=",iflag_fisrtilp_qsat, ". Might as well stop here."
509 ! stop
510 ! endif
511 
512  IF (cpartiel) THEN
513 
514  ! print*,'Dans partiel k=',k
515  !
516  ! Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
517  ! nuageuse a partir des PDF de Sandrine Bony.
518  ! rneb : fraction nuageuse
519  ! zqn : eau totale dans le nuage
520  ! zcond : eau condensee moyenne dans la maille.
521  ! on prend en compte le réchauffement qui diminue la partie
522  ! condensee
523  !
524  ! Version avec les raqts
525 
526  if (iflag_pdf.eq.0) then
527 
528  do i=1,klon
529  zdelq = min(ratqs(i,k),0.99) * zq(i)
530  rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
531  zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
532  enddo
533 
534  else
535  !
536  ! Version avec les nouvelles PDFs.
537  do i=1,klon
538  if(zq(i).lt.1.e-15) then
539  ncoreczq=ncoreczq+1
540  zq(i)=1.e-15
541  endif
542  enddo
543 
544  if (iflag_cld_th>=5) then
545 
546  call cloudth(klon,klev,k,ztv, &
547  zq,zqta,fraca, &
548  qcloud,ctot,zpspsk,paprs,ztla,zthl, &
549  ratqs,zqs,t)
550 
551  do i=1,klon
552  rneb(i,k)=ctot(i,k)
553  zqn(i)=qcloud(i)
554  enddo
555 
556  endif
557 
558  if (iflag_cld_th <= 4) then
559  lognormale = .true.
560  elseif (iflag_cld_th >= 6) then
561  ! lognormale en l'absence des thermiques
562  lognormale = fraca(:,k) < 1e-10
563  else
564  ! Dans le cas iflag_cld_th=5, on prend systématiquement la
565  ! bi-gaussienne
566  lognormale = .false.
567  end if
568 
569 !CR: variation de qsat avec T en présence de glace ou non
570 !initialisations
571  do i=1,klon
572  dt(i) = 0.
573  n_i(i)=0
574  tbef(i)=zt(i)
575  qlbef(i)=0.
576  enddo
577 
578 
579 !Boucle iterative: ATTENTION, l'option -1 n'est plus activable ici
580  if (iflag_fisrtilp_qsat.ge.0) then
581  do iter=1,iflag_fisrtilp_qsat+1
582 
583  do i=1,klon
584 ! do while ((abs(DT(i)).gt.DDT0).or.(n_i(i).eq.0))
585  convergence(i)=abs(dt(i)).gt.ddt0
586  if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
587  tbef(i)=tbef(i)+dt(i)
588  if (.not.ice_thermo) then
589  zdelta = max(0.,sign(1.,rtt-tbef(i)))
590  else
591  if (iflag_t_glace.eq.0) then
592  zdelta = max(0.,sign(1.,t_glace_min_old-tbef(i)))
593  else if (iflag_t_glace.eq.1) then
594  zdelta = max(0.,sign(1.,t_glace_min-tbef(i)))
595  endif
596  endif
597  zcvm5 = r5les*rlvtt*(1.-zdelta) + r5ies*rlstt*zdelta
598  zcvm5 = zcvm5 /rcpd/(1.0+rvtmp2*zq(i))
599  zqs(i) = r2es*foeew(tbef(i),zdelta)/pplay(i,k)
600  zqs(i) = min(0.5,zqs(i))
601  zcor = 1./(1.-retv*zqs(i))
602  zqs(i) = zqs(i)*zcor
603  zdqs(i) = foede(tbef(i),zdelta,zcvm5,zqs(i),zcor)
604  zpdf_sig(i)=ratqs(i,k)*zq(i)
605  zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
606  zpdf_delta(i)=log(zq(i)/zqs(i))
607  zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
608  zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
609  zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
610  zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
611  zpdf_e1(i)=1.-erf(zpdf_e1(i))
612  zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
613  zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
614  zpdf_e2(i)=1.-erf(zpdf_e2(i))
615 
616  if (zpdf_e1(i).lt.1.e-10) then
617  rneb(i,k)=0.
618  zqn(i)=zqs(i)
619  else
620  rneb(i,k)=0.5*zpdf_e1(i)
621  zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
622  endif
623 
624  endif !convergence
625  enddo ! boucle en i
626 
627  if (.not. ice_thermo) then
628 
629  do i=1,klon
630  if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
631 
632  qlbef(i)=max(0.,zqn(i)-zqs(i))
633  num=-tbef(i)+zt(i)+rneb(i,k)*rlvtt/rcpd/(1.0+rvtmp2*zq(i))*qlbef(i)
634  denom=1.+rneb(i,k)*zdqs(i)
635  dt(i)=num/denom
636  n_i(i)=n_i(i)+1
637  endif
638  enddo
639 
640  else
641 
642 !calcul de la fraction de glace
643 !CR: on utilise la nouvelle fonction de JBM pour l ancien calcul
644 ! zfice(i) = icefrac_lsc(Tbef(i), t_glace_min, &
645 ! t_glace_max, exposant_glace)
646 ! zfice(i) = 1.0 - (Tbef(i)-ztglace) / (RTT-ztglace)
647 ! zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
648 ! zfice(i) = zfice(i)**nexpo
649  if (iflag_t_glace.eq.1) then
650  CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
651  endif
652 
653  do i=1,klon
654  if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
655 
656  if (iflag_t_glace.eq.0) then
657  zfice(i) = 1.0 - (tbef(i)-t_glace_min_old) / (rtt-t_glace_min_old)
658  zfice(i) = min(max(zfice(i),0.0),1.0)
659  zfice(i) = zfice(i)**exposant_glace_old
660  dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) / (t_glace_min_old - rtt)
661  endif
662 
663  if (iflag_t_glace.eq.1) then
664  dzfice(i)= exposant_glace * zfice(i)**(exposant_glace-1) / (t_glace_min - t_glace_max)
665  endif
666 
667  if ((zfice(i).eq.0).or.(zfice(i).eq.1)) then
668  dzfice(i)=0.
669  endif
670 
671  if (zfice(i).lt.1) then
672  cste=rlvtt
673  else
674  cste=rlstt
675  endif
676 
677  qlbef(i)=max(0.,zqn(i)-zqs(i))
678  num=-tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*rlvtt+zfice(i)*rlstt)/rcpd/(1.0+rvtmp2*zq(i))*qlbef(i)
679  denom=1.+rneb(i,k)*((1-zfice(i))*rlvtt+zfice(i)*rlstt)/cste*zdqs(i) &
680  -(rlstt-rlvtt)/rcpd/(1.0+rvtmp2*zq(i))*rneb(i,k)*qlbef(i)*dzfice(i)
681  dt(i)=num/denom
682  n_i(i)=n_i(i)+1
683 
684  endif ! fin convergence
685  enddo ! fin boucle i
686 
687  endif !ice_thermo
688 
689 ! endif
690 ! enddo
691 
692 
693  enddo
694  endif
695 
696 
697  endif ! iflag_pdf
698 
699 
700 ! if (iflag_fisrtilp_qsat.eq.-1) then
701 !CR: ATTENTION option fausse mais a existe: pour la re-activer, prendre iflag_fisrtilp_qsat=0 et activer les lignes suivantes:
702  IF (1.eq.0) THEN
703  DO i=1,klon
704  IF (rneb(i,k) .LE. 0.0) THEN
705  zqn(i) = 0.0
706  rneb(i,k) = 0.0
707  zcond(i) = 0.0
708  rhcl(i,k)=zq(i)/zqs(i)
709  ELSE IF (rneb(i,k) .GE. 1.0) THEN
710  zqn(i) = zq(i)
711  rneb(i,k) = 1.0
712  zcond(i) = max(0.0,zqn(i)-zqs(i))/(1+zdqs(i))
713  rhcl(i,k)=1.0
714  ELSE
715  zcond(i) = max(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1+zdqs(i))
716  rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
717  ENDIF
718  ENDDO
719  ENDIF
720 
721 ! ELSE
722 
723  DO i=1,klon
724  IF (rneb(i,k) .LE. 0.0) THEN
725  zqn(i) = 0.0
726  rneb(i,k) = 0.0
727  zcond(i) = 0.0
728  rhcl(i,k)=zq(i)/zqs(i)
729  ELSE IF (rneb(i,k) .GE. 1.0) THEN
730  zqn(i) = zq(i)
731  rneb(i,k) = 1.0
732  zcond(i) = max(0.0,zqn(i)-zqs(i))
733  rhcl(i,k)=1.0
734  ELSE
735  zcond(i) = max(0.0,zqn(i)-zqs(i))*rneb(i,k)
736  rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
737  ENDIF
738  ENDDO
739 
740 
741 ! ENDIF
742 
743  ! do i=1,klon
744  ! IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
745  ! IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
746  ! rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
747  !c zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
748  !c On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
749  !c la convection.
750  !c ATTENTION !!! Il va falloir verifier tout ca.
751  ! zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
752  !c print*,'ZDQS ',zdqs(i)
753  !c--Olivier
754  ! rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
755  ! IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
756  ! IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
757  !c--fin
758  ! ENDDO
759  ELSE
760  DO i = 1, klon
761  IF (zq(i).GT.zqs(i)) THEN
762  rneb(i,k) = 1.0
763  ELSE
764  rneb(i,k) = 0.0
765  ENDIF
766  zcond(i) = max(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
767  ENDDO
768  ENDIF
769  !
770  DO i = 1, klon
771  zq(i) = zq(i) - zcond(i)
772  ! zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
773  ENDDO
774 !AJ<
775  IF (.NOT. ice_thermo) THEN
776  if (iflag_fisrtilp_qsat.lt.1) then
777  DO i = 1, klon
778  zt(i) = zt(i) + zcond(i) * rlvtt/rcpd/(1.0+rvtmp2*zq(i))
779  ENDDO
780  else if (iflag_fisrtilp_qsat.gt.0) then
781  DO i= 1, klon
782  zt(i) = zt(i) + zcond(i) * rlvtt/rcpd/(1.0+rvtmp2*(zq(i)+zcond(i)))
783  ENDDO
784  endif
785  ELSE
786  if (iflag_t_glace.eq.1) then
787  CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
788  endif
789  if (iflag_fisrtilp_qsat.lt.1) then
790  DO i = 1, klon
791 ! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
792 ! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
793 ! t_glace_max, exposant_glace)
794  if (iflag_t_glace.eq.0) then
795  zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (rtt-t_glace_min_old)
796  zfice(i) = min(max(zfice(i),0.0),1.0)
797  zfice(i) = zfice(i)**exposant_glace_old
798  endif
799  zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * rlvtt/rcpd/(1.0+rvtmp2*zq(i)) &
800  +zfice(i)*zcond(i) * rlstt/rcpd/(1.0+rvtmp2*zq(i))
801  ENDDO
802  else
803  DO i=1, klon
804 ! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
805 ! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
806 ! t_glace_max, exposant_glace)
807  if (iflag_t_glace.eq.0) then
808  zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (rtt-t_glace_min_old)
809  zfice(i) = min(max(zfice(i),0.0),1.0)
810  zfice(i) = zfice(i)**exposant_glace_old
811  endif
812  zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * rlvtt/rcpd/(1.0+rvtmp2*(zq(i)+zcond(i))) &
813  +zfice(i)*zcond(i) * rlstt/rcpd/(1.0+rvtmp2*(zq(i)+zcond(i)))
814  ENDDO
815  endif
816 ! print*,zt(i),zrfl(i),zifl(i),'temp1'
817  ENDIF
819  !
820  ! Partager l'eau condensee en precipitation et eau liquide nuageuse
821  !
822  DO i = 1, klon
823  IF (rneb(i,k).GT.0.0) THEN
824  zoliq(i) = zcond(i)
825  zrho(i) = pplay(i,k) / zt(i) / rd
826  zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*rg)
827  ENDIF
828  ENDDO
829 !AJ<
830  IF (.NOT. ice_thermo) THEN
831  IF (iflag_t_glace.EQ.0) THEN
832  DO i = 1, klon
833  IF (rneb(i,k).GT.0.0) THEN
834  zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.13-t_glace_min_old)
835  zfice(i) = min(max(zfice(i),0.0),1.0)
836  zfice(i) = zfice(i)**exposant_glace_old
837 ! zfice(i) = zfice(i)**nexpo
838  !! zfice(i)=0.
839  ENDIF
840  ENDDO
841  ELSE ! of IF (iflag_t_glace.EQ.0)
842  CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
843 ! DO i = 1, klon
844 ! IF (rneb(i,k).GT.0.0) THEN
845 ! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
846 ! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
847 ! t_glace_max, exposant_glace)
848 ! ENDIF
849 ! ENDDO
850  ENDIF
851  ENDIF
852  DO i = 1, klon
853  IF (rneb(i,k).GT.0.0) THEN
854  zneb(i) = max(rneb(i,k), seuil_neb)
855  ! zt(i) = zt(i)+zcond(i)*zfice(i)*RLMLT/RCPD/(1.0+RVTMP2*zq(i))
856 ! print*,zt(i),'fractionglace'
858  radliq(i,k) = zoliq(i)/REAL(ninter+1)
859  ENDIF
860  ENDDO
861  !
862  DO n = 1, ninter
863  DO i = 1, klon
864  IF (rneb(i,k).GT.0.0) THEN
865  zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
866  ! Initialization of zpluie and zice:
867  zpluie=0
868  zice=0
869  IF (zneb(i).EQ.seuil_neb) THEN
870  ztot = 0.0
871  ELSE
872  ! quantite d'eau a eliminer: zchau
873  ! meme chose pour la glace: zfroi
874  if (ptconv(i,k)) then
875  zcl =cld_lc_con
876  zct =1./cld_tau_con
877  zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
878  *fallvc(zrhol(i)) * zfice(i)
879  else
880  zcl =cld_lc_lsc
881  zct =1./cld_tau_lsc
882  zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
883  *fallvs(zrhol(i)) * zfice(i)
884  endif
885  zchau = zct *dtime/REAL(ninter) * zoliq(i) &
886  *(1.0-exp(-(zoliq(i)/zneb(i)/zcl )**2)) *(1.-zfice(i))
887 !AJ<
888  IF (.NOT. ice_thermo) THEN
889  ztot = zchau + zfroi
890  ELSE
891  zpluie = min(max(zchau,0.0),zoliq(i)*(1.-zfice(i)))
892  zice = min(max(zfroi,0.0),zoliq(i)*zfice(i))
893  ztot = zpluie + zice
894  ENDIF
896  ztot = max(ztot ,0.0)
897  ENDIF
898  ztot = min(ztot,zoliq(i))
899 !AJ<
900  ! zoliqp = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie , 0.0)
901  ! zoliqi = MAX(zoliq(i)*zfice(i)-1.*zice , 0.0)
902  zoliqp(i) = max(zoliq(i)*(1.-zfice(i))-1.*zpluie , 0.0)
903  zoliqi(i) = max(zoliq(i)*zfice(i)-1.*zice , 0.0)
904  zoliq(i) = max(zoliq(i)-ztot , 0.0)
906  radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
907  ENDIF
908  ENDDO
909  ENDDO
910  !
911  IF (.NOT. ice_thermo) THEN
912  DO i = 1, klon
913  IF (rneb(i,k).GT.0.0) THEN
914  d_ql(i,k) = zoliq(i)
915  zrfl(i) = zrfl(i)+ max(zcond(i)-zoliq(i),0.0) &
916  * (paprs(i,k)-paprs(i,k+1))/(rg*dtime)
917  ENDIF
918  ENDDO
919  ELSE
920  DO i = 1, klon
921  IF (rneb(i,k).GT.0.0) THEN
922 !CR on prend en compte la phase glace
923  if (.not.ice_thermo) then
924  d_ql(i,k) = zoliq(i)
925  d_qi(i,k) = 0.
926  else
927  d_ql(i,k) = (1-zfice(i))*zoliq(i)
928  d_qi(i,k) = zfice(i)*zoliq(i)
929  endif
930 !AJ<
931  zrfl(i) = zrfl(i)+ max(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
932  *(paprs(i,k)-paprs(i,k+1))/(rg*dtime)
933  zifl(i) = zifl(i)+ max(zcond(i)*zfice(i)-zoliqi(i),0.0) &
934  *(paprs(i,k)-paprs(i,k+1))/(rg*dtime)
935  ! zrfl(i) = zrfl(i)+ zpluie &
936  ! *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
937  ! zifl(i) = zifl(i)+ zice &
938  ! *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
939 
940  ENDIF
941  ENDDO
942  ENDIF
943 
944 !CR: la fonte est faite au debut
945 ! IF (ice_thermo) THEN
946 ! DO i = 1, klon
947 ! zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
948 ! zmelt = MIN(MAX(zmelt,0.),1.)
949 ! zrfl(i)=zrfl(i)+zmelt*zifl(i)
950 ! zifl(i)=zifl(i)*(1.-zmelt)
951 ! print*,zt(i),'octavio1'
952 ! zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
953 ! *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
954 ! print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
955 ! ENDDO
956 ! ENDIF
957 
958 
959  IF (.NOT. ice_thermo) THEN
960  DO i = 1, klon
961  IF (zt(i).LT.rtt) THEN
962  psfl(i,k)=zrfl(i)
963  ELSE
964  prfl(i,k)=zrfl(i)
965  ENDIF
966  ENDDO
967  ELSE
968  ! JAM*************************************************
969  ! Revoir partie ci-dessous: à quoi servent psfl et prfl?
970  ! *****************************************************
971 
972  DO i = 1, klon
973  ! IF (zt(i).LT.RTT) THEN
974  psfl(i,k)=zifl(i)
975  ! ELSE
976  prfl(i,k)=zrfl(i)
977  ! ENDIF
979  ENDDO
980  ENDIF
981  !
982  !
983  ! Calculer les tendances de q et de t:
984  !
985  DO i = 1, klon
986  d_q(i,k) = zq(i) - q(i,k)
987  d_t(i,k) = zt(i) - t(i,k)
988  ENDDO
989  !
990  !AA--------------- Calcul du lessivage stratiforme -------------
991 
992  DO i = 1,klon
993  !
994  if(zcond(i).gt.zoliq(i)+1.e-10) then
995  beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
996  else
997  beta(i,k) = 0.
998  endif
999  zprec_cond(i) = max(zcond(i)-zoliq(i),0.0) &
1000  * (paprs(i,k)-paprs(i,k+1))/rg
1001  IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
1002  !AA lessivage nucleation LMD5 dans la couche elle-meme
1003  IF (iflag_t_glace.EQ.0) THEN
1004  if (t(i,k) .GE. t_glace_min_old) THEN
1005  zalpha_tr = a_tr_sca(3)
1006  else
1007  zalpha_tr = a_tr_sca(4)
1008  endif
1009  ELSE ! of IF (iflag_t_glace.EQ.0)
1010  if (t(i,k) .GE. t_glace_min) THEN
1011  zalpha_tr = a_tr_sca(3)
1012  else
1013  zalpha_tr = a_tr_sca(4)
1014  endif
1015  ENDIF
1016  zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))
1017  pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
1018  frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
1019  !
1020  ! nucleation avec un facteur -1 au lieu de -0.5
1021  zfrac_lessi = 1. - exp(-zprec_cond(i)/zneb(i))
1022  pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
1023  ENDIF
1024  !
1025  ENDDO ! boucle sur i
1026  !
1027  !AA Lessivage par impaction dans les couches en-dessous
1028  DO kk = k-1, 1, -1
1029  DO i = 1, klon
1030  IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
1031  IF (iflag_t_glace.EQ.0) THEN
1032  if (t(i,kk) .GE. t_glace_min_old) THEN
1033  zalpha_tr = a_tr_sca(1)
1034  else
1035  zalpha_tr = a_tr_sca(2)
1036  endif
1037  ELSE ! of IF (iflag_t_glace.EQ.0)
1038  if (t(i,kk) .GE. t_glace_min) THEN
1039  zalpha_tr = a_tr_sca(1)
1040  else
1041  zalpha_tr = a_tr_sca(2)
1042  endif
1043  ENDIF
1044  zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))
1045  pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
1046  frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
1047  ENDIF
1048  ENDDO
1049  ENDDO
1050  !
1051  !AA----------------------------------------------------------
1052  ! FIN DE BOUCLE SUR K
1053  end DO
1054  !
1055  !AA-----------------------------------------------------------
1056  !
1057  ! Pluie ou neige au sol selon la temperature de la 1ere couche
1058  !
1059 !CR: si la thermo de la glace est active, on calcule zifl directement
1060  IF (.NOT.ice_thermo) THEN
1061  DO i = 1, klon
1062  IF ((t(i,1)+d_t(i,1)) .LT. rtt) THEN
1063 !AJ<
1064 ! snow(i) = zrfl(i)
1065  snow(i) = zrfl(i)+zifl(i)
1067  zlh_solid(i) = rlstt-rlvtt
1068  ELSE
1069  rain(i) = zrfl(i)
1070  zlh_solid(i) = 0.
1071  ENDIF
1072  ENDDO
1073 
1074  ELSE
1075  DO i = 1, klon
1076  snow(i) = zifl(i)
1077  rain(i) = zrfl(i)
1078  ENDDO
1079 
1080  ENDIF
1081  !
1082  ! For energy conservation : when snow is present, the solification
1083  ! latent heat is considered.
1084 !CR: si thermo de la glace, neige deja prise en compte
1085  IF (.not.ice_thermo) THEN
1086  DO k = 1, klev
1087  DO i = 1, klon
1088  zcpair=rcpd*(1.0+rvtmp2*(q(i,k)+d_q(i,k)))
1089  zmair=(paprs(i,k)-paprs(i,k+1))/rg
1090  zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
1091  d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair)
1092  END DO
1093  END DO
1094  ENDIF
1095  !
1096 
1097  if (ncoreczq>0) then
1098  WRITE(lunout,*)'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
1099  endif
1100 
1101 END SUBROUTINE fisrtilp
subroutine fisrtilp(dtime, paprs, pplay, t, q, ptconv, ratqs, d_t, d_q, d_ql, d_qi, rneb, radliq, rain, snow, pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, beta, prfl, psfl, rhcl, zqta, fraca, ztv, zpspsk, ztla, zthl, iflag_cld_th, iflag_ice_thermo)
Definition: fisrtilp.F90:12
subroutine icefrac_lsc(np, temp, sig, icefrac)
!$Id t_glace_min REAL exposant_glace REAL rei_max REAL coefw_cld_cv REAL tmax_fonte_cv INTEGER iflag_cld_cv common nuagecom && t_glace_min
Definition: nuage.h:4
integer, save klon
Definition: dimphy.F90:3
!$Id!Parameters for nlm real spfac!IM cf epmax real ptcrit real omtrain real dttrig real beta
Definition: cv30param.h:5
!$Id t_glace_max
Definition: nuage.h:4
integer, save klev
Definition: dimphy.F90:7
!$Id t_glace_min REAL exposant_glace REAL rei_max REAL coefw_cld_cv REAL tmax_fonte_cv INTEGER iflag_cld_cv common nuagecom exposant_glace
Definition: nuage.h:4
!$Id t_glace_min REAL exposant_glace REAL rei_max REAL coefw_cld_cv REAL tmax_fonte_cv INTEGER iflag_t_glace
Definition: nuage.h:4
!$Id cld_lc_con REAL cld_tau_con REAL ffallv_lsc
Definition: fisrtilp.h:10
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL false
Definition: calcul_STDlev.h:26
!$Id cld_lc_con REAL cld_tau_lsc
Definition: fisrtilp.h:10
integer, save klevm1
Definition: dimphy.F90:9
!FH On elimine toutes les clefs physiques dans la dynamique prt_level
!$Header!integer nvarmx parameter(nfmx=10, imx=200, jmx=150, lmx=200, nvarmx=1000) real xd(imx
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL true
subroutine cloudth(ngrid, klev, ind2, ztv, po, zqta, fraca, qcloud, ctot, zpspsk, paprs, ztla, zthl, ratqs, zqs, t)
Definition: cloudth.F90:5
Definition: dimphy.F90:1
!$Header!gestion des impressions de sorties et de débogage la sortie standard prt_level COMMON comprint lunout
Definition: iniprint.h:7
real rg
Definition: comcstphy.h:1