Context Navigation

verttransform.f90 @ f9ce123

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie

Last change on this file since f9ce123 was 8a65cb0, checked in by Espen Sollum ATMOS <espen@…>, 9 years ago
Added code, makefile for dev branch
Property mode set to `100644`
File size: 22.2 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	!*****************************************************************************
52	! *
53	! Variables: *
54	! nx,ny,nz field dimensions in x,y and z direction *
55	! clouds(0:nxmax,0:nymax,0:nzmax,numwfmem) cloud field for wet deposition *
56	! uu(0:nxmax,0:nymax,nzmax,numwfmem) wind components in x-direction [m/s]*
57	! vv(0:nxmax,0:nymax,nzmax,numwfmem) wind components in y-direction [m/s]*
58	! ww(0:nxmax,0:nymax,nzmax,numwfmem) wind components in z-direction *
59	! [deltaeta/s] *
60	! tt(0:nxmax,0:nymax,nzmax,numwfmem) temperature [K] *
61	! pv(0:nxmax,0:nymax,nzmax,numwfmem) potential voriticity (pvu) *
62	! ps(0:nxmax,0:nymax,numwfmem) surface pressure [Pa] *
63	! *
64	!*****************************************************************************
65
66	use par_mod
67	use com_mod
68	use cmapf_mod, only: cc2gll
69
70	implicit none
71
72	integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
73	integer :: rain_cloud_above,kz_inv
74	! integer :: icloudtop !PS
75	real :: f_qvsat,pressure
76	! real :: rh,lsp,convp
77	real :: rh,lsp,convp,prec,rhmin,precmin
78	real :: rhoh(nuvzmax),pinmconv(nzmax)
79	real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
80	real :: xlon,ylat,xlonr,dzdx,dzdy
81	real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf
82	real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
83	real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
84	real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
85	real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
86	real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
87	real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
88	logical lconvectprec
89	real,parameter :: const=r_air/ga
90	parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics
91
92	logical :: init = .true.
93
94	!hg
95	integer :: cloud_ver,cloud_min, cloud_max
96	real :: cloud_col_wat, cloud_water
97	!hg temporary variables for testing
98	real :: rcw(0:nxmax-1,0:nymax-1)
99	real :: rpc(0:nxmax-1,0:nymax-1)
100	!hg
101
102	!*************************************************************************
103	! If verttransform is called the first time, initialize heights of the *
104	! z levels in meter. The heights are the heights of model levels, where *
105	! u,v,T and qv are given. *
106	!*************************************************************************
107
108	if (init) then
109
110	! Search for a point with high surface pressure (i.e. not above significant topography)
111	! Then, use this point to construct a reference z profile, to be used at all times
112	!**************************************************************************************
113
114	do jy=0,nymin1
115	do ix=0,nxmin1
116	if (ps(ix,jy,1,n).gt.100000.) then
117	ixm=ix
118	jym=jy
119	goto 3
120	endif
121	end do
122	end do
123	3 continue
124
125
126	tvold=tt2(ixm,jym,1,n)(1.+0.378ew(td2(ixm,jym,1,n))/ &
127	ps(ixm,jym,1,n))
128	pold=ps(ixm,jym,1,n)
129	height(1)=0.
130
131	do kz=2,nuvz
132	pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
133	tv=tth(ixm,jym,kz,n)(1.+0.608qvh(ixm,jym,kz,n))
134
135	if (abs(tv-tvold).gt.0.2) then
136	height(kz)=height(kz-1)+constlog(pold/pint) &
137	(tv-tvold)/log(tv/tvold)
138	else
139	height(kz)=height(kz-1)+constlog(pold/pint)tv
140	endif
141
142	tvold=tv
143	pold=pint
144	end do
145
146
147	! Determine highest levels that can be within PBL
148	!************************************************
149
150	do kz=1,nz
151	if (height(kz).gt.hmixmax) then
152	nmixz=kz
153	goto 9
154	endif
155	end do
156	9 continue
157
158	! Do not repeat initialization of the Cartesian z grid
159	!*****************************************************
160
161	init=.false.
162
163	endif
164
165
166	! Loop over the whole grid
167	!*************************
168
169	do jy=0,nymin1
170	do ix=0,nxmin1
171	tvold=tt2(ix,jy,1,n)(1.+0.378ew(td2(ix,jy,1,n))/ps(ix,jy,1,n))
172	pold=ps(ix,jy,1,n)
173	uvwzlev(ix,jy,1)=0.
174	wzlev(1)=0.
175	rhoh(1)=pold/(r_air*tvold)
176
177
178	! Compute heights of eta levels
179	!******************************
180
181	do kz=2,nuvz
182	pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)
183	tv=tth(ix,jy,kz,n)(1.+0.608qvh(ix,jy,kz,n))
184	rhoh(kz)=pint/(r_air*tv)
185
186	if (abs(tv-tvold).gt.0.2) then
187	uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+constlog(pold/pint) &
188	(tv-tvold)/log(tv/tvold)
189	else
190	uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+constlog(pold/pint)tv
191	endif
192
193	tvold=tv
194	pold=pint
195	end do
196
197
198	do kz=2,nwz-1
199	wzlev(kz)=(uvwzlev(ix,jy,kz+1)+uvwzlev(ix,jy,kz))/2.
200	end do
201	wzlev(nwz)=wzlev(nwz-1)+uvwzlev(ix,jy,nuvz)-uvwzlev(ix,jy,nuvz-1)
202
203	! pinmconv=(h2-h1)/(p2-p1)
204
205	pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
206	((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
207	(aknew(1)+bknew(1)*ps(ix,jy,1,n)))
208	do kz=2,nz-1
209	pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
210	((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
211	(aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
212	end do
213	pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
214	((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
215	(aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
216
217	! Levels, where u,v,t and q are given
218	!************************************
219
220	uu(ix,jy,1,n)=uuh(ix,jy,1)
221	vv(ix,jy,1,n)=vvh(ix,jy,1)
222	tt(ix,jy,1,n)=tth(ix,jy,1,n)
223	qv(ix,jy,1,n)=qvh(ix,jy,1,n)
224	!hg adding the cloud water
225	clwc(ix,jy,1,n)=clwch(ix,jy,1,n)
226	ciwc(ix,jy,1,n)=ciwch(ix,jy,1,n)
227	!hg
228	pv(ix,jy,1,n)=pvh(ix,jy,1)
229	rho(ix,jy,1,n)=rhoh(1)
230	uu(ix,jy,nz,n)=uuh(ix,jy,nuvz)
231	vv(ix,jy,nz,n)=vvh(ix,jy,nuvz)
232	tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
233	qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
234
235	!hg adding the cloud water
236	clwc(ix,jy,nz,n)=clwch(ix,jy,nuvz,n)
237	ciwc(ix,jy,nz,n)=ciwch(ix,jy,nuvz,n)
238	!hg
239	pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
240	rho(ix,jy,nz,n)=rhoh(nuvz)
241	kmin=2
242	do iz=2,nz-1
243	do kz=kmin,nuvz
244	if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
245	uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
246	vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
247	tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
248	qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
249	!hg adding the cloud water
250	clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
251	ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
252	!hg
253	pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
254	rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
255	goto 30
256	endif
257	if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
258	(height(iz).le.uvwzlev(ix,jy,kz))) then
259	dz1=height(iz)-uvwzlev(ix,jy,kz-1)
260	dz2=uvwzlev(ix,jy,kz)-height(iz)
261	dz=dz1+dz2
262	uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)dz2+uuh(ix,jy,kz)dz1)/dz
263	vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)dz2+vvh(ix,jy,kz)dz1)/dz
264	tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)dz2+tth(ix,jy,kz,n)dz1)/dz
265	qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)dz2+qvh(ix,jy,kz,n)dz1)/dz
266	!hg adding the cloud water
267	clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)dz2+clwch(ix,jy,kz,n)dz1)/dz
268	ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)dz2+ciwch(ix,jy,kz,n)dz1)/dz
269	!hg
270
271	pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)dz2+pvh(ix,jy,kz)dz1)/dz
272	rho(ix,jy,iz,n)=(rhoh(kz-1)dz2+rhoh(kz)dz1)/dz
273	kmin=kz
274	goto 30
275	endif
276	end do
277	30 continue
278	end do
279
280
281	! Levels, where w is given
282	!*************************
283
284	ww(ix,jy,1,n)=wwh(ix,jy,1)*pinmconv(1)
285	ww(ix,jy,nz,n)=wwh(ix,jy,nwz)*pinmconv(nz)
286	kmin=2
287	do iz=2,nz
288	do kz=kmin,nwz
289	if ((height(iz).gt.wzlev(kz-1)).and. &
290	(height(iz).le.wzlev(kz))) then
291	dz1=height(iz)-wzlev(kz-1)
292	dz2=wzlev(kz)-height(iz)
293	dz=dz1+dz2
294	ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)pinmconv(kz-1)dz2 &
295	+wwh(ix,jy,kz)pinmconv(kz)dz1)/dz
296	kmin=kz
297	goto 40
298	endif
299	end do
300	40 continue
301	end do
302
303	! Compute density gradients at intermediate levels
304	!*************************************************
305
306	drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
307	(height(2)-height(1))
308	do kz=2,nz-1
309	drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
310	(height(kz+1)-height(kz-1))
311	end do
312	drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
313
314	end do
315	end do
316
317
318	!****************************************************************
319	! Compute slope of eta levels in windward direction and resulting
320	! vertical wind correction
321	!****************************************************************
322
323	do jy=1,ny-2
324	cosf=cos((real(jy)dy+ylat0)pi180)
325	do ix=1,nx-2
326
327	kmin=2
328	do iz=2,nz-1
329
330	ui=uu(ix,jy,iz,n)*dxconst/cosf
331	vi=vv(ix,jy,iz,n)*dyconst
332
333	do kz=kmin,nz
334	if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
335	(height(iz).le.uvwzlev(ix,jy,kz))) then
336	dz1=height(iz)-uvwzlev(ix,jy,kz-1)
337	dz2=uvwzlev(ix,jy,kz)-height(iz)
338	dz=dz1+dz2
339	kl=kz-1
340	klp=kz
341	kmin=kz
342	goto 47
343	endif
344	end do
345
346	47 ix1=ix-1
347	jy1=jy-1
348	ixp=ix+1
349	jyp=jy+1
350
351	dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
352	dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
353	dzdx=(dzdx1dz2+dzdx2dz1)/dz
354
355	dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
356	dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
357	dzdy=(dzdy1dz2+dzdy2dz1)/dz
358
359	ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdxui+dzdyvi)
360
361	end do
362
363	end do
364	end do
365
366
367	! If north pole is in the domain, calculate wind velocities in polar
368	! stereographic coordinates
369	!*******************************************************************
370
371	if (nglobal) then
372	do jy=int(switchnorthg)-2,nymin1
373	ylat=ylat0+real(jy)*dy
374	do ix=0,nxmin1
375	xlon=xlon0+real(ix)*dx
376	do iz=1,nz
377	call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
378	vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
379	end do
380	end do
381	end do
382
383
384	do iz=1,nz
385
386	! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
387	!
388	! AMSnauffer Nov 18 2004 Added check for case vv=0
389	!
390	xlon=xlon0+real(nx/2-1)*dx
391	xlonr=xlon*pi/180.
392	ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)2+vv(nx/2-1,nymin1,iz,n)2)
393	if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
394	ddpol=atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr
395	else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
396	ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
397	vv(nx/2-1,nymin1,iz,n))-xlonr
398	else
399	ddpol=pi/2-xlonr
400	endif
401	if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
402	if(ddpol.gt.2.0pi) ddpol=ddpol-2.0pi
403
404	! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
405	xlon=180.0
406	xlonr=xlon*pi/180.
407	ylat=90.0
408	uuaux=-ffpol*sin(xlonr+ddpol)
409	vvaux=-ffpol*cos(xlonr+ddpol)
410	call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
411	jy=nymin1
412	do ix=0,nxmin1
413	uupol(ix,jy,iz,n)=uupolaux
414	vvpol(ix,jy,iz,n)=vvpolaux
415	end do
416	end do
417
418
419	! Fix: Set W at pole to the zonally averaged W of the next equator-
420	! ward parallel of latitude
421
422	do iz=1,nz
423	wdummy=0.
424	jy=ny-2
425	do ix=0,nxmin1
426	wdummy=wdummy+ww(ix,jy,iz,n)
427	end do
428	wdummy=wdummy/real(nx)
429	jy=nymin1
430	do ix=0,nxmin1
431	ww(ix,jy,iz,n)=wdummy
432	end do
433	end do
434
435	endif
436
437
438	! If south pole is in the domain, calculate wind velocities in polar
439	! stereographic coordinates
440	!*******************************************************************
441
442	if (sglobal) then
443	do jy=0,int(switchsouthg)+3
444	ylat=ylat0+real(jy)*dy
445	do ix=0,nxmin1
446	xlon=xlon0+real(ix)*dx
447	do iz=1,nz
448	call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
449	vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
450	end do
451	end do
452	end do
453
454	do iz=1,nz
455
456	! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
457	!
458	! AMSnauffer Nov 18 2004 Added check for case vv=0
459	!
460	xlon=xlon0+real(nx/2-1)*dx
461	xlonr=xlon*pi/180.
462	ffpol=sqrt(uu(nx/2-1,0,iz,n)2+vv(nx/2-1,0,iz,n)2)
463	if (vv(nx/2-1,0,iz,n).lt.0.) then
464	ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
465	else if (vv(nx/2-1,0,iz,n).gt.0.) then
466	ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
467	else
468	ddpol=pi/2-xlonr
469	endif
470	if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
471	if(ddpol.gt.2.0pi) ddpol=ddpol-2.0pi
472
473	! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
474	xlon=180.0
475	xlonr=xlon*pi/180.
476	ylat=-90.0
477	uuaux=+ffpol*sin(xlonr-ddpol)
478	vvaux=-ffpol*cos(xlonr-ddpol)
479	call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
480
481	jy=0
482	do ix=0,nxmin1
483	uupol(ix,jy,iz,n)=uupolaux
484	vvpol(ix,jy,iz,n)=vvpolaux
485	end do
486	end do
487
488
489	! Fix: Set W at pole to the zonally averaged W of the next equator-
490	! ward parallel of latitude
491
492	do iz=1,nz
493	wdummy=0.
494	jy=1
495	do ix=0,nxmin1
496	wdummy=wdummy+ww(ix,jy,iz,n)
497	end do
498	wdummy=wdummy/real(nx)
499	jy=0
500	do ix=0,nxmin1
501	ww(ix,jy,iz,n)=wdummy
502	end do
503	end do
504	endif
505
506
507	!write (,) 'initializing clouds, n:',n,nymin1,nxmin1,nz
508	! create a cloud and rainout/washout field, clouds occur where rh>80%
509	! total cloudheight is stored at level 0
510
511	if (readclouds) write(,) 'using cloud water from ECMWF'
512	if (.not.readclouds) write(,) 'using cloud water from parameterization'
513
514	rcw(:,:)=0
515	rpc(:,:)=0
516
517	do jy=0,nymin1
518	do ix=0,nxmin1
519	! OLD METHOD
520	rain_cloud_above=0
521	lsp=lsprec(ix,jy,1,n)
522	convp=convprec(ix,jy,1,n)
523	cloudsh(ix,jy,n)=0
524	do kz_inv=1,nz-1
525	kz=nz-kz_inv+1
526	pressure=rho(ix,jy,kz,n)r_airtt(ix,jy,kz,n)
527	rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
528	clouds(ix,jy,kz,n)=0
529	if (rh.gt.0.8) then ! in cloud
530	if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
531	rain_cloud_above=1
532	cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
533	height(kz)-height(kz-1)
534	if (lsp.ge.convp) then
535	clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
536	else
537	clouds(ix,jy,kz,n)=2 ! convp dominated rainout
538	endif
539	else ! no precipitation
540	clouds(ix,jy,kz,n)=1 ! cloud
541	endif
542	else ! no cloud
543	if (rain_cloud_above.eq.1) then ! scavenging
544	if (lsp.ge.convp) then
545	clouds(ix,jy,kz,n)=5 ! lsp dominated washout
546	else
547	clouds(ix,jy,kz,n)=4 ! convp dominated washout
548	endif
549	endif
550	endif
551	end do
552	!END OLD METHOD
553
554	! PS 2012
555	! lsp=lsprec(ix,jy,1,n)
556	! convp=convprec(ix,jy,1,n)
557	! prec=lsp+convp
558	! if (lsp.gt.convp) then ! prectype='lsp'
559	! lconvectprec = .false.
560	! else ! prectype='cp '
561	! lconvectprec = .true.
562	! endif
563	!HG METHOD
564	!readclouds =.true.
565	! if (readclouds) then
566	!hg added APR 2014 Cloud Water=clwc(ix,jy,kz,n) Cloud Ice=ciwc(ix,jy,kz,n)
567	!hg Use the cloud water variables to determine existence of clouds. This makes the PS code obsolete
568	! cloud_min=99999
569	! cloud_max=-1
570	! cloud_col_wat=0
571
572	! do kz=1, nz
573	! !clw & ciw are given in kg/kg
574	! cloud_water=clwc(ix,jy,kz,n)+ciwc(ix,jy,kz,n)
575	! if (cloud_water .gt. 0) then
576	! cloud_min=min(nint(height(kz)),cloud_min) !hg needs reset each grid
577	! cloud_max=max(nint(height(kz)),cloud_max) !hg needs reset each grid
578	! cloud_col_wat=cloud_col_wat+cloud_water !hg needs reset each grid
579	! endif
580	! cloud_ver=max(0,cloud_max-cloud_min)
581
582	! if (clwc(ix,jy,kz,n).gt.0 .or. ciwc(ix,jy,kz,n).gt.0) &
583	! !write(,) 'WATER',clwc(ix,jy,kz,n), 'ICE',ciwc(ix,jy,kz,n),'RH',rh,'KZ',kz &
584	! write(,) 'WATER',cloud_water,'RH',rh,'PREC',prec,'HEIGHT',height(kz) &
585	! enddo
586	! if ( cloud_min .ne. 99999 .and. cloud_max .ne. -1) write(,) 'CB', cloud_min, ' CT',cloud_max
587	! if (prec .gt. 0) write(,) 'PREC',prec,'Cloud Bot',cloud_min,'Cloud Top',cloud_max, 'Cloud Vert. ext',cloud_ver &
588	! ,'COLUMN cloud water',cloud_col_wat,'Conevctive' ,lconvectprec
589	! icloudbot(ix,jy,n)=cloud_min
590	! icloudthck(ix,jy,n)=cloud_ver
591	! rcw(ix,jy)=cloud_col_wat
592	! rpc(ix,jy)=prec
593	!write(,) 'Using clouds from ECMWF' !hg END Henrik Code
594	!END HG METHOD
595
596
597
598	! else ! windfields does not contain cloud data
599	! rhmin = 0.90 ! standard condition for presence of clouds
600	!PS note that original by Sabine Eckhart was 80%
601	!PS however, for T<-20 C we consider saturation over ice
602	!PS so I think 90% should be enough
603	! icloudbot(ix,jy,n)=icmv
604	! icloudtop=icmv ! this is just a local variable
605	!98 do kz=1,nz
606	! pressure=rho(ix,jy,kz,n)r_airtt(ix,jy,kz,n)
607	! rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
608	!ps if (prec.gt.0.01) print*,'relhum',prec,kz,rh,height(kz)
609	! if (rh .gt. rhmin) then
610	! if (icloudbot(ix,jy,n) .eq. icmv) then
611	! icloudbot(ix,jy,n)=nint(height(kz))
612	! endif
613	! icloudtop=nint(height(kz)) ! use int to save memory
614	! endif
615	! enddo
616	!PS try to get a cloud thicker than 50 m
617	!PS if there is at least .01 mm/h - changed to 0.002 and put into
618	!PS parameter precpmin
619	! if ((icloudbot(ix,jy,n) .eq. icmv .or. &
620	! icloudtop-icloudbot(ix,jy,n) .lt. 50) .and. &
621	! prec .gt. precmin) then
622	! rhmin = rhmin - 0.05
623	! if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum.
624	! endif
625
626	!PS is based on looking at a limited set of comparison data
627	! if (lconvectprec .and. icloudtop .lt. 6000 .and. &
628	! prec .gt. precmin) then
629	!
630	! if (convp .lt. 0.1) then
631	! icloudbot(ix,jy,n) = 500
632	! icloudtop = 8000
633	! else
634	! icloudbot(ix,jy,n) = 0
635	! icloudtop = 10000
636	! endif
637	! endif
638	! if (icloudtop .ne. icmv) then
639	! icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n)
640	! else
641	! icloudthck(ix,jy,n) = icmv
642	! endif
643	!PS get rid of too thin clouds
644	! if (icloudthck(ix,jy,n) .lt. 50) then
645	! icloudbot(ix,jy,n)=icmv
646	! icloudthck(ix,jy,n)=icmv
647	! endif
648	!hg__________________________________
649	! rcw(ix,jy)=2E-7prec*0.36
650	! rpc(ix,jy)=prec
651	!hg end______________________________
652
653	! endif !hg read clouds
654
655	end do
656	end do
657
658	end subroutine verttransform

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

Download in other formats:

Original Format