verttransform.f90

Go to the documentation of this file.

!**********************************************************************
! 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
  !*****************************************************************************
  !  Changes Arnold, D. and Morton, D. (2015):                                 *
  !   -- description of local and common variables                             *
  !*****************************************************************************

  use par_mod
  use com_mod
  use cmapf_mod, only: cc2gll

  implicit none

  !***********************************************************************
  ! Subroutine Parameters:                                               *
  !    input                                                             *
  ! n                   temporal index for meteorological fields (1 to 3)*  
  ! uuh,vvh, wwh        wind components in ECMWF model levels            *  
  ! pvh                 potential vorticity in ECMWF model levels        *
  integer :: n
  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
  real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
  !***********************************************************************

  !***********************************************************************
  ! Local variables                                                      *
  !                                                                      *
  !  ew                   subroutine/function to calculate saturation    *  
  !                         water vaport for a given air temperature     *
  ! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition     *
  ! rain_cloud_above [0,1]  whether there is a raining cloud or not      *  
  ! kz_inv  inverted indez for kz                                        *
  ! f_qvsat  Saturation water vapor specific humidity (kg/kg)            *
  ! pressure                                                             *
  ! rh       relative humidity                                           *    
  ! lsp      large scale precip at one grid cell                         *
  ! convp    convective precip at one grid cell                          *
  ! uvzlev    height of the eta half-levels                              *
  ! rhoh     density in model levels                                     *
  ! pinmconv conversion factor                                           *
  ! pint     pressure on model levels (using akz, bkz, ps)               *
  ! tv       virtual temperature                                         *
  ! tvold,pold      dummy variables to keep previous value               *
  ! dz1, dz2, dz  differences heights model levels, pressure levels      *    
  ! ui, vi                                                               *   
  ! xlon,xlat                                                            *
  ! xlonr         xlon*pi/180.                                           * 
  ! dzdx,dzdy     slope corrections                                      *
  ! dzdx1,dzdx2,dzdy1   slope corrections in each direction              *
  ! uuaux,vvaux,uupolaux,vvpolaux  auxiliary variables for polar sterero.*
  ! ddpol,ffpol  for special treatment of poles                          *
  ! wdummy                                                               *
  ! wzlev                                                                *
  ! uvwzlev                                                              *

  integer :: rain_cloud_above,kz_inv
  real :: f_qvsat,pressure
  real :: rh,lsp,convp
  real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
  real :: xlon,ylat,xlonr,dzdx,dzdy
  real :: dzdx1,dzdx2,dzdy1,dzdy2
  real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
  real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)

  ! Other variables:
  ! ix,jy,kz,kzmin     loop control indices in each direction            *
  ! kl,klp
  ! ix1,jy1,ixp,jyp,ixm,jym 
  integer :: ix,jy,kz,iz,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
  !***********************************************************************

  !***********************************************************************
  ! Local constants                                                      *
  real,parameter :: const=r_air/ga
  !***********************************************************************

  !***********************************************************************
  ! Global variables                                                     *
  !     from com_mod                                                     *
  ! nxmin1,nymin1           nx-1, ny-1, respectively                     *
  ! tt2                  2 meter temperature                             *
  !  ps                   surface pressure                               * 
  ! td2                  2 meter dew point                               *
  ! height                  heights of all levels                        *
  ! akm, bkm  coeffs.  which regulate vertical discretization of ecmwf   *
  ! akz, bkz  model discretization coeffizients at the centre of layers  *
  ! nuvz,nwz                vertical dimension of original ECMWF data    *
  ! nx,ny,nz                actual dimensions of wind fields in x,y and z*
  ! aknew,bknew model discretization coeffs. at the interpolated levels  *
  ! uu,vv,ww [m/2]       wind components in x,y and z direction          *
  ! qv                   specific humidity data                          *
  ! pv (pvu)             potential vorticity                             *
  ! rho [kg/m3]          air density                                     *
  ! drhodz [kg/m2]       vertical air density gradient                   *
  ! dxconst,dyconst         auxiliary variables for utransform,vtransform*
  ! ylat0                 geographical latitude of lower left grid point *
  ! xlon0                 geographical longitude of lower left grid point*
  ! uupol,vvpol [m/s]   wind components in polar stereographic projection*
  ! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition     *    
  ! lsprec [mm/h]        large scale total precipitation                 *
  ! convprec [mm/h]      convective precipitation                        *
  ! cloudsh                                                              *
  !                                                                      *
  !***********************************************************************

