Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

verttransform_nests.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 11 years ago
sources for flexwrf v3.1
File size: 20.8 KB

Line
1	!***********************************************************************
2	!* Copyright 2012,2013 *
3	!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine, *
4	!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
5	!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso, *
6	!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann, *
7	!* *
8	!* This file is part of FLEXPART WRF *
9	!* *
10	!* FLEXPART is free software: you can redistribute it and/or modify *
11	!* it under the terms of the GNU General Public License as published by*
12	!* the Free Software Foundation, either version 3 of the License, or *
13	!* (at your option) any later version. *
14	!* *
15	!* FLEXPART is distributed in the hope that it will be useful, *
16	!* but WITHOUT ANY WARRANTY; without even the implied warranty of *
17	!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
18	!* GNU General Public License for more details. *
19	!* *
20	!* You should have received a copy of the GNU General Public License *
21	!* along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
22	!***********************************************************************
23	subroutine verttransform_nests(n,uuhn,vvhn,wwhn,pvhn,divhn)
24	! i i i i i
25	!*******************************************************************************
26	! *
27	! Note: This is the FLEXPART_WRF version of subroutine verttransform_nests. *
28	! The computational grid is the WRF x-y grid rather than lat-lon. *
29	! *
30	! This subroutine transforms temperature, dew point temperature and *
31	! wind components from eta to meter coordinates. *
32	! The vertical wind component is transformed from Pa/s to m/s using *
33	! the conversion factor pinmconv. *
34	! In addition, this routine calculates vertical density gradients *
35	! needed for the parameterization of the turbulent velocities. *
36	! It is similar to verttransform, but makes the transformations for *
37	! the nested grids. *
38	! *
39	! Author: A. Stohl, G. Wotawa *
40	! *
41	! 12 August 1996 *
42	! Update: 16 January 1998 *
43	! *
44	! Major update: 17 February 1999 *
45	! by G. Wotawa *
46	! *
47	! - Vertical levels for u, v and w are put together *
48	! - Slope correction for vertical velocity: Modification of calculation *
49	! procedure *
50	! *
51	! Changes, Bernd C. Krueger, Feb. 2001: (marked "C-cv") *
52	! Variables tthn and qvhn (on eta coordinates) from common block *
53	! *
54	! 16 Nov 2005, R. Easter - changes for FLEXPART_WRF *
55	! 17 Nov 2005 - R. Easter - terrain correction applied to ww. There are *
56	! now 3 options, controlled by "method_w_terrain_correction" *
57	! *
58	! 11 June 2007 -- convert TKEhn to tken
59	! 25 June 2007 -- convert ptthn to pttn
60	! Jan 2012, J Brioude: modified to handle different wind options and openmp
61	!*******************************************************************************
62	! *
63	! Variables: *
64	! nxn,nyn,nuvz,nwz field dimensions in x,y and z direction *
65	! uun wind components in x-direction [m/s] *
66	! vvn wind components in y-direction [m/s] *
67	! wwn wind components in z-direction [deltaeta/s] *
68	! ttn temperature [K] *
69	! pvn potential vorticity (pvu) *
70	! psn surface pressure [Pa] *
71	! *
72	!*******************************************************************************
73
74	use par_mod
75	use com_mod
76
77	! include 'includepar'
78	! include 'includecom'
79	implicit none
80	integer :: ix,jy,kz,iz,n,l,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
81	integer :: icloudtop
82	real :: rh,lsp,convp,prec,rhmin
83	integer :: method_z_compute,aa,dimx,dimy
84	real :: uvzlev(nuvzmax),wzlev(nwzmax),rhoh(nuvzmax),pinmconv(nzmax)
85	real :: uvwzlev(0:nxmaxn-1,0:nymaxn-1,nzmax)
86	real :: ew,pint,tv,tvold,pold,const,dz1,dz2,dz,ui,vi
87	real :: dzdx,dzdy
88	real :: dzdx1,dzdx2,dzdy1,dzdy2
89	real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
90	real :: divn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
91	real(kind=4) :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
92	real(kind=4) :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
93	real(kind=4) :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
94	real(kind=4) :: divhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
95	real :: wwhn_svaa(nwzmax),u,v
96	parameter(const=r_air/ga)
97	! integer :: rain_cloud_above,kz_inv
98
99	real :: f_qvsat,pressure
100	! real :: rh,lsp,convp
101	real,parameter :: precmin = 0.002
102
103	! CDA
104	logical :: lconvectprec = .true.
105
106
107	! set method_z_compute
108	method_z_compute = 10
109
110
111	! Loop over all nests
112	!********************
113
114	do l=1,numbnests
115	dimy=nyn(l)-1
116	dimx=nxn(l)-1
117	! print*,'start omp '
118	! Loop over the whole grid
119	!*************************
120	!!!$OMP PARALLEL DEFAULT(SHARED) &
121	!!!$OMP PRIVATE(ix,jy,ixm,jym,tvold,pold,pint,tv,rhoh,uvzlev,wzlev, &
122	!!!$OMP uvwzlev,pinmconv,kz,iz,kmin,dz1,dz2,dz,ix1,jy1,ixp,jyp, &
123	!!!$OMP dzdy,dzdx,aa,u,v )
124	!$OMP DO
125	do jy=0,dimy
126	do ix=0,dimx
127
128	tvold=tt2n(ix,jy,1,n,l)(1.+0.378ew(td2n(ix,jy,1,n,l))/ &
129	psn(ix,jy,1,n,l))
130	pold=psn(ix,jy,1,n,l)
131	uvzlev(1)=0.
132	wzlev(1)=0.
133	rhoh(1)=pold/(r_air*tvold)
134
135
136	! Compute heights of eta levels
137	!******************************
138
139	do kz=2,nuvz
140	! FLEXPART_WRF - pphn hold pressure
141	! pint=akz(kz)+bkz(kz)*psn(ix,jy,1,n,l)
142	pint=pphn(ix,jy,kz,n,l)
143	tv=tthn(ix,jy,kz,n,l)(1.+0.608qvhn(ix,jy,kz,n,l))
144	rhoh(kz)=pint/(r_air*tv)
145
146	if (abs(tv-tvold).gt.0.2) then
147	uvzlev(kz)=uvzlev(kz-1)+constlog(pold/pint) &
148	(tv-tvold)/log(tv/tvold)
149	else
150	uvzlev(kz)=uvzlev(kz-1)+constlog(pold/pint)tv
151	endif
152
153	tvold=tv
154	pold=pint
155	end do
156
157
158	! print*,'etap 1',ix,jy
159
160	do kz=2,nwz-1
161	wzlev(kz)=(uvzlev(kz+1)+uvzlev(kz))/2.
162	end do
163	wzlev(nwz)=wzlev(nwz-1)+ &
164	uvzlev(nuvz)-uvzlev(nuvz-1)
165
166	! FLEXPART_WRF - get uvzlev & wzlev from zzh
167	if (method_z_compute .eq. 10) then
168	do kz = 2, nuvz
169	if ((add_sfc_level .eq. 1) .and. (kz .eq. 2)) then
170	uvzlev(kz) = 0.5*(zzhn(ix,jy,3,n,l) + &
171	zzhn(ix,jy,1,n,l)) &
172	- zzhn(ix,jy,1,n,l)
173	else
174	uvzlev(kz) = 0.5*(zzhn(ix,jy,kz+1,n,l) + &
175	zzhn(ix,jy,kz ,n,l)) &
176	- zzhn(ix,jy,1,n,l)
177	end if
178	end do
179	do kz = 2, nwz
180	wzlev(kz) = zzhn(ix,jy,kz+add_sfc_level,n,l) &
181	- zzhn(ix,jy,1,n,l)
182	end do
183	end if
184
185	! print*,'etap 2',ix,jy
186	! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
187	! upon the transformation to z levels. In order to save computer memory, this is
188	! not done anymore in the standard version. However, this option can still be
189	! switched on by replacing the following lines with those below, that are
190	! currently commented out.
191	! Note that one change is also necessary in gridcheck.f,
192	! and three changes in verttransform.f
193	!
194	! *** NOTE -- the doubled vertical resolution has not been tested in FLEXPART_WRF
195	!*******************************************************************************
196	uvwzlev(ix,jy,1)=0.0
197	do kz=2,nuvz
198	uvwzlev(ix,jy,kz)=uvzlev(kz)
199	enddo
200	! Switch on following lines to use doubled vertical resolution
201	! Switch off the three lines above.
202	!
203	! *** NOTE -- the doubled vertical resolution has not been tested in FLEXPART_WRF
204	!*************************************************************
205	!22 uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
206	! do 23 kz=2,nwz
207	!23 uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
208	! End doubled vertical resolution
209
210	! pinmconv=(h2-h1)/(p2-p1)
211	!
212	! in flexpart_ecmwf, pinmconv is used to convert etadot to w
213	! in FLEXPART_WRF, vertical velocity is already m/s, so pinmconv=1.0
214	!
215	! pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/
216	! + ((aknew(2)+bknew(2)*psn(ix,jy,1,n,l))-
217	! + (aknew(1)+bknew(1)*psn(ix,jy,1,n,l)))
218	if (wind_option.eq.0) then
219	pinmconv(1)=1.0
220	do kz=2,nz-1
221	pinmconv(kz)=1.0
222	enddo
223	pinmconv(nz)=1.0
224	elseif (wind_option.eq.1) then
225	! pinmconv(1)=(uvzlev(2)-uvzlev(1)) &
226	! /(eta_u_wrf(1)-1.)
227	pinmconv(1)=(wzlev(2)-0.) &
228	/(eta_w_wrf(2)-1.)
229	do kz=2,nz-1
230	! pinmconv(kz)=(uvzlev(kz)-uvzlev(kz-1)) &
231	! /(eta_u_wrf(kz)-eta_u_wrf(kz-1))
232	pinmconv(kz)=(wzlev(kz+1)-wzlev(kz-1)) &
233	/(eta_w_wrf(kz+1)-eta_w_wrf(kz-1))
234	enddo
235	pinmconv(nwz)=pinmconv(nwz-1)
236	endif
237
238
239	! print*,'etap 3',ix,jy
240
241	! Levels, where u,v,t and q are given
242	!************************************
243
244	uun(ix,jy,1,n,l)=uuhn(ix,jy,1,l)
245	vvn(ix,jy,1,n,l)=vvhn(ix,jy,1,l)
246	divn(ix,jy,1,l)=divhn(ix,jy,1,l)
247	ttn(ix,jy,1,n,l)=tthn(ix,jy,1,n,l)
248	qvn(ix,jy,1,n,l)=qvhn(ix,jy,1,n,l)
249	pvn(ix,jy,1,n,l)=pvhn(ix,jy,1,l)
250	rhon(ix,jy,1,n,l)=rhoh(1)
251	uun(ix,jy,nz,n,l)=uuhn(ix,jy,nuvz,l)
252	vvn(ix,jy,nz,n,l)=vvhn(ix,jy,nuvz,l)
253	ttn(ix,jy,nz,n,l)=tthn(ix,jy,nuvz,n,l)
254	qvn(ix,jy,nz,n,l)=qvhn(ix,jy,nuvz,n,l)
255	pvn(ix,jy,nz,n,l)=pvhn(ix,jy,nuvz,l)
256	rhon(ix,jy,nz,n,l)=rhoh(nuvz)
257	tken(ix,jy,1,n,l)=tkehn(ix,jy,1,n,l)
258	tken(ix,jy,nz,n,l)=tkehn(ix,jy,nuvz,n,l)
259	pttn(ix,jy,1,n,l)=ptthn(ix,jy,1,n,l)
260	pttn(ix,jy,nz,n,l)=ptthn(ix,jy,nuvz,n,l)
261
262
263	! print*,'etap 3.5',ix,jy
264	kmin=2
265	do iz=2,nz-1
266	do kz=kmin,nuvz
267	if(heightmid(iz).gt.uvzlev(nuvz)) then
268	divn(ix,jy,iz,l)=divn(ix,jy,nz,l)
269	goto 230
270	endif
271	!! print*,'etap 3.7',kz,iz,heightmid(iz),uvzlev(kz-1),uvzlev(kz)
272	if ((heightmid(iz).gt.uvzlev(kz-1)).and. &
273	(heightmid(iz).le.uvzlev(kz))) then
274	dz1=heightmid(iz)-uvzlev(kz-1)
275	dz2=uvzlev(kz)-heightmid(iz)
276	dz=dz1+dz2
277	divn(ix,jy,iz,l)=(divhn(ix,jy,kz-1,l)dz2+divhn(ix,jy,kz,l)dz1)/dz
278	kmin=kz
279	goto 230
280	endif
281	end do
282	230 continue
283	end do
284
285	! print*,'etap 4',ix,jy
286
287	kmin=2
288	do iz=2,nz-1
289	do kz=kmin,nuvz
290	if(height(iz).gt.uvzlev(nuvz)) then
291	uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l)
292	vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l)
293	ttn(ix,jy,iz,n,l)=ttn(ix,jy,nz,n,l)
294	qvn(ix,jy,iz,n,l)=qvn(ix,jy,nz,n,l)
295	pvn(ix,jy,iz,n,l)=pvn(ix,jy,nz,n,l)
296	rhon(ix,jy,iz,n,l)=rhon(ix,jy,nz,n,l)
297	tken(ix,jy,iz,n,l)=tken(ix,jy,nz,n,l)
298	pttn(ix,jy,iz,n,l)=pttn(ix,jy,nz,n,l)
299	goto 30
300	endif
301	if ((height(iz).gt.uvzlev(kz-1)).and. &
302	(height(iz).le.uvzlev(kz))) then
303	dz1=height(iz)-uvzlev(kz-1)
304	dz2=uvzlev(kz)-height(iz)
305	dz=dz1+dz2
306	uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+ &
307	uuhn(ix,jy,kz,l)*dz1)/dz
308	vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+ &
309	vvhn(ix,jy,kz,l)*dz1)/dz
310	ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2+ &
311	tthn(ix,jy,kz,n,l)*dz1)/dz
312	qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2+ &
313	qvhn(ix,jy,kz,n,l)*dz1)/dz
314	pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+ &
315	pvhn(ix,jy,kz,l)*dz1)/dz
316	rhon(ix,jy,iz,n,l)=(rhoh(kz-1)dz2+rhoh(kz)dz1)/dz
317	tken(ix,jy,iz,n,l)=(tkehn(ix,jy,kz-1,n,l)*dz2+ &
318	tkehn(ix,jy,kz,n,l)*dz1)/dz
319	pttn(ix,jy,iz,n,l)=(ptthn(ix,jy,kz-1,n,l)*dz2+ &
320	ptthn(ix,jy,kz,n,l)*dz1)/dz
321
322
323	kmin=kz
324	goto 30
325	endif
326	end do
327	30 continue
328	end do
329
330	! print*,'continue to ww in nests'
331	! Levels, where w is given
332	!*************************
333
334	if (method_w_terrain_correction .eq. 20) then
335	! apply w correction assuming that the WRF w is "absolute w";
336	! apply it here to wwh; set wwh=0 at iz=1
337	ix1 = max( ix-1, 0 )
338	jy1 = max( jy-1, 0 )
339	ixp = min( ix+1, nxn(l)-1 )
340	jyp = min( jy+1, nyn(l)-1 )
341	if (wind_option.eq.0) then
342	dzdx=(oron(ixp,jy,l)-oron(ix1,jy,l))/(dxn(l)(ixp-ix1)m_xn(ix,jy,1,l))
343	dzdy=(oron(ix,jyp,l)-oron(ix,jy1,l))/(dyn(l)(jyp-jy1)m_yn(ix,jy,1,l))
344
345	do kz = 1, nwz-1
346	! wwhn_svaa(kz) = wwhn(ix,jy,kz,l)
347	wwhn(ix,jy,kz,l) = wwhn(ix,jy,kz,l) &
348	- (uuhn(ix,jy,kz,l)dzdx + vvhn(ix,jy,kz,l)dzdy)
349	! if (kz .eq. 1) wwhn(ix,jy,kz,l) = 0.0
350	end do
351	elseif (wind_option.ge.1) then
352	do kz = 2, nwz-1
353	! wwhn_svaa(kz) = wwhn(ix,jy,kz,l)
354	dzdx=(zzhn(ixp,jy,kz+add_sfc_level,n,l) - zzhn(ix1,jy,kz+add_sfc_level,n,l) &
355	-zzhn(ixp,jy,1,n,l) &
356	+zzhn(ix1,jy,1,n,l)) &
357	/(dxn(l)(ixp-ix1)m_xn(ix,jy,1,l))
358	dzdy=(zzhn(ix,jyp,kz+add_sfc_level,n,l) - zzhn(ix,jy1,kz+add_sfc_level,n,l) &
359	-zzhn(ix,jyp,1,n,l) &
360	+zzhn(ix,jy1,1,n,l)) &
361	/(dyn(l)(jyp-jy1)m_yn(ix,jy,1,l))
362
363	dzdx=(zzhn(ixp,jy,kz+add_sfc_level,n,l) - zzhn(ix1,jy,kz+add_sfc_level,n,l) &
364	-zzhn(ixp,jy,1,n,l)+zzhn(ix1,jy,1,n,l))/(dxn(l)(ixp-ix1)m_xn(ix,jy,1,l))
365	dzdy=(zzhn(ix,jyp,kz+add_sfc_level,n,l) - zzhn(ix,jy1,kz+add_sfc_level,n,l) &
366	-zzhn(ix,jyp,1,n,l)+zzhn(ix,jy1,1,n,l))/(dyn(l)(jyp-jy1)m_yn(ix,jy,1,l))
367	u=0.5*(uuhn(ix,jy,kz+add_sfc_level,l)+uuhn(ix,jy,kz-1+add_sfc_level,l))
368	v=0.5*(vvhn(ix,jy,kz+add_sfc_level,l)+vvhn(ix,jy,kz-1+add_sfc_level,l))
369
370	wwhn(ix,jy,kz,l) = wwhn(ix,jy,kz,l)*pinmconv(kz) &
371	! + (uuhn(ix,jy,kz,l)dzdx + vvhn(ix,jy,kz,l)dzdy) ! variation of geopot on sigma is necessary
372	+ (udzdx + vdzdy) ! variation of geopot on sigma is necessary
373	if (kz .eq. 1) wwhn(ix,jy,kz,l) = wwhn(ix,jy,kz,l)*pinmconv(kz)
374
375	! if (kz .eq. 1) wwhn(ix,jy,kz,l) = 0.0
376	end do
377	endif
378	if (wind_option.eq.-1) then
379	!! ww(ix,jy,1,n)=wwh(ix,jy,1)
380	wwn(ix,jy,1,n,l)=0.
381	do iz=2,nz
382	wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz-1,n,l)-(height(iz)-height(iz-1))* &
383	divn(ix,jy,iz-1,l)
384	enddo
385	else
386	! print*,'converting ww in nest'
387	wwn(ix,jy,1,n,l)=wwhn(ix,jy,1,l)
388	wwn(ix,jy,nz,n,l)=wwhn(ix,jy,nwz,l)
389	kmin=2
390	do iz=2,nz
391	do kz=kmin,nwz
392	if ((height(iz).gt.wzlev(kz-1)).and. &
393	(height(iz).le.wzlev(kz))) then
394	dz1=height(iz)-wzlev(kz-1)
395	dz2=wzlev(kz)-height(iz)
396	dz=dz1+dz2
397	wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*dz2+ &
398	wwhn(ix,jy,kz,l)*dz1)/dz
399	kmin=kz
400	goto 40
401	endif
402	end do
403	40 continue
404	end do
405	endif
406
407	! if (method_w_terrain_correction .eq. 20) then
408	! do kz = 1, nwz
409	! wwhn(ix,jy,kz,l) = wwhn_svaa(kz)
410	! end do
411	! end if
412	end if
413
414	! Compute density gradients at intermediate levels
415	!*************************************************
416
417	drhodzn(ix,jy,1,n,l)=(rhon(ix,jy,2,n,l)-rhon(ix,jy,1,n,l))/ &
418	(height(2)-height(1))
419	do kz=2,nz-1
420	drhodzn(ix,jy,kz,n,l)=(rhon(ix,jy,kz+1,n,l)- &
421	rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
422	end do
423	drhodzn(ix,jy,nz,n,l)=drhodzn(ix,jy,nz-1,n,l)
424
425	end do
426	end do
427	!!!$OMP END DO
428	!!!$OMP END PARALLEL
429
430	! print*,'end of ww, now clouds, nests'
431	!****************************************************************
432	! Compute slope of eta levels in windward direction and resulting
433	! vertical wind correction
434	!
435	! See notes in verttransform.f about the w correction done here.
436	!****************************************************************
437	!write (,) 'initializing clouds, n:',n,nymin1,nxmin1,nz^M
438	! create a cloud and rainout/washout field, clouds occur where rh>80%^M
439	! total cloudheight is stored at level 0^M
440
441	do 100 jy=0,nyn(l)-1
442	do 100 ix=0,nxn(l)-1
443	! rain_cloud_above=0
444	lsp=lsprecn(ix,jy,1,n,l)
445	convp=convprecn(ix,jy,1,n,l)
446
447	! cloudsh(ix,jy,n)=0
448
449	prec=lsp+convp
450	if (lsp.gt.convp) then ! prectype='lsp'
451	lconvectprec = .false.
452	else ! prectype='cp '
453	lconvectprec = .true.
454	endif
455	rhmin = 0.90 ! standard condition for presence of clouds
456
457	!CPS note that original by Sabine Eckhart was 80%
458	!CPS however, for T<-20 C we consider saturation over ice
459	!CPS so I think 90% should be enough
460
461	icloudbotn(ix,jy,n,l)=icmv
462	icloudtop=icmv ! this is just a local variable
463	98 do kz=1,nz
464	pressure=rhon(ix,jy,kz,n,l)r_airttn(ix,jy,kz,n,l)
465	rh=qvn(ix,jy,kz,n,l)/f_qvsat(pressure,ttn(ix,jy,kz,n,l))
466	!cps if (prec.gt.0.01) print*,'relhum',prec,kz,rh,height(kz)
467	if (rh .gt. rhmin) then
468	if (icloudbotn(ix,jy,n,l) .eq. icmv) then
469	icloudbotn(ix,jy,n,l)=nint(height(kz))
470	endif
471	icloudtop=nint(height(kz)) ! use int to save memory
472	endif
473	enddo
474
475	!CPS try to get a cloud thicker than 50 m
476	!CPS if there is at least .01 mm/h - changed to 0.002 and put into
477	!CPS parameter precpmin
478	if ((icloudbotn(ix,jy,n,l) .eq. icmv .or. &
479	icloudtop-icloudbotn(ix,jy,n,l) .lt. 50) .and. &
480	prec .gt. precmin) then
481	rhmin = rhmin - 0.05
482	if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum.
483	endif
484	!CPS implement a rough fix for badly represented convection
485	!CPS is based on looking at a limited set of comparison data
486	if (lconvectprec .and. icloudtop .lt. 6000 .and. &
487	prec .gt. precmin) then
488	if (convp .lt. 0.1) then
489	icloudbotn(ix,jy,n,l) = 500
490	icloudtop = 8000
491	else
492	icloudbotn(ix,jy,n,l) = 0
493	icloudtop = 10000
494	endif
495	endif
496	if (icloudtop .ne. icmv) then
497	icloudthckn(ix,jy,n,l) = icloudtop-icloudbotn(ix,jy,n,l)
498	else
499	icloudthckn(ix,jy,n,l) = icmv
500	endif
501	!CPS get rid of too thin clouds
502	if (icloudthckn(ix,jy,n,l) .lt. 50) then
503	icloudbotn(ix,jy,n,l)=icmv
504	icloudthckn(ix,jy,n,l)=icmv
505	endif
506
507	100 continue
508	enddo ! nests
509
510	return
511	end

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format