source: LMDZ5/trunk/libf/phylmd/conema3.F90 @ 2306

Last change on this file since 2306 was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

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