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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 47 - (hide annotations)
Fri Jul 1 15:00:48 2011 UTC (12 years, 10 months ago) by guez
Original Path: trunk/libf/phylmd/CV3_routines/cv3_undilute1.f90
File size: 9353 byte(s)
Split "thermcell.f" and "cv3_routines.f".
Removed copies of files that are now in "L_util".
Moved "mva9" and "diagetpq" to their own files.
Unified variable names across procedures.

1 guez 47
2     SUBROUTINE cv3_undilute1(len,nd,t,q,qs,gz,plcl,p,nk,icb &
3     ,tp,tvp,clw,icbs)
4     use cvparam3
5     use cvthermo
6     implicit none
7    
8     !----------------------------------------------------------------
9     ! Equivalent de TLIFT entre NK et ICB+1 inclus
10     !
11     ! Differences with convect4:
12     ! - specify plcl in input
13     ! - icbs is the first level above LCL (may differ from icb)
14     ! - in the iterations, used x(icbs) instead x(icb)
15     ! - many minor differences in the iterations
16     ! - tvp is computed in only one time
17     ! - icbs: first level above Plcl (IMIN de TLIFT) in output
18     ! - if icbs=icb, compute also tp(icb+1),tvp(icb+1) & clw(icb+1)
19     !----------------------------------------------------------------
20    
21    
22     ! inputs:
23     integer len, nd
24     integer nk(len), icb(len)
25     real t(len,nd), q(len,nd), qs(len,nd), gz(len,nd)
26     real p(len,nd)
27     real plcl(len) ! convect3
28    
29     ! outputs:
30     real tp(len,nd), tvp(len,nd), clw(len,nd)
31    
32     ! local variables:
33     integer i, k
34     integer icb1(len), icbs(len), icbsmax2 ! convect3
35     real tg, qg, alv, s, ahg, tc, denom, es, rg
36     real ah0(len), cpp(len)
37     real tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len)
38     real qsicb(len) ! convect3
39     real cpinv(len) ! convect3
40    
41     !-------------------------------------------------------------------
42     ! --- Calculates the lifted parcel virtual temperature at nk,
43     ! --- the actual temperature, and the adiabatic
44     ! --- liquid water content. The procedure is to solve the equation.
45     ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
46     !-------------------------------------------------------------------
47    
48     do 320 i=1,len
49     tnk(i)=t(i,nk(i))
50     qnk(i)=q(i,nk(i))
51     gznk(i)=gz(i,nk(i))
52     ! ori ticb(i)=t(i,icb(i))
53     ! ori gzicb(i)=gz(i,icb(i))
54     320 continue
55     !
56     ! *** Calculate certain parcel quantities, including static energy ***
57     !
58     do 330 i=1,len
59     ah0(i)=(cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) &
60     +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15))+gznk(i)
61     cpp(i)=cpd*(1.-qnk(i))+qnk(i)*cpv
62     cpinv(i)=1./cpp(i)
63     330 continue
64     !
65     ! *** Calculate lifted parcel quantities below cloud base ***
66     !
67     do i=1,len !convect3
68     icb1(i)=MAX(icb(i),2) !convect3
69     icb1(i)=MIN(icb(i),nl) !convect3
70     ! if icb is below LCL, start loop at ICB+1:
71     ! (icbs est le premier niveau au-dessus du LCL)
72     icbs(i)=icb1(i) !convect3
73     if (plcl(i).lt.p(i,icb1(i))) then
74     icbs(i)=MIN(icbs(i)+1,nl) !convect3
75     endif
76     enddo !convect3
77    
78     do i=1,len !convect3
79     ticb(i)=t(i,icbs(i)) !convect3
80     gzicb(i)=gz(i,icbs(i)) !convect3
81     qsicb(i)=qs(i,icbs(i)) !convect3
82     enddo !convect3
83    
84     !
85     ! Re-compute icbsmax (icbsmax2): !convect3
86     ! !convect3
87     icbsmax2=2 !convect3
88     do 310 i=1,len !convect3
89     icbsmax2=max(icbsmax2,icbs(i)) !convect3
90     310 continue !convect3
91    
92     ! initialization outputs:
93    
94     do k=1,icbsmax2 ! convect3
95     do i=1,len ! convect3
96     tp(i,k) = 0.0 ! convect3
97     tvp(i,k) = 0.0 ! convect3
98     clw(i,k) = 0.0 ! convect3
99     enddo ! convect3
100     enddo ! convect3
101    
102     ! tp and tvp below cloud base:
103    
104     do 350 k=minorig,icbsmax2-1
105     do 340 i=1,len
106     tp(i,k)=tnk(i)-(gz(i,k)-gznk(i))*cpinv(i)
107     tvp(i,k)=tp(i,k)*(1.+qnk(i)/eps-qnk(i)) !whole thing (convect3)
108     340 continue
109     350 continue
110     !
111     ! *** Find lifted parcel quantities above cloud base ***
112     !
113     do 360 i=1,len
114     tg=ticb(i)
115     ! ori qg=qs(i,icb(i))
116     qg=qsicb(i) ! convect3
117     !debug alv=lv0-clmcpv*(ticb(i)-t0)
118     alv=lv0-clmcpv*(ticb(i)-273.15)
119     !
120     ! First iteration.
121     !
122     ! ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
123     s=cpd*(1.-qnk(i))+cl*qnk(i) &
124     +alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3
125     s=1./s
126     ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
127     ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3
128     tg=tg+s*(ah0(i)-ahg)
129     ! ori tg=max(tg,35.0)
130     !debug tc=tg-t0
131     tc=tg-273.15
132     denom=243.5+tc
133     denom=MAX(denom,1.0) ! convect3
134     ! ori if(tc.ge.0.0)then
135     es=6.112*exp(17.67*tc/denom)
136     ! ori else
137     ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
138     ! ori endif
139     ! ori qg=eps*es/(p(i,icb(i))-es*(1.-eps))
140     qg=eps*es/(p(i,icbs(i))-es*(1.-eps))
141     !
142     ! Second iteration.
143     !
144    
145     ! ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
146     ! ori s=1./s
147     ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
148     ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3
149     tg=tg+s*(ah0(i)-ahg)
150     ! ori tg=max(tg,35.0)
151     !debug tc=tg-t0
152     tc=tg-273.15
153     denom=243.5+tc
154     denom=MAX(denom,1.0) ! convect3
155     ! ori if(tc.ge.0.0)then
156     es=6.112*exp(17.67*tc/denom)
157     ! ori else
158     ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
159     ! ori end if
160     ! ori qg=eps*es/(p(i,icb(i))-es*(1.-eps))
161     qg=eps*es/(p(i,icbs(i))-es*(1.-eps))
162    
163     alv=lv0-clmcpv*(ticb(i)-273.15)
164    
165     ! ori c approximation here:
166     ! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
167     ! ori & -gz(i,icb(i))-alv*qg)/cpd
168    
169     ! convect3: no approximation:
170     tp(i,icbs(i))=(ah0(i)-gz(i,icbs(i))-alv*qg) &
171     /(cpd+(cl-cpd)*qnk(i))
172    
173     ! ori clw(i,icb(i))=qnk(i)-qg
174     ! ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
175     clw(i,icbs(i))=qnk(i)-qg
176     clw(i,icbs(i))=max(0.0,clw(i,icbs(i)))
177    
178     rg=qg/(1.-qnk(i))
179     ! ori tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
180     ! convect3: (qg utilise au lieu du vrai mixing ratio rg)
181     tvp(i,icbs(i))=tp(i,icbs(i))*(1.+qg/eps-qnk(i)) !whole thing
182    
183     360 continue
184     !
185     ! ori do 380 k=minorig,icbsmax2
186     ! ori do 370 i=1,len
187     ! ori tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)
188     ! ori 370 continue
189     ! ori 380 continue
190     !
191    
192     ! -- The following is only for convect3:
193     !
194     ! * icbs is the first level above the LCL:
195     ! if plcl<p(icb), then icbs=icb+1
196     ! if plcl>p(icb), then icbs=icb
197     !
198     ! * the routine above computes tvp from minorig to icbs (included).
199     !
200     ! * to compute buoybase (in cv3_trigger.F), both tvp(icb) and tvp(icb+1)
201     ! must be known. This is the case if icbs=icb+1, but not if icbs=icb.
202     !
203     ! * therefore, in the case icbs=icb, we compute tvp at level icb+1
204     ! (tvp at other levels will be computed in cv3_undilute2.F)
205     !
206    
207     do i=1,len
208     ticb(i)=t(i,icb(i)+1)
209     gzicb(i)=gz(i,icb(i)+1)
210     qsicb(i)=qs(i,icb(i)+1)
211     enddo
212    
213     do 460 i=1,len
214     tg=ticb(i)
215     qg=qsicb(i) ! convect3
216     !debug alv=lv0-clmcpv*(ticb(i)-t0)
217     alv=lv0-clmcpv*(ticb(i)-273.15)
218     !
219     ! First iteration.
220     !
221     ! ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
222     s=cpd*(1.-qnk(i))+cl*qnk(i) &
223     +alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3
224     s=1./s
225     ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
226     ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3
227     tg=tg+s*(ah0(i)-ahg)
228     ! ori tg=max(tg,35.0)
229     !debug tc=tg-t0
230     tc=tg-273.15
231     denom=243.5+tc
232     denom=MAX(denom,1.0) ! convect3
233     ! ori if(tc.ge.0.0)then
234     es=6.112*exp(17.67*tc/denom)
235     ! ori else
236     ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
237     ! ori endif
238     ! ori qg=eps*es/(p(i,icb(i))-es*(1.-eps))
239     qg=eps*es/(p(i,icb(i)+1)-es*(1.-eps))
240     !
241     ! Second iteration.
242     !
243    
244     ! ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
245     ! ori s=1./s
246     ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
247     ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3
248     tg=tg+s*(ah0(i)-ahg)
249     ! ori tg=max(tg,35.0)
250     !debug tc=tg-t0
251     tc=tg-273.15
252     denom=243.5+tc
253     denom=MAX(denom,1.0) ! convect3
254     ! ori if(tc.ge.0.0)then
255     es=6.112*exp(17.67*tc/denom)
256     ! ori else
257     ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
258     ! ori end if
259     ! ori qg=eps*es/(p(i,icb(i))-es*(1.-eps))
260     qg=eps*es/(p(i,icb(i)+1)-es*(1.-eps))
261    
262     alv=lv0-clmcpv*(ticb(i)-273.15)
263    
264     ! ori c approximation here:
265     ! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
266     ! ori & -gz(i,icb(i))-alv*qg)/cpd
267    
268     ! convect3: no approximation:
269     tp(i,icb(i)+1)=(ah0(i)-gz(i,icb(i)+1)-alv*qg) &
270     /(cpd+(cl-cpd)*qnk(i))
271    
272     ! ori clw(i,icb(i))=qnk(i)-qg
273     ! ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
274     clw(i,icb(i)+1)=qnk(i)-qg
275     clw(i,icb(i)+1)=max(0.0,clw(i,icb(i)+1))
276    
277     rg=qg/(1.-qnk(i))
278     ! ori tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
279     ! convect3: (qg utilise au lieu du vrai mixing ratio rg)
280     tvp(i,icb(i)+1)=tp(i,icb(i)+1)*(1.+qg/eps-qnk(i)) !whole thing
281    
282     460 continue
283    
284     return
285     end

  ViewVC Help
Powered by ViewVC 1.1.21