GCC Code Coverage Report

Directory: ./
File: phys/thermcell_plume_6A.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 165 361 45.7%
Branches: 182 396 46.0%
Line Branch Exec Source
1 !
2 ! $Id: thermcell_plume_6A.F90 3451 2019-01-27 11:07:30Z fhourdin $
3 !
4 3914853 SUBROUTINE thermcell_plume_6A(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, &
5 480 & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
6 480 & lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
7 480 & ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
8 & ,lev_out,lunout1,igout)
9 ! & ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
10 !--------------------------------------------------------------------------
11 !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
12 !--------------------------------------------------------------------------
13 USE IOIPSL, ONLY : getin
14 USE ioipsl_getin_p_mod, ONLY : getin_p
15
16 USE print_control_mod, ONLY: prt_level
17 IMPLICIT NONE
18
19 !
20 ! $Header$
21 !
22 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
23 ! veillez � n'utiliser que des ! pour les commentaires
24 ! et � bien positionner les & des lignes de continuation
25 ! (les placer en colonne 6 et en colonne 73)
26 !
27 !
28 ! A1.0 Fundamental constants
29 REAL RPI,RCLUM,RHPLA,RKBOL,RNAVO
30 ! A1.1 Astronomical constants
31 REAL RDAY,REA,REPSM,RSIYEA,RSIDAY,ROMEGA
32 ! A1.1.bis Constantes concernant l'orbite de la Terre:
33 REAL R_ecc, R_peri, R_incl
34 ! A1.2 Geoide
35 REAL RA,RG,R1SA
36 ! A1.3 Radiation
37 ! REAL RSIGMA,RI0
38 REAL RSIGMA
39 ! A1.4 Thermodynamic gas phase
40 REAL RMO3,RMCO2,RMC,RMCH4,RMN2O,RMCFC11,RMCFC12
41 REAL R,RMD,RMV,RD,RV,RCPD,RCPV,RCVD,RCVV
42 REAL RKAPPA,RETV, eps_w
43 ! A1.5,6 Thermodynamic liquid,solid phases
44 REAL RCW,RCS
45 ! A1.7 Thermodynamic transition of phase
46 REAL RLVTT,RLSTT,RLMLT,RTT,RATM
47 ! A1.8 Curve of saturation
48 REAL RESTT,RALPW,RBETW,RGAMW,RALPS,RBETS,RGAMS
49 REAL RALPD,RBETD,RGAMD
50 !
51 COMMON/YOMCST/RPI ,RCLUM ,RHPLA ,RKBOL ,RNAVO &
52 & ,RDAY ,REA ,REPSM ,RSIYEA,RSIDAY,ROMEGA &
53 & ,R_ecc, R_peri, R_incl &
54 & ,RA ,RG ,R1SA &
55 & ,RSIGMA &
56 & ,R ,RMD ,RMV ,RD ,RV ,RCPD &
57 & ,RMO3 ,RMCO2 ,RMC ,RMCH4 ,RMN2O ,RMCFC11 ,RMCFC12 &
58 & ,RCPV ,RCVD ,RCVV ,RKAPPA,RETV, eps_w &
59 & ,RCW ,RCS &
60 & ,RLVTT ,RLSTT ,RLMLT ,RTT ,RATM &
61 & ,RESTT ,RALPW ,RBETW ,RGAMW ,RALPS ,RBETS ,RGAMS &
62 & ,RALPD ,RBETD ,RGAMD
63 ! ------------------------------------------------------------------
64 !$OMP THREADPRIVATE(/YOMCST/)
65 !
66 ! $Id: YOETHF.h 2799 2017-02-24 18:50:33Z jyg $
67 !
68 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
69 ! veillez n'utiliser que des ! pour les commentaires
70 ! et bien positionner les & des lignes de continuation
71 ! (les placer en colonne 6 et en colonne 73)
72 !
73 !* COMMON *YOETHF* DERIVED CONSTANTS SPECIFIC TO ECMWF THERMODYNAMICS
74 !
75 ! *R__ES* *CONSTANTS USED FOR COMPUTATION OF SATURATION
76 ! MIXING RATIO OVER LIQUID WATER(*R_LES*) OR
77 ! ICE(*R_IES*).
78 ! *RVTMP2* *RVTMP2=RCPV/RCPD-1.
79 ! *RHOH2O* *DENSITY OF LIQUID WATER. (RATM/100.)
80 !
81 REAL R2ES, R3LES, R3IES, R4LES, R4IES, R5LES, R5IES
82 REAL RVTMP2, RHOH2O
83 REAL R5ALVCP,R5ALSCP,RALVDCP,RALSDCP,RALFDCP,RTWAT,RTBER,RTBERCU
84 REAL RTICE,RTICECU,RTWAT_RTICE_R,RTWAT_RTICECU_R,RKOOP1,RKOOP2
85 LOGICAL OK_BAD_ECMWF_THERMO ! If TRUE, then variables set by rrtm/suphec.F90
86 ! If FALSE, then variables set by suphel.F90
87 COMMON /YOETHF/R2ES, R3LES, R3IES, R4LES, R4IES, R5LES, R5IES, &
88 & RVTMP2, RHOH2O, &
89 & R5ALVCP,R5ALSCP,RALVDCP,RALSDCP, &
90 & RALFDCP,RTWAT,RTBER,RTBERCU, &
91 & RTICE,RTICECU,RTWAT_RTICE_R,RTWAT_RTICECU_R,RKOOP1,&
92 & RKOOP2, &
93 & OK_BAD_ECMWF_THERMO
94
95 !$OMP THREADPRIVATE(/YOETHF/)
96 !
97 ! $Header$
98 !
99 !
100 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
101 ! veillez n'utiliser que des ! pour les commentaires
102 ! et bien positionner les & des lignes de continuation
103 ! (les placer en colonne 6 et en colonne 73)
104 !
105 ! ------------------------------------------------------------------
106 ! This COMDECK includes the Thermodynamical functions for the cy39
107 ! ECMWF Physics package.
108 ! Consistent with YOMCST Basic physics constants, assuming the
109 ! partial pressure of water vapour is given by a first order
110 ! Taylor expansion of Qs(T) w.r.t. to Temperature, using constants
111 ! in YOETHF
112 ! ------------------------------------------------------------------
113 REAL PTARG, PDELARG, P5ARG, PQSARG, PCOARG
114 REAL FOEEW, FOEDE, qsats, qsatl, dqsats, dqsatl
115 LOGICAL thermcep
116 PARAMETER (thermcep=.TRUE.)
117 !
118 FOEEW ( PTARG,PDELARG ) = EXP ( &
119 & (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT) &
120 & / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) )
121 !
122 FOEDE ( PTARG,PDELARG,P5ARG,PQSARG,PCOARG ) = PQSARG*PCOARG*P5ARG &
123 & / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG))**2
124 !
125 qsats(ptarg) = 100.0 * 0.622 * 10.0 &
126 & ** (2.07023 - 0.00320991 * ptarg &
127 & - 2484.896 / ptarg + 3.56654 * LOG10(ptarg))
128 qsatl(ptarg) = 100.0 * 0.622 * 10.0 &
129 & ** (23.8319 - 2948.964 / ptarg &
130 & - 5.028 * LOG10(ptarg) &
131 & - 29810.16 * EXP( - 0.0699382 * ptarg) &
132 & + 25.21935 * EXP( - 2999.924 / ptarg))
133 !
134 dqsats(ptarg,pqsarg) = RLVTT/RCPD*pqsarg * (3.56654/ptarg &
135 & +2484.896*LOG(10.)/ptarg**2 &
136 & -0.00320991*LOG(10.))
137 dqsatl(ptarg,pqsarg) = RLVTT/RCPD*pqsarg*LOG(10.)* &
138 & (2948.964/ptarg**2-5.028/LOG(10.)/ptarg &
139 & +25.21935*2999.924/ptarg**2*EXP(-2999.924/ptarg) &
140 & +29810.16*0.0699382*EXP(-0.0699382*ptarg))
141 integer :: iflag_thermals,nsplit_thermals
142
143 !!! nrlmd le 10/04/2012
144 integer :: iflag_trig_bl,iflag_clos_bl
145 integer :: tau_trig_shallow,tau_trig_deep
146 real :: s_trig
147 !!! fin nrlmd le 10/04/2012
148
149 real,parameter :: r_aspect_thermals=2.,l_mix_thermals=30.
150 real :: alp_bl_k
151 real :: tau_thermals,fact_thermals_ed_dz
152 integer,parameter :: w2di_thermals=0
153 integer :: isplit
154
155 integer :: iflag_coupl,iflag_clos,iflag_wake
156 integer :: iflag_thermals_ed,iflag_thermals_optflux,iflag_thermals_closure
157
158 common/ctherm1/iflag_thermals,nsplit_thermals,iflag_thermals_closure
159 common/ctherm2/tau_thermals,alp_bl_k,fact_thermals_ed_dz
160 common/ctherm4/iflag_coupl,iflag_clos,iflag_wake
161 common/ctherm5/iflag_thermals_ed,iflag_thermals_optflux
162
163 !!! nrlmd le 10/04/2012
164 common/ctherm6/iflag_trig_bl,iflag_clos_bl
165 common/ctherm7/tau_trig_shallow,tau_trig_deep
166 common/ctherm8/s_trig
167 !!! fin nrlmd le 10/04/2012
168
169 !$OMP THREADPRIVATE(/ctherm1/,/ctherm2/,/ctherm4/,/ctherm5/)
170 !$OMP THREADPRIVATE(/ctherm6/,/ctherm7/,/ctherm8/)
171
172 INTEGER itap
173 INTEGER lunout1,igout
174 INTEGER ngrid,klev
175 REAL ptimestep
176 REAL ztv(ngrid,klev)
177 REAL zthl(ngrid,klev)
178 REAL po(ngrid,klev)
179 REAL zl(ngrid,klev)
180 REAL rhobarz(ngrid,klev)
181 REAL zlev(ngrid,klev+1)
182 REAL pplev(ngrid,klev+1)
183 REAL pphi(ngrid,klev)
184 REAL zpspsk(ngrid,klev)
185 REAL alim_star(ngrid,klev)
186 REAL f0(ngrid)
187 INTEGER lalim(ngrid)
188 integer lev_out ! niveau pour les print
189 integer nbpb
190
191 real alim_star_tot(ngrid)
192
193 REAL ztva(ngrid,klev)
194 REAL ztla(ngrid,klev)
195 REAL zqla(ngrid,klev)
196 REAL zqta(ngrid,klev)
197 REAL zha(ngrid,klev)
198
199 REAL detr_star(ngrid,klev)
200 REAL coefc
201 REAL entr_star(ngrid,klev)
202 960 REAL detr(ngrid,klev)
203 960 REAL entr(ngrid,klev)
204
205 REAL csc(ngrid,klev)
206
207 REAL zw2(ngrid,klev+1)
208 REAL w_est(ngrid,klev+1)
209 REAL f_star(ngrid,klev+1)
210 960 REAL wa_moy(ngrid,klev+1)
211
212 REAL ztva_est(ngrid,klev)
213 960 REAL ztv_est(ngrid,klev)
214 960 REAL zqla_est(ngrid,klev)
215 REAL zqsatth(ngrid,klev)
216 960 REAL zta_est(ngrid,klev)
217 960 REAL ztemp(ngrid),zqsat(ngrid)
218 REAL zdw2,zdw2bis
219 REAL zw2modif
220 REAL zw2fact,zw2factbis
221 960 REAL zeps(ngrid,klev)
222
223 REAL linter(ngrid)
224 INTEGER lmix(ngrid)
225 INTEGER lmix_bis(ngrid)
226 960 REAL wmaxa(ngrid)
227
228 INTEGER ig,l,k,lt,it,lm
229
230 960 real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
231 960 real zbuoyjam(ngrid,klev),zdqtjam(ngrid,klev)
232 real zbuoybis,zdz2,zdz3,lmel,entrbis,zdzbis
233 960 real d_temp(ngrid)
234 real ztv1,ztv2,factinv,zinv,zlmel
235 real zlmelup,zlmeldwn,zlt,zltdwn,zltup
236 real atv1,atv2,btv1,btv2
237 real ztv_est1,ztv_est2
238 real zcor,zdelta,zcvm5,qlbef
239 real zbetalpha, coefzlmel
240 real eps
241 REAL REPS,RLvCp,DDT0
242 PARAMETER (DDT0=.01)
243 logical Zsat
244 960 LOGICAL active(ngrid),activetmp(ngrid)
245 REAL fact_gamma,fact_gamma2,fact_epsilon2
246
247 REAL, SAVE :: fact_epsilon=0.002
248 REAL, SAVE :: betalpha=0.9
249 REAL, SAVE :: afact=2./3.
250 REAL, SAVE :: fact_shell=1.
251 REAL,SAVE :: detr_min=1.e-5
252 REAL,SAVE :: entr_min=1.e-5
253 REAL,SAVE :: detr_q_coef=0.012
254 REAL,SAVE :: detr_q_power=0.5
255 REAL,SAVE :: mix0=0.
256 INTEGER,SAVE :: thermals_flag_alim=0
257
258 !$OMP THREADPRIVATE(fact_epsilon, betalpha, afact, fact_shell)
259 !$OMP THREADPRIVATE(detr_min, entr_min, detr_q_coef, detr_q_power)
260 !$OMP THREADPRIVATE( mix0, thermals_flag_alim)
261
262 LOGICAL, SAVE :: first=.true.
263 !$OMP THREADPRIVATE(first)
264
265
266 REAL c2(ngrid,klev)
267
268
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 480 times.
 480 if (ngrid==1) print*,'THERMCELL PLUME MODIFIE 2014/07/11'
