Context Navigation

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

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

Download in other formats:

Original Format