GCC Code Coverage Report

Directory: ./
File: phys/fisrtilp.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 383 585 65.5%
Branches: 224 392 57.1%
Line Branch Exec Source
1 ! $Id: fisrtilp.F90 3992 2021-10-21 14:46:18Z musat $
2 !
3 !
4 508518467 SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
5 480 d_t, d_q, d_ql, d_qi, rneb, radliq, rain, snow, &
6 480 pfrac_impa, pfrac_nucl, pfrac_1nucl, &
7 frac_impa, frac_nucl, beta, &
8 prfl, psfl, rhcl, zqta, fraca, &
9 ztv, zpspsk, ztla, zthl, iflag_cld_th, &
10 iflag_ice_thermo)
11
12 !
13 USE dimphy
14 USE icefrac_lsc_mod ! compute ice fraction (JBM 3/14)
15 USE print_control_mod, ONLY: prt_level, lunout
16 USE cloudth_mod
17 USE ioipsl_getin_p_mod, ONLY : getin_p
18 USE phys_local_var_mod, ONLY: ql_seri,qs_seri
19 USE phys_local_var_mod, ONLY: rneblsvol
20 ! flag to include modifications to ensure energy conservation (if flag >0)
21 USE add_phys_tend_mod, only : fl_cor_ebil
22 IMPLICIT none
23 !======================================================================
24 ! Auteur(s): Z.X. Li (LMD/CNRS)
25 ! Date: le 20 mars 1995
26 ! Objet: condensation et precipitation stratiforme.
27 ! schema de nuage
28 ! Fusion de fisrt (physique sursaturation, P. LeVan K. Laval)
29 ! et ilp (il pleut, L. Li)
30 ! Principales parties:
31 ! P0> Thermalisation des precipitations venant de la couche du dessus
32 ! P1> Evaporation de la precipitation (qui vient du niveau k+1)
33 ! P2> Formation du nuage (en k)
34 ! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
35 ! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour
36 ! les valeurs de T et Q initiales
37 ! P2.A.2> Prise en compte du couplage entre eau condensee et T.
38 ! P2.A.3> Calcul des valeures finales associees a la formation des nuages
39 ! P2.B> Nuage "tout ou rien"
40 ! P2.C> Prise en compte de la Chaleur latente apres formation nuage
41 ! P3> Formation de la precipitation (en k)
42 !======================================================================
43 ! JLD:
44 ! * Routine probablement fausse (au moins incoherente) si thermcep = .false.
45 ! * fl_cor_ebil doit etre > 0 ;
46 ! fl_cor_ebil= 0 pour reproduire anciens bugs
47 !======================================================================
48 include "YOMCST.h"
49 include "fisrtilp.h"
50 include "nuage.h" ! JBM (3/14)
51
52 !
53 ! Principaux inputs:
54 !
55 REAL, INTENT(IN) :: dtime ! intervalle du temps (s)
56 REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! pression a inter-couche
57 REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! pression au milieu de couche
58 REAL, DIMENSION(klon,klev), INTENT(IN) :: t ! temperature (K)
59 REAL, DIMENSION(klon,klev), INTENT(IN) :: q ! humidite specifique (kg/kg)
60 LOGICAL, DIMENSION(klon,klev), INTENT(IN) :: ptconv ! points ou le schema de conv. prof. est actif
61 INTEGER, INTENT(IN) :: iflag_cld_th
62 INTEGER, INTENT(IN) :: iflag_ice_thermo
63 !
64 ! Inputs lies aux thermiques
65 !
66 REAL, DIMENSION(klon,klev), INTENT(IN) :: ztv
67 REAL, DIMENSION(klon,klev), INTENT(IN) :: zqta, fraca
68 REAL, DIMENSION(klon,klev), INTENT(IN) :: zpspsk, ztla
69 REAL, DIMENSION(klon,klev), INTENT(IN) :: zthl
70 !
71 ! Input/output
72 REAL, DIMENSION(klon,klev), INTENT(INOUT):: ratqs ! determine la largeur de distribution de vapeur
73 !
74 ! Principaux outputs:
75 !
76 REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t ! incrementation de la temperature (K)
77 REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q ! incrementation de la vapeur d'eau
78 REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_ql ! incrementation de l'eau liquide
79 REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_qi ! incrementation de l'eau glace
80 REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneb ! fraction nuageuse
81 REAL, DIMENSION(klon,klev), INTENT(OUT) :: radliq ! eau liquide utilisee dans rayonnements
82 REAL, DIMENSION(klon,klev), INTENT(OUT) :: rhcl ! humidite relative en ciel clair
83 REAL, DIMENSION(klon), INTENT(OUT) :: rain
84 REAL, DIMENSION(klon), INTENT(OUT) :: snow
85 REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: prfl
86 REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: psfl
87
88 !AA
89 ! Coeffients de fraction lessivee : pour OFF-LINE
90 !
91 REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfrac_nucl
92 REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfrac_1nucl
93 REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfrac_impa
94 !
95 ! Fraction d'aerosols lessivee par impaction et par nucleation
96 ! POur ON-LINE
97 !
98 REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_impa
99 REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_nucl
100 !AA
101 ! --------------------------------------------------------------------------------
102 !
103 ! Options du programme:
104 !
105 REAL, SAVE :: seuil_neb=0.001 ! un nuage existe vraiment au-dela
106 !$OMP THREADPRIVATE(seuil_neb)
107
108 !<LTP
109 REAL smallestreal
110 REAL, SAVE :: rain_int_min=0.001 !intensité locale minimum pour la pluie avant diminution de la fraction précipitante associée = 0.001 mm/s
111 !>LTP
112 !$OMP THREADPRIVATE(rain_int_min)
113
114
115 INTEGER ninter ! sous-intervals pour la precipitation
116 PARAMETER (ninter=5)
117 INTEGER,SAVE :: iflag_evap_prec=1 ! evaporation de la pluie
118 !$OMP THREADPRIVATE(iflag_evap_prec)
119 !
120 LOGICAL cpartiel ! condensation partielle
121 PARAMETER (cpartiel=.TRUE.)
122 REAL t_coup
123 PARAMETER (t_coup=234.0)
124 REAL DDT0
125 PARAMETER (DDT0=.01)
126 REAL ztfondue
127 PARAMETER (ztfondue=278.15)
128 ! --------------------------------------------------------------------------------
129 !
130 ! Variables locales:
131 !
132 INTEGER i, k, n, kk
133 INTEGER,save::itap=0
134 !$OMP THREADPRIVATE(itap)
135
136 REAL qsl, qsi
137 real zct ,zcl
138 INTEGER ncoreczq
139 960 REAL ctot(klon,klev)
140 960 REAL ctot_vol(klon,klev)
141 960 REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
142 960 REAL zdqsdT_raw(klon)
143 960 REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
144
145 960 logical lognormale(klon)
146 logical ice_thermo
147 960 LOGICAL convergence(klon)
148 960 INTEGER n_i(klon), iter
149 REAL cste
150
151 960 real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
152 960 real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
153 real erf
154 960 REAL qcloud(klon)
155
156 960 REAL zrfl(klon), zrfln(klon), zqev, zqevt
157 !<LTP
158 960 REAL zrflclr(klon), zrflcld(klon)
159 960 REAL d_zrfl_clr_cld(klon), d_zifl_clr_cld(klon)
160 960 REAL d_zrfl_cld_clr(klon), d_zifl_cld_clr(klon)
161 !>LTP
162
163 960 REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
164 !<LTP
165 960 REAL ziflclr(klon), ziflcld(klon)
166 !>LTP
167 960 REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
168 960 REAL zoliqp(klon), zoliqi(klon)
169 960 REAL zt(klon)
170 ! JBM (3/14) nexpo is replaced by exposant_glace
171 ! REAL nexpo ! exponentiel pour glace/eau
172 ! INTEGER, PARAMETER :: nexpo=6
173 INTEGER exposant_glace_old
174 REAL t_glace_min_old
175 960 REAL zdz(klon),zrho(klon),ztot , zrhol(klon)
176 960 REAL zchau ,zfroi ,zfice(klon),zneb(klon),znebprecip(klon)
177 !<LTP
178 960 REAL znebprecipclr(klon), znebprecipcld(klon)
179 960 REAL tot_zneb(klon), tot_znebn(klon), d_tot_zneb(klon)
180 960 REAL d_znebprecip_clr_cld(klon), d_znebprecip_cld_clr(klon)
181 !>LTP
182
183 REAL zmelt, zpluie, zice
184 960 REAL dzfice(klon)
185 REAL zsolid
186 !!!!
187 ! Variables pour Bergeron
188 REAL zcp, coef1, DeltaT, Deltaq, Deltaqprecl
189 960 REAL zqpreci(klon), zqprecl(klon)
190 ! Variable pour conservation enegie des precipitations
191 960 REAL zmqc(klon)
192 !
193 LOGICAL appel1er
194 SAVE appel1er
195 !$OMP THREADPRIVATE(appel1er)
196 !
197 ! iflag_oldbug_fisrtilp=0 enleve le BUG par JYG : tglace_min -> tglace_max
198 ! iflag_oldbug_fisrtilp=1 ajoute le BUG
199 INTEGER,SAVE :: iflag_oldbug_fisrtilp=0 !=0 sans bug
200 !$OMP THREADPRIVATE(iflag_oldbug_fisrtilp)
201 !---------------------------------------------------------------
202 !
203 !AA Variables traceurs:
204 !AA Provisoire !!! Parametres alpha du lessivage
205 !AA A priori on a 4 scavenging # possibles
206 !
207 REAL a_tr_sca(4)
208 save a_tr_sca
209 !$OMP THREADPRIVATE(a_tr_sca)
210 !
211 ! Variables intermediaires
212 !
213 REAL zalpha_tr
214 REAL zfrac_lessi
215 960 REAL zprec_cond(klon)
216 !AA
217 ! RomP >>> 15 nov 2012
218 REAL beta(klon,klev) ! taux de conversion de l'eau cond
219 ! RomP <<<
220 960 REAL zmair(klon), zcpair, zcpeau
221 ! Pour la conversion eau-neige
222 960 REAL zlh_solid(klon), zm_solid
223 !---------------------------------------------------------------
224 !
225 ! Fonctions en ligne:
226 !
227 REAL fallvs,fallvc ! Vitesse de chute pour cristaux de glace
228 ! (Heymsfield & Donner, 1990)
229 REAL zzz
230
231 include "YOETHF.h"
232 include "FCTTRE.h"
233 fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
234 fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc
235 !
236 DATA appel1er /.TRUE./
237 !ym
238 !CR: pour iflag_ice_thermo=2, on active que la convection
239 ! ice_thermo = iflag_ice_thermo .GE. 1
240
241
242 480 itap=itap+1
243
2/2
✓ Branch 0 taken 477120 times.
✓ Branch 1 taken 480 times.
 477600 znebprecip(:)=0.
