source: branches/FP_AI/src/verttransform_nests.f90 @ 23

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

start tracking test environment directory FP_AI

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