Context Navigation

verttransform.f90 @ 37

Last change on this file since 37 was 37, checked in by pesei, 9 years ago
Wet dep quick fix and other small changes. Wet depo quick fix not final yet.
File size: 18.2 KB

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

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

Download in other formats:

Original Format