244
245 !<LTP
246 smallestreal=1.e-9
247
2/2
✓ Branch 0 taken 477120 times.
✓ Branch 1 taken 480 times.
 477600 znebprecipclr(:)=0.
248
2/2
✓ Branch 0 taken 477120 times.
✓ Branch 1 taken 480 times.
 477600 znebprecipcld(:)=0.
249 !>LTP
250
251 480 ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
252 zdelq=0.0
253
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 ctot_vol(1:klon,1:klev)=0.0
254
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 rneblsvol(1:klon,1:klev)=0.0
255
256
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 480 times.
 480 if (prt_level>9)write(lunout,*)'NUAGES4 A. JAM'
257
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 479 times.
 480 IF (appel1er) THEN
258 1 CALL getin_p('iflag_oldbug_fisrtilp',iflag_oldbug_fisrtilp)
259 1 CALL getin_p('iflag_evap_prec',iflag_evap_prec)
260 1 CALL getin_p('seuil_neb',seuil_neb)
261 !<LTP
262 1 CALL getin_p('rain_int_min',rain_int_min)
263 !>LTP
264 1 write(lunout,*)' iflag_oldbug_fisrtilp =',iflag_oldbug_fisrtilp
265 !
266 1 WRITE(lunout,*) 'fisrtilp, ninter:', ninter
267 1 WRITE(lunout,*) 'fisrtilp, iflag_evap_prec:', iflag_evap_prec
268 !<LTP
269 1 WRITE(lunout,*) 'fisrtilp, rain_int_min:', rain_int_min
270 !>LTP
271 1 WRITE(lunout,*) 'fisrtilp, cpartiel:', cpartiel
272 1 WRITE(lunout,*) 'FISRTILP VERSION LUDO'
273
274
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
275 1 WRITE(lunout,*) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
276 1 WRITE(lunout,*) 'Je prefere un sous-intervalle de 6 minutes'
277 ! CALL abort
278 ENDIF
279 1 appel1er = .FALSE.
280 !
281 !AA initialiation provisoire
282 1 a_tr_sca(1) = -0.5
283 1 a_tr_sca(2) = -0.5
284 1 a_tr_sca(3) = -0.5
285 1 a_tr_sca(4) = -0.5
286 !
287 !AA Initialisation a 1 des coefs des fractions lessivees
288 !
289 !cdir collapse
290
2/2
✓ Branch 0 taken 39 times.
✓ Branch 1 taken 1 times.
 40 DO k = 1, klev
291
2/2
✓ Branch 0 taken 38766 times.
✓ Branch 1 taken 39 times.
 38806 DO i = 1, klon
292 38766 pfrac_nucl(i,k)=1.
293 38766 pfrac_1nucl(i,k)=1.
294 38766 pfrac_impa(i,k)=1.
295 38805 beta(i,k)=0. !RomP initialisation
296 ENDDO
297 ENDDO
298
299 ENDIF ! test sur appel1er
300 !
301 !MAf Initialisation a 0 de zoliq
302 ! DO i = 1, klon
303 ! zoliq(i)=0.
304 ! ENDDO
305 ! Determiner les nuages froids par leur temperature
306 ! nexpo regle la raideur de la transition eau liquide / eau glace.
307 !
308 !CR: on est oblige de definir des valeurs fisrt car les valeurs de newmicro ne sont pas les memes par defaut
309
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 480 times.
 480 IF (iflag_t_glace.EQ.0) THEN
310 ! ztglace = RTT - 15.0
311 ✗ t_glace_min_old = RTT - 15.0
312 !AJ<
313 ✗ IF (ice_thermo) THEN
314 ! nexpo = 2
315 exposant_glace_old = 2
316 ELSE
317 ! nexpo = 6
318 exposant_glace_old = 6
319 ENDIF
320
321 ENDIF
322
323 !! RLVTT = 2.501e6 ! pas de redefinition des constantes physiques (jyg)
324 !! RLSTT = 2.834e6 ! pas de redefinition des constantes physiques (jyg)
325 !>AJ
326 !cc nexpo = 1
327 !
328 ! Initialiser les sorties:
329 !
330 !cdir collapse
331
2/2
✓ Branch 0 taken 19200 times.
✓ Branch 1 taken 480 times.
 19680 DO k = 1, klev+1
332
2/2
✓ Branch 0 taken 19084800 times.
✓ Branch 1 taken 19200 times.
 19104480 DO i = 1, klon
333 19084800 prfl(i,k) = 0.0
334 19104000 psfl(i,k) = 0.0
335 ENDDO
336 ENDDO
337
338 !cdir collapse
339
340
2/2
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
 19200 DO k = 1, klev
341
2/2
✓ Branch 0 taken 18607680 times.
✓ Branch 1 taken 18720 times.
 18626880 DO i = 1, klon
342 18607680 d_t(i,k) = 0.0
343 18607680 d_q(i,k) = 0.0
344 18607680 d_ql(i,k) = 0.0
345 18607680 d_qi(i,k) = 0.0
346 18607680 rneb(i,k) = 0.0
347 18607680 radliq(i,k) = 0.0
348 18607680 frac_nucl(i,k) = 1.
349 18626400 frac_impa(i,k) = 1.
350 ENDDO
351 ENDDO
352
2/2
✓ Branch 0 taken 477120 times.
✓ Branch 1 taken 480 times.
 477600 DO i = 1, klon
353 477120 rain(i) = 0.0
354 477120 snow(i) = 0.0
355 477120 zoliq(i)=0.
356 ! ENDDO
357 !
358 ! Initialiser le flux de precipitation a zero
359 !
360 ! DO i = 1, klon
361 477120 zrfl(i) = 0.0
362 477120 zifl(i) = 0.0
363 !<LTP
364 477120 zrflclr(i) = 0.0
365 477120 ziflclr(i) = 0.0
366 477120 zrflcld(i) = 0.0
367 477120 ziflcld(i) = 0.0
368 477120 tot_zneb(i) = 0.0
369 477120 tot_znebn(i) = 0.0
370 477120 d_tot_zneb(i) = 0.0
371 !>LTP
372
373 477600 zneb(i) = seuil_neb
374 ENDDO
375 !
376 !
377 !AA Pour plus de securite
378
379 zalpha_tr = 0.
380 zfrac_lessi = 0.
381
382 !AA==================================================================
383 !
384 480 ncoreczq=0
385 ! BOUCLE VERTICALE (DU HAUT VERS LE BAS)
386 !
387
2/2
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
 19200 DO k = klev, 1, -1
388 !
389 !AA===============================================================
390 !
391 ! Initialisation temperature et vapeur
392
2/2
✓ Branch 0 taken 18607680 times.
✓ Branch 1 taken 18720 times.
 18626400 DO i = 1, klon
393 18607680 zt(i)=t(i,k)
394 18626400 zq(i)=q(i,k)
395 ENDDO
396 !
397 ! ----------------------------------------------------------------
398 ! P0> Thermalisation des precipitations venant de la couche du dessus
399 ! ----------------------------------------------------------------
400 ! Calculer la varition de temp. de l'air du a la chaleur sensible
401 ! transporter par la pluie. On thermalise la pluie avec l'air de la couche.
402 ! Cette quantite de pluie qui est thermalisee, et devra continue a l'etre lors
403 ! des differentes transformations thermodynamiques. Cette masse d'eau doit
404 ! donc etre ajoute a l'humidite de la couche lorsque l'on calcule la variation
405 ! de l'enthalpie de la couche avec la temperature
406 ! Variables calculees ou modifiees:
407 ! - zt: temperature de la cocuhe
408 ! - zmqc: masse de precip qui doit etre thermalisee
409 !
410
2/2
✓ Branch 0 taken 18240 times.
✓ Branch 1 taken 480 times.
 18720 IF(k.LE.klevm1) THEN
411
2/2
✓ Branch 0 taken 18130560 times.
✓ Branch 1 taken 18240 times.
 18148800 DO i = 1, klon
412 !IM
413 18130560 zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
414 ! il n'y a pas encore d'eau liquide ni glace dans la maiille, donc zq suffit
415 18130560 zcpair=RCPD*(1.0+RVTMP2*zq(i))
416 18130560 zcpeau=RCPD*RVTMP2
417
1/2
✓ Branch 0 taken 18130560 times.
✗ Branch 1 not taken.
 18148800 if (fl_cor_ebil .GT. 0) then
418 ! zmqc: masse de precip qui doit etre thermalisee avec l'air de la couche atm
419 ! pour s'assurer que la precip arrivant au sol aura bien la temperature de la
420 ! derniere couche
421 18130560 zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair(i)
422 ! t(i,k+1)+d_t(i,k+1): nouvelle temp de la couche au dessus
423 zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) &
424 18130560 / (zcpair + zmqc(i)*zcpeau)
425 else ! si on maintient les anciennes erreurs
426 zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau &
427 + zmair(i)*zcpair*zt(i) ) &
428 ✗ / (zmair(i)*zcpair + zrfl(i)*dtime*zcpeau)
429 end if
430 ENDDO
431 ELSE ! IF(k.LE.klevm1)
432
2/2
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 477120 times.
 477600 DO i = 1, klon
433 477120 zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
434 477600 zmqc(i) = 0.
435 ENDDO
436 ENDIF ! end IF(k.LE.klevm1)
437 !
438 ! ----------------------------------------------------------------
439 ! P1> Calcul de l'evaporation de la precipitation
440 ! ----------------------------------------------------------------
441 ! On evapore une partie des precipitations venant de la maille du dessus.
442 ! On calcule l'evaporation et la sublimation des precipitations, jusqu'a
443 ! ce que la fraction de cette couche qui est sous le nuage soit saturee.
444 ! Variables calculees ou modifiees:
445 ! - zrfl et zifl: flux de precip liquide et glace
446 ! - zq, zt: humidite et temperature de la cocuhe
447 ! - zmqc: masse de precip qui doit etre thermalisee
448 !
449
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18720 times.
 18720 IF (iflag_evap_prec>=1) THEN
