verttransform_ecmwf.f90 32.6 KB
Newer Older
1
! **********************************************************************
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
! 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/>.   *
20
! **********************************************************************
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
!                         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.           *
!                                                                            *
!     Author: A. Stohl, G. Wotawa                                            *
!                                                                            *
!     12 August 1996                                                         *
!     Update: 16 January 1998                                                *
!                                                                            *
38
39
40
!                                                                            *
!*****************************************************************************
!  CHANGES                                                                   *
41
42
43
44
45
46
47
!     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                                                            *
!                                                                            *
48
!  Bernd C. Krueger, Feb. 2001:
49
!   Variables tth and qvh (on eta coordinates) from common block
50
51
52
53
54
55
!                                                                            *
! Sabine Eckhardt, March 2007:
!   added the variable cloud for use with scavenging - descr. in com_mod
!                                                                            *
! Who? When?                                                                 *
!   Unified ECMWF and GFS builds
56
57
! Marian Harustak, 12.5.2017 
!     - Renamed from verttransform to verttransform_ecmwf
58
59
60
61
62
63
64
65
66
67
68
69
70
71
!                                                                            *
! 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:                                                *
!   - fix bug of ticket:140 (find reference position for vertical grid)      *
!   - put back SAVE attribute for INIT, just to be safe                      *
!   - minor changes, most of them just cosmetics                             *
!   for details see changelog.txt                                            *
!                                                                            *
! ****************************************************************************

72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
!*****************************************************************************
!                                                                            *
! Variables:                                                                 *
! nx,ny,nz                        field dimensions in x,y and z direction    *
! clouds(0:nxmax,0:nymax,0:nzmax,numwfmem) cloud field for wet deposition    *
! uu(0:nxmax,0:nymax,nzmax,numwfmem)     wind components in x-direction [m/s]*
! vv(0:nxmax,0:nymax,nzmax,numwfmem)     wind components in y-direction [m/s]*
! ww(0:nxmax,0:nymax,nzmax,numwfmem)     wind components in z-direction      *
!                                          [deltaeta/s]                      *
! tt(0:nxmax,0:nymax,nzmax,numwfmem)     temperature [K]                     *
! pv(0:nxmax,0:nymax,nzmax,numwfmem)     potential voriticity (pvu)          *
! ps(0:nxmax,0:nymax,numwfmem)           surface pressure [Pa]               *
!                                                                            *
!*****************************************************************************

  use par_mod
  use com_mod
  use cmapf_mod, only: cc2gll

  implicit none

  real,intent(in),dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: uuh,vvh,pvh
  real,intent(in),dimension(0:nxmax-1,0:nymax-1,nwzmax) :: wwh

  real,dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: rhoh,uvzlev,wzlev
  real,dimension(0:nxmax-1,0:nymax-1,nzmax) :: pinmconv
98
99
100
101
102
103
104
!>   for reference profile (PS)  
  real ::  tvoldref, poldref, pintref, psmean, psstd
  integer :: ixyref(2)
  integer, parameter :: ioldref = 1 ! PS 2018-06-13: set to 0 if you 
!!   want old method of searching reference location for the vertical grid
!!   1 for new method (options for other methods 2,... in the future)

105
  real,dimension(0:nymax-1) :: cosf
106
107
  real,dimension(0:nxmax-1,0:nymax-1) ::  tvold,pold,pint,tv
!! automatic arrays introduced in v10 by ?? to achieve better loop order (PS)
108
109
110

  integer,dimension(0:nxmax-1,0:nymax-1) :: rain_cloud_above,idx

111
  integer :: ix,jy,kz,iz,n,kmin,ix1,jy1,ixp,jyp,ixref,jyref,kz_inv
112
113
114
115
116
117
118
119
  real :: f_qvsat,pressure,rh,lsp,convp,cloudh_min,prec
  real :: ew,dz1,dz2,dz
  real :: xlon,ylat,xlonr,dzdx,dzdy
  real :: dzdx1,dzdx2,dzdy1,dzdy2
  real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
  real,parameter :: const=r_air/ga
  real,parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagnostics

120
  logical, save :: init = .true. ! PS 2018-06-13: add back save attribute
121

