GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: phylmd/conema3.F90 Lines: 0 137 0.0 %
Date: 2023-06-30 12:51:15 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