source: LMDZ6/branches/contrails/libf/phylmd/conema3.f90 @ 5452

Last change on this file since 5452 was 5400, checked in by evignon, 3 weeks ago

ajout de omp threadprivate manquants

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.5 KB
Line 
1
2! $Id: conema3.f90 5400 2024-12-10 10:35:28Z aborella $
3
4SUBROUTINE 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  USE yomcst_mod_h
12  USE conema3_mod_h
13  USE yoethf_mod_h
14  IMPLICIT NONE
15  ! ======================================================================
16  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
17  ! Objet: schema de convection de Emanuel (1991) interface
18  ! Mai 1998: Interface modifiee pour implementation dans LMDZ
19  ! ======================================================================
20  ! Arguments:
21  ! dtime---input-R-pas d'integration (s)
22  ! paprs---input-R-pression inter-couches (Pa)
23  ! pplay---input-R-pression au milieu des couches (Pa)
24  ! t-------input-R-temperature (K)
25  ! q-------input-R-humidite specifique (kg/kg)
26  ! u-------input-R-vitesse du vent zonal (m/s)
27  ! v-------input-R-vitesse duvent meridien (m/s)
28  ! tra-----input-R-tableau de rapport de melange des traceurs
29  ! work*: input et output: deux variables de travail,
30  ! on peut les mettre a 0 au debut
31
32  ! d_t-----output-R-increment de la temperature
33  ! d_q-----output-R-increment de la vapeur d'eau
34  ! d_u-----output-R-increment de la vitesse zonale
35  ! d_v-----output-R-increment de la vitesse meridienne
36  ! d_tra---output-R-increment du contenu en traceurs
37  ! rain----output-R-la pluie (mm/s)
38  ! snow----output-R-la neige (mm/s)
39  ! kbas----output-R-bas du nuage (integer)
40  ! ktop----output-R-haut du nuage (integer)
41  ! upwd----output-R-saturated updraft mass flux (kg/m**2/s)
42  ! dnwd----output-R-saturated downdraft mass flux (kg/m**2/s)
43  ! dnwdbis-output-R-unsaturated downdraft mass flux (kg/m**2/s)
44  ! bas-----output-R-bas du nuage (real)
45  ! top-----output-R-haut du nuage (real)
46  ! Ma------output-R-flux ascendant non dilue (kg/m**2/s)
47  ! cape----output-R-CAPE
48  ! tvp-----output-R-virtual temperature of the lifted parcel
49  ! rflag---output-R-flag sur le fonctionnement de convect
50  ! pbase---output-R-pression a la base du nuage (Pa)
51  ! bbase---output-R-buoyancy a la base du nuage (K)
52  ! dtvpdt1-output-R-derivative of parcel virtual temp wrt T1
53  ! dtvpdq1-output-R-derivative of parcel virtual temp wrt Q1
54  ! dplcldt-output-R-derivative of the PCP pressure wrt T1
55  ! dplcldr-output-R-derivative of the PCP pressure wrt Q1
56  ! ======================================================================
57
58  INTEGER i, l, m, itra
59  INTEGER ntra ! if no tracer transport
60    ! is needed, set ntra = 1 (or 0)
61  REAL dtime
62
63  REAL d_t2(klon, klev), d_q2(klon, klev) ! sbl
64  REAL d_u2(klon, klev), d_v2(klon, klev) ! sbl
65  REAL em_d_t2(klev), em_d_q2(klev) ! sbl
66  REAL em_d_u2(klev), em_d_v2(klev) ! sbl
67
68  REAL paprs(klon, klev+1), pplay(klon, klev)
69  REAL t(klon, klev), q(klon, klev), d_t(klon, klev), d_q(klon, klev)
70  REAL u(klon, klev), v(klon, klev), tra(klon, klev, ntra)
71  REAL d_u(klon, klev), d_v(klon, klev), d_tra(klon, klev, ntra)
72  REAL work1(klon, klev), work2(klon, klev)
73  REAL upwd(klon, klev), dnwd(klon, klev), dnwdbis(klon, klev)
74  REAL rain(klon)
75  REAL snow(klon)
76  REAL cape(klon), tvp(klon, klev), rflag(klon)
77  REAL pbase(klon), bbase(klon)
78  REAL dtvpdt1(klon, klev), dtvpdq1(klon, klev)
79  REAL dplcldt(klon), dplcldr(klon)
80  INTEGER kbas(klon), ktop(klon)
81
82  REAL wd(klon)
83  REAL qcond_incld(klon, klev)
84
85  LOGICAL, SAVE :: first = .TRUE.
86  !$OMP THREADPRIVATE(first)
87
88  ! ym      REAL em_t(klev)
89  REAL, ALLOCATABLE, SAVE :: em_t(:)
90  !$OMP THREADPRIVATE(em_t)
91  ! ym      REAL em_q(klev)
92  REAL, ALLOCATABLE, SAVE :: em_q(:)
93  !$OMP THREADPRIVATE(em_q)
94  ! ym      REAL em_qs(klev)
95  REAL, ALLOCATABLE, SAVE :: em_qs(:)
96  !$OMP THREADPRIVATE(em_qs)
97  ! ym      REAL em_u(klev), em_v(klev), em_tra(klev,nbtr)
98  REAL, ALLOCATABLE, SAVE :: em_u(:), em_v(:), em_tra(:, :)
99  !$OMP THREADPRIVATE(em_u,em_v,em_tra)
100  ! ym      REAL em_ph(klev+1), em_p(klev)
101  REAL, ALLOCATABLE, SAVE :: em_ph(:), em_p(:)
102  !$OMP THREADPRIVATE(em_ph,em_p)
103  ! ym      REAL em_work1(klev), em_work2(klev)
104  REAL, ALLOCATABLE, SAVE :: em_work1(:), em_work2(:)
105  !$OMP THREADPRIVATE(em_work1,em_work2)
106  ! ym      REAL em_precip, em_d_t(klev), em_d_q(klev)
107  REAL, SAVE :: em_precip
108  !$OMP THREADPRIVATE(em_precip)
109  REAL, ALLOCATABLE, SAVE :: em_d_t(:), em_d_q(:)
110  !$OMP THREADPRIVATE(em_d_t,em_d_q)
111  ! ym      REAL em_d_u(klev), em_d_v(klev), em_d_tra(klev,nbtr)
112  REAL, ALLOCATABLE, SAVE :: em_d_u(:), em_d_v(:), em_d_tra(:, :)
113  !$OMP THREADPRIVATE(em_d_u,em_d_v,em_d_tra)
114  ! ym      REAL em_upwd(klev), em_dnwd(klev), em_dnwdbis(klev)
115  REAL, ALLOCATABLE, SAVE :: em_upwd(:), em_dnwd(:), em_dnwdbis(:)
116  !$OMP THREADPRIVATE(em_upwd,em_dnwd,em_dnwdbis)
117  REAL em_dtvpdt1(klev), em_dtvpdq1(klev)
118  REAL em_dplcldt, em_dplcldr
119  ! ym      SAVE em_t,em_q, em_qs, em_ph, em_p, em_work1, em_work2
120  ! ym      SAVE em_u,em_v, em_tra
121  ! ym      SAVE em_d_u,em_d_v, em_d_tra
122  ! ym      SAVE em_precip, em_d_t, em_d_q, em_upwd, em_dnwd, em_dnwdbis
123
124  INTEGER em_bas, em_top
125  SAVE em_bas, em_top
126  !$OMP THREADPRIVATE(em_bas,em_top)
127  REAL em_wd
128  REAL em_qcond(klev)
129  REAL em_qcondc(klev)
130
131  REAL zx_t, zx_qs, zdelta, zcor
132  INTEGER iflag
133  REAL sigsum
134  ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
135  ! VARIABLES A SORTIR
136  ! ccccccccccccccccccccccccccccccccccccccccccccccccc
137
138  ! ym      REAL emmip(klev) !variation de flux ascnon dilue i et i+1
139  REAL, ALLOCATABLE, SAVE :: emmip(:)
140  !$OMP THREADPRIVATE(emmip)
141  ! ym      SAVE emmip
142  ! ym      real emMke(klev)
143  REAL, ALLOCATABLE, SAVE :: emmke(:)
144  !$OMP THREADPRIVATE(emmke)
145  ! ym      save emMke
146  REAL top
147  REAL bas
148  ! ym      real emMa(klev)
149  REAL, ALLOCATABLE, SAVE :: emma(:)
150  !$OMP THREADPRIVATE(emma)
151  ! ym      save emMa
152  REAL ma(klon, klev)
153  REAL ment(klev, klev)
154  REAL qent(klev, klev)
155  REAL tps(klev), tls(klev)
156  REAL sij(klev, klev)
157  REAL em_cape, em_tvp(klev)
158  REAL em_pbase, em_bbase
159  INTEGER iw, j, k, ix, iy
160
161  ! -- sb: pour schema nuages:
162
163  INTEGER iflagcon
164  INTEGER em_ifc(klev)
165
166  REAL em_pradj
167  REAL em_cldf(klev), em_cldq(klev)
168  REAL em_ftadj(klev), em_fradj(klev)
169
170  INTEGER ifc(klon, klev)
171  REAL pradj(klon)
172  REAL cldf(klon, klev), cldq(klon, klev)
173  REAL ftadj(klon, klev), fqadj(klon, klev)
174
175  ! sb --
176
177  ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
178
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)
2611792 FORMAT ('sig avant convect ', /, 10(1X,E13.5))
262
263    ! sb d      write(*,1793) (em_work2(m),m=1,klev)
2641793 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
411END SUBROUTINE conema3
412
Note: See TracBrowser for help on using the repository browser.