122
123
124
125
126
127
128
129
130
131
132
133

  !ZHG SEP 2014 tests  
  ! integer :: cloud_ver,cloud_min, cloud_max 
  ! integer ::teller(5), convpteller=0, lspteller=0
  ! real :: cloud_col_wat, cloud_water
  !ZHG 2015 temporary variables for testing
  ! real :: rcw(0:nxmax-1,0:nymax-1)
  ! real :: rpc(0:nxmax-1,0:nymax-1)
  ! character(len=60) :: zhgpath='/xnilu_wrk/flex_wrk/zhg/'
  ! character(len=60) :: fnameA,fnameB,fnameC,fnameD,fnameE,fnameF,fnameG,fnameH
  ! CHARACTER(LEN=3)  :: aspec
  ! integer :: virr=0
134
135
  !real :: tot_cloud_h
  !real :: dbg_height(nzmax) 
136
137
138
139
140
141
!ZHG

!*************************************************************************
! If verttransform is called the first time, initialize heights of the   *
! z levels in meter. The heights are the heights of model levels, where  *
! u,v,T and qv are given, and of the interfaces, where w is given. So,   *
142
! the vertical rESOlution in the z system is doubled. As reference point,*
143
144
145
146
147
! the lower left corner of the grid is used.                             *
! Unlike in the eta system, no difference between heights for u,v and    *
! heights for w exists.                                                  *
!*************************************************************************

148
!ESO measure CPU time
149
150
151
152
!  call mpif_mtime('verttransform',0)

  if (init) then

153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
!! Search for a point with high surface pressure (i.e. not above significant
!!  topography) Then, use this point to construct a reference z profile, 
!!  to be used at all times
! *****************************************************************************

    if (ioldref .eq. 0) then ! old reference grid point
       do jy=0,nymin1
         do ix=0,nxmin1
           if (ps(ix,jy,1,n).gt.1000.e2) then
             ixref=ix
             jyref=jy
             goto 3
           endif
         end do
       end do
3     continue

      print*,'oldheights at' ,ixref,jyref,ps(ixref,jyref,1,n)
    else ! new reference grid point
!     PS: the old version fails if the pressure is <=1000 hPa in the whole
!     domain. Let us find a good replacement, not just a quick fix.
!     Search point near to mean pressure after excluding mountains

      psmean = sum( ps(:,:,1,n) ) / (nx*ny)
      print*,'height: fg psmean',psmean
      psstd = sqrt(sum( (ps(:,:,1,n) - psmean)**2 ) / (nx*ny))

!>    new psmean using only points within $\plusminus\sigma$  
!      psmean = sum( ps(:,:,1,n), abs(ps(:,:,1,n)-psmean) < psstd ) / &
!        count(abs(ps(:,:,1,n)-psmean) < psstd)

!>    new psmean using only points with $p\gt \overbar{p}+\sigma_p$  
!>    (reject mountains, accept valleys)
      psmean = sum( ps(:,:,1,n), ps(:,:,1,n) > psmean - psstd ) / &
        count(ps(:,:,1,n) > psmean - psstd)
      print*,'height: std, new psmean',psstd,psmean
      ixyref = minloc( abs( ps(:,:,1,n) - psmean ) )
      ixref = ixyref(1)
      jyref = ixyref(2)
      print*,'newheights at' ,ixref,jyref,ps(ixref,jyref,1,n)
    endif  

    tvoldref=tt2(ixref,jyref,1,n)* &
      ( 1. + 0.378*ew(td2(ixref,jyref,1,n) ) / ps(ixref,jyref,1,n))
    poldref=ps(ixref,jyref,1,n)
198
199
200
201
    height(1)=0.

    do kz=2,nuvz

202
203
204
205
206
207
208
      pintref = akz(kz)+bkz(kz)*ps(ixref,jyref,1,n)
      tv = tth(ixref,jyref,kz,n)*(1.+0.608*qvh(ixref,jyref,kz,n))

      if (abs(tv(ixref,jyref)-tvoldref) .gt. 0.2) then
        height(kz) = height(kz-1) +      &
          const*log( poldref/pintref ) * &
          ( tv(ixref,jyref) - tvoldref ) / log( tv(ixref,jyref) / tvoldref )
209
      else
210
211
        height(kz) = height(kz-1) + &
          const*log( poldref/pintref ) * tv(ixref,jyref)
212
213
      endif

214
215
216
      tvoldref = tv(ixref,jyref)
      poldref = pintref
      print*,'height=',kz,height(kz),tvoldref,poldref
217

218
    end do
219

220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235

