source: flexpart.git/src/verttransform_nests.f90 @ 02095e3

devrelease-10univie
Last change on this file since 02095e3 was 72d3a5a, checked in by Espen Sollum ATMOS <eso@…>, 2 years ago

Bugfix in verttransform (call to function ew)

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