source: LMDZ5/branches/testing/libf/phylmd/conema3.F90 @ 5425

Last change on this file since 5425 was 2408, checked in by Laurent Fairhead, 9 years ago

Merged trunk changes r2298:2396 into testing branch

  • 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 2408 2015-12-14 10:43:09Z jyg $
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  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)
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.