450
2/2
✓ Branch 0 taken 18607680 times.
✓ Branch 1 taken 18720 times.
 18626400 DO i = 1, klon
451 ! S'il y a des precipitations
452
2/2
✓ Branch 0 taken 4842693 times.
✓ Branch 1 taken 13764987 times.
 18626400 IF (zrfl(i)+zifl(i).GT.0.) THEN
453 ! Calcul du qsat
454 IF (thermcep) THEN
455 4842693 zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
456 4842693 zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
457 4842693 zqs(i)=MIN(0.5,zqs(i))
458 4842693 zcor=1./(1.-RETV*zqs(i))
459 4842693 zqs(i)=zqs(i)*zcor
460 ELSE
461 IF (zt(i) .LT. t_coup) THEN
462 zqs(i) = qsats(zt(i)) / pplay(i,k)
463 ELSE
464 zqs(i) = qsatl(zt(i)) / pplay(i,k)
465 ENDIF
466 ENDIF
467 ENDIF ! (zrfl(i)+zifl(i).GT.0.)
468 ENDDO
469 !AJ<
470
471
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18720 times.
 18720 IF (.NOT. ice_thermo) THEN
472 ✗ DO i = 1, klon
473 ! S'il y a des precipitations
474 ✗ IF (zrfl(i)+zifl(i).GT.0.) THEN
475 ! Evap max pour ne pas saturer la fraction sous le nuage
476 ! Evap max jusqu'à atteindre la saturation dans la partie
477 ! de la maille qui est sous le nuage de la couche du dessus
478 !!! On ne tient compte de cette fraction que sous une seule
479 !!! couche sous le nuage
480 ✗ zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
481 ! Ajout de la prise en compte des precip a thermiser
482 ! avec petite reecriture
483 ✗ if (fl_cor_ebil .GT. 0) then ! nouveau
484 ! Calcul de l'evaporation du flux de precip herite
485 ! d'au-dessus
486 zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
487 ✗ * zmair(i)/pplay(i,k)*zt(i)*RD
488 ✗ zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) * dtime/zmair(i)
489
490 ! Seuil pour ne pas saturer la fraction sous le nuage
491 ✗ zqev = MIN (zqev, zqevt)
492 ! Nouveau flux de precip
493 ✗ zrfln(i) = zrfl(i) - zqev*zmair(i)/dtime
494 ! Aucun flux liquide pour T < t_coup, on reevapore tout.
495 ✗ IF (zt(i) .LT. t_coup.and.reevap_ice) THEN
496 ✗ zrfln(i)=0.
497 ✗ zqev = (zrfl(i)-zrfln(i))/zmair(i)*dtime
498 END IF
499 ! Nouvelle vapeur
500 ✗ zq(i) = zq(i) + zqev
501 ✗ zmqc(i) = zmqc(i)-zqev
502 ! Nouvelle temperature (chaleur latente)
503 zt(i) = zt(i) - zqev &
504 ✗ * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
505 !!JLD debut de partie a supprimer a terme
506 else ! if (fl_cor_ebil .GT. 0)
507 ! Calcul de l'evaporation du flux de precip herite
508 ! d'au-dessus
509 zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
510 ✗ * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
511 zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
512 ✗ * RG*dtime/(paprs(i,k)-paprs(i,k+1))
513 ! Seuil pour ne pas saturer la fraction sous le nuage
514 ✗ zqev = MIN (zqev, zqevt)
515 ! Nouveau flux de precip
516 zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
517 ✗ /RG/dtime
518 ! Aucun flux liquide pour T < t_coup
519 ✗ IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
520 ! Nouvelle vapeur
521 zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
522 ✗ * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
523 ! Nouvelle temperature (chaleur latente)
524 zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
525 * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
526 ✗ * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
527 end if ! if (fl_cor_ebil .GT. 0)
528 !!JLD fin de partie a supprimer a terme
529 ✗ zrfl(i) = zrfln(i)
530 ✗ zifl(i) = 0.
531 ENDIF ! (zrfl(i)+zifl(i).GT.0.)
532 ENDDO
533 !
534 ELSE ! (.NOT. ice_thermo)
535 ! ================================
536 ! Avec thermodynamique de la glace
537 ! ================================
538
2/2
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 18607680 times.
 18626400 DO i = 1, klon
539
540
541 !AJ<
542 ! S'il y a des precipitations
543
2/2
✓ Branch 0 taken 4842693 times.
✓ Branch 1 taken 13764987 times.
 18626400 IF (zrfl(i)+zifl(i).GT.0.) THEN
544
545 !LTP<
546 !On ne tient compte que du flux de précipitation en ciel clair dans le calcul de l'évaporation.
547
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4842693 times.
 4842693 IF (iflag_evap_prec==4) THEN
548 ✗ zrfl(i) = zrflclr(i)
549 ✗ zifl(i) = ziflclr(i)
550 ENDIF
551
552 !>LTP
553
554
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4842693 times.
 4842693 IF (iflag_evap_prec==1) THEN
555 ✗ znebprecip(i)=zneb(i)
556 ELSE
557 4842693 znebprecip(i)=MAX(zneb(i),znebprecip(i))
558 ENDIF
559
560
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4842693 times.
 4842693 IF (iflag_evap_prec==4) THEN
561 ! Evap max pour ne pas saturer toute la maille
562 ✗ zqev0 = MAX (0.0, zqs(i)-zq(i))
563 ELSE
564 ! Evap max pour ne pas saturer la fraction sous le nuage
565 4842693 zqev0 = MAX (0.0, (zqs(i)-zq(i))*znebprecip(i) )
566 ENDIF
567
568 !JAM
569 ! On differencie qsat pour l'eau et la glace
570 ! Si zdelta=1. --> glace
571 ! Si zdelta=0. --> eau liquide
572
573 ! Calcul du qsat par rapport a l'eau liquide
574 4842693 qsl= R2ES*FOEEW(zt(i),0.)/pplay(i,k)
575 4842693 qsl= MIN(0.5,qsl)
576 4842693 zcor= 1./(1.-RETV*qsl)
577 4842693 qsl= qsl*zcor
578
579 ! Calcul de l'evaporation du flux de precip venant du dessus
580 ! Formulation en racine du flux de precip
581 ! (Klemp & Wilhelmson, 1978; Sundqvist, 1988)
582
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4842693 times.
 4842693 IF (iflag_evap_prec==3) THEN
583 zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl) &
584 *SQRT(zrfl(i)/max(1.e-4,znebprecip(i))) &
585 ✗ *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
586 !<LTP
587
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4842693 times.
 4842693 ELSE IF (iflag_evap_prec==4) THEN
588 zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl) &
589 *SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) &
590 ✗ *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
591 !>LTP
592 ELSE
593 zqevt = 1.*coef_eva*(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
594 4842693 *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
595 ENDIF
596
597
598 zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
599 4842693 *RG*dtime/(paprs(i,k)-paprs(i,k+1))
600
601 ! Calcul du qsat par rapport a la glace
602 4842693 qsi= R2ES*FOEEW(zt(i),1.)/pplay(i,k)
603 4842693 qsi= MIN(0.5,qsi)
604 4842693 zcor= 1./(1.-RETV*qsi)
605 4842693 qsi= qsi*zcor
606
607 ! Calcul de la sublimation du flux de precip solide herite
608 ! d'au-dessus
609
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4842693 times.
 4842693 IF (iflag_evap_prec==3) THEN
610 zqevti = znebprecip(i)*coef_eva*(1.0-zq(i)/qsi) &
611 *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
612 ✗ *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
613 !<LTP
614
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4842693 times.
 4842693 ELSE IF (iflag_evap_prec==4) THEN
615 zqevti = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsi) &
616 *SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) &
617 ✗ *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
618 !>LTP
619 ELSE
620 zqevti = 1.*coef_eva*(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
621 4842693 *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
622 ENDIF
623 zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
624 4842693 *RG*dtime/(paprs(i,k)-paprs(i,k+1))
625
626
627 !JAM
628 ! Limitation de l'evaporation. On s'assure qu'on ne sature pas
629 ! la fraction de la couche sous le nuage sinon on repartit zqev0
630 ! en conservant la proportion liquide / glace
631
632
2/2
✓ Branch 0 taken 432804 times.
✓ Branch 1 taken 4409889 times.
 4842693 IF (zqevt+zqevti.GT.zqev0) THEN
633 432804 zqev=zqev0*zqevt/(zqevt+zqevti)
634 432804 zqevi=zqev0*zqevti/(zqevt+zqevti)
635 ELSE
636 !JLD je ne comprends pas les lignes ci-dessous. On repartit les precips
637 ! liquides et solides meme si on ne sature pas la couche.
638 ! A mon avis, le test est inutile, et il faudrait tout remplacer par:
639 ! zqev=zqevt
640 ! zqevi=zqevti
641
2/2
✓ Branch 0 taken 3578587 times.
✓ Branch 1 taken 831302 times.
 4409889 IF (zqevt+zqevti.GT.0.) THEN
642 3578587 zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
643 3578587 zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
644 ELSE
645 zqev=0.
646 zqevi=0.
647 ENDIF
648 ENDIF
649
650 ! Nouveaux flux de precip liquide et solide
651 zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
652 4842693 /RG/dtime)
653 zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
654 4842693 /RG/dtime)
655
656 ! Mise a jour de la vapeur, temperature et flux de precip
657 zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
658 4842693 * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
659
1/2
✓ Branch 0 taken 4842693 times.
✗ Branch 1 not taken.
 4842693 if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips
660 zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
661 4842693 * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
662 zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
663 * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
664 * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
665 + (zifln(i)-zifl(i)) &
666 * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
667 4842693 * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
668 else ! sans correction thermalisation des precips
669 zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
670 * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
671 * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
672 + (zifln(i)-zifl(i)) &
673 * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
674 ✗ * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
675 end if
676 ! Nouvelles vaeleurs des precips liquides et solides
677 4842693 zrfl(i) = zrfln(i)
678 4842693 zifl(i) = zifln(i)
679
680 !<LTP
681
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4842693 times.
 4842693 IF (iflag_evap_prec==4) THEN
