/[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 189 - (hide annotations)
Tue Mar 29 15:20:23 2016 UTC (8 years, 1 month ago) by guez
File size: 9508 byte(s)
There was a function gr_phy_write_3d in dyn3d and a function
gr_phy_write_2d in module grid_change. Moved them into a new module
gr_phy_write_m under a generic interface gr_phy_write. Replaced calls
to gr_fi_ecrit by calls to gr_phy_write.

Removed arguments len, nloc and nd of cv30_compress.

Removed arguments wd and wd1 of cv30_uncompress, wd of cv30_yield, wd
of concvl, wd1 of cv_driver. Was just filled with 0. Removed option
ok_gust in physiq, never used.

In cv30_unsat, cv30_yield and cv_driver, we only need to define b to
level nl - 1.

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

  ViewVC Help
Powered by ViewVC 1.1.21