Context Navigation

verttransform_grib2nc4_gfs.f90 @ 94106e2

FPv9.3.2grib2nc4_repair

Last change on this file since 94106e2 was 8624a75, checked in by Don Morton <Don.Morton@…>, 7 years ago

Enhancements to FPv9.3.2

Documented in Ticket #182 (as well as CTBTO ticket 357)

Property mode set to 100644

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

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

Download in other formats:

Original Format