269 Zsat=.false.
270 ! Initialisation
271
272 480 RLvCp = RLVTT/RCPD
273
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 479 times.
 480 IF (first) THEN
274
275 1 CALL getin_p('thermals_fact_epsilon',fact_epsilon)
276 1 CALL getin_p('thermals_betalpha',betalpha)
277 1 CALL getin_p('thermals_afact',afact)
278 1 CALL getin_p('thermals_fact_shell',fact_shell)
279 1 CALL getin_p('thermals_detr_min',detr_min)
280 1 CALL getin_p('thermals_entr_min',entr_min)
281 1 CALL getin_p('thermals_detr_q_coef',detr_q_coef)
282 1 CALL getin_p('thermals_detr_q_power',detr_q_power)
283 1 CALL getin_p('thermals_mix0',mix0)
284 1 CALL getin_p('thermals_flag_alim',thermals_flag_alim)
285
286
287 1 first=.false.
288 ENDIF
289
290 480 zbetalpha=betalpha/(1.+betalpha)
291
292
293 ! Initialisations des variables r?elles
294 if (1==1) then
295
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 ztva(:,:)=ztv(:,:)
296
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 ztva_est(:,:)=ztva(:,:)
297
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 ztv_est(:,:)=ztv(:,:)
298
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 ztla(:,:)=zthl(:,:)
299
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 zqta(:,:)=po(:,:)
300
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 zqla(:,:)=0.
301
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 zha(:,:) = ztva(:,:)
302 else
303 ztva(:,:)=0.
304 ztv_est(:,:)=0.
305 ztva_est(:,:)=0.
306 ztla(:,:)=0.
307 zqta(:,:)=0.
308 zha(:,:) =0.
309 endif
310
311
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 zqla_est(:,:)=0.
312
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 zqsatth(:,:)=0.
313
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 zqla(:,:)=0.
314
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 detr_star(:,:)=0.
315
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 entr_star(:,:)=0.
316
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 alim_star(:,:)=0.
317
2/2
✓ Branch 0 taken 477120 times.
✓ Branch 1 taken 480 times.
 477600 alim_star_tot(:)=0.
318
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 csc(:,:)=0.
319
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 detr(:,:)=0.
320
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 entr(:,:)=0.
321
4/4
✓ Branch 0 taken 19200 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 19084800 times.
✓ Branch 3 taken 19200 times.
 19104480 zw2(:,:)=0.