! Determine highest levels that can be within PBL
!************************************************

    do kz=1,nz
      if (height(kz).gt.hmixmax) then
        nmixz=kz
        goto 9
      endif
    end do
9   continue

! Do not repeat initialization of the Cartesian z grid
!*****************************************************
    init=.false.

236
  endif ! init block (vertical grid construction)
237
238
239
240
241


! Loop over the whole grid
!*************************

242
  tvold(:,:)=tt2(:,:,1,n) * (1.+0.378*ew( td2(:,:,1,n) ) / ps(:,:,1,n))
243
244
245
246
247
248
249
250
251
252
  pold=ps(:,:,1,n)
  uvzlev(:,:,1)=0.
  wzlev(:,:,1)=0.
  rhoh(:,:,1)=pold/(r_air*tvold)


! Compute heights of eta levels
!******************************

  do kz=2,nuvz
253
254
255
    pint(:,:)=akz(kz)+bkz(kz)*ps(:,:,1,n)
    tv(:,:)=tth(:,:,kz,n)*(1.+0.608*qvh(:,:,kz,n))
    rhoh(:,:,kz)=pint(:,:)/(r_air*tv(:,:))
256

257
258
259
    where (abs(tv(:,:)-tvold(:,:)).gt.0.2) 
      uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+const*log(pold(:,:)/pint(:,:))* &
           (tv(:,:)-tvold(:,:))/log(tv(:,:)/tvold(:,:))
260
    elsewhere
261
      uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+const*log(pold(:,:)/pint(:,:))*tv(:,:)
262
263
    endwhere

264
265
    tvold(:,:)=tv(:,:)
    pold(:,:)=pint(:,:)
266
267
268
269
270
271
272
273
274
275
276
277
  end do


  do kz=2,nwz-1
    wzlev(:,:,kz)=(uvzlev(:,:,kz+1)+uvzlev(:,:,kz))/2.
  end do
  wzlev(:,:,nwz)=wzlev(:,:,nwz-1)+ &
       uvzlev(:,:,nuvz)-uvzlev(:,:,nuvz-1)

! pinmconv=(h2-h1)/(p2-p1)

  pinmconv(:,:,1)=(uvzlev(:,:,2))/ &
278
279
    ((aknew(2)+bknew(2)*ps(:,:,1,n))- &
     (aknew(1)+bknew(1)*ps(:,:,1,n)))
280
281
  do kz=2,nz-1
    pinmconv(:,:,kz)=(uvzlev(:,:,kz+1)-uvzlev(:,:,kz-1))/ &
282
283
      ((aknew(kz+1)+bknew(kz+1)*ps(:,:,1,n))- &
       (aknew(kz-1)+bknew(kz-1)*ps(:,:,1,n)))
284
285
286
  end do
  pinmconv(:,:,nz)=(uvzlev(:,:,nz)-uvzlev(:,:,nz-1))/ &
       ((aknew(nz)+bknew(nz)*ps(:,:,1,n))- &
287
        (aknew(nz-1)+bknew(nz-1)*ps(:,:,1,n)))
288

289
290
! Levels where u,v,t and q are given
!***********************************
291
292
293
294
295

  uu(:,:,1,n)=uuh(:,:,1)
  vv(:,:,1,n)=vvh(:,:,1)
  tt(:,:,1,n)=tth(:,:,1,n)
  qv(:,:,1,n)=qvh(:,:,1,n)
296
!HG adding the cloud water 
297
298
299
300
  if (readclouds) then
    clwc(:,:,1,n)=clwch(:,:,1,n)
    if (.not.sumclouds) ciwc(:,:,1,n)=ciwch(:,:,1,n)
  end if
301
!HG 
302
303
304
305
306
307
308
  pv(:,:,1,n)=pvh(:,:,1)
  rho(:,:,1,n)=rhoh(:,:,1)

  uu(:,:,nz,n)=uuh(:,:,nuvz)
  vv(:,:,nz,n)=vvh(:,:,nuvz)
  tt(:,:,nz,n)=tth(:,:,nuvz,n)
  qv(:,:,nz,n)=qvh(:,:,nuvz,n)
309
!HG adding the cloud water
310
311
  if (readclouds) then
    clwc(:,:,nz,n)=clwch(:,:,nuvz,n)
312
    if (.not. sumclouds) ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n)
313
  end if
