verttransform_nests.f90 24.5 KB
Newer Older
Matthias Langer's avatar
 
Matthias Langer committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
!                                                                     *
! This file is part of FLEXPART.                                      *
!                                                                     *
! FLEXPART is free software: you can redistribute it and/or modify    *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or   *
! (at your option) any later version.                                 *
!                                                                     *
! FLEXPART is distributed in the hope that it will be useful,         *
! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
! GNU General Public License for more details.                        *
!                                                                     *
! You should have received a copy of the GNU General Public License   *
! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
!**********************************************************************

subroutine verttransform_nests(n,uuhn,vvhn,wwhn,pvhn)
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
!                            i   i    i    i   i
!*****************************************************************************
!                                                                            *
!     This subroutine transforms temperature, dew point temperature and      *
!     wind components from eta to meter coordinates.                         *
!     The vertical wind component is transformed from Pa/s to m/s using      *
!     the conversion factor pinmconv.                                        *
!     In addition, this routine calculates vertical density gradients        *
!     needed for the parameterization of the turbulent velocities.           *
!     It is similar to verttransform, but makes the transformations for      *
!     the nested grids.                                                      *
!                                                                            *
!     Author: A. Stohl, G. Wotawa                                            *
!                                                                            *
!     12 August 1996                                                         *
!     Update: 16 January 1998                                                *
!                                                                            *
40
41
42
!                                                                            *
!*****************************************************************************
!  CHANGES                                                                   *
43
44
45
46
47
48
49
50
51
52
!     Major update: 17 February 1999                                         *
!     by G. Wotawa                                                           *
!                                                                            *
!     - Vertical levels for u, v and w are put together                      *
!     - Slope correction for vertical velocity: Modification of calculation  *
!       procedure                                                            *
!                                                                            *
!*****************************************************************************
!  Changes, Bernd C. Krueger, Feb. 2001:       (marked "C-cv")
!   Variables tthn and qvhn (on eta coordinates) from common block
53
54
55
56
57
58
59
60
! Sabine Eckhardt, March 2007:
!   added the variable cloud for use with scavenging - descr. in com_mod
!                                                                            *
! Who? When?                                                                 *
!   Unified ECMWF and GFS builds
! Marian Harustak, 12.5.2017 
!     - Renamed from verttransform to verttransform_ecmwf
!                                                                            *
61
62
! ESO, 2016
! -note that divide-by-zero occurs when nxmaxn,nymaxn etc. are larger than 
63
64
65
66
67
68
69
70
71
72
73
74
75
!  the actual field dimensions /fixed, PS 2018/
!                                                                            *
! Don Morton, 2017-05-30: 
!   modification of a bug in ew. Don Morton (CTBTO project)                  *
!                                                                            *
!  undocumented modifications by NILU for v10                                *
!                                                                            *
!  Petra Seibert, 2018-06-13:                                                *
!   - insert proper boundaries for implied loops in array expressions        *s
!   - minor changes, most of them just cosmetics                             *
!   for details see changelog.txt                                            *
!                                                                            *
! ****************************************************************************
76
77
78
79
80
81
82
83
84
85
86
!                                                                            *
! Variables:                                                                 *
! nxn,nyn,nuvz,nwz                field dimensions in x,y and z direction    *
! uun                             wind components in x-direction [m/s]       *
! vvn                             wind components in y-direction [m/s]       *
! wwn                             wind components in z-direction [deltaeta/s]*
! ttn                             temperature [K]                            *
! pvn                             potential vorticity (pvu)                  *
! psn                             surface pressure [Pa]                      *
!                                                                            *
!*****************************************************************************
Matthias Langer's avatar
 
Matthias Langer committed
87
88
89
90
91
92

  use par_mod
  use com_mod

  implicit none

93
94
  real,intent(in),dimension(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests) :: &
    uuhn,vvhn,pvhn
95
96
97
98
  real,intent(in),dimension(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests) :: wwhn

  real,dimension(0:nxmaxn-1,0:nymaxn-1,nuvzmax) :: rhohn,uvzlev,wzlev
  real,dimension(0:nxmaxn-1,0:nymaxn-1,nzmax) :: pinmconv
