GCC Code Coverage Report


Directory: ./
File: phys/conema3.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 137 0.0%
Branches: 0 156 0.0%

Line Branch Exec Source
1
2 ! $Id: conema3.F90 2346 2015-08-21 15:13:46Z emillour $
3
4 SUBROUTINE conema3(dtime, paprs, pplay, t, q, u, v, tra, ntra, work1, work2, &
5 d_t, d_q, d_u, d_v, d_tra, rain, snow, kbas, ktop, upwd, dnwd, dnwdbis, &
6 bas, top, ma, cape, tvp, rflag, pbase, bbase, dtvpdt1, dtvpdq1, dplcldt, &
7 dplcldr, qcond_incld)
8
9 USE dimphy
10 USE infotrac_phy, ONLY: nbtr
11 IMPLICIT NONE
12 ! ======================================================================
13 ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
14 ! Objet: schema de convection de Emanuel (1991) interface
15 ! Mai 1998: Interface modifiee pour implementation dans LMDZ
16 ! ======================================================================
17 ! Arguments:
18 ! dtime---input-R-pas d'integration (s)
19 ! paprs---input-R-pression inter-couches (Pa)
20 ! pplay---input-R-pression au milieu des couches (Pa)
21 ! t-------input-R-temperature (K)
22 ! q-------input-R-humidite specifique (kg/kg)
23 ! u-------input-R-vitesse du vent zonal (m/s)
24 ! v-------input-R-vitesse duvent meridien (m/s)
25 ! tra-----input-R-tableau de rapport de melange des traceurs
26 ! work*: input et output: deux variables de travail,
27 ! on peut les mettre a 0 au debut
28
29 ! d_t-----output-R-increment de la temperature
30 ! d_q-----output-R-increment de la vapeur d'eau
31 ! d_u-----output-R-increment de la vitesse zonale
32 ! d_v-----output-R-increment de la vitesse meridienne
33 ! d_tra---output-R-increment du contenu en traceurs
34 ! rain----output-R-la pluie (mm/s)
35 ! snow----output-R-la neige (mm/s)
36 ! kbas----output-R-bas du nuage (integer)
37 ! ktop----output-R-haut du nuage (integer)
38 ! upwd----output-R-saturated updraft mass flux (kg/m**2/s)
39 ! dnwd----output-R-saturated downdraft mass flux (kg/m**2/s)
40 ! dnwdbis-output-R-unsaturated downdraft mass flux (kg/m**2/s)
41 ! bas-----output-R-bas du nuage (real)
42 ! top-----output-R-haut du nuage (real)
43 ! Ma------output-R-flux ascendant non dilue (kg/m**2/s)
44 ! cape----output-R-CAPE
45 ! tvp-----output-R-virtual temperature of the lifted parcel
46 ! rflag---output-R-flag sur le fonctionnement de convect
47 ! pbase---output-R-pression a la base du nuage (Pa)
48 ! bbase---output-R-buoyancy a la base du nuage (K)
49 ! dtvpdt1-output-R-derivative of parcel virtual temp wrt T1
50 ! dtvpdq1-output-R-derivative of parcel virtual temp wrt Q1
51 ! dplcldt-output-R-derivative of the PCP pressure wrt T1
52 ! dplcldr-output-R-derivative of the PCP pressure wrt Q1
53 ! ======================================================================
54
55 include "conema3.h"
56 INTEGER i, l, m, itra
57 INTEGER ntra ! if no tracer transport
58 ! is needed, set ntra = 1 (or 0)
59 REAL dtime
60
61 REAL d_t2(klon, klev), d_q2(klon, klev) ! sbl
62 REAL d_u2(klon, klev), d_v2(klon, klev) ! sbl
63 REAL em_d_t2(klev), em_d_q2(klev) ! sbl
64 REAL em_d_u2(klev), em_d_v2(klev) ! sbl
65
66 REAL paprs(klon, klev+1), pplay(klon, klev)
67 REAL t(klon, klev), q(klon, klev), d_t(klon, klev), d_q(klon, klev)
68 REAL u(klon, klev), v(klon, klev), tra(klon, klev, ntra)
69 REAL d_u(klon, klev), d_v(klon, klev), d_tra(klon, klev, ntra)
70 REAL work1(klon, klev), work2(klon, klev)
71 REAL upwd(klon, klev), dnwd(klon, klev), dnwdbis(klon, klev)
72 REAL rain(klon)
73 REAL snow(klon)
74 REAL cape(klon), tvp(klon, klev), rflag(klon)
75 REAL pbase(klon), bbase(klon)
76 REAL dtvpdt1(klon, klev), dtvpdq1(klon, klev)
77 REAL dplcldt(klon), dplcldr(klon)
78 INTEGER kbas(klon), ktop(klon)
79
80 REAL wd(klon)
81 REAL qcond_incld(klon, klev)
82
83 LOGICAL, SAVE :: first = .TRUE.
84 !$OMP THREADPRIVATE(first)
85
86 ! ym REAL em_t(klev)
87 REAL, ALLOCATABLE, SAVE :: em_t(:)
88 !$OMP THREADPRIVATE(em_t)
89 ! ym REAL em_q(klev)
90 REAL, ALLOCATABLE, SAVE :: em_q(:)
91 !$OMP THREADPRIVATE(em_q)
92 ! ym REAL em_qs(klev)
93 REAL, ALLOCATABLE, SAVE :: em_qs(:)
94 !$OMP THREADPRIVATE(em_qs)
95 ! ym REAL em_u(klev), em_v(klev), em_tra(klev,nbtr)
96 REAL, ALLOCATABLE, SAVE :: em_u(:), em_v(:), em_tra(:, :)
97 !$OMP THREADPRIVATE(em_u,em_v,em_tra)
98 ! ym REAL em_ph(klev+1), em_p(klev)
99 REAL, ALLOCATABLE, SAVE :: em_ph(:), em_p(:)
100 !$OMP THREADPRIVATE(em_ph,em_p)
101 ! ym REAL em_work1(klev), em_work2(klev)
102 REAL, ALLOCATABLE, SAVE :: em_work1(:), em_work2(:)
103 !$OMP THREADPRIVATE(em_work1,em_work2)
104 ! ym REAL em_precip, em_d_t(klev), em_d_q(klev)
105 REAL, SAVE :: em_precip
106 !$OMP THREADPRIVATE(em_precip)
107 REAL, ALLOCATABLE, SAVE :: em_d_t(:), em_d_q(:)
108 !$OMP THREADPRIVATE(em_d_t,em_d_q)
109 ! ym REAL em_d_u(klev), em_d_v(klev), em_d_tra(klev,nbtr)
110 REAL, ALLOCATABLE, SAVE :: em_d_u(:), em_d_v(:), em_d_tra(:, :)
111 !$OMP THREADPRIVATE(em_d_u,em_d_v,em_d_tra)
112 ! ym REAL em_upwd(klev), em_dnwd(klev), em_dnwdbis(klev)
113 REAL, ALLOCATABLE, SAVE :: em_upwd(:), em_dnwd(:), em_dnwdbis(:)
114 !$OMP THREADPRIVATE(em_upwd,em_dnwd,em_dnwdbis)
115 REAL em_dtvpdt1(klev), em_dtvpdq1(klev)
116 REAL em_dplcldt, em_dplcldr
117 ! ym SAVE em_t,em_q, em_qs, em_ph, em_p, em_work1, em_work2
118 ! ym SAVE em_u,em_v, em_tra
119 ! ym SAVE em_d_u,em_d_v, em_d_tra
120 ! ym SAVE em_precip, em_d_t, em_d_q, em_upwd, em_dnwd, em_dnwdbis
121
122 INTEGER em_bas, em_top
123 SAVE em_bas, em_top
124 !$OMP THREADPRIVATE(em_bas,em_top)
125 REAL em_wd
126 REAL em_qcond(klev)
127 REAL em_qcondc(klev)
128
129 REAL zx_t, zx_qs, zdelta, zcor
130 INTEGER iflag
131 REAL sigsum
132 ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
133 ! VARIABLES A SORTIR
134 ! ccccccccccccccccccccccccccccccccccccccccccccccccc
135
136 ! ym REAL emmip(klev) !variation de flux ascnon dilue i et i+1
137 REAL, ALLOCATABLE, SAVE :: emmip(:)
138 !$OMP THREADPRIVATE(emmip)
139 ! ym SAVE emmip
140 ! ym real emMke(klev)
141 REAL, ALLOCATABLE, SAVE :: emmke(:)
142 !$OMP THREADPRIVATE(emMke)
143 ! ym save emMke
144 REAL top
145 REAL bas
146 ! ym real emMa(klev)
147 REAL, ALLOCATABLE, SAVE :: emma(:)
148 !$OMP THREADPRIVATE(emMa)
149 ! ym save emMa
150 REAL ma(klon, klev)
151 REAL ment(klev, klev)
152 REAL qent(klev, klev)
153 REAL tps(klev), tls(klev)
154 REAL sij(klev, klev)
155 REAL em_cape, em_tvp(klev)
156 REAL em_pbase, em_bbase
157 INTEGER iw, j, k, ix, iy
158
159 ! -- sb: pour schema nuages:
160
161 INTEGER iflagcon
162 INTEGER em_ifc(klev)
163
164 REAL em_pradj
165 REAL em_cldf(klev), em_cldq(klev)
166 REAL em_ftadj(klev), em_fradj(klev)
167
168 INTEGER ifc(klon, klev)
169 REAL pradj(klon)
170 REAL cldf(klon, klev), cldq(klon, klev)
171 REAL ftadj(klon, klev), fqadj(klon, klev)
172
173 ! sb --
174
175 ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
176
177 include "YOMCST.h"
178 include "YOETHF.h"
179 include "FCTTRE.h"
180
181 IF (first) THEN
182
183 ALLOCATE (em_t(klev))
184 ALLOCATE (em_q(klev))
185 ALLOCATE (em_qs(klev))
186 ALLOCATE (em_u(klev), em_v(klev), em_tra(klev,nbtr))
187 ALLOCATE (em_ph(klev+1), em_p(klev))
188 ALLOCATE (em_work1(klev), em_work2(klev))
189 ALLOCATE (em_d_t(klev), em_d_q(klev))
190 ALLOCATE (em_d_u(klev), em_d_v(klev), em_d_tra(klev,nbtr))
191 ALLOCATE (em_upwd(klev), em_dnwd(klev), em_dnwdbis(klev))
192 ALLOCATE (emmip(klev))
193 ALLOCATE (emmke(klev))
194 ALLOCATE (emma(klev))
195
196 first = .FALSE.
197 END IF
198
199 qcond_incld(:, :) = 0.
200
201 ! @$$ print*,'debut conema'
202
203 DO i = 1, klon
204 DO l = 1, klev + 1
205 em_ph(l) = paprs(i, l)/100.0
206 END DO
207
208 DO l = 1, klev
209 em_p(l) = pplay(i, l)/100.0
210 em_t(l) = t(i, l)
211 em_q(l) = q(i, l)
212 em_u(l) = u(i, l)
213 em_v(l) = v(i, l)
214 DO itra = 1, ntra
215 em_tra(l, itra) = tra(i, l, itra)
216 END DO
217 ! @$$ print*,'em_t',em_t
218 ! @$$ print*,'em_q',em_q
219 ! @$$ print*,'em_qs',em_qs
220 ! @$$ print*,'em_u',em_u
221 ! @$$ print*,'em_v',em_v
222 ! @$$ print*,'em_tra',em_tra
223 ! @$$ print*,'em_p',em_p
224
225
226
227 zx_t = em_t(l)
228 zdelta = max(0., sign(1.,rtt-zx_t))
229 zx_qs = r2es*foeew(zx_t, zdelta)/em_p(l)/100.0
230 zx_qs = min(0.5, zx_qs)
231 ! @$$ print*,'zx_qs',zx_qs
232 zcor = 1./(1.-retv*zx_qs)
233 zx_qs = zx_qs*zcor
234 em_qs(l) = zx_qs
235 ! @$$ print*,'em_qs',em_qs
236
237 em_work1(l) = work1(i, l)
238 em_work2(l) = work2(i, l)
239 emmke(l) = 0
240 ! emMa(l)=0
241 ! Ma(i,l)=0
242
243 em_dtvpdt1(l) = 0.
244 em_dtvpdq1(l) = 0.
245 dtvpdt1(i, l) = 0.
246 dtvpdq1(i, l) = 0.
247 END DO
248
249 em_dplcldt = 0.
250 em_dplcldr = 0.
251 rain(i) = 0.0
252 snow(i) = 0.0
253 kbas(i) = 1
254 ktop(i) = 1
255 ! ajout SB:
256 bas = 1
257 top = 1
258
259
260 ! sb3d write(*,1792) (em_work1(m),m=1,klev)
261 1792 FORMAT ('sig avant convect ', /, 10(1X,E13.5))
262
263 ! sb d write(*,1793) (em_work2(m),m=1,klev)
264 1793 FORMAT ('w avant convect ', /, 10(1X,E13.5))
265
266 ! @$$ print*,'avant convect'
267 ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
268
269
270 ! print*,'avant convect i=',i
271 CALL convect3(dtime, epmax, ok_adj_ema, em_t, em_q, em_qs, em_u, em_v, &
272 em_tra, em_p, em_ph, klev, klev+1, klev-1, ntra, dtime, iflag, em_d_t, &
273 em_d_q, em_d_u, em_d_v, em_d_tra, em_precip, em_bas, em_top, em_upwd, &
274 em_dnwd, em_dnwdbis, em_work1, em_work2, emmip, emmke, emma, ment, &
275 qent, tps, tls, sij, em_cape, em_tvp, em_pbase, em_bbase, em_dtvpdt1, &
276 em_dtvpdq1, em_dplcldt, em_dplcldr, & ! sbl
277 em_d_t2, em_d_q2, em_d_u2, em_d_v2, em_wd, em_qcond, em_qcondc) !sbl
278 ! print*,'apres convect '
279
280 ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
281
282 ! -- sb: Appel schema statistique de nuages couple a la convection
283 ! (Bony et Emanuel 2001):
284
285 ! -- creer cvthermo.h qui contiendra les cstes thermo de LMDZ:
286
287 iflagcon = 3
288 ! CALL cv_thermo(iflagcon)
289
290 ! -- appel schema de nuages:
291
292 ! CALL CLOUDS_SUB_LS(klev,em_q,em_qs,em_t
293 ! i ,em_p,em_ph,dtime,em_qcondc
294 ! o ,em_cldf,em_cldq,em_pradj,em_ftadj,em_fradj,em_ifc)
295
296 DO k = 1, klev
297 cldf(i, k) = em_cldf(k) ! cloud fraction (0-1)
298 cldq(i, k) = em_cldq(k) ! in-cloud water content (kg/kg)
299 ftadj(i, k) = em_ftadj(k) ! (dT/dt)_{LS adj} (K/s)
300 fqadj(i, k) = em_fradj(k) ! (dq/dt)_{LS adj} (kg/kg/s)
301 ifc(i, k) = em_ifc(k) ! flag convergence clouds_gno (1 ou 2)
302 END DO
303 pradj(i) = em_pradj ! precip from LS supersat adj (mm/day)
304
305 ! sb --
306
307 ! SB:
308 IF (iflag/=1 .AND. iflag/=4) THEN
309 em_cape = 0.
310 DO l = 1, klev
311 em_upwd(l) = 0.
312 em_dnwd(l) = 0.
313 em_dnwdbis(l) = 0.
314 emma(l) = 0.
315 em_tvp(l) = 0.
316 END DO
317 END IF
318 ! fin SB
319
320 ! If sig has been set to zero, then set Ma to zero
321
322 sigsum = 0.
323 DO k = 1, klev
324 sigsum = sigsum + em_work1(k)
325 END DO
326 IF (sigsum==0.0) THEN
327 DO k = 1, klev
328 emma(k) = 0.
329 END DO
330 END IF
331
332 ! sb3d print*,'i, iflag=',i,iflag
333
334 ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
335
336 ! SORTIE DES ICB ET INB
337 ! en fait inb et icb correspondent au niveau ou se trouve
338 ! le nuage,le numero d'interface
339 ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
340
341 ! modif SB:
342 IF (iflag==1 .OR. iflag==4) THEN
343 top = em_top
344 bas = em_bas
345 kbas(i) = em_bas
346 ktop(i) = em_top
347 END IF
348
349 pbase(i) = em_pbase
350 bbase(i) = em_bbase
351 rain(i) = em_precip/86400.0
352 snow(i) = 0.0
353 cape(i) = em_cape
354 wd(i) = em_wd
355 rflag(i) = real(iflag)
356 ! SB kbas(i) = em_bas
357 ! SB ktop(i) = em_top
358 dplcldt(i) = em_dplcldt
359 dplcldr(i) = em_dplcldr
360 DO l = 1, klev
361 d_t2(i, l) = dtime*em_d_t2(l)
362 d_q2(i, l) = dtime*em_d_q2(l)
363 d_u2(i, l) = dtime*em_d_u2(l)
364 d_v2(i, l) = dtime*em_d_v2(l)
365
366 d_t(i, l) = dtime*em_d_t(l)
367 d_q(i, l) = dtime*em_d_q(l)
368 d_u(i, l) = dtime*em_d_u(l)
369 d_v(i, l) = dtime*em_d_v(l)
370 DO itra = 1, ntra
371 d_tra(i, l, itra) = dtime*em_d_tra(l, itra)
372 END DO
373 upwd(i, l) = em_upwd(l)
374 dnwd(i, l) = em_dnwd(l)
375 dnwdbis(i, l) = em_dnwdbis(l)
376 work1(i, l) = em_work1(l)
377 work2(i, l) = em_work2(l)
378 ma(i, l) = emma(l)
379 tvp(i, l) = em_tvp(l)
380 dtvpdt1(i, l) = em_dtvpdt1(l)
381 dtvpdq1(i, l) = em_dtvpdq1(l)
382
383 IF (iflag_clw==0) THEN
384 qcond_incld(i, l) = em_qcondc(l)
385 ELSE IF (iflag_clw==1) THEN
386 qcond_incld(i, l) = em_qcond(l)
387 END IF
388 END DO
389 END DO
390
391 ! On calcule une eau liquide diagnostique en fonction de la
392 ! precip.
393 IF (iflag_clw==2) THEN
394 DO l = 1, klev
395 DO i = 1, klon
396 IF (ktop(i)-kbas(i)>0 .AND. l>=kbas(i) .AND. l<=ktop(i)) THEN
397 qcond_incld(i, l) = rain(i)*8.E4 & ! s *(pplay(i,l
398 ! )-paprs(i,ktop(i)+1))
399 /(pplay(i,kbas(i))-pplay(i,ktop(i)))
400 ! s **2
401 ELSE
402 qcond_incld(i, l) = 0.
403 END IF
404 END DO
405 PRINT *, 'l=', l, ', qcond_incld=', qcond_incld(1, l)
406 END DO
407 END IF
408
409
410 RETURN
411 END SUBROUTINE conema3
412
413