!-----------------------------------------------------------------------------



  logical :: init = .true.


  !*************************************************************************
  ! 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.                                                  *
  !*************************************************************************


  !      do 897 kz=1,nuvz
  !       write (*,*) 'akz: ',akz(kz),'bkz',bkz(kz)
  !897     continue

  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=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
         ps(ixm,jym,1,n))
    pold=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))


  ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
  ! upon the transformation to z levels. In order to save computer memory, this is
  ! not done anymore in the standard version. However, this option can still be
  ! switched on by replacing the following lines with those below, that are
  ! currently commented out.
  ! Note that two more changes are necessary in this subroutine below.
  ! One change is also necessary in gridcheck.f, and another one in verttransform_nests.
  !*****************************************************************************

      if (abs(tv-tvold).gt.0.2) then
        height(kz)= &
             height(kz-1)+const*log(pold/pint)* &
             (tv-tvold)/log(tv/tvold)
      else
        height(kz)=height(kz-1)+ &
             const*log(pold/pint)*tv
      endif

  ! Switch on following lines to use doubled vertical resolution
  !*************************************************************
  !    if (abs(tv-tvold).gt.0.2) then
  !      height((kz-1)*2)=
  !    +      height(max((kz-2)*2,1))+const*log(pold/pint)*
  !    +      (tv-tvold)/log(tv/tvold)
  !    else
  !      height((kz-1)*2)=height(max((kz-2)*2,1))+
  !    +      const*log(pold/pint)*tv
  !    endif
  ! End doubled vertical resolution

      tvold=tv
      pold=pint
    end do


  ! Switch on following lines to use doubled vertical resolution
  !*************************************************************
  !  do 7 kz=3,nz-1,2
  !    height(kz)=0.5*(height(kz-1)+height(kz+1))
  !  height(nz)=height(nz-1)+height(nz-1)-height(nz-2)
  ! End doubled vertical resolution


  ! 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.

  endif


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

  do jy=0,nymin1
    do ix=0,nxmin1
      tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
           ps(ix,jy,1,n))
      pold=ps(ix,jy,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(ix,jy,1,n)
        tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n))
        rhoh(kz)=pint/(r_air*tv)

        if (abs(tv-tvold).gt.0.2) then
          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
               (tv-tvold)/log(tv/tvold)
        else
          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
        endif

        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)

      uvwzlev(ix,jy,1)=0.0
      do kz=2,nuvz
        uvwzlev(ix,jy,kz)=uvzlev(kz)
      end do

  ! Switch on following lines to use doubled vertical resolution
  ! Switch off the three lines above.
  !*************************************************************
  !22          uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
  !     do 23 kz=2,nwz
  !23          uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
  ! End doubled vertical resolution

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

      pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
           ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
           (aknew(1)+bknew(1)*ps(ix,jy,1,n)))
      do kz=2,nz-1
        pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
             ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
             (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
      end do
      pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
           ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
           (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))

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

      uu(ix,jy,1,n)=uuh(ix,jy,1)
      vv(ix,jy,1,n)=vvh(ix,jy,1)
      tt(ix,jy,1,n)=tth(ix,jy,1,n)
      qv(ix,jy,1,n)=qvh(ix,jy,1,n)
      pv(ix,jy,1,n)=pvh(ix,jy,1)
      rho(ix,jy,1,n)=rhoh(1)
      uu(ix,jy,nz,n)=uuh(ix,jy,nuvz)
      vv(ix,jy,nz,n)=vvh(ix,jy,nuvz)
      tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
      qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
      pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
      rho(ix,jy,nz,n)=rhoh(nuvz)
      kmin=2
      do iz=2,nz-1
        do kz=kmin,nuvz
          if(height(iz).gt.uvzlev(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)
            pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
            rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
            goto 30
          endif
          if ((height(iz).gt.uvzlev(kz-1)).and. &
               (height(iz).le.uvzlev(kz))) then
           dz1=height(iz)-uvzlev(kz-1)
           dz2=uvzlev(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
           pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
           rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
           kmin=kz
           goto 30
          endif
        end do
30      continue
      end do


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

      ww(ix,jy,1,n)=wwh(ix,jy,1)*pinmconv(1)
      ww(ix,jy,nz,n)=wwh(ix,jy,nwz)*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
           ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 &
                +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz
           kmin=kz
           goto 40
          endif
        end do
40      continue
      end do

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

      drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
           (height(2)-height(1))
      do kz=2,nz-1
        drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
             (height(kz+1)-height(kz-1))
      end do
      drhodz(ix,jy,nz,n)=drhodz(ix,jy,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
    do ix=1,nx-2

      kmin=2
      do iz=2,nz-1

        ui=uu(ix,jy,iz,n)*dxconst/cos((real(jy)*dy+ylat0)*pi180)
        vi=vv(ix,jy,iz,n)*dyconst

        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

        ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdx*ui+dzdy*vi)

      end do

    end do
  end do


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

  if (nglobal) then
    do jy=int(switchnorthg)-2,nymin1
      ylat=ylat0+real(jy)*dy
      do ix=0,nxmin1
        xlon=xlon0+real(ix)*dx
        do iz=1,nz
          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 jy=0,int(switchsouthg)+3
      ylat=ylat0+real(jy)*dy
      do ix=0,nxmin1
        xlon=xlon0+real(ix)*dx
        do iz=1,nz
          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


  !write (*,*) 'initializing clouds, n:',n,nymin1,nxmin1,nz
  !   create a cloud and rainout/washout field, clouds occur where rh>80%
  !   total cloudheight is stored at level 0
  do jy=0,nymin1
    do ix=0,nxmin1
      rain_cloud_above=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=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.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 do
  end do

  !do 102 kz=1,nuvz
  !write(an,'(i02)') kz+10
  !write(*,*) nuvz,nymin1,nxmin1,'--',an,'--'
  !open(4,file='/nilu_wrk2/sec/cloudtest/cloud'//an,form='formatted')
  !do 101 jy=0,nymin1
  !    write(4,*) (clouds(ix,jy,kz,n),ix=1,nxmin1)
  !101   continue
  ! close(4)
  !102   continue

  ! open(4,file='/nilu_wrk2/sec/cloudtest/height',form='formatted')
  ! do 103 jy=0,nymin1
  !     write (4,*)
  !+       (height(kz),kz=1,nuvz)
  !103   continue
  ! close(4)

  !open(4,file='/nilu_wrk2/sec/cloudtest/p',form='formatted')
  ! do 104 jy=0,nymin1
  !     write (4,*)
  !+       (r_air*tt(ix,jy,1,n)*rho(ix,jy,1,n),ix=1,nxmin1)
  !104   continue
  ! close(4)
end subroutine verttransform_ecmwf