314
!HG
315
316
317
318
319
  pv(:,:,nz,n)=pvh(:,:,nuvz)
  rho(:,:,nz,n)=rhoh(:,:,nuvz)

  kmin=2
  idx=kmin
320
  iz_loop: do iz=2,nz-1
321
322
    do jy=0,nymin1
      do ix=0,nxmin1
323
324
        if (height(iz).gt.uvzlev(ix,jy,nuvz)) then
        
325
326
327
328
          uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
          vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
          tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
          qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
329
!HG adding the cloud water
330
331
          if (readclouds) then
            clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
332
            if (.not. sumclouds) ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
333
          end if
334
!HG
335
336
          pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
          rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
337

338
        else
339
340
341
342

          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
343
              idx(ix,jy)=kz
344
              exit kz_loop
345
            endif
346
347
          enddo kz_loop

348
349
350
351
352
        endif
      enddo
    enddo
    do jy=0,nymin1
      do ix=0,nxmin1
353
354
        if (height(iz).le.uvzlev(ix,jy,nuvz)) then
        
355
356
357
358
359
360
361
362
363
364
          kz=idx(ix,jy)
          dz1=height(iz)-uvzlev(ix,jy,kz-1)
          dz2=uvzlev(ix,jy,kz)-height(iz)
          dz=dz1+dz2
          uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
          vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
          tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
               +tth(ix,jy,kz,n)*dz1)/dz
          qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
               +qvh(ix,jy,kz,n)*dz1)/dz
365
!HG adding the cloud water
366
367
          if (readclouds) then
            clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2+clwch(ix,jy,kz,n)*dz1)/dz
368
369
            if (.not. sumclouds) ciwc(ix,jy,iz,n)= &
              (ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz
370
          end if
371
!HG
372
373
          pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
          rho(ix,jy,iz,n)=(rhoh(ix,jy,kz-1)*dz2+rhoh(ix,jy,kz)*dz1)/dz
374

375
376
377
        endif
      enddo
    enddo
378
  enddo iz_loop
379
380


381
! Levels where w is given
382
383
384
385
386
387
!*************************

  ww(:,:,1,n)=wwh(:,:,1)*pinmconv(:,:,1)
  ww(:,:,nz,n)=wwh(:,:,nwz)*pinmconv(:,:,nz)
  kmin=2
  idx=kmin
388
  iz_loop2: do iz=2,nz
389
390
    do jy=0,nymin1
      do ix=0,nxmin1
391
392
393
        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
394
            idx(ix,jy)=kz
395
            exit kz_loop2
396
          endif
397
        enddo kz_loop2
398
399
400
401
402
403
404
405
406
      enddo
    enddo
    do jy=0,nymin1
      do ix=0,nxmin1
        kz=idx(ix,jy)
        dz1=height(iz)-wzlev(ix,jy,kz-1)
        dz2=wzlev(ix,jy,kz)-height(iz)
        dz=dz1+dz2
        ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(ix,jy,kz-1)*dz2 &
407
                      + wwh(ix,jy,kz)  *pinmconv(ix,jy,kz)*dz1)/dz
408
409
      enddo
    enddo
410
  enddo iz_loop2
411
412
413
414

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

415
  drhodz(:,:,1,n)=(rho(:,:,2,n)-rho(:,:,1,n))/(height(2)-height(1))
416
417
  do kz=2,nz-1
    drhodz(:,:,kz,n)=(rho(:,:,kz+1,n)-rho(:,:,kz-1,n))/ &
418
      (height(kz+1)-height(kz-1))
419
420
421
422
423
424
425
426
427
428
  end do
  drhodz(:,:,nz,n)=drhodz(:,:,nz-1,n)


!****************************************************************
! Compute slope of eta levels in windward direction and resulting
! vertical wind correction
!****************************************************************

  do jy=1,ny-2
429
    cosf(jy) = 1. / cos( ( real(jy)*dy + ylat0 )*pi180 )
430
431
432
433
  enddo

  kmin=2
  idx=kmin
434
  iz_loop3: do iz=2,nz-1
435
436
437
    do jy=1,ny-2
      do ix=1,nx-2

438
439
440
        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
441
            idx(ix,jy)=kz
442
            exit kz_loop3
443
          endif
444
        enddo kz_loop3
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
      enddo
    enddo

    do jy=1,ny-2
      do ix=1,nx-2
        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.
460
        dzdx2=(uvzlev(ixp,jy,kz)  -uvzlev(ix1,jy,kz)  )/2.