682 ✗ zrflclr(i) = zrfl(i)
683 ✗ ziflclr(i) = zifl(i)
684 ✗ IF(zrflclr(i) + ziflclr(i) .LE. 0) THEN
685 ✗ znebprecipclr(i) = 0.
686 ENDIF
687 ✗ zrfl(i) = zrflclr(i) + zrflcld(i)
688 ✗ zifl(i) = ziflclr(i) + ziflcld(i)
689 ENDIF
690 !>LTP
691
692
693 ! print*,'REEVAP ',itap,k,znebprecip(1),zqev0,zqev,zqevi,zrfl(1)
694
695 !CR ATTENTION: deplacement de la fonte de la glace
696 !jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
697 !!! zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2 !!!!!!!!! jyg
698 !jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
699 4842693 zmelt = ((zt(i)-273.15)/(ztfondue-273.15)) ! jyg
700 ! fraction de la precip solide qui est fondue
701 4842693 zmelt = MIN(MAX(zmelt,0.),1.)
702 ! Fusion de la glace
703 !<LTP
704
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4842693 times.
 4842693 IF (iflag_evap_prec==4) THEN
705 ✗ zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)
706 ✗ zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)
707 ✗ zrfl(i)=zrflclr(i)+zrflcld(i)
708 !>LTP
709 ELSE
710 4842693 zrfl(i)=zrfl(i)+zmelt*zifl(i)
711 ENDIF
712
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4842693 times.
 4842693 if (fl_cor_ebil .LE. 0) then
713 ! the following line should not be here. Indeed, if zifl is modified
714 ! now, zifl(i)*zmelt is no more the amount of ice that has melt
715 ! and therefore the change in temperature computed below is wrong
716 ✗ zifl(i)=zifl(i)*(1.-zmelt)
717 end if
718 ! Chaleur latente de fusion
719
1/2
✓ Branch 0 taken 4842693 times.
✗ Branch 1 not taken.
 4842693 if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips
720 zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
721 4842693 *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
722 else ! sans correction thermalisation des precips
723 zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
724 ✗ *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
725 end if
726
1/2
✓ Branch 0 taken 4842693 times.
✗ Branch 1 not taken.
 4842693 if (fl_cor_ebil .GT. 0) then ! correction bug, deplacement ligne precedente
727 !<LTP
728
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4842693 times.
 4842693 IF (iflag_evap_prec==4) THEN
729 ✗ ziflclr(i)=ziflclr(i)*(1.-zmelt)
730 ✗ ziflcld(i)=ziflcld(i)*(1.-zmelt)
731 ✗ zifl(i)=ziflclr(i)+ziflcld(i)
732 !>LTP
733 ELSE
734 4842693 zifl(i)=zifl(i)*(1.-zmelt)
735 ENDIF
736 end if
737
738 ELSE
739 ! Si on n'a plus de pluies, on reinitialise a 0 la farcion
740 ! sous nuageuse utilisee pour la pluie.
741 13764987 znebprecip(i)=0.
742 ENDIF ! (zrfl(i)+zifl(i).GT.0.)
743 ENDDO
744
745 ENDIF ! (.NOT. ice_thermo)
746
747 ! ----------------------------------------------------------------
748 ! Fin evaporation de la precipitation
749 ! ----------------------------------------------------------------
750 ENDIF ! (iflag_evap_prec>=1)
751 !
752 ! Calculer Qs et L/Cp*dQs/dT:
753 !
754 IF (thermcep) THEN
755
2/2
✓ Branch 0 taken 18607680 times.
✓ Branch 1 taken 18720 times.
 18626400 DO i = 1, klon
756 18607680 zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
757 18607680 zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
758
1/2
✓ Branch 0 taken 18607680 times.
✗ Branch 1 not taken.
 18607680 if (fl_cor_ebil .GT. 0) then ! nouveau
759 18607680 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
760 else
761 ✗ zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
762 endif
763 18607680 zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
764 18607680 zqs(i) = MIN(0.5,zqs(i))
765 18607680 zcor = 1./(1.-RETV*zqs(i))
766 18607680 zqs(i) = zqs(i)*zcor
767 18607680 zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
768 zdqsdT_raw(i) = zdqs(i)* &
769 18626400 & RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
770 ENDDO
771 ELSE
772 DO i = 1, klon
773 IF (zt(i).LT.t_coup) THEN
774 zqs(i) = qsats(zt(i))/pplay(i,k)
775 zdqs(i) = dqsats(zt(i),zqs(i))
776 ELSE
777 zqs(i) = qsatl(zt(i))/pplay(i,k)
778 zdqs(i) = dqsatl(zt(i),zqs(i))
779 ENDIF
780 ENDDO
781 ENDIF
782 !
783 ! Determiner la condensation partielle et calculer la quantite
784 ! de l'eau condensee:
785 !
786 !verification de la valeur de iflag_fisrtilp_qsat pour iflag_ice_thermo=1
787 ! if ((iflag_ice_thermo.eq.1).and.(iflag_fisrtilp_qsat.ne.0)) then
788 ! write(*,*) " iflag_ice_thermo==1 requires iflag_fisrtilp_qsat==0", &
789 ! " but iflag_fisrtilp_qsat=",iflag_fisrtilp_qsat, ". Might as well stop here."
790 ! stop
791 ! endif
792
793 ! ----------------------------------------------------------------
794 ! P2> Formation du nuage
795 ! ----------------------------------------------------------------
796 ! Variables calculees:
797 ! rneb : fraction nuageuse
798 ! zcond : eau condensee moyenne dans la maille.
799 ! rhcl: humidite relative ciel-clair
800 ! zt : temperature de la maille
801 ! ----------------------------------------------------------------
802 !
803 IF (cpartiel) THEN
804 ! -------------------------
805 ! P2.A> Nuage fractionnaire
806 ! -------------------------
807 !
808 ! Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
809 ! nuageuse a partir des PDF de Sandrine Bony.
810 ! rneb : fraction nuageuse
811 ! zqn : eau totale dans le nuage
812 ! zcond : eau condensee moyenne dans la maille.
813 ! on prend en compte le réchauffement qui diminue la partie
814 ! condensee
815 !
816 ! Version avec les raqts
817
818 ! ----------------------------------------------------------------
819 ! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
820 ! ----------------------------------------------------------------
821
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18720 times.
 18720 if (iflag_pdf.eq.0) then
822
823 ! version creneau de (Li, 1998)
824 ✗ do i=1,klon
825 ✗ zdelq = min(ratqs(i,k),0.99) * zq(i)
826 ✗ rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
827 ✗ zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
828 enddo
829
830 else ! if (iflag_pdf.eq.0)
831 ! ----------------------------------------------------------------
832 ! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour
833 ! les valeurs de T et Q initiales
834 ! ----------------------------------------------------------------
835
2/2
✓ Branch 0 taken 18607680 times.
✓ Branch 1 taken 18720 times.
 18626400 do i=1,klon
836
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18607680 times.
 18626400 if(zq(i).lt.1.e-15) then
837 ✗ ncoreczq=ncoreczq+1
838 ✗ zq(i)=1.e-15
839 endif
840 enddo
841
842
1/2
✓ Branch 0 taken 18720 times.
✗ Branch 1 not taken.
 18720 if (iflag_cld_th>=5) then
843
844
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18720 times.
 18720 if (iflag_cloudth_vert<=2) then
845 call cloudth(klon,klev,k,ztv, &
846 zq,zqta,fraca, &
847 qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
848 ✗ ratqs,zqs,t)
849
1/2
✓ Branch 0 taken 18720 times.
✗ Branch 1 not taken.
 18720 elseif (iflag_cloudth_vert>=3 .and. iflag_cloudth_vert<=5) then
850 call cloudth_v3(klon,klev,k,ztv, &
851 zq,zqta,fraca, &
852 qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
853 18720 ratqs,zqs,t)
854 !----------------------------------
855 !Version these Jean Jouhaud, Decembre 2018
856 !----------------------------------
857 ✗ elseif (iflag_cloudth_vert==6) then
858 call cloudth_v6(klon,klev,k,ztv, &
859 zq,zqta,fraca, &
860 qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
861 ✗ ratqs,zqs,t)
862
863 endif
864
2/2
✓ Branch 0 taken 18607680 times.
✓ Branch 1 taken 18720 times.
 18626400 do i=1,klon
865 18607680 rneb(i,k)=ctot(i,k)
866 18607680 rneblsvol(i,k)=ctot_vol(i,k)
867 18626400 zqn(i)=qcloud(i)
868 enddo
869
870 endif
871
872
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18720 times.
 18720 if (iflag_cld_th <= 4) then
873 ✗ lognormale = .true.
874
1/2
✓ Branch 0 taken 18720 times.
✗ Branch 1 not taken.
 18720 elseif (iflag_cld_th >= 6) then
875 ! lognormale en l'absence des thermiques
876
2/2
✓ Branch 0 taken 18607680 times.
✓ Branch 1 taken 18720 times.
 18626400 lognormale = fraca(:,k) < 1e-10
877 else
878 ! Dans le cas iflag_cld_th=5, on prend systématiquement la
879 ! bi-gaussienne
880 ✗ lognormale = .false.
881 end if
882
883 !CR: variation de qsat avec T en presence de glace ou non
884 !initialisations
885
2/2
✓ Branch 0 taken 18607680 times.
✓ Branch 1 taken 18720 times.
 18626400 do i=1,klon
886 18607680 DT(i) = 0.
887 18607680 n_i(i)=0
888 18607680 Tbef(i)=zt(i)
889 18626400 qlbef(i)=0.
890 enddo
891
892 ! ----------------------------------------------------------------
893 ! P2.A.2> Prise en compte du couplage entre eau condensee et T.
894 ! Calcul des grandeurs nuageuses en tenant compte de l'effet de
895 ! la condensation sur T, et donc sur qsat et sur les grandeurs nuageuses
896 ! qui en dependent. Ce changement de temperature est provisoire, et
897 ! la valeur definitive sera calcule plus tard.
898 ! Variables calculees:
899 ! rneb : nebulosite
900 ! zcond: eau condensee en moyenne dans la maille
901 ! note JLD: si on n'a pas de pdf lognormale, ce qui se passe ne me semble
902 ! pas clair, il n'y a probablement pas de prise en compte de l'effet de
903 ! T sur qsat
904 ! ----------------------------------------------------------------
905
906 !Boucle iterative: ATTENTION, l'option -1 n'est plus activable ici
907
1/2
✓ Branch 0 taken 18720 times.
✗ Branch 1 not taken.
 18720 if (iflag_fisrtilp_qsat.ge.0) then