99

100
  real,dimension(0:nymaxn-1) :: cosf
101
102
  real,dimension(0:nxmaxn-1,0:nymaxn-1) :: tvold,pold,pint,tv
!! automatic arrays introduced in v10 by ?? to achieve better loop order (PS)
103
104
105
106

  integer,dimension(0:nxmaxn-1,0:nymaxn-1) :: rain_cloud_above, idx

  integer :: ix,jy,kz,iz,n,l,kmin,kl,klp,ix1,jy1,ixp,jyp,kz_inv
107
  integer :: nxx, nyy ! max of nest where we are working
108
109
110
  real :: f_qvsat,pressure,rh,lsp,convp,cloudh_min,prec

  real :: ew,dz1,dz2,dz
Matthias Langer's avatar
 
Matthias Langer committed
111
112
113
  real :: dzdx,dzdy
  real :: dzdx1,dzdx2,dzdy1,dzdy2
  real,parameter :: const=r_air/ga
114
  real :: tot_cloud_h
Matthias Langer's avatar
 
Matthias Langer committed
115

116
!  real,parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagn.
Matthias Langer's avatar
 
Matthias Langer committed
117

118
119
! Loop over all nests
!********************
Matthias Langer's avatar
 
Matthias Langer committed
120

121
122
123
  l_loop: do l=1,numbnests
    nxx=nxn(l)-1 ! shorthand
    nyy=nyn(l)-1 ! shorthand
Matthias Langer's avatar
 
Matthias Langer committed
124

