source: LMDZ4/trunk/libf/phy_IPCC_AR4/conema3.F @ 932

Last change on this file since 932 was 868, checked in by Laurent Fairhead, 17 years ago

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