Index: src/verttransform_ecmwf.f90
===================================================================
--- src/verttransform_ecmwf.f90	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
+++ src/verttransform_ecmwf.f90	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
@@ -0,0 +1,866 @@
+!**********************************************************************
+! 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_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                                                *
+!                                                                            *
+!     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:
+!   Variables tth and qvh (on eta coordinates) from common block
+!
+! Sabine Eckhardt, March 2007
+! added the variable cloud for use with scavenging - descr. in com_mod
+!
+! Unified ECMWF and GFS builds
+! Marian Harustak, 12.5.2017 
+!     - Renamed from verttransform to verttransform_ecmwf
+!
+!*****************************************************************************
+!                                                                            *
+! 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
+!  use mpi_mod
+
+  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
+  real,dimension(0:nxmax-1,0:nymax-1) ::  tvold,pold,pint,tv
+  real,dimension(0:nymax-1) :: cosf
+
+  integer,dimension(0:nxmax-1,0:nymax-1) :: rain_cloud_above,idx
+
+  integer :: ix,jy,kz,iz,n,kmin,ix1,jy1,ixp,jyp,ixm,jym,kz_inv
+  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
+
+  logical :: init = .true.
+
+  !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
+  real :: tot_cloud_h
+  real :: dbg_height(nzmax) 
+!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,   *
+! the vertical resolution in the z system is doubled. As reference point,*
+! 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.                                                  *
+!*************************************************************************
+
+
+!eso measure CPU time
+!  call mpif_mtime('verttransform',0)
+
+  if (init) then
+
+! 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
+!*****************************************************************************
+
+    do jy=0,nymin1
+      do ix=0,nxmin1
+        if (ps(ix,jy,1,n).gt.100000.) then
+          ixm=ix
+          jym=jy
+          goto 3
+        endif
+      end do
+    end do
+3   continue
+
+
+    tvold(ixm,jym)=tt2(ixm,jym,1,n)*(1.+0.378*ew*(td2(ixm,jym,1,n))/ &
+         ps(ixm,jym,1,n))
+    pold(ixm,jym)=ps(ixm,jym,1,n)
+    height(1)=0.
+
+    do kz=2,nuvz
+      pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
+      tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
+
+      if (abs(tv(ixm,jym)-tvold(ixm,jym)).gt.0.2) then
+        height(kz)= &
+             height(kz-1)+const*log(pold(ixm,jym)/pint(ixm,jym))* &
+             (tv(ixm,jym)-tvold(ixm,jym))/log(tv(ixm,jym)/tvold(ixm,jym))
+      else
+        height(kz)=height(kz-1)+ &
+             const*log(pold(ixm,jym)/pint(ixm,jym))*tv(ixm,jym)
+      endif
+
+      tvold(ixm,jym)=tv(ixm,jym)
+      pold(ixm,jym)=pint(ixm,jym)
+    end do
+
+
+! 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.
+
+    dbg_height = height
+
+  endif
+
+
+! Loop over the whole grid
+!*************************
+
+
+  do jy=0,nymin1
+    do ix=0,nxmin1
+      tvold(ix,jy)=tt2(ix,jy,1,n)*(1.+0.378*ew*(td2(ix,jy,1,n))/ &
+           ps(ix,jy,1,n))
+    enddo
+  enddo
+  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
+    pint=akz(kz)+bkz(kz)*ps(:,:,1,n)
+    tv=tth(:,:,kz,n)*(1.+0.608*qvh(:,:,kz,n))
+    rhoh(:,:,kz)=pint(:,:)/(r_air*tv)
+
+    where (abs(tv-tvold).gt.0.2) 
+      uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+const*log(pold/pint)* &
+           (tv-tvold)/log(tv/tvold)
+    elsewhere
+      uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+const*log(pold/pint)*tv
+    endwhere
+
+    tvold=tv
+    pold=pint
+  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))/ &
+       ((aknew(2)+bknew(2)*ps(:,:,1,n))- &
+       (aknew(1)+bknew(1)*ps(:,:,1,n)))
+  do kz=2,nz-1
+    pinmconv(:,:,kz)=(uvzlev(:,:,kz+1)-uvzlev(:,:,kz-1))/ &
+         ((aknew(kz+1)+bknew(kz+1)*ps(:,:,1,n))- &
+         (aknew(kz-1)+bknew(kz-1)*ps(:,:,1,n)))
+  end do
+  pinmconv(:,:,nz)=(uvzlev(:,:,nz)-uvzlev(:,:,nz-1))/ &
+       ((aknew(nz)+bknew(nz)*ps(:,:,1,n))- &
+       (aknew(nz-1)+bknew(nz-1)*ps(:,:,1,n)))
+
+! Levels, where u,v,t and q are given
+!************************************
+
+
+  uu(:,:,1,n)=uuh(:,:,1)
+  vv(:,:,1,n)=vvh(:,:,1)
+  tt(:,:,1,n)=tth(:,:,1,n)
+  qv(:,:,1,n)=qvh(:,:,1,n)
+!hg adding the cloud water 
+  if (readclouds) then
+    clwc(:,:,1,n)=clwch(:,:,1,n)
+    if (.not.sumclouds) ciwc(:,:,1,n)=ciwch(:,:,1,n)
+  end if
+!hg 
+  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)
+!hg adding the cloud water
+  if (readclouds) then
+    clwc(:,:,nz,n)=clwch(:,:,nuvz,n)
+    if (.not.sumclouds) ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n)
+  end if
+!hg
+  pv(:,:,nz,n)=pvh(:,:,nuvz)
+  rho(:,:,nz,n)=rhoh(:,:,nuvz)
+
+
+  kmin=2
+  idx=kmin
+  do iz=2,nz-1
+    do jy=0,nymin1
+      do ix=0,nxmin1
+        if(height(iz).gt.uvzlev(ix,jy,nuvz)) then
+          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)
+!hg adding the cloud water
+          if (readclouds) then
+            clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
+            if (.not.sumclouds) ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
+          end if
+!hg
+          pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
+          rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
+        else
+          innuvz: 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
+              idx(ix,jy)=kz
+              exit innuvz
+            endif
+          enddo innuvz
+        endif
+      enddo
+    enddo
+    do jy=0,nymin1
+      do ix=0,nxmin1
+        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)
+          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
+!hg adding the cloud water
+          if (readclouds) then
+            clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2+clwch(ix,jy,kz,n)*dz1)/dz
+            if (.not.sumclouds) &
+                 &ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz
+          end if
+!hg
+          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
+        endif
+      enddo
+    enddo
+  enddo
+
+
+! Levels, where w is given
+!*************************
+
+  ww(:,:,1,n)=wwh(:,:,1)*pinmconv(:,:,1)
+  ww(:,:,nz,n)=wwh(:,:,nwz)*pinmconv(:,:,nz)
+  kmin=2
+  idx=kmin
+  do iz=2,nz
+    do jy=0,nymin1
+      do ix=0,nxmin1
+        inn:         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
+            idx(ix,jy)=kz
+            exit inn
+          endif
+        enddo inn
+      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 &
+             +wwh(ix,jy,kz)*pinmconv(ix,jy,kz)*dz1)/dz
+      enddo
+    enddo
+  enddo
+
+! Compute density gradients at intermediate levels
+!*************************************************
+
+  drhodz(:,:,1,n)=(rho(:,:,2,n)-rho(:,:,1,n))/ &
+       (height(2)-height(1))
+  do kz=2,nz-1
+    drhodz(:,:,kz,n)=(rho(:,:,kz+1,n)-rho(:,:,kz-1,n))/ &
+         (height(kz+1)-height(kz-1))
+  end do
+  drhodz(:,:,nz,n)=drhodz(:,:,nz-1,n)
+
+!    end do
+!  end do
+
+
+!****************************************************************
+! Compute slope of eta levels in windward direction and resulting
+! vertical wind correction
+!****************************************************************
+
+  do jy=1,ny-2
+    cosf(jy)=1./cos((real(jy)*dy+ylat0)*pi180)
+  enddo
+
+  kmin=2
+  idx=kmin
+  do iz=2,nz-1
+    do jy=1,ny-2
+      do ix=1,nx-2
+
+        inneta: 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
+            idx(ix,jy)=kz
+            exit inneta
+          endif
+        enddo inneta
+      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.
+        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
+
+        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)
+
+      end do
+
+    end do
+  end do
+
+! If north pole is in the domain, calculate wind velocities in polar
+! stereographic coordinates
+!*******************************************************************
+
+  if (nglobal) then
+    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), &
+               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
+               vvpol(ix,jy,iz,n))
+        end do
+      end do
+    end do
+
+
+    do iz=1,nz
+
+! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
+!
+!   AMSnauffer Nov 18 2004 Added check for case vv=0
+!
+      xlon=xlon0+real(nx/2-1)*dx
+      xlonr=xlon*pi/180.
+      ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
+           vv(nx/2-1,nymin1,iz,n)**2)
+      if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
+        ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
+             vv(nx/2-1,nymin1,iz,n))-xlonr
+      else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
+        ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
+             vv(nx/2-1,nymin1,iz,n))-xlonr
+      else
+        ddpol=pi/2-xlonr
+      endif
+      if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
+      if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
+
+! 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)
+      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
+           vvpolaux)
+
+      jy=nymin1
+      do ix=0,nxmin1
+        uupol(ix,jy,iz,n)=uupolaux
+        vvpol(ix,jy,iz,n)=vvpolaux
+      end do
+    end do
+
+
+! Fix: Set W at pole to the zonally averaged W of the next equator-
+! ward parallel of latitude
+
+    do iz=1,nz
+      wdummy=0.
+      jy=ny-2
+      do ix=0,nxmin1
+        wdummy=wdummy+ww(ix,jy,iz,n)
+      end do
+      wdummy=wdummy/real(nx)
+      jy=nymin1
+      do ix=0,nxmin1
+        ww(ix,jy,iz,n)=wdummy
+      end do
+    end do
+
+  endif
+
+
+! If south pole is in the domain, calculate wind velocities in polar
+! stereographic coordinates
+!*******************************************************************
+
+  if (sglobal) then
+    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), &
+               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
+               vvpol(ix,jy,iz,n))
+        end do
+      end do
+    end do
+
+    do iz=1,nz
+
+! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
+!
+!   AMSnauffer Nov 18 2004 Added check for case vv=0
+!
+      xlon=xlon0+real(nx/2-1)*dx
+      xlonr=xlon*pi/180.
+      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
+           vv(nx/2-1,0,iz,n)**2)
+      if (vv(nx/2-1,0,iz,n).lt.0.) then
+        ddpol=atan(uu(nx/2-1,0,iz,n)/ &
+             vv(nx/2-1,0,iz,n))+xlonr
+      else if (vv(nx/2-1,0,iz,n).gt.0.) then
+        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
+             vv(nx/2-1,0,iz,n))+xlonr
+      else
+        ddpol=pi/2-xlonr
+      endif
+      if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
+      if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
+
+! 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)
+      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
+           vvpolaux)
+
+      jy=0
+      do ix=0,nxmin1
+        uupol(ix,jy,iz,n)=uupolaux
+        vvpol(ix,jy,iz,n)=vvpolaux
+      end do
+    end do
+
+
+! Fix: Set W at pole to the zonally averaged W of the next equator-
+! ward parallel of latitude
+
+    do iz=1,nz
+      wdummy=0.
+      jy=1
+      do ix=0,nxmin1
+        wdummy=wdummy+ww(ix,jy,iz,n)
+      end do
+      wdummy=wdummy/real(nx)
+      jy=0
+      do ix=0,nxmin1
+        ww(ix,jy,iz,n)=wdummy
+      end do
+    end do
+  endif
+
+
+!***********************************************************************************  
+  if (readclouds) then !HG METHOD
+! The method is loops all grids vertically and constructs the 3D matrix for clouds
+! Cloud top and cloud bottom gid cells are assigned as well as the total column 
+! cloud water. For precipitating grids, the type and whether it is in or below 
+! cloud scavenging are assigned with numbers 2-5 (following the old metod).
+! Distinction is done for lsp and convp though they are treated the same in regards
+! to scavenging. Also clouds that are not precipitating are defined which may be 
+! to include future cloud processing by non-precipitating-clouds. 
+!***********************************************************************************
+    write(*,*) 'Global ECMWF fields: using cloud water'
+    clw(:,:,:,n)=0.0
+!    icloud_stats(:,:,:,n)=0.0
+    ctwc(:,:,n)=0.0
+    clouds(:,:,:,n)=0
+! If water/ice are read separately into clwc and ciwc, store sum in clwc
+    if (.not.sumclouds) then 
+      clwc(:,:,:,n) = clwc(:,:,:,n) + ciwc(:,:,:,n)
+    end if
+    do jy=0,nymin1
+      do ix=0,nxmin1
+        lsp=lsprec(ix,jy,1,n)
+        convp=convprec(ix,jy,1,n)
+        prec=lsp+convp
+        tot_cloud_h=0
+! Find clouds in the vertical
+        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 
+            clw(ix,jy,kz,n)=(clwc(ix,jy,kz,n)*rho(ix,jy,kz,n))*(height(kz+1)-height(kz))
+            tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz)) 
+            
+!            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
+        if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
+
+          do kz=nz,1,-1 !go Bottom up!
+            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
+            elseif((clw(ix,jy,kz,n).le.0) .and. (cloudh_min.ge.height(kz))) then ! is below cloud
+              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
+              clouds(ix,jy,kz,n)=0
+            endif !clw>0
+          end do !nz
+        endif ! precipitation
+      end do
+    end do
+
+! eso: copy the relevant data to clw4 to reduce amount of communicated data for MPI
+!    ctwc(:,:,n) = icloud_stats(:,:,4,n)
+
+!**************************************************************************
+  else       ! use old definitions
+!**************************************************************************
+!   create a cloud and rainout/washout field, clouds occur where rh>80%
+!   total cloudheight is stored at level 0
+    write(*,*) 'Global fields: using cloud water from Parameterization'
+    do jy=0,nymin1
+      do ix=0,nxmin1
+! OLD METHOD
+        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
+          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
+          if (rh.gt.0.8) then ! in cloud
+            if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
+              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
+          else ! no cloud
+            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
+          endif
+        end do
+!END OLD METHOD
+      end do
+    end do
+  endif !readclouds
+
+
+     !********* 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
+!hg__________________________________
+!           rcw(ix,jy)=2E-7*prec**0.36
+!           rpc(ix,jy)=prec
+!hg end______________________________
+
+!      endif !hg read clouds
+
+
+
+
+!eso measure CPU time
+!  call mpif_mtime('verttransform',1)
+
+!eso print out the same measure as in Leo's routine
+    ! 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
+