125
! Loop over the whole grid (partly implicitly
126
!*************************
127
! initialise the automatic arrays
Matthias Langer's avatar
 
Matthias Langer committed
128

129
130
131
132
133
134
  tvold(0:nxx,0:nyy)=tt2n(0:nxx,0:nyy,1,n,l) * &
    (1.+0.378*ew( td2n(0:nxx,0:nyy,1,n,l) ) / psn(0:nxx,0:nyy,1,n,l))
  pold(0:nxx,0:nyy)=psn(0:nxx,0:nyy,1,n,l)
  uvzlev(0:nxx,0:nyy,1)=0.
  wzlev(0:nxx,0:nyy,1)=0.
  rhohn(0:nxx,0:nyy,1)=pold(0:nxx,0:nyy)/(r_air*tvold(0:nxx,0:nyy))
Matthias Langer's avatar
 
Matthias Langer committed
135
136


137
138
! Compute heights of eta levels
!******************************
Matthias Langer's avatar
 
Matthias Langer committed
139

140
    do kz=2,nuvz
141
142
143
144
145
146
147
148
149
150
      pint(:,:)=akz(kz)+bkz(kz)*psn(0:nxx,0:nyy,1,n,l)
      tv(0:nxx,0:nyy)=tthn(0:nxx,0:nyy,kz,n,l)* &
        (1.+0.608*qvhn(0:nxx,0:nyy,kz,n,l))
      rhohn(0:nxx,0:nyy,kz)=pint(0:nxx,0:nyy)/(r_air*tv(0:nxx,0:nyy))

      where (abs(tv(0:nxx,0:nyy)-tvold(0:nxx,0:nyy)).gt.0.2) 
        uvzlev(0:nxx,0:nyy,kz)=uvzlev(0:nxx,0:nyy,kz-1) + &
          const*log( pold(0:nxx,0:nyy)/pint(0:nxx,0:nyy) ) * &
          ( tv(0:nxx,0:nyy) - tvold(0:nxx,0:nyy) ) / &
          log( tv(0:nxx,0:nyy)/tvold(0:nxx,0:nyy) )
151
      elsewhere
152
153
154
        uvzlev(0:nxx,0:nyy,kz)=uvzlev(0:nxx,0:nyy,kz-1) + &
           const*log( pold(0:nxx,0:nyy)/pint(0:nxx,0:nyy) ) *  &
           tv(0:nxx,0:nyy)
155
      endwhere
Matthias Langer's avatar
 
Matthias Langer committed
156

157
158
      tvold(0:nxx,0:nyy)=tv(0:nxx,0:nyy)
      pold(0:nxx,0:nyy)=pint(0:nxx,0:nyy)
Matthias Langer's avatar
 
Matthias Langer committed
159

160
    end do
Matthias Langer's avatar
 
Matthias Langer committed
161

162
    do kz=2,nwz-1
163
      wzlev(0:nxx,0:nyy,kz)=(uvzlev(0:nxx,0:nyy,kz+1)+uvzlev(0:nxx,0:nyy,kz))/2.
164
    end do
165
166
    wzlev(0:nxx,0:nyy,nwz)=wzlev(0:nxx,0:nyy,nwz-1)+ &
         uvzlev(0:nxx,0:nyy,nuvz)-uvzlev(0:nxx,0:nyy,nuvz-1)
Matthias Langer's avatar
 
Matthias Langer committed
167
168


169
170
171
    pinmconv(0:nxx,0:nyy,1)=(uvzlev(0:nxx,0:nyy,2))/ &
      ((aknew(2)+bknew(2)*psn(0:nxx,0:nyy,1,n,l))- &
       (aknew(1)+bknew(1)*psn(0:nxx,0:nyy,1,n,l)))
172
    do kz=2,nz-1
173
174
175
176
      pinmconv(0:nxx,0:nyy,kz)= &
        (uvzlev(0:nxx,0:nyy,kz+1)-uvzlev(0:nxx,0:nyy,kz-1))/ &
        ((aknew(kz+1)+bknew(kz+1)*psn(0:nxx,0:nyy,1,n,l))- &
         (aknew(kz-1)+bknew(kz-1)*psn(0:nxx,0:nyy,1,n,l)))
177
    end do
178
179
180
181
    pinmconv(0:nxx,0:nyy,nz)= & 
      (uvzlev(0:nxx,0:nyy,nz)-uvzlev(0:nxx,0:nyy,nz-1))/ &
      ((aknew(nz)+bknew(nz)*psn(0:nxx,0:nyy,1,n,l))- &
       (aknew(nz-1)+bknew(nz-1)*psn(0:nxx,0:nyy,1,n,l)))
182

183
! Levels where u,v,t and q are given
184
185
!************************************

186
187
188
189
    uun(0:nxx,0:nyy,1,n,l)=uuhn(0:nxx,0:nyy,1,l)
    vvn(0:nxx,0:nyy,1,n,l)=vvhn(0:nxx,0:nyy,1,l)
    ttn(0:nxx,0:nyy,1,n,l)=tthn(0:nxx,0:nyy,1,n,l)
    qvn(0:nxx,0:nyy,1,n,l)=qvhn(0:nxx,0:nyy,1,n,l)
190
    if (readclouds_nest(l)) then
191
192
193
      clwcn(0:nxx,0:nyy,1,n,l)=clwchn(0:nxx,0:nyy,1,n,l)
      if (.not. sumclouds_nest(l)) &
        ciwcn(0:nxx,0:nyy,1,n,l)=ciwchn(0:nxx,0:nyy,1,n,l)
194
    end if
195
196
    pvn(0:nxx,0:nyy,1,n,l)=pvhn(0:nxx,0:nyy,1,l)
    rhon(0:nxx,0:nyy,1,n,l)=rhohn(0:nxx,0:nyy,1)
197

198
199
200
201
    uun(0:nxx,0:nyy,nz,n,l)=uuhn(0:nxx,0:nyy,nuvz,l)
    vvn(0:nxx,0:nyy,nz,n,l)=vvhn(0:nxx,0:nyy,nuvz,l)
    ttn(0:nxx,0:nyy,nz,n,l)=tthn(0:nxx,0:nyy,nuvz,n,l)
    qvn(0:nxx,0:nyy,nz,n,l)=qvhn(0:nxx,0:nyy,nuvz,n,l)
202
    if (readclouds_nest(l)) then
203
204
205
      clwcn(0:nxx,0:nyy,nz,n,l)=clwchn(0:nxx,0:nyy,nuvz,n,l)
      if (.not. sumclouds_nest(l)) &
        ciwcn(0:nxx,0:nyy,nz,n,l)=ciwchn(0:nxx,0:nyy,nuvz,n,l)
206
    end if
207
208
    pvn(0:nxx,0:nyy,nz,n,l)=pvhn(0:nxx,0:nyy,nuvz,l)
    rhon(0:nxx,0:nyy,nz,n,l)=rhohn(0:nxx,0:nyy,nuvz)
209
210
211

    kmin=2
    idx=kmin
212
213
214
215
216
217
    iz_loop: do iz=2,nz-1

      do jy=0,nyy
        do ix=0,nxx
          if( height(iz).gt.uvzlev(ix,jy,nuvz)) then

Matthias Langer's avatar
 
Matthias Langer committed
218
219
220
221
            uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l)
            vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l)
            ttn(ix,jy,iz,n,l)=ttn(ix,jy,nz,n,l)
            qvn(ix,jy,iz,n,l)=qvn(ix,jy,nz,n,l)