908 ! Iteration pour condensation avec variation de qsat(T)
909 ! -----------------------------------------------------
910
2/2
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 93600 times.
 112320 do iter=1,iflag_fisrtilp_qsat+1
911
912
2/2
✓ Branch 0 taken 93038400 times.
✓ Branch 1 taken 93600 times.
 93132000 do i=1,klon
913 ! do while ((abs(DT(i)).gt.DDT0).or.(n_i(i).eq.0))
914 ! !! convergence = .true. tant que l'on n'a pas converge !!
915 ! ------------------------------
916 93038400 convergence(i)=abs(DT(i)).gt.DDT0
917
6/6
✓ Branch 0 taken 91462338 times.
✓ Branch 1 taken 1576062 times.
✓ Branch 2 taken 23826200 times.
✓ Branch 3 taken 67636138 times.
✓ Branch 4 taken 18879112 times.
✓ Branch 5 taken 6523150 times.
 93132000 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
918 ! si on n'a pas converge
919 !
920 ! P2.A.2.1> Calcul de la fraction nuageuse et de la quantite d'eau condensee
921 ! ---------------------------------------------------------------
922 ! Variables calculees:
923 ! rneb : nebulosite
924 ! zqn : eau condensee, dans le nuage (in cloud water content)
925 ! zcond: eau condensee en moyenne dans la maille
926 ! rhcl: humidite relative ciel-clair
927 !
928 18879112 Tbef(i)=Tbef(i)+DT(i) ! nouvelle temperature
929
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18879112 times.
 18879112 if (.not.ice_thermo) then
930 ✗ zdelta = MAX(0.,SIGN(1.,RTT-Tbef(i)))
931 else
932
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18879112 times.
 18879112 if (iflag_t_glace.eq.0) then
933 ✗ zdelta = MAX(0.,SIGN(1.,t_glace_min_old-Tbef(i)))
934
1/2
✓ Branch 0 taken 18879112 times.
✗ Branch 1 not taken.
 18879112 else if (iflag_t_glace.ge.1) then
935
1/2
✓ Branch 0 taken 18879112 times.
✗ Branch 1 not taken.
 18879112 if (iflag_oldbug_fisrtilp.EQ.0) then
936 18879112 zdelta = MAX(0.,SIGN(1.,t_glace_max-Tbef(i)))
937 else
938 !avec bug : zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
939 ✗ zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
940 endif
941 endif
942 endif
943 ! Calcul de rneb, qzn et zcond pour les PDF lognormales
944 18879112 zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
945
1/2
✓ Branch 0 taken 18879112 times.
✗ Branch 1 not taken.
 18879112 if (fl_cor_ebil .GT. 0) then
946 18879112 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
947 else
948 ✗ zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
949 end if
950 18879112 zqs(i) = R2ES*FOEEW(Tbef(i),zdelta)/pplay(i,k)
951 18879112 zqs(i) = MIN(0.5,zqs(i))
952 18879112 zcor = 1./(1.-RETV*zqs(i))
953 18879112 zqs(i) = zqs(i)*zcor
954 18879112 zdqs(i) = FOEDE(Tbef(i),zdelta,zcvm5,zqs(i),zcor)
955 18879112 zpdf_sig(i)=ratqs(i,k)*zq(i)
956 18879112 zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
957 18879112 zpdf_delta(i)=log(zq(i)/zqs(i))
958 18879112 zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
959 18879112 zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
960 18879112 zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
961 18879112 zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
962 18879112 zpdf_e1(i)=1.-erf(zpdf_e1(i))
963 18879112 zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
964 18879112 zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
965 18879112 zpdf_e2(i)=1.-erf(zpdf_e2(i))
966
967
2/2
✓ Branch 0 taken 11958421 times.
✓ Branch 1 taken 6920691 times.
 18879112 if (zpdf_e1(i).lt.1.e-10) then
968 11958421 rneb(i,k)=0.
969 11958421 zqn(i)=zqs(i)
970 else
971 6920691 rneb(i,k)=0.5*zpdf_e1(i)
972 6920691 zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
973 endif
974
975 ! If vertical heterogeneity, change fraction by volume as well
976
1/2
✓ Branch 0 taken 18879112 times.
✗ Branch 1 not taken.
 18879112 if (iflag_cloudth_vert>=3) then
977 18879112 ctot_vol(i,k)=rneb(i,k)
978 18879112 rneblsvol(i,k)=ctot_vol(i,k)
979 endif
980
981 endif !convergence
982
983 enddo ! boucle en i
984
985 ! P2.A.2.2> Calcul APPROCHE de la variation de temperature DT
986 ! due a la condensation.
987 ! ---------------------------------------------------------------
988 ! Variables calculees:
989 ! DT : variation de temperature due a la condensation
990
991
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 93600 times.
 112320 if (.not. ice_thermo) then
992 ! --------------------------
993 ✗ do i=1,klon
994 ✗ if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
995
996 ✗ qlbef(i)=max(0.,zqn(i)-zqs(i))
997 ✗ if (fl_cor_ebil .GT. 0) then
998 ✗ num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
999 else
1000 ✗ num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
1001 end if
1002 ✗ denom=1.+rneb(i,k)*zdqs(i)
1003 ✗ DT(i)=num/denom
1004 ✗ n_i(i)=n_i(i)+1
1005 endif
1006 enddo
1007
1008 else ! if (.not. ice_thermo)
1009 ! --------------------------
1010
1/2
✓ Branch 0 taken 93600 times.
✗ Branch 1 not taken.
 93600 if (iflag_t_glace.ge.1) then
1011
3/4
✓ Branch 0 taken 93600 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 93038400 times.
✓ Branch 3 taken 93600 times.
 93132000 CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
1012 endif
1013
1014
2/2
✓ Branch 0 taken 93600 times.
✓ Branch 1 taken 93038400 times.
 93132000 do i=1,klon
1015
6/6
✓ Branch 0 taken 91462338 times.
✓ Branch 1 taken 1576062 times.
✓ Branch 2 taken 23826200 times.
✓ Branch 3 taken 67636138 times.
✓ Branch 4 taken 18879112 times.
✓ Branch 5 taken 6523150 times.
 93132000 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
1016
1017
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18879112 times.
 18879112 if (iflag_t_glace.eq.0) then
1018 ✗ zfice(i) = 1.0 - (Tbef(i)-t_glace_min_old) / (RTT-t_glace_min_old)
1019 ✗ zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1020 ✗ zfice(i) = zfice(i)**exposant_glace_old
1021 dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) &
1022 ✗ & / (t_glace_min_old - RTT)
1023 endif
1024
1025
3/4
✓ Branch 0 taken 18879112 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 16884652 times.
✓ Branch 3 taken 1994460 times.
 18879112 if (iflag_t_glace.ge.1.and.zfice(i)>0.) then
1026 dzfice(i)= exposant_glace * zfice(i)**(exposant_glace-1) &
1027 16884652 & / (t_glace_min - t_glace_max)
1028 endif
1029
1030
4/4
✓ Branch 0 taken 16884652 times.
✓ Branch 1 taken 1994460 times.
✓ Branch 2 taken 11794904 times.
✓ Branch 3 taken 5089748 times.
 18879112 if ((zfice(i).eq.0).or.(zfice(i).eq.1)) then
1031 13789364 dzfice(i)=0.
1032 endif
1033
1034
2/2
✓ Branch 0 taken 7084208 times.
✓ Branch 1 taken 11794904 times.
 18879112 if (zfice(i).lt.1) then
1035 7084208 cste=RLVTT
1036 else
1037 11794904 cste=RLSTT
1038 endif
1039
1040 18879112 qlbef(i)=max(0.,zqn(i)-zqs(i))
1041
1/2
✓ Branch 0 taken 18879112 times.
✗ Branch 1 not taken.
 18879112 if (fl_cor_ebil .GT. 0) then
1042 num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
1043 18879112 & +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
1044 denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
1045 -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k) &
1046 18879112 & *qlbef(i)*dzfice(i)
1047 else
1048 num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
1049 ✗ & +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
1050 denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
1051 ✗ -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))*rneb(i,k)*qlbef(i)*dzfice(i)
1052 end if
1053 18879112 DT(i)=num/denom
1054 18879112 n_i(i)=n_i(i)+1
1055
1056 endif ! fin convergence
1057 enddo ! fin boucle i
1058
1059 endif !ice_thermo
1060
1061 enddo ! iter=1,iflag_fisrtilp_qsat+1
1062 ! Fin d'iteration pour condensation avec variation de qsat(T)
1063 ! -----------------------------------------------------------
1064 endif ! if (iflag_fisrtilp_qsat.ge.0)
1065 ! ----------------------------------------------------------------
1066 ! Fin de P2.A.2> la prise en compte du couplage entre eau condensee et T
1067 ! ----------------------------------------------------------------
1068
1069 endif ! iflag_pdf
1070
1071 ! if (iflag_fisrtilp_qsat.eq.-1) then
1072 !------------------------------------------
1073 !CR: ATTENTION option fausse mais a existe:
1074 ! pour la re-activer, prendre iflag_fisrtilp_qsat=0 et
1075 ! activer les lignes suivantes:
1076 IF (1.eq.0) THEN
1077 DO i=1,klon
1078 IF (rneb(i,k) .LE. 0.0) THEN
1079 zqn(i) = 0.0
1080 rneb(i,k) = 0.0
1081 zcond(i) = 0.0
1082 rhcl(i,k)=zq(i)/zqs(i)
1083 ELSE IF (rneb(i,k) .GE. 1.0) THEN
1084 zqn(i) = zq(i)
1085 rneb(i,k) = 1.0
1086 zcond(i) = MAX(0.0,zqn(i)-zqs(i))/(1+zdqs(i))
1087 rhcl(i,k)=1.0
1088 ELSE
1089 zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1+zdqs(i))
1090 rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
1091 ENDIF
1092 ENDDO
1093 ENDIF
1094 !------------------------------------------
1095
1096 ! ELSE
1097 ! ----------------------------------------------------------------
1098 ! P2.A.3> Calcul des valeures finales associees a la formation des nuages
1099 ! Variables calculees:
1100 ! rneb : nebulosite
1101 ! zcond: eau condensee en moyenne dans la maille
1102 ! zq : eau vapeur dans la maille
1103 ! zt : temperature de la maille
1104 ! rhcl: humidite relative ciel-clair
1105 ! ----------------------------------------------------------------
1106 !
1107 ! Bornage de l'eau in-cloud (zqn) et de la fraction nuageuse (rneb)
1108 ! Calcule de l'eau condensee moyenne dans la maille (zcond),
1109 ! et de l'humidite relative ciel-clair (rhcl)
1110
2/2
✓ Branch 0 taken 18607680 times.
✓ Branch 1 taken 18720 times.
 18626400 DO i=1,klon