322
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 zbuoy(:,:)=0.
323
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 zbuoyjam(:,:)=0.
324
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 gamma(:,:)=0.
325
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 18607680 times.
✓ Branch 3 taken 18720 times.
 18626880 zeps(:,:)=0.
326
4/4
✓ Branch 0 taken 19200 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 19084800 times.
✓ Branch 3 taken 19200 times.
 19104480 w_est(:,:)=0.
327
4/4
✓ Branch 0 taken 19200 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 19084800 times.
✓ Branch 3 taken 19200 times.
 19104480 f_star(:,:)=0.
328
4/4
✓ Branch 0 taken 19200 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 19084800 times.
✓ Branch 3 taken 19200 times.
 19104480 wa_moy(:,:)=0.
329
2/2
✓ Branch 0 taken 477120 times.
✓ Branch 1 taken 480 times.
 477600 linter(:)=1.
330 ! linter(:)=1.
331 ! Initialisation des variables entieres
332
2/2
✓ Branch 0 taken 477120 times.
✓ Branch 1 taken 480 times.
 477600 lmix(:)=1
333
2/2
✓ Branch 0 taken 477120 times.
✓ Branch 1 taken 480 times.
 477600 lmix_bis(:)=2
334
2/2
✓ Branch 0 taken 477120 times.
✓ Branch 1 taken 480 times.
 477600 wmaxa(:)=0.
335
336
337 !-------------------------------------------------------------------------
338 ! On ne considere comme actif que les colonnes dont les deux premieres
339 ! couches sont instables.
340 !-------------------------------------------------------------------------
341
342
2/2
✓ Branch 0 taken 477120 times.
✓ Branch 1 taken 480 times.
 477600 active(:)=ztv(:,1)>ztv(:,2)
343
2/2
✓ Branch 0 taken 477120 times.
✓ Branch 1 taken 480 times.
 477600 d_temp(:)=0. ! Pour activer un contraste de temperature a la base
344 ! du panache
345 ! Cet appel pourrait être fait avant thermcell_plume dans thermcell_main
346 480 CALL thermcell_alim(thermals_flag_alim,ngrid,klev,ztv,d_temp,zlev,alim_star,lalim)
347
348 !------------------------------------------------------------------------------
349 ! Calcul dans la premiere couche
350 ! On decide dans cette version que le thermique n'est actif que si la premiere
351 ! couche est instable.
352 ! Pourrait etre change si on veut que le thermiques puisse se d??clencher
353 ! dans une couche l>1
354 !------------------------------------------------------------------------------
355
2/2
✓ Branch 0 taken 477120 times.
✓ Branch 1 taken 480 times.
 477600 do ig=1,ngrid
356 ! Le panache va prendre au debut les caracteristiques de l'air contenu
357 ! dans cette couche.
358
2/2
✓ Branch 0 taken 248157 times.
✓ Branch 1 taken 228963 times.
 477600 if (active(ig)) then
359 248157 ztla(ig,1)=zthl(ig,1)
360 248157 zqta(ig,1)=po(ig,1)
361 248157 zqla(ig,1)=zl(ig,1)
362 !cr: attention, prise en compte de f*(1)=1
363 248157 f_star(ig,2)=alim_star(ig,1)
364 zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
365 & *(zlev(ig,2)-zlev(ig,1)) &
366 248157 & *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
367 248157 w_est(ig,2)=zw2(ig,2)
368 endif
369 enddo
370 !
371
372 !==============================================================================
373 !boucle de calcul de la vitesse verticale dans le thermique
374 !==============================================================================
375
2/2
✓ Branch 0 taken 17760 times.
✓ Branch 1 taken 480 times.
 18240 do l=2,klev-1
376 !==============================================================================
377
378
379 ! On decide si le thermique est encore actif ou non
380 ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
381
2/2
✓ Branch 0 taken 17653440 times.
✓ Branch 1 taken 17760 times.
 17671200 do ig=1,ngrid
382 active(ig)=active(ig) &
383 & .and. zw2(ig,l)>1.e-10 &
384
5/6
✓ Branch 0 taken 1552788 times.
✓ Branch 1 taken 16100652 times.
✓ Branch 2 taken 1304631 times.
✓ Branch 3 taken 248157 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1304631 times.
 34020009 & .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
385 enddo
386
387
388
389 !---------------------------------------------------------------------------
390 ! calcul des proprietes thermodynamiques et de la vitesse de la couche l
391 ! sans tenir compte du detrainement et de l'entrainement dans cette
392 ! couche
393 ! C'est a dire qu'on suppose
394 ! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
395 ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
396 ! avant) a l'alimentation pour avoir un calcul plus propre
397 !---------------------------------------------------------------------------
398
399
2/2
✓ Branch 0 taken 17653440 times.
✓ Branch 1 taken 17760 times.
 17671200 ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
400 17760 call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
401
2/2
✓ Branch 0 taken 17653440 times.
✓ Branch 1 taken 17760 times.
 17671200 do ig=1,ngrid
402 ! print*,'active',active(ig),ig,l
403
2/2
✓ Branch 0 taken 1304631 times.
✓ Branch 1 taken 16348809 times.
 17671200 if(active(ig)) then
404 1304631 zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
405 1304631 ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
406 1304631 zta_est(ig,l)=ztva_est(ig,l)
407 1304631 ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
408 ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1) &
409 1304631 & -zqla_est(ig,l))-zqla_est(ig,l))
410
411
412 !Modif AJAM
413
414 1304631 zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
415 1304631 zdz=zlev(ig,l+1)-zlev(ig,l)
416 1304631 lmel=fact_thermals_ed_dz*zlev(ig,l)
417 ! lmel=0.09*zlev(ig,l)
418 1304631 zlmel=zlev(ig,l)+lmel
419 1304631 zlmelup=zlmel+(zdz/2)
420 1304631 zlmeldwn=zlmel-(zdz/2)
421
422 lt=l+1
423 zlt=zlev(ig,lt)
424 1304631 zdz3=zlev(ig,lt+1)-zlt
425 1304631 zltdwn=zlt-zdz3/2
426 1304631 zltup=zlt+zdz3/2
427
428 !=========================================================================
429 ! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus
430 !=========================================================================
431
432 !--------------------------------------------------
433
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1304631 times.
 1304631 if (iflag_thermals_ed.lt.8) then