222
223
224
!hg adding the cloud water
            if (readclouds_nest(l)) then
              clwcn(ix,jy,iz,n,l)=clwcn(ix,jy,nz,n,l)
225
226
              if (.not. sumclouds_nest(l))  &
                  ciwcn(ix,jy,iz,n,l)=ciwcn(ix,jy,nz,n,l)
227
228
            end if
!hg
Matthias Langer's avatar
 
Matthias Langer committed
229
230
            pvn(ix,jy,iz,n,l)=pvn(ix,jy,nz,n,l)
            rhon(ix,jy,iz,n,l)=rhon(ix,jy,nz,n,l)
231

232
          else
233
234
235
236

            kz_loop: do kz=idx(ix,jy),nuvz
              if (idx(ix,jy) .le. kz .and. height(iz).gt.uvzlev(ix,jy,kz-1) &
                  .and. height(iz).le.uvzlev(ix,jy,kz)) then
237
                idx(ix,jy)=kz
238
                exit kz_loop
239
              endif
240
241
            enddo kz_loop
            
Matthias Langer's avatar
 
Matthias Langer committed
242
          endif
243
244
        enddo
      enddo
245
246
247

      do jy=0,nyy
        do ix=0,nxx
248
249
250
251
          if(height(iz).le.uvzlev(ix,jy,nuvz)) then
            kz=idx(ix,jy)
            dz1=height(iz)-uvzlev(ix,jy,kz-1)
            dz2=uvzlev(ix,jy,kz)-height(iz)
252
            dz=dz1+dz2
253
254
255
256
257
258
259
260
            uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+uuhn(ix,jy,kz,l)*dz1)/dz
            vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+vvhn(ix,jy,kz,l)*dz1)/dz
            ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2 &
                 +tthn(ix,jy,kz,n,l)*dz1)/dz
            qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2 &
                 +qvhn(ix,jy,kz,n,l)*dz1)/dz
!hg adding the cloud water
            if (readclouds_nest(l)) then
261
262
263
264
              clwcn(ix,jy,iz,n,l)=(clwchn(ix,jy,kz-1,n,l)*dz2 + &
                clwchn(ix,jy,kz,n,l)*dz1)/dz
              if (.not. sumclouds_nest(l)) ciwcn(ix,jy,iz,n,l)= &
                  (ciwchn(ix,jy,kz-1,n,l)*dz2+ciwchn(ix,jy,kz,n,l)*dz1)/dz
265
266
267
268
            end if
!hg
            pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+pvhn(ix,jy,kz,l)*dz1)/dz
            rhon(ix,jy,iz,n,l)=(rhohn(ix,jy,kz-1)*dz2+rhohn(ix,jy,kz)*dz1)/dz
269
            
Matthias Langer's avatar
 
Matthias Langer committed
270
          endif
271
272
        enddo
      enddo
273
274

    enddo iz_loop
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314



