CTBTO FLEXPART WO4 (2015-10-15)
 All Classes Files Functions Variables
verttransform_nests.f90
Go to the documentation of this file.
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.378*ew(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.608*qvhn(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)+const*log(pold/pint)* &
115  (tv-tvold)/log(tv/tvold)
116  else
117  uvzlev(kz)=uvzlev(kz-1)+const*log(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)*dxconst*xresoln(l)/ &
267  cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
268  vi=vvn(ix,jy,iz,n,l)*dyconst*yresoln(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=(dzdx1*dz2+dzdx2*dz1)/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=(dzdy1*dz2+dzdy2*dz1)/dz
295 
296  wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l)+(dzdx*ui+dzdy*vi)
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_air*ttn(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
ew
real function ew(x)
Definition: ew.f90:22
verttransform_nests
subroutine verttransform_nests(n, uuhn, vvhn, wwhn, pvhn)
Definition: verttransform_nests.f90:22
par_mod
Definition: par_mod.F90:34
f_qvsat
real function f_qvsat(p, t)
Definition: qvsat.f90:32
com_mod
Definition: com_mod.f90:13
Generated on Thu Oct 15 2015 06:20:58 for CTBTO FLEXPART WO4 (2015-10-15) by   doxygen 1.8.6