1111
2/2
✓ Branch 0 taken 12695925 times.
✓ Branch 1 taken 5911755 times.
 18626400 IF (rneb(i,k) .LE. 0.0) THEN
1112 12695925 zqn(i) = 0.0
1113 12695925 rneb(i,k) = 0.0
1114 12695925 zcond(i) = 0.0
1115 12695925 rhcl(i,k)=zq(i)/zqs(i)
1116
2/2
✓ Branch 0 taken 48234 times.
✓ Branch 1 taken 5863521 times.
 5911755 ELSE IF (rneb(i,k) .GE. 1.0) THEN
1117 48234 zqn(i) = zq(i)
1118 48234 rneb(i,k) = 1.0
1119 48234 zcond(i) = MAX(0.0,zqn(i)-zqs(i))
1120 48234 rhcl(i,k)=1.0
1121 ELSE
1122 5863521 zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
1123 5863521 rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
1124 ENDIF
1125 ENDDO
1126
1127
1128 ! If vertical heterogeneity, change fraction by volume as well
1129
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18720 times.
 18720 if (iflag_cloudth_vert>=3) then
1130
2/2
✓ Branch 0 taken 18607680 times.
✓ Branch 1 taken 18720 times.
 18626400 ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
1131
2/2
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 18607680 times.
 18626400 rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
1132 endif
1133
1134 ! ENDIF
1135
1136 ELSE ! de IF (cpartiel)
1137 ! -------------------------
1138 ! P2.B> Nuage "tout ou rien"
1139 ! -------------------------
1140 ! note JLD: attention, rhcl non calcule. Ca peut avoir des consequences?
1141 DO i = 1, klon
1142 IF (zq(i).GT.zqs(i)) THEN
1143 rneb(i,k) = 1.0
1144 ELSE
1145 rneb(i,k) = 0.0
1146 ENDIF
1147 zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
1148 ENDDO
1149 ENDIF ! de IF (cpartiel)
1150 !
1151 ! Mise a jour vapeur d'eau
1152 ! -------------------------
1153
2/2
✓ Branch 0 taken 18607680 times.
✓ Branch 1 taken 18720 times.
 18626400 DO i = 1, klon
1154 18626400 zq(i) = zq(i) - zcond(i)
1155 ! zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
1156 ENDDO
1157 !AJ<
1158 ! ------------------------------------
1159 ! P2.C> Prise en compte de la Chaleur latente apres formation nuage
1160 ! -------------------------------------
1161 ! Variable calcule:
1162 ! zt : temperature de la maille
1163 !
1164
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18720 times.
 18720 IF (.NOT. ice_thermo) THEN
1165 ✗ if (iflag_fisrtilp_qsat.lt.1) then
1166 ✗ DO i = 1, klon
1167 ✗ zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
1168 ENDDO
1169 else if (iflag_fisrtilp_qsat.gt.0) then
1170 ✗ DO i= 1, klon
1171 ✗ if (fl_cor_ebil .GT. 0) then
1172 ✗ zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1173 else
1174 ✗ zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
1175 end if
1176 ENDDO
1177 endif
1178 ELSE
1179
1/2
✓ Branch 0 taken 18720 times.
✗ Branch 1 not taken.
 18720 if (iflag_t_glace.ge.1) then
1180
3/4
✓ Branch 0 taken 18720 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626400 CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
1181 endif
1182
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18720 times.
 18720 if (iflag_fisrtilp_qsat.lt.1) then
1183 ✗ DO i = 1, klon
1184 ! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
1185 ! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
1186 ! t_glace_max, exposant_glace)
1187 ✗ if (iflag_t_glace.eq.0) then
1188 ✗ zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
1189 ✗ zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1190 ✗ zfice(i) = zfice(i)**exposant_glace_old
1191 endif
1192 zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
1193 ✗ +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
1194 ENDDO
1195 else
1196
2/2
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 18607680 times.
 18626400 DO i=1, klon
1197 ! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
1198 ! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
1199 ! t_glace_max, exposant_glace)
1200
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18607680 times.
 18607680 if (iflag_t_glace.eq.0) then
1201 ✗ zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
1202 ✗ zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1203 ✗ zfice(i) = zfice(i)**exposant_glace_old
1204 endif
1205
1/2
✓ Branch 0 taken 18607680 times.
✗ Branch 1 not taken.
 18626400 if (fl_cor_ebil .GT. 0) then
1206 zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
1207 & * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
1208 18607680 +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1209 else
1210 zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
1211 ✗ +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
1212 end if
1213 ENDDO
1214 endif
1215 ! print*,zt(i),zrfl(i),zifl(i),'temp1'
1216 ENDIF
1217 !>AJ
1218
1219 ! ----------------------------------------------------------------
1220 ! P3> Formation des precipitations
1221 ! ----------------------------------------------------------------
1222 !
1223 ! Partager l'eau condensee en precipitation et eau liquide nuageuse
1224 !
1225
1226 !<LTP
1227
1228
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18720 times.
 18720 IF (iflag_evap_prec==4) THEN
1229 !Partitionnement des precipitations venant du dessus en précipitations nuageuses
1230 !et précipitations ciel clair
1231
1232 !0) Calculate tot_zneb, la fraction nuageuse totale au-dessus du nuage
1233 !en supposant un recouvrement maximum aléatoire (voir Jakob and Klein, 2000)
1234
1235 ✗ DO i=1, klon
1236 tot_znebn(i) = 1 - (1-tot_zneb(i))*(1 - max(rneb(i,k),zneb(i))) &
1237 ✗ /(1-min(zneb(i),1-smallestreal))
1238 ✗ d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i)
1239 ✗ tot_zneb(i) = tot_znebn(i)
1240
1241
1242 !1) Cloudy to clear air
1243 ✗ d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i,k),znebprecipcld(i))
1244 ✗ IF (znebprecipcld(i) .GT. 0) THEN
1245 ✗ d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i)
1246 ✗ d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i)
1247 ELSE
1248 ✗ d_zrfl_cld_clr(i) = 0.
1249 ✗ d_zifl_cld_clr(i) = 0.
1250 ENDIF
1251
1252 !2) Clear to cloudy air
1253 d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i,k) &
1254 ✗ - d_tot_zneb(i) - zneb(i)))
1255 ✗ IF (znebprecipclr(i) .GT. 0) THEN
1256 ✗ d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i)
1257 ✗ d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i)
1258 ELSE
1259 ✗ d_zrfl_clr_cld(i) = 0.
1260 ✗ d_zifl_clr_cld(i) = 0.
1261 ENDIF
1262
1263 !Update variables
1264 ✗ znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i)
1265 ✗ znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i)
1266 ✗ zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i)
1267 ✗ ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i)
1268 ✗ zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i)
1269 ✗ ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i)
1270
1271 ENDDO
1272 ENDIF
1273
1274 !>LTP
1275
1276
1277
1278 ! Initialisation de zoliq (eau condensee moyenne dans la maille)
1279
2/2
✓ Branch 0 taken 18607680 times.
✓ Branch 1 taken 18720 times.
 18626400 DO i = 1, klon
1280
2/2
✓ Branch 0 taken 5911755 times.
✓ Branch 1 taken 12695925 times.
 18626400 IF (rneb(i,k).GT.0.0) THEN
1281 5911755 zoliq(i) = zcond(i)
1282 5911755 zrho(i) = pplay(i,k) / zt(i) / RD
1283 5911755 zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
1284 ENDIF
1285 ENDDO
1286 !AJ<
1287
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18720 times.
 18720 IF (.NOT. ice_thermo) THEN
1288 ✗ IF (iflag_t_glace.EQ.0) THEN
1289 ✗ DO i = 1, klon
1290 ✗ IF (rneb(i,k).GT.0.0) THEN
1291 ✗ zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.13-t_glace_min_old)
1292 ✗ zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1293 ✗ zfice(i) = zfice(i)**exposant_glace_old
1294 ! zfice(i) = zfice(i)**nexpo
1295 !! zfice(i)=0.
1296 ENDIF
1297 ENDDO
1298 ELSE ! of IF (iflag_t_glace.EQ.0)
1299 ✗ CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
1300 ! DO i = 1, klon
1301 ! IF (rneb(i,k).GT.0.0) THEN
1302 ! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
1303 ! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
1304 ! t_glace_max, exposant_glace)
1305 ! ENDIF
1306 ! ENDDO
1307 ENDIF
1308 ENDIF
1309
1310 ! Calcul de radliq (eau condensee pour le rayonnement)
1311 ! Iteration pour realiser une moyenne de l'eau nuageuse lors de la precip
1312 ! Remarque: ce n'est donc pas l'eau restante en fin de precip mais une
1313 ! eau moyenne restante dans le nuage sur la duree du pas de temps qui est
1314 ! transmise au rayonnement;
1315 ! ----------------------------------------------------------------
1316
2/2
✓ Branch 0 taken 18607680 times.
✓ Branch 1 taken 18720 times.
 18626400 DO i = 1, klon
1317
2/2
✓ Branch 0 taken 5911755 times.
✓ Branch 1 taken 12695925 times.
 18626400 IF (rneb(i,k).GT.0.0) THEN
1318 5911755 zneb(i) = MAX(rneb(i,k), seuil_neb)
1319 ! zt(i) = zt(i)+zcond(i)*zfice(i)*RLMLT/RCPD/(1.0+RVTMP2*zq(i))
1320 ! print*,zt(i),'fractionglace'
1321 !>AJ
1322 5911755 radliq(i,k) = zoliq(i)/REAL(ninter+1)
1323 ENDIF
1324 ENDDO
1325 !
1326
2/2
✓ Branch 0 taken 93600 times.
✓ Branch 1 taken 18720 times.
 112320 DO n = 1, ninter