!         do kz=kmin,nuvz
!           if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
!             uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l)
!             vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l)
!             ttn(ix,jy,iz,n,l)=ttn(ix,jy,nz,n,l)
!             qvn(ix,jy,iz,n,l)=qvn(ix,jy,nz,n,l)
!             pvn(ix,jy,iz,n,l)=pvn(ix,jy,nz,n,l)
!             rhon(ix,jy,iz,n,l)=rhon(ix,jy,nz,n,l)
!             goto 30
!           endif
!           if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
!           (height(iz).le.uvwzlev(ix,jy,kz))) then
!            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
!            dz2=uvwzlev(ix,jy,kz)-height(iz)
!            dz=dz1+dz2
!            uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+ &
!            uuhn(ix,jy,kz,l)*dz1)/dz
!            vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+ &
!            vvhn(ix,jy,kz,l)*dz1)/dz
!            ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2+ &
!            tthn(ix,jy,kz,n,l)*dz1)/dz
!            qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2+ &
!            qvhn(ix,jy,kz,n,l)*dz1)/dz
!            pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+ &
!            pvhn(ix,jy,kz,l)*dz1)/dz
!            rhon(ix,jy,iz,n,l)=(rhohn(kz-1)*dz2+rhohn(kz)*dz1)/dz
!            kmin=kz
!            goto 30
!           endif
!         end do
! 30      continue
!       end do


! Levels, where w is given
!*************************

315
316
    wwn(0:nxx,0:nyy,1,n,l)=wwhn(0:nxx,0:nyy,1,l)*pinmconv(0:nxx,0:nyy,1)
    wwn(0:nxx,0:nyy,nz,n,l)=wwhn(0:nxx,0:nyy,nwz,l)*pinmconv(0:nxx,0:nyy,nz)
317
318
    kmin=2
    idx=kmin
319
320
321
322
323
324
    iz_loop2: do iz=2,nz
      do jy=0,nyy
        do ix=0,nxx
          kz_loop2:   do kz=idx(ix,jy),nwz
            if (idx(ix,jy) .le. kz .and. height(iz).gt.wzlev(ix,jy,kz-1) &
                .and. height(iz).le.wzlev(ix,jy,kz)) then
325
              idx(ix,jy)=kz
326
              exit kz_loop2
327
            endif
328
          enddo kz_loop2
329
330
        enddo
      enddo
331
332
      do jy=0,nyy
        do ix=0,nxx
333
334
335
336
337
338
339
340
          kz=idx(ix,jy)
          dz1=height(iz)-wzlev(ix,jy,kz-1)
          dz2=wzlev(ix,jy,kz)-height(iz)
          dz=dz1+dz2
          wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(ix,jy,kz-1)*dz2 &
               +wwhn(ix,jy,kz,l)*pinmconv(ix,jy,kz)*dz1)/dz
        enddo
      enddo
341
    enddo iz_loop2
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364

!       wwn(ix,jy,1,n,l)=wwhn(ix,jy,1,l)*pinmconv(1)
!       wwn(ix,jy,nz,n,l)=wwhn(ix,jy,nwz,l)*pinmconv(nz)
!       kmin=2
!       do iz=2,nz
!         do kz=kmin,nwz
!           if ((height(iz).gt.wzlev(kz-1)).and. &
!           (height(iz).le.wzlev(kz))) then
!             dz1=height(iz)-wzlev(kz-1)
!             dz2=wzlev(kz)-height(iz)
!             dz=dz1+dz2
!             wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(kz-1)*dz2 &
!             +wwhn(ix,jy,kz,l)*pinmconv(kz)*dz1)/dz
!             kmin=kz
!             goto 40
!           endif
!         end do
! 40      continue
!       end do

! Compute density gradients at intermediate levels
!*************************************************

365
366
    drhodzn(0:nxx,0:nyy,1,n,l)= &
      ( rhon(0:nxx,0:nyy,2,n,l)-rhon(0:nxx,0:nyy,1,n,l) )/(height(2)-height(1))
367
    do kz=2,nz-1
368
369
370
      drhodzn(0:nxx,0:nyy,kz,n,l)= &
        ( rhon(0:nxx,0:nyy,kz+1,n,l)-rhon(0:nxx,0:nyy,kz-1,n,l) ) / &
        ( height(kz+1)-height(kz-1) )
