Context Navigation

verttransform_gfs.f90 @ db91eb7

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

Last change on this file since db91eb7 was db91eb7, checked in by Sabine <sabine.eckhardt@…>, 5 years ago
reading clouds in NCEP
Property mode set to `100644`
File size: 21.4 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 routine from verttransform to verttransform_gfs
53	!
54	!*****************************************************************************
55	! *
56	! Variables: *
57	! nx,ny,nz field dimensions in x,y and z direction *
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	! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition *
65	! *
66	!*****************************************************************************
67
68	use par_mod
69	use com_mod
70	use cmapf_mod
71
72	implicit none
73
74	integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
75	integer :: rain_cloud_above,kz_inv
76	real :: f_qvsat,pressure
77	real :: rh,lsp,cloudh_min,convp,prec
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	real,parameter :: const=r_air/ga
89
90	! NCEP version
91	integer :: llev, i
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. *
100	!*************************************************************************
101
102	if (init) then
103
104	! Search for a point with high surface pressure (i.e. not above significant topography)
105	! Then, use this point to construct a reference z profile, to be used at all times
106	!*****************************************************************************
107
108	do jy=0,nymin1
109	do ix=0,nxmin1
110	if (ps(ix,jy,1,n).gt.100000.) then
111	ixm=ix
112	jym=jy
113	goto 3
114	endif
115	end do
116	end do
117	3 continue
118
119
120	tvold=tt2(ixm,jym,1,n)(1.+0.378ew(td2(ixm,jym,1,n))/ &
121	ps(ixm,jym,1,n))
122	pold=ps(ixm,jym,1,n)
123	height(1)=0.
124
125	do kz=2,nuvz
126	pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
127	tv=tth(ixm,jym,kz,n)(1.+0.608qvh(ixm,jym,kz,n))
128
129	if (abs(tv-tvold).gt.0.2) then
130	height(kz)=height(kz-1)+constlog(pold/pint) &
131	(tv-tvold)/log(tv/tvold)
132	else
133	height(kz)=height(kz-1)+constlog(pold/pint)tv
134	endif
135
136	tvold=tv
137	pold=pint
138	end do
139
140
141	! Determine highest levels that can be within PBL
142	!************************************************
143
144	do kz=1,nz
145	if (height(kz).gt.hmixmax) then
146	nmixz=kz
147	goto 9
148	endif
149	end do
150	9 continue
151
152	! Do not repeat initialization of the Cartesian z grid
153	!*****************************************************
154
155	init=.false.
156
157	endif
158
159
160	! Loop over the whole grid
161	!*************************
162
163	do jy=0,nymin1
164	do ix=0,nxmin1
165
166	! NCEP version: find first level above ground
167	llev = 0
168	do i=1,nuvz
169	if (ps(ix,jy,1,n).lt.akz(i)) llev=i
170	end do
171	llev = llev+1
172	if (llev.gt.nuvz-2) llev = nuvz-2
173	! if (llev.eq.nuvz-2) write(,) 'verttransform
174	! +WARNING: LLEV eq NUZV-2'
175	! NCEP version
176
177
178	! compute height of pressure levels above ground
179	!***********************************************
180
181	tvold=tth(ix,jy,llev,n)(1.+0.608qvh(ix,jy,llev,n))
182	pold=akz(llev)
183	wzlev(llev)=0.
184	uvwzlev(ix,jy,llev)=0.
185	rhoh(llev)=pold/(r_air*tvold)
186
187	do kz=llev+1,nuvz
188	pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)
189	tv=tth(ix,jy,kz,n)(1.+0.608qvh(ix,jy,kz,n))
190	rhoh(kz)=pint/(r_air*tv)
191
192	if (abs(tv-tvold).gt.0.2) then
193	uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+constlog(pold/pint) &
194	(tv-tvold)/log(tv/tvold)
195	else
196	uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+constlog(pold/pint)tv
197	endif
198	wzlev(kz)=uvwzlev(ix,jy,kz)
199
200	tvold=tv
201	pold=pint
202	end do
203
204	! pinmconv=(h2-h1)/(p2-p1)
205
206	pinmconv(llev)=(uvwzlev(ix,jy,llev+1)-uvwzlev(ix,jy,llev))/ &
207	((aknew(llev+1)+bknew(llev+1)*ps(ix,jy,1,n))- &
208	(aknew(llev)+bknew(llev)*ps(ix,jy,1,n)))
209	do kz=llev+1,nz-1
210	pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
211	((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
212	(aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
213	end do
214	pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
215	((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
216	(aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
217
218
219	! Levels, where u,v,t and q are given
220	!************************************
221
222	uu(ix,jy,1,n)=uuh(ix,jy,llev)
223	vv(ix,jy,1,n)=vvh(ix,jy,llev)
224	tt(ix,jy,1,n)=tth(ix,jy,llev,n)
225	qv(ix,jy,1,n)=qvh(ix,jy,llev,n)
226	! IP & SEC, 201812 add clouds
227	if (readclouds) then
228	clwc(ix,jy,1,n)=clwch(ix,jy,llev,n)
229	endif
230	pv(ix,jy,1,n)=pvh(ix,jy,llev)
231	rho(ix,jy,1,n)=rhoh(llev)
232	pplev(ix,jy,1,n)=akz(llev)
233	uu(ix,jy,nz,n)=uuh(ix,jy,nuvz)
234	vv(ix,jy,nz,n)=vvh(ix,jy,nuvz)
235	tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
236	qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
237	! IP & SEC, 201812 add clouds
238	if (readclouds) then
239	clwc(ix,jy,nz,n)=clwch(ix,jy,nuvz,n)
240	endif
241	pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
242	rho(ix,jy,nz,n)=rhoh(nuvz)
243	pplev(ix,jy,nz,n)=akz(nuvz)
244	kmin=llev+1
245	do iz=2,nz-1
246	do kz=kmin,nuvz
247	if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
248	uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
249	vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
250	tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
251	qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
252	! IP & SEC, 201812 add clouds
253	if (readclouds) then
254	clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
255	endif
256	pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
257	rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
258	pplev(ix,jy,iz,n)=pplev(ix,jy,nz,n)
259	goto 30
260	endif
261	if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
262	(height(iz).le.uvwzlev(ix,jy,kz))) then
263	dz1=height(iz)-uvwzlev(ix,jy,kz-1)
264	dz2=uvwzlev(ix,jy,kz)-height(iz)
265	dz=dz1+dz2
266	uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)dz2+uuh(ix,jy,kz)dz1)/dz
267	vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)dz2+vvh(ix,jy,kz)dz1)/dz
268	tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
269	+tth(ix,jy,kz,n)*dz1)/dz
270	qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
271	+qvh(ix,jy,kz,n)*dz1)/dz
272	! IP & SEC, 201812 add clouds
273	if (readclouds) then
274	clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2 &
275	+clwch(ix,jy,kz,n)*dz1)/dz
276	endif
277	pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)dz2+pvh(ix,jy,kz)dz1)/dz
278	rho(ix,jy,iz,n)=(rhoh(kz-1)dz2+rhoh(kz)dz1)/dz
279	pplev(ix,jy,iz,n)=(akz(kz-1)dz2+akz(kz)dz1)/dz
280	endif
281	end do
282	30 continue
283	end do
284
285
286	! Levels, where w is given
287	!*************************
288
289	ww(ix,jy,1,n)=wwh(ix,jy,llev)*pinmconv(llev)
290	ww(ix,jy,nz,n)=wwh(ix,jy,nwz)*pinmconv(nz)
291	kmin=llev+1
292	do iz=2,nz
293	do kz=kmin,nwz
294	if ((height(iz).gt.wzlev(kz-1)).and. &
295	(height(iz).le.wzlev(kz))) then
296	dz1=height(iz)-wzlev(kz-1)
297	dz2=wzlev(kz)-height(iz)
298	dz=dz1+dz2
299	ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)pinmconv(kz-1)dz2 &
300	+wwh(ix,jy,kz)pinmconv(kz)dz1)/dz
301	endif
302	end do
303	end do
304
305
306	! Compute density gradients at intermediate levels
307	!*************************************************
308
309	drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
310	(height(2)-height(1))
311	do kz=2,nz-1
312	drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
313	(height(kz+1)-height(kz-1))
314	end do
315	drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
316
317	end do
318	end do
319
320
321	!****************************************************************
322	! Compute slope of eta levels in windward direction and resulting
323	! vertical wind correction
324	!****************************************************************
325
326	do jy=1,ny-2
327	cosf=cos((real(jy)dy+ylat0)pi180)
328	do ix=1,nx-2
329
330	! NCEP version: find first level above ground
331	llev = 0
332	do i=1,nuvz
333	if (ps(ix,jy,1,n).lt.akz(i)) llev=i
334	end do
335	llev = llev+1
336	if (llev.gt.nuvz-2) llev = nuvz-2
337	! if (llev.eq.nuvz-2) write(,) 'verttransform
338	! +WARNING: LLEV eq NUZV-2'
339	! NCEP version
340
341	kmin=llev+1
342	do iz=2,nz-1
343
344	ui=uu(ix,jy,iz,n)*dxconst/cosf
345	vi=vv(ix,jy,iz,n)*dyconst
346
347	do kz=kmin,nz
348	if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
349	(height(iz).le.uvwzlev(ix,jy,kz))) then
350	dz1=height(iz)-uvwzlev(ix,jy,kz-1)
351	dz2=uvwzlev(ix,jy,kz)-height(iz)
352	dz=dz1+dz2
353	kl=kz-1
354	klp=kz
355	goto 47
356	endif
357	end do
358
359	47 ix1=ix-1
360	jy1=jy-1
361	ixp=ix+1
362	jyp=jy+1
363
364	dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
365	dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
366	dzdx=(dzdx1dz2+dzdx2dz1)/dz
367
368	dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
369	dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
370	dzdy=(dzdy1dz2+dzdy2dz1)/dz
371
372	ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdxui+dzdyvi)
373
374	end do
375
376	end do
377	end do
378
379
380	! If north pole is in the domain, calculate wind velocities in polar
381	! stereographic coordinates
382	!*******************************************************************
383
384	if (nglobal) then
385	do jy=int(switchnorthg)-2,nymin1
386	ylat=ylat0+real(jy)*dy
387	do ix=0,nxmin1
388	xlon=xlon0+real(ix)*dx
389	do iz=1,nz
390	call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
391	vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
392	vvpol(ix,jy,iz,n))
393	end do
394	end do
395	end do
396
397
398	do iz=1,nz
399
400	! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
401	xlon=xlon0+real(nx/2-1)*dx
402	xlonr=xlon*pi/180.
403	ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)2+vv(nx/2-1,nymin1,iz,n)2)
404	if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
405	ddpol=atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr
406	elseif (vv(nx/2-1,nymin1,iz,n).gt.0.) then
407	ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
408	vv(nx/2-1,nymin1,iz,n))-xlonr
409	else
410	ddpol=pi/2-xlonr
411	endif
412	if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
413	if(ddpol.gt.2.0pi) ddpol=ddpol-2.0pi
414
415	! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
416	xlon=180.0
417	xlonr=xlon*pi/180.
418	ylat=90.0
419	uuaux=-ffpol*sin(xlonr+ddpol)
420	vvaux=-ffpol*cos(xlonr+ddpol)
421	call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
422	jy=nymin1
423	do ix=0,nxmin1
424	uupol(ix,jy,iz,n)=uupolaux
425	vvpol(ix,jy,iz,n)=vvpolaux
426	end do
427	end do
428
429
430	! Fix: Set W at pole to the zonally averaged W of the next equator-
431	! ward parallel of latitude
432
433	do iz=1,nz
434	wdummy=0.
435	jy=ny-2
436	do ix=0,nxmin1
437	wdummy=wdummy+ww(ix,jy,iz,n)
438	end do
439	wdummy=wdummy/real(nx)
440	jy=nymin1
441	do ix=0,nxmin1
442	ww(ix,jy,iz,n)=wdummy
443	end do
444	end do
445
446	endif
447
448
449	! If south pole is in the domain, calculate wind velocities in polar
450	! stereographic coordinates
451	!*******************************************************************
452
453	if (sglobal) then
454	do jy=0,int(switchsouthg)+3
455	ylat=ylat0+real(jy)*dy
456	do ix=0,nxmin1
457	xlon=xlon0+real(ix)*dx
458	do iz=1,nz
459	call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
460	vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
461	end do
462	end do
463	end do
464
465	do iz=1,nz
466
467	! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
468	xlon=xlon0+real(nx/2-1)*dx
469	xlonr=xlon*pi/180.
470	ffpol=sqrt(uu(nx/2-1,0,iz,n)2+vv(nx/2-1,0,iz,n)2)
471	if(vv(nx/2-1,0,iz,n).lt.0.) then
472	ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
473	elseif (vv(nx/2-1,0,iz,n).gt.0.) then
474	ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))-xlonr
475	else
476	ddpol=pi/2-xlonr
477	endif
478	if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
479	if(ddpol.gt.2.0pi) ddpol=ddpol-2.0pi
480
481	! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
482	xlon=180.0
483	xlonr=xlon*pi/180.
484	ylat=-90.0
485	uuaux=+ffpol*sin(xlonr-ddpol)
486	vvaux=-ffpol*cos(xlonr-ddpol)
487	call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
488
489	jy=0
490	do ix=0,nxmin1
491	uupol(ix,jy,iz,n)=uupolaux
492	vvpol(ix,jy,iz,n)=vvpolaux
493	end do
494	end do
495
496
497	! Fix: Set W at pole to the zonally averaged W of the next equator-
498	! ward parallel of latitude
499
500	do iz=1,nz
501	wdummy=0.
502	jy=1
503	do ix=0,nxmin1
504	wdummy=wdummy+ww(ix,jy,iz,n)
505	end do
506	wdummy=wdummy/real(nx)
507	jy=0
508	do ix=0,nxmin1
509	ww(ix,jy,iz,n)=wdummy
510	end do
511	end do
512	endif
513
514
515
516	!***********************************************************************************
517	! IP & SEC, 201812 GFS clouds read
518	if (readclouds) then
519	! The method is loops all grids vertically and constructs the 3D matrix for clouds
520	! Cloud top and cloud bottom gid cells are assigned as well as the total column
521	! cloud water. For precipitating grids, the type and whether it is in or below
522	! cloud scavenging are assigned with numbers 2-5 (following the old metod).
523	! Distinction is done for lsp and convp though they are treated the same in regards
524	! to scavenging. Also clouds that are not precipitating are defined which may be
525	! to include future cloud processing by non-precipitating-clouds.
526	!***********************************************************************************
527	write(,) 'Global NCEP fields: using cloud water'
528	clw(:,:,:,n)=0.0
529	ctwc(:,:,n)=0.0
530	clouds(:,:,:,n)=0
531	! If water/ice are read separately into clwc and ciwc, store sum in clwc
532	do jy=0,nymin1
533	do ix=0,nxmin1
534	lsp=lsprec(ix,jy,1,n)
535	convp=convprec(ix,jy,1,n)
536	prec=lsp+convp
537	! Find clouds in the vertical
538	do kz=1, nz-1 !go from top to bottom
539	if (clwc(ix,jy,kz,n).gt.0) then
540	! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3
541	clw(ix,jy,kz,n)=(clwc(ix,jy,kz,n)rho(ix,jy,kz,n))(height(kz+1)-height(kz))
542	ctwc(ix,jy,n) = ctwc(ix,jy,n)+clw(ix,jy,kz,n)
543	cloudh_min=min(height(kz+1),height(kz))
544	endif
545	end do
546
547	! If Precipitation. Define removal type in the vertical
548	if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
549
550	do kz=nz,1,-1 !go Bottom up!
551	if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud
552	cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1)
553	clouds(ix,jy,kz,n)=1 ! is a cloud
554	if (lsp.ge.convp) then
555	clouds(ix,jy,kz,n)=3 ! lsp in-cloud
556	else
557	clouds(ix,jy,kz,n)=2 ! convp in-cloud
558	endif ! convective or large scale
559	elseif((clw(ix,jy,kz,n).le.0) .and. (cloudh_min.ge.height(kz))) then ! is below cloud
560	if (lsp.ge.convp) then
561	clouds(ix,jy,kz,n)=5 ! lsp dominated washout
562	else
563	clouds(ix,jy,kz,n)=4 ! convp dominated washout
564	endif ! convective or large scale
565	endif
566
567	if (height(kz).ge. 19000) then ! set a max height for removal
568	clouds(ix,jy,kz,n)=0
569	endif !clw>0
570	end do !nz
571	endif ! precipitation
572	end do
573	end do
574	else
575	write(,) 'Global NCEP fields: using cloud water from Parameterization'
576	! write (,) 'initializing clouds, n:',n,nymin1,nxmin1,nz
577	! create a cloud and rainout/washout field, clouds occur where rh>80%
578	! total cloudheight is stored at level 0
579	do jy=0,nymin1
580	do ix=0,nxmin1
581	rain_cloud_above=0
582	lsp=lsprec(ix,jy,1,n)
583	convp=convprec(ix,jy,1,n)
584	cloudsh(ix,jy,n)=0
585	do kz_inv=1,nz-1
586	kz=nz-kz_inv+1
587	pressure=rho(ix,jy,kz,n)r_airtt(ix,jy,kz,n)
588	rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
589	clouds(ix,jy,kz,n)=0
590	if (rh.gt.0.8) then ! in cloud
591	if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
592	rain_cloud_above=1
593	cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1)
594	if (lsp.ge.convp) then
595	clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
596	else
597	clouds(ix,jy,kz,n)=2 ! convp dominated rainout
598	endif
599	else ! no precipitation
600	clouds(ix,jy,kz,n)=1 ! cloud
601	endif
602	else ! no cloud
603	if (rain_cloud_above.eq.1) then ! scavenging
604	if (lsp.ge.convp) then
605	clouds(ix,jy,kz,n)=5 ! lsp dominated washout
606	else
607	clouds(ix,jy,kz,n)=4 ! convp dominated washout
608	endif
609	endif
610	endif
611	end do
612	end do
613	end do
614	endif ! IP & SEC 201812, GFS clouds read
615
616
617	end subroutine verttransform_gfs

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

Download in other formats:

Original Format