461
462
463
        dzdx=(dzdx1*dz2+dzdx2*dz1)/dz

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

467
468
        ww(ix,jy,iz,n)=ww(ix,jy,iz,n) +( dzdx*uu(ix,jy,iz,n)*dxconst*cosf(jy) &
          + dzdy*vv(ix,jy,iz,n)*dyconst)
469

470
471
      enddo
    enddo
472

473
  enddo iz_loop3
474
475
476
477
478
479

! If north pole is in the domain, calculate wind velocities in polar
! stereographic coordinates
!*******************************************************************

  if (nglobal) then
480
    
481
482
483
484
485
486
    do iz=1,nz
      do jy=int(switchnorthg)-2,nymin1
        ylat=ylat0+real(jy)*dy
        do ix=0,nxmin1
          xlon=xlon0+real(ix)*dx
          call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
487
488
489
490
               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), vvpol(ix,jy,iz,n))
        enddo
      enddo
    enddo
491
492
493
494


    do iz=1,nz

495
!   CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
496
!
497
498
!   AMS nauffer Nov 18 2004 Added check for case vv=0

499
500
      xlon=xlon0+real(nx/2-1)*dx
      xlonr=xlon*pi/180.
501
      ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+vv(nx/2-1,nymin1,iz,n)**2)
502
      if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
503
        ddpol=atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr
504
      else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
505
        ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr
506
      else
507
        ddpol=pi/2.-xlonr
508
      endif
509
510
      if (ddpol.lt.0.) ddpol=2.0*pi+ddpol
      if (ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
511
512
513
514
515
516
517

! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
      xlon=180.0
      xlonr=xlon*pi/180.
      ylat=90.0
      uuaux=-ffpol*sin(xlonr+ddpol)
      vvaux=-ffpol*cos(xlonr+ddpol)
518
      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, vvpolaux)
519
520
521
522
523

      jy=nymin1
      do ix=0,nxmin1
        uupol(ix,jy,iz,n)=uupolaux
        vvpol(ix,jy,iz,n)=vvpolaux
524
525
      enddo
    enddo
526

527
528
! Fix: Set W (vertical wind) at pole to the zonally averaged W of the next 
! equator-ward parallel 
529
530
531
532
533
534

    do iz=1,nz
      wdummy=0.
      jy=ny-2
      do ix=0,nxmin1
        wdummy=wdummy+ww(ix,jy,iz,n)
535
      enddo
536
537
538
539
      wdummy=wdummy/real(nx)
      jy=nymin1
      do ix=0,nxmin1
        ww(ix,jy,iz,n)=wdummy
540
541
      enddo
    enddo
542
543
544
545
546
547
548
549
550

  endif


! If south pole is in the domain, calculate wind velocities in polar
! stereographic coordinates
!*******************************************************************

  if (sglobal) then
551
  
552
553
554
555
556
557
    do iz=1,nz
      do jy=0,int(switchsouthg)+3
        ylat=ylat0+real(jy)*dy
        do ix=0,nxmin1
          xlon=xlon0+real(ix)*dx
          call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
558
559
560
561
               vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
        enddo
      enddo
    enddo
562
563
564

    do iz=1,nz

565
!   CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
566
!
567
568
!   AMS nauffer Nov 18 2004 Added check for case vv=0

569
570
      xlon=xlon0+real(nx/2-1)*dx
      xlonr=xlon*pi/180.
571
      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+vv(nx/2-1,0,iz,n)**2)
572
      if (vv(nx/2-1,0,iz,n).lt.0.) then
573
        ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
574
      else if (vv(nx/2-1,0,iz,n).gt.0.) then
575
        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
576
      else
577
        ddpol=pi/2.-xlonr
578
      endif
