/[lmdze]/trunk/Sources/phylmd/CV30_routines/cv30_undilute1.f
ViewVC logotype

Contents of /trunk/Sources/phylmd/CV30_routines/cv30_undilute1.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 196 - (show annotations)
Mon May 23 13:50:39 2016 UTC (8 years ago) by guez
File size: 6326 byte(s)
Removed argument icbmax of cv30_feed, not used in cv_driver (not used
in LMDZ either).

Clearer to have iflag1 = 0 in cv30_feed than in cv_driver. Clearer to
have iflag1 = 42 in cv30_uncompress than in cv_driver.

Removed argument iflag of cv30_compress. Since there is iflag1 = 42 in
cv30_uncompress, iflag1 that comes out of cv_driver is disconnected
from iflag1 computed by cv30_feed and cv30_trigger.

Removed some references to convect3 and convect4 in comments. This
program derives from the convect3 version, we do not need to know
about other versions.

Bug fix in cv30_undilute1: icbs1 was not made >= 2.

1 module cv30_undilute1_m
2
3 implicit none
4
5 contains
6
7 SUBROUTINE cv30_undilute1(t1, q1, qs1, gz1, plcl1, p1, nk1, icb1, tp1, &
8 tvp1, clw1, icbs1)
9
10 ! UNDILUTE (ADIABATIC) UPDRAFT / 1st part
11 ! (up through ICB1 + 1)
12 ! Calculates the lifted parcel virtual temperature at nk1, the
13 ! actual temperature, and the adiabatic liquid water content.
14
15 ! Equivalent de TLIFT entre NK1 et ICB1+1 inclus
16
17 ! Differences with convect4:
18 ! - icbs1 is the first level above LCL (may differ from icb1)
19 ! - in the iterations, used x(icbs1) instead x(icb1)
20 ! - tvp1 is computed in only one time
21 ! - icbs1: first level above Plcl1 (IMIN de TLIFT) in output
22 ! - if icbs1=icb1, compute also tp1(icb1+1), tvp1(icb1+1) & clw1(icb1+1)
23
24 use cv30_param_m, only: minorig, nl
25 use cv_thermo_m, only: cl, clmcpv, cpd, cpv, eps, lv0, rrv
26 USE dimphy, ONLY: klev, klon
27
28 ! inputs:
29 integer, intent(in):: nk1(klon), icb1(klon)
30 real, intent(in):: t1(klon, klev)
31 real, intent(in):: q1(klon, klev), qs1(klon, klev), gz1(klon, klev)
32 real, intent(in):: p1(klon, klev)
33 real, intent(in):: plcl1(klon)
34
35 ! outputs:
36 real tp1(klon, klev), tvp1(klon, klev), clw1(klon, klev)
37
38 ! local variables:
39 integer i, k
40 integer icbs1(klon), icbsmax2
41 real tg, qg, alv, s, ahg, tc, denom, es
42 real ah0(klon), cpp(klon)
43 real tnk(klon), qnk(klon), gznk(klon), ticb(klon), gzicb(klon)
44 real qsicb(klon)
45 real cpinv(klon)
46
47 !-------------------------------------------------------------------
48
49 ! Calculates the lifted parcel virtual temperature at nk1,
50 ! the actual temperature, and the adiabatic
51 ! liquid water content. The procedure is to solve the equation.
52 ! cp*tp1+L*qp+phi=cp*tnk+L*qnk+gznk.
53
54 do i=1, klon
55 tnk(i)=t1(i, nk1(i))
56 qnk(i)=q1(i, nk1(i))
57 gznk(i)=gz1(i, nk1(i))
58 end do
59
60 ! *** Calculate certain parcel quantities, including static energy ***
61
62 do i=1, klon
63 ah0(i)=(cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) &
64 +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15))+gznk(i)
65 cpp(i)=cpd*(1.-qnk(i))+qnk(i)*cpv
66 cpinv(i)=1./cpp(i)
67 end do
68
69 ! *** Calculate lifted parcel quantities below cloud base ***
70
71 do i=1, klon
72 ! if icb1 is below LCL, start loop at ICB1+1:
73 ! (icbs1 est le premier niveau au-dessus du LCL)
74 icbs1(i)=MIN(max(icb1(i), 2), nl)
75 if (plcl1(i) < p1(i, icbs1(i))) then
76 icbs1(i)=MIN(icbs1(i)+1, nl)
77 endif
78 enddo
79
80 do i=1, klon
81 ticb(i)=t1(i, icbs1(i))
82 gzicb(i)=gz1(i, icbs1(i))
83 qsicb(i)=qs1(i, icbs1(i))
84 enddo
85
86 ! Re-compute icbsmax (icbsmax2):
87 icbsmax2=2
88 do i=1, klon
89 icbsmax2=max(icbsmax2, icbs1(i))
90 end do
91
92 ! initialization outputs:
93
94 do k=1, icbsmax2
95 do i=1, klon
96 tp1(i, k) = 0.0
97 tvp1(i, k) = 0.0
98 clw1(i, k) = 0.0
99 enddo
100 enddo
101
102 ! tp1 and tvp1 below cloud base:
103
104 do k=minorig, icbsmax2-1
105 do i=1, klon
106 tp1(i, k)=tnk(i)-(gz1(i, k)-gznk(i))*cpinv(i)
107 tvp1(i, k)=tp1(i, k)*(1.+qnk(i)/eps-qnk(i))
108 end do
109 end do
110
111 ! *** Find lifted parcel quantities above cloud base ***
112
113 do i=1, klon
114 tg=ticb(i)
115 qg=qsicb(i)
116 !debug alv=lv0-clmcpv*(ticb(i)-t0)
117 alv=lv0-clmcpv*(ticb(i)-273.15)
118
119 ! First iteration.
120
121 s=cpd*(1.-qnk(i))+cl*qnk(i) &
122 +alv*alv*qg/(rrv*ticb(i)*ticb(i))
123 s=1./s
124
125 ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i)
126 tg=tg+s*(ah0(i)-ahg)
127
128 !debug tc=tg-t0
129 tc=tg-273.15
130 denom=243.5+tc
131 denom=MAX(denom, 1.0)
132
133 es=6.112*exp(17.67*tc/denom)
134 qg=eps*es/(p1(i, icbs1(i))-es*(1.-eps))
135
136 ! Second iteration.
137
138 ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i)
139 tg=tg+s*(ah0(i)-ahg)
140
141 !debug tc=tg-t0
142 tc=tg-273.15
143 denom=243.5+tc
144 denom=MAX(denom, 1.0)
145
146 es=6.112*exp(17.67*tc/denom)
147
148 qg=eps*es/(p1(i, icbs1(i))-es*(1.-eps))
149
150 alv=lv0-clmcpv*(ticb(i)-273.15)
151
152 ! no approximation:
153 tp1(i, icbs1(i))=(ah0(i)-gz1(i, icbs1(i))-alv*qg) &
154 /(cpd+(cl-cpd)*qnk(i))
155
156 clw1(i, icbs1(i))=qnk(i)-qg
157 clw1(i, icbs1(i))=max(0.0, clw1(i, icbs1(i)))
158
159 ! (qg utilise au lieu du vrai mixing ratio rg)
160 tvp1(i, icbs1(i))=tp1(i, icbs1(i))*(1.+qg/eps-qnk(i))
161
162 end do
163
164 ! * icbs1 is the first level above the LCL:
165 ! if plcl1<p1(icb1), then icbs1=icb1+1
166 ! if plcl1>p1(icb1), then icbs1=icb1
167
168 ! * the routine above computes tvp1 from minorig to icbs1 (included).
169
170 ! * to compute buoybase (in cv30_trigger.F), both tvp1(icb1) and
171 ! tvp1(icb1+1) must be known. This is the case if icbs1=icb1+1,
172 ! but not if icbs1=icb1.
173
174 ! * therefore, in the case icbs1=icb1, we compute tvp1 at level icb1+1
175 ! (tvp1 at other levels will be computed in cv30_undilute2.F)
176
177 do i=1, klon
178 ticb(i)=t1(i, icb1(i)+1)
179 gzicb(i)=gz1(i, icb1(i)+1)
180 qsicb(i)=qs1(i, icb1(i)+1)
181 enddo
182
183 do i=1, klon
184 tg=ticb(i)
185 qg=qsicb(i)
186 !debug alv=lv0-clmcpv*(ticb(i)-t0)
187 alv=lv0-clmcpv*(ticb(i)-273.15)
188
189 ! First iteration.
190
191 s=cpd*(1.-qnk(i))+cl*qnk(i) &
192 +alv*alv*qg/(rrv*ticb(i)*ticb(i))
193 s=1./s
194
195 ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i)
196 tg=tg+s*(ah0(i)-ahg)
197
198 !debug tc=tg-t0
199 tc=tg-273.15
200 denom=243.5+tc
201 denom=MAX(denom, 1.0)
202
203 es=6.112*exp(17.67*tc/denom)
204
205 qg=eps*es/(p1(i, icb1(i)+1)-es*(1.-eps))
206
207 ! Second iteration.
208
209 ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i)
210 tg=tg+s*(ah0(i)-ahg)
211
212 !debug tc=tg-t0
213 tc=tg-273.15
214 denom=243.5+tc
215 denom=MAX(denom, 1.0)
216
217 es=6.112*exp(17.67*tc/denom)
218
219 qg=eps*es/(p1(i, icb1(i)+1)-es*(1.-eps))
220
221 alv=lv0-clmcpv*(ticb(i)-273.15)
222
223 ! no approximation:
224 tp1(i, icb1(i)+1)=(ah0(i)-gz1(i, icb1(i)+1)-alv*qg) &
225 /(cpd+(cl-cpd)*qnk(i))
226
227 clw1(i, icb1(i)+1)=qnk(i)-qg
228 clw1(i, icb1(i)+1)=max(0.0, clw1(i, icb1(i)+1))
229
230 ! (qg utilise au lieu du vrai mixing ratio rg)
231 tvp1(i, icb1(i)+1)=tp1(i, icb1(i)+1)*(1.+qg/eps-qnk(i)) !whole thing
232 end do
233
234 end SUBROUTINE cv30_undilute1
235
236 end module cv30_undilute1_m

  ViewVC Help
Powered by ViewVC 1.1.21