Matthias Langer's avatar
 
Matthias Langer committed
371
    end do
372
    drhodzn(0:nxx,0:nyy,nz,n,l)=drhodzn(0:nxx,0:nyy,nz-1,n,l)
Matthias Langer's avatar
 
Matthias Langer committed
373
374


375
376
377
378
379
380
381
! drhodzn(ix,jy,1,n,l)=(rhon(ix,jy,2,n,l)-rhon(ix,jy,1,n,l))/ &
! (height(2)-height(1))
! do kz=2,nz-1
!   drhodzn(ix,jy,kz,n,l)=(rhon(ix,jy,kz+1,n,l)- &
!   rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
! end do
! drhodzn(ix,jy,nz,n,l)=drhodzn(ix,jy,nz-1,n,l)
Matthias Langer's avatar
 
Matthias Langer committed
382
383
384



385
386
!   end do
! end do
Matthias Langer's avatar
 
Matthias Langer committed
387
388


389
390
391
392
!****************************************************************
! Compute slope of eta levels in windward direction and resulting
! vertical wind correction
!****************************************************************
Matthias Langer's avatar
 
Matthias Langer committed
393

394
    do jy=1,nyy-1
395
!    cosf=cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
396
      cosf(jy) = 1. / cos( ( real(jy)*dyn(l) + ylat0n(l) ) * pi180 )
397
    end do
Matthias Langer's avatar
 
Matthias Langer committed
398

399
400
    kmin=2
    idx=kmin
401
402
403
    iz_loop3: do iz=2,nz-1
      do jy=1,nyy-1
        do ix=1,nxx-1
404

405
406
407
          kz_loop3: do kz=idx(ix,jy),nz
            if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)) &
                .and. (height(iz).le.uvzlev(ix,jy,kz))) then
408
              idx(ix,jy)=kz
409
              exit kz_loop3
410
            endif
411
          enddo kz_loop3
412
413
414
        enddo
      enddo

415
416
      do jy=1,nyy-1
        do ix=1,nxx-1
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
          kz=idx(ix,jy)
          dz1=height(iz)-uvzlev(ix,jy,kz-1)
          dz2=uvzlev(ix,jy,kz)-height(iz)
          dz=dz1+dz2
          ix1=ix-1
          jy1=jy-1
          ixp=ix+1
          jyp=jy+1

          dzdx1=(uvzlev(ixp,jy,kz-1)-uvzlev(ix1,jy,kz-1))/2.
          dzdx2=(uvzlev(ixp,jy,kz)-uvzlev(ix1,jy,kz))/2.
          dzdx=(dzdx1*dz2+dzdx2*dz1)/dz

          dzdy1=(uvzlev(ix,jyp,kz-1)-uvzlev(ix,jy1,kz-1))/2.
          dzdy2=(uvzlev(ix,jyp,kz)-uvzlev(ix,jy1,kz))/2.
          dzdy=(dzdy1*dz2+dzdy2*dz1)/dz

434
435
436
          wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l) + &
            ( dzdx*uun(ix,jy,iz,n,l)*dxconst*xresoln(l)*cosf(jy) + &
              dzdy*vvn(ix,jy,iz,n,l)*dyconst*yresoln(l) )
Matthias Langer's avatar
 
Matthias Langer committed
437

438
        end do
Matthias Langer's avatar
 
Matthias Langer committed
439
440

      end do
441
    end do iz_loop3
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489


!   do jy=1,nyn(l)-2
!     do ix=1,nxn(l)-2
!       kmin=2
!       do iz=2,nz-1

!         ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/cosf(jy)
!         vi=vvn(ix,jy,iz,n,l)*dyconst*yresoln(l)

!         do kz=kmin,nz
!           if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
!                (height(iz).le.uvwzlev(ix,jy,kz))) then
!             dz1=height(iz)-uvwzlev(ix,jy,kz-1)
!             dz2=uvwzlev(ix,jy,kz)-height(iz)
!             dz=dz1+dz2
!             kl=kz-1
!             klp=kz
!             kmin=kz
!             goto 47
!           endif
!         end do