579
580
      if (ddpol.lt.0.) ddpol=2.0*pi+ddpol
      if (ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
581
582
583
584
585
586
587

! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
      xlon=180.0
      xlonr=xlon*pi/180.
      ylat=-90.0
      uuaux=+ffpol*sin(xlonr-ddpol)
      vvaux=-ffpol*cos(xlonr-ddpol)
588
      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
589
590
591
592
593

      jy=0
      do ix=0,nxmin1
        uupol(ix,jy,iz,n)=uupolaux
        vvpol(ix,jy,iz,n)=vvpolaux
594
595
      enddo
    enddo
596

597
598
! Fix: Set W (vertical wind) at pole to the zonally averaged W of the next 
! equator-ward parallel 
599
600
601
602
603
604

    do iz=1,nz
      wdummy=0.
      jy=1
      do ix=0,nxmin1
        wdummy=wdummy+ww(ix,jy,iz,n)
605
      enddo
606
607
608
609
      wdummy=wdummy/real(nx)
      jy=0
      do ix=0,nxmin1
        ww(ix,jy,iz,n)=wdummy
610
611
      enddo
    enddo
612
613
614
  endif


615
616
617
618
619
620
!****************************************************************************** 
  if (readclouds) 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
621
! cloud scavenging are assigned with numbers 2-5 (following the old metod).
622
623
624
625
626
627
! 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(*,*) 'Global ECMWF fields: using cloud water'
628
629
630
631
    clw(:,:,:,n)=0.0
    ctwc(:,:,n)=0.0
    clouds(:,:,:,n)=0
! If water/ice are read separately into clwc and ciwc, store sum in clwc
632
    if (.not. sumclouds) then 
633
      clwc(:,:,:,n) = clwc(:,:,:,n) + ciwc(:,:,:,n)
634
    endif
635
636
637
638
    do jy=0,nymin1
      do ix=0,nxmin1
        lsp=lsprec(ix,jy,1,n)
        convp=convprec(ix,jy,1,n)
639
        prec=lsp+convp  ! Note PS: prec is not used currently
640
!        tot_cloud_h=0
641
! Find clouds in the vertical
642
!! Note PS: bad loop order. 
643
644
645
        do kz=1, nz-1 !go from top to bottom
          if (clwc(ix,jy,kz,n).gt.0) then      
! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3 
646
647
            clw(ix,jy,kz,n)=(clwc(ix,jy,kz,n)*rho(ix,jy,kz,n))* &
              (height(kz+1)-height(kz))
648
!            tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz)) 
649
650
651
652
653
654
655
656
657
658
659
660
661
662
            
!            icloud_stats(ix,jy,4,n)= icloud_stats(ix,jy,4,n)+clw(ix,jy,kz,n)          ! Column cloud water [m3/m3]
            ctwc(ix,jy,n) = ctwc(ix,jy,n)+clw(ix,jy,kz,n)
!            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
663
664
        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
665
666

          do kz=nz,1,-1 !go Bottom up!
667
!! Note PS: bad loop order
668
669
670
671
672
673
674
675
            if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud
              cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1) 
              clouds(ix,jy,kz,n)=1                               ! is a cloud
              if (lsp.ge.convp) then
                clouds(ix,jy,kz,n)=3                            ! lsp in-cloud
              else
                clouds(ix,jy,kz,n)=2                             ! convp in-cloud
              endif                                              ! convective or large scale
676
            elseif( clw(ix,jy,kz,n).le.0 .and. cloudh_min.ge.height(kz)) then ! is below cloud
677
678
679
680
681
682
683
684
              if (lsp.ge.convp) then
                clouds(ix,jy,kz,n)=5                             ! lsp dominated washout
              else
                clouds(ix,jy,kz,n)=4                             ! convp dominated washout
              endif                                              ! convective or large scale 
            endif

            if (height(kz).ge. 19000) then                        ! set a max height for removal
685
!! Note PS: such hardcoded limits would better be parameters defined in par_mod
686
              clouds(ix,jy,kz,n)=0
687
688
            endif ! clw>0
          enddo ! kz
689
690
691
692
        endif ! precipitation
      end do
    end do

693
! ESO: copy the relevant data to clw4 to reduce amount of communicated data for MPI
694
695
696
697
698
!    ctwc(:,:,n) = icloud_stats(:,:,4,n)

!**************************************************************************
  else       ! use old definitions
!**************************************************************************
699

700
701
!   create a cloud and rainout/washout field, clouds occur where rh>80%
!   total cloudheight is stored at level 0
702
703

!PS    write(*,*) 'Global fields: using cloud water from Parameterization'
704
705
706
707
708
709
710
    do jy=0,nymin1
      do ix=0,nxmin1
        rain_cloud_above(ix,jy)=0
        lsp=lsprec(ix,jy,1,n)
        convp=convprec(ix,jy,1,n)
        cloudsh(ix,jy,n)=0
        do kz_inv=1,nz-1
711
!! Note PS: bad loop order. 
712
713
714
715
          kz=nz-kz_inv+1
          pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
          rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
          clouds(ix,jy,kz,n)=0