434 !--------------------------------------------------
435 !AJ052014: J'ai remplac?? la boucle do par un do while
436 ! afin de faire moins de calcul dans la boucle
437 !--------------------------------------------------
438 ✗ do while (zlmelup.gt.zltup)
439 ✗ lt=lt+1
440 ✗ zlt=zlev(ig,lt)
441 ✗ zdz3=zlev(ig,lt+1)-zlt
442 ✗ zltdwn=zlt-zdz3/2
443 ✗ zltup=zlt+zdz3/2
444 enddo
445 !--------------------------------------------------
446 !AJ052014: Si iflag_thermals_ed<8 (par ex 6), alors
447 ! on cherche o?? se trouve l'altitude d'inversion
448 ! en calculant ztv1 (interpolation de la valeur de
449 ! theta au niveau lt en utilisant les niveaux lt-1 et
450 ! lt-2) et ztv2 (interpolation avec les niveaux lt+1
451 ! et lt+2). Si theta r??ellement calcul??e au niveau lt
452 ! comprise entre ztv1 et ztv2, alors il y a inversion
453 ! et on calcule son altitude zinv en supposant que ztv(lt)
454 ! est une combinaison lineaire de ztv1 et ztv2.
455 ! Ensuite, on calcule la flottabilite en comparant
456 ! la temperature de la couche l a celle de l'air situe
457 ! l+lmel plus haut, ce qui necessite de savoir quel fraction
458 ! de cet air est au-dessus ou en-dessous de l'inversion
459 !--------------------------------------------------
460 ✗ atv1=(ztv(ig,lt-1)-ztv(ig,lt-2))/(zlev(ig,lt-1)-zlev(ig,lt-2))
461 btv1=(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) &
462 ✗ & /(zlev(ig,lt-1)-zlev(ig,lt-2))
463 ✗ atv2=(ztv(ig,lt+2)-ztv(ig,lt+1))/(zlev(ig,lt+2)-zlev(ig,lt+1))
464 btv2=(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) &
465 ✗ & /(zlev(ig,lt+2)-zlev(ig,lt+1))
466
467 ✗ ztv1=atv1*zlt+btv1
468 ✗ ztv2=atv2*zlt+btv2
469
470 ✗ if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
471
472 !--------------------------------------------------
473 !AJ052014: D??calage de zinv qui est entre le haut
474 ! et le bas de la couche lt
475 !--------------------------------------------------
476 ✗ factinv=(ztv2-ztv(ig,lt))/(ztv2-ztv1)
477 ✗ zinv=zltdwn+zdz3*factinv
478
479
480 ✗ if (zlmeldwn.ge.zinv) then
481 ✗ ztv_est(ig,l)=atv2*zlmel+btv2
482 zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) &
483 ✗ & +(1.-fact_shell)*zbuoy(ig,l)
484 ✗ elseif (zlmelup.ge.zinv) then
485 ✗ ztv_est2=atv2*0.5*(zlmelup+zinv)+btv2
486 ✗ ztv_est1=atv1*0.5*(zinv+zlmeldwn)+btv1
487 ✗ ztv_est(ig,l)=((zlmelup-zinv)/zdz)*ztv_est2+((zinv-zlmeldwn)/zdz)*ztv_est1
488
489 zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zinv)/zdz)*(ztva_est(ig,l)- &
490 & ztv_est2)/ztv_est2+((zinv-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
491 ✗ & ztv_est1)/ztv_est1)+(1.-fact_shell)*zbuoy(ig,l)
492
493 else
494 ✗ ztv_est(ig,l)=atv1*zlmel+btv1
495 zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) &
496 ✗ & +(1.-fact_shell)*zbuoy(ig,l)
497 endif
498
499 else ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
500
501 ✗ if (zlmeldwn.gt.zltdwn) then
502 zbuoyjam(ig,l)=fact_shell*RG*((ztva_est(ig,l)- &
503 ✗ & ztv(ig,lt))/ztv(ig,lt))+(1.-fact_shell)*zbuoy(ig,l)
504 else
505 zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- &
506 & ztv(ig,lt))/ztv(ig,lt)+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
507 ✗ & ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
508
509 endif
510
511 ! zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- &
512 ! & ztv1)/ztv1+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
513 ! & ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
514 ! zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- &
515 ! & po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- &
516 ! & po(ig,lt-1))/po(ig,lt-1))
517 endif ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
518
519 else ! if (iflag_thermals_ed.lt.8) then
520 lt=l+1
521 zlt=zlev(ig,lt)
522 zdz2=zlev(ig,lt)-zlev(ig,l)
523
524
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1304631 times.
 1304631 do while (lmel.gt.zdz2)
525 ✗ lt=lt+1
526 ✗ zlt=zlev(ig,lt)
527 ✗ zdz2=zlt-zlev(ig,l)
528 enddo
529 1304631 zdz3=zlev(ig,lt+1)-zlt
530 1304631 zltdwn=zlev(ig,lt)-zdz3/2
531 zlmelup=zlmel+(zdz/2)
532 1304631 coefzlmel=Min(1.,(zlmelup-zltdwn)/zdz)
533 zbuoyjam(ig,l)=1.*RG*(coefzlmel*(ztva_est(ig,l)- &
534 & ztv(ig,lt))/ztv(ig,lt)+(1.-coefzlmel)*(ztva_est(ig,l)- &
535 1304631 & ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
536 endif ! if (iflag_thermals_ed.lt.8) then
537
538 !------------------------------------------------
539 !AJAM:nouveau calcul de w?
540 !------------------------------------------------
541 zdz=zlev(ig,l+1)-zlev(ig,l)
542 zdzbis=zlev(ig,l)-zlev(ig,l-1)
543 zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
544
545 1304631 zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
546 zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
547 1304631 zdw2=afact*zbuoy(ig,l)/fact_epsilon
548 1304631 zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
549 ! zdw2bis=0.5*(zdw2+zdw2bis)
550 lm=Max(1,l-2)
551 ! zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) &
552 ! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
553 ! zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) &
554 ! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
555 ! w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
556 ! w_est(ig,l+1)=(zdz/zdzbis)*Max(0.0001,exp(-zw2fact)* &
557 ! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
558 ! & Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
559 ! w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))*zdw2+w_est(ig,l)*exp(-zw2fact))
560
561 !--------------------------------------------------
562 !AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l)
563 !--------------------------------------------------
564
1/2
✓ Branch 0 taken 1304631 times.
✗ Branch 1 not taken.
 1304631 if (iflag_thermals_ed==8) then
565 ! Ancienne version
566 ! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
567 ! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
568 ! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
569
570 1304631 w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
571
572 ! Nouvelle version Arnaud
573 else
574 ! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
575 ! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
576 ! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
577
578 ✗ w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
579
580 ! w_est(ig,l+1)=Max(0.0001,(zdz/(zdzbis+zdz))*(exp(-zw2fact)* &
581 ! & (w_est(ig,l)-zdw2bis)+zdw2)+(zdzbis/(zdzbis+zdz))* &
582 ! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2bis))
583
584
585
586 ! w_est(ig,l+1)=Max(0.0001,(w_est(ig,l)+zdw2bis*zw2fact)*exp(-zw2fact))
587
588 ! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ &
589 ! & (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis))
590
591 ! w_est(ig,l+1)=Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
592
593 endif
594
595
596
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1304631 times.
 1304631 if (iflag_thermals_ed<6) then
597 ✗ zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
598 ! fact_epsilon=0.0005/(zalpha+0.025)**0.5
599 ! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
600 ✗ fact_epsilon=0.0002/(zalpha+0.1)
601 ✗ zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
602 zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
603 ✗ zdw2=afact*zbuoy(ig,l)/fact_epsilon
604 ✗ zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
605 ! w_est(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
606
607 ! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
608 ! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
609 ! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
610
611 ✗ w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
612
613
614 endif
615 !--------------------------------------------------
616 !AJ052014: J'ai comment? ce if plus n?cessaire puisqu'
617 !on fait max(0.0001,.....)
618 !--------------------------------------------------
619
620 ! if (w_est(ig,l+1).lt.0.) then
621 ! w_est(ig,l+1)=zw2(ig,l)
622 ! w_est(ig,l+1)=0.0001
623 ! endif
624
625 endif
626 enddo
627
628
629 !-------------------------------------------------
630 !calcul des taux d'entrainement et de detrainement
631 !-------------------------------------------------
632
633
2/2
✓ Branch 0 taken 17653440 times.
✓ Branch 1 taken 17760 times.
 17671200 do ig=1,ngrid
634
2/2
✓ Branch 0 taken 1304631 times.
✓ Branch 1 taken 16348809 times.
 17671200 if (active(ig)) then
635
636 ! zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
637 1304631 zw2m=w_est(ig,l+1)
638 ! zw2m=zw2(ig,l)
639 1304631 zdz=zlev(ig,l+1)-zlev(ig,l)
640 1304631 zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
641 ! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
642 zbuoybis=zbuoy(ig,l)
643 1304631 zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
644 1304631 zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
645
646
647 ! entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0., &
648 ! & afact*zbuoybis/zw2m - fact_epsilon )
649
650 ! entr_star(ig,l)=MAX(0.,f_star(ig,l)*zdz*zbetalpha* &
651 ! & afact*zbuoybis/zw2m - fact_epsilon )
652
653
654
655 ! zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
656
657 !=========================================================================
658 ! 4. Calcul de l'entrainement et du detrainement
659 !=========================================================================
660
661 ! entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0., &
662 ! & afact*zbuoyjam(ig,l)/zw2m - fact_epsilon )
663 ! entrbis=entr_star(ig,l)
664
665
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1304631 times.
 1304631 if (iflag_thermals_ed.lt.6) then
