source: LMDZ4/branches/LMDZ4_par_0/libf/phylmd/conema3.F @ 5456

Last change on this file since 5456 was 634, checked in by Laurent Fairhead, 20 years ago

Modifications faites à la physique pour la rendre parallele YM
Une branche de travail LMDZ4_par_0 a été créée provisoirement afin de tester
les modifs pleinement avant leurs inclusions dans le tronc principal
LF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.5 KB
RevLine 
[524]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
[634]11      USE dimphy
[524]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"
[634]57cym#include "dimphy.h"
[524]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
[634]87      LOGICAL,SAVE :: first=.true.
88cym      REAL em_t(klev)
89      REAL,ALLOCATABLE,SAVE :: em_t(:)
90cym      REAL em_q(klev)
91      REAL,ALLOCATABLE,SAVE :: em_q(:)
92cym      REAL em_qs(klev)
93      REAL,ALLOCATABLE,SAVE :: em_qs(:) 
94cym      REAL em_u(klev), em_v(klev), em_tra(klev,ntrac)
95      REAL,ALLOCATABLE,SAVE :: em_u(:),em_v(:),em_tra(:,:)
96cym      REAL em_ph(klev+1), em_p(klev)
97      REAL,ALLOCATABLE,SAVE ::em_ph(:),em_p(:)
98cym      REAL em_work1(klev), em_work2(klev)
99      REAL,ALLOCATABLE,SAVE ::em_work1(:),em_work2(:)
100cym      REAL em_precip, em_d_t(klev), em_d_q(klev)
101      REAL,SAVE :: em_precip
102      REAL,ALLOCATABLE,SAVE :: em_d_t(:),em_d_q(:)
103cym      REAL em_d_u(klev), em_d_v(klev), em_d_tra(klev,ntrac)
104      REAL,ALLOCATABLE,SAVE ::em_d_u(:),em_d_v(:),em_d_tra(:,:)
105cym      REAL em_upwd(klev), em_dnwd(klev), em_dnwdbis(klev)
106      REAL,ALLOCATABLE,SAVE :: em_upwd(:),em_dnwd(:),em_dnwdbis(:)
[524]107      REAL em_dtvpdt1(klev), em_dtvpdq1(klev)
108      REAL em_dplcldt, em_dplcldr
[634]109cym      SAVE em_t,em_q, em_qs, em_ph, em_p, em_work1, em_work2
110cym      SAVE em_u,em_v, em_tra
111cym      SAVE em_d_u,em_d_v, em_d_tra
112cym      SAVE em_precip, em_d_t, em_d_q, em_upwd, em_dnwd, em_dnwdbis
113
[524]114      INTEGER em_bas, em_top
115      SAVE em_bas, em_top
116
117      REAL em_wd
118      REAL em_qcond(klev)
119      REAL em_qcondc(klev)
120c
121      REAL zx_t, zx_qs, zdelta, zcor
122      INTEGER iflag
123      REAL sigsum
124ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
125c     VARIABLES A SORTIR
126cccccccccccccccccccccccccccccccccccccccccccccccccc
127 
[634]128cym      REAL emmip(klev) !variation de flux ascnon dilue i et i+1
129      REAL,ALLOCATABLE,SAVE ::emmip(:)
130cym      SAVE emmip
131cym      real emMke(klev)
132      REAL,ALLOCATABLE,SAVE ::emMke(:)
133cym      save emMke
[524]134      real top
135      real bas
[634]136cym      real emMa(klev)
137      REAL,ALLOCATABLE,SAVE ::emMa(:)
138cym      save emMa
[524]139      real Ma(klon,klev)
140      real Ment(klev,klev)
141      real Qent(klev,klev)
142      real TPS(klev),TLS(klev)
143      real SIJ(klev,klev)
144      real em_CAPE, em_TVP(klev)
145      real em_pbase, em_bbase
146      integer iw,j,k,ix,iy
147
148c -- sb: pour schema nuages:
149
150       integer iflagcon
151       integer em_ifc(klev)
152     
153       real em_pradj
154       real em_cldf(klev), em_cldq(klev)
155       real em_ftadj(klev), em_fradj(klev)
156
157       integer ifc(klon,klev)
158       real pradj(klon)
159       real cldf(klon,klev), cldq(klon,klev)
160       real ftadj(klon,klev), fqadj(klon,klev)
161
162c sb --
163
164ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
165c
166#include "YOMCST.h"
167#include "YOETHF.h"
168#include "FCTTRE.h"
[634]169     
170      if (first) then
171 
172        allocate(em_t(klev))
173        allocate(em_q(klev))
174        allocate(em_qs(klev))
175        allocate(em_u(klev), em_v(klev), em_tra(klev,ntrac))
176        allocate(em_ph(klev+1), em_p(klev))
177        allocate(em_work1(klev), em_work2(klev))
178        allocate(em_d_t(klev), em_d_q(klev))
179        allocate(em_d_u(klev), em_d_v(klev), em_d_tra(klev,ntrac))
180        allocate(em_upwd(klev), em_dnwd(klev), em_dnwdbis(klev))
181        allocate(emmip(klev))
182        allocate(emMke(klev))
183        allocate(emMa(klev))
184 
185        first=.false.
186      endif
187 
[524]188      qcond_incld(:,:) = 0.
189c
190c$$$      print*,'debut conema'
191
192      DO 999 i = 1, klon
193      DO l = 1, klev+1
194         em_ph(l) = paprs(i,l) / 100.0
195      ENDDO
196c
197      DO l = 1, klev
198         em_p(l) = pplay(i,l) / 100.0
199         em_t(l) = t(i,l)
200         em_q(l) = q(i,l)
201         em_u(l) = u(i,l)
202         em_v(l) = v(i,l)
203         do itra = 1, ntra
204          em_tra(l,itra) = tra(i,l,itra)
205         enddo
206c$$$      print*,'em_t',em_t
207c$$$      print*,'em_q',em_q
208c$$$      print*,'em_qs',em_qs
209c$$$      print*,'em_u',em_u
210c$$$      print*,'em_v',em_v
211c$$$      print*,'em_tra',em_tra
212c$$$      print*,'em_p',em_p
213
214 
215c
216         zx_t = em_t(l)
217         zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
218         zx_qs= r2es * FOEEW(zx_t,zdelta)/em_p(l)/100.0
219         zx_qs=MIN(0.5,zx_qs)
220c$$$       print*,'zx_qs',zx_qs
221         zcor=1./(1.-retv*zx_qs)
222         zx_qs=zx_qs*zcor
223         em_qs(l) = zx_qs
224c$$$      print*,'em_qs',em_qs
225c
226         em_work1(l) = work1(i,l)
227         em_work2(l) = work2(i,l)
228         emMke(l)=0
229c        emMa(l)=0
230c        Ma(i,l)=0
231     
232         em_dtvpdt1(l) = 0.
233         em_dtvpdq1(l) = 0.
234         dtvpdt1(i,l) = 0.
235         dtvpdq1(i,l) = 0.
236      ENDDO
237c
238      em_dplcldt = 0.
239      em_dplcldr = 0.
240      rain(i) = 0.0
241      snow(i) = 0.0
242      kbas(i) = 1
243      ktop(i) = 1
244c ajout SB:
245      bas = 1
246      top = 1
247 
248 
249c sb3d      write(*,1792) (em_work1(m),m=1,klev)
2501792  format('sig avant convect ',/,10(1X,E13.5))
251c
252c sb d      write(*,1793) (em_work2(m),m=1,klev)
2531793  format('w avant convect ',/,10(1X,E13.5))
254 
255c$$$      print*,'avant convect'
256ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
257c
258
259c     print*,'avant convect i=',i
260      CALL convect3(dtime,epmax,ok_adj_ema,
261     .              em_t, em_q, em_qs,em_u ,em_v ,
262     .              em_tra, em_p, em_ph,
263     .              klev, klev+1, klev-1,ntra, dtime, iflag,
264     .              em_d_t, em_d_q,em_d_u,em_d_v,
265     .              em_d_tra, em_precip,
266     .              em_bas, em_top,em_upwd, em_dnwd, em_dnwdbis,
267     .              em_work1, em_work2,emmip,emMke,emMa,Ment,
268     .  Qent,TPS,TLS,SIJ,em_CAPE,em_TVP,em_pbase,em_bbase,
269     .  em_dtvpdt1,em_dtvpdq1,em_dplcldt,em_dplcldr, ! sbl
270     .  em_d_t2,em_d_q2,em_d_u2,em_d_v2,em_wd,em_qcond,em_qcondc)!sbl
271c     print*,'apres convect '
272c
273ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
274c
275c -- sb: Appel schema statistique de nuages couple a la convection
276c (Bony et Emanuel 2001):
277
278c -- creer cvthermo.h qui contiendra les cstes thermo de LMDZ:
279
280        iflagcon = 3
281c       CALL cv_thermo(iflagcon)
282
283c -- appel schema de nuages:
284
285c       CALL CLOUDS_SUB_LS(klev,em_q,em_qs,em_t
286c    i          ,em_p,em_ph,dtime,em_qcondc
287c    o          ,em_cldf,em_cldq,em_pradj,em_ftadj,em_fradj,em_ifc)
288
289        do k = 1, klev
290         cldf(i,k)  = em_cldf(k)  ! cloud fraction (0-1)
291         cldq(i,k)  = em_cldq(k)  ! in-cloud water content (kg/kg)
292         ftadj(i,k) = em_ftadj(k) ! (dT/dt)_{LS adj} (K/s)
293         fqadj(i,k) = em_fradj(k) ! (dq/dt)_{LS adj} (kg/kg/s)
294         ifc(i,k)   = em_ifc(k)   ! flag convergence clouds_gno (1 ou 2)
295        enddo
296        pradj(i) = em_pradj       ! precip from LS supersat adj (mm/day)
297
298c sb --
299c
300c SB:
301      if (iflag.ne.1 .and. iflag.ne.4) then
302         em_CAPE = 0.
303      do l = 1, klev
304         em_upwd(l) = 0.
305         em_dnwd(l) = 0.
306         em_dnwdbis(l) = 0.
307         emMa(l) = 0.
308         em_TVP(l) = 0.
309      enddo
310      endif
311c fin SB
312c
313c  If sig has been set to zero, then set Ma to zero
314c
315      sigsum = 0.
316      do k = 1,klev
317        sigsum = sigsum + em_work1(k)
318      enddo
319      if (sigsum .eq. 0.0) then
320        do k = 1,klev
321          emMa(k) = 0.
322        enddo
323      endif
324c
325c sb3d       print*,'i, iflag=',i,iflag
326c
327ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
328c
329c       SORTIE DES ICB ET INB
330c       en fait inb et icb correspondent au niveau ou se trouve
331c       le nuage,le numero d'interface
332cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
333 
334c modif SB:
335      if (iflag.EQ.1 .or. iflag.EQ.4) then
336       top=em_top
337       bas=em_bas
338       kbas(i) = em_bas
339       ktop(i) = em_top
340      endif
341 
342      pbase(i) = em_pbase
343      bbase(i) = em_bbase
344      rain(i) = em_precip/ 86400.0
345      snow(i) = 0.0
346      cape(i) = em_CAPE
347      wd(i) = em_wd
348      rflag(i) = float(iflag)
349c SB      kbas(i) = em_bas
350c SB      ktop(i) = em_top
351      dplcldt(i) = em_dplcldt
352      dplcldr(i) = em_dplcldr
353      DO l = 1, klev
354         d_t2(i,l) = dtime * em_d_t2(l)
355         d_q2(i,l) = dtime * em_d_q2(l)
356         d_u2(i,l) = dtime * em_d_u2(l)
357         d_v2(i,l) = dtime * em_d_v2(l)
358
359         d_t(i,l) = dtime * em_d_t(l)
360         d_q(i,l) = dtime * em_d_q(l)
361         d_u(i,l) = dtime * em_d_u(l)
362         d_v(i,l) = dtime * em_d_v(l)
363         do itra = 1, ntra
364         d_tra(i,l,itra) = dtime * em_d_tra(l,itra)
365         enddo
366         upwd(i,l) = em_upwd(l)
367         dnwd(i,l) = em_dnwd(l)
368         dnwdbis(i,l) = em_dnwdbis(l)
369         work1(i,l) = em_work1(l)
370         work2(i,l) = em_work2(l)
371         Ma(i,l)=emMa(l)
372         tvp(i,l)=em_TVP(l)
373         dtvpdt1(i,l) = em_dtvpdt1(l)
374         dtvpdq1(i,l) = em_dtvpdq1(l)
375
376         if (iflag_clw.eq.0) then
377            qcond_incld(i,l) = em_qcondc(l)
378         else if (iflag_clw.eq.1) then
379            qcond_incld(i,l) = em_qcond(l)
380         endif
381      ENDDO
382  999 CONTINUE
383
384c   On calcule une eau liquide diagnostique en fonction de la
385c  precip.
386      if ( iflag_clw.eq.2 ) then
387      do l=1,klev
388         do i=1,klon
389            if (ktop(i)-kbas(i).gt.0.and.
390     s         l.ge.kbas(i).and.l.le.ktop(i)) then
391               qcond_incld(i,l)=rain(i)*8.e4
392c    s         *(pplay(i,l      )-paprs(i,ktop(i)+1))
393     s         /(pplay(i,kbas(i))-pplay(i,ktop(i)))
394c    s         **2
395            else
396               qcond_incld(i,l)=0.
397            endif
398         enddo
399         print*,'l=',l,',   qcond_incld=',qcond_incld(1,l)
400      enddo
401      endif
402 
403
404      RETURN
405      END
406
Note: See TracBrowser for help on using the repository browser.