My Project
 All Classes Files Functions Variables Macros
tlift.F
Go to the documentation of this file.
1 !
2 ! $Header$
3 !
4  SUBROUTINE tlift(P,T,RR,RS,GZ,PLCL,ICB,NK,
5  . tvp,tpk,clw,nd,nl,
6  . dtvpdt1,dtvpdq1)
7 c
8 c argument nk ajoute(jyg) = niveau de depart de la
9 c convection
10 c
11  parameter(na=60)
12  REAL gz(nd),tpk(nd),clw(nd)
13  REAL t(nd),rr(nd),rs(nd),tvp(nd),p(nd)
14  REAL dtvpdt1(nd),dtvpdq1(nd) ! Derivatives of parcel virtual
15 c temperature wrt t1 and q1
16 c
17  REAL clw_new(na),qi(na)
18  REAL dtpdt1(na),dtpdq1(na) ! Derivatives of parcel temperature
19 c wrt t1 and q1
20 
21 c
22  LOGICAL ice_conv
23 c
24 c *** assign values of thermodynamic constants ***
25 c
26 c sb cpd=1005.7
27 c sb cpv=1870.0
28 c sb cl=4190.0
29 c sb cpvmcl=2320.0
30 c sb rv=461.5
31 c sb rd=287.04
32 c sb eps=rd/rv
33 c sb alv0=2.501e6
34 ccccccccccccccccccccccc
35 c constantes coherentes avec le modele du centre europeen
36 c sb rd = 1000.0 * 1.380658e-23 * 6.0221367e+23 / 28.9644
37 c sb rv = 1000.0 * 1.380658e-23 * 6.0221367e+23 / 18.0153
38 c sb cpd = 3.5 * rd
39 c sb cpv = 4.0 * rv
40 c sb cl = 4218.0
41 c sb ci=2090.0
42 c sb cpvmcl=cl-cpv
43 c sb clmci=cl-ci
44 c sb eps=rd/rv
45 c sb alv0=2.5008e+06
46 c sb alf0=3.34e+05
47 
48 cccccccccccc
49 c on utilise les constantes thermo du centre europeen: (sb)
50 c
51 #include "YOMCST.h"
52  gravity = rg !sb: Pr que gravite ne devienne pas humidite!
53 c
54  cpd = rcpd
55  cpv = rcpv
56  cl = rcw
57  ci = rcs
58  cpvmcl = cl-cpv
59  clmci = cl-ci
60  eps = rd/rv
61  alv0 = rlvtt
62  alf0 = rlmlt ! (ALF0 = RLSTT-RLVTT)
63 c
64 cccccccccccccccccccccc
65 c
66 c *** calculate certain parcel quantities, including static energy ***
67 c
68  icb1=max(icb,2)
69  icb1=min(icb,nl)
70 c
71 cjyg1
72 cc cpp=cpd*(1.-rr(1))+rr(1)*cpv
73  cpp=cpd*(1.-rr(nk))+rr(nk)*cpv
74 cjyg2
75  cpinv=1./cpp
76 cjyg1
77 c icb may be below condensation level
78 ccc DO 100 i=1,icb1-1
79 ccc tpk(i)=t(1)-gz(i)*cpinv
80 ccc tvp(i)=tpk(i)*(1.+rr(1)/eps)
81  DO 50 i=1,icb1
82  clw(i)=0.0
83 50 CONTINUE
84 c
85  DO 100 i=nk,icb1
86  tpk(i)=t(nk)-(gz(i) - gz(nk))*cpinv
87 cjyg1
88 ccc tvp(i)=tpk(i)*(1.+rr(nk)/eps)
89  tvp(i)=tpk(i)*(1.+rr(nk)/eps-rr(nk))
90 cjyg2
91  dtvpdt1(i) = 1.+rr(nk)/eps-rr(nk)
92  dtvpdq1(i) = tpk(i)*(1./eps-1.)
93 c
94 cjyg2
95 
96  100 CONTINUE
97 
98 c
99 c *** find lifted parcel temperature and mixing ratio ***
100 c
101 cjyg1
102 cc ah0=(cpd*(1.-rr(1))+cl*rr(1))*t(1)
103 cc $ +rr(1)*(alv0-cpvmcl*(t(1)-273.15))
104  ah0=(cpd*(1.-rr(nk))+cl*rr(nk))*t(nk)
105  $ +rr(nk)*(alv0-cpvmcl*(t(nk)-273.15)) + gz(nk)
106 cjyg2
107 c
108 cjyg1
109  imin = icb1
110 c If icb is below lcl, start loop at icb+1
111  IF (plcl .LT. p(icb1)) imin = min(imin+1,nl)
112 c
113 ccc DO 300 i=icb1,nl
114  DO 300 i=imin,nl
115 cjyg2
116  alv=alv0-cpvmcl*(t(i)-273.15)
117  alf=alf0+clmci*(t(i)-273.15)
118 
119  rg=rs(i)
120  tg=t(i)
121 c s=cpd+alv*alv*rg/(rv*t(i)*t(i))
122 cjyg1
123 cc s=cpd*(1.-rr(1))+cl*rr(1)+alv*alv*rg/(rv*t(i)*t(i))
124  s=cpd*(1.-rr(nk))+cl*rr(nk)+alv*alv*rg/(rv*t(i)*t(i))
125 cjyg2
126  s=1./s
127 
128  DO 200 j=1,2
129 cjyg1
130 cc ahg=cpd*tg+(cl-cpd)*rr(1)*tg+alv*rg+gz(i)
131  ahg=cpd*tg+(cl-cpd)*rr(nk)*tg+alv*rg+gz(i)
132 cjyg2
133  tg=tg+s*(ah0-ahg)
134  tc=tg-273.15
135  denom=243.5+tc
136  denom=max(denom,1.0)
137 c
138 c formule de bolton pour psat
139 c
140  es=6.112*exp(17.67*tc/denom)
141  rg=eps*es/(p(i)-es*(1.-eps))
142 
143 
144  200 CONTINUE
145 
146 cjyg1
147 cc tpk(i)=(ah0-gz(i)-alv*rg)/(cpd+(cl-cpd)*rr(1))
148  tpk(i)=(ah0-gz(i)-alv*rg)/(cpd+(cl-cpd)*rr(nk))
149 cjyg2
150 c tpk(i)=(ah0-gz(i)-alv*rg-(cl-cpd)*t(i)*rr(1))/cpd
151 
152 cjyg1
153 cc clw(i)=rr(1)-rg
154  clw(i)=rr(nk)-rg
155 cjyg2
156  clw(i)=max(0.0,clw(i))
157 cjyg1
158 ccc tvp(i)=tpk(i)*(1.+rg/eps)
159  tvp(i)=tpk(i)*(1.+rg/eps-rr(nk))
160 cjyg2
161 c
162 cjyg1 derivatives
163 c
164  dtpdt1(i) = cpd*s
165  dtpdq1(i) = alv*s
166 c
167  dtvpdt1(i) = dtpdt1(i)*(1. + rg/eps -
168  . rr(nk) + alv*rg/(rd*tpk(i)) )
169  dtvpdq1(i) = dtpdq1(i)*(1. + rg/eps -
170  . rr(nk) + alv*rg/(rd*tpk(i)) ) - tpk(i)
171 c
172 cjyg2
173 
174  300 CONTINUE
175 c
176  ice_conv = .false.
177 
178  IF (ice_conv) THEN
179 c
180 cjam
181 c rajout de la procedure icefrac
182 c
183 c sb CALL icefrac(t,clw,clw_new,qi,nd,nl)
184 
185  DO 400 i=icb1,nl
186  IF (t(i).LT.263.15) THEN
187  tg=tpk(i)
188  tc=tpk(i)-273.15
189  denom=243.5+tc
190  es=6.112*exp(17.67*tc/denom)
191  alv=alv0-cpvmcl*(t(i)-273.15)
192  alf=alf0+clmci*(t(i)-273.15)
193 
194  DO j=1,4
195  esi=exp(23.33086-(6111.72784/tpk(i))+0.15215*log(tpk(i)))
196  qsat_new=eps*esi/(p(i)-esi*(1.-eps))
197 ccc snew= cpd*(1.-rr(1))+cl*rr(1)+alv*alv*qsat_new/(rv*tpk(i)*tpk(i))
198  snew= cpd*(1.-rr(nk))+cl*rr(nk)
199  . +alv*alv*qsat_new/(rv*tpk(i)*tpk(i))
200 c
201  snew=1./snew
202  tpk(i)=tg+(alf*qi(i)+alv*rg*(1.-(esi/es)))*snew
203 c@$$ print*,'################################'
204 c@$$ print*,tpk(i)
205 c@$$ print*,(alf*qi(i)+alv*rg*(1.-(esi/es)))*snew
206  ENDDO
207 ccc clw(i)=rr(1)-qsat_new
208  clw(i)=rr(nk)-qsat_new
209  clw(i)=max(0.0,clw(i))
210 cjyg1
211 ccc tvp(i)=tpk(i)*(1.+qsat_new/eps)
212  tvp(i)=tpk(i)*(1.+qsat_new/eps-rr(nk))
213 cjyg2
214  ELSE
215  CONTINUE
216  ENDIF
217 
218  400 CONTINUE
219 c
220  ENDIF
221 c
222 
223 ******************************************************
224 ** bk : rajout de la temperature des ascendances
225 ** non dilues au niveau klev = nd
226 ** posons le environ egal a celui de klev-1
227 ********************************************************
228 
229  tpk(nl+1)=tpk(nl)
230 
231 *******************************************************
232 
233  rg = gravity ! RG redevient la gravite de YOMCST (sb)
234 
235 
236  RETURN
237  END
238 
239 
240 
241 
242 
243 
244 
245