71 integer :: ix,jy,kz,iz,n,l,kmin,kl,klp,ix1,jy1,ixp,jyp
72 integer :: rain_cloud_above,kz_inv
73 real ::
f_qvsat,pressure,rh,lsp,convp
74 real :: uvzlev(nuvzmax),wzlev(nwzmax),rhoh(nuvzmax),pinmconv(nzmax)
75 real :: uvwzlev(0:nxmaxn-1,0:nymaxn-1,nzmax)
76 real ::
ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
78 real :: dzdx1,dzdx2,dzdy1,dzdy2
79 real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
80 real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
81 real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
82 real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
83 real,
parameter :: const=r_air/ga
97 tvold=tt2n(ix,jy,1,n,l)*(1.+0.378*
ew(td2n(ix,jy,1,n,l))/ &
102 rhoh(1)=pold/(r_air*tvold)
109 pint=akz(kz)+bkz(kz)*psn(ix,jy,1,n,l)
110 tv=tthn(ix,jy,kz,n,l)*(1.+0.608*qvhn(ix,jy,kz,n,l))
111 rhoh(kz)=pint/(r_air*tv)
113 if (abs(tv-tvold).gt.0.2)
then
114 uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
115 (tv-tvold)/log(tv/tvold)
117 uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
126 wzlev(kz)=(uvzlev(kz+1)+uvzlev(kz))/2.
128 wzlev(nwz)=wzlev(nwz-1)+ &
129 uvzlev(nuvz)-uvzlev(nuvz-1)
141 uvwzlev(ix,jy,kz)=uvzlev(kz)
154 pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
155 ((aknew(2)+bknew(2)*psn(ix,jy,1,n,l))- &
156 (aknew(1)+bknew(1)*psn(ix,jy,1,n,l)))
158 pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
159 ((aknew(kz+1)+bknew(kz+1)*psn(ix,jy,1,n,l))- &
160 (aknew(kz-1)+bknew(kz-1)*psn(ix,jy,1,n,l)))
162 pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
163 ((aknew(nz)+bknew(nz)*psn(ix,jy,1,n,l))- &
164 (aknew(nz-1)+bknew(nz-1)*psn(ix,jy,1,n,l)))
170 uun(ix,jy,1,n,l)=uuhn(ix,jy,1,l)
171 vvn(ix,jy,1,n,l)=vvhn(ix,jy,1,l)
172 ttn(ix,jy,1,n,l)=tthn(ix,jy,1,n,l)
173 qvn(ix,jy,1,n,l)=qvhn(ix,jy,1,n,l)
174 pvn(ix,jy,1,n,l)=pvhn(ix,jy,1,l)
175 rhon(ix,jy,1,n,l)=rhoh(1)
176 uun(ix,jy,nz,n,l)=uuhn(ix,jy,nuvz,l)
177 vvn(ix,jy,nz,n,l)=vvhn(ix,jy,nuvz,l)
178 ttn(ix,jy,nz,n,l)=tthn(ix,jy,nuvz,n,l)
179 qvn(ix,jy,nz,n,l)=qvhn(ix,jy,nuvz,n,l)
180 pvn(ix,jy,nz,n,l)=pvhn(ix,jy,nuvz,l)
181 rhon(ix,jy,nz,n,l)=rhoh(nuvz)
185 if(height(iz).gt.uvzlev(nuvz))
then
186 uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l)
187 vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l)
188 ttn(ix,jy,iz,n,l)=ttn(ix,jy,nz,n,l)
189 qvn(ix,jy,iz,n,l)=qvn(ix,jy,nz,n,l)
190 pvn(ix,jy,iz,n,l)=pvn(ix,jy,nz,n,l)
191 rhon(ix,jy,iz,n,l)=rhon(ix,jy,nz,n,l)
194 if ((height(iz).gt.uvzlev(kz-1)).and. &
195 (height(iz).le.uvzlev(kz)))
then
196 dz1=height(iz)-uvzlev(kz-1)
197 dz2=uvzlev(kz)-height(iz)
199 uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+ &
200 uuhn(ix,jy,kz,l)*dz1)/dz
201 vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+ &
202 vvhn(ix,jy,kz,l)*dz1)/dz
203 ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2+ &
204 tthn(ix,jy,kz,n,l)*dz1)/dz
205 qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2+ &
206 qvhn(ix,jy,kz,n,l)*dz1)/dz
207 pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+ &
208 pvhn(ix,jy,kz,l)*dz1)/dz
209 rhon(ix,jy,iz,n,l)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
221 wwn(ix,jy,1,n,l)=wwhn(ix,jy,1,l)*pinmconv(1)
222 wwn(ix,jy,nz,n,l)=wwhn(ix,jy,nwz,l)*pinmconv(nz)
226 if ((height(iz).gt.wzlev(kz-1)).and. &
227 (height(iz).le.wzlev(kz)))
then
228 dz1=height(iz)-wzlev(kz-1)
229 dz2=wzlev(kz)-height(iz)
231 wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(kz-1)*dz2 &
232 +wwhn(ix,jy,kz,l)*pinmconv(kz)*dz1)/dz
243 drhodzn(ix,jy,1,n,l)=(rhon(ix,jy,2,n,l)-rhon(ix,jy,1,n,l))/ &
244 (height(2)-height(1))
246 drhodzn(ix,jy,kz,n,l)=(rhon(ix,jy,kz+1,n,l)- &
247 rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
249 drhodzn(ix,jy,nz,n,l)=drhodzn(ix,jy,nz-1,n,l)
266 ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/ &
267 cos((
real(jy)*dyn(l)+ylat0n(l))*pi180)
268 vi=vvn(ix,jy,iz,n,l)*dyconst*yresoln(l)
271 if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
272 (height(iz).le.uvwzlev(ix,jy,kz)))
then
273 dz1=height(iz)-uvwzlev(ix,jy,kz-1)
274 dz2=uvwzlev(ix,jy,kz)-height(iz)
288 dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
289 dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
290 dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
292 dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
293 dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
294 dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
296 wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l)+(dzdx*ui+dzdy*vi)
309 lsp=lsprecn(ix,jy,1,n,l)
310 convp=convprecn(ix,jy,1,n,l)
311 cloudsnh(ix,jy,n,l)=0
314 pressure=rhon(ix,jy,kz,n,l)*r_air*ttn(ix,jy,kz,n,l)
315 rh=qvn(ix,jy,kz,n,l)/
f_qvsat(pressure,ttn(ix,jy,kz,n,l))
316 cloudsn(ix,jy,kz,n,l)=0
318 if ((lsp.gt.0.01).or.(convp.gt.0.01))
then
320 cloudsnh(ix,jy,n,l)=cloudsnh(ix,jy,n,l)+ &
321 height(kz)-height(kz-1)
322 if (lsp.ge.convp)
then
323 cloudsn(ix,jy,kz,n,l)=3
325 cloudsn(ix,jy,kz,n,l)=2
328 cloudsn(ix,jy,kz,n,l)=1
331 if (rain_cloud_above.eq.1)
then
332 if (lsp.ge.convp)
then
333 cloudsn(ix,jy,kz,n,l)=5
335 cloudsn(ix,jy,kz,n,l)=4
real function f_qvsat(p, t)