1327
2/2
✓ Branch 0 taken 93038400 times.
✓ Branch 1 taken 93600 times.
 93150720 DO i = 1, klon
1328
2/2
✓ Branch 0 taken 29558775 times.
✓ Branch 1 taken 63479625 times.
 93132000 IF (rneb(i,k).GT.0.0) THEN
1329 29558775 zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
1330 ! Initialization of zpluie and zice:
1331 zpluie=0
1332 zice=0
1333
2/2
✓ Branch 0 taken 20836225 times.
✓ Branch 1 taken 8722550 times.
 29558775 IF (zneb(i).EQ.seuil_neb) THEN
1334 ztot = 0.0
1335 ELSE
1336 ! quantite d'eau a eliminer: zchau (Sundqvist, 1978)
1337 ! meme chose pour la glace: zfroi (Zender & Kiehl, 1997)
1338
2/2
✓ Branch 0 taken 2522735 times.
✓ Branch 1 taken 18313490 times.
 20836225 if (ptconv(i,k)) then
1339 2522735 zcl =cld_lc_con
1340 2522735 zct =1./cld_tau_con
1341 zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
1342 2522735 *fallvc(zrhol(i)) * zfice(i)
1343 else
1344 18313490 zcl =cld_lc_lsc
1345 18313490 zct =1./cld_tau_lsc
1346 zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
1347 18313490 *fallvs(zrhol(i)) * zfice(i)
1348 endif
1349
1350 ! si l'heterogeneite verticale est active, on utilise
1351 ! la fraction volumique "vraie" plutot que la fraction
1352 ! surfacique modifiee, qui est plus grande et reduit
1353 ! sinon l'eau in-cloud de facon artificielle
1354
2/4
✓ Branch 0 taken 20836225 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 20836225 times.
 20836225 if ((iflag_cloudth_vert>=3).AND.(iflag_rain_incloud_vol==1)) then
1355 zchau = zct *dtime/REAL(ninter) * zoliq(i) &
1356 ✗ *(1.0-EXP(-(zoliq(i)/ctot_vol(i,k)/zcl )**2)) *(1.-zfice(i))
1357 else
1358 zchau = zct *dtime/REAL(ninter) * zoliq(i) &
1359 20836225 *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl )**2)) *(1.-zfice(i))
1360 endif
1361 !AJ<
1362
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 20836225 times.
 20836225 IF (.NOT. ice_thermo) THEN
1363 ✗ ztot = zchau + zfroi
1364 ELSE
1365 20836225 zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
1366 20836225 zice = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i))
1367 20836225 ztot = zpluie + zice
1368 ENDIF
1369 !>AJ
1370 20836225 ztot = MAX(ztot ,0.0)
1371 ENDIF
1372 29558775 ztot = MIN(ztot,zoliq(i))
1373 !AJ<
1374 ! zoliqp = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie , 0.0)
1375 ! zoliqi = MAX(zoliq(i)*zfice(i)-1.*zice , 0.0)
1376 !JLD : les 2 variables zoliqp et zoliqi crorresponent a des pseudo precip
1377 ! temporaires et ne doivent pas etre calcule (alors qu'elles le sont
1378 ! si iflag_bergeron <> 2
1379 ! A SUPPRIMER A TERME
1380 29558775 zoliqp(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie , 0.0)
1381 29558775 zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zice , 0.0)
1382 29558775 zoliq(i) = MAX(zoliq(i)-ztot , 0.0)
1383 !>AJ
1384 29558775 radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
1385 ENDIF
1386 ENDDO ! i = 1,klon
1387 ENDDO ! n = 1,ninter
1388
1389 ! ----------------------------------------------------------------
1390 !
1391
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18720 times.
 18720 IF (.NOT. ice_thermo) THEN
1392 ✗ DO i = 1, klon
1393 ✗ IF (rneb(i,k).GT.0.0) THEN
1394 ✗ d_ql(i,k) = zoliq(i)
1395 zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
1396 ✗ * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1397 ENDIF
1398 ENDDO
1399 ELSE
1400 !
1401 !CR&JYG<
1402 ! On prend en compte l'effet Bergeron dans les flux de precipitation :
1403 ! Si T < 0 C, alors les precipitations liquides sont converties en glace, ce qui
1404 ! provoque un accroissement de temperature DeltaT. L'effet de DeltaT sur le condensat
1405 ! et les precipitations est grossierement pris en compte en linearisant les equations
1406 ! et en approximant le processus de precipitation liquide par un processus a seuil.
1407 ! On fait l'hypothese que le condensat nuageux n'est pas modifié dans cette opération.
1408 ! Le condensat precipitant liquide est supprime (dans la limite DeltaT<273-T).
1409 ! Le condensat precipitant solide est augmente.
1410 ! La vapeur d'eau est augmentee.
1411 !
1412
1/2
✓ Branch 0 taken 18720 times.
✗ Branch 1 not taken.
 18720 IF ((iflag_bergeron .EQ. 2)) THEN
1413
2/2
✓ Branch 0 taken 18607680 times.
✓ Branch 1 taken 18720 times.
 18626400 DO i = 1, klon
1414
2/2
✓ Branch 0 taken 5911755 times.
✓ Branch 1 taken 12695925 times.
 18626400 IF (rneb(i,k) .GT. 0.0) THEN
1415 5911755 zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
1416 5911755 zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
1417
1/2
✓ Branch 0 taken 5911755 times.
✗ Branch 1 not taken.
 5911755 if (fl_cor_ebil .GT. 0) then
