Index: src/verttransform.f90
===================================================================
--- src/verttransform.f90	(revision 0f20c31f55c80f501e55c69aad5a2dd4081628f1)
+++ src/verttransform.f90	(revision 4d456391fa86b9dc93fe32bac700d3c926f265cb)
@@ -21,5 +21,5 @@
 
 subroutine verttransform(n,uuh,vvh,wwh,pvh)
-!                      i  i   i   i   i
+!                         i  i   i   i   i
 !*****************************************************************************
 !                                                                            *
@@ -67,17 +67,22 @@
   use com_mod
   use cmapf_mod, only: cc2gll
+! eso TODO:
+! only used for timing of CPU measurement. Remove this (and calls to mpif_mtime below)
+! as this routine should not be dependent on MPI
+!  use mpi_mod
+! :TODO
 
   implicit none
 
   integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
-  integer :: rain_cloud_above,kz_inv
-! integer :: icloudtop !PS
-  real :: f_qvsat,pressure
-! real :: rh,lsp,convp
-  real :: rh,lsp,convp,prec,rhmin,precmin
-  real :: rhoh(nuvzmax),pinmconv(nzmax)
-  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
+  integer :: rain_cloud_above(0:nxmax-1,0:nymax-1),kz_inv,idx(0:nxmax-1,0:nymax-1)
+  real :: f_qvsat,pressure,cosf(0:nymax-1)
+  real :: rh,lsp,convp,tim,tim2,rhmin,precmin,prec
+  real :: uvzlev(0:nxmax-1,0:nymax-1,nuvzmax),rhoh(0:nxmax-1,0:nymax-1,nuvzmax)
+  real :: pinmconv(0:nxmax-1,0:nymax-1,nzmax)
+  real :: ew,pint(0:nxmax-1,0:nymax-1),tv(0:nxmax-1,0:nymax-1)
+  real :: tvold(0:nxmax-1,0:nymax-1),pold(0:nxmax-1,0:nymax-1),dz1,dz2,dz,ui,vi
   real :: xlon,ylat,xlonr,dzdx,dzdy
-  real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf
+  real :: dzdx1,dzdx2,dzdy1,dzdy2
   real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
   real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
@@ -85,7 +90,7 @@
   real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
   real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
-  real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
-  logical lconvectprec
+  real :: wzlev(0:nxmax-1,0:nymax-1,nwzmax)
   real,parameter :: const=r_air/ga
+!  integer:: ihr, imin, isec, i100th,ihr2, imin2, isec2, i100th2
   parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics
 
@@ -103,6 +108,14 @@
 ! 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.                                                *
+! 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
@@ -110,5 +123,5 @@
 ! 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
@@ -124,7 +137,7 @@
 
 
-    tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
+    tvold(ixm,jym)=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)
+    pold(ixm,jym)=ps(ixm,jym,1,n)
     height(1)=0.
 
@@ -133,13 +146,15 @@
       tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
 
-      if (abs(tv-tvold).gt.0.2) then
-        height(kz)=height(kz-1)+const*log(pold/pint)* &
-             (tv-tvold)/log(tv/tvold)
+      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/pint)*tv
+        height(kz)=height(kz-1)+ &
+             const*log(pold(ixm,jym)/pint(ixm,jym))*tv(ixm,jym)
       endif
 
-      tvold=tv
-      pold=pint
+      tvold(ixm,jym)=tv(ixm,jym)
+      pold(ixm,jym)=pint(ixm,jym)
     end do
 
@@ -167,11 +182,15 @@
 !*************************
 
+
   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)
-      uvwzlev(ix,jy,1)=0.
-      wzlev(1)=0.
-      rhoh(1)=pold/(r_air*tvold)
+      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)
 
 
@@ -179,102 +198,118 @@
 !******************************
 
