Context Navigation

verttransform_gfs.f90 @ f251e57

univie

Last change on this file since f251e57 was f251e57, checked in by pesei <petra seibert at univie ac at>, 6 years ago
fix ticket:140 (ref position for vert. grid) in verttransform_gfs, minor changes in both verttransforms
Property mode set to `100644`
File size: 20.3 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_gfs(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	! CHANGE 17/11/2005 Caroline Forster, NCEP GFS version *
41	! *
42	! - Vertical levels for u, v and w are put together *
43	! - Slope correction for vertical velocity: Modification of calculation *
44	! procedure *
45	! *
46	!*****************************************************************************
47	! Changes, Bernd C. Krueger, Feb. 2001:
48	! Variables tth and qvh (on eta coordinates) from common block
49	!
50	! Unified ECMWF and GFS builds
51	! Marian Harustak, 12.5.2017
52	! - Renamed from verttransform to verttransform_ecmwf
53	!
54	! undocumented modifications by NILU for v10 *
55	! *
56	! Petra Seibert, 2018-06-13: *
57	! - fix bug of ticket:140 (find reference position for vertical grid) *
58	! - put back SAVE attribute for INIT, just to be safe *
59	! - minor changes, most of them just cosmetics *
60	! for details see changelog.txt *
61	! *
62	!*****************************************************************************
63	! *
64	! Variables: *
65	! nx,ny,nz field dimensions in x,y and z direction *
66	! uu(0:nxmax,0:nymax,nzmax,2) wind components in x-direction [m/s] *
67	! vv(0:nxmax,0:nymax,nzmax,2) wind components in y-direction [m/s] *
68	! ww(0:nxmax,0:nymax,nzmax,2) wind components in z-direction [deltaeta/s]*
69	! tt(0:nxmax,0:nymax,nzmax,2) temperature [K] *
70	! pv(0:nxmax,0:nymax,nzmax,2) potential voriticity (pvu) *
71	! ps(0:nxmax,0:nymax,2) surface pressure [Pa] *
72	! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition *
73	! *
74	!*****************************************************************************
75
76	use par_mod
77	use com_mod
78	use cmapf_mod
79
80	implicit none
81
82	integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp
83	integer :: rain_cloud_above,kz_inv
84	real :: f_qvsat,pressure
85	real :: rh,lsp,convp
86	real :: rhoh(nuvzmax),pinmconv(nzmax)
87	real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
88
89	!> for reference profile (PS)
90	real :: tvoldref, poldref, pintref, psmean, psstd
91	integer :: ixyref(2), ixref,jyref
92	integer, parameter :: ioldref = 1 ! PS 2018-06-13: set to 0 if you
93	!! want old method of searching reference location for the vertical grid
94	!! 1 for new method (options for other methods 2,... in the future)
95
96	real :: xlon,ylat,xlonr,dzdx,dzdy
97	real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf
98	real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
99	real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
100	real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
101	real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
102	real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
103	real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
104	real,parameter :: const=r_air/ga
105
106	logical, save :: init = .true. ! PS 2018-06-13: add back save attribute
107
108	!> GFS version
109	integer :: llev, i
110
111
112
113	!*************************************************************************
114	! If verttransform is called the first time, initialize heights of the *
115	! z levels in meter. The heights are the heights of model levels, where *
116	! u,v,T and qv are given. *
117	!*************************************************************************
118
119	if (init) then
120
121	! Search for a point with high surface pressure (i.e. not above significant topography)
122	! Then, use this point to construct a reference z profile, to be used at all times
123	!*****************************************************************************
124
125	if (ioldref .eq. 0) then ! old reference grid point
126	do jy=0,nymin1
127	do ix=0,nxmin1
128	if (ps(ix,jy,1,n).gt.1000.e2) then
129	ixref=ix
130	jyref=jy
131	goto 3
132	endif
133	end do
134	end do
135	3 continue
136	! print*,'oldheights at' ,ixref,jyref,ps(ixref,jyref,1,n)
137	else ! new reference grid point
138	! PS: the old version fails if the pressure is <=1000 hPa in the whole
139	! domain. Let us find a good replacement, not just a quick fix.
140	! Search point near to mean pressure after excluding mountains
141
142	psmean = sum( ps(:,:,1,n) ) / (nx*ny)
143	! print*,'height: fg psmean',psmean
144	psstd = sqrt(sum( (ps(:,:,1,n) - psmean)*2 ) / (nxny))
145
146	!> new psmean using only points within $\plusminus\sigma$
147	! psmean = sum( ps(:,:,1,n), abs(ps(:,:,1,n)-psmean) < psstd ) / &
148	! count(abs(ps(:,:,1,n)-psmean) < psstd)
149
150	!> new psmean using only points with $p\gt \overbar{p}+\sigma_p$
151	!> (reject mountains, accept valleys)
152	psmean = sum( ps(:,:,1,n), ps(:,:,1,n) > psmean - psstd ) / &
153	count(ps(:,:,1,n) > psmean - psstd)
154	! print*,'height: std, new psmean',psstd,psmean
155	ixyref = minloc( abs( ps(:,:,1,n) - psmean ) )
156	ixref = ixyref(1)
157	jyref = ixyref(2)
158	! print*,'newheights at' ,ixref,jyref,ps(ixref,jyref,1,n)
159	endif
160
161	tvoldref=tt2(ixref,jyref,1,n)* &
162	( 1. + 0.378*ew(td2(ixref,jyref,1,n) ) / ps(ixref,jyref,1,n))
163	poldref=ps(ixref,jyref,1,n)
164	height(1)=0.
165	kz=1
166	! print*,'height=',kz,height(kz),tvoldref,poldref
167
168	do kz=2,nuvz
169
170	pintref = akz(kz)
171	! Note that for GFS data, the akz variable contains the input level
172	! pressure values. bkz is zero (I removed all terms with bkz that
173	! were erroneously copied from verttransform_ecmwf). [PS, June 2018]
174
175	tv = tth(ixref,jyref,kz,n)(1.+0.608qvh(ixref,jyref,kz,n))
176
177	if (abs(tv-tvoldref) .gt. 0.2) then
178	height(kz)=height(kz-1) + &
179	const * log(poldref/pintref) * (tv-tvoldref) / log(tv/tvoldref)
180	else
181	height(kz) = height(kz-1) + constlog(poldref/pintref)tv
182	endif
183
184	tvoldref=tv
185	poldref=pintref
186	! print*,'height=',kz,height(kz),tvoldref,poldref
187	end do
188
189
190	! Determine highest levels that can be within PBL
191	!************************************************
192
193	do kz=1,nz
194	if (height(kz).gt.hmixmax) then
195	nmixz=kz
196	goto 9
197	endif
198	end do
199	9 continue
200
201	! Do not repeat initialization of the Cartesian z grid
202	!*****************************************************
203
204	init=.false.
205
206	endif ! init block (vertical grid construction)
207
208
209	! Loop over the whole grid
210	!*************************
211
212	do jy=0,nymin1
213	do ix=0,nxmin1
214
215	! NCEP version: find first level above ground
216	llev = 0
217	do i=1,nuvz
218	if (ps(ix,jy,1,n).lt.akz(i)) llev=i
219	end do
220	llev = llev+1
221	if (llev.gt.nuvz-2) llev = nuvz-2
222	! if (llev.eq.nuvz-2) write(,) 'verttransform
223	! +WARNING: LLEV eq NUZV-2'
224	! NCEP version
225
226
227	! compute height of pressure levels above ground
228	!***********************************************
229
230	tvold=tth(ix,jy,llev,n)(1.+0.608qvh(ix,jy,llev,n))
231	pold=akz(llev)
232	wzlev(llev)=0.
233	uvwzlev(ix,jy,llev)=0.
234	rhoh(llev)=pold/(r_air*tvold)
235
236	do kz=llev+1,nuvz
237	pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)
238	tv=tth(ix,jy,kz,n)(1.+0.608qvh(ix,jy,kz,n))
239	rhoh(kz)=pint/(r_air*tv)
240
241	if (abs(tv-tvold).gt.0.2) then
242	uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+constlog(pold/pint) &
243	(tv-tvold)/log(tv/tvold)
244	else
245	uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+constlog(pold/pint)tv
246	endif
247	wzlev(kz)=uvwzlev(ix,jy,kz)
248
249	tvold=tv
250	pold=pint
251	end do
252
253	! pinmconv=(h2-h1)/(p2-p1)
254
255	pinmconv(llev)=(uvwzlev(ix,jy,llev+1)-uvwzlev(ix,jy,llev))/ &
256	((aknew(llev+1)+bknew(llev+1)*ps(ix,jy,1,n))- &
257	(aknew(llev)+bknew(llev)*ps(ix,jy,1,n)))
258	do kz=llev+1,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
268	! Levels, where u,v,t and q are given
269	!************************************
270
271	uu(ix,jy,1,n)=uuh(ix,jy,llev)
272	vv(ix,jy,1,n)=vvh(ix,jy,llev)
273	tt(ix,jy,1,n)=tth(ix,jy,llev,n)
274	qv(ix,jy,1,n)=qvh(ix,jy,llev,n)
275	pv(ix,jy,1,n)=pvh(ix,jy,llev)
276	rho(ix,jy,1,n)=rhoh(llev)
277	pplev(ix,jy,1,n)=akz(llev)
278	uu(ix,jy,nz,n)=uuh(ix,jy,nuvz)
279	vv(ix,jy,nz,n)=vvh(ix,jy,nuvz)
280	tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
281	qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
282	pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
283	rho(ix,jy,nz,n)=rhoh(nuvz)
284	pplev(ix,jy,nz,n)=akz(nuvz)
285	kmin=llev+1
286	do iz=2,nz-1
287	do kz=kmin,nuvz
288	if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
289	uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
290	vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
291	tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
292	qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
293	pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
294	rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
295	pplev(ix,jy,iz,n)=pplev(ix,jy,nz,n)
296	goto 30
297	endif
298	if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
299	(height(iz).le.uvwzlev(ix,jy,kz))) then
300	dz1=height(iz)-uvwzlev(ix,jy,kz-1)
301	dz2=uvwzlev(ix,jy,kz)-height(iz)
302	dz=dz1+dz2
303	uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)dz2+uuh(ix,jy,kz)dz1)/dz
304	vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)dz2+vvh(ix,jy,kz)dz1)/dz
305	tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
306	+tth(ix,jy,kz,n)*dz1)/dz
307	qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
308	+qvh(ix,jy,kz,n)*dz1)/dz
309	pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)dz2+pvh(ix,jy,kz)dz1)/dz
310	rho(ix,jy,iz,n)=(rhoh(kz-1)dz2+rhoh(kz)dz1)/dz
311	pplev(ix,jy,iz,n)=(akz(kz-1)dz2+akz(kz)dz1)/dz
312	endif
313	end do
314	30 continue
315	end do
316
317
318	! Levels, where w is given
319	!*************************
320
321	ww(ix,jy,1,n)=wwh(ix,jy,llev)*pinmconv(llev)
322	ww(ix,jy,nz,n)=wwh(ix,jy,nwz)*pinmconv(nz)
323	kmin=llev+1
324	do iz=2,nz
325	do kz=kmin,nwz
326	if ((height(iz).gt.wzlev(kz-1)).and. &
327	(height(iz).le.wzlev(kz))) then
328	dz1=height(iz)-wzlev(kz-1)
329	dz2=wzlev(kz)-height(iz)
330	dz=dz1+dz2
331	ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)pinmconv(kz-1)dz2 &
332	+wwh(ix,jy,kz)pinmconv(kz)dz1)/dz
333	endif
334	end do
335	end do
336
337
338	! Compute density gradients at intermediate levels
339	!*************************************************
340
341	drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
342	(height(2)-height(1))
343	do kz=2,nz-1
344	drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
345	(height(kz+1)-height(kz-1))
346	end do
347	drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
348
349	end do
350	end do
351
352
353	!****************************************************************
354	! Compute slope of eta levels in windward direction and resulting
355	! vertical wind correction
356	!****************************************************************
357
358	do jy=1,ny-2
359	cosf=cos((real(jy)dy+ylat0)pi180)
360	do ix=1,nx-2
361
362	! NCEP version: find first level above ground
363	llev = 0
364	do i=1,nuvz
365	if (ps(ix,jy,1,n).lt.akz(i)) llev=i
366	end do
367	llev = llev+1
368	if (llev.gt.nuvz-2) llev = nuvz-2
369	! if (llev.eq.nuvz-2) write(,) 'verttransform
370	! +WARNING: LLEV eq NUZV-2'
371	! NCEP version
372
373	kmin=llev+1
374	do iz=2,nz-1
375
376	ui=uu(ix,jy,iz,n)*dxconst/cosf
377	vi=vv(ix,jy,iz,n)*dyconst
378
379	do kz=kmin,nz
380	if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
381	(height(iz).le.uvwzlev(ix,jy,kz))) then
382	dz1=height(iz)-uvwzlev(ix,jy,kz-1)
383	dz2=uvwzlev(ix,jy,kz)-height(iz)
384	dz=dz1+dz2
385	kl=kz-1
386	klp=kz
387	goto 47
388	endif
389	end do
390
391	47 ix1=ix-1
392	jy1=jy-1
393	ixp=ix+1
394	jyp=jy+1
395
396	dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
397	dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
398	dzdx=(dzdx1dz2+dzdx2dz1)/dz
399
400	dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
401	dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
402	dzdy=(dzdy1dz2+dzdy2dz1)/dz
403
404	ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdxui+dzdyvi)
405
406	end do
407
408	end do
409	end do
410
411
412	! If north pole is in the domain, calculate wind velocities in polar
413	! stereographic coordinates
414	!*******************************************************************
415
416	if (nglobal) then
417	do jy=int(switchnorthg)-2,nymin1
418	ylat=ylat0+real(jy)*dy
419	do ix=0,nxmin1
420	xlon=xlon0+real(ix)*dx
421	do iz=1,nz
422	call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
423	vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
424	vvpol(ix,jy,iz,n))
425	end do
426	end do
427	end do
428
429
430	do iz=1,nz
431
432	! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
433	xlon=xlon0+real(nx/2-1)*dx
434	xlonr=xlon*pi/180.
435	ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)2+vv(nx/2-1,nymin1,iz,n)2)
436	if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
437	ddpol=atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr
438	elseif (vv(nx/2-1,nymin1,iz,n).gt.0.) then
439	ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
440	vv(nx/2-1,nymin1,iz,n))-xlonr
441	else
442	ddpol=pi/2-xlonr
443	endif
444	if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
445	if(ddpol.gt.2.0pi) ddpol=ddpol-2.0pi
446
447	! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
448	xlon=180.0
449	xlonr=xlon*pi/180.
450	ylat=90.0
451	uuaux=-ffpol*sin(xlonr+ddpol)
452	vvaux=-ffpol*cos(xlonr+ddpol)
453	call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
454	jy=nymin1
455	do ix=0,nxmin1
456	uupol(ix,jy,iz,n)=uupolaux
457	vvpol(ix,jy,iz,n)=vvpolaux
458	end do
459	end do
460
461
462	! Fix: Set W at pole to the zonally averaged W of the next equator-
463	! ward parallel of latitude
464
465	do iz=1,nz
466	wdummy=0.
467	jy=ny-2
468	do ix=0,nxmin1
469	wdummy=wdummy+ww(ix,jy,iz,n)
470	end do
471	wdummy=wdummy/real(nx)
472	jy=nymin1
473	do ix=0,nxmin1
474	ww(ix,jy,iz,n)=wdummy
475	end do
476	end do
477
478	endif
479
480
481	! If south pole is in the domain, calculate wind velocities in polar
482	! stereographic coordinates
483	!*******************************************************************
484
485	if (sglobal) then
486	do jy=0,int(switchsouthg)+3
487	ylat=ylat0+real(jy)*dy
488	do ix=0,nxmin1
489	xlon=xlon0+real(ix)*dx
490	do iz=1,nz
491	call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
492	vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
493	end do
494	end do
495	end do
496
497	do iz=1,nz
498
499	! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
500	xlon=xlon0+real(nx/2-1)*dx
501	xlonr=xlon*pi/180.
502	ffpol=sqrt(uu(nx/2-1,0,iz,n)2+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)/vv(nx/2-1,0,iz,n))+xlonr
505	elseif (vv(nx/2-1,0,iz,n).gt.0.) then
506	ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))-xlonr
507	else
508	ddpol=pi/2-xlonr
509	endif
510	if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
511	if(ddpol.gt.2.0pi) ddpol=ddpol-2.0pi
512
513	! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
514	xlon=180.0
515	xlonr=xlon*pi/180.
516	ylat=-90.0
517	uuaux=+ffpol*sin(xlonr-ddpol)
518	vvaux=-ffpol*cos(xlonr-ddpol)
519	call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
520
521	jy=0
522	do ix=0,nxmin1
523	uupol(ix,jy,iz,n)=uupolaux
524	vvpol(ix,jy,iz,n)=vvpolaux
525	end do
526	end do
527
528
529	! Fix: Set W at pole to the zonally averaged W of the next equator-
530	! ward parallel of latitude
531
532	do iz=1,nz
533	wdummy=0.
534	jy=1
535	do ix=0,nxmin1
536	wdummy=wdummy+ww(ix,jy,iz,n)
537	end do
538	wdummy=wdummy/real(nx)
539	jy=0
540	do ix=0,nxmin1
541	ww(ix,jy,iz,n)=wdummy
542	end do
543	end do
544	endif
545
546
547	! write (,) 'initializing clouds, n:',n,nymin1,nxmin1,nz
548	! create a cloud and rainout/washout field, clouds occur where rh>80%
549	! total cloudheight is stored at level 0
550	do jy=0,nymin1
551	do ix=0,nxmin1
552	rain_cloud_above=0
553	lsp=lsprec(ix,jy,1,n)
554	convp=convprec(ix,jy,1,n)
555	cloudsh(ix,jy,n)=0
556	do kz_inv=1,nz-1
557	kz=nz-kz_inv+1
558	pressure=rho(ix,jy,kz,n)r_airtt(ix,jy,kz,n)
559	rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
560	clouds(ix,jy,kz,n)=0
561	if (rh.gt.0.8) then ! in cloud
562	if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
563	rain_cloud_above=1
564	cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1)
565	if (lsp.ge.convp) then
566	clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
567	else
568	clouds(ix,jy,kz,n)=2 ! convp dominated rainout
569	endif
570	else ! no precipitation
571	clouds(ix,jy,kz,n)=1 ! cloud
572	endif
573	else ! no cloud
574	if (rain_cloud_above.eq.1) then ! scavenging
575	if (lsp.ge.convp) then
576	clouds(ix,jy,kz,n)=5 ! lsp dominated washout
577	else
578	clouds(ix,jy,kz,n)=4 ! convp dominated washout
579	endif
580	endif
581	endif
582	end do
583	end do
584	end do
585
586
587	end subroutine verttransform_gfs

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

Download in other formats:

Original Format