716
717
718
719
720
721

          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
722
723
724
725
726
727
728
729
730
731
732
              rain_cloud_above(ix,jy)=1
              cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
                   height(kz)-height(kz-1)
              if (lsp.ge.convp) then
                clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
              else
                clouds(ix,jy,kz,n)=2 ! convp dominated rainout
              endif
            else ! no precipitation
              clouds(ix,jy,kz,n)=1 ! cloud
            endif
733

734
          else ! no cloud
735

736
737
738
739
740
741
742
            if (rain_cloud_above(ix,jy).eq.1) then ! scavenging
              if (lsp.ge.convp) then
                clouds(ix,jy,kz,n)=5 ! lsp dominated washout
              else
                clouds(ix,jy,kz,n)=4 ! convp dominated washout
              endif
            endif
743

744
745
          endif
        end do
746

747
748
      end do
    end do
749
  endif ! END OLD METHOD
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892


     !********* TEST ***************
     ! WRITE OUT SOME TEST VARIABLES
     !********* TEST ************'**
!teller(:)=0
!virr=virr+1
!WRITE(aspec, '(i3.3)'), virr

!if (readclouds) then
!fnameH=trim(zhgpath)//trim(aspec)//'Vertical_placement.txt'
!else
!fnameH=trim(zhgpath)//trim(aspec)//'Vertical_placement_old.txt'
!endif
!
!OPEN(UNIT=118, FILE=fnameH,FORM='FORMATTED',STATUS = 'UNKNOWN')
!do kz_inv=1,nz-1
!  kz=nz-kz_inv+1
!  !kz=91
!  do jy=0,nymin1
!     do ix=0,nxmin1
!          if (clouds(ix,jy,kz,n).eq.1) teller(1)=teller(1)+1 ! no precipitation cloud
!          if (clouds(ix,jy,kz,n).eq.2) teller(2)=teller(2)+1 ! convp dominated rainout
!          if (clouds(ix,jy,kz,n).eq.3) teller(3)=teller(3)+1 ! lsp dominated rainout
!          if (clouds(ix,jy,kz,n).eq.4) teller(4)=teller(4)+1 ! convp dominated washout
!          if (clouds(ix,jy,kz,n).eq.5) teller(5)=teller(5)+1 ! lsp dominated washout
!          
!        !  write(*,*) height(kz),teller
!     end do
!  end do
!  write(118,*) height(kz),teller
!  teller(:)=0
!end do
!teller(:)=0
!write(*,*) teller 
!write(*,*) aspec
!
!fnameA=trim(zhgpath)//trim(aspec)//'cloudV.txt'
!fnameB=trim(zhgpath)//trim(aspec)//'cloudT.txt'
!fnameC=trim(zhgpath)//trim(aspec)//'cloudB.txt'
!fnameD=trim(zhgpath)//trim(aspec)//'cloudW.txt'
!fnameE=trim(zhgpath)//trim(aspec)//'old_cloudV.txt'
!fnameF=trim(zhgpath)//trim(aspec)//'lsp.txt'
!fnameG=trim(zhgpath)//trim(aspec)//'convp.txt'
!if (readclouds) then
!OPEN(UNIT=111, FILE=fnameA,FORM='FORMATTED',STATUS = 'UNKNOWN')
!OPEN(UNIT=112, FILE=fnameB,FORM='FORMATTED',STATUS = 'UNKNOWN')
!OPEN(UNIT=113, FILE=fnameC,FORM='FORMATTED',STATUS = 'UNKNOWN')
!OPEN(UNIT=114, FILE=fnameD,FORM='FORMATTED',STATUS = 'UNKNOWN')
!else
!OPEN(UNIT=115, FILE=fnameE,FORM='FORMATTED',STATUS = 'UNKNOWN')
!OPEN(UNIT=116, FILE=fnameF,FORM='FORMATTED',STATUS = 'UNKNOWN')
!OPEN(UNIT=117, FILE=fnameG,FORM='FORMATTED',STATUS = 'UNKNOWN')
!endif
!
!do ix=0,nxmin1
!if (readclouds) then
!write(111,*) (icloud_stats(ix,jy,1,n),jy=0,nymin1)
!write(112,*) (icloud_stats(ix,jy,2,n),jy=0,nymin1)
!write(113,*) (icloud_stats(ix,jy,3,n),jy=0,nymin1)
!write(114,*) (icloud_stats(ix,jy,4,n),jy=0,nymin1)
!else
!write(115,*) (cloudsh(ix,jy,n),jy=0,nymin1)    !integer
!write(116,*) (lsprec(ix,jy,1,n),jy=0,nymin1)   !7.83691406E-02 
!write(117,*) (convprec(ix,jy,1,n),jy=0,nymin1) !5.38330078E-02
!endif
!end do
!
!if (readclouds) then
!CLOSE(111)
!CLOSE(112)
!CLOSE(113)
!CLOSE(114)
!else
!CLOSE(115)
!CLOSE(116)
!CLOSE(117)
!endif
!
!END ********* TEST *************** END
! WRITE OUT SOME TEST VARIABLES
!END ********* TEST *************** END


