Context Navigation

verttransform.f90 @ 20

Last change on this file since 20 was 20, checked in by igpis, 11 years ago
move version 9.1.8 form branches to trunk. Contributions from HSO, saeck, pesei, NIK, RT, XKF, IP and others
File size: 22.8 KB

Line
1	!**********************************************************************
2	! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
3	! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
4	! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
5	! *
6	! This file is part of FLEXPART. *
7	! *
8	! FLEXPART is free software: you can redistribute it and/or modify *
9	! it under the terms of the GNU General Public License as published by*
10	! the Free Software Foundation, either version 3 of the License, or *
11	! (at your option) any later version. *
12	! *
13	! FLEXPART is distributed in the hope that it will be useful, *
14	! but WITHOUT ANY WARRANTY; without even the implied warranty of *
15	! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
16	! GNU General Public License for more details. *
17	! *
18	! You should have received a copy of the GNU General Public License *
19	! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
20	!**********************************************************************
21
22	subroutine verttransform(n,uuh,vvh,wwh,pvh)
23	! i i i i i
24	!*****************************************************************************
25	! *
26	! This subroutine transforms temperature, dew point temperature and *
27	! wind components from eta to meter coordinates. *
28	! The vertical wind component is transformed from Pa/s to m/s using *
29	! the conversion factor pinmconv. *
30	! In addition, this routine calculates vertical density gradients *
31	! needed for the parameterization of the turbulent velocities. *
32	! *
33	! Author: A. Stohl, G. Wotawa *
34	! *
35	! 12 August 1996 *
36	! Update: 16 January 1998 *
37	! *
38	! Major update: 17 February 1999 *
39	! by G. Wotawa *
40	! *
41	! - Vertical levels for u, v and w are put together *
42	! - Slope correction for vertical velocity: Modification of calculation *
43	! procedure *
44	! *
45	!*****************************************************************************
46	! Changes, Bernd C. Krueger, Feb. 2001:
47	! Variables tth and qvh (on eta coordinates) from common block
48	!*****************************************************************************
49	! Sabine Eckhardt, March 2007
50	! added the variable cloud for use with scavenging - descr. in com_mod
51	! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification
52	! note that also other subroutines are affected by the fix
53	!*****************************************************************************
54	! *
55	! Variables: *
56	! nx,ny,nz field dimensions in x,y and z direction *
57	! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition *
58	! uu(0:nxmax,0:nymax,nzmax,2) wind components in x-direction [m/s] *
59	! vv(0:nxmax,0:nymax,nzmax,2) wind components in y-direction [m/s] *
60	! ww(0:nxmax,0:nymax,nzmax,2) wind components in z-direction [deltaeta/s]*
61	! tt(0:nxmax,0:nymax,nzmax,2) temperature [K] *
62	! pv(0:nxmax,0:nymax,nzmax,2) potential voriticity (pvu) *
63	! ps(0:nxmax,0:nymax,2) surface pressure [Pa] *
64	! *
65	!*****************************************************************************
66
67	use par_mod
68	use com_mod
69	use cmapf_mod, only: cc2gll
70
71	implicit none
72
73	integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
74	integer :: rain_cloud_above,kz_inv !SE
75	integer icloudtop !PS
76	real :: f_qvsat,pressure
77	!real :: rh,lsp,convp
78	real :: rh,lsp,convp,prec,rhmin
79	real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
80	real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
81	real :: xlon,ylat,xlonr,dzdx,dzdy
82	real :: dzdx1,dzdx2,dzdy1,dzdy2, precmin
83	real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
84	real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
85	real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
86	real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
87	real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
88	real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
89	logical lconvectprec
90	real,parameter :: const=r_air/ga
91	parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics
92
93	logical :: init = .true.
94
95
96	!*************************************************************************
97	! If verttransform is called the first time, initialize heights of the *
98	! z levels in meter. The heights are the heights of model levels, where *
99	! u,v,T and qv are given, and of the interfaces, where w is given. So, *
100	! the vertical resolution in the z system is doubled. As reference point,*
101	! the lower left corner of the grid is used. *
102	! Unlike in the eta system, no difference between heights for u,v and *
103	! heights for w exists. *
104	!*************************************************************************
105
106
107	! do 897 kz=1,nuvz
108	! write (,) 'akz: ',akz(kz),'bkz',bkz(kz)
109	!897 continue
110
111	if (init) then
112
113	! Search for a point with high surface pressure (i.e. not above significant topography)
114	! Then, use this point to construct a reference z profile, to be used at all times
115	!*****************************************************************************
116
117	do jy=0,nymin1
118	do ix=0,nxmin1
119	if (ps(ix,jy,1,n).gt.100000.) then
120	ixm=ix
121	jym=jy
122	goto 3
123	endif
124	end do
125	end do
126	3 continue
127
128
129	tvold=tt2(ixm,jym,1,n)(1.+0.378ew(td2(ixm,jym,1,n))/ &
130	ps(ixm,jym,1,n))
131	pold=ps(ixm,jym,1,n)
132	height(1)=0.
133
134	do kz=2,nuvz
135	pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
136	tv=tth(ixm,jym,kz,n)(1.+0.608qvh(ixm,jym,kz,n))
137
138
139	! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
140	! upon the transformation to z levels. In order to save computer memory, this is
141	! not done anymore in the standard version. However, this option can still be
142	! switched on by replacing the following lines with those below, that are
143	! currently commented out.
144	! Note that two more changes are necessary in this subroutine below.
145	! One change is also necessary in gridcheck.f, and another one in verttransform_nests.
146	!*****************************************************************************
147
148	if (abs(tv-tvold).gt.0.2) then
149	height(kz)= &
150	height(kz-1)+constlog(pold/pint) &
151	(tv-tvold)/log(tv/tvold)
152	else
153	height(kz)=height(kz-1)+ &
154	constlog(pold/pint)tv
155	endif
156
157	! Switch on following lines to use doubled vertical resolution
158	!*************************************************************
159	! if (abs(tv-tvold).gt.0.2) then
160	! height((kz-1)*2)=
161	! + height(max((kz-2)2,1))+constlog(pold/pint)*
162	! + (tv-tvold)/log(tv/tvold)
163	! else
164	! height((kz-1)2)=height(max((kz-2)2,1))+
165	! + constlog(pold/pint)tv
166	! endif
167	! End doubled vertical resolution
168
169	tvold=tv
170	pold=pint
171	end do
172
173
174	! Switch on following lines to use doubled vertical resolution
175	!*************************************************************
176	! do 7 kz=3,nz-1,2
177	! height(kz)=0.5*(height(kz-1)+height(kz+1))
178	! height(nz)=height(nz-1)+height(nz-1)-height(nz-2)
179	! End doubled vertical resolution
180
181
182	! Determine highest levels that can be within PBL
183	!************************************************
184
185	do kz=1,nz
186	if (height(kz).gt.hmixmax) then
187	nmixz=kz
188	goto 9
189	endif
190	end do
191	9 continue
192
193	! Do not repeat initialization of the Cartesian z grid
194	!*****************************************************
195
196	init=.false.
197
198	endif
199
200
201	! Loop over the whole grid
202	!*************************
203
204	do jy=0,nymin1
205	do ix=0,nxmin1
206	tvold=tt2(ix,jy,1,n)(1.+0.378ew(td2(ix,jy,1,n))/ &
207	ps(ix,jy,1,n))
208	pold=ps(ix,jy,1,n)
209	uvzlev(1)=0.
210	wzlev(1)=0.
211	rhoh(1)=pold/(r_air*tvold)
212
213
214	! Compute heights of eta levels
215	!******************************
216
217	do kz=2,nuvz
218	pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)
219	tv=tth(ix,jy,kz,n)(1.+0.608qvh(ix,jy,kz,n))
220	rhoh(kz)=pint/(r_air*tv)
221
222	if (abs(tv-tvold).gt.0.2) then
223	uvzlev(kz)=uvzlev(kz-1)+constlog(pold/pint) &
224	(tv-tvold)/log(tv/tvold)
225	else
226	uvzlev(kz)=uvzlev(kz-1)+constlog(pold/pint)tv
227	endif
228
229	tvold=tv
230	pold=pint
231	end do
232
233
234	do kz=2,nwz-1
235	wzlev(kz)=(uvzlev(kz+1)+uvzlev(kz))/2.
236	end do
237	wzlev(nwz)=wzlev(nwz-1)+ &
238	uvzlev(nuvz)-uvzlev(nuvz-1)
239
240	uvwzlev(ix,jy,1)=0.0
241	do kz=2,nuvz
242	uvwzlev(ix,jy,kz)=uvzlev(kz)
243	end do
244
245	! Switch on following lines to use doubled vertical resolution
246	! Switch off the three lines above.
247	!*************************************************************
248	!22 uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
249	! do 23 kz=2,nwz
250	!23 uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
251	! End doubled vertical resolution
252
253	! pinmconv=(h2-h1)/(p2-p1)
254
255	pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
256	((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
257	(aknew(1)+bknew(1)*ps(ix,jy,1,n)))
258	do kz=2,nz-1
259	pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
260	((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
261	(aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
262	end do
263	pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
264	((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
265	(aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
266
267	! Levels, where u,v,t and q are given
268	!************************************
269
270	uu(ix,jy,1,n)=uuh(ix,jy,1)
271	vv(ix,jy,1,n)=vvh(ix,jy,1)
272	tt(ix,jy,1,n)=tth(ix,jy,1,n)
273	qv(ix,jy,1,n)=qvh(ix,jy,1,n)
274	pv(ix,jy,1,n)=pvh(ix,jy,1)
275	rho(ix,jy,1,n)=rhoh(1)
276	uu(ix,jy,nz,n)=uuh(ix,jy,nuvz)
277	vv(ix,jy,nz,n)=vvh(ix,jy,nuvz)
278	tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
279	qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
280	pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
281	rho(ix,jy,nz,n)=rhoh(nuvz)
282	kmin=2
283	do iz=2,nz-1
284	do kz=kmin,nuvz
285	if(height(iz).gt.uvzlev(nuvz)) then
286	uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
287	vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
288	tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
289	qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
290	pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
291	rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
292	goto 30
293	endif
294	if ((height(iz).gt.uvzlev(kz-1)).and. &
295	(height(iz).le.uvzlev(kz))) then
296	dz1=height(iz)-uvzlev(kz-1)
297	dz2=uvzlev(kz)-height(iz)
298	dz=dz1+dz2
299	uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)dz2+uuh(ix,jy,kz)dz1)/dz
300	vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)dz2+vvh(ix,jy,kz)dz1)/dz
301	tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
302	+tth(ix,jy,kz,n)*dz1)/dz
303	qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
304	+qvh(ix,jy,kz,n)*dz1)/dz
305	pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)dz2+pvh(ix,jy,kz)dz1)/dz
306	rho(ix,jy,iz,n)=(rhoh(kz-1)dz2+rhoh(kz)dz1)/dz
307	kmin=kz
308	goto 30
309	endif
310	end do
311	30 continue
312	end do
313
314
315	! Levels, where w is given
316	!*************************
317
318	ww(ix,jy,1,n)=wwh(ix,jy,1)*pinmconv(1)
319	ww(ix,jy,nz,n)=wwh(ix,jy,nwz)*pinmconv(nz)
320	kmin=2
321	do iz=2,nz
322	do kz=kmin,nwz
323	if ((height(iz).gt.wzlev(kz-1)).and. &
324	(height(iz).le.wzlev(kz))) then
325	dz1=height(iz)-wzlev(kz-1)
326	dz2=wzlev(kz)-height(iz)
327	dz=dz1+dz2
328	ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)pinmconv(kz-1)dz2 &
329	+wwh(ix,jy,kz)pinmconv(kz)dz1)/dz
330	kmin=kz
331	goto 40
332	endif
333	end do
334	40 continue
335	end do
336
337	! Compute density gradients at intermediate levels
338	!*************************************************
339
340	drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
341	(height(2)-height(1))
342	do kz=2,nz-1
343	drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
344	(height(kz+1)-height(kz-1))
345	end do
346	drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
347
348	end do
349	end do
350
351
352	!****************************************************************
353	! Compute slope of eta levels in windward direction and resulting
354	! vertical wind correction
355	!****************************************************************
356
357	do jy=1,ny-2
358	do ix=1,nx-2
359
360	kmin=2
361	do iz=2,nz-1
362
363	ui=uu(ix,jy,iz,n)dxconst/cos((real(jy)dy+ylat0)*pi180)
364	vi=vv(ix,jy,iz,n)*dyconst
365
366	do kz=kmin,nz
367	if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
368	(height(iz).le.uvwzlev(ix,jy,kz))) then
369	dz1=height(iz)-uvwzlev(ix,jy,kz-1)
370	dz2=uvwzlev(ix,jy,kz)-height(iz)
371	dz=dz1+dz2
372	kl=kz-1
373	klp=kz
374	kmin=kz
375	goto 47
376	endif
377	end do
378
379	47 ix1=ix-1
380	jy1=jy-1
381	ixp=ix+1
382	jyp=jy+1
383
384	dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
385	dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
386	dzdx=(dzdx1dz2+dzdx2dz1)/dz
387
388	dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
389	dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
390	dzdy=(dzdy1dz2+dzdy2dz1)/dz
391
392	ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdxui+dzdyvi)
393
394	end do
395
396	end do
397	end do
398
399
400	! If north pole is in the domain, calculate wind velocities in polar
401	! stereographic coordinates
402	!*******************************************************************
403
404	if (nglobal) then
405	do jy=int(switchnorthg)-2,nymin1
406	ylat=ylat0+real(jy)*dy
407	do ix=0,nxmin1
408	xlon=xlon0+real(ix)*dx
409	do iz=1,nz
410	call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
411	vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
412	vvpol(ix,jy,iz,n))
413	end do
414	end do
415	end do
416
417
418	do iz=1,nz
419
420	! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
421	!
422	! AMSnauffer Nov 18 2004 Added check for case vv=0
423	!
424	xlon=xlon0+real(nx/2-1)*dx
425	xlonr=xlon*pi/180.
426	ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
427	vv(nx/2-1,nymin1,iz,n)**2)
428	if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
429	ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
430	vv(nx/2-1,nymin1,iz,n))-xlonr
431	else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
432	ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
433	vv(nx/2-1,nymin1,iz,n))-xlonr
434	else
435	ddpol=pi/2-xlonr
436	endif
437	if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
438	if(ddpol.gt.2.0pi) ddpol=ddpol-2.0pi
439
440	! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
441	xlon=180.0
442	xlonr=xlon*pi/180.
443	ylat=90.0
444	uuaux=-ffpol*sin(xlonr+ddpol)
445	vvaux=-ffpol*cos(xlonr+ddpol)
446	call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
447	vvpolaux)
448
449	jy=nymin1
450	do ix=0,nxmin1
451	uupol(ix,jy,iz,n)=uupolaux
452	vvpol(ix,jy,iz,n)=vvpolaux
453	end do
454	end do
455
456
457	! Fix: Set W at pole to the zonally averaged W of the next equator-
458	! ward parallel of latitude
459
460	do iz=1,nz
461	wdummy=0.
462	jy=ny-2
463	do ix=0,nxmin1
464	wdummy=wdummy+ww(ix,jy,iz,n)
465	end do
466	wdummy=wdummy/real(nx)
467	jy=nymin1
468	do ix=0,nxmin1
469	ww(ix,jy,iz,n)=wdummy
470	end do
471	end do
472
473	endif
474
475
476	! If south pole is in the domain, calculate wind velocities in polar
477	! stereographic coordinates
478	!*******************************************************************
479
480	if (sglobal) then
481	do jy=0,int(switchsouthg)+3
482	ylat=ylat0+real(jy)*dy
483	do ix=0,nxmin1
484	xlon=xlon0+real(ix)*dx
485	do iz=1,nz
486	call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
487	vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
488	vvpol(ix,jy,iz,n))
489	end do
490	end do
491	end do
492
493	do iz=1,nz
494
495	! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
496	!
497	! AMSnauffer Nov 18 2004 Added check for case vv=0
498	!
499	xlon=xlon0+real(nx/2-1)*dx
500	xlonr=xlon*pi/180.
501	ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
502	vv(nx/2-1,0,iz,n)**2)
503	if (vv(nx/2-1,0,iz,n).lt.0.) then
504	ddpol=atan(uu(nx/2-1,0,iz,n)/ &
505	vv(nx/2-1,0,iz,n))+xlonr
506	else if (vv(nx/2-1,0,iz,n).gt.0.) then
507	ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
508	vv(nx/2-1,0,iz,n))+xlonr
509	else
510	ddpol=pi/2-xlonr
511	endif
512	if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
513	if(ddpol.gt.2.0pi) ddpol=ddpol-2.0pi
514
515	! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
516	xlon=180.0
517	xlonr=xlon*pi/180.
518	ylat=-90.0
519	uuaux=+ffpol*sin(xlonr-ddpol)
520	vvaux=-ffpol*cos(xlonr-ddpol)
521	call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
522	vvpolaux)
523
524	jy=0
525	do ix=0,nxmin1
526	uupol(ix,jy,iz,n)=uupolaux
527	vvpol(ix,jy,iz,n)=vvpolaux
528	end do
529	end do
530
531
532	! Fix: Set W at pole to the zonally averaged W of the next equator-
533	! ward parallel of latitude
534
535	do iz=1,nz
536	wdummy=0.
537	jy=1
538	do ix=0,nxmin1
539	wdummy=wdummy+ww(ix,jy,iz,n)
540	end do
541	wdummy=wdummy/real(nx)
542	jy=0
543	do ix=0,nxmin1
544	ww(ix,jy,iz,n)=wdummy
545	end do
546	end do
547	endif
548
549
550	!write (,) 'initializing clouds, n:',n,nymin1,nxmin1,nz
551	! create a cloud and rainout/washout field, clouds occur where rh>80%
552	! total cloudheight is stored at level 0
553
554
555
556	do jy=0,nymin1
557	do ix=0,nxmin1
558
559
560
561	! rain_cloud_above=0
562	! lsp=lsprec(ix,jy,1,n)
563	! convp=convprec(ix,jy,1,n)
564	! cloudsh(ix,jy,n)=0
565	! do kz_inv=1,nz-1
566	! kz=nz-kz_inv+1
567	! pressure=rho(ix,jy,kz,n)r_airtt(ix,jy,kz,n)
568	! rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
569	! clouds(ix,jy,kz,n)=0
570	! if (rh.gt.0.8) then ! in cloud
571	! if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
572	! rain_cloud_above=1
573	! cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
574	! height(kz)-height(kz-1)
575	! if (lsp.ge.convp) then
576	! clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
577	! else
578	! clouds(ix,jy,kz,n)=2 ! convp dominated rainout
579	! endif
580	! else ! no precipitation
581	! clouds(ix,jy,kz,n)=1 ! cloud
582	! endif
583	! else ! no cloud
584	! if (rain_cloud_above.eq.1) then ! scavenging
585	! if (lsp.ge.convp) then
586	! clouds(ix,jy,kz,n)=5 ! lsp dominated washout
587	! else
588	! clouds(ix,jy,kz,n)=4 ! convp dominated washout
589	! endif
590	! endif
591	! endif
592	! end do
593
594
595	! PS 3012
596
597	lsp=lsprec(ix,jy,1,n)
598	convp=convprec(ix,jy,1,n)
599	prec=lsp+convp
600	if (lsp.gt.convp) then ! prectype='lsp'
601	lconvectprec = .false.
602	else ! prectype='cp '
603	lconvectprec = .true.
604	endif
605	rhmin = 0.90 ! standard condition for presence of clouds
606	!PS note that original by Sabine Eckhart was 80%
607	!PS however, for T<-20 C we consider saturation over ice
608	!PS so I think 90% should be enough
609	icloudbot(ix,jy,n)=icmv
610	icloudtop=icmv ! this is just a local variable
611	98 do kz=1,nz
612	pressure=rho(ix,jy,kz,n)r_airtt(ix,jy,kz,n)
613	rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
614	!ps if (prec.gt.0.01) print*,'relhum',prec,kz,rh,height(kz)
615	if (rh .gt. rhmin) then
616	if (icloudbot(ix,jy,n) .eq. icmv) then
617	icloudbot(ix,jy,n)=nint(height(kz))
618	endif
619	icloudtop=nint(height(kz)) ! use int to save memory
620	endif
621	enddo
622
623	!PS try to get a cloud thicker than 50 m
624	!PS if there is at least .01 mm/h - changed to 0.002 and put into
625	!PS parameter precpmin
626	if ((icloudbot(ix,jy,n) .eq. icmv .or. &
627	icloudtop-icloudbot(ix,jy,n) .lt. 50) .and. &
628	prec .gt. precmin) then
629	rhmin = rhmin - 0.05
630	if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum.
631	endif
632	!PS implement a rough fix for badly represented convection
633	!PS is based on looking at a limited set of comparison data
634	if (lconvectprec .and. icloudtop .lt. 6000 .and. &
635	prec .gt. precmin) then
636
637	if (convp .lt. 0.1) then
638	icloudbot(ix,jy,n) = 500
639	icloudtop = 8000
640	else
641	icloudbot(ix,jy,n) = 0
642	icloudtop = 10000
643	endif
644	endif
645	if (icloudtop .ne. icmv) then
646	icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n)
647	else
648	icloudthck(ix,jy,n) = icmv
649	endif
650	!PS get rid of too thin clouds
651	if (icloudthck(ix,jy,n) .lt. 50) then
652	icloudbot(ix,jy,n)=icmv
653	icloudthck(ix,jy,n)=icmv
654	endif
655
656
657	end do
658	end do
659
660	!do 102 kz=1,nuvz
661	!write(an,'(i02)') kz+10
662	!write(,) nuvz,nymin1,nxmin1,'--',an,'--'
663	!open(4,file='/nilu_wrk2/sec/cloudtest/cloud'//an,form='formatted')
664	!do 101 jy=0,nymin1
665	! write(4,*) (clouds(ix,jy,kz,n),ix=1,nxmin1)
666	!101 continue
667	! close(4)
668	!102 continue
669
670	! open(4,file='/nilu_wrk2/sec/cloudtest/height',form='formatted')
671	! do 103 jy=0,nymin1
672	! write (4,*)
673	!+ (height(kz),kz=1,nuvz)
674	!103 continue
675	! close(4)
676
677	!open(4,file='/nilu_wrk2/sec/cloudtest/p',form='formatted')
678	! do 104 jy=0,nymin1
679	! write (4,*)
680	!+ (r_airtt(ix,jy,1,n)rho(ix,jy,1,n),ix=1,nxmin1)
681	!104 continue
682	! close(4)
683
684
685	end subroutine verttransform

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

Download in other formats:

Original Format