666 ✗ fact_epsilon=0.0002/(zalpha+0.1)
667 endif
668
669
670
671 detr_star(ig,l)=f_star(ig,l)*zdz &
672 & *( mix0 * 0.1 / (zalpha+0.001) &
673 & + MAX(detr_min, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m &
674 1304631 & + detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power))
675
676 ! detr_star(ig,l)=(zdz/zdzbis)*detr_star(ig,l)+ &
677 ! & ((zdzbis-zdz)/zdzbis)*detr_star(ig,l-1)
678
679 zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
680
681 entr_star(ig,l)=f_star(ig,l)*zdz* ( &
682 & mix0 * 0.1 / (zalpha+0.001) &
683 & + zbetalpha*MAX(entr_min, &
684 1304631 & afact*zbuoyjam(ig,l)/zw2m - fact_epsilon))
685
686
687 ! entr_star(ig,l)=f_star(ig,l)*zdz* ( &
688 ! & mix0 * 0.1 / (zalpha+0.001) &
689 ! & + MAX(entr_min, &
690 ! & zbetalpha*afact*zbuoyjam(ig,l)/zw2m - fact_epsilon + &
691 ! & detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power))
692
693
694 ! entr_star(ig,l)=(zdz/zdzbis)*entr_star(ig,l)+ &
695 ! & ((zdzbis-zdz)/zdzbis)*entr_star(ig,l-1)
696
697 ! entr_star(ig,l)=Max(0.,f_star(ig,l)*zdz*zbetalpha* &
698 ! & afact*zbuoy(ig,l)/zw2m &
699 ! & - 1.*fact_epsilon)
700
701
702 ! En dessous de lalim, on prend le max de alim_star et entr_star pour
703 ! alim_star et 0 sinon
704
2/2
✓ Branch 0 taken 359864 times.
✓ Branch 1 taken 944767 times.
 1304631 if (l.lt.lalim(ig)) then
705 359864 alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
706 359864 entr_star(ig,l)=0.
707 endif
708 ! if (l.lt.lalim(ig).and.alim_star(ig,l)>alim_star(ig,l-1)) then
709 ! alim_star(ig,l)=entrbis
710 ! endif
711
712 ! print*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l)
713 ! Calcul du flux montant normalise
714 f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
715 1304631 & -detr_star(ig,l)
716
717 endif
718 enddo
719
720
721 !============================================================================
722 ! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique
723 !===========================================================================
724
725
6/6
✓ Branch 0 taken 17653440 times.
✓ Branch 1 taken 17760 times.
✓ Branch 2 taken 1304631 times.
✓ Branch 3 taken 16348809 times.
✓ Branch 4 taken 248157 times.
✓ Branch 5 taken 1056474 times.
 17671200 activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
726
2/2
✓ Branch 0 taken 17653440 times.
✓ Branch 1 taken 17760 times.
 17671200 do ig=1,ngrid
727
2/2
✓ Branch 0 taken 1056474 times.
✓ Branch 1 taken 16596966 times.
 17671200 if (activetmp(ig)) then
728 Zsat=.false.
729 ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
730 & (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
731 1056474 & /(f_star(ig,l+1)+detr_star(ig,l))
732 zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
733 & (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
734 1056474 & /(f_star(ig,l+1)+detr_star(ig,l))
735
736 endif
737 enddo
738
739
2/2
✓ Branch 0 taken 17653440 times.
✓ Branch 1 taken 17760 times.
 17671200 ztemp(:)=zpspsk(:,l)*ztla(:,l)
740 17760 call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
741
2/2
✓ Branch 0 taken 17653440 times.
✓ Branch 1 taken 17760 times.
 17671200 do ig=1,ngrid
742
2/2
✓ Branch 0 taken 1056474 times.
✓ Branch 1 taken 16596966 times.
 17671200 if (activetmp(ig)) then
743 ! on ecrit de maniere conservative (sat ou non)
744 ! T = Tl +Lv/Cp ql
745 1056474 zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
746 1056474 ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
747 1056474 ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
748 !on rajoute le calcul de zha pour diagnostiques (temp potentielle)
749 1056474 zha(ig,l) = ztva(ig,l)
750 ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l) &
751 1056474 & -zqla(ig,l))-zqla(ig,l))
752 1056474 zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
753 1056474 zdz=zlev(ig,l+1)-zlev(ig,l)
754 1056474 zdzbis=zlev(ig,l)-zlev(ig,l-1)
755 1056474 zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
756 !!!!!!! fact_epsilon=0.002
757 1056474 zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
758 zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
759 1056474 zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
760 1056474 zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
761 ! zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) &
762 ! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
763 ! lm=Max(1,l-2)
764 ! zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) &
765 ! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
766 ! zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
767 ! zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ &
768 ! & (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis))
769 ! zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
770 ! zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
771 ! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
772 ! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
773
1/2
✓ Branch 0 taken 1056474 times.
✗ Branch 1 not taken.
 1056474 if (iflag_thermals_ed==8) then
774 1056474 zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
775 else
776 ✗ zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
777 endif
778 ! zw2(ig,l+1)=Max(0.0001,(zdz/(zdz+zdzbis))*(exp(-zw2fact)* &
779 ! & (zw2(ig,l)-zdw2)+zdw2bis)+(zdzbis/(zdz+zdzbis))* &
780 ! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2bis))
781
782
783
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1056474 times.
 1056474 if (iflag_thermals_ed.lt.6) then
784 ✗ zalpha=f0(ig)*f_star(ig,l)/sqrt(zw2(ig,l+1))/rhobarz(ig,l)
785 ! fact_epsilon=0.0005/(zalpha+0.025)**0.5
786 ! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
787 ✗ fact_epsilon=0.0002/(zalpha+0.1)**1
788 ✗ zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
789 zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
790 ✗ zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
791 ✗ zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
792
793 ! zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
794 ! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
795 ! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
796 ! zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
797 ✗ zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
798
799 endif
800
801
802 endif
803 enddo
804
805
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 17760 times.
 17760 if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
806 !
807 !===========================================================================
808 ! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
809 !===========================================================================
810
811 17760 nbpb=0
812
2/2
✓ Branch 0 taken 17653440 times.
✓ Branch 1 taken 17760 times.
 17671200 do ig=1,ngrid
813
3/4
✓ Branch 0 taken 1056474 times.
✓ Branch 1 taken 16596966 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1056474 times.
 17653440 if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
814 ! stop'On tombe sur le cas particulier de thermcell_dry'
815 ! print*,'On tombe sur le cas particulier de thermcell_plume'
816 ✗ nbpb=nbpb+1
817 ✗ zw2(ig,l+1)=0.
818 ✗ linter(ig)=l+1
819 endif
820
821
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 17653440 times.
 17653440 if (zw2(ig,l+1).lt.0.) then
822 linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
823 ✗ & -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
824 ✗ zw2(ig,l+1)=0.
825 !+CR:04/05/12:correction calcul linter pour calcul de zmax continu
826
2/2
✓ Branch 0 taken 248157 times.
✓ Branch 1 taken 17405283 times.
 17653440 elseif (f_star(ig,l+1).lt.0.) then
827 linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) &
828 248157 & -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
829 248157 zw2(ig,l+1)=0.
830 !fin CR:04/05/12
831 endif
832
833 17653440 wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
834
835
2/2
✓ Branch 0 taken 1030476 times.
✓ Branch 1 taken 16622964 times.
 17671200 if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
836 ! lmix est le niveau de la couche ou w (wa_moy) est maximum
837 !on rajoute le calcul de lmix_bis
838
2/2
✓ Branch 0 taken 655258 times.
✓ Branch 1 taken 375218 times.
 1030476 if (zqla(ig,l).lt.1.e-10) then
839 655258 lmix_bis(ig)=l+1
840 endif
841 1030476 lmix(ig)=l+1
842 1030476 wmaxa(ig)=wa_moy(ig,l+1)
843 endif
844 enddo
845
846
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 17760 times.
 18240 if (nbpb>0) then
847 ✗ print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
848 endif
849
850 !=========================================================================
851 ! FIN DE LA BOUCLE VERTICALE
852 enddo
853 !=========================================================================
854
855 !on recalcule alim_star_tot
856
2/2
✓ Branch 0 taken 477120 times.
✓ Branch 1 taken 480 times.
 477600 do ig=1,ngrid
