source: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_conema3.f90 @ 5157

Last change on this file since 5157 was 5153, checked in by abarral, 7 weeks ago

Revert FCTTRE to INCLUDE to assess impact of inlining

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