Context Navigation

verttransform_gfs.f90 @ b398fb6

FPv9.3.1FPv9.3.1b_testingFPv9.3.2fp9.3.1-20161214-nc4grib2nc4_repair

Last change on this file since b398fb6 was 496c607, checked in by Don Morton <Don.Morton@…>, 8 years ago

Initial commit of FPv9.3.1

Currently, this is a clone of snapshot FPv9.3.0

Property mode set to 100644

File size: 24.5 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	! Changes Arnold, D. and Morton, D. (2015): *
51	! -- description of local and common variables *
52	!*****************************************************************************
53
54	use par_mod
55	use com_mod
56	use cmapf_mod
57
58	implicit none
59	!***********************************************************************
60	! Subroutine Parameters: *
61	! input *
62	! n temporal index for meteorological fields (1 to 3)*
63	! uuh,vvh, wwh wind components in ECMWF model levels *
64	! pvh potential vorticity in ECMWF model levels *
65	integer :: n
66	real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
67	real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
68	real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
69	real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
70	!***********************************************************************
71
72	!***********************************************************************
73	! Local variables *
74	! *
75	! ew subroutine/function to calculate saturation *
76	! water vaport for a given air temperature *
77	! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition *
78	! rain_cloud_above [0,1] whether there is a raining cloud or not *
79	! kz_inv inverted indez for kz *
80	! f_qvsat Saturation water vapor specific humidity (kg/kg) *
81	! pressure *
82	! rh relative humidity *
83	! lsp large scale precip at one grid cell *
84	! convp convective precip at one grid cell *
85	! uvzlev height of the eta half-levels *
86	! rhoh density in model levels *
87	! pinmconv conversion factor *
88	! pint pressure on model levels (using akz, bkz, ps) *
89	! tv virtual temperature *
90	! tvold,pold dummy variables to keep previous value *
91	! dz1, dz2, dz differences heights model levels, pressure levels *
92	! ui, vi *
93	! xlon,xlat *
94	! xlonr xlonpi/180.
95	! dzdx,dzdy slope corrections *
96	! dzdx1,dzdx2,dzdy1 slope corrections in each direction *
97	! uuaux,vvaux,uupolaux,vvpolaux auxiliary variables for polar sterero.*
98	! ddpol,ffpol for special treatment of poles *
99	! wdummy *
100	! wzlev *
101	! uvwzlev *
102
103	integer :: rain_cloud_above,kz_inv
104	real :: f_qvsat,pressure
105	real :: rh,lsp,convp
106	real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
107	real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
108	real :: xlon,ylat,xlonr,dzdx,dzdy
109	real :: dzdx1,dzdx2,dzdy1,dzdy2
110	real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
111	real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
112
113	! NCEP version
114	integer :: llev, i
115
116	! Other variables:
117	! ix,jy,kz,kzmin loop control indices in each direction *
118	! kl,klp
119	! ix1,jy1,ixp,jyp,ixm,jym
120	integer :: ix,jy,kz,iz,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
121	!***********************************************************************
122
123	!***********************************************************************
124	! Local constants *
125	real,parameter :: const=r_air/ga
126	!***********************************************************************
127
128	!***********************************************************************
129	! Global variables *
130	! from par_mod and com_mod *
131	! nxmin1,nymin1 nx-1, ny-1, respectively *
132	! tt2 2 meter temperature *
133	! ps surface pressure *
134	! td2 2 meter dew point *
135	! height heights of all levels *
136	! akm, bkm coeffs. which regulate vertical discretization of ecmwf *
137	! akz, bkz model discretization coeffizients at the centre of layers *
138	! tv virtual temperature *
139	! nuvz,nwz vertical dimension of original ECMWF data *
140	! nx,ny,nz actual dimensions of wind fields in x,y and z*
141	! aknew,bknew model discretization coeffs. at the interpolated levels *
142	! uu,vv,ww [m/2] wind components in x,y and z direction *
143	! qv specific humidity data *
144	! pv (pvu) potential vorticity *
145	! rho [kg/m3] air density *
146	! drhodz [kg/m2] vertical air density gradient *
147	! dxconst,dyconst auxiliary variables for utransform,vtransform*
148	! ylat0 geographical latitude of lower left grid point *
149	! xlon0 geographical longitude of lower left grid point*
150	! uupol,vvpol [m/s] wind components in polar stereographic projection*
151	! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition *
152	! lsprec [mm/h] large scale total precipitation *
153	! convprec [mm/h] convective precipitation *
154	! cloudsh *
155	! *
156	!***********************************************************************
157
158	!-----------------------------------------------------------------------------
159
160
161	logical :: init = .true.
162
163
164	!*************************************************************************
165	! If verttransform is called the first time, initialize heights of the *
166	! z levels in meter. The heights are the heights of model levels, where *
167	! u,v,T and qv are given, and of the interfaces, where w is given. So, *
168	! the vertical resolution in the z system is doubled. As reference point,*
169	! the lower left corner of the grid is used. *
170	! Unlike in the eta system, no difference between heights for u,v and *
171	! heights for w exists. *
172	!*************************************************************************
173
174	if (init) then
175
176	! Search for a point with high surface pressure (i.e. not above significant topography)
177	! Then, use this point to construct a reference z profile, to be used at all times
178	!*****************************************************************************
179
180	do jy=0,nymin1
181	do ix=0,nxmin1
182	if (ps(ix,jy,1,n).gt.100000.) then
183	ixm=ix
184	jym=jy
185	goto 3
186	endif
187	end do
188	end do
189	3 continue
190
191
192	tvold=tt2(ixm,jym,1,n)(1.+0.378ew(td2(ixm,jym,1,n))/ &
193	ps(ixm,jym,1,n))
194	pold=ps(ixm,jym,1,n)
195	height(1)=0.
196
197	do kz=2,nuvz
198	pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
199	tv=tth(ixm,jym,kz,n)(1.+0.608qvh(ixm,jym,kz,n))
200
201
202	! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
203	! upon the transformation to z levels. In order to save computer memory, this is
204	! not done anymore in the standard version. However, this option can still be
205	! switched on by replacing the following lines with those below, that are
206	! currently commented out.
207	! Note that two more changes are necessary in this subroutine below.
208	! One change is also necessary in gridcheck.f, and another one in verttransform_nests.
209	!*****************************************************************************
210
211	if (abs(tv-tvold).gt.0.2) then
212	height(kz)= &
213	height(kz-1)+constlog(pold/pint) &
214	(tv-tvold)/log(tv/tvold)
215	else
216	height(kz)=height(kz-1)+ &
217	constlog(pold/pint)tv
218	endif
219
220	! Switch on following lines to use doubled vertical resolution
221	!*************************************************************
222	! if (abs(tv-tvold).gt.0.2) then
223	! height((kz-1)*2)=
224	! + height(max((kz-2)2,1))+constlog(pold/pint)*
225	! + (tv-tvold)/log(tv/tvold)
226	! else
227	! height((kz-1)2)=height(max((kz-2)2,1))+
228	! + constlog(pold/pint)tv
229	! endif
230	! End doubled vertical resolution
231
232	tvold=tv
233	pold=pint
234	end do
235
236
237	! Switch on following lines to use doubled vertical resolution
238	!*************************************************************
239	! do 7 kz=3,nz-1,2
240	! height(kz)=0.5*(height(kz-1)+height(kz+1))
241	! height(nz)=height(nz-1)+height(nz-1)-height(nz-2)
242	! End doubled vertical resolution
243
244
245	! Determine highest levels that can be within PBL
246	!************************************************
247
248	do kz=1,nz
249	if (height(kz).gt.hmixmax) then
250	nmixz=kz
251	goto 9
252	endif
253	end do
254	9 continue
255
256	! Do not repeat initialization of the Cartesian z grid
257	!*****************************************************
258
259	init=.false.
260
261	endif
262
263
264	! Loop over the whole grid
265	!*************************
266
267	do jy=0,nymin1
268	do ix=0,nxmin1
269
270	! NCEP version: find first level above ground
271	llev = 0
272	do i=1,nuvz
273	if (ps(ix,jy,1,n).lt.akz(i)) llev=i
274	end do
275	llev = llev+1
276	if (llev.gt.nuvz-2) llev = nuvz-2
277	! if (llev.eq.nuvz-2) write(,) 'verttransform
278	! +WARNING: LLEV eq NUZV-2'
279	! NCEP version
280
281
282	! compute height of pressure levels above ground
283	!***********************************************
284
285	tvold=tth(ix,jy,llev,n)(1.+0.608qvh(ix,jy,llev,n))
286	pold=akz(llev)
287	uvzlev(llev)=0.
288	wzlev(llev)=0.
289	uvwzlev(ix,jy,llev)=0.
290	rhoh(llev)=pold/(r_air*tvold)
291
292	do kz=llev+1,nuvz
293	pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)
294	tv=tth(ix,jy,kz,n)(1.+0.608qvh(ix,jy,kz,n))
295	rhoh(kz)=pint/(r_air*tv)
296
297	if (abs(tv-tvold).gt.0.2) then
298	uvzlev(kz)=uvzlev(kz-1)+constlog(pold/pint) &
299	(tv-tvold)/log(tv/tvold)
300	else
301	uvzlev(kz)=uvzlev(kz-1)+constlog(pold/pint)tv
302	endif
303	wzlev(kz)=uvzlev(kz)
304	uvwzlev(ix,jy,kz)=uvzlev(kz)
305
306	tvold=tv
307	pold=pint
308	end do
309
310
311	! Switch on following lines to use doubled vertical resolution
312	! Switch off the three lines above.
313	!*************************************************************
314	!22 uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
315	! do 23 kz=2,nwz
316	!23 uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
317	! End doubled vertical resolution
318
319	! pinmconv=(h2-h1)/(p2-p1)
320
321	pinmconv(llev)=(uvwzlev(ix,jy,llev+1)-uvwzlev(ix,jy,llev))/ &
322	((aknew(llev+1)+bknew(llev+1)*ps(ix,jy,1,n))- &
323	(aknew(llev)+bknew(llev)*ps(ix,jy,1,n)))
324	do kz=llev+1,nz-1
325	pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
326	((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
327	(aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
328	end do
329	pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
330	((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
331	(aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
332
333
334	! Levels, where u,v,t and q are given
335	!************************************
336
337	uu(ix,jy,1,n)=uuh(ix,jy,llev)
338	vv(ix,jy,1,n)=vvh(ix,jy,llev)
339	tt(ix,jy,1,n)=tth(ix,jy,llev,n)
340	qv(ix,jy,1,n)=qvh(ix,jy,llev,n)
341	pv(ix,jy,1,n)=pvh(ix,jy,llev)
342	rho(ix,jy,1,n)=rhoh(llev)
343	pplev(ix,jy,1,n)=akz(llev)
344	uu(ix,jy,nz,n)=uuh(ix,jy,nuvz)
345	vv(ix,jy,nz,n)=vvh(ix,jy,nuvz)
346	tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
347	qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
348	pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
349	rho(ix,jy,nz,n)=rhoh(nuvz)
350	pplev(ix,jy,nz,n)=akz(nuvz)
351	kmin=llev+1
352	do iz=2,nz-1
353	do kz=kmin,nuvz
354	if(height(iz).gt.uvzlev(nuvz)) then
355	uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
356	vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
357	tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
358	qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
359	pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
360	rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
361	pplev(ix,jy,iz,n)=pplev(ix,jy,nz,n)
362	goto 30
363	endif
364	if ((height(iz).gt.uvzlev(kz-1)).and. &
365	(height(iz).le.uvzlev(kz))) then
366	dz1=height(iz)-uvzlev(kz-1)
367	dz2=uvzlev(kz)-height(iz)
368	dz=dz1+dz2
369	uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)dz2+uuh(ix,jy,kz)dz1)/dz
370	vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)dz2+vvh(ix,jy,kz)dz1)/dz
371	tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
372	+tth(ix,jy,kz,n)*dz1)/dz
373	qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
374	+qvh(ix,jy,kz,n)*dz1)/dz
375	pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)dz2+pvh(ix,jy,kz)dz1)/dz
376	rho(ix,jy,iz,n)=(rhoh(kz-1)dz2+rhoh(kz)dz1)/dz
377	pplev(ix,jy,iz,n)=(akz(kz-1)dz2+akz(kz)dz1)/dz
378	endif
379	end do
380	30 continue
381	end do
382
383
384	! Levels, where w is given
385	!*************************
386
387	ww(ix,jy,1,n)=wwh(ix,jy,llev)*pinmconv(llev)
388	ww(ix,jy,nz,n)=wwh(ix,jy,nwz)*pinmconv(nz)
389	kmin=llev+1
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	ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)pinmconv(kz-1)dz2 &
398	+wwh(ix,jy,kz)pinmconv(kz)dz1)/dz
399
400	endif
401	end do
402	end do
403
404
405	! Compute density gradients at intermediate levels
406	!*************************************************
407
408	drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
409	(height(2)-height(1))
410	do kz=2,nz-1
411	drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
412	(height(kz+1)-height(kz-1))
413	end do
414	drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
415
416	end do
417	end do
418
419
420	!****************************************************************
421	! Compute slope of eta levels in windward direction and resulting
422	! vertical wind correction
423	!****************************************************************
424
425	do jy=1,ny-2
426	do ix=1,nx-2
427
428	! NCEP version: find first level above ground
429	llev = 0
430	do i=1,nuvz
431	if (ps(ix,jy,1,n).lt.akz(i)) llev=i
432	end do
433	llev = llev+1
434	if (llev.gt.nuvz-2) llev = nuvz-2
435	! if (llev.eq.nuvz-2) write(,) 'verttransform
436	! +WARNING: LLEV eq NUZV-2'
437	! NCEP version
438
439	kmin=llev+1
440	do iz=2,nz-1
441
442	ui=uu(ix,jy,iz,n)dxconst/cos((real(jy)dy+ylat0)*pi180)
443	vi=vv(ix,jy,iz,n)*dyconst
444
445	do kz=kmin,nz
446	if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
447	(height(iz).le.uvwzlev(ix,jy,kz))) then
448	dz1=height(iz)-uvwzlev(ix,jy,kz-1)
449	dz2=uvwzlev(ix,jy,kz)-height(iz)
450	dz=dz1+dz2
451	kl=kz-1
452	klp=kz
453	goto 47
454	endif
455	end do
456
457	47 ix1=ix-1
458	jy1=jy-1
459	ixp=ix+1
460	jyp=jy+1
461
462	dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
463	dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
464	dzdx=(dzdx1dz2+dzdx2dz1)/dz
465
466	dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
467	dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
468	dzdy=(dzdy1dz2+dzdy2dz1)/dz
469
470	ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdxui+dzdyvi)
471
472	end do
473
474	end do
475	end do
476
477
478	! If north pole is in the domain, calculate wind velocities in polar
479	! stereographic coordinates
480	!*******************************************************************
481
482	if (nglobal) then
483	do jy=int(switchnorthg)-2,nymin1
484	ylat=ylat0+real(jy)*dy
485	do ix=0,nxmin1
486	xlon=xlon0+real(ix)*dx
487	do iz=1,nz
488	call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
489	vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
490	vvpol(ix,jy,iz,n))
491	end do
492	end do
493	end do
494
495
496	do iz=1,nz
497
498	! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
499	xlon=xlon0+real(nx/2-1)*dx
500	xlonr=xlon*pi/180.
501	ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
502	vv(nx/2-1,nymin1,iz,n)**2)
503	if(vv(nx/2-1,nymin1,iz,n).lt.0.) then
504	ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
505	vv(nx/2-1,nymin1,iz,n))-xlonr
506	elseif (vv(nx/2-1,nymin1,iz,n).gt.0.) then
507	ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
508	vv(nx/2-1,nymin1,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=nymin1
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=ny-2
538	do ix=0,nxmin1
539	wdummy=wdummy+ww(ix,jy,iz,n)
540	end do
541	wdummy=wdummy/real(nx)
542	jy=nymin1
543	do ix=0,nxmin1
544	ww(ix,jy,iz,n)=wdummy
545	end do
546	end do
547
548	endif
549
550
551	! If south pole is in the domain, calculate wind velocities in polar
552	! stereographic coordinates
553	!*******************************************************************
554
555	if (sglobal) then
556	do jy=0,int(switchsouthg)+3
557	ylat=ylat0+real(jy)*dy
558	do ix=0,nxmin1
559	xlon=xlon0+real(ix)*dx
560	do iz=1,nz
561	call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
562	vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
563	vvpol(ix,jy,iz,n))
564	end do
565	end do
566	end do
567
568	do iz=1,nz
569
570	! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
571	xlon=xlon0+real(nx/2-1)*dx
572	xlonr=xlon*pi/180.
573	ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
574	vv(nx/2-1,0,iz,n)**2)
575	if(vv(nx/2-1,0,iz,n).lt.0.) then
576	ddpol=atan(uu(nx/2-1,0,iz,n)/ &
577	vv(nx/2-1,0,iz,n))+xlonr
578	elseif (vv(nx/2-1,0,iz,n).gt.0.) then
579	ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
580	vv(nx/2-1,0,iz,n))-xlonr
581	else
582	ddpol=pi/2-xlonr
583	endif
584	if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
585	if(ddpol.gt.2.0pi) ddpol=ddpol-2.0pi
586
587	! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
588	xlon=180.0
589	xlonr=xlon*pi/180.
590	ylat=-90.0
591	uuaux=+ffpol*sin(xlonr-ddpol)
592	vvaux=-ffpol*cos(xlonr-ddpol)
593	call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
594	vvpolaux)
595
596	jy=0
597	do ix=0,nxmin1
598	uupol(ix,jy,iz,n)=uupolaux
599	vvpol(ix,jy,iz,n)=vvpolaux
600	end do
601	end do
602
603
604	! Fix: Set W at pole to the zonally averaged W of the next equator-
605	! ward parallel of latitude
606
607	do iz=1,nz
608	wdummy=0.
609	jy=1
610	do ix=0,nxmin1
611	wdummy=wdummy+ww(ix,jy,iz,n)
612	end do
613	wdummy=wdummy/real(nx)
614	jy=0
615	do ix=0,nxmin1
616	ww(ix,jy,iz,n)=wdummy
617	end do
618	end do
619	endif
620
621
622	! write (,) 'initializing clouds, n:',n,nymin1,nxmin1,nz
623	! create a cloud and rainout/washout field, clouds occur where rh>80%
624	! total cloudheight is stored at level 0
625	do jy=0,nymin1
626	do ix=0,nxmin1
627	rain_cloud_above=0
628	lsp=lsprec(ix,jy,1,n)
629	convp=convprec(ix,jy,1,n)
630	cloudsh(ix,jy,n)=0
631	do kz_inv=1,nz-1
632	kz=nz-kz_inv+1
633	pressure=rho(ix,jy,kz,n)r_airtt(ix,jy,kz,n)
634	rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
635	clouds(ix,jy,kz,n)=0
636	if (rh.gt.0.8) then ! in cloud
637	if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
638	rain_cloud_above=1
639	cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
640	height(kz)-height(kz-1)
641	if (lsp.ge.convp) then
642	clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
643	else
644	clouds(ix,jy,kz,n)=2 ! convp dominated rainout
645	endif
646	else ! no precipitation
647	clouds(ix,jy,kz,n)=1 ! cloud
648	endif
649	else ! no cloud
650	if (rain_cloud_above.eq.1) then ! scavenging
651	if (lsp.ge.convp) then
652	clouds(ix,jy,kz,n)=5 ! lsp dominated washout
653	else
654	clouds(ix,jy,kz,n)=4 ! convp dominated washout
655	endif
656	endif
657	endif
658	end do
659	end do
660	end do
661
662
663	end subroutine verttransform_gfs

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

Download in other formats:

Original Format