-      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
-          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
-               (tv-tvold)/log(tv/tvold)
-        else
-          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
-        endif
-
-        tvold=tv
-        pold=pint
-      end do
-
-
-      do kz=2,nwz-1
-        wzlev(kz)=(uvwzlev(ix,jy,kz+1)+uvwzlev(ix,jy,kz))/2.
-      end do
-      wzlev(nwz)=wzlev(nwz-1)+uvwzlev(ix,jy,nuvz)-uvwzlev(ix,jy,nuvz-1)
+  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)=(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)))
+  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(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)
+
+  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 
-      clwc(ix,jy,1,n)=clwch(ix,jy,1,n)
-      ciwc(ix,jy,1,n)=ciwch(ix,jy,1,n)   
+  clwc(:,:,1,n)=clwch(:,:,1,n)
+  ciwc(:,:,1,n)=ciwch(:,:,1,n)   
 !hg 
-      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(:,:,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
-      clwc(ix,jy,nz,n)=clwch(ix,jy,nuvz,n)
-      ciwc(ix,jy,nz,n)=ciwch(ix,jy,nuvz,n)
+  clwc(:,:,nz,n)=clwch(:,:,nuvz,n)
+  ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n)
 !hg
-      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.uvwzlev(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)
+  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
-            clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
-            ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
+          clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
+          ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
 !hg
-            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.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
-            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)=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
-            clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2+clwch(ix,jy,kz,n)*dz1)/dz
-            ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz
+          clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2+clwch(ix,jy,kz,n)*dz1)/dz
+          ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz
 !hg
-
-            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
+          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
 
 
@@ -282,36 +317,45 @@
 !*************************
 
-      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
+  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
-        end do
-40      continue
-      end do
+        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(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
+  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
 
 
@@ -322,40 +366,43 @@
 
   do jy=1,ny-2
-    cosf=cos((real(jy)*dy+ylat0)*pi180)
-    do ix=1,nx-2
-
-      kmin=2
-      do iz=2,nz-1
-
-        ui=uu(ix,jy,iz,n)*dxconst/cosf
-        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
+    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
-        end do
-
-47      ix1=ix-1
+        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=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
-        dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
+        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=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
-        dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
+        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*ui+dzdy*vi)
+        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
@@ -363,5 +410,4 @@
     end do
   end do
-
 
 ! If north pole is in the domain, calculate wind velocities in polar
@@ -370,11 +416,12 @@
 
   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
+    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))
+               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
+               vvpol(ix,jy,iz,n))
         end do
       end do
@@ -390,7 +437,9 @@
       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)
+      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
+        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)/ &
@@ -408,5 +457,7 @@
       uuaux=-ffpol*sin(xlonr+ddpol)
       vvaux=-ffpol*cos(xlonr+ddpol)
-      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
+      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
+           vvpolaux)
+
       jy=nymin1
       do ix=0,nxmin1
@@ -441,11 +492,12 @@
 
   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
+    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))
+               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
+               vvpol(ix,jy,iz,n))
         end do
       end do
@@ -460,9 +512,12 @@
       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)
+      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
+        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
+        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
+             vv(nx/2-1,0,iz,n))+xlonr
       else
         ddpol=pi/2-xlonr
@@ -477,5 +532,6 @@
       uuaux=+ffpol*sin(xlonr-ddpol)
       vvaux=-ffpol*cos(xlonr-ddpol)
-      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
+      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
+           vvpolaux)
 
       jy=0
@@ -505,8 +561,4 @@
 
 
-!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
-
   if (readclouds) write(*,*) 'using cloud water from ECMWF'
 !  if (.not.readclouds) write(*,*) 'using cloud water from parameterization'
@@ -518,5 +570,5 @@
     do ix=0,nxmin1
 ! OLD METHOD
-      rain_cloud_above=0
+      rain_cloud_above(ix,jy)=0
       lsp=lsprec(ix,jy,1,n)
       convp=convprec(ix,jy,1,n)
@@ -529,5 +581,5 @@
         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
+            rain_cloud_above(ix,jy)=1
             cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
                  height(kz)-height(kz-1)
@@ -541,5 +593,5 @@
           endif
         else ! no cloud
