Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

verttransform.f90 @ 23

Last change on this file since 23 was 23, checked in by igpis, 10 years ago
start tracking test environment directory FP_AI
File size: 17.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(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	! *
53	! Variables: *
54	! nx,ny,nz field dimensions in x,y and z direction *
55	! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition *
56	! uu(0:nxmax,0:nymax,nzmax,2) wind components in x-direction [m/s] *
57	! vv(0:nxmax,0:nymax,nzmax,2) wind components in y-direction [m/s] *
58	! ww(0:nxmax,0:nymax,nzmax,2) wind components in z-direction [deltaeta/s]*
59	! tt(0:nxmax,0:nymax,nzmax,2) temperature [K] *
60	! pv(0:nxmax,0:nymax,nzmax,2) potential voriticity (pvu) *
61	! ps(0:nxmax,0:nymax,2) surface pressure [Pa] *
62	! *
63	!*****************************************************************************
64
65	use par_mod
66	use com_mod
67	use cmapf_mod, only: cc2gll
68
69	implicit none
70
71	integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
72	integer :: rain_cloud_above,kz_inv
73	real :: f_qvsat,pressure
74	real :: rh,lsp,convp
75	real :: rhoh(nuvzmax),pinmconv(nzmax)
76	real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
77	real :: xlon,ylat,xlonr,dzdx,dzdy
78	real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf
79	real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
80	real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
81	real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
82	real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
83	real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
84	real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
85	real,parameter :: const=r_air/ga
86
87	logical :: init = .true.
88
89
90	!*************************************************************************
91	! If verttransform is called the first time, initialize heights of the *
92	! z levels in meter. The heights are the heights of model levels, where *
93	! u,v,T and qv are given. *
94	!*************************************************************************
95
96	if (init) then
97
98	! Search for a point with high surface pressure (i.e. not above significant topography)
99	! Then, use this point to construct a reference z profile, to be used at all times
100	!**************************************************************************************
101
102	do jy=0,nymin1
103	do ix=0,nxmin1
104	if (ps(ix,jy,1,n).gt.100000.) then
105	ixm=ix
106	jym=jy
107	goto 3
108	endif
109	end do
110	end do
111	3 continue
112
113
114	tvold=tt2(ixm,jym,1,n)(1.+0.378ew(td2(ixm,jym,1,n))/ &
115	ps(ixm,jym,1,n))
116	pold=ps(ixm,jym,1,n)
117	height(1)=0.
118
119	do kz=2,nuvz
120	pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
121	tv=tth(ixm,jym,kz,n)(1.+0.608qvh(ixm,jym,kz,n))
122
123	if (abs(tv-tvold).gt.0.2) then
124	height(kz)=height(kz-1)+constlog(pold/pint) &
125	(tv-tvold)/log(tv/tvold)
126	else
127	height(kz)=height(kz-1)+constlog(pold/pint)tv
128	endif
129
130	tvold=tv
131	pold=pint
132	end do
133
134
135	! Determine highest levels that can be within PBL
136	!************************************************
137
138	do kz=1,nz
139	if (height(kz).gt.hmixmax) then
140	nmixz=kz
141	goto 9
142	endif
143	end do
144	9 continue
145
146	! Do not repeat initialization of the Cartesian z grid
147	!*****************************************************
148
149	init=.false.
150
151	endif
152
153
154	! Loop over the whole grid
155	!*************************
156
157	do jy=0,nymin1
158	do ix=0,nxmin1
159	tvold=tt2(ix,jy,1,n)(1.+0.378ew(td2(ix,jy,1,n))/ps(ix,jy,1,n))
160	pold=ps(ix,jy,1,n)
161	uvwzlev(ix,jy,1)=0.
162	wzlev(1)=0.
163	rhoh(1)=pold/(r_air*tvold)
164
165
166	! Compute heights of eta levels
167	!******************************
168
169	do kz=2,nuvz
170	pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)
171	tv=tth(ix,jy,kz,n)(1.+0.608qvh(ix,jy,kz,n))
172	rhoh(kz)=pint/(r_air*tv)
173
174	if (abs(tv-tvold).gt.0.2) then
175	uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+constlog(pold/pint) &
176	(tv-tvold)/log(tv/tvold)
177	else
178	uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+constlog(pold/pint)tv
179	endif
180
181	tvold=tv
182	pold=pint
183	end do
184
185
186	do kz=2,nwz-1
187	wzlev(kz)=(uvwzlev(ix,jy,kz+1)+uvwzlev(ix,jy,kz))/2.
188	end do
189	wzlev(nwz)=wzlev(nwz-1)+uvwzlev(ix,jy,nuvz)-uvwzlev(ix,jy,nuvz-1)
190
191	! pinmconv=(h2-h1)/(p2-p1)
192
193	pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
194	((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
195	(aknew(1)+bknew(1)*ps(ix,jy,1,n)))
196	do kz=2,nz-1
197	pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
198	((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
199	(aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
200	end do
201	pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
202	((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
203	(aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
204
205	! Levels, where u,v,t and q are given
206	!************************************
207
208	uu(ix,jy,1,n)=uuh(ix,jy,1)
209	vv(ix,jy,1,n)=vvh(ix,jy,1)
210	tt(ix,jy,1,n)=tth(ix,jy,1,n)
211	qv(ix,jy,1,n)=qvh(ix,jy,1,n)
212	pv(ix,jy,1,n)=pvh(ix,jy,1)
213	rho(ix,jy,1,n)=rhoh(1)
214	uu(ix,jy,nz,n)=uuh(ix,jy,nuvz)
215	vv(ix,jy,nz,n)=vvh(ix,jy,nuvz)
216	tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
217	qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
218	pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
219	rho(ix,jy,nz,n)=rhoh(nuvz)
220	kmin=2
221	do iz=2,nz-1
222	do kz=kmin,nuvz
223	if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
224	uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
225	vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
226	tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
227	qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
228	pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
229	rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
230	goto 30
231	endif
232	if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
233	(height(iz).le.uvwzlev(ix,jy,kz))) then
234	dz1=height(iz)-uvwzlev(ix,jy,kz-1)
235	dz2=uvwzlev(ix,jy,kz)-height(iz)
236	dz=dz1+dz2
237	uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)dz2+uuh(ix,jy,kz)dz1)/dz
238	vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)dz2+vvh(ix,jy,kz)dz1)/dz
239	tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
240	+tth(ix,jy,kz,n)*dz1)/dz
241	qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
242	+qvh(ix,jy,kz,n)*dz1)/dz
243	pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)dz2+pvh(ix,jy,kz)dz1)/dz
244	rho(ix,jy,iz,n)=(rhoh(kz-1)dz2+rhoh(kz)dz1)/dz
245	kmin=kz
246	goto 30
247	endif
248	end do
249	30 continue
250	end do
251
252
253	! Levels, where w is given
254	!*************************
255
256	ww(ix,jy,1,n)=wwh(ix,jy,1)*pinmconv(1)
257	ww(ix,jy,nz,n)=wwh(ix,jy,nwz)*pinmconv(nz)
258	kmin=2
259	do iz=2,nz
260	do kz=kmin,nwz
261	if ((height(iz).gt.wzlev(kz-1)).and. &
262	(height(iz).le.wzlev(kz))) then
263	dz1=height(iz)-wzlev(kz-1)
264	dz2=wzlev(kz)-height(iz)
265	dz=dz1+dz2
266	ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)pinmconv(kz-1)dz2 &
267	+wwh(ix,jy,kz)pinmconv(kz)dz1)/dz
268	kmin=kz
269	goto 40
270	endif
271	end do
272	40 continue
273	end do
274
275	! Compute density gradients at intermediate levels
276	!*************************************************
277
278	drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
279	(height(2)-height(1))
280	do kz=2,nz-1
281	drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
282	(height(kz+1)-height(kz-1))
283	end do
284	drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
285
286	end do
287	end do
288
289
290	!****************************************************************
291	! Compute slope of eta levels in windward direction and resulting
292	! vertical wind correction
293	!****************************************************************
294
295	do jy=1,ny-2
296	cosf=cos((real(jy)dy+ylat0)pi180)
297	do ix=1,nx-2
298
299	kmin=2
300	do iz=2,nz-1
301
302	ui=uu(ix,jy,iz,n)*dxconst/cosf
303	vi=vv(ix,jy,iz,n)*dyconst
304
305	do kz=kmin,nz
306	if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
307	(height(iz).le.uvwzlev(ix,jy,kz))) then
308	dz1=height(iz)-uvwzlev(ix,jy,kz-1)
309	dz2=uvwzlev(ix,jy,kz)-height(iz)
310	dz=dz1+dz2
311	kl=kz-1
312	klp=kz
313	kmin=kz
314	goto 47
315	endif
316	end do
317
318	47 ix1=ix-1
319	jy1=jy-1
320	ixp=ix+1
321	jyp=jy+1
322
323	dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
324	dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
325	dzdx=(dzdx1dz2+dzdx2dz1)/dz
326
327	dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
328	dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
329	dzdy=(dzdy1dz2+dzdy2dz1)/dz
330
331	ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdxui+dzdyvi)
332
333	end do
334
335	end do
336	end do
337
338
339	! If north pole is in the domain, calculate wind velocities in polar
340	! stereographic coordinates
341	!*******************************************************************
342
343	if (nglobal) then
344	do jy=int(switchnorthg)-2,nymin1
345	ylat=ylat0+real(jy)*dy
346	do ix=0,nxmin1
347	xlon=xlon0+real(ix)*dx
348	do iz=1,nz
349	call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
350	vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
351	end do
352	end do
353	end do
354
355
356	do iz=1,nz
357
358	! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
359	!
360	! AMSnauffer Nov 18 2004 Added check for case vv=0
361	!
362	xlon=xlon0+real(nx/2-1)*dx
363	xlonr=xlon*pi/180.
364	ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)2+vv(nx/2-1,nymin1,iz,n)2)
365	if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
366	ddpol=atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr
367	else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
368	ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
369	vv(nx/2-1,nymin1,iz,n))-xlonr
370	else
371	ddpol=pi/2-xlonr
372	endif
373	if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
374	if(ddpol.gt.2.0pi) ddpol=ddpol-2.0pi
375
376	! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
377	xlon=180.0
378	xlonr=xlon*pi/180.
379	ylat=90.0
380	uuaux=-ffpol*sin(xlonr+ddpol)
381	vvaux=-ffpol*cos(xlonr+ddpol)
382	call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
383	jy=nymin1
384	do ix=0,nxmin1
385	uupol(ix,jy,iz,n)=uupolaux
386	vvpol(ix,jy,iz,n)=vvpolaux
387	end do
388	end do
389
390
391	! Fix: Set W at pole to the zonally averaged W of the next equator-
392	! ward parallel of latitude
393
394	do iz=1,nz
395	wdummy=0.
396	jy=ny-2
397	do ix=0,nxmin1
398	wdummy=wdummy+ww(ix,jy,iz,n)
399	end do
400	wdummy=wdummy/real(nx)
401	jy=nymin1
402	do ix=0,nxmin1
403	ww(ix,jy,iz,n)=wdummy
404	end do
405	end do
406
407	endif
408
409
410	! If south pole is in the domain, calculate wind velocities in polar
411	! stereographic coordinates
412	!*******************************************************************
413
414	if (sglobal) then
415	do jy=0,int(switchsouthg)+3
416	ylat=ylat0+real(jy)*dy
417	do ix=0,nxmin1
418	xlon=xlon0+real(ix)*dx
419	do iz=1,nz
420	call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
421	vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
422	end do
423	end do
424	end do
425
426	do iz=1,nz
427
428	! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
429	!
430	! AMSnauffer Nov 18 2004 Added check for case vv=0
431	!
432	xlon=xlon0+real(nx/2-1)*dx
433	xlonr=xlon*pi/180.
434	ffpol=sqrt(uu(nx/2-1,0,iz,n)2+vv(nx/2-1,0,iz,n)2)
435	if (vv(nx/2-1,0,iz,n).lt.0.) then
436	ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
437	else if (vv(nx/2-1,0,iz,n).gt.0.) then
438	ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
439	else
440	ddpol=pi/2-xlonr
441	endif
442	if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
443	if(ddpol.gt.2.0pi) ddpol=ddpol-2.0pi
444
445	! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
446	xlon=180.0
447	xlonr=xlon*pi/180.
448	ylat=-90.0
449	uuaux=+ffpol*sin(xlonr-ddpol)
450	vvaux=-ffpol*cos(xlonr-ddpol)
451	call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
452
453	jy=0
454	do ix=0,nxmin1
455	uupol(ix,jy,iz,n)=uupolaux
456	vvpol(ix,jy,iz,n)=vvpolaux
457	end do
458	end do
459
460
461	! Fix: Set W at pole to the zonally averaged W of the next equator-
462	! ward parallel of latitude
463
464	do iz=1,nz
465	wdummy=0.
466	jy=1
467	do ix=0,nxmin1
468	wdummy=wdummy+ww(ix,jy,iz,n)
469	end do
470	wdummy=wdummy/real(nx)
471	jy=0
472	do ix=0,nxmin1
473	ww(ix,jy,iz,n)=wdummy
474	end do
475	end do
476	endif
477
478
479	!write (,) 'initializing clouds, n:',n,nymin1,nxmin1,nz
480	! create a cloud and rainout/washout field, clouds occur where rh>80%
481	! total cloudheight is stored at level 0
482	do jy=0,nymin1
483	do ix=0,nxmin1
484	rain_cloud_above=0
485	lsp=lsprec(ix,jy,1,n)
486	convp=convprec(ix,jy,1,n)
487	cloudsh(ix,jy,n)=0
488	do kz_inv=1,nz-1
489	kz=nz-kz_inv+1
490	pressure=rho(ix,jy,kz,n)r_airtt(ix,jy,kz,n)
491	rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
492	clouds(ix,jy,kz,n)=0
493	if (rh.gt.0.8) then ! in cloud
494	if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
495	rain_cloud_above=1
496	cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
497	height(kz)-height(kz-1)
498	if (lsp.ge.convp) then
499	clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
500	else
501	clouds(ix,jy,kz,n)=2 ! convp dominated rainout
502	endif
503	else ! no precipitation
504	clouds(ix,jy,kz,n)=1 ! cloud
505	endif
506	else ! no cloud
507	if (rain_cloud_above.eq.1) then ! scavenging
508	if (lsp.ge.convp) then
509	clouds(ix,jy,kz,n)=5 ! lsp dominated washout
510	else
511	clouds(ix,jy,kz,n)=4 ! convp dominated washout
512	endif
513	endif
514	endif
515	end do
516	end do
517	end do
518
519	end subroutine verttransform

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

Download in other formats:

Original Format