Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

verttransform.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 11 years ago
sources for flexwrf v3.1
File size: 32.2 KB

Line
1	!***********************************************************************
2	!* Copyright 2012,2013 *
3	!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine, *
4	!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
5	!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso, *
6	!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann, *
7	!* *
8	!* This file is part of FLEXPART WRF *
9	!* *
10	!* FLEXPART is free software: you can redistribute it and/or modify *
11	!* it under the terms of the GNU General Public License as published by*
12	!* the Free Software Foundation, either version 3 of the License, or *
13	!* (at your option) any later version. *
14	!* *
15	!* FLEXPART is distributed in the hope that it will be useful, *
16	!* but WITHOUT ANY WARRANTY; without even the implied warranty of *
17	!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
18	!* GNU General Public License for more details. *
19	!* *
20	!* You should have received a copy of the GNU General Public License *
21	!* along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
22	!***********************************************************************
23
24	subroutine verttransform(n,uuh,vvh,wwh,pvh,divh)
25	! i i i i i
26	!*******************************************************************************
27	! *
28	! This subroutine transforms temperature, dew point temperature and *
29	! wind components from eta to meter coordinates. *
30	! The vertical wind component is transformed from Pa/s to m/s using *
31	! the conversion factor pinmconv. *
32	! In addition, this routine calculates vertical density gradients *
33	! needed for the parameterization of the turbulent velocities. *
34	! *
35	! Note: This is the FLEXPART_WRF version of subroutine assignland. *
36	! The computational grid is the WRF x-y grid rather than lat-lon. *
37	! *
38	! Author: A. Stohl, G. Wotawa *
39	! *
40	! 12 August 1996 *
41	! Update: 16 January 1998 *
42	! *
43	! Major update: 17 February 1999 *
44	! by G. Wotawa *
45	! *
46	! - Vertical levels for u, v and w are put together *
47	! - Slope correction for vertical velocity: Modification of calculation *
48	! procedure *
49	! *
50	! Changes, Bernd C. Krueger, Feb. 2001: *
51	! Variables tth and qvh (on eta coordinates) from common block *
52	! *
53	! Oct-Nov 2005 - R. Easter - conversion to wrf *
54	! 17 Nov 2005 - R. Easter - terrain correction applied to ww. There are *
55	! now 3 options, controlled by "method_w_terrain_correction" *
56	! *
57	! 11 June 2007, conversion of tkeh to tke
58	! 25 June 2007 conversion of ptth to ptt
59	! Jan 2012, J Brioude: modified to handle different wind options and openmp
60	!*******************************************************************************
61	! *
62	! Variables: *
63	! nx,ny,nz field dimensions in x,y and z direction *
64	! uu(0:nxmax,0:nymax,nzmax,2) wind components in x-direction [m/s] *
65	! vv(0:nxmax,0:nymax,nzmax,2) wind components in y-direction [m/s] *
66	! ww(0:nxmax,0:nymax,nzmax,2) wind components in z-direction [deltaeta/s] *
67	! tt(0:nxmax,0:nymax,nzmax,2) temperature [K] *
68	! pv(0:nxmax,0:nymax,nzmax,2) potential voriticity (pvu) *
69	! ps(0:nxmax,0:nymax,2) surface pressure [Pa] *
70	! *
71	!*******************************************************************************
72
73	use par_mod
74	use com_mod
75	! include 'includepar'
76	! include 'includecom'
77
78	implicit none
79
80	integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
81	! CDA
82	integer :: icloudtop
83
84	integer :: method_z_compute,aa
85	real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
86	real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
87	real :: xlon,ylat,xlonr,dzdx,dzdy
88	real :: dzdx1,dzdx2,dzdy1,dzdy2
89	real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
90	real(kind=4) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
91	real(kind=4) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
92	real(kind=4) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
93	real(kind=4) :: divh(0:nxmax-1,0:nymax-1,nuvzmax)
94
95	! real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
96	! real :: divh(0:nxmax-1,0:nymax-1,nuvzmax)
97	real :: div(0:nxmax-1,0:nymax-1,nuvzmax)
98	! real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
99	real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
100	! real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
101	real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
102	real :: wwh_svaa(nwzmax), wtc_stat(4,nzmax),u,v
103	real,parameter :: const=r_air/ga
104	! CDA cloud commented
105	! integer :: rain_cloud_above,kz_inv
106
107	integer :: kz_inv
108	real :: f_qvsat,pressure
109	! CDA some new declarations and mods
110	! real :: rh,lsp,convp
111	real :: rh,lsp,convp,prec,rhmin
112	real,parameter :: precmin = 0.002
113
114
115
116	logical :: init = .true.
117	! CDA
118	logical :: lconvectprec = .true.
119
120
121
122	! set method_w_terrain_correction & method_z_compute
123	method_w_terrain_correction = 20
124	method_z_compute = 10
125	aa=0
126	do iz = 1, nz
127	do ix = 1, 4
128	wtc_stat(ix,iz) = 0.0
129	end do
130	end do
131
132
133	!*************************************************************************
134	! If verttransform is called the first time, initialize heights of the *
135	! z levels in meter. The heights are the heights of model levels, where *
136	! u,v,T and qv are given, and of the interfaces, where w is given. So, *
137	! the vertical resolution in the z system is doubled. As reference point,*
138	! the lower left corner of the grid is used. *
139	! Unlike in the eta system, no difference between heights for u,v and *
140	! heights for w exists. *
141	!*************************************************************************
142
143	if (init) then
144
145	! Search for a point with high surface pressure (i.e. not above significant topography)
146	! Then, use this point to construct a reference z profile, to be used at all times
147	!
148	! FLEXPART_WRF - use grid point with highest surface pressure
149	!**************************************************************************************
150
151	pint = -1.0
152	ixm = -999888777
153	jym = -999888777
154	do jy=0,nymin1
155	do ix=0,nxmin1
156	! if (ps(ix,jy,1,n).gt.100000.) then
157	if (ps(ix,jy,1,n).gt.pint) then
158	pint = ps(ix,jy,1,n)
159	ixm=ix
160	jym=jy
161	! goto 3
162	endif
163	enddo
164	enddo
165	3 continue
166	! write(*,'(/a,2i4,1pe11.2)')
167	! & 'verttransform -- ixm,jym,ps() =', ixm, jym, pint
168
169
170	tvold=tt2(ixm,jym,1,n)(1.+0.378ew(td2(ixm,jym,1,n))/ &
171	ps(ixm,jym,1,n))
172	pold=ps(ixm,jym,1,n)
173	height(1)=0.
174
175	do kz=2,nuvz
176	! use pressure from wrf met file
177	! pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
178	pint=pph(ixm,jym,kz,n)
179	tv=tth(ixm,jym,kz,n)(1.+0.608qvh(ixm,jym,kz,n))
180
181
182	! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
183	! upon the transformation to z levels. In order to save computer memory, this is
184	! not done anymore in the standard version. However, this option can still be
185	! switched on by replacing the following lines with those below, that are
186	! currently commented out.
187	! Note that two more changes are necessary in this subroutine below.
188	! One change is also necessary in gridcheck.f, and another one in verttransform_nests.
189	!*************************************************************************************
190
191	if (abs(tv-tvold).gt.0.2) then
192	height(kz)= &
193	height(kz-1)+constlog(pold/pint) &
194	(tv-tvold)/log(tv/tvold)
195	else
196	height(kz)=height(kz-1)+ &
197	constlog(pold/pint)tv
198	endif
199
200	!
201	! *** NOTE -- the doubled vertical resolution has not been tested in FLEXPART_WRF
202	!
203	! Switch on following lines to use doubled vertical resolution
204	!*************************************************************
205	! if (abs(tv-tvold).gt.0.2) then
206	! height((kz-1)*2)=
207	! + height(max((kz-2)2,1))+constlog(pold/pint)*
208	! + (tv-tvold)/log(tv/tvold)
209	! else
210	! height((kz-1)2)=height(max((kz-2)2,1))+
211	! + constlog(pold/pint)tv
212	! endif
213	! End doubled vertical resolution
214
215	! FLEXPART_WRF - get height from zzh
216	if (method_z_compute .eq. 10) then
217	if ((add_sfc_level .eq. 1) .and. (kz .eq. 2)) then
218	height(kz) = 0.5*(zzh(ixm,jym, 3,n)+zzh(ixm,jym, 1,n)) &
219	- zzh(ixm,jym,1,n)
220	else
221	height(kz) = 0.5*(zzh(ixm,jym,kz+1,n)+zzh(ixm,jym,kz,n)) &
222	- zzh(ixm,jym,1,n)
223	end if
224	end if
225
226	tvold=tv
227	pold=pint
228	enddo
229	do kz=1,nz-1
230	heightmid(kz)=0.5*(height(kz)+height(kz+1))
231	enddo
232	!
233	! *** NOTE -- the doubled vertical resolution has not been tested in FLEXPART_WRF
234	!
235	! Switch on following lines to use doubled vertical resolution
236	!*************************************************************
237	! do 7 kz=3,nz-1,2
238	! height(kz)=0.5*(height(kz-1)+height(kz+1))
239	! height(nz)=height(nz-1)+height(nz-1)-height(nz-2)
240	! End doubled vertical resolution
241
242
243	! Determine highest levels that can be within PBL
244	!************************************************
245
246	do kz=1,nz
247	if (height(kz).gt.hmixmax) then
248	nmixz=kz
249	goto 9
250	endif
251	end do
252	9 continue
253
254	! Do not repeat initialization of the Cartesian z grid
255	!*****************************************************
256
257	init=.false.
258
259	endif
260
261
262	! Loop over the whole grid
263	!*************************
264
265	!!!$OMP PARALLEL DEFAULT(SHARED) &
266	!!!$OMP PRIVATE(ix,jy,ixm,jym,tvold,pold,pint,tv,rhoh,uvzlev,wzlev, &
267	!!!$OMP uvwzlev,pinmconv,kz,iz,kmin,dz1,dz2,dz,ix1,jy1,ixp,jyp, &
268	!!!$OMP dzdy,dzdx,aa,u,v,wwh_svaa )
269	!!!$OMP DO
270	do jy=0,nymin1
271	do ix=0,nxmin1
272	tvold=tt2(ix,jy,1,n)(1.+0.378ew(td2(ix,jy,1,n))/ &
273	ps(ix,jy,1,n))
274	pold=ps(ix,jy,1,n)
275	uvzlev(1)=0.
276	wzlev(1)=0.
277	rhoh(1)=pold/(r_air*tvold)
278
279
280	! Compute heights of eta levels
281	!******************************
282
283	do kz=2,nuvz
284	! use pressure from wrf met file
285	! pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)
286	pint=pph(ix,jy,kz,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
306	wzlev(nwz)=wzlev(nwz-1)+ &
307	uvzlev(nuvz)-uvzlev(nuvz-1)
308
309	! FLEXPART_WRF - get uvzlev & wzlev from zzh
310	if (method_z_compute .eq. 10) then
311	do kz = 2, nuvz
312	if ((add_sfc_level .eq. 1) .and. (kz .eq. 2)) then
313	uvzlev(kz) = 0.5*(zzh(ix,jy, 3,n) + zzh(ix,jy, 1,n)) &
314	- zzh(ix,jy,1,n)
315	else
316	uvzlev(kz) = 0.5*(zzh(ix,jy,kz+1,n) + zzh(ix,jy,kz,n)) &
317	- zzh(ix,jy,1,n)
318	end if
319	end do
320	do kz = 2, nwz
321	wzlev(kz) = zzh(ix,jy,kz+add_sfc_level,n) &
322	- zzh(ix,jy,1,n)
323	end do
324	end if
325
326	uvwzlev(ix,jy,1)=0.0
327	do kz=2,nuvz
328	uvwzlev(ix,jy,kz)=uvzlev(kz)
329	end do
330
331	! if ((ix .eq. ixm) .and. (jy .eq. jym)) then
332	! write(*,'(/a)')
333	! & 'kz, height, uvzlev, wzlev, zzh-zzh(1) at ixm,jym (in km)'
334	! write(,'(i3,4f8.3)') (kz, height(kz)1.0e-3,
335	! & uvzlev(kz)1.0e-3, wzlev(kz)1.0e-3,
336	! & (zzh(ix,jy,kz,n)-zzh(ix,jy,1,n))*1.0e-3, kz=nz,1,-1)
337	! ixm = -9
338	! end if
339
340	! Switch on following lines to use doubled vertical resolution
341	! Switch off the three lines above.
342	!*************************************************************
343	!22 uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
344	! do 23 kz=2,nwz
345	!23 uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
346	! End doubled vertical resolution
347
348	! pinmconv=(h2-h1)/(p2-p1)
349	!
350	! in flexpart_ecmwf, pinmconv is used to convert etadot to w
351	! in FLEXPART_WRF, vertical velocity is already m/s, so pinmconv=1.0
352
353	if (wind_option.le.0) then
354	pinmconv(1)=1.0
355	do kz=2,nz-1
356	pinmconv(kz)=1.0
357	enddo
358	pinmconv(nz)=1.0
359	elseif (wind_option.ge.1) then
360
361	! pinmconv(1)=(uvzlev(1+add_sfc_level)-0.) &
362	! /(eta_u_wrf(1)-1.)
363	pinmconv(1)=(wzlev(2)-0.) &
364	/(eta_w_wrf(2)-1.)
365	do kz=2,nz-1
366
367	! pinmconv(kz)=(uvzlev(kz+add_sfc_level)-uvzlev(kz-1+add_sfc_level)) &
368	! /(eta_u_wrf(kz)-eta_u_wrf(kz-1))
369	! /(pph(ix,jy,kz+add_sfc_level,n)-pph(ix,jy,kz-1+add_sfc_level,n)) &
370	! *(pph(ix,jy,1,n)-pph(ix,jy,nz,n))
371	! *(ps(ix,jy,1,n)-p_top_wrf)
372
373	pinmconv(kz)=(wzlev(kz+1)-wzlev(kz-1)) &
374	/(eta_w_wrf(kz+1)-eta_w_wrf(kz-1))
375	enddo
376
377	pinmconv(nwz)=pinmconv(nwz-1) !
378	endif
379	! Levels, where u,v,t and q are given
380	!************************************
381
382	uu(ix,jy,1,n)=uuh(ix,jy,1)
383	vv(ix,jy,1,n)=vvh(ix,jy,1)
384	div(ix,jy,1)=divh(ix,jy,1)
385	tt(ix,jy,1,n)=tth(ix,jy,1,n)
386	qv(ix,jy,1,n)=qvh(ix,jy,1,n)
387	pv(ix,jy,1,n)=pvh(ix,jy,1)
388	rho(ix,jy,1,n)=rhoh(1)
389	uu(ix,jy,nz,n)=uuh(ix,jy,nuvz)
390	vv(ix,jy,nz,n)=vvh(ix,jy,nuvz)
391	tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
392	qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
393	pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
394	rho(ix,jy,nz,n)=rhoh(nuvz)
395	tke(ix,jy,1,n)=tkeh(ix,jy,1,n)
396	tke(ix,jy,nz,n)=tkeh(ix,jy,nuvz,n)
397	ptt(ix,jy,1,n)=ptth(ix,jy,1,n)
398	ptt(ix,jy,nz,n)=ptth(ix,jy,nuvz,n)
399
400
401	kmin=2
402	do iz=2,nz-1
403	do kz=kmin,nuvz
404	if(heightmid(iz).gt.uvzlev(nuvz)) then
405	div(ix,jy,iz)=div(ix,jy,nz)
406	goto 230
407	endif
408	if ((heightmid(iz).gt.uvzlev(kz-1)).and. &
409	(heightmid(iz).le.uvzlev(kz))) then
410	dz1=heightmid(iz)-uvzlev(kz-1)
411	dz2=uvzlev(kz)-heightmid(iz)
412	dz=dz1+dz2
413	div(ix,jy,iz)=(divh(ix,jy,kz-1)dz2+divh(ix,jy,kz)dz1)/dz
414	kmin=kz
415	goto 230
416	endif
417	end do
418	230 continue
419	end do
420
421	kmin=2
422	do iz=2,nz-1
423	do kz=kmin,nuvz
424	if(height(iz).gt.uvzlev(nuvz)) then
425	uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
426	vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
427	tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
428	qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
429	pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
430	rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
431	tke(ix,jy,iz,n)=tke(ix,jy,nz,n)
432	ptt(ix,jy,iz,n)=ptt(ix,jy,nz,n)
433
434	goto 30
435	endif
436	if ((height(iz).gt.uvzlev(kz-1)).and. &
437	(height(iz).le.uvzlev(kz))) then
438	dz1=height(iz)-uvzlev(kz-1)
439	dz2=uvzlev(kz)-height(iz)
440	dz=dz1+dz2
441	uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)dz2+uuh(ix,jy,kz)dz1)/dz
442	vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)dz2+vvh(ix,jy,kz)dz1)/dz
443	tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
444	+tth(ix,jy,kz,n)*dz1)/dz
445	qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
446	+qvh(ix,jy,kz,n)*dz1)/dz
447	pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)dz2+pvh(ix,jy,kz)dz1)/dz
448	rho(ix,jy,iz,n)=(rhoh(kz-1)dz2+rhoh(kz)dz1)/dz
449	tke(ix,jy,iz,n)=(tkeh(ix,jy,kz-1,n)*dz2 &
450	+tkeh(ix,jy,kz,n)*dz1)/dz
451	ptt(ix,jy,iz,n)=(ptth(ix,jy,kz-1,n)*dz2 &
452	+ptth(ix,jy,kz,n)*dz1)/dz
453
454
455	kmin=kz
456	goto 30
457	endif
458	end do
459	30 continue
460	end do
461
462
463	! Levels, where w is given
464	!*************************
465
466	! ww(ix,jy,1,n)=wwh(ix,jy,1)*pinmconv(1)
467	! ww(ix,jy,nz,n)=wwh(ix,jy,nwz)*pinmconv(nz)
468	! kmin=2
469	! do iz=2,nz
470	! do kz=kmin,nwz
471	! if ((height(iz).gt.wzlev(kz-1)).and. &
472	! (height(iz).le.wzlev(kz))) then
473	! dz1=height(iz)-wzlev(kz-1)
474	! dz2=wzlev(kz)-height(iz)
475	! dz=dz1+dz2
476	!! ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)dz2pinmconv(kz-1)
477	!! + +wwh(ix,jy,kz)dz1pinmconv(kz))/dz
478	! kmin=kz
479	! goto 40
480	! endif
481	! end do
482	!40 continue
483	! end do
484
485	if (method_w_terrain_correction .eq. 20) then
486	! apply w correction assuming that the WRF w is "absolute w";
487	! apply it here to wwh; set wwh=0 at iz=1
488	! do iz = 1, nz
489	! wtc_stat(1,iz) = wtc_stat(1,iz) + ww(ix,jy,iz,n)
490	! wtc_stat(2,iz) = wtc_stat(2,iz) + abs(ww(ix,jy,iz,n))
491	! end do
492
493	! if ((ix.eq.0) .and. (jy.eq.0)) write(,)
494	! & 'verttransform doing method_w_terrain_correction =',
495	! & method_w_terrain_correction
496	ix1 = max( ix-1, 0 )
497	jy1 = max( jy-1, 0 )
498	ixp = min( ix+1, nx-1 )
499	jyp = min( jy+1, ny-1 )
500	if (wind_option.eq.0) then
501	dzdx=(oro(ixp,jy) - oro(ix1,jy))/(dx(ixp-ix1)m_x(ix,jy,1))
502	dzdy=(oro(ix,jyp) - oro(ix,jy1))/(dy(jyp-jy1)m_y(ix,jy,1))
503	do kz = 1, nwz-1
504	wwh_svaa(kz) = wwh(ix,jy,kz)
505	wwh(ix,jy,kz) = wwh(ix,jy,kz)*pinmconv(kz) &
506	! wwh(ix,jy,kz) =
507	- (uuh(ix,jy,kz)dzdx + vvh(ix,jy,kz)dzdy) !this is correct. term of variation of geopot not necessary
508
509	if (kz .eq. 1) wwh(ix,jy,kz) = 0.0
510	aa=aa+1
511	end do
512	elseif (wind_option.ge.1) then
513	do kz = 2, nwz-1
514	wwh_svaa(kz) = wwh(ix,jy,kz)
515	! dzdx=(zzh(ixp,jy,kz,n) - zzh(ix1,jy,kz,n))
516	! + /(dx*(ixp-ix1))
517	! dzdy=(zzh(ix,jyp,kz,n) - zzh(ix,jy1,kz,n))
518	! + /(dy*(jyp-jy1))
519	! dzdx=(zzh(ixp,jy,kz,n) - zzh(ix1,jy,kz,n)-zzh(ixp,jy,1,n) &
520	! +zzh(ix1,jy,1,n)) &
521	! /(dx(ixp-ix1)m_x(ix,jy,1))
522	! dzdy=(zzh(ix,jyp,kz,n) - zzh(ix,jy1,kz,n)-zzh(ix,jyp,1,n) &
523	! +zzh(ix,jy1,1,n)) &
524	! /(dy(jyp-jy1)m_y(ix,jy,1))
525
526	dzdx=(zzh(ixp,jy,kz+add_sfc_level,n) - zzh(ix1,jy,kz+add_sfc_level,n) &
527	-zzh(ixp,jy,1,n)+zzh(ix1,jy,1,n))/(dx(ixp-ix1)m_x(ix,jy,1))
528	dzdy=(zzh(ix,jyp,kz+add_sfc_level,n) - zzh(ix,jy1,kz+add_sfc_level,n) &
529	-zzh(ix,jyp,1,n)+zzh(ix,jy1,1,n))/(dy(jyp-jy1)m_y(ix,jy,1))
530	u=0.5*(uuh(ix,jy,kz+add_sfc_level)+uuh(ix,jy,kz-1+add_sfc_level))
531	v=0.5*(vvh(ix,jy,kz+add_sfc_level)+vvh(ix,jy,kz-1+add_sfc_level))
532	wwh(ix,jy,kz) = wwh(ix,jy,kz)*pinmconv(kz) &
533	+ (udzdx + vdzdy) ! variation of geopot on sigma is necessary
534
535	! wwh(ix,jy,kz) = wwh(ix,jy,kz)*pinmconv(kz) &
536	! + (uuh(ix,jy,kz)dzdx + vvh(ix,jy,kz)dzdy) ! variation of geopot on sigma is necessary
537	! if (kz .eq. 1) wwh(ix,jy,kz) = 0.0
538	if (kz .eq. 1) wwh(ix,jy,kz) = wwh(ix,jy,kz)*pinmconv(kz)
539	! aa=aa+1
540	end do
541	endif
542	if (wind_option.eq.-1) then
543	! ww(ix,jy,1,n)=wwh(ix,jy,1)
544	ww(ix,jy,1,n)=0.
545	do iz=2,nz
546	ww(ix,jy,iz,n)=ww(ix,jy,iz-1,n)-(height(iz)-height(iz-1))* &
547	div(ix,jy,iz-1)
548	enddo
549	else
550
551	ww(ix,jy,1,n)=wwh(ix,jy,1)
552	ww(ix,jy,nz,n)=wwh(ix,jy,nwz)
553	kmin=2
554	do iz=2,nz
555	do kz=kmin,nwz
556	if ((height(iz).gt.wzlev(kz-1)).and. &
557	(height(iz).le.wzlev(kz))) then
558	dz1=height(iz)-wzlev(kz-1)
559	dz2=wzlev(kz)-height(iz)
560	dz=dz1+dz2
561	ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*dz2 &
562	+wwh(ix,jy,kz)*dz1) &
563	/dz
564	kmin=kz
565	goto 4000
566	endif
567	end do
568	4000 continue
569	end do
570	endif
571
572	! do kz = 1, nwz
573	! wwh(ix,jy,kz) = wwh_svaa(kz)
574	! end do
575
576	! do iz = 1, nz
577	! wtc_stat(3,iz) = wtc_stat(3,iz) + ww(ix,jy,iz,n)
578	! wtc_stat(4,iz) = wtc_stat(4,iz) + abs(ww(ix,jy,iz,n))
579	! end do
580	end if
581
582	! Compute density gradients at intermediate levels
583	!*************************************************
584
585	drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
586	(height(2)-height(1))
587	do kz=2,nz-1
588	drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
589	(height(kz+1)-height(kz-1))
590	end do
591
592	drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
593
594	end do
595	end do
596	!!!$OMP END DO
597	!!!$OMP END PARALLEL
598
599	!****************************************************************
600	! Compute slope of eta levels in windward direction and resulting
601	! vertical wind correction
602	!
603	! The ECMWF model uses a hybrid-pressure vertical coordinate, "eta"
604	! The "eta" coordinate transitions from terrain-following near
605	! the surface to constant pressure in the stratosphere.
606	! The vertical velocities in the ECMWF grib files are "eta_dot"
607	! FLEXPART uses a "height above ground" vertical coordinate
608	! which we will call "hag".
609	! The vertical velocity is uses (in ww array) is "hag_dot".
610	! Converting from eta_dot to hag_dot involves
611	! >> multiplying by pinmconv = [d(hag)/d(eta)]
612	! >> adding a term that accounts for the fact that
613	! "eta" varies on constant "hag" surfaces.
614	! This term is [ud(hag)/dx + vd(hag)/dy], with the
615	! partial derivatives taken with "eta" being constant
616	!
617	! The WRF model uses a similar (to ECMWF) vertical coordinate.
618	! HOWEVER, the vertical velocities in the WRF output files
619	! are the "true/absolute w" in m/s. (Is this true?)
620	! Converting from "absolute w" to hag_dot involves
621	! adding a term that accounts for the fact that
622	! "absolute z" varies on constant "hag" surfaces.
623	! This term is [- ud(oro)/dx - vd(oro)/dy]
624	!
625	! The FLEXPART code did not apply the terrain corrections
626	! at jy=0 & ny-1; ix=0 & nx-1; iz=1 & nz.
627	! FLEXPART_WRF applies the correction at all grid points
628	!****************************************************************
629
630
631	! If north pole is in the domain, calculate wind velocities in polar
632	! stereographic coordinates
633	!*******************************************************************
634
635	if (nglobal) then
636	write(,)
637	write(,) '* stopping in verttransform *'
638	write(,) ' the nglobal code section should not be active'
639	write(,)
640	stop
641	! do 74 jy=int(switchnorthg)-2,nymin1
642	! ylat=ylat0+real(jy)*dy
643	! do 74 ix=0,nxmin1
644	! xlon=xlon0+real(ix)*dx
645	! do 74 iz=1,nz
646	!74 call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n),
647	! + vv(ix,jy,iz,n),uupol(ix,jy,iz,n),
648	! + vvpol(ix,jy,iz,n))
649	!
650	!
651	! do 76 iz=1,nz
652	!
653	!* CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
654	! xlon=xlon0+real(nx/2-1)*dx
655	! xlonr=xlon*pi/180.
656	! ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+
657	! & vv(nx/2-1,nymin1,iz,n)**2)
658	! if(vv(nx/2-1,nymin1,iz,n).lt.0.) then
659	! ddpol=atan(uu(nx/2-1,nymin1,iz,n)/
660	! & vv(nx/2-1,nymin1,iz,n))-xlonr
661	! else
662	! ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/
663	! & vv(nx/2-1,nymin1,iz,n))-xlonr
664	! endif
665	! if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
666	! if(ddpol.gt.2.0pi) ddpol=ddpol-2.0pi
667	!
668	!* CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
669	! xlon=180.0
670	! xlonr=xlon*pi/180.
671	! ylat=90.0
672	! uuaux=-ffpol*sin(xlonr+ddpol)
673	! vvaux=-ffpol*cos(xlonr+ddpol)
674	! call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,
675	! + vvpolaux)
676	!
677	! jy=nymin1
678	! do 76 ix=0,nxmin1
679	! uupol(ix,jy,iz,n)=uupolaux
680	! vvpol(ix,jy,iz,n)=vvpolaux
681	!76 continue
682	!
683	!
684	!* Fix: Set W at pole to the zonally averaged W of the next equator-
685	!* ward parallel of latitude
686	!
687	! do 85 iz=1,nz
688	! wdummy=0.
689	! jy=ny-2
690	! do 80 ix=0,nxmin1
691	!80 wdummy=wdummy+ww(ix,jy,iz,n)
692	! wdummy=wdummy/real(nx)
693	! jy=nymin1
694	! do 85 ix=0,nxmin1
695	!85 ww(ix,jy,iz,n)=wdummy
696
697	endif
698
699
700	! If south pole is in the domain, calculate wind velocities in polar
701	! stereographic coordinates
702	!*******************************************************************
703
704	if (sglobal) then
705	write(,)
706	write(,) '* stopping in verttransform *'
707	write(,) ' the sglobal code section should not be active'
708	write(,)
709	stop
710	! do 77 jy=0,int(switchsouthg)+3
711	! ylat=ylat0+real(jy)*dy
712	! do 77 ix=0,nxmin1
713	! xlon=xlon0+real(ix)*dx
714	! do 77 iz=1,nz
715	!77 call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n),
716	! + vv(ix,jy,iz,n),uupol(ix,jy,iz,n),
717	! + vvpol(ix,jy,iz,n))
718	!
719	! do 79 iz=1,nz
720	!
721	!* CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
722	! xlon=xlon0+real(nx/2-1)*dx
723	! xlonr=xlon*pi/180.
724	! ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+
725	! & vv(nx/2-1,0,iz,n)**2)
726	! if(vv(nx/2-1,0,iz,n).lt.0.) then
727	! ddpol=atan(uu(nx/2-1,0,iz,n)/
728	! & vv(nx/2-1,0,iz,n))+xlonr
729	! else
730	! ddpol=pi+atan(uu(nx/2-1,0,iz,n)/
731	! & vv(nx/2-1,0,iz,n))+xlonr
732	! endif
733	! if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
734	! if(ddpol.gt.2.0pi) ddpol=ddpol-2.0pi
735	!
736	!* CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
737	! xlon=180.0
738	! xlonr=xlon*pi/180.
739	! ylat=-90.0
740	! uuaux=+ffpol*sin(xlonr-ddpol)
741	! vvaux=-ffpol*cos(xlonr-ddpol)
742	! call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,
743	! + vvpolaux)
744	!
745	! jy=0
746	! do 79 ix=0,nxmin1
747	! uupol(ix,jy,iz,n)=uupolaux
748	!79 vvpol(ix,jy,iz,n)=vvpolaux
749	!
750	!
751	!* Fix: Set W at pole to the zonally averaged W of the next equator-
752	!* ward parallel of latitude
753	!
754	! do 95 iz=1,nz
755	! wdummy=0.
756	! jy=1
757	! do 90 ix=0,nxmin1
758	!90 wdummy=wdummy+ww(ix,jy,iz,n)
759	! wdummy=wdummy/real(nx)
760	! jy=0
761	! do 95 ix=0,nxmin1
762	!95 ww(ix,jy,iz,n)=wdummy
763	endif
764
765	!write (,) 'initializing clouds, n:',n,nymin1,nxmin1,nz^M
766	! create a cloud and rainout/washout field, clouds occur where rh>80%^M
767	! total cloudheight is stored at level 0^M
768	do 100 jy=0,nymin1
769	do 100 ix=0,nxmin1
770	! rain_cloud_above=0
771	lsp=lsprec(ix,jy,1,n)
772	convp=convprec(ix,jy,1,n)
773	! cloudsh(ix,jy,n)=0
774
775	prec=lsp+convp
776	if (lsp.gt.convp) then ! prectype='lsp'
777	lconvectprec = .false.
778	else ! prectype='cp '
779	lconvectprec = .true.
780	endif
781	rhmin = 0.90 ! standard condition for presence of clouds
782
783	!CPS note that original by Sabine Eckhart was 80%
784	!CPS however, for T<-20 C we consider saturation over ice
785	!CPS so I think 90% should be enough
786
787
788	icloudbot(ix,jy,n)=icmv
789	icloudtop=icmv ! this is just a local variable
790	98 do kz=1,nz
791	pressure=rho(ix,jy,kz,n)r_airtt(ix,jy,kz,n)
792	rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
793	!cps if (prec.gt.0.01) print*,'relhum',prec,kz,rh,height(kz)
794	if (rh .gt. rhmin) then
795	if (icloudbot(ix,jy,n) .eq. icmv) then
796	icloudbot(ix,jy,n)=nint(height(kz))
797	endif
798	icloudtop=nint(height(kz)) ! use int to save memory
799	endif
800	enddo
801
802
803	!CPS try to get a cloud thicker than 50 m
804	!CPS if there is at least .01 mm/h - changed to 0.002 and put into
805	!CPS parameter precpmin
806	if ((icloudbot(ix,jy,n) .eq. icmv .or. &
807	icloudtop-icloudbot(ix,jy,n) .lt. 50) .and. &
808	prec .gt. precmin) then
809	rhmin = rhmin - 0.05
810	if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum.
811	endif
812	!CPS implement a rough fix for badly represented convection
813	!CPS is based on looking at a limited set of comparison data
814	if (lconvectprec .and. icloudtop .lt. 6000 .and. &
815	prec .gt. precmin) then
816	if (convp .lt. 0.1) then
817	icloudbot(ix,jy,n) = 500
818	icloudtop = 8000
819	else
820	icloudbot(ix,jy,n) = 0
821	icloudtop = 10000
822	endif
823	endif
824	if (icloudtop .ne. icmv) then
825	icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n)
826	else
827	icloudthck(ix,jy,n) = icmv
828	endif
829	!CPS get rid of too thin clouds
830	if (icloudthck(ix,jy,n) .lt. 50) then
831	icloudbot(ix,jy,n)=icmv
832	icloudthck(ix,jy,n)=icmv
833	endif
834
835	100 continue
836
837
838
839
840
841
842	! do kz_inv=1,nz-1
843	! kz=nz-kz_inv+1
844	! pressure=rho(ix,jy,kz,n)r_airtt(ix,jy,kz,n)
845	! rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
846	! clouds(ix,jy,kz,n)=0
847	! if (rh.gt.0.8) then ! in cloud
848	! if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
849	! rain_cloud_above=1
850	! cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
851	! height(kz)-height(kz-1)
852	! if (lsp.ge.convp) then
853	! clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
854	! else
855	! clouds(ix,jy,kz,n)=2 ! convp dominated rainout
856	! endif
857	! else ! no precipitation
858	! clouds(ix,jy,kz,n)=1 ! cloud
859	! endif
860	! else ! no cloud
861	! if (rain_cloud_above.eq.1) then ! scavenging
862	! if (lsp.ge.convp) then
863	! clouds(ix,jy,kz,n)=5 ! lsp dominated washout
864	! else
865	! clouds(ix,jy,kz,n)=4 ! convp dominated washout
866	! endif
867	! endif
868	! endif
869	! end do
870	! end do
871	! end do
872
873
874	end subroutine verttransform

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

Download in other formats:

Original Format