source: trunk/src/verttransform_nests.f90 @ 28

Last change on this file since 28 was 24, checked in by igpis, 10 years ago

version 9.2 beta. Changes from HH, AST, MC, NIK, IP. Changes in vert transform. New SPECIES input includes scavenging coefficients

File size: 12.4 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
22subroutine 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 :: 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,cosf
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      uvwzlev(ix,jy,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          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
115          (tv-tvold)/log(tv/tvold)
116        else
117          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,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)=(uvwzlev(ix,jy,kz+1)+uvwzlev(ix,jy,kz))/2.
127      end do
128      wzlev(nwz)=wzlev(nwz-1)+uvwzlev(ix,jy,nuvz)-uvwzlev(ix,jy,nuvz-1)
129
130  ! pinmconv=(h2-h1)/(p2-p1)
131
132      pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
133      ((aknew(2)+bknew(2)*psn(ix,jy,1,n,l))- &
134      (aknew(1)+bknew(1)*psn(ix,jy,1,n,l)))
135      do kz=2,nz-1
136        pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
137        ((aknew(kz+1)+bknew(kz+1)*psn(ix,jy,1,n,l))- &
138        (aknew(kz-1)+bknew(kz-1)*psn(ix,jy,1,n,l)))
139      end do
140      pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
141      ((aknew(nz)+bknew(nz)*psn(ix,jy,1,n,l))- &
142      (aknew(nz-1)+bknew(nz-1)*psn(ix,jy,1,n,l)))
143
144
145  ! Levels, where u,v,t and q are given
146  !************************************
147
148      uun(ix,jy,1,n,l)=uuhn(ix,jy,1,l)
149      vvn(ix,jy,1,n,l)=vvhn(ix,jy,1,l)
150      ttn(ix,jy,1,n,l)=tthn(ix,jy,1,n,l)
151      qvn(ix,jy,1,n,l)=qvhn(ix,jy,1,n,l)
152      pvn(ix,jy,1,n,l)=pvhn(ix,jy,1,l)
153      rhon(ix,jy,1,n,l)=rhoh(1)
154      uun(ix,jy,nz,n,l)=uuhn(ix,jy,nuvz,l)
155      vvn(ix,jy,nz,n,l)=vvhn(ix,jy,nuvz,l)
156      ttn(ix,jy,nz,n,l)=tthn(ix,jy,nuvz,n,l)
157      qvn(ix,jy,nz,n,l)=qvhn(ix,jy,nuvz,n,l)
158      pvn(ix,jy,nz,n,l)=pvhn(ix,jy,nuvz,l)
159      rhon(ix,jy,nz,n,l)=rhoh(nuvz)
160      kmin=2
161      do iz=2,nz-1
162        do kz=kmin,nuvz
163          if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
164            uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l)
165            vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l)
166            ttn(ix,jy,iz,n,l)=ttn(ix,jy,nz,n,l)
167            qvn(ix,jy,iz,n,l)=qvn(ix,jy,nz,n,l)
168            pvn(ix,jy,iz,n,l)=pvn(ix,jy,nz,n,l)
169            rhon(ix,jy,iz,n,l)=rhon(ix,jy,nz,n,l)
170            goto 30
171          endif
172          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
173          (height(iz).le.uvwzlev(ix,jy,kz))) then
174           dz1=height(iz)-uvwzlev(ix,jy,kz-1)
175           dz2=uvwzlev(ix,jy,kz)-height(iz)
176           dz=dz1+dz2
177           uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+ &
178           uuhn(ix,jy,kz,l)*dz1)/dz
179           vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+ &
180           vvhn(ix,jy,kz,l)*dz1)/dz
181           ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2+ &
182           tthn(ix,jy,kz,n,l)*dz1)/dz
183           qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2+ &
184           qvhn(ix,jy,kz,n,l)*dz1)/dz
185           pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+ &
186           pvhn(ix,jy,kz,l)*dz1)/dz
187           rhon(ix,jy,iz,n,l)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
188           kmin=kz
189           goto 30
190          endif
191        end do
19230      continue
193      end do
194
195
196  ! Levels, where w is given
197  !*************************
198
199      wwn(ix,jy,1,n,l)=wwhn(ix,jy,1,l)*pinmconv(1)
200      wwn(ix,jy,nz,n,l)=wwhn(ix,jy,nwz,l)*pinmconv(nz)
201      kmin=2
202      do iz=2,nz
203        do kz=kmin,nwz
204          if ((height(iz).gt.wzlev(kz-1)).and. &
205          (height(iz).le.wzlev(kz))) then
206            dz1=height(iz)-wzlev(kz-1)
207            dz2=wzlev(kz)-height(iz)
208            dz=dz1+dz2
209            wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(kz-1)*dz2 &
210            +wwhn(ix,jy,kz,l)*pinmconv(kz)*dz1)/dz
211            kmin=kz
212            goto 40
213          endif
214        end do
21540      continue
216      end do
217
218  ! Compute density gradients at intermediate levels
219  !*************************************************
220
221      drhodzn(ix,jy,1,n,l)=(rhon(ix,jy,2,n,l)-rhon(ix,jy,1,n,l))/ &
222      (height(2)-height(1))
223      do kz=2,nz-1
224        drhodzn(ix,jy,kz,n,l)=(rhon(ix,jy,kz+1,n,l)- &
225        rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
226      end do
227      drhodzn(ix,jy,nz,n,l)=drhodzn(ix,jy,nz-1,n,l)
228
229    end do
230  end do
231
232
233  !****************************************************************
234  ! Compute slope of eta levels in windward direction and resulting
235  ! vertical wind correction
236  !****************************************************************
237
238  do jy=1,nyn(l)-2
239    cosf=cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
240    do ix=1,nxn(l)-2
241
242      kmin=2
243      do iz=2,nz-1
244
245        ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/cosf
246        vi=vvn(ix,jy,iz,n,l)*dyconst*yresoln(l)
247
248        do kz=kmin,nz
249          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
250               (height(iz).le.uvwzlev(ix,jy,kz))) then
251            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
252            dz2=uvwzlev(ix,jy,kz)-height(iz)
253            dz=dz1+dz2
254            kl=kz-1
255            klp=kz
256            kmin=kz
257            goto 47
258          endif
259        end do
260
26147      ix1=ix-1
262        jy1=jy-1
263        ixp=ix+1
264        jyp=jy+1
265
266        dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
267        dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
268        dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
269
270        dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
271        dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
272        dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
273
274        wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l)+(dzdx*ui+dzdy*vi)
275
276      end do
277
278    end do
279  end do
280
281
282  !write (*,*) 'initializing nested cloudsn, n:',n
283  !   create a cloud and rainout/washout field, cloudsn occur where rh>80%
284  do jy=0,nyn(l)-1
285    do ix=0,nxn(l)-1
286      rain_cloud_above=0
287      lsp=lsprecn(ix,jy,1,n,l)
288      convp=convprecn(ix,jy,1,n,l)
289      cloudsnh(ix,jy,n,l)=0
290      do kz_inv=1,nz-1
291         kz=nz-kz_inv+1
292         pressure=rhon(ix,jy,kz,n,l)*r_air*ttn(ix,jy,kz,n,l)
293         rh=qvn(ix,jy,kz,n,l)/f_qvsat(pressure,ttn(ix,jy,kz,n,l))
294         cloudsn(ix,jy,kz,n,l)=0
295         if (rh.gt.0.8) then ! in cloud
296            if ((lsp.gt.0.01).or.(convp.gt.0.01)) then
297               rain_cloud_above=1
298               cloudsnh(ix,jy,n,l)=cloudsnh(ix,jy,n,l)+ &
299               height(kz)-height(kz-1)
300               if (lsp.ge.convp) then
301                  cloudsn(ix,jy,kz,n,l)=3 ! lsp dominated rainout
302               else
303                  cloudsn(ix,jy,kz,n,l)=2 ! convp dominated rainout
304               endif
305            else ! no precipitation
306                  cloudsn(ix,jy,kz,n,l)=1 ! cloud
307            endif
308         else ! no cloud
309            if (rain_cloud_above.eq.1) then ! scavenging
310               if (lsp.ge.convp) then
311                  cloudsn(ix,jy,kz,n,l)=5 ! lsp dominated washout
312               else
313                  cloudsn(ix,jy,kz,n,l)=4 ! convp dominated washout
314               endif
315            endif
316         endif
317      end do
318    end do
319  end do
320
321  end do
322
323end subroutine verttransform_nests
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG