source: LMDZ6/trunk/libf/phylmd/conema3_mod.f90

Last change on this file was 6058, checked in by fhourdin, 5 weeks ago

Travail pour la replayisation de la convection

Reunion de tous les anciens common devenus modules, dans lmdz_cv_ini.
Déplacement de presque toutes les routines d'initialisation dans lmdz_cv_ini.
Encapsulage de certains sous-programmes dans des modules.
Suppression de programmes inutilisés (cv3_crit et cv3_incp)
Reste :

  • à sortir des routines d'initialisation "_pre" de cv_driver et

cva_driver

  • à passer le variables argunement en intent(in/out/inout).

La convergence numérique a été testée pour
iflag_con=3/30/4
en 3D parallèle.
La compilation de la version isotopique fonctionne.

  • 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.7 KB
Line 
1MODULE conema3_mod
2
3CONTAINS
4
5SUBROUTINE conema3(dtime, paprs, pplay, t, q, u, v, tra, ntra, work1, work2, &
6    d_t, d_q, d_u, d_v, d_tra, rain, snow, kbas, ktop, upwd, dnwd, dnwdbis, &
7    bas, top, ma, cape, tvp, rflag, pbase, bbase, dtvpdt1, dtvpdq1, dplcldt, &
8    dplcldr, qcond_incld)
9
10  USE dimphy
11  USE infotrac_phy, ONLY: nbtr
12! Replayisation USE yomcst_mod_h
13  USE lmdz_cv_ini, ONLY : retv,lv0,rtt,cpd
14
15! Replayisation USE conema3_mod_h
16  USE lmdz_cv_ini, ONLY : epmax,iflag_clw,ok_adj_ema
17  USE yoethf_mod_h
18  IMPLICIT NONE
19  ! ======================================================================
20  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
21  ! Objet: schema de convection de Emanuel (1991) interface
22  ! Mai 1998: Interface modifiee pour implementation dans LMDZ
23  ! ======================================================================
24  ! Arguments:
25  ! dtime---input-R-pas d'integration (s)
26  ! paprs---input-R-pression inter-couches (Pa)
27  ! pplay---input-R-pression au milieu des couches (Pa)
28  ! t-------input-R-temperature (K)
29  ! q-------input-R-humidite specifique (kg/kg)
30  ! u-------input-R-vitesse du vent zonal (m/s)
31  ! v-------input-R-vitesse duvent meridien (m/s)
32  ! tra-----input-R-tableau de rapport de melange des traceurs
33  ! work*: input et output: deux variables de travail,
34  ! on peut les mettre a 0 au debut
35
36  ! d_t-----output-R-increment de la temperature
37  ! d_q-----output-R-increment de la vapeur d'eau
38  ! d_u-----output-R-increment de la vitesse zonale
39  ! d_v-----output-R-increment de la vitesse meridienne
40  ! d_tra---output-R-increment du contenu en traceurs
41  ! rain----output-R-la pluie (mm/s)
42  ! snow----output-R-la neige (mm/s)
43  ! kbas----output-R-bas du nuage (integer)
44  ! ktop----output-R-haut du nuage (integer)
45  ! upwd----output-R-saturated updraft mass flux (kg/m**2/s)
46  ! dnwd----output-R-saturated downdraft mass flux (kg/m**2/s)
47  ! dnwdbis-output-R-unsaturated downdraft mass flux (kg/m**2/s)
48  ! bas-----output-R-bas du nuage (real)
49  ! top-----output-R-haut du nuage (real)
50  ! Ma------output-R-flux ascendant non dilue (kg/m**2/s)
51  ! cape----output-R-CAPE
52  ! tvp-----output-R-virtual temperature of the lifted parcel
53  ! rflag---output-R-flag sur le fonctionnement de convect
54  ! pbase---output-R-pression a la base du nuage (Pa)
55  ! bbase---output-R-buoyancy a la base du nuage (K)
56  ! dtvpdt1-output-R-derivative of parcel virtual temp wrt T1
57  ! dtvpdq1-output-R-derivative of parcel virtual temp wrt Q1
58  ! dplcldt-output-R-derivative of the PCP pressure wrt T1
59  ! dplcldr-output-R-derivative of the PCP pressure wrt Q1
60  ! ======================================================================
61
62  INTEGER i, l, m, itra
63  INTEGER ntra ! if no tracer transport
64    ! is needed, set ntra = 1 (or 0)
65  REAL dtime
66
67  REAL d_t2(klon, klev), d_q2(klon, klev) ! sbl
68  REAL d_u2(klon, klev), d_v2(klon, klev) ! sbl
69  REAL em_d_t2(klev), em_d_q2(klev) ! sbl
70  REAL em_d_u2(klev), em_d_v2(klev) ! sbl
71
72  REAL paprs(klon, klev+1), pplay(klon, klev)
73  REAL t(klon, klev), q(klon, klev), d_t(klon, klev), d_q(klon, klev)
74  REAL u(klon, klev), v(klon, klev), tra(klon, klev, ntra)
75  REAL d_u(klon, klev), d_v(klon, klev), d_tra(klon, klev, ntra)
76  REAL work1(klon, klev), work2(klon, klev)
77  REAL upwd(klon, klev), dnwd(klon, klev), dnwdbis(klon, klev)
78  REAL rain(klon)
79  REAL snow(klon)
80  REAL cape(klon), tvp(klon, klev), rflag(klon)
81  REAL pbase(klon), bbase(klon)
82  REAL dtvpdt1(klon, klev), dtvpdq1(klon, klev)
83  REAL dplcldt(klon), dplcldr(klon)
84  INTEGER kbas(klon), ktop(klon)
85
86  REAL wd(klon)
87  REAL qcond_incld(klon, klev)
88
89  LOGICAL, SAVE :: first = .TRUE.
90  !$OMP THREADPRIVATE(first)
91
92  ! ym      REAL em_t(klev)
93  REAL, ALLOCATABLE, SAVE :: em_t(:)
94  !$OMP THREADPRIVATE(em_t)
95  ! ym      REAL em_q(klev)
96  REAL, ALLOCATABLE, SAVE :: em_q(:)
97  !$OMP THREADPRIVATE(em_q)
98  ! ym      REAL em_qs(klev)
99  REAL, ALLOCATABLE, SAVE :: em_qs(:)
100  !$OMP THREADPRIVATE(em_qs)
101  ! ym      REAL em_u(klev), em_v(klev), em_tra(klev,nbtr)
102  REAL, ALLOCATABLE, SAVE :: em_u(:), em_v(:), em_tra(:, :)
103  !$OMP THREADPRIVATE(em_u,em_v,em_tra)
104  ! ym      REAL em_ph(klev+1), em_p(klev)
105  REAL, ALLOCATABLE, SAVE :: em_ph(:), em_p(:)
106  !$OMP THREADPRIVATE(em_ph,em_p)
107  ! ym      REAL em_work1(klev), em_work2(klev)
108  REAL, ALLOCATABLE, SAVE :: em_work1(:), em_work2(:)
109  !$OMP THREADPRIVATE(em_work1,em_work2)
110  ! ym      REAL em_precip, em_d_t(klev), em_d_q(klev)
111  REAL, SAVE :: em_precip
112  !$OMP THREADPRIVATE(em_precip)
113  REAL, ALLOCATABLE, SAVE :: em_d_t(:), em_d_q(:)
114  !$OMP THREADPRIVATE(em_d_t,em_d_q)
115  ! ym      REAL em_d_u(klev), em_d_v(klev), em_d_tra(klev,nbtr)
116  REAL, ALLOCATABLE, SAVE :: em_d_u(:), em_d_v(:), em_d_tra(:, :)
117  !$OMP THREADPRIVATE(em_d_u,em_d_v,em_d_tra)
118  ! ym      REAL em_upwd(klev), em_dnwd(klev), em_dnwdbis(klev)
119  REAL, ALLOCATABLE, SAVE :: em_upwd(:), em_dnwd(:), em_dnwdbis(:)
120  !$OMP THREADPRIVATE(em_upwd,em_dnwd,em_dnwdbis)
121  REAL em_dtvpdt1(klev), em_dtvpdq1(klev)
122  REAL em_dplcldt, em_dplcldr
123  ! ym      SAVE em_t,em_q, em_qs, em_ph, em_p, em_work1, em_work2
124  ! ym      SAVE em_u,em_v, em_tra
125  ! ym      SAVE em_d_u,em_d_v, em_d_tra
126  ! ym      SAVE em_precip, em_d_t, em_d_q, em_upwd, em_dnwd, em_dnwdbis
127
128  INTEGER em_bas, em_top
129  SAVE em_bas, em_top
130  !$OMP THREADPRIVATE(em_bas,em_top)
131  REAL em_wd
132  REAL em_qcond(klev)
133  REAL em_qcondc(klev)
134
135  REAL zx_t, zx_qs, zdelta, zcor
136  INTEGER iflag
137  REAL sigsum
138  ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
139  ! VARIABLES A SORTIR
140  ! ccccccccccccccccccccccccccccccccccccccccccccccccc
141
142  ! ym      REAL emmip(klev) !variation de flux ascnon dilue i et i+1
143  REAL, ALLOCATABLE, SAVE :: emmip(:)
144  !$OMP THREADPRIVATE(emmip)
145  ! ym      SAVE emmip
146  ! ym      real emMke(klev)
147  REAL, ALLOCATABLE, SAVE :: emmke(:)
148  !$OMP THREADPRIVATE(emmke)
149  ! ym      save emMke
150  REAL top
151  REAL bas
152  ! ym      real emMa(klev)
153  REAL, ALLOCATABLE, SAVE :: emma(:)
154  !$OMP THREADPRIVATE(emma)
155  ! ym      save emMa
156  REAL ma(klon, klev)
157  REAL ment(klev, klev)
158  REAL qent(klev, klev)
159  REAL tps(klev), tls(klev)
160  REAL sij(klev, klev)
161  REAL em_cape, em_tvp(klev)
162  REAL em_pbase, em_bbase
163  INTEGER iw, j, k, ix, iy
164
165  ! -- sb: pour schema nuages:
166
167  INTEGER iflagcon
168  INTEGER em_ifc(klev)
169
170  REAL em_pradj
171  REAL em_cldf(klev), em_cldq(klev)
172  REAL em_ftadj(klev), em_fradj(klev)
173
174  INTEGER ifc(klon, klev)
175  REAL pradj(klon)
176  REAL cldf(klon, klev), cldq(klon, klev)
177  REAL ftadj(klon, klev), fqadj(klon, klev)
178
179  include "FCTTRE.h"
180
181  REAL ::  rcpd,rlvtt
182
183  ! FH 2026/01/27 historiques : cohabitation de constantes cv et yomcst a nettoyer
184  ! rcpd used in FCTTHRE.h
185  rcpd=cpd
186  rlvtt=lv0
187
188
189  ! sb --
190
191
192  IF (first) THEN
193
194    ALLOCATE (em_t(klev))
195    ALLOCATE (em_q(klev))
196    ALLOCATE (em_qs(klev))
197    ALLOCATE (em_u(klev), em_v(klev), em_tra(klev,nbtr))
198    ALLOCATE (em_ph(klev+1), em_p(klev))
199    ALLOCATE (em_work1(klev), em_work2(klev))
200    ALLOCATE (em_d_t(klev), em_d_q(klev))
201    ALLOCATE (em_d_u(klev), em_d_v(klev), em_d_tra(klev,nbtr))
202    ALLOCATE (em_upwd(klev), em_dnwd(klev), em_dnwdbis(klev))
203    ALLOCATE (emmip(klev))
204    ALLOCATE (emmke(klev))
205    ALLOCATE (emma(klev))
206
207    first = .FALSE.
208  END IF
209
210  qcond_incld(:, :) = 0.
211
212  ! @$$      print*,'debut conema'
213
214  DO i = 1, klon
215    DO l = 1, klev + 1
216      em_ph(l) = paprs(i, l)/100.0
217    END DO
218
219    DO l = 1, klev
220      em_p(l) = pplay(i, l)/100.0
221      em_t(l) = t(i, l)
222      em_q(l) = q(i, l)
223      em_u(l) = u(i, l)
224      em_v(l) = v(i, l)
225      DO itra = 1, ntra
226        em_tra(l, itra) = tra(i, l, itra)
227      END DO
228      ! @$$      print*,'em_t',em_t
229      ! @$$      print*,'em_q',em_q
230      ! @$$      print*,'em_qs',em_qs
231      ! @$$      print*,'em_u',em_u
232      ! @$$      print*,'em_v',em_v
233      ! @$$      print*,'em_tra',em_tra
234      ! @$$      print*,'em_p',em_p
235
236
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)
2721792 FORMAT ('sig avant convect ', /, 10(1X,E13.5))
273
274    ! sb d      write(*,1793) (em_work2(m),m=1,klev)
2751793 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
421  RETURN
422END SUBROUTINE conema3
423
424END MODULE conema3_mod
Note: See TracBrowser for help on using the repository browser.