-          if (rain_cloud_above.eq.1) then ! scavenging
+          if (rain_cloud_above(ix,jy).eq.1) then ! scavenging
             if (lsp.ge.convp) then
               clouds(ix,jy,kz,n)=5 ! lsp dominated washout
@@ -553,6 +605,6 @@
 
 ! PS 2012
-!          lsp=lsprec(ix,jy,1,n)
-!          convp=convprec(ix,jy,1,n)
+!      lsp=lsprec(ix,jy,1,n)
+!      convp=convprec(ix,jy,1,n)
 !          prec=lsp+convp
 !          if (lsp.gt.convp) then !  prectype='lsp'
@@ -560,5 +612,5 @@
 !          else ! prectype='cp '
 !            lconvectprec = .true.
-!          endif
+!           endif
 !HG METHOD
 !readclouds =.true.
@@ -572,23 +624,16 @@
 !        do kz=1, nz
 !         !clw & ciw are given in kg/kg  
-!         cloud_water=clwc(ix,jy,kz,n)+ciwc(ix,jy,kz,n)
-!         if (cloud_water .gt. 0) then
-!          cloud_min=min(nint(height(kz)),cloud_min) !hg needs reset each grid
-!          cloud_max=max(nint(height(kz)),cloud_max) !hg needs reset each grid
-!          cloud_col_wat=cloud_col_wat+cloud_water !hg needs reset each grid 
-!         endif
-!         cloud_ver=max(0,cloud_max-cloud_min)
-
-! if (clwc(ix,jy,kz,n).gt.0 .or.  ciwc(ix,jy,kz,n).gt.0) &
-!     !write(*,*) 'WATER',clwc(ix,jy,kz,n), 'ICE',ciwc(ix,jy,kz,n),'RH',rh,'KZ',kz &
-!     write(*,*) 'WATER',cloud_water,'RH',rh,'PREC',prec,'HEIGHT',height(kz) &
-!       enddo
-!       if ( cloud_min .ne. 99999 .and. cloud_max .ne. -1) write(*,*) 'CB', cloud_min, '   CT',cloud_max
-! if (prec .gt. 0) write(*,*) 'PREC',prec,'Cloud Bot',cloud_min,'Cloud Top',cloud_max, 'Cloud Vert. ext',cloud_ver &
-!                             ,'COLUMN cloud water',cloud_col_wat,'Conevctive' ,lconvectprec
-!       icloudbot(ix,jy,n)=cloud_min
-!       icloudthck(ix,jy,n)=cloud_ver
-!       rcw(ix,jy)=cloud_col_wat
-!       rpc(ix,jy)=prec
+          ! cloud_water=clwc(ix,jy,kz,n)+ciwc(ix,jy,kz,n)
+          ! if (cloud_water .gt. 0) then
+          !   cloud_min=min(nint(height(kz)),cloud_min) !hg needs reset each grid
+          !   cloud_max=max(nint(height(kz)),cloud_max) !hg needs reset each grid
+          !   cloud_col_wat=cloud_col_wat+cloud_water !hg needs reset each grid 
+          ! endif
+          ! cloud_ver=max(0,cloud_max-cloud_min)
+
+          ! icloudbot(ix,jy,n)=cloud_min
+          ! icloudthck(ix,jy,n)=cloud_ver
+          ! rcw(ix,jy)=cloud_col_wat
+          ! rpc(ix,jy)=prec
 !write(*,*) 'Using clouds from ECMWF' !hg END Henrik Code
 !END HG METHOD
@@ -613,5 +658,5 @@
 !              icloudtop=nint(height(kz)) ! use int to save memory
 !            endif
-!          enddo
+    ! 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
@@ -622,5 +667,5 @@
 !            rhmin = rhmin - 0.05
 !            if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum.
-!          endif
+!   end if
 
 !PS is based on looking at a limited set of comparison data
@@ -656,3 +701,39 @@
   end do
 
-end subroutine verttransform
+!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)
+
+!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
+
