source: trunk/libf/phylmd/conema3.F @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago

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