1418 5911755 zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1419 5911755 coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i)
1420 ! Calcul de DT si toute les precips liquides congelent
1421 5911755 DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
1422 ! On ne veut pas que T devienne superieur a la temp. de congelation.
1423 ! donc que Delta > RTT-zt(i
1424 5911755 DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
1425 5911755 zt(i) = zt(i) + DeltaT
1426 ! Eau vaporisee du fait de l'augmentation de T
1427 5911755 Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT
1428 ! on reajoute cette eau vaporise a la vapeur et on l'enleve des precips
1429 5911755 zq(i) = zq(i) + Deltaq
1430 ! Les 3 max si dessous prtotegent uniquement des erreurs d'arrondies
1431 5911755 zcond(i) = max( zcond(i)- Deltaq, 0. )
1432 ! precip liquide qui congele ou qui s'evapore
1433 5911755 Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
1434 5911755 zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
1435 ! bilan eau glacee
1436 5911755 zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
1437 else ! if (fl_cor_ebil .GT. 0)
1438 ! ancien calcul
1439 ✗ zcp=RCPD*(1.0+RVTMP2*(zq(i)+zcond(i)))
1440 ✗ coef1 = RLMLT*zdqs(i)/RLVTT
1441 ✗ DeltaT = max( min( RTT-zt(i), RLMLT*zqprecl(i)/zcp/(1.+coef1) ) , 0.)
1442 ✗ zqpreci(i) = zqpreci(i) + zcp/RLMLT*DeltaT
1443 ✗ zqprecl(i) = max( zqprecl(i) - zcp/RLMLT*(1.+coef1)*DeltaT, 0. )
1444 ✗ zcond(i) = max( zcond(i) - zcp/RLVTT*zdqs(i)*DeltaT, 0. )
1445 ✗ zq(i) = zq(i) + zcp/RLVTT*zdqs(i)*DeltaT
1446 ✗ zt(i) = zt(i) + DeltaT
1447 end if ! if (fl_cor_ebil .GT. 0)
1448 ENDIF ! rneb(i,k) .GT. 0.0
1449 ENDDO
1450
2/2
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 18607680 times.
 18626400 DO i = 1, klon
1451
2/2
✓ Branch 0 taken 5911755 times.
✓ Branch 1 taken 12695925 times.
 18626400 IF (rneb(i,k).GT.0.0) THEN
1452 5911755 d_ql(i,k) = (1-zfice(i))*zoliq(i)
1453 5911755 d_qi(i,k) = zfice(i)*zoliq(i)
1454 !<LTP
1455
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5911755 times.
 5911755 IF (iflag_evap_prec == 4) THEN
1456 zrflcld(i) = zrflcld(i)+zqprecl(i) &
1457 ✗ *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1458 ziflcld(i) = ziflcld(i)+ zqpreci(i) &
1459 ✗ *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1460 ✗ znebprecipcld(i) = rneb(i,k)
1461 ✗ zrfl(i) = zrflcld(i) + zrflclr(i)
1462 ✗ zifl(i) = ziflcld(i) + ziflclr(i)
1463 !>LTP
1464 ELSE
1465 zrfl(i) = zrfl(i)+ zqprecl(i) &
1466 5911755 *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1467 zifl(i) = zifl(i)+ zqpreci(i) &
1468 5911755 *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1469
1470 ENDIF !iflag_evap_prec==4
1471
1472 ENDIF
1473 ENDDO
1474 !!
1475 ELSE ! iflag_bergeron
1476 !>CR&JYG
1477 !!
1478 ✗ DO i = 1, klon
1479 ✗ IF (rneb(i,k).GT.0.0) THEN
1480 !CR on prend en compte la phase glace
1481 !JLD inutile car on ne passe jamais ici si .not.ice_thermo
1482 ! if (.not.ice_thermo) then
1483 ! d_ql(i,k) = zoliq(i)
1484 ! d_qi(i,k) = 0.
1485 ! else
1486 ✗ d_ql(i,k) = (1-zfice(i))*zoliq(i)
1487 ✗ d_qi(i,k) = zfice(i)*zoliq(i)
1488 ! endif
1489 !<LTP
1490 ✗ IF (iflag_evap_prec == 4) THEN
1491 zrflcld(i) = zrflcld(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
1492 ✗ *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1493 ziflcld(i) = ziflcld(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
1494 ✗ *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1495 ✗ znebprecipcld(i) = rneb(i,k)
1496 ✗ zrfl(i) = zrflcld(i) + zrflclr(i)
1497 ✗ zifl(i) = ziflcld(i) + ziflclr(i)
1498 !>LTP
1499 ELSE
1500 !AJ<
1501 zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
1502 ✗ *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1503 zifl(i) = zifl(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
1504 ✗ *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1505 ! zrfl(i) = zrfl(i)+ zpluie &
1506 ! *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1507 ! zifl(i) = zifl(i)+ zice &
1508 ! *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1509 ENDIF !iflag_evap_prec == 4
1510
1511 !CR : on prend en compte l'effet Bergeron dans les flux de precipitation
1512 ✗ IF ((iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)) THEN
1513 !<LTP
1514 ✗ IF (iflag_evap_prec == 4) THEN
1515 ✗ zsolid = zrfl(i)
1516 ✗ ziflclr(i) = ziflclr(i) +zrflclr(i)
1517 ✗ ziflcld(i) = ziflcld(i) +zrflcld(i)
1518 ✗ zifl(i) = ziflclr(i)+ziflcld(i)
1519 ✗ zrflcld(i)=0.
1520 ✗ zrflclr(i)=0.
1521 ✗ zrfl(i) = zrflclr(i)+zrflcld(i)
1522 !>LTP
1523 ELSE
1524 ✗ zsolid = zrfl(i)
1525 ✗ zifl(i) = zifl(i)+zrfl(i)
1526 ✗ zrfl(i) = 0.
1527 ENDIF!iflag_evap_prec==4
1528
1529 ✗ if (fl_cor_ebil .GT. 0) then
1530 zt(i)=zt(i)+zsolid*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
1531 ✗ *(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
1532 else
1533 zt(i)=zt(i)+zsolid*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
1534 ✗ *(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))
1535 end if
1536 ENDIF ! (iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)
1537 !RC
1538
1539 ENDIF ! rneb(i,k).GT.0.0
1540 ENDDO
1541
1542 ENDIF ! iflag_bergeron .EQ. 2
1543 ENDIF ! .NOT. ice_thermo
1544
1545 !CR: la fonte est faite au debut
1546 ! IF (ice_thermo) THEN
1547 ! DO i = 1, klon
1548 ! zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
1549 ! zmelt = MIN(MAX(zmelt,0.),1.)
1550 ! zrfl(i)=zrfl(i)+zmelt*zifl(i)
1551 ! zifl(i)=zifl(i)*(1.-zmelt)
1552 ! print*,zt(i),'octavio1'
1553 ! zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
1554 ! *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
1555 ! print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
1556 ! ENDDO
1557 ! ENDIF
1558
1559
1560 !<LTP
1561
1562 !Limitation de la fraction surfacique couverte par les précipitations lorsque l'intensité locale du flux de précipitation descend en
1563 !dessous de rain_int_min
1564
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18720 times.
 18720 IF (iflag_evap_prec==4) THEN
1565 ✗ DO i=1, klon
1566 ✗ IF (zrflclr(i) + ziflclr(i) .GT. 0 ) THEN
1567 ✗ znebprecipclr(i) = min(znebprecipclr(i), max(zrflclr(i)/(znebprecipclr(i)*rain_int_min), ziflclr(i)/(znebprecipclr(i)*rain_int_min)))
1568 ELSE
1569 ✗ znebprecipclr(i)=0.
1570 ENDIF
1571
1572 ✗ IF (zrflcld(i) + ziflcld(i) .GT. 0 ) THEN
1573 ✗ znebprecipcld(i) = min(znebprecipcld(i), max(zrflcld(i)/(znebprecipcld(i)*rain_int_min), ziflcld(i)/(znebprecipcld(i)*rain_int_min)))
1574 ELSE
1575 ✗ znebprecipcld(i)=0.
1576 ENDIF
1577 ENDDO
1578 ENDIf
1579
1580 !>LTP
1581
1582
1583
1584
1585
1586
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18720 times.
 18720 IF (.NOT. ice_thermo) THEN
1587 ✗ DO i = 1, klon
1588 ✗ IF (zt(i).LT.RTT) THEN
1589 ✗ psfl(i,k)=zrfl(i)
1590 ELSE
1591 ✗ prfl(i,k)=zrfl(i)
1592 ENDIF
1593 ENDDO
1594 ELSE
1595 ! JAM*************************************************
1596 ! Revoir partie ci-dessous: a quoi servent psfl et prfl?
1597 ! *****************************************************
1598
1599
2/2
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 18607680 times.
 18626400 DO i = 1, klon
1600 ! IF (zt(i).LT.RTT) THEN
1601 18607680 psfl(i,k)=zifl(i)
1602 ! ELSE
1603 18626400 prfl(i,k)=zrfl(i)
1604 ! ENDIF
1605 !>AJ
1606 ENDDO
1607 ENDIF
1608 ! ----------------------------------------------------------------
1609 ! Fin de formation des precipitations
1610 ! ----------------------------------------------------------------
1611 !
1612 ! Calculer les tendances de q et de t:
1613 !
1614
2/2
✓ Branch 0 taken 18607680 times.
✓ Branch 1 taken 18720 times.
 18626400 DO i = 1, klon
1615 18607680 d_q(i,k) = zq(i) - q(i,k)
1616 18626400 d_t(i,k) = zt(i) - t(i,k)
1617 ENDDO
1618 !
1619 !AA--------------- Calcul du lessivage stratiforme -------------
1620
1621
2/2
✓ Branch 0 taken 18607680 times.
✓ Branch 1 taken 18720 times.
 18626400 DO i = 1,klon
1622 !
1623
2/2
✓ Branch 0 taken 4104068 times.
✓ Branch 1 taken 14503612 times.
 18607680 if(zcond(i).gt.zoliq(i)+1.e-10) then
1624 4104068 beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1625 else
1626 14503612 beta(i,k) = 0.
1627 endif
1628 zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
1629 18607680 * (paprs(i,k)-paprs(i,k+1))/RG
1630
4/4
✓ Branch 0 taken 5911755 times.
✓ Branch 1 taken 12695925 times.
✓ Branch 2 taken 4167059 times.
✓ Branch 3 taken 1744696 times.
 18626400 IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
1631 !AA lessivage nucleation LMD5 dans la couche elle-meme
1632
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4167059 times.
 4167059 IF (iflag_t_glace.EQ.0) THEN
1633 ✗ if (t(i,k) .GE. t_glace_min_old) THEN
1634 ✗ zalpha_tr = a_tr_sca(3)
1635 else
1636 ✗ zalpha_tr = a_tr_sca(4)
1637 endif
1638 ELSE ! of IF (iflag_t_glace.EQ.0)
1639
2/2
✓ Branch 0 taken 1958503 times.
✓ Branch 1 taken 2208556 times.
 4167059 if (t(i,k) .GE. t_glace_min) THEN
1640 1958503 zalpha_tr = a_tr_sca(3)
1641 else
1642 2208556 zalpha_tr = a_tr_sca(4)
1643 endif
1644 ENDIF
1645 4167059 zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1646 4167059 pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
1647 4167059 frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
1648 !
1649 ! nucleation avec un facteur -1 au lieu de -0.5
1650 4167059 zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
1651 4167059 pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
1652 ENDIF
1653 !
1654 ENDDO ! boucle sur i
1655 !
1656 !AA Lessivage par impaction dans les couches en-dessous
1657
2/2
✓ Branch 0 taken 355680 times.
✓ Branch 1 taken 18720 times.
 374880 DO kk = k-1, 1, -1
1658
2/2
✓ Branch 0 taken 353545920 times.
✓ Branch 1 taken 355680 times.
 353920320 DO i = 1, klon
1659
4/4
✓ Branch 0 taken 71560548 times.
✓ Branch 1 taken 281985372 times.
✓ Branch 2 taken 43822325 times.
✓ Branch 3 taken 27738223 times.
 353901600 IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
1660
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 43822325 times.
 43822325 IF (iflag_t_glace.EQ.0) THEN
1661 ✗ if (t(i,kk) .GE. t_glace_min_old) THEN
1662 ✗ zalpha_tr = a_tr_sca(1)
1663 else
1664 ✗ zalpha_tr = a_tr_sca(2)
1665 endif
1666 ELSE ! of IF (iflag_t_glace.EQ.0)
1667
2/2
✓ Branch 0 taken 35399627 times.
✓ Branch 1 taken 8422698 times.
 43822325 if (t(i,kk) .GE. t_glace_min) THEN
1668 35399627 zalpha_tr = a_tr_sca(1)
1669 else
1670 8422698 zalpha_tr = a_tr_sca(2)
1671 endif
1672 ENDIF
1673 43822325 zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1674 43822325 pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
1675 43822325 frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
1676 ENDIF
1677 ENDDO
1678 ENDDO
1679 !
1680 !AA===============================================================
1681 ! FIN DE LA BOUCLE VERTICALE
1682 end DO
1683 !
1684 !AA==================================================================
1685 !
1686 ! Pluie ou neige au sol selon la temperature de la 1ere couche
1687 !
1688 !CR: si la thermo de la glace est active, on calcule zifl directement
1689
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 480 times.
 480 IF (.NOT.ice_thermo) THEN
1690 ✗ DO i = 1, klon
1691 ✗ IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
1692 !AJ<
1693 ! snow(i) = zrfl(i)
1694 ✗ snow(i) = zrfl(i)+zifl(i)
1695 !>AJ
1696 ✗ zlh_solid(i) = RLSTT-RLVTT
1697 ELSE
1698 ✗ rain(i) = zrfl(i)
1699 ✗ zlh_solid(i) = 0.
1700 ENDIF
1701 ENDDO
1702
1703 ELSE
1704
2/2
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 477120 times.
 477600 DO i = 1, klon
1705 477120 snow(i) = zifl(i)
1706 477600 rain(i) = zrfl(i)
1707 ENDDO
1708
1709 ENDIF
1710 !
1711 ! For energy conservation : when snow is present, the solification
1712 ! latent heat is considered.
1713 !CR: si thermo de la glace, neige deja prise en compte
1714
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 480 times.
 480 IF (.not.ice_thermo) THEN
1715 ✗ DO k = 1, klev
1716 ✗ DO i = 1, klon
1717 ✗ zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
1718 ✗ zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
1719 ✗ zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
1720 ✗ d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair(i))
1721 END DO
1722 END DO
1723 ENDIF
1724 !
1725
1726
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 480 times.
 480 if (ncoreczq>0) then
1727 ✗ WRITE(lunout,*)'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
1728 endif
1729
1730 480 END SUBROUTINE fisrtilp
1731