! 47      ix1=ix-1
!         jy1=jy-1
!         ixp=ix+1
!         jyp=jy+1

!         dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
!         dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
!         dzdx=(dzdx1*dz2+dzdx2*dz1)/dz

!         dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
!         dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
!         dzdy=(dzdy1*dz2+dzdy2*dz1)/dz

!         wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l)+(dzdx*ui+dzdy*vi)

!       end do

!     end do
!   end do


!write (*,*) 'initializing nested cloudsn, n:',n
!   create a cloud and rainout/washout field, cloudsn occur where rh>80%


490
491
492
493
494
495
!****************************************************************************** 
    if (readclouds_nest(l)) then !HG METHOD 
    
! Loops over all grid cells vertically and construct the 3D matrix for clouds
! Cloud top and cloud bottom grid cells are assigned as well as the total column
! cloud water. For precipitating columns, the type and whether it is in or below
496
! cloud scavenging are assigned with numbers 2-5 (following the old metod).
497
498
499
500
501
502
503
504
505
506
! A distinction is made between lsp and convp though they are treated the equally 
! with regard to scavenging. Also, clouds that are not precipitating are defined which 
! may be used in the future for cloud processing by non-precipitating-clouds. 
!*******************************************************************************

!PS      write(*,*) 'Nested ECMWF fields: using cloud water'
      clwn(0:nxx,0:nyy,:,n,l)=0.
!    icloud_stats(0:nxx,0:nyy,:,n)=0.
      ctwcn(0:nxx,0:nyy,n,l)=0.
      cloudsn(0:nxx,0:nyy,:,n,l)=0
507
! If water/ice are read separately into clwc and ciwc, store sum in clwcn
508
509
510
      if (.not. sumclouds_nest(l)) then 
        clwcn(0:nxx,0:nyy,:,n,l) = clwcn(0:nxx,0:nyy,:,n,l) + &
                                   ciwcn(0:nxx,0:nyy,:,n,l)
511
      end if
512
513
      do jy=0,nyy
        do ix=0,nxx
514
515
516
517
518
          lsp=lsprecn(ix,jy,1,n,l)
          convp=convprecn(ix,jy,1,n,l)
          prec=lsp+convp
          tot_cloud_h=0
! Find clouds in the vertical
519
!! Note PS: bad loop order. 
520
521
          do kz=1, nz-1 !go from top to bottom
            if (clwcn(ix,jy,kz,n,l).gt.0) then      
522
523
524
! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3
              clwn(ix,jy,kz,n,l)=(clwcn(ix,jy,kz,n,l)*rhon(ix,jy,kz,n,l)) * &
                (height(kz+1)-height(kz)) 
525
              tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz)) 
526
              ctwcn(ix,jy,n,l) = ctwcn(ix,jy,n,l)+clwn(ix,jy,kz,n,l)
527
528
529
530
531
532
533
534
535
536
537
538
!            icloud_stats(ix,jy,4,n)= icloud_stats(ix,jy,4,n)+clw(ix,jy,kz,n)          ! Column cloud water [m3/m3]
!           icloud_stats(ix,jy,3,n)= min(height(kz+1),height(kz))                     ! Cloud BOT height stats      [m]
              cloudh_min=min(height(kz+1),height(kz))
!ZHG 2015 extra for testing
!            clh(ix,jy,kz,n)=height(kz+1)-height(kz)
!            icloud_stats(ix,jy,1,n)=icloud_stats(ix,jy,1,n)+(height(kz+1)-height(kz)) ! Cloud total vertical extent [m]
!            icloud_stats(ix,jy,2,n)= max(icloud_stats(ix,jy,2,n),height(kz))          ! Cloud TOP height            [m]
!ZHG
            endif
          end do

! If Precipitation. Define removal type in the vertical
539
540
        if (lsp.gt.0.01 .or. convp.gt.0.01) then ! cloud and precipitation
!! Note PS: such hardcoded limits would better be parameters defined in par_mod
541
542

            do kz=nz,1,-1 !go Bottom up!
