source: LMDZ.3.3/trunk/libf/phylmd/conema3.F @ 481

Last change on this file since 481 was 207, checked in by lmdz, 24 years ago

petit detail
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.7 KB
Line 
1c $Header$
2      SUBROUTINE conema (dtime,paprs,pplay,t,q,u,v,tra,ntra,
3     .             work1,work2,d_t,d_q,d_u,d_v,d_tra,
4     .             rain, snow, kbas, ktop,
5     .             upwd,dnwd,dnwdbis,bas,top,Ma,cape,tvp,rflag,
6     .             pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr)
7
8      IMPLICIT none
9c======================================================================
10c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
11c Objet: schema de convection de Emanuel (1991) interface
12c Mai 1998: Interface modifiee pour implementation dans LMDZ
13c======================================================================
14c Arguments:
15c dtime---input-R-pas d'integration (s)
16c paprs---input-R-pression inter-couches (Pa)
17c pplay---input-R-pression au milieu des couches (Pa)
18c t-------input-R-temperature (K)
19c q-------input-R-humidite specifique (kg/kg)
20c u-------input-R-vitesse du vent zonal (m/s)
21c v-------input-R-vitesse duvent meridien (m/s)
22c tra-----input-R-tableau de rapport de melange des traceurs
23c work*: input et output: deux variables de travail,
24c                            on peut les mettre a 0 au debut
25c
26C d_t-----output-R-increment de la temperature
27c d_q-----output-R-increment de la vapeur d'eau
28c d_u-----output-R-increment de la vitesse zonale
29c d_v-----output-R-increment de la vitesse meridienne
30c d_tra---output-R-increment du contenu en traceurs
31c rain----output-R-la pluie (mm/s)
32c snow----output-R-la neige (mm/s)
33c kbas----output-R-bas du nuage (integer)
34c ktop----output-R-haut du nuage (integer)
35c upwd----output-R-saturated updraft mass flux (kg/m**2/s)
36c dnwd----output-R-saturated downdraft mass flux (kg/m**2/s)
37c dnwdbis-output-R-unsaturated downdraft mass flux (kg/m**2/s)
38c bas-----output-R-bas du nuage (real)
39c top-----output-R-haut du nuage (real)
40c Ma------output-R-flux ascendant non dilue (kg/m**2/s)
41c cape----output-R-CAPE
42c tvp-----output-R-virtual temperature of the lifted parcel
43c rflag---output-R-flag sur le fonctionnement de convect
44c pbase---output-R-pression a la base du nuage (Pa)
45c bbase---output-R-buoyancy a la base du nuage (K)
46c dtvpdt1-output-R-derivative of parcel virtual temp wrt T1
47c dtvpdq1-output-R-derivative of parcel virtual temp wrt Q1
48c dplcldt-output-R-derivative of the PCP pressure wrt T1
49c dplcldr-output-R-derivative of the PCP pressure wrt Q1
50c======================================================================
51c
52#include "dimensions.h"
53#include "dimphy.h"
54      INTEGER i, l,m,itra
55      INTEGER ntra,ntrac !number of tracers; if no tracer transport
56                         ! is needed, set ntra = 1 (or 0)
57      PARAMETER (ntrac=nqmx-2)
58      REAL dtime
59c
60
61c
62      REAL paprs(klon,klev+1), pplay(klon,klev)
63      REAL t(klon,klev), q(klon,klev), d_t(klon,klev), d_q(klon,klev)
64      REAL u(klon,klev), v(klon,klev), tra(klon,klev,ntra)
65      REAL d_u(klon,klev), d_v(klon,klev), d_tra(klon,klev,ntra)
66      REAL work1(klon,klev), work2(klon,klev)
67      REAL upwd(klon,klev), dnwd(klon,klev), dnwdbis(klon,klev)
68      REAL rain(klon)
69      REAL snow(klon)
70      REAL cape(klon), tvp(klon,klev), rflag(klon)
71      REAL pbase(klon), bbase(klon)
72      REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
73      REAL dplcldt(klon), dplcldr(klon)
74      INTEGER kbas(klon), ktop(klon)
75c
76      REAL em_t(klev)
77      REAL em_q(klev)
78      REAL em_qs(klev)
79      REAL em_u(klev), em_v(klev), em_tra(klev,ntrac)
80      REAL em_ph(klev+1), em_p(klev)
81      REAL em_work1(klev), em_work2(klev)
82      REAL em_precip, em_d_t(klev), em_d_q(klev)
83      REAL em_d_u(klev), em_d_v(klev), em_d_tra(klev,ntrac)
84      REAL em_upwd(klev), em_dnwd(klev), em_dnwdbis(klev)
85      REAL em_dtvpdt1(klev), em_dtvpdq1(klev)
86      REAL em_dplcldt, em_dplcldr
87      SAVE em_t,em_q, em_qs, em_ph, em_p, em_work1, em_work2
88      SAVE em_u,em_v, em_tra
89      SAVE em_d_u,em_d_v, em_d_tra
90      SAVE em_precip, em_d_t, em_d_q, em_upwd, em_dnwd, em_dnwdbis
91      INTEGER em_bas, em_top
92      SAVE em_bas, em_top
93c
94      REAL zx_t, zx_qs, zdelta, zcor
95      INTEGER iflag
96      REAL sigsum
97ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
98c     VARIABLES A SORTIR
99cccccccccccccccccccccccccccccccccccccccccccccccccc
100 
101      REAL emmip(klev) !variation de flux ascnon dilue i et i+1
102      SAVE emmip
103      real emMke(klev)
104      save emMke
105      real top
106      real bas
107      real emMa(klev)
108      save emMa
109      real Ma(klon,klev)
110      real Ment(klev,klev)
111      real Qent(klev,klev)
112      real TPS(klev),TLS(klev)
113      real SIJ(klev,klev)
114      real em_CAPE, em_TVP(klev)
115      real em_pbase, em_bbase
116      integer iw,j,k,ix,iy
117ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
118c
119#include "YOMCST.h"
120#include "YOETHF.h"
121#include "FCTTRE.h"
122 
123c
124c$$$      print*,'debut conema'
125
126      DO 999 i = 1, klon
127      DO l = 1, klev+1
128         em_ph(l) = paprs(i,l) / 100.0
129      ENDDO
130c
131      DO l = 1, klev
132         em_p(l) = pplay(i,l) / 100.0
133         em_t(l) = t(i,l)
134         em_q(l) = q(i,l)
135         em_u(l) = u(i,l)
136         em_v(l) = v(i,l)
137         do itra = 1, ntra
138          em_tra(l,itra) = tra(i,l,itra)
139         enddo
140c$$$      print*,'em_t',em_t
141c$$$      print*,'em_q',em_q
142c$$$      print*,'em_qs',em_qs
143c$$$      print*,'em_u',em_u
144c$$$      print*,'em_v',em_v
145c$$$      print*,'em_tra',em_tra
146c$$$      print*,'em_p',em_p
147
148 
149c
150         zx_t = em_t(l)
151         zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
152         zx_qs= r2es * FOEEW(zx_t,zdelta)/em_p(l)/100.0
153         zx_qs=MIN(0.5,zx_qs)
154c$$$       print*,'zx_qs',zx_qs
155         zcor=1./(1.-retv*zx_qs)
156         zx_qs=zx_qs*zcor
157         em_qs(l) = zx_qs
158c$$$      print*,'em_qs',em_qs
159c
160         em_work1(l) = work1(i,l)
161         em_work2(l) = work2(i,l)
162         emMke(l)=0
163c        emMa(l)=0
164c        Ma(i,l)=0
165     
166         em_dtvpdt1(l) = 0.
167         em_dtvpdq1(l) = 0.
168         dtvpdt1(i,l) = 0.
169         dtvpdq1(i,l) = 0.
170      ENDDO
171c
172      em_dplcldt = 0.
173      em_dplcldr = 0.
174      rain(i) = 0.0
175      snow(i) = 0.0
176      kbas(i) = 1
177      ktop(i) = 1
178c ajout SB:
179      bas = 1
180      top = 1
181 
182 
183c sb3d      write(*,1792) (em_work1(m),m=1,klev)
1841792  format('sig avant convect ',/,10(1X,E13.5))
185c
186c sb d      write(*,1793) (em_work2(m),m=1,klev)
1871793  format('w avant convect ',/,10(1X,E13.5))
188 
189c$$$      print*,'avant convect'
190ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
191c
192      CALL convect(dtime, em_t, em_q, em_qs,em_u ,em_v ,
193     .              em_tra, em_p, em_ph,
194     .              klev, klev+1, klev-1,ntra, dtime, iflag,
195     .              em_d_t, em_d_q,em_d_u,em_d_v,
196     .              em_d_tra, em_precip,
197     .              em_bas, em_top,em_upwd, em_dnwd, em_dnwdbis,
198     .              em_work1, em_work2,emmip,emMke,emMa,Ment,
199c SB 11sept98     .              Qent,TPS,TLS,SIJ)
200c 19oct98     .              Qent,TPS,TLS,SIJ,em_CAPE,em_TVP)
201     .  Qent,TPS,TLS,SIJ,em_CAPE,em_TVP,em_pbase,em_bbase,
202     .  em_dtvpdt1,em_dtvpdq1,em_dplcldt,em_dplcldr)
203c
204ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
205c
206c SB:
207      if (iflag.ne.1 .and. iflag.ne.4) then
208         em_CAPE = 0.
209      do l = 1, klev
210         em_upwd(l) = 0.
211         em_dnwd(l) = 0.
212         em_dnwdbis(l) = 0.
213         emMa(l) = 0.
214         em_TVP(l) = 0.
215      enddo
216      endif
217c fin SB
218c
219c  If sig has been set to zero, then set Ma to zero
220c
221      sigsum = 0.
222      do k = 1,klev
223        sigsum = sigsum + em_work1(k)
224      enddo
225      if (sigsum .eq. 0.0) then
226        do k = 1,klev
227          emMa(k) = 0.
228        enddo
229      endif
230c
231c sb3d       print*,'i, iflag=',i,iflag
232c
233ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
234c
235c       SORTIE DES ICB ET INB
236c       en fait inb et icb correspondent au niveau ou se trouve
237c       le nuage,le numero d'interface
238cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
239 
240c modif SB:
241      if (iflag.EQ.1 .or. iflag.EQ.4) then
242       top=em_top
243       bas=em_bas
244       kbas(i) = em_bas
245       ktop(i) = em_top
246      endif
247 
248      pbase(i) = em_pbase
249      bbase(i) = em_bbase
250      rain(i) = em_precip/ 86400.0
251      snow(i) = 0.0
252      cape(i) = em_CAPE
253      rflag(i) = float(iflag)
254c SB      kbas(i) = em_bas
255c SB      ktop(i) = em_top
256      dplcldt(i) = em_dplcldt
257      dplcldr(i) = em_dplcldr
258      DO l = 1, klev
259         d_t(i,l) = dtime * em_d_t(l)
260         d_q(i,l) = dtime * em_d_q(l)
261         d_u(i,l) = dtime * em_d_u(l)
262         d_v(i,l) = dtime * em_d_v(l)
263         do itra = 1, ntra
264         d_tra(i,l,itra) = dtime * em_d_tra(l,itra)
265         enddo
266         upwd(i,l) = em_upwd(l)
267         dnwd(i,l) = em_dnwd(l)
268         dnwdbis(i,l) = em_dnwdbis(l)
269         work1(i,l) = em_work1(l)
270         work2(i,l) = em_work2(l)
271         Ma(i,l)=emMa(l)
272         tvp(i,l)=em_TVP(l)
273         dtvpdt1(i,l) = em_dtvpdt1(l)
274         dtvpdq1(i,l) = em_dtvpdq1(l)
275      ENDDO
276  999 CONTINUE
277 
278
279      RETURN
280      END
281
Note: See TracBrowser for help on using the repository browser.