source: flexpart.git/src/verttransform_nests.f90 @ 7999df47

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 7999df47 was db712a8, checked in by Espen Sollum ATMOS <eso@…>, 8 years ago

Completed handling of nested wind fields with cloud water (for wet deposition).

  • Property mode set to 100644
File size: 22.7 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! ESO, 2016
55! -note that divide-by-zero occurs when nxmaxn,nymaxn etc. are larger than
56!  the actual field dimensions
57!*****************************************************************************
58!                                                                            *
59! Variables:                                                                 *
60! nxn,nyn,nuvz,nwz                field dimensions in x,y and z direction    *
61! uun                             wind components in x-direction [m/s]       *
62! vvn                             wind components in y-direction [m/s]       *
63! wwn                             wind components in z-direction [deltaeta/s]*
64! ttn                             temperature [K]                            *
65! pvn                             potential vorticity (pvu)                  *
66! psn                             surface pressure [Pa]                      *
67!                                                                            *
68!*****************************************************************************
69
70  use par_mod
71  use com_mod
72
73  implicit none
74
75  real,intent(in),dimension(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests) :: uuhn,vvhn,pvhn
76  real,intent(in),dimension(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests) :: wwhn
77
78  real,dimension(0:nxmaxn-1,0:nymaxn-1,nuvzmax) :: rhohn,uvzlev,wzlev
79  real,dimension(0:nxmaxn-1,0:nymaxn-1,nzmax) :: pinmconv
80  real,dimension(0:nxmaxn-1,0:nymaxn-1) :: tvold,pold,pint,tv
81  real,dimension(0:nymaxn-1) :: cosf
82
83  integer,dimension(0:nxmaxn-1,0:nymaxn-1) :: rain_cloud_above, idx
84
85  integer :: ix,jy,kz,iz,n,l,kmin,kl,klp,ix1,jy1,ixp,jyp,kz_inv
86  real :: f_qvsat,pressure,rh,lsp,convp,cloudh_min,prec
87
88  real :: ew,dz1,dz2,dz
89  real :: dzdx,dzdy
90  real :: dzdx1,dzdx2,dzdy1,dzdy2
91  real,parameter :: const=r_air/ga
92  real :: tot_cloud_h
93  integer :: nxm1, nym1
94
95!  real,parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagnostics
96
97! Loop over all nests
98!********************
99
100  do l=1,numbnests
101    nxm1=nxn(l)-1
102    nym1=nyn(l)-1
103
104! Loop over the whole grid
105!*************************
106
107    do jy=0,nyn(l)-1
108      do ix=0,nxn(l)-1
109        tvold(ix,jy)=tt2n(ix,jy,1,n,l)*(1.+0.378*ew*(td2n(ix,jy,1,n,l))/ &
110             psn(ix,jy,1,n,l))
111      end do
112    end do
113
114    pold(0:nxm1,0:nym1)=psn(0:nxm1,0:nym1,1,n,l)
115    uvzlev(0:nxm1,0:nym1,1)=0.
116    wzlev(0:nxm1,0:nym1,1)=0.
117    rhohn(0:nxm1,0:nym1,1)=pold(0:nxm1,0:nym1)/(r_air*tvold(0:nxm1,0:nym1))
118
119! Compute heights of eta levels
120!******************************
121
122    do kz=2,nuvz
123      pint(0:nxm1,0:nym1)=akz(kz)+bkz(kz)*psn(0:nxm1,0:nym1,1,n,l)
124      tv(0:nxm1,0:nym1)=tthn(0:nxm1,0:nym1,kz,n,l)*(1.+0.608*qvhn(0:nxm1,0:nym1,kz,n,l))
125      rhohn(0:nxm1,0:nym1,kz)=pint(0:nxm1,0:nym1)/(r_air*tv(0:nxm1,0:nym1))
126
127      where (abs(tv(0:nxm1,0:nym1)-tvold(0:nxm1,0:nym1)).gt.0.2)
128        uvzlev(0:nxm1,0:nym1,kz)=uvzlev(0:nxm1,0:nym1,kz-1)+const*&
129             &log(pold(0:nxm1,0:nym1)/pint(0:nxm1,0:nym1))* &
130             (tv(0:nxm1,0:nym1)-tvold(0:nxm1,0:nym1))/&
131             &log(tv(0:nxm1,0:nym1)/tvold(0:nxm1,0:nym1))
132      elsewhere
133        uvzlev(0:nxm1,0:nym1,kz)=uvzlev(0:nxm1,0:nym1,kz-1)+const*&
134             &log(pold(0:nxm1,0:nym1)/pint(0:nxm1,0:nym1))*tv(0:nxm1,0:nym1)
135      endwhere
136
137      tvold(0:nxm1,0:nym1)=tv(0:nxm1,0:nym1)
138      pold(0:nxm1,0:nym1)=pint(0:nxm1,0:nym1)
139
140    end do
141
142    do kz=2,nwz-1
143      wzlev(0:nxm1,0:nym1,kz)=(uvzlev(0:nxm1,0:nym1,kz+1)+uvzlev(0:nxm1,0:nym1,kz))/2.
144    end do
145    wzlev(0:nxm1,0:nym1,nwz)=wzlev(0:nxm1,0:nym1,nwz-1)+ &
146         uvzlev(0:nxm1,0:nym1,nuvz)-uvzlev(0:nxm1,0:nym1,nuvz-1)
147
148
149    pinmconv(0:nxm1,0:nym1,1)=(uvzlev(0:nxm1,0:nym1,2))/ &
150         ((aknew(2)+bknew(2)*psn(0:nxm1,0:nym1,1,n,l))- &
151         (aknew(1)+bknew(1)*psn(0:nxm1,0:nym1,1,n,l)))
152    do kz=2,nz-1
153      pinmconv(0:nxm1,0:nym1,kz)=(uvzlev(0:nxm1,0:nym1,kz+1)-uvzlev(0:nxm1,0:nym1,kz-1))/ &
154           ((aknew(kz+1)+bknew(kz+1)*psn(0:nxm1,0:nym1,1,n,l))- &
155           (aknew(kz-1)+bknew(kz-1)*psn(0:nxm1,0:nym1,1,n,l)))
156    end do
157    pinmconv(0:nxm1,0:nym1,nz)=(uvzlev(0:nxm1,0:nym1,nz)-uvzlev(0:nxm1,0:nym1,nz-1))/ &
158         ((aknew(nz)+bknew(nz)*psn(0:nxm1,0:nym1,1,n,l))- &
159         (aknew(nz-1)+bknew(nz-1)*psn(0:nxm1,0:nym1,1,n,l)))
160
161! Levels, where u,v,t and q are given
162!************************************
163
164    uun(0:nxm1,0:nym1,1,n,l)=uuhn(0:nxm1,0:nym1,1,l)
165    vvn(0:nxm1,0:nym1,1,n,l)=vvhn(0:nxm1,0:nym1,1,l)
166    ttn(0:nxm1,0:nym1,1,n,l)=tthn(0:nxm1,0:nym1,1,n,l)
167    qvn(0:nxm1,0:nym1,1,n,l)=qvhn(0:nxm1,0:nym1,1,n,l)
168    if (readclouds_nest(l)) then
169      clwcn(0:nxm1,0:nym1,1,n,l)=clwchn(0:nxm1,0:nym1,1,n,l)
170      if (.not.sumclouds_nest(l)) ciwcn(0:nxm1,0:nym1,1,n,l)=ciwchn(0:nxm1,0:nym1,1,n,l)
171    end if
172    pvn(0:nxm1,0:nym1,1,n,l)=pvhn(0:nxm1,0:nym1,1,l)
173    rhon(0:nxm1,0:nym1,1,n,l)=rhohn(0:nxm1,0:nym1,1)
174
175    uun(0:nxm1,0:nym1,nz,n,l)=uuhn(0:nxm1,0:nym1,nuvz,l)
176    vvn(0:nxm1,0:nym1,nz,n,l)=vvhn(0:nxm1,0:nym1,nuvz,l)
177    ttn(0:nxm1,0:nym1,nz,n,l)=tthn(0:nxm1,0:nym1,nuvz,n,l)
178    qvn(0:nxm1,0:nym1,nz,n,l)=qvhn(0:nxm1,0:nym1,nuvz,n,l)
179    if (readclouds_nest(l)) then
180      clwcn(0:nxm1,0:nym1,nz,n,l)=clwchn(0:nxm1,0:nym1,nuvz,n,l)
181      if (.not.sumclouds_nest(l)) ciwcn(0:nxm1,0:nym1,nz,n,l)=ciwchn(0:nxm1,0:nym1,nuvz,n,l)
182    end if
183    pvn(0:nxm1,0:nym1,nz,n,l)=pvhn(0:nxm1,0:nym1,nuvz,l)
184    rhon(0:nxm1,0:nym1,nz,n,l)=rhohn(0:nxm1,0:nym1,nuvz)
185
186
187    kmin=2
188    idx=kmin
189    do iz=2,nz-1
190      do jy=0,nyn(l)-1
191        do ix=0,nxn(l)-1
192          if(height(iz).gt.uvzlev(ix,jy,nuvz)) then
193            uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l)
194            vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l)
195            ttn(ix,jy,iz,n,l)=ttn(ix,jy,nz,n,l)
196            qvn(ix,jy,iz,n,l)=qvn(ix,jy,nz,n,l)
197!hg adding the cloud water
198            if (readclouds_nest(l)) then
199              clwcn(ix,jy,iz,n,l)=clwcn(ix,jy,nz,n,l)
200              if (.not.sumclouds_nest(l)) ciwcn(ix,jy,iz,n,l)=ciwcn(ix,jy,nz,n,l)
201            end if
202!hg
203            pvn(ix,jy,iz,n,l)=pvn(ix,jy,nz,n,l)
204            rhon(ix,jy,iz,n,l)=rhon(ix,jy,nz,n,l)
205          else
206            innuvz: do kz=idx(ix,jy),nuvz
207              if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. &
208                   (height(iz).le.uvzlev(ix,jy,kz))) then
209                idx(ix,jy)=kz
210                exit innuvz
211              endif
212            enddo innuvz
213          endif
214        enddo
215      enddo
216      do jy=0,nyn(l)-1
217        do ix=0,nxn(l)-1
218          if(height(iz).le.uvzlev(ix,jy,nuvz)) then
219            kz=idx(ix,jy)
220            dz1=height(iz)-uvzlev(ix,jy,kz-1)
221            dz2=uvzlev(ix,jy,kz)-height(iz)
222            dz=dz1+dz2
223            uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+uuhn(ix,jy,kz,l)*dz1)/dz
224            vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+vvhn(ix,jy,kz,l)*dz1)/dz
225            ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2 &
226                 +tthn(ix,jy,kz,n,l)*dz1)/dz
227            qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2 &
228                 +qvhn(ix,jy,kz,n,l)*dz1)/dz
229!hg adding the cloud water
230            if (readclouds_nest(l)) then
231              clwcn(ix,jy,iz,n,l)=(clwchn(ix,jy,kz-1,n,l)*dz2+clwchn(ix,jy,kz,n,l)*dz1)/dz
232              if (.not.sumclouds_nest(l)) &
233                   &ciwcn(ix,jy,iz,n,l)=(ciwchn(ix,jy,kz-1,n,l)*dz2+ciwchn(ix,jy,kz,n,l)*dz1)/dz
234            end if
235!hg
236            pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+pvhn(ix,jy,kz,l)*dz1)/dz
237            rhon(ix,jy,iz,n,l)=(rhohn(ix,jy,kz-1)*dz2+rhohn(ix,jy,kz)*dz1)/dz
238          endif
239        enddo
240      enddo
241    enddo
242
243
244
245!         do kz=kmin,nuvz
246!           if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
247!             uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l)
248!             vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l)
249!             ttn(ix,jy,iz,n,l)=ttn(ix,jy,nz,n,l)
250!             qvn(ix,jy,iz,n,l)=qvn(ix,jy,nz,n,l)
251!             pvn(ix,jy,iz,n,l)=pvn(ix,jy,nz,n,l)
252!             rhon(ix,jy,iz,n,l)=rhon(ix,jy,nz,n,l)
253!             goto 30
254!           endif
255!           if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
256!           (height(iz).le.uvwzlev(ix,jy,kz))) then
257!            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
258!            dz2=uvwzlev(ix,jy,kz)-height(iz)
259!            dz=dz1+dz2
260!            uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+ &
261!            uuhn(ix,jy,kz,l)*dz1)/dz
262!            vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+ &
263!            vvhn(ix,jy,kz,l)*dz1)/dz
264!            ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2+ &
265!            tthn(ix,jy,kz,n,l)*dz1)/dz
266!            qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2+ &
267!            qvhn(ix,jy,kz,n,l)*dz1)/dz
268!            pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+ &
269!            pvhn(ix,jy,kz,l)*dz1)/dz
270!            rhon(ix,jy,iz,n,l)=(rhohn(kz-1)*dz2+rhohn(kz)*dz1)/dz
271!            kmin=kz
272!            goto 30
273!           endif
274!         end do
275! 30      continue
276!       end do
277
278
279! Levels, where w is given
280!*************************
281
282    wwn(0:nxm1,0:nym1,1,n,l)=wwhn(0:nxm1,0:nym1,1,l)*pinmconv(0:nxm1,0:nym1,1)
283    wwn(0:nxm1,0:nym1,nz,n,l)=wwhn(0:nxm1,0:nym1,nwz,l)*pinmconv(0:nxm1,0:nym1,nz)
284    kmin=2
285    idx=kmin
286    do iz=2,nz
287      do jy=0,nym1
288        do ix=0,nxm1
289          inn:         do kz=idx(ix,jy),nwz
290            if(idx(ix,jy) .le. kz .and. height(iz).gt.wzlev(ix,jy,kz-1).and. &
291                 height(iz).le.wzlev(ix,jy,kz)) then
292              idx(ix,jy)=kz
293              exit inn
294            endif
295          enddo inn
296        enddo
297      enddo
298      do jy=0,nym1
299        do ix=0,nxm1
300          kz=idx(ix,jy)
301          dz1=height(iz)-wzlev(ix,jy,kz-1)
302          dz2=wzlev(ix,jy,kz)-height(iz)
303          dz=dz1+dz2
304          wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(ix,jy,kz-1)*dz2 &
305               +wwhn(ix,jy,kz,l)*pinmconv(ix,jy,kz)*dz1)/dz
306        enddo
307      enddo
308    enddo
309
310!       wwn(ix,jy,1,n,l)=wwhn(ix,jy,1,l)*pinmconv(1)
311!       wwn(ix,jy,nz,n,l)=wwhn(ix,jy,nwz,l)*pinmconv(nz)
312!       kmin=2
313!       do iz=2,nz
314!         do kz=kmin,nwz
315!           if ((height(iz).gt.wzlev(kz-1)).and. &
316!           (height(iz).le.wzlev(kz))) then
317!             dz1=height(iz)-wzlev(kz-1)
318!             dz2=wzlev(kz)-height(iz)
319!             dz=dz1+dz2
320!             wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(kz-1)*dz2 &
321!             +wwhn(ix,jy,kz,l)*pinmconv(kz)*dz1)/dz
322!             kmin=kz
323!             goto 40
324!           endif
325!         end do
326! 40      continue
327!       end do
328
329! Compute density gradients at intermediate levels
330!*************************************************
331
332    drhodzn(0:nxm1,0:nym1,1,n,l)=(rhon(0:nxm1,0:nym1,2,n,l)-rhon(0:nxm1,0:nym1,1,n,l))/ &
333         (height(2)-height(1))
334    do kz=2,nz-1
335      drhodzn(0:nxm1,0:nym1,kz,n,l)=(rhon(0:nxm1,0:nym1,kz+1,n,l)-rhon(0:nxm1,0:nym1,kz-1,n,l))/ &
336           (height(kz+1)-height(kz-1))
337    end do
338    drhodzn(0:nxm1,0:nym1,nz,n,l)=drhodzn(0:nxm1,0:nym1,nz-1,n,l)
339
340
341! drhodzn(ix,jy,1,n,l)=(rhon(ix,jy,2,n,l)-rhon(ix,jy,1,n,l))/ &
342! (height(2)-height(1))
343! do kz=2,nz-1
344!   drhodzn(ix,jy,kz,n,l)=(rhon(ix,jy,kz+1,n,l)- &
345!   rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
346! end do
347! drhodzn(ix,jy,nz,n,l)=drhodzn(ix,jy,nz-1,n,l)
348
349
350
351!   end do
352! end do
353
354
355!****************************************************************
356! Compute slope of eta levels in windward direction and resulting
357! vertical wind correction
358!****************************************************************
359
360    do jy=1,nyn(l)-2
361!    cosf=cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
362      cosf(jy)=1./cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
363    end do
364
365    kmin=2
366    idx=kmin
367    do iz=2,nz-1
368      do jy=1,nyn(l)-2
369        do ix=1,nxn(l)-2
370
371          inneta: do kz=idx(ix,jy),nz
372            if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. &
373                 (height(iz).le.uvzlev(ix,jy,kz))) then
374              idx(ix,jy)=kz
375              exit inneta
376            endif
377          enddo inneta
378        enddo
379      enddo
380
381      do jy=1,nyn(l)-2
382        do ix=1,nxn(l)-2
383          kz=idx(ix,jy)
384          dz1=height(iz)-uvzlev(ix,jy,kz-1)
385          dz2=uvzlev(ix,jy,kz)-height(iz)
386          dz=dz1+dz2
387          ix1=ix-1
388          jy1=jy-1
389          ixp=ix+1
390          jyp=jy+1
391
392          dzdx1=(uvzlev(ixp,jy,kz-1)-uvzlev(ix1,jy,kz-1))/2.
393          dzdx2=(uvzlev(ixp,jy,kz)-uvzlev(ix1,jy,kz))/2.
394          dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
395
396          dzdy1=(uvzlev(ix,jyp,kz-1)-uvzlev(ix,jy1,kz-1))/2.
397          dzdy2=(uvzlev(ix,jyp,kz)-uvzlev(ix,jy1,kz))/2.
398          dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
399
400          wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l)+(dzdx*uun(ix,jy,iz,n,l)*dxconst*xresoln(l)*cosf(jy)+&
401               &dzdy*vvn(ix,jy,iz,n,l)*dyconst*yresoln(l))
402
403        end do
404
405      end do
406    end do
407
408
409!   do jy=1,nyn(l)-2
410!     do ix=1,nxn(l)-2
411!       kmin=2
412!       do iz=2,nz-1
413
414!         ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/cosf(jy)
415!         vi=vvn(ix,jy,iz,n,l)*dyconst*yresoln(l)
416
417!         do kz=kmin,nz
418!           if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
419!                (height(iz).le.uvwzlev(ix,jy,kz))) then
420!             dz1=height(iz)-uvwzlev(ix,jy,kz-1)
421!             dz2=uvwzlev(ix,jy,kz)-height(iz)
422!             dz=dz1+dz2
423!             kl=kz-1
424!             klp=kz
425!             kmin=kz
426!             goto 47
427!           endif
428!         end do
429
430! 47      ix1=ix-1
431!         jy1=jy-1
432!         ixp=ix+1
433!         jyp=jy+1
434
435!         dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
436!         dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
437!         dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
438
439!         dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
440!         dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
441!         dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
442
443!         wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l)+(dzdx*ui+dzdy*vi)
444
445!       end do
446
447!     end do
448!   end do
449
450
451!write (*,*) 'initializing nested cloudsn, n:',n
452!   create a cloud and rainout/washout field, cloudsn occur where rh>80%
453
454
455!*********************************************************************************** 
456    if (readclouds_nest(l)) then !HG METHOD
457! The method is loops all grids vertically and constructs the 3D matrix for clouds
458! Cloud top and cloud bottom gid cells are assigned as well as the total column
459! cloud water. For precipitating grids, the type and whether it is in or below
460! cloud scavenging are assigned with numbers 2-5 (following the old metod).
461! Distinction is done for lsp and convp though they are treated the same in regards
462! to scavenging. Also clouds that are not precipitating are defined which may be
463! to include future cloud processing by non-precipitating-clouds.
464!***********************************************************************************
465      write(*,*) 'Nested ECMWF fields: using cloud water'
466      clwn(0:nxm1,0:nym1,:,n,l)=0.0
467!    icloud_stats(0:nxm1,0:nym1,:,n)=0.0
468      clw4n(0:nxm1,0:nym1,n,l)=0.0
469      cloudsn(0:nxm1,0:nym1,:,n,l)=0
470! If water/ice are read separately into clwc and ciwc, store sum in clwcn
471      if (.not.sumclouds_nest(l)) then
472        clwcn(0:nxm1,0:nym1,:,n,l) = clwcn(0:nxm1,0:nym1,:,n,l) + ciwcn(0:nxm1,0:nym1,:,n,l)
473      end if
474      do jy=0,nym1
475        do ix=0,nxm1
476          lsp=lsprecn(ix,jy,1,n,l)
477          convp=convprecn(ix,jy,1,n,l)
478          prec=lsp+convp
479          tot_cloud_h=0
480! Find clouds in the vertical
481          do kz=1, nz-1 !go from top to bottom
482            if (clwcn(ix,jy,kz,n,l).gt.0) then     
483! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3
484              clwn(ix,jy,kz,n,l)=(clwcn(ix,jy,kz,n,l)*rhon(ix,jy,kz,n,l))*(height(kz+1)-height(kz))
485              tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz))
486              clw4n(ix,jy,n,l) = clw4n(ix,jy,n,l)+clwn(ix,jy,kz,n,l)
487!            icloud_stats(ix,jy,4,n)= icloud_stats(ix,jy,4,n)+clw(ix,jy,kz,n)          ! Column cloud water [m3/m3]
488!           icloud_stats(ix,jy,3,n)= min(height(kz+1),height(kz))                     ! Cloud BOT height stats      [m]
489              cloudh_min=min(height(kz+1),height(kz))
490!ZHG 2015 extra for testing
491!            clh(ix,jy,kz,n)=height(kz+1)-height(kz)
492!            icloud_stats(ix,jy,1,n)=icloud_stats(ix,jy,1,n)+(height(kz+1)-height(kz)) ! Cloud total vertical extent [m]
493!            icloud_stats(ix,jy,2,n)= max(icloud_stats(ix,jy,2,n),height(kz))          ! Cloud TOP height            [m]
494!ZHG
495            endif
496          end do
497
498! If Precipitation. Define removal type in the vertical
499          if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
500
501            do kz=nz,1,-1 !go Bottom up!
502              if (clwn(ix,jy,kz,n,l).gt.0.0) then ! is in cloud
503                cloudshn(ix,jy,n,l)=cloudshn(ix,jy,n,l)+height(kz)-height(kz-1)
504                cloudsn(ix,jy,kz,n,l)=1                               ! is a cloud
505                if (lsp.ge.convp) then
506                  cloudsn(ix,jy,kz,n,l)=3                            ! lsp in-cloud
507                else
508                  cloudsn(ix,jy,kz,n,l)=2                             ! convp in-cloud
509                endif                                              ! convective or large scale
510              else if((clwn(ix,jy,kz,n,l).le.0.0) .and. (cloudh_min.ge.height(kz))) then ! is below cloud
511                if (lsp.ge.convp) then
512                  cloudsn(ix,jy,kz,n,l)=5                             ! lsp dominated washout
513                else
514                  cloudsn(ix,jy,kz,n,l)=4                             ! convp dominated washout
515                endif                                              ! convective or large scale
516              endif
517
518              if (height(kz).ge. 19000) then                        ! set a max height for removal
519                cloudsn(ix,jy,kz,n,l)=0
520              endif !clw>0
521            end do !nz
522          endif ! precipitation
523        end do
524      end do
525!********************************************************************
526    else ! old method:
527!********************************************************************
528      write(*,*) 'Nested fields: using cloud water from Parameterization'
529      do jy=0,nyn(l)-1
530        do ix=0,nxn(l)-1
531          rain_cloud_above(ix,jy)=0
532          lsp=lsprecn(ix,jy,1,n,l)
533          convp=convprecn(ix,jy,1,n,l)
534          cloudshn(ix,jy,n,l)=0
535          do kz_inv=1,nz-1
536            kz=nz-kz_inv+1
537            pressure=rhon(ix,jy,kz,n,l)*r_air*ttn(ix,jy,kz,n,l)
538            rh=qvn(ix,jy,kz,n,l)/f_qvsat(pressure,ttn(ix,jy,kz,n,l))
539            cloudsn(ix,jy,kz,n,l)=0
540            if (rh.gt.0.8) then ! in cloud
541              if ((lsp.gt.0.01).or.(convp.gt.0.01)) then
542                rain_cloud_above(ix,jy)=1
543                cloudshn(ix,jy,n,l)=cloudshn(ix,jy,n,l)+ &
544                     height(kz)-height(kz-1)
545                if (lsp.ge.convp) then
546                  cloudsn(ix,jy,kz,n,l)=3 ! lsp dominated rainout
547                else
548                  cloudsn(ix,jy,kz,n,l)=2 ! convp dominated rainout
549                endif
550              else ! no precipitation
551                cloudsn(ix,jy,kz,n,l)=1 ! cloud
552              endif
553            else ! no cloud
554              if (rain_cloud_above(ix,jy).eq.1) then ! scavenging
555                if (lsp.ge.convp) then
556                  cloudsn(ix,jy,kz,n,l)=5 ! lsp dominated washout
557                else
558                  cloudsn(ix,jy,kz,n,l)=4 ! convp dominated washout
559                endif
560              endif
561            endif
562          end do
563        end do
564      end do
565    end if
566
567  end do ! end loop over nests
568
569end subroutine verttransform_nests
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG