source: flexpart.git/src/verttransform_nests.f90

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file was 92fab65, checked in by Ignacio Pisso <ip@…>, 4 years ago

add SPDX-License-Identifier to all .f90 files

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