Context Navigation

verttransform_nests.f90

Last change on this file was 13, checked in by saeck, 11 years ago
update to wetdepo.f90
File size: 13.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_nests(n,uuhn,vvhn,wwhn,pvhn)
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	! It is similar to verttransform, but makes the transformations for *
33	! the nested grids. *
34	! *
35	! Author: A. Stohl, G. Wotawa *
36	! *
37	! 12 August 1996 *
38	! Update: 16 January 1998 *
39	! *
40	! Major update: 17 February 1999 *
41	! by G. Wotawa *
42	! *
43	! - Vertical levels for u, v and w are put together *
44	! - Slope correction for vertical velocity: Modification of calculation *
45	! procedure *
46	! *
47	!*****************************************************************************
48	! Changes, Bernd C. Krueger, Feb. 2001: (marked "C-cv")
49	! Variables tthn and qvhn (on eta coordinates) from common block
50	!*****************************************************************************
51	! Sabine Eckhardt, March 2007
52	! add the variable cloud for use with scavenging - descr. in com_mod
53	!*****************************************************************************
54	! *
55	! Variables: *
56	! nxn,nyn,nuvz,nwz field dimensions in x,y and z direction *
57	! uun wind components in x-direction [m/s] *
58	! vvn wind components in y-direction [m/s] *
59	! wwn wind components in z-direction [deltaeta/s]*
60	! ttn temperature [K] *
61	! pvn potential vorticity (pvu) *
62	! psn surface pressure [Pa] *
63	! *
64	!*****************************************************************************
65
66	use par_mod
67	use com_mod
68
69	implicit none
70
71	integer :: ix,jy,kz,iz,n,l,kmin,kl,klp,ix1,jy1,ixp,jyp
72	integer :: rain_cloud_above,kz_inv
73	real :: f_qvsat,pressure,rh,lsp,convp
74	real :: uvzlev(nuvzmax),wzlev(nwzmax),rhoh(nuvzmax),pinmconv(nzmax)
75	real :: uvwzlev(0:nxmaxn-1,0:nymaxn-1,nzmax)
76	real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
77	real :: dzdx,dzdy
78	real :: dzdx1,dzdx2,dzdy1,dzdy2
79	real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
80	real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
81	real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
82	real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
83	real,parameter :: const=r_air/ga
84
85
86	! Loop over all nests
87	!********************
88
89	do l=1,numbnests
90
91	! Loop over the whole grid
92	!*************************
93
94	do jy=0,nyn(l)-1
95	do ix=0,nxn(l)-1
96
97	tvold=tt2n(ix,jy,1,n,l)(1.+0.378ew(td2n(ix,jy,1,n,l))/ &
98	psn(ix,jy,1,n,l))
99	pold=psn(ix,jy,1,n,l)
100	uvzlev(1)=0.
101	wzlev(1)=0.
102	rhoh(1)=pold/(r_air*tvold)
103
104
105	! Compute heights of eta levels
106	!******************************
107
108	do kz=2,nuvz
109	pint=akz(kz)+bkz(kz)*psn(ix,jy,1,n,l)
110	tv=tthn(ix,jy,kz,n,l)(1.+0.608qvhn(ix,jy,kz,n,l))
111	rhoh(kz)=pint/(r_air*tv)
112
113	if (abs(tv-tvold).gt.0.2) then
114	uvzlev(kz)=uvzlev(kz-1)+constlog(pold/pint) &
115	(tv-tvold)/log(tv/tvold)
116	else
117	uvzlev(kz)=uvzlev(kz-1)+constlog(pold/pint)tv
118	endif
119
120	tvold=tv
121	pold=pint
122	end do
123
124
125	do kz=2,nwz-1
126	wzlev(kz)=(uvzlev(kz+1)+uvzlev(kz))/2.
127	end do
128	wzlev(nwz)=wzlev(nwz-1)+ &
129	uvzlev(nuvz)-uvzlev(nuvz-1)
130
131	! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
132	! upon the transformation to z levels. In order to save computer memory, this is
133	! not done anymore in the standard version. However, this option can still be
134	! switched on by replacing the following lines with those below, that are
135	! currently commented out.
136	! Note that one change is also necessary in gridcheck.f,
137	! and three changes in verttransform.f
138	!*****************************************************************************
139	uvwzlev(ix,jy,1)=0.0
140	do kz=2,nuvz
141	uvwzlev(ix,jy,kz)=uvzlev(kz)
142	end do
143
144	! Switch on following lines to use doubled vertical resolution
145	! Switch off the three lines above.
146	!*************************************************************
147	!22 uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
148	! do 23 kz=2,nwz
149	!23 uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
150	! End doubled vertical resolution
151
152	! pinmconv=(h2-h1)/(p2-p1)
153
154	pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
155	((aknew(2)+bknew(2)*psn(ix,jy,1,n,l))- &
156	(aknew(1)+bknew(1)*psn(ix,jy,1,n,l)))
157	do kz=2,nz-1
158	pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
159	((aknew(kz+1)+bknew(kz+1)*psn(ix,jy,1,n,l))- &
160	(aknew(kz-1)+bknew(kz-1)*psn(ix,jy,1,n,l)))
161	end do
162	pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
163	((aknew(nz)+bknew(nz)*psn(ix,jy,1,n,l))- &
164	(aknew(nz-1)+bknew(nz-1)*psn(ix,jy,1,n,l)))
165
166
167	! Levels, where u,v,t and q are given
168	!************************************
169
170	uun(ix,jy,1,n,l)=uuhn(ix,jy,1,l)
171	vvn(ix,jy,1,n,l)=vvhn(ix,jy,1,l)
172	ttn(ix,jy,1,n,l)=tthn(ix,jy,1,n,l)
173	qvn(ix,jy,1,n,l)=qvhn(ix,jy,1,n,l)
174	pvn(ix,jy,1,n,l)=pvhn(ix,jy,1,l)
175	rhon(ix,jy,1,n,l)=rhoh(1)
176	uun(ix,jy,nz,n,l)=uuhn(ix,jy,nuvz,l)
177	vvn(ix,jy,nz,n,l)=vvhn(ix,jy,nuvz,l)
178	ttn(ix,jy,nz,n,l)=tthn(ix,jy,nuvz,n,l)
179	qvn(ix,jy,nz,n,l)=qvhn(ix,jy,nuvz,n,l)
180	pvn(ix,jy,nz,n,l)=pvhn(ix,jy,nuvz,l)
181	rhon(ix,jy,nz,n,l)=rhoh(nuvz)
182	kmin=2
183	do iz=2,nz-1
184	do kz=kmin,nuvz
185	if(height(iz).gt.uvzlev(nuvz)) then
186	uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l)
187	vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l)
188	ttn(ix,jy,iz,n,l)=ttn(ix,jy,nz,n,l)
189	qvn(ix,jy,iz,n,l)=qvn(ix,jy,nz,n,l)
190	pvn(ix,jy,iz,n,l)=pvn(ix,jy,nz,n,l)
191	rhon(ix,jy,iz,n,l)=rhon(ix,jy,nz,n,l)
192	goto 30
193	endif
194	if ((height(iz).gt.uvzlev(kz-1)).and. &
195	(height(iz).le.uvzlev(kz))) then
196	dz1=height(iz)-uvzlev(kz-1)
197	dz2=uvzlev(kz)-height(iz)
198	dz=dz1+dz2
199	uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+ &
200	uuhn(ix,jy,kz,l)*dz1)/dz
201	vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+ &
202	vvhn(ix,jy,kz,l)*dz1)/dz
203	ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2+ &
204	tthn(ix,jy,kz,n,l)*dz1)/dz
205	qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2+ &
206	qvhn(ix,jy,kz,n,l)*dz1)/dz
207	pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+ &
208	pvhn(ix,jy,kz,l)*dz1)/dz
209	rhon(ix,jy,iz,n,l)=(rhoh(kz-1)dz2+rhoh(kz)dz1)/dz
210	kmin=kz
211	goto 30
212	endif
213	end do
214	30 continue
215	end do
216
217
218	! Levels, where w is given
219	!*************************
220
221	wwn(ix,jy,1,n,l)=wwhn(ix,jy,1,l)*pinmconv(1)
222	wwn(ix,jy,nz,n,l)=wwhn(ix,jy,nwz,l)*pinmconv(nz)
223	kmin=2
224	do iz=2,nz
225	do kz=kmin,nwz
226	if ((height(iz).gt.wzlev(kz-1)).and. &
227	(height(iz).le.wzlev(kz))) then
228	dz1=height(iz)-wzlev(kz-1)
229	dz2=wzlev(kz)-height(iz)
230	dz=dz1+dz2
231	wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)pinmconv(kz-1)dz2 &
232	+wwhn(ix,jy,kz,l)pinmconv(kz)dz1)/dz
233	kmin=kz
234	goto 40
235	endif
236	end do
237	40 continue
238	end do
239
240	! Compute density gradients at intermediate levels
241	!*************************************************
242
243	drhodzn(ix,jy,1,n,l)=(rhon(ix,jy,2,n,l)-rhon(ix,jy,1,n,l))/ &
244	(height(2)-height(1))
245	do kz=2,nz-1
246	drhodzn(ix,jy,kz,n,l)=(rhon(ix,jy,kz+1,n,l)- &
247	rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
248	end do
249	drhodzn(ix,jy,nz,n,l)=drhodzn(ix,jy,nz-1,n,l)
250
251	end do
252	end do
253
254
255	!****************************************************************
256	! Compute slope of eta levels in windward direction and resulting
257	! vertical wind correction
258	!****************************************************************
259
260	do jy=1,nyn(l)-2
261	do ix=1,nxn(l)-2
262
263	kmin=2
264	do iz=2,nz-1
265
266	ui=uun(ix,jy,iz,n,l)dxconstxresoln(l)/ &
267	cos((real(jy)dyn(l)+ylat0n(l))pi180)
268	vi=vvn(ix,jy,iz,n,l)dyconstyresoln(l)
269
270	do kz=kmin,nz
271	if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
272	(height(iz).le.uvwzlev(ix,jy,kz))) then
273	dz1=height(iz)-uvwzlev(ix,jy,kz-1)
274	dz2=uvwzlev(ix,jy,kz)-height(iz)
275	dz=dz1+dz2
276	kl=kz-1
277	klp=kz
278	kmin=kz
279	goto 47
280	endif
281	end do
282
283	47 ix1=ix-1
284	jy1=jy-1
285	ixp=ix+1
286	jyp=jy+1
287
288	dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
289	dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
290	dzdx=(dzdx1dz2+dzdx2dz1)/dz
291
292	dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
293	dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
294	dzdy=(dzdy1dz2+dzdy2dz1)/dz
295
296	wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l)+(dzdxui+dzdyvi)
297
298	end do
299
300	end do
301	end do
302
303
304	!write (,) 'initializing nested cloudsn, n:',n
305	! create a cloud and rainout/washout field, cloudsn occur where rh>80%
306	do jy=0,nyn(l)-1
307	do ix=0,nxn(l)-1
308	rain_cloud_above=0
309	lsp=lsprecn(ix,jy,1,n,l)
310	convp=convprecn(ix,jy,1,n,l)
311	cloudsnh(ix,jy,n,l)=0
312	do kz_inv=1,nz-1
313	kz=nz-kz_inv+1
314	pressure=rhon(ix,jy,kz,n,l)r_airttn(ix,jy,kz,n,l)
315	rh=qvn(ix,jy,kz,n,l)/f_qvsat(pressure,ttn(ix,jy,kz,n,l))
316	cloudsn(ix,jy,kz,n,l)=0
317	if (rh.gt.0.8) then ! in cloud
318	if ((lsp.gt.0.01).or.(convp.gt.0.01)) then
319	rain_cloud_above=1
320	cloudsnh(ix,jy,n,l)=cloudsnh(ix,jy,n,l)+ &
321	height(kz)-height(kz-1)
322	if (lsp.ge.convp) then
323	cloudsn(ix,jy,kz,n,l)=3 ! lsp dominated rainout
324	else
325	cloudsn(ix,jy,kz,n,l)=2 ! convp dominated rainout
326	endif
327	else ! no precipitation
328	cloudsn(ix,jy,kz,n,l)=1 ! cloud
329	endif
330	else ! no cloud
331	if (rain_cloud_above.eq.1) then ! scavenging
332	if (lsp.ge.convp) then
333	cloudsn(ix,jy,kz,n,l)=5 ! lsp dominated washout
334	else
335	cloudsn(ix,jy,kz,n,l)=4 ! convp dominated washout
336	endif
337	endif
338	endif
339	end do
340	end do
341	end do
342
343	end do
344
345	end subroutine verttransform_nests

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

Download in other formats:

Original Format