! PS 2012
!      lsp=lsprec(ix,jy,1,n)
!      convp=convprec(ix,jy,1,n)
!          prec=lsp+convp
!          if (lsp.gt.convp) then !  prectype='lsp'
!            lconvectprec = .false.
!          else ! prectype='cp'
!            lconvectprec = .true.
!           endif
!      else ! windfields does not contain cloud data 
!          rhmin = 0.90 ! standard condition for presence of clouds
!PS       note that original by Sabine Eckhart was 80%
!PS       however, for T<-20 C we consider saturation over ice
!PS       so I think 90% should be enough          
!          icloudbot(ix,jy,n)=icmv
!          icloudtop=icmv ! this is just a local variable
!98        do kz=1,nz
!            pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
!            rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
!ps            if (prec.gt.0.01) print*,'relhum',prec,kz,rh,height(kz)
!            if (rh .gt. rhmin) then
!              if (icloudbot(ix,jy,n) .eq. icmv) then
!                icloudbot(ix,jy,n)=nint(height(kz))
!              endif
!              icloudtop=nint(height(kz)) ! use int to save memory
!            endif
    ! end do
!PS try to get a cloud thicker than 50 m 
!PS if there is at least .01 mm/h  - changed to 0.002 and put into
!PS parameter precpmin        
!          if ((icloudbot(ix,jy,n) .eq. icmv .or. &
!              icloudtop-icloudbot(ix,jy,n) .lt. 50) .and. &
!              prec .gt. precmin) then
!            rhmin = rhmin - 0.05
!            if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum.
!   end if

!PS is based on looking at a limited set of comparison data
!          if (lconvectprec .and. icloudtop .lt. 6000 .and. &
!             prec .gt. precmin) then 
!
!            if (convp .lt. 0.1) then
!              icloudbot(ix,jy,n) = 500
!              icloudtop =         8000
!            else
!              icloudbot(ix,jy,n) = 0
!              icloudtop =      10000
!            endif
!          endif
!          if (icloudtop .ne. icmv) then
!            icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n)
!          else
!            icloudthck(ix,jy,n) = icmv
!          endif
!PS  get rid of too thin clouds      
!          if (icloudthck(ix,jy,n) .lt. 50) then
!            icloudbot(ix,jy,n)=icmv
!            icloudthck(ix,jy,n)=icmv
!          endif
893
!HG__________________________________
894
895
!           rcw(ix,jy)=2E-7*prec**0.36
!           rpc(ix,jy)=prec
896
!HG end______________________________
897

898
!      endif !HG read clouds
899
900
901
902




903
!ESO measure CPU time
904
905
!  call mpif_mtime('verttransform',1)

906
!ESO print out the same measure as in Leo's routine
907
908
909
910
911
912
913
914
915
    ! write(*,*) 'verttransform: ', &
    !      sum(tt(:,:,:,n)*tt(:,:,:,n)), &
    !      sum(uu(:,:,:,n)*uu(:,:,:,n)),sum(vv(:,:,:,n)*vv(:,:,:,n)),&
    !      sum(qv(:,:,:,n)*qv(:,:,:,n)),sum(pv(:,:,:,n)*pv(:,:,:,n)),&
    !      sum(rho(:,:,:,n)*rho(:,:,:,n)),sum(drhodz(:,:,:,n)*drhodz(:,:,:,n)),&
    !      sum(ww(:,:,:,n)*ww(:,:,:,n)), &
    !      sum(clouds(:,:,:,n)), sum(cloudsh(:,:,n)),sum(idx),sum(pinmconv)
end subroutine verttransform_ecmwf