543
544
545
546
!! Note PS: bad loop order
             if (clwn(ix,jy,kz,n,l) .gt. 0.) then ! is in cloud
               cloudshn(ix,jy,n,l)=cloudshn(ix,jy,n,l)+height(kz)-height(kz-1) 
               cloudsn(ix,jy,kz,n,l)=1                             ! is a cloud
547
                if (lsp.ge.convp) then
548
                  cloudsn(ix,jy,kz,n,l)=3                        ! lsp in-cloud
549
                else
550
551
552
553
554
555
                  cloudsn(ix,jy,kz,n,l)=2                      ! convp in-cloud
                endif                               ! convective or large scale
              elseif ( clwn(ix,jy,kz,n,l) .le. 0. .and. &
                  cloudh_min .ge. height(kz) ) then ! is below cloud
                if (lsp .ge. convp) then
                  cloudsn(ix,jy,kz,n,l)=5              ! lsp dominated washout
556
                else
557
558
                  cloudsn(ix,jy,kz,n,l)=4            ! convp dominated washout
                endif                             ! convective or large scale 
559
560
              endif

561
              if (height(kz).ge. 19000) then ! set a max height for removal
562
                cloudsn(ix,jy,kz,n,l)=0
563
564
              endif ! clw>0
            end do ! kz
565
566
567
          endif ! precipitation
        end do
      end do
568
      
569
570
571
!********************************************************************
    else ! old method:
!********************************************************************
572
573
574
575
 
!PS      write(*,*) 'Nested fields: using cloud water from Parameterization'
      do jy=0,nyy
        do ix=0,nxx
576
577
578
579
580
          rain_cloud_above(ix,jy)=0
          lsp=lsprecn(ix,jy,1,n,l)
          convp=convprecn(ix,jy,1,n,l)
          cloudshn(ix,jy,n,l)=0
          do kz_inv=1,nz-1
581
!! Note PS: bad loop order. 
582
583
584
585
            kz=nz-kz_inv+1
            pressure=rhon(ix,jy,kz,n,l)*r_air*ttn(ix,jy,kz,n,l)
            rh=qvn(ix,jy,kz,n,l)/f_qvsat(pressure,ttn(ix,jy,kz,n,l))
            cloudsn(ix,jy,kz,n,l)=0
586
587
588
589
590
591

            if (rh .gt. 0.8) then ! in cloud
!! Note PS: such hardcoded limits would better be parameters defined in par_mod

              if (lsp.gt.0.01 .or. convp.gt.0.01) then ! cloud and precipitation
!! Note PS: such hardcoded limits would better be parameters defined in par_mod
592
593
594
595
                rain_cloud_above(ix,jy)=1
                cloudshn(ix,jy,n,l)=cloudshn(ix,jy,n,l)+ &
                     height(kz)-height(kz-1)
                if (lsp.ge.convp) then
Matthias Langer's avatar
 
Matthias Langer committed
596
                  cloudsn(ix,jy,kz,n,l)=3 ! lsp dominated rainout
597
                else
Matthias Langer's avatar
 
Matthias Langer committed
598
                  cloudsn(ix,jy,kz,n,l)=2 ! convp dominated rainout
599
600
601
602
                endif
              else ! no precipitation
                cloudsn(ix,jy,kz,n,l)=1 ! cloud
              endif
603

604
            else ! no cloud
605

606
607
              if (rain_cloud_above(ix,jy).eq.1) then ! scavenging
                if (lsp.ge.convp) then
Matthias Langer's avatar
 
Matthias Langer committed
608
                  cloudsn(ix,jy,kz,n,l)=5 ! lsp dominated washout
609
                else
Matthias Langer's avatar
 
Matthias Langer committed
610
                  cloudsn(ix,jy,kz,n,l)=4 ! convp dominated washout
611
612
                endif
              endif
Matthias Langer's avatar
 
Matthias Langer committed
613

614
615
616
617
618
619
620
621
622
            endif
          end do ! kz
          
        end do ! ix
      end do ! jy
      
    end if! old method:

  end do l_loop! end loop over nests
Matthias Langer's avatar
 
Matthias Langer committed
623
624

end subroutine verttransform_nests