857 477600 alim_star_tot(ig)=0.
858 enddo
859
2/2
✓ Branch 0 taken 477120 times.
✓ Branch 1 taken 480 times.
 477600 do ig=1,ngrid
860
2/2
✓ Branch 0 taken 609661 times.
✓ Branch 1 taken 477120 times.
 1087261 do l=1,lalim(ig)-1
861 1086781 alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
862 enddo
863 enddo
864
865
866
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 480 times.
 480 if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
867
868
869
870 480 return
871 end
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
908 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
909 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
910 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
911 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
912 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
913 ✗ SUBROUTINE thermcell_plume_5B(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, &
914 ✗ & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
915 ✗ & lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
916 ✗ & ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
917 & ,lev_out,lunout1,igout)
918 !& ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
919
920 !--------------------------------------------------------------------------
921 !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
922 ! Version conforme a l'article de Rio et al. 2010.
923 ! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
924 !--------------------------------------------------------------------------
925
926 USE print_control_mod, ONLY: prt_level
927 IMPLICIT NONE
928
929 !
930 ! $Header$
931 !
932 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
933 ! veillez � n'utiliser que des ! pour les commentaires
934 ! et � bien positionner les & des lignes de continuation
935 ! (les placer en colonne 6 et en colonne 73)
936 !
937 !
938 ! A1.0 Fundamental constants
939 REAL RPI,RCLUM,RHPLA,RKBOL,RNAVO
940 ! A1.1 Astronomical constants
941 REAL RDAY,REA,REPSM,RSIYEA,RSIDAY,ROMEGA
942 ! A1.1.bis Constantes concernant l'orbite de la Terre:
943 REAL R_ecc, R_peri, R_incl
944 ! A1.2 Geoide
945 REAL RA,RG,R1SA
946 ! A1.3 Radiation
947 ! REAL RSIGMA,RI0
948 REAL RSIGMA
949 ! A1.4 Thermodynamic gas phase
950 REAL RMO3,RMCO2,RMC,RMCH4,RMN2O,RMCFC11,RMCFC12
951 REAL R,RMD,RMV,RD,RV,RCPD,RCPV,RCVD,RCVV
952 REAL RKAPPA,RETV, eps_w
953 ! A1.5,6 Thermodynamic liquid,solid phases
954 REAL RCW,RCS
955 ! A1.7 Thermodynamic transition of phase
956 REAL RLVTT,RLSTT,RLMLT,RTT,RATM
957 ! A1.8 Curve of saturation
958 REAL RESTT,RALPW,RBETW,RGAMW,RALPS,RBETS,RGAMS
959 REAL RALPD,RBETD,RGAMD
960 !
961 COMMON/YOMCST/RPI ,RCLUM ,RHPLA ,RKBOL ,RNAVO &
962 & ,RDAY ,REA ,REPSM ,RSIYEA,RSIDAY,ROMEGA &
963 & ,R_ecc, R_peri, R_incl &
964 & ,RA ,RG ,R1SA &
965 & ,RSIGMA &
966 & ,R ,RMD ,RMV ,RD ,RV ,RCPD &
967 & ,RMO3 ,RMCO2 ,RMC ,RMCH4 ,RMN2O ,RMCFC11 ,RMCFC12 &
968 & ,RCPV ,RCVD ,RCVV ,RKAPPA,RETV, eps_w &
969 & ,RCW ,RCS &
970 & ,RLVTT ,RLSTT ,RLMLT ,RTT ,RATM &
971 & ,RESTT ,RALPW ,RBETW ,RGAMW ,RALPS ,RBETS ,RGAMS &
972 & ,RALPD ,RBETD ,RGAMD
973 ! ------------------------------------------------------------------
974 !$OMP THREADPRIVATE(/YOMCST/)
975 !
976 ! $Id: YOETHF.h 2799 2017-02-24 18:50:33Z jyg $
977 !
978 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
979 ! veillez n'utiliser que des ! pour les commentaires
980 ! et bien positionner les & des lignes de continuation
981 ! (les placer en colonne 6 et en colonne 73)
982 !
983 !* COMMON *YOETHF* DERIVED CONSTANTS SPECIFIC TO ECMWF THERMODYNAMICS
984 !
985 ! *R__ES* *CONSTANTS USED FOR COMPUTATION OF SATURATION
986 ! MIXING RATIO OVER LIQUID WATER(*R_LES*) OR
987 ! ICE(*R_IES*).
988 ! *RVTMP2* *RVTMP2=RCPV/RCPD-1.
989 ! *RHOH2O* *DENSITY OF LIQUID WATER. (RATM/100.)
990 !
991 REAL R2ES, R3LES, R3IES, R4LES, R4IES, R5LES, R5IES
992 REAL RVTMP2, RHOH2O
993 REAL R5ALVCP,R5ALSCP,RALVDCP,RALSDCP,RALFDCP,RTWAT,RTBER,RTBERCU
994 REAL RTICE,RTICECU,RTWAT_RTICE_R,RTWAT_RTICECU_R,RKOOP1,RKOOP2
995 LOGICAL OK_BAD_ECMWF_THERMO ! If TRUE, then variables set by rrtm/suphec.F90
996 ! If FALSE, then variables set by suphel.F90
997 COMMON /YOETHF/R2ES, R3LES, R3IES, R4LES, R4IES, R5LES, R5IES, &
998 & RVTMP2, RHOH2O, &
999 & R5ALVCP,R5ALSCP,RALVDCP,RALSDCP, &
1000 & RALFDCP,RTWAT,RTBER,RTBERCU, &
1001 & RTICE,RTICECU,RTWAT_RTICE_R,RTWAT_RTICECU_R,RKOOP1,&
1002 & RKOOP2, &
1003 & OK_BAD_ECMWF_THERMO
1004
1005 !$OMP THREADPRIVATE(/YOETHF/)
1006 !
1007 ! $Header$
1008 !
1009 !
1010 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
1011 ! veillez n'utiliser que des ! pour les commentaires
1012 ! et bien positionner les & des lignes de continuation
1013 ! (les placer en colonne 6 et en colonne 73)
1014 !
1015 ! ------------------------------------------------------------------
1016 ! This COMDECK includes the Thermodynamical functions for the cy39
1017 ! ECMWF Physics package.
1018 ! Consistent with YOMCST Basic physics constants, assuming the
1019 ! partial pressure of water vapour is given by a first order
1020 ! Taylor expansion of Qs(T) w.r.t. to Temperature, using constants
1021 ! in YOETHF
1022 ! ------------------------------------------------------------------
1023 REAL PTARG, PDELARG, P5ARG, PQSARG, PCOARG
1024 REAL FOEEW, FOEDE, qsats, qsatl, dqsats, dqsatl
1025 LOGICAL thermcep
1026 PARAMETER (thermcep=.TRUE.)
1027 !
1028 FOEEW ( PTARG,PDELARG ) = EXP ( &
1029 & (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT) &
1030 & / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) )
1031 !
1032 FOEDE ( PTARG,PDELARG,P5ARG,PQSARG,PCOARG ) = PQSARG*PCOARG*P5ARG &
1033 & / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG))**2
1034 !
1035 qsats(ptarg) = 100.0 * 0.622 * 10.0 &
1036 & ** (2.07023 - 0.00320991 * ptarg &
1037 & - 2484.896 / ptarg + 3.56654 * LOG10(ptarg))
1038 qsatl(ptarg) = 100.0 * 0.622 * 10.0 &
1039 & ** (23.8319 - 2948.964 / ptarg &
1040 & - 5.028 * LOG10(ptarg) &
1041 & - 29810.16 * EXP( - 0.0699382 * ptarg) &
1042 & + 25.21935 * EXP( - 2999.924 / ptarg))
1043 !
1044 dqsats(ptarg,pqsarg) = RLVTT/RCPD*pqsarg * (3.56654/ptarg &
1045 & +2484.896*LOG(10.)/ptarg**2 &
1046 & -0.00320991*LOG(10.))
1047 dqsatl(ptarg,pqsarg) = RLVTT/RCPD*pqsarg*LOG(10.)* &
1048 & (2948.964/ptarg**2-5.028/LOG(10.)/ptarg &
1049 & +25.21935*2999.924/ptarg**2*EXP(-2999.924/ptarg) &
1050 & +29810.16*0.0699382*EXP(-0.0699382*ptarg))
1051 integer :: iflag_thermals,nsplit_thermals
1052
1053 !!! nrlmd le 10/04/2012
1054 integer :: iflag_trig_bl,iflag_clos_bl
1055 integer :: tau_trig_shallow,tau_trig_deep
1056 real :: s_trig
1057 !!! fin nrlmd le 10/04/2012
1058
1059 real,parameter :: r_aspect_thermals=2.,l_mix_thermals=30.
1060 real :: alp_bl_k
1061 real :: tau_thermals,fact_thermals_ed_dz
1062 integer,parameter :: w2di_thermals=0
1063 integer :: isplit
1064
1065 integer :: iflag_coupl,iflag_clos,iflag_wake
1066 integer :: iflag_thermals_ed,iflag_thermals_optflux,iflag_thermals_closure
1067
1068 common/ctherm1/iflag_thermals,nsplit_thermals,iflag_thermals_closure
1069 common/ctherm2/tau_thermals,alp_bl_k,fact_thermals_ed_dz
1070 common/ctherm4/iflag_coupl,iflag_clos,iflag_wake
1071 common/ctherm5/iflag_thermals_ed,iflag_thermals_optflux
1072
1073 !!! nrlmd le 10/04/2012
1074 common/ctherm6/iflag_trig_bl,iflag_clos_bl
1075 common/ctherm7/tau_trig_shallow,tau_trig_deep
1076 common/ctherm8/s_trig
1077 !!! fin nrlmd le 10/04/2012
1078
1079 !$OMP THREADPRIVATE(/ctherm1/,/ctherm2/,/ctherm4/,/ctherm5/)
1080 !$OMP THREADPRIVATE(/ctherm6/,/ctherm7/,/ctherm8/)
1081
1082 INTEGER itap
1083 INTEGER lunout1,igout
1084 INTEGER ngrid,klev
1085 REAL ptimestep
1086 REAL ztv(ngrid,klev)
1087 REAL zthl(ngrid,klev)
1088 REAL po(ngrid,klev)
1089 REAL zl(ngrid,klev)
1090 REAL rhobarz(ngrid,klev)
1091 REAL zlev(ngrid,klev+1)
1092 REAL pplev(ngrid,klev+1)
1093 REAL pphi(ngrid,klev)
1094 REAL zpspsk(ngrid,klev)
1095 REAL alim_star(ngrid,klev)
1096 REAL f0(ngrid)
1097 INTEGER lalim(ngrid)
1098 integer lev_out ! niveau pour les print
1099 integer nbpb
1100
1101 real alim_star_tot(ngrid)
1102
1103 REAL ztva(ngrid,klev)
1104 REAL ztla(ngrid,klev)
1105 REAL zqla(ngrid,klev)
1106 REAL zqta(ngrid,klev)
1107 REAL zha(ngrid,klev)
1108
1109 REAL detr_star(ngrid,klev)
1110 REAL coefc
1111 REAL entr_star(ngrid,klev)
1112 ✗ REAL detr(ngrid,klev)
1113 ✗ REAL entr(ngrid,klev)
1114
1115 REAL csc(ngrid,klev)
1116
1117 REAL zw2(ngrid,klev+1)
1118 REAL w_est(ngrid,klev+1)
1119 REAL f_star(ngrid,klev+1)
1120 ✗ REAL wa_moy(ngrid,klev+1)
1121
1122 REAL ztva_est(ngrid,klev)
1123 ✗ REAL zqla_est(ngrid,klev)
1124 REAL zqsatth(ngrid,klev)
1125 ✗ REAL zta_est(ngrid,klev)
1126 ✗ REAL zbuoyjam(ngrid,klev)
1127 ✗ REAL ztemp(ngrid),zqsat(ngrid)
1128 REAL zdw2
1129 REAL zw2modif
1130 REAL zw2fact
1131 ✗ REAL zeps(ngrid,klev)
1132
1133 REAL linter(ngrid)
1134 INTEGER lmix(ngrid)
1135 INTEGER lmix_bis(ngrid)
1136 ✗ REAL wmaxa(ngrid)
1137
1138 INTEGER ig,l,k
1139
1140 ✗ real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
1141 real zbuoybis
1142 real zcor,zdelta,zcvm5,qlbef,zdz2
1143 real betalpha,zbetalpha
1144 real eps, afact
1145 REAL REPS,RLvCp,DDT0
1146 PARAMETER (DDT0=.01)
1147 logical Zsat
1148 ✗ LOGICAL active(ngrid),activetmp(ngrid)
1149 REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
1150 REAL c2(ngrid,klev)
1151 Zsat=.false.
1152 ! Initialisation
1153
1154 ✗ RLvCp = RLVTT/RCPD
1155 fact_epsilon=0.002
1156 betalpha=0.9
1157 afact=2./3.
1158
1159 zbetalpha=betalpha/(1.+betalpha)
1160
1161
1162 ! Initialisations des variables reeles
1163 if (1==1) then
1164 ✗ ztva(:,:)=ztv(:,:)
1165 ✗ ztva_est(:,:)=ztva(:,:)
1166 ✗ ztla(:,:)=zthl(:,:)
1167 ✗ zqta(:,:)=po(:,:)
1168 ✗ zha(:,:) = ztva(:,:)
1169 else
1170 ztva(:,:)=0.
1171 ztva_est(:,:)=0.
1172 ztla(:,:)=0.
1173 zqta(:,:)=0.
1174 zha(:,:) =0.
1175 endif
1176
1177 ✗ zqla_est(:,:)=0.
1178 ✗ zqsatth(:,:)=0.
1179 ✗ zqla(:,:)=0.
1180 ✗ detr_star(:,:)=0.
1181 ✗ entr_star(:,:)=0.
1182 ✗ alim_star(:,:)=0.
1183 ✗ alim_star_tot(:)=0.
1184 ✗ csc(:,:)=0.
1185 ✗ detr(:,:)=0.
1186 ✗ entr(:,:)=0.
1187 ✗ zw2(:,:)=0.
1188 ✗ zbuoy(:,:)=0.
1189 ✗ zbuoyjam(:,:)=0.
1190 ✗ gamma(:,:)=0.
1191 ✗ zeps(:,:)=0.
1192 ✗ w_est(:,:)=0.
1193 ✗ f_star(:,:)=0.
1194 ✗ wa_moy(:,:)=0.
1195 ✗ linter(:)=1.
1196 ! linter(:)=1.
1197 ! Initialisation des variables entieres
1198 ✗ lmix(:)=1
1199 ✗ lmix_bis(:)=2
1200 ✗ wmaxa(:)=0.
1201 ✗ lalim(:)=1
1202
1203
1204 !-------------------------------------------------------------------------
1205 ! On ne considere comme actif que les colonnes dont les deux premieres
1206 ! couches sont instables.
1207 !-------------------------------------------------------------------------
1208 ✗ active(:)=ztv(:,1)>ztv(:,2)
1209
1210 !-------------------------------------------------------------------------
1211 ! Definition de l'alimentation a l'origine dans thermcell_init
1212 !-------------------------------------------------------------------------
1213 ✗ do l=1,klev-1
1214 ✗ do ig=1,ngrid
1215 ✗ if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
1216 alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
1217 ✗ & *sqrt(zlev(ig,l+1))
1218 ✗ lalim(ig)=l+1
1219 ✗ alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
1220 endif
1221 enddo
1222 enddo
1223 ✗ do l=1,klev
1224 ✗ do ig=1,ngrid
1225 ✗ if (alim_star_tot(ig) > 1.e-10 ) then
1226 ✗ alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
1227 endif
1228 enddo
1229 enddo
1230 ✗ alim_star_tot(:)=1.
1231
1232
1233
1234 !------------------------------------------------------------------------------
1235 ! Calcul dans la premiere couche
1236 ! On decide dans cette version que le thermique n'est actif que si la premiere
1237 ! couche est instable.
1238 ! Pourrait etre change si on veut que le thermiques puisse se d??clencher
1239 ! dans une couche l>1
1240 !------------------------------------------------------------------------------
1241 ✗ do ig=1,ngrid
1242 ! Le panache va prendre au debut les caracteristiques de l'air contenu
1243 ! dans cette couche.
1244 ✗ if (active(ig)) then
1245 ✗ ztla(ig,1)=zthl(ig,1)
1246 ✗ zqta(ig,1)=po(ig,1)
1247 ✗ zqla(ig,1)=zl(ig,1)
1248 !cr: attention, prise en compte de f*(1)=1
1249 ✗ f_star(ig,2)=alim_star(ig,1)
1250 zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
1251 & *(zlev(ig,2)-zlev(ig,1)) &
1252 ✗ & *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
1253 ✗ w_est(ig,2)=zw2(ig,2)
1254 endif
1255 enddo
1256 !
1257
1258 !==============================================================================
1259 !boucle de calcul de la vitesse verticale dans le thermique
1260 !==============================================================================
1261 ✗ do l=2,klev-1
1262 !==============================================================================
1263
1264
1265 ! On decide si le thermique est encore actif ou non
1266 ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
1267 ✗ do ig=1,ngrid
1268 active(ig)=active(ig) &
1269 & .and. zw2(ig,l)>1.e-10 &
1270 ✗ & .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
1271 enddo
1272
1273
1274
1275 !---------------------------------------------------------------------------
1276 ! calcul des proprietes thermodynamiques et de la vitesse de la couche l
1277 ! sans tenir compte du detrainement et de l'entrainement dans cette
1278 ! couche
1279 ! C'est a dire qu'on suppose
1280 ! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
1281 ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
1282 ! avant) a l'alimentation pour avoir un calcul plus propre
1283 !---------------------------------------------------------------------------
1284
1285 ✗ ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
1286 ✗ call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
1287
1288 ✗ do ig=1,ngrid
1289 ! print*,'active',active(ig),ig,l
1290 ✗ if(active(ig)) then
1291 ✗ zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
1292 ✗ ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
1293 ✗ zta_est(ig,l)=ztva_est(ig,l)
1294 ✗ ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
1295 ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1) &
1296 ✗ & -zqla_est(ig,l))-zqla_est(ig,l))
1297
1298 !------------------------------------------------
1299 !AJAM:nouveau calcul de w?
1300 !------------------------------------------------
1301 ✗ zdz=zlev(ig,l+1)-zlev(ig,l)
1302 ✗ zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
1303
1304 ✗ zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
1305 ✗ zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
1306 ✗ w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
1307
1308
1309 ✗ if (w_est(ig,l+1).lt.0.) then
1310 ✗ w_est(ig,l+1)=zw2(ig,l)
1311 endif
1312 endif
1313 enddo
1314
1315
1316 !-------------------------------------------------
1317 !calcul des taux d'entrainement et de detrainement
1318 !-------------------------------------------------
1319
1320 ✗ do ig=1,ngrid
1321 ✗ if (active(ig)) then
1322
1323 ✗ zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
1324 zw2m=w_est(ig,l+1)
1325 ✗ zdz=zlev(ig,l+1)-zlev(ig,l)
1326 ✗ zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
1327 ! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
1328 zbuoybis=zbuoy(ig,l)
1329 ✗ zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
1330 ✗ zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
1331
1332
1333 entr_star(ig,l)=f_star(ig,l)*zdz* zbetalpha*MAX(0., &
1334 ✗ & afact*zbuoybis/zw2m - fact_epsilon )
1335
1336
1337 detr_star(ig,l)=f_star(ig,l)*zdz &
1338 & *MAX(1.e-3, -afact*zbetalpha*zbuoy(ig,l)/zw2m &
1339 ✗ & + 0.012*(zdqt(ig,l)/zw2m)**0.5 )
1340
1341 ! En dessous de lalim, on prend le max de alim_star et entr_star pour
1342 ! alim_star et 0 sinon
1343 ✗ if (l.lt.lalim(ig)) then
1344 ✗ alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
1345 ✗ entr_star(ig,l)=0.
1346 endif
1347
1348 ! Calcul du flux montant normalise
1349 f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
1350 ✗ & -detr_star(ig,l)
1351
1352 endif
1353 enddo
1354
1355
1356 !----------------------------------------------------------------------------
1357 !calcul de la vitesse verticale en melangeant Tl et qt du thermique
1358 !---------------------------------------------------------------------------
1359 ✗ activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
1360 ✗ do ig=1,ngrid
1361 ✗ if (activetmp(ig)) then
1362 Zsat=.false.
1363 ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
1364 & (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
1365 ✗ & /(f_star(ig,l+1)+detr_star(ig,l))
1366 zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
1367 & (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
1368 ✗ & /(f_star(ig,l+1)+detr_star(ig,l))
1369
1370 endif
1371 enddo
1372
1373 ✗ ztemp(:)=zpspsk(:,l)*ztla(:,l)
1374 ✗ call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
1375
1376 ✗ do ig=1,ngrid
1377 ✗ if (activetmp(ig)) then
1378 ! on ecrit de maniere conservative (sat ou non)
1379 ! T = Tl +Lv/Cp ql
1380 ✗ zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
1381 ✗ ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
1382 ✗ ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
1383 !on rajoute le calcul de zha pour diagnostiques (temp potentielle)
1384 ✗ zha(ig,l) = ztva(ig,l)
1385 ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l) &
1386 ✗ & -zqla(ig,l))-zqla(ig,l))
1387 ✗ zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
1388 ✗ zdz=zlev(ig,l+1)-zlev(ig,l)
1389 ✗ zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
1390
1391 ✗ zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
1392 ✗ zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
1393 ✗ zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
1394 endif
1395 enddo
1396
1397 ✗ if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
1398 !
1399 !---------------------------------------------------------------------------
1400 !initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
1401 !---------------------------------------------------------------------------
1402
1403 ✗ nbpb=0
1404 ✗ do ig=1,ngrid
1405 ✗ if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
1406 ! stop'On tombe sur le cas particulier de thermcell_dry'
1407 ! print*,'On tombe sur le cas particulier de thermcell_plume'
1408 ✗ nbpb=nbpb+1
1409 ✗ zw2(ig,l+1)=0.
1410 ✗ linter(ig)=l+1
1411 endif
1412
1413 ✗ if (zw2(ig,l+1).lt.0.) then
1414 linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
1415 ✗ & -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
1416 ✗ zw2(ig,l+1)=0.
1417 ✗ elseif (f_star(ig,l+1).lt.0.) then
1418 linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) &
1419 ✗ & -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
1420 ! print*,"linter plume", linter(ig)
1421 ✗ zw2(ig,l+1)=0.
1422 endif
1423
1424 ✗ wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
1425
1426 ✗ if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
1427 ! lmix est le niveau de la couche ou w (wa_moy) est maximum
1428 !on rajoute le calcul de lmix_bis
1429 ✗ if (zqla(ig,l).lt.1.e-10) then
1430 ✗ lmix_bis(ig)=l+1
1431 endif
1432 ✗ lmix(ig)=l+1
1433 ✗ wmaxa(ig)=wa_moy(ig,l+1)
1434 endif
1435 enddo
1436
1437 ✗ if (nbpb>0) then
1438 ✗ print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
1439 endif
1440
1441 !=========================================================================
1442 ! FIN DE LA BOUCLE VERTICALE
1443 enddo
1444 !=========================================================================
1445
1446 !on recalcule alim_star_tot
1447 ✗ do ig=1,ngrid
1448 ✗ alim_star_tot(ig)=0.
1449 enddo
1450 ✗ do ig=1,ngrid
1451 ✗ do l=1,lalim(ig)-1
1452 ✗ alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
1453 enddo
1454 enddo
1455
1456
1457 ✗ if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
1458
1459
1460
1461 ✗ return
1462 end
1463