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

Last change on this file since 5151 was 5144, checked in by abarral, 7 weeks ago

Put YOMCST.h into modules

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