Index: src/FLEXPART.f90
===================================================================
--- src/FLEXPART.f90	(revision b0434e1968915cd3d4a8afcbf2cc275acecb9ff7)
+++ src/FLEXPART.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -315,9 +315,4 @@
   endif
 
-  !! testing !!
-  ! open(999,file=trim(path(1))//'OH_FIELDS/jscalar_50N.txt',action='write',status='new')
-  ! open(998,file=trim(path(1))//'OH_FIELDS/jscalar_50S.txt',action='write',status='new')
-
-
   ! Write basic information on the simulation to a file "header"
   ! and open files that are to be kept open throughout the simulation
Index: src/FLEXPART_MPI.f90
===================================================================
--- src/FLEXPART_MPI.f90	(revision b0434e1968915cd3d4a8afcbf2cc275acecb9ff7)
+++ src/FLEXPART_MPI.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -348,9 +348,4 @@
   endif
 
-  !! testing !!
-  ! open(999,file=trim(path(1))//'OH_FIELDS/jscalar_50N.txt',action='write',status='new')
-  ! open(998,file=trim(path(1))//'OH_FIELDS/jscalar_50S.txt',action='write',status='new')
-
-
   ! Write basic information on the simulation to a file "header"
   ! and open files that are to be kept open throughout the simulation
Index: src/calcpar_nests.f90
===================================================================
--- src/calcpar_nests.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/calcpar_nests.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -103,4 +103,5 @@
       ustarn(ix,jy,1,n,l)=scalev(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), &
            td2n(ix,jy,1,n,l),surfstrn(ix,jy,1,n,l))
+      if (ustarn(ix,jy,1,n,l).le.1.e-8) ustarn(ix,jy,1,n,l)=1.e-8
 
   ! 2) Calculation of inverse Obukhov length scale
@@ -135,5 +136,5 @@
         subsceff=min(excessoron(ix,jy,l),hmixplus)
       else
-        subsceff=0
+        subsceff=0.0
       endif
   !
@@ -149,6 +150,7 @@
 
       if (DRYDEP) then
-        z0(4)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga
-        z0(9)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga
+        ! z0(4)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga
+        ! z0(9)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga
+        z0(7)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga
 
   ! Calculate relative humidity at surface
@@ -228,6 +230,8 @@
   end do
 
-
-    call calcpv_nests(l,n,uuhn,vvhn,pvhn)
+  ! Calculation of potential vorticity on 3-d grid
+  !***********************************************
+
+  call calcpv_nests(l,n,uuhn,vvhn,pvhn)
 
   end do
Index: src/calcpv.f90
===================================================================
--- src/calcpv.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/calcpv.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -70,8 +70,6 @@
   enddo
 
-! :DBG: halts with SIGFPE if compiling with -ffpe-trap
-!  where(ppml.le.0) ppml=10000.0
 !  ppmk(:,:,1:nuvz)=(100000./ppml(:,:,1:nuvz))**kappa
-  ppmk=(100000./ppml)**kappa
+  ppmk(0:nxmin1,0:nymin1,1:nuvz)=(100000./ppml(0:nxmin1,0:nymin1,1:nuvz))**kappa
 
   do jy=0,nymin1
Index: src/calcpv_nests.f90
===================================================================
--- src/calcpv_nests.f90	(revision 4fbe7a5609fe43bdc3940768d8f4d949237cc44a)
+++ src/calcpv_nests.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -68,5 +68,6 @@
     enddo
   enddo
-  ppmk=(100000./ppml)**kappa
+!  ppmk=(100000./ppml)**kappa
+  ppmk(0:nxn(l)-1,0:nyn(l)-1,1:nuvz)=(100000./ppml(0:nxn(l)-1,0:nyn(l)-1,1:nuvz))**kappa
 
   do jy=0,nyn(l)-1
Index: src/com_mod.f90
===================================================================
--- src/com_mod.f90	(revision b0434e1968915cd3d4a8afcbf2cc275acecb9ff7)
+++ src/com_mod.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -134,5 +134,7 @@
   logical :: readclouds
 !ESO DEC 2015 whether or not both clwc and ciwc are present (if so they are summed)
-  logical :: sumclouds 
+  logical :: sumclouds
+
+  logical,dimension(maxnests) :: readclouds_nest, sumclouds_nest
   
 
@@ -370,5 +372,4 @@
 !ZHG Sep 2015  
    real :: icloud_stats(0:nxmax-1,0:nymax-1,5,numwfmem)
-   real :: clh(0:nxmax-1,0:nymax-1,nuvzmax,numwfmem)
    real :: clw4(0:nxmax-1,0:nymax-1,numwfmem) ! eso: =icloud_stats(:,:,4,:)
 
@@ -487,6 +488,7 @@
 
   real,allocatable,dimension(:,:,:,:,:) :: uun, vvn, wwn, ttn, qvn, pvn,&
-       & rhon, drhodzn, tthn, qvhn
-  integer,allocatable,dimension(:,:,:,:) :: cloudsnh
+       & rhon, drhodzn, tthn, qvhn, clwcn, ciwcn, clwn, clwchn, ciwchn
+  real,allocatable,dimension(:,:,:,:) :: clw4n
+  integer,allocatable,dimension(:,:,:,:) :: cloudshn
   integer(kind=1),allocatable,dimension(:,:,:,:,:) :: cloudsn
 
@@ -788,11 +790,27 @@
     allocate(qvn(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests))
     allocate(pvn(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests))
-    allocate(cloudsn(0:nxmaxn-1,0:nymaxn-1,0:nzmax,numwfmem,numbnests))
-    allocate(cloudsnh(0:nxmaxn-1,0:nymaxn-1,numwfmem,numbnests))
+    allocate(clwcn(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests))
+    allocate(ciwcn(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests))
+    allocate(clwn(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests))
+
+    allocate(cloudsn(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests))
+    allocate(cloudshn(0:nxmaxn-1,0:nymaxn-1,numwfmem,numbnests))
     allocate(rhon(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests))
     allocate(drhodzn(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests))
     allocate(tthn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,numwfmem,numbnests))
     allocate(qvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,numwfmem,numbnests))
-
+    allocate(clwchn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,numwfmem,numbnests))
+    allocate(ciwchn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,numwfmem,numbnests))
+    allocate(clw4n(0:nxmax-1,0:nymax-1,numwfmem,numbnests))
+
+!    clw4n(:,:,:,:)=0.
+    clwcn(:,:,:,:,:)=0.
+    ciwcn(:,:,:,:,:)=0.
+    clwchn(:,:,:,:,:)=0.
+    ciwchn(:,:,:,:,:)=0.
+    ! clwn(:,:,:,:,:)=0.
+
+    ! cloudsn(:,:,:,:,:)=0
+    ! cloudshn(:,:,:,:)=0
     
   end subroutine com_mod_allocate_nests
Index: src/concoutput.f90
===================================================================
--- src/concoutput.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
+++ src/concoutput.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -626,22 +626,24 @@
 !*************************
 
-  do ks=1,nspec
-    do kp=1,maxpointspec_act
-      do i=1,numreceptor
-        creceptor(i,ks)=0.
-      end do
-      do jy=0,numygrid-1
-        do ix=0,numxgrid-1
-          do l=1,nclassunc
-            do nage=1,nageclass
-              do kz=1,numzgrid
-                gridunc(ix,jy,kz,ks,kp,l,nage)=0.
-              end do
-            end do
-          end do
-        end do
-      end do
-    end do
-  end do
+  ! do ks=1,nspec
+  !   do kp=1,maxpointspec_act
+  !     do i=1,numreceptor
+  !       creceptor(i,ks)=0.
+  !     end do
+  !     do jy=0,numygrid-1
+  !       do ix=0,numxgrid-1
+  !         do l=1,nclassunc
+  !           do nage=1,nageclass
+  !             do kz=1,numzgrid
+  !               gridunc(ix,jy,kz,ks,kp,l,nage)=0.
+  !             end do
+  !           end do
+  !         end do
+  !       end do
+  !     end do
+  !   end do
+  ! end do
+  creceptor(:,:)=0.
+  gridunc(:,:,:,:,:,:,:)=0.
 
 
Index: src/concoutput_nest.f90
===================================================================
--- src/concoutput_nest.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
+++ src/concoutput_nest.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -545,22 +545,24 @@
   !*************************
 
-  do ks=1,nspec
-  do kp=1,maxpointspec_act
-    do i=1,numreceptor
-      creceptor(i,ks)=0.
-    end do
-    do jy=0,numygridn-1
-      do ix=0,numxgridn-1
-        do l=1,nclassunc
-          do nage=1,nageclass
-            do kz=1,numzgrid
-              griduncn(ix,jy,kz,ks,kp,l,nage)=0.
-            end do
-          end do
-        end do
-      end do
-    end do
-  end do
-  end do
+  ! do ks=1,nspec
+  ! do kp=1,maxpointspec_act
+  !   do i=1,numreceptor
+  !     creceptor(i,ks)=0.
+  !   end do
+  !   do jy=0,numygridn-1
+  !     do ix=0,numxgridn-1
+  !       do l=1,nclassunc
+  !         do nage=1,nageclass
+  !           do kz=1,numzgrid
+  !             griduncn(ix,jy,kz,ks,kp,l,nage)=0.
+  !           end do
+  !         end do
+  !       end do
+  !     end do
+  !   end do
+  ! end do
+  ! end do
+  creceptor(:,:)=0.
+  griduncn(:,:,:,:,:,:,:)=0.
 
 
Index: src/ecmwf_mod.f90
===================================================================
--- src/ecmwf_mod.f90	(revision b0434e1968915cd3d4a8afcbf2cc275acecb9ff7)
+++ src/ecmwf_mod.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -41,6 +41,9 @@
   !*********************************************
 
-  integer,parameter :: nxmax=361,nymax=181,nuvzmax=152,nwzmax=152,nzmax=152 !ECMWF new 
+!  integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF new 
+  integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new 
   integer,parameter :: nxshift=359
+!  integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !test BUG 
+!  integer,parameter :: nxshift=0
 
 ! 
@@ -49,5 +52,6 @@
   !*********************************************
 
-  integer,parameter :: maxnests=1,nxmaxn=361,nymaxn=181 !ECMWF
+  integer,parameter :: maxnests=1,nxmaxn=361,nymaxn=181
+!  integer,parameter :: maxnests=1,nxmaxn=86,nymaxn=31
 
 
Index: src/gridcheck.f90
===================================================================
--- src/gridcheck.f90	(revision b0434e1968915cd3d4a8afcbf2cc275acecb9ff7)
+++ src/gridcheck.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -194,6 +194,4 @@
   elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
     isec1(6)=246         ! indicatorOfParameter
-    ! readclouds=.true.
-    ! sumclouds=.false.
   elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
     isec1(6)=247         ! indicatorOfParameter
@@ -202,6 +200,4 @@
   elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc
     isec1(6)=201031      ! indicatorOfParameter
-    ! readclouds=.true.
-    ! sumclouds=.true.
   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
     isec1(6)=134         ! indicatorOfParameter
@@ -366,4 +362,20 @@
   endif ! gotGrid
 
+  if (nx.gt.nxmax) then
+   write(*,*) 'FLEXPART error: Too many grid points in x direction.'
+    write(*,*) 'Reduce resolution of wind fields.'
+    write(*,*) 'Or change parameter settings in file par_mod.'
+    write(*,*) nx,nxmax
+    stop
+  endif
+
+  if (ny.gt.nymax) then
+   write(*,*) 'FLEXPART error: Too many grid points in y direction.'
+    write(*,*) 'Reduce resolution of wind fields.'
+    write(*,*) 'Or change parameter settings in file par_mod.'
+    write(*,*) ny,nymax
+    stop
+  endif
+
   k=isec1(8)
   if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
@@ -410,20 +422,4 @@
   nwz =iwmax
   if(nuvz.eq.nlev_ec) nwz=nlev_ec+1
-
-  if (nx.gt.nxmax) then
-   write(*,*) 'FLEXPART error: Too many grid points in x direction.'
-    write(*,*) 'Reduce resolution of wind fields.'
-    write(*,*) 'Or change parameter settings in file par_mod.'
-    write(*,*) nx,nxmax
-    stop
-  endif
-
-  if (ny.gt.nymax) then
-   write(*,*) 'FLEXPART error: Too many grid points in y direction.'
-    write(*,*) 'Reduce resolution of wind fields.'
-    write(*,*) 'Or change parameter settings in file par_mod.'
-    write(*,*) ny,nymax
-    stop
-  endif
 
   if (nuvz+1.gt.nuvzmax) then
Index: src/gridcheck_nests.f90
===================================================================
--- src/gridcheck_nests.f90	(revision b0434e1968915cd3d4a8afcbf2cc275acecb9ff7)
+++ src/gridcheck_nests.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -172,6 +172,4 @@
   elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
     isec1(6)=246         ! indicatorOfParameter
-    ! readclouds=.true.
-    ! sumclouds=.false.
   elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
     isec1(6)=247         ! indicatorOfParameter
@@ -180,12 +178,10 @@
   elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc
     isec1(6)=201031      ! indicatorOfParameter
-    ! readclouds=.true.
-    ! sumclouds=.true.
   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
     isec1(6)=134         ! indicatorOfParameter
   elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
     isec1(6)=135         ! indicatorOfParameter
-  elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot !added bymc to make it consistent with new gridchek.f90
-    isec1(6)=135         ! indicatorOfParameter    !                     ! " " " " " " " " " " " " " " " " " " " " " " " "  " " " 
+  elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot !added bymc to make it consistent with new gridcheck.f90
+    isec1(6)=135         ! indicatorOfParameter    !
   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
     isec1(6)=151         ! indicatorOfParameter
@@ -258,4 +254,22 @@
   endif ! ifield
 
+  if (nxn(l).gt.nxmaxn) then
+  write(*,*) 'FLEXPART error: Too many grid points in x direction.'
+  write(*,*) 'Reduce resolution of wind fields (file GRIDSPEC)'
+  write(*,*) 'for nesting level ',l
+  write(*,*) 'Or change parameter settings in file par_mod.'
+  write(*,*) nxn(l),nxmaxn
+  stop
+  endif
+
+  if (nyn(l).gt.nymaxn) then
+  write(*,*) 'FLEXPART error: Too many grid points in y direction.'
+  write(*,*) 'Reduce resolution of wind fields (file GRIDSPEC)'
+  write(*,*) 'for nesting level ',l
+  write(*,*) 'Or change parameter settings in file par_mod.'
+  write(*,*) nyn(l),nymaxn
+  stop
+  endif
+
   !HSO  get the second part of the grid dimensions only from GRiB1 messages
   if (isec1(6) .eq. 167 .and. (gotGrib.eq.0)) then !added by mc to make it consistent with new gridchek.f90 note that gotGrid must be changed in gotGrib!!
@@ -333,22 +347,4 @@
   if(nuvzn.eq.nlev_ec) nwzn=nlev_ecn+1
 
-  if (nxn(l).gt.nxmaxn) then
-  write(*,*) 'FLEXPART error: Too many grid points in x direction.'
-  write(*,*) 'Reduce resolution of wind fields (file GRIDSPEC)'
-  write(*,*) 'for nesting level ',l
-  write(*,*) 'Or change parameter settings in file par_mod.'
-  write(*,*) nxn(l),nxmaxn
-  stop
-  endif
-
-  if (nyn(l).gt.nymaxn) then
-  write(*,*) 'FLEXPART error: Too many grid points in y direction.'
-  write(*,*) 'Reduce resolution of wind fields (file GRIDSPEC)'
-  write(*,*) 'for nesting level ',l
-  write(*,*) 'Or change parameter settings in file par_mod.'
-  write(*,*) nyn(l),nymaxn
-  stop
-  endif
-
   if ((nuvzn.gt.nuvzmax).or.(nwzn.gt.nwzmax)) then
   write(*,*) 'FLEXPART error: Nested wind fields have too many'// &
Index: src/initialize.f90
===================================================================
--- src/initialize.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/initialize.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -99,5 +99,5 @@
   jyp=jy+1
 
-  h=max(hmix(ix ,jy ,1,memind(1)), &
+  h=max(hmix(ix ,jy,1,memind(1)), &
        hmix(ixp,jy ,1,memind(1)), &
        hmix(ix ,jyp,1,memind(1)), &
@@ -157,13 +157,13 @@
     wp=rannumb(nrand+2)
     if (.not.turbswitch) then     ! modified by mc
+      wp=wp*sigw
+    else if (cblflag.eq.1) then   ! modified by mc
+      if(-h/ol.gt.5) then
+!if (ol.lt.0.) then
+!if (ol.gt.0.) then !by mc : only for test correct is lt.0
+        call initialize_cbl_vel(idummy,zt,ust,wst,h,sigw,wp,ol)
+      else
         wp=wp*sigw
-    else if (cblflag.eq.1) then   ! modified by mc
-        if(-h/ol.gt.5) then
-        !if (ol.lt.0.) then
-        !if (ol.gt.0.) then !by mc : only for test correct is lt.0
-            call initialize_cbl_vel(idummy,zt,ust,wst,h,sigw,wp,ol)
-        else
-            wp=wp*sigw
-        end if            
+      end if
     end if
 
Index: src/interpol_rain_nests.f90
===================================================================
--- src/interpol_rain_nests.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/interpol_rain_nests.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -88,8 +88,12 @@
   !***********************************************************
 
-  if (xt.ge.(real(nxn(ngrid)-1)-0.0001)) &
-       xt=real(nxn(ngrid)-1)-0.0001
-  if (yt.ge.(real(nyn(ngrid)-1)-0.0001)) &
-       yt=real(nyn(ngrid)-1)-0.0001
+  ! if (xt.ge.(real(nxn(ngrid)-1)-0.0001)) &
+  !      xt=real(nxn(ngrid)-1)-0.0001
+  ! if (yt.ge.(real(nyn(ngrid)-1)-0.0001)) &
+  !      yt=real(nyn(ngrid)-1)-0.0001
+
+! ESO make it consistent with interpol_rain
+  if (xt.ge.(real(nxn(ngrid)-1))) xt=real(nxn(ngrid)-1)-0.00001
+  if (yt.ge.(real(nyn(ngrid)-1))) yt=real(nyn(ngrid)-1)-0.00001
 
 
@@ -105,4 +109,5 @@
   ix=int(xt)
   jy=int(yt)
+
   ixp=ix+1
   jyp=jy+1
Index: src/makefile
===================================================================
--- src/makefile	(revision b0434e1968915cd3d4a8afcbf2cc275acecb9ff7)
+++ src/makefile	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -388,5 +388,5 @@
 readdepo.o: com_mod.o par_mod.o
 readlanduse.o: com_mod.o par_mod.o
-readlanduse_int1.o: com_mod.o par_mod.o
+#readlanduse_int1.o: com_mod.o par_mod.o
 readOHfield.o: com_mod.o oh_mod.o par_mod.o
 readoutgrid.o: com_mod.o outg_mod.o par_mod.o
Index: src/mpi_mod.f90
===================================================================
--- src/mpi_mod.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
+++ src/mpi_mod.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -860,6 +860,4 @@
 !   step
 !
-! TODO
-!   GFS version
 !
 !*******************************************************************************
@@ -922,5 +920,4 @@
 
 ! The non-reader processes need to know if cloud water were read.
-! TODO: only at first step or always?
     call MPI_Bcast(readclouds,1,MPI_LOGICAL,id_read,MPI_COMM_WORLD,mp_ierr)
     if (mp_ierr /= 0) goto 600 
@@ -1017,7 +1014,4 @@
     if (mp_ierr /= 0) goto 600 
 
-    call MPI_Bcast(z0,numclass,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
-    if (mp_ierr /= 0) goto 600 
-
     if (mp_measure_time) call mpif_mtime('commtime',1)
 
@@ -1045,6 +1039,4 @@
 !   step
 !
-! TODO
-!   Transfer cloud water/ice if and when available for nested
 !
 !***********************************************************************
@@ -1060,9 +1052,9 @@
 
 ! Common array sizes used for communications
-    integer :: d3_size1 = nxmaxn*nymaxn*nzmax*maxnests
-    integer :: d3_size2 = nxmaxn*nymaxn*nuvzmax*maxnests
-    integer :: d2_size1 = nxmaxn*nymaxn*maxnests
-    integer :: d2_size2 = nxmaxn*nymaxn*maxspec*maxnests
-    integer :: d2_size3 = nxmaxn*nymaxn*maxnests
+    integer :: d3_size1 = nxmaxn*nymaxn*nzmax
+    integer :: d3_size2 = nxmaxn*nymaxn*nuvzmax
+    integer :: d2_size1 = nxmaxn*nymaxn
+    integer :: d2_size2 = nxmaxn*nymaxn*maxspec
+    integer :: d2_size3 = nxmaxn*nymaxn
 
     integer :: d3s1,d3s2,d2s1,d2s2
@@ -1106,4 +1098,8 @@
 !**********************************************************************
 
+! The non-reader processes need to know if cloud water were read.
+    call MPI_Bcast(readclouds_nest,maxnests,MPI_LOGICAL,id_read,MPI_COMM_WORLD,mp_ierr)
+    if (mp_ierr /= 0) goto 600 
+
 ! Static fields/variables sent only at startup
     if (first_call) then
@@ -1120,6 +1116,6 @@
 
 ! MPI prefers contiguous arrays for sending (else a buffer is created),
-! hence the loop
-
+! hence the loop over nests
+!**********************************************************************
     do i=1, numbnests 
 ! 3D fields
@@ -1145,8 +1141,16 @@
       if (mp_ierr /= 0) goto 600 
       call MPI_Bcast(cloudsn(:,:,:,li:ui,i),d3s1,MPI_INTEGER1,id_read,MPI_COMM_WORLD,mp_ierr)
-      if (mp_ierr /= 0) goto 600 
+      if (mp_ierr /= 0) goto 600
+
+! cloud water/ice:
+    if (readclouds_nest(i)) then
+      ! call MPI_Bcast(icloud_stats(:,:,:,li:ui),d2s1*5,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
+      ! if (mp_ierr /= 0) goto 600
+      call MPI_Bcast(clw4n(:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
+      if (mp_ierr /= 0) goto 600
+    end if
 
 ! 2D fields
-      call MPI_Bcast(cloudsnh(:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
+      call MPI_Bcast(cloudshn(:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
       if (mp_ierr /= 0) goto 600
       call MPI_Bcast(vdepn(:,:,:,li:ui,i),d2s2,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
@@ -1206,5 +1210,4 @@
 !   mind    -- index where to place new fields
 !
-! TODO
 !
 !*******************************************************************************
@@ -1392,5 +1395,4 @@
 !   memstat -- input, used to resolve windfield index being received
 !
-! TODO
 !
 !*******************************************************************************
@@ -2019,5 +2021,4 @@
 ! In the implementation with 3 fields, the processes may have posted
 ! MPI_Irecv requests that should be cancelled here
-!! TODO:
 ! if (.not.lmp_sync) then
 !   r=mp_pid*nvar_async
@@ -2105,5 +2106,4 @@
 !   
 !
-! TODO
 !
 !*******************************************************************************
@@ -2132,5 +2132,5 @@
     uu(:,:,:,li:ui) = 10.0
     vv(:,:,:,li:ui) = 0.0
-    uupol(:,:,:,li:ui) = 10.0 ! TODO check if ok
+    uupol(:,:,:,li:ui) = 10.0
     vvpol(:,:,:,li:ui)=0.0
     ww(:,:,:,li:ui)=0.
@@ -2164,5 +2164,4 @@
     tropopause(:,:,:,li:ui)=10000.
     oli(:,:,:,li:ui)=0.01
-    z0=1.0
  
   end subroutine set_fields_synthetic
Index: src/netcdf_output_mod.f90
===================================================================
--- src/netcdf_output_mod.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
+++ src/netcdf_output_mod.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -79,4 +79,9 @@
   implicit none
 
+  private
+
+  public :: writeheader_netcdf,concoutput_surf_nest_netcdf,concoutput_netcdf,&
+       &concoutput_nest_netcdf,concoutput_surf_netcdf
+
 !  include 'netcdf.inc'
 
@@ -97,5 +102,5 @@
   real,parameter :: eps=nxmax/3.e5
 
-  private:: writemetadata, output_units, nf90_err
+!  private:: writemetadata, output_units, nf90_err
 
   ! switch output of release point information on/off
@@ -628,8 +633,16 @@
 
   ! volume
-  call nf90_err(nf90_put_var(ncid, volID, volume(:,:,:)))
+  if (lnest) then
+    call nf90_err(nf90_put_var(ncid, volID, volumen(:,:,:)))
+  else
+    call nf90_err(nf90_put_var(ncid, volID, volume(:,:,:)))
+  end if
 
   ! area
-  call nf90_err(nf90_put_var(ncid, areaID, area(:,:)))
+  if (lnest) then
+    call nf90_err(nf90_put_var(ncid, areaID, arean(:,:)))
+  else
+    call nf90_err(nf90_put_var(ncid, areaID, area(:,:)))
+  end if
 
   if (write_releases.eqv..true.) then
Index: src/par_mod.f90
===================================================================
--- src/par_mod.f90	(revision b0434e1968915cd3d4a8afcbf2cc275acecb9ff7)
+++ src/par_mod.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -139,7 +139,5 @@
   !integer,parameter :: nxmax=361,nymax=181,nuvzmax=152,nwzmax=152,nzmax=152 !ECMWF new 
   !integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF
-  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=26,nwzmax=26,nzmax=26
   !integer,parameter :: nxmax=721,nymax=361,nuvzmax=64,nwzmax=64,nzmax=64
-  !integer,parameter :: nxmax=1201,nymax=235,nuvzmax=58,nwzmax=58,nzmax=58
 
 !  integer,parameter :: nxshift=359 ! for ECMWF
@@ -216,7 +214,7 @@
   !**************************************************
 
-  integer,parameter :: maxpart=40000000
+  integer,parameter :: maxpart=400000
   integer,parameter :: maxspec=6
-  integer,parameter :: minmass=0.0 !0.0001
+  real,parameter :: minmass=0.0 !0.0001
 
   ! maxpart                 Maximum number of particles
Index: src/readavailable.f90
===================================================================
--- src/readavailable.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
+++ src/readavailable.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -283,5 +283,5 @@
   stop
 
-999   write(*,*) ' #### FLEXPART MODEL ERROR! AVAILABLE IILE #### '
+999   write(*,*) ' #### FLEXPART MODEL ERROR! AVAILABLE FILE #### '
   write(*,'(a)') '     '//path(4)(1:length(4))
   write(*,*) ' #### CANNOT BE OPENED           #### '
Index: src/readwind_nests.f90
===================================================================
--- src/readwind_nests.f90	(revision b0434e1968915cd3d4a8afcbf2cc275acecb9ff7)
+++ src/readwind_nests.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -347,20 +347,18 @@
 ! -add check for if one of clwc/ciwc missing (error),
 !    also if all 3 cw fields present, use qc and disregard the others
-! -use same flags readclouds/sumclouds as in mother grid? this assumes
-!    that both the nested and mother grids contain CW in same format
         if(isec1(6).eq.246) then  !! CLWC  Cloud liquid water content [kg/kg]
-          clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
-          ! readclouds=.true.
-          ! sumclouds=.false.
+          clwchn(i,j,nlev_ec-k+2,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
+          readclouds_nest(l)=.true.
+          sumclouds_nest(l)=.false.
         endif
         if(isec1(6).eq.247) then  !! CIWC  Cloud ice water content
-          ciwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
+          ciwchn(i,j,nlev_ec-k+2,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
         endif
 !ZHG end
 !ESO read qc (=clwc+ciwc)
         if(isec1(6).eq.201031) then  !! QC  Cloud liquid water content [kg/kg]
-          clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
-          ! readclouds=.true.
-          ! sumclouds=.true.
+          clwchn(i,j,nlev_ec-k+2,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
+          readclouds_nest(l)=.true.
+          sumclouds_nest(l)=.true.
         endif
 
Index: src/timemanager_mpi.f90
===================================================================
--- src/timemanager_mpi.f90	(revision 6a678e3fe206788e9704a724ec54e6d3dcc908e7)
+++ src/timemanager_mpi.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -237,10 +237,10 @@
     if (lmp_sync.and.lmp_use_reader.and.memstat.gt.0) then
       call mpif_gf_send_vars(memstat)
-      call mpif_gf_send_vars_nest(memstat)
+      if (numbnests>0) call mpif_gf_send_vars_nest(memstat)
 ! Version 2  (lmp_sync=.false., see below) is also used whenever 2 new fields are
 ! read (as at first time step), in which case async send/recv is impossible.
     else if (.not.lmp_sync.and.lmp_use_reader.and.memstat.ge.32) then
       call mpif_gf_send_vars(memstat)
-      call mpif_gf_send_vars_nest(memstat)
+      if (numbnests>0) call mpif_gf_send_vars_nest(memstat)
     end if
 
@@ -260,7 +260,4 @@
 ! Issued at start of each new field period. 
       if (memstat.ne.0.and.memstat.lt.32.and.lmp_use_reader) then
-! TODO: z0(7) changes with time, so should be dimension (numclass,2) to
-! allow transfer of the future value in the background
-        call MPI_Bcast(z0,numclass,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
         call mpif_gf_request
       end if
Index: src/verttransform.f90
===================================================================
--- src/verttransform.f90	(revision b0434e1968915cd3d4a8afcbf2cc275acecb9ff7)
+++ src/verttransform.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -67,46 +67,42 @@
   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(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,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,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 :: 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)
-  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
+  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
+  ! 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)
+  ! 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
+  ! CHARACTER(LEN=3)  :: aspec
+  ! integer :: virr=0
   real :: tot_cloud_h
+  real :: dbg_height(nzmax) 
 !ZHG
 
@@ -143,5 +139,5 @@
 
 
-    tvold(ixm,jym)=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(ixm,jym)=ps(ixm,jym,1,n)
@@ -182,4 +178,6 @@
     init=.false.
 
+    dbg_height = height
+
   endif
 
@@ -191,5 +189,5 @@
   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))/ &
+      tvold(ix,jy)=tt2(ix,jy,1,n)*(1.+0.378*ew*(td2(ix,jy,1,n))/ &
            ps(ix,jy,1,n))
     enddo
@@ -250,17 +248,21 @@
   qv(:,:,1,n)=qvh(:,:,1,n)
 !hg adding the cloud water 
-  clwc(:,:,1,n)=clwch(:,:,1,n)
-  if (.not.sumclouds) ciwc(:,:,1,n)=ciwch(:,:,1,n)   
+  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
-  clwc(:,:,nz,n)=clwch(:,:,nuvz,n)
-  if (.not.sumclouds) ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n)
+  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)
@@ -279,6 +281,8 @@
           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)
-          if (.not.sumclouds) ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
+          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)
@@ -309,7 +313,9 @@
                +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
-          if (.not.sumclouds) &
-               &ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz
+          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
@@ -601,5 +607,5 @@
             icloud_stats(ix,jy,3,n)= min(height(kz+1),height(kz))                     ! Cloud BOT height stats      [m]
 !ZHG 2015 extra for testing
-            clh(ix,jy,kz,n)=height(kz+1)-height(kz)
+!            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]
Index: src/verttransform_nests.f90
===================================================================
--- src/verttransform_nests.f90	(revision b0434e1968915cd3d4a8afcbf2cc275acecb9ff7)
+++ src/verttransform_nests.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -21,46 +21,50 @@
 
 subroutine verttransform_nests(n,uuhn,vvhn,wwhn,pvhn)
-  !                            i   i    i    i   i
-  !*****************************************************************************
-  !                                                                            *
-  !     This subroutine transforms temperature, dew point temperature and      *
-  !     wind components from eta to meter coordinates.                         *
-  !     The vertical wind component is transformed from Pa/s to m/s using      *
-  !     the conversion factor pinmconv.                                        *
-  !     In addition, this routine calculates vertical density gradients        *
-  !     needed for the parameterization of the turbulent velocities.           *
-  !     It is similar to verttransform, but makes the transformations for      *
-  !     the nested grids.                                                      *
-  !                                                                            *
-  !     Author: A. Stohl, G. Wotawa                                            *
-  !                                                                            *
-  !     12 August 1996                                                         *
-  !     Update: 16 January 1998                                                *
-  !                                                                            *
-  !     Major update: 17 February 1999                                         *
-  !     by G. Wotawa                                                           *
-  !                                                                            *
-  !     - Vertical levels for u, v and w are put together                      *
-  !     - Slope correction for vertical velocity: Modification of calculation  *
-  !       procedure                                                            *
-  !                                                                            *
-  !*****************************************************************************
-  !  Changes, Bernd C. Krueger, Feb. 2001:       (marked "C-cv")
-  !   Variables tthn and qvhn (on eta coordinates) from common block
-  !*****************************************************************************
-  ! Sabine Eckhardt, March 2007
-  ! add the variable cloud for use with scavenging - descr. in com_mod
-  !*****************************************************************************
-  !                                                                            *
-  ! Variables:                                                                 *
-  ! nxn,nyn,nuvz,nwz                field dimensions in x,y and z direction    *
-  ! uun                             wind components in x-direction [m/s]       *
-  ! vvn                             wind components in y-direction [m/s]       *
-  ! wwn                             wind components in z-direction [deltaeta/s]*
-  ! ttn                             temperature [K]                            *
-  ! pvn                             potential vorticity (pvu)                  *
-  ! psn                             surface pressure [Pa]                      *
-  !                                                                            *
-  !*****************************************************************************
+!                            i   i    i    i   i
+!*****************************************************************************
+!                                                                            *
+!     This subroutine transforms temperature, dew point temperature and      *
+!     wind components from eta to meter coordinates.                         *
+!     The vertical wind component is transformed from Pa/s to m/s using      *
+!     the conversion factor pinmconv.                                        *
+!     In addition, this routine calculates vertical density gradients        *
+!     needed for the parameterization of the turbulent velocities.           *
+!     It is similar to verttransform, but makes the transformations for      *
+!     the nested grids.                                                      *
+!                                                                            *
+!     Author: A. Stohl, G. Wotawa                                            *
+!                                                                            *
+!     12 August 1996                                                         *
+!     Update: 16 January 1998                                                *
+!                                                                            *
+!     Major update: 17 February 1999                                         *
+!     by G. Wotawa                                                           *
+!                                                                            *
+!     - Vertical levels for u, v and w are put together                      *
+!     - Slope correction for vertical velocity: Modification of calculation  *
+!       procedure                                                            *
+!                                                                            *
+!*****************************************************************************
+!  Changes, Bernd C. Krueger, Feb. 2001:       (marked "C-cv")
+!   Variables tthn and qvhn (on eta coordinates) from common block
+!*****************************************************************************
+! Sabine Eckhardt, March 2007
+! add the variable cloud for use with scavenging - descr. in com_mod
+!*****************************************************************************
+! ESO, 2016
+! -note that divide-by-zero occurs when nxmaxn,nymaxn etc. are larger than 
+!  the actual field dimensions
+!*****************************************************************************
+!                                                                            *
+! Variables:                                                                 *
+! nxn,nyn,nuvz,nwz                field dimensions in x,y and z direction    *
+! uun                             wind components in x-direction [m/s]       *
+! vvn                             wind components in y-direction [m/s]       *
+! wwn                             wind components in z-direction [deltaeta/s]*
+! ttn                             temperature [K]                            *
+! pvn                             potential vorticity (pvu)                  *
+! psn                             surface pressure [Pa]                      *
+!                                                                            *
+!*****************************************************************************
 
   use par_mod
@@ -69,256 +73,497 @@
   implicit none
 
-  integer :: ix,jy,kz,iz,n,l,kmin,kl,klp,ix1,jy1,ixp,jyp
-  integer :: rain_cloud_above,kz_inv
-  real :: f_qvsat,pressure,rh,lsp,convp
-  real :: wzlev(nwzmax),rhoh(nuvzmax),pinmconv(nzmax)
-  real :: uvwzlev(0:nxmaxn-1,0:nymaxn-1,nzmax)
-  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi,cosf
+  real,intent(in),dimension(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests) :: uuhn,vvhn,pvhn
+  real,intent(in),dimension(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests) :: wwhn
+
+  real,dimension(0:nxmaxn-1,0:nymaxn-1,nuvzmax) :: rhohn,uvzlev,wzlev
+  real,dimension(0:nxmaxn-1,0:nymaxn-1,nzmax) :: pinmconv
+  real,dimension(0:nxmaxn-1,0:nymaxn-1) :: tvold,pold,pint,tv
+  real,dimension(0:nymaxn-1) :: cosf
+
+  integer,dimension(0:nxmaxn-1,0:nymaxn-1) :: rain_cloud_above, idx
+
+  integer :: ix,jy,kz,iz,n,l,kmin,kl,klp,ix1,jy1,ixp,jyp,kz_inv
+  real :: f_qvsat,pressure,rh,lsp,convp,cloudh_min,prec
+
+  real :: ew,dz1,dz2,dz
   real :: dzdx,dzdy
   real :: dzdx1,dzdx2,dzdy1,dzdy2
-  real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
-  real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
-  real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
-  real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
   real,parameter :: const=r_air/ga
-
-
-  ! Loop over all nests
-  !********************
+  real :: tot_cloud_h
+  integer :: nxm1, nym1
+
+!  real,parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagnostics
+
+! Loop over all nests
+!********************
 
   do l=1,numbnests
-
-  ! Loop over the whole grid
-  !*************************
-
-  do jy=0,nyn(l)-1
-    do ix=0,nxn(l)-1
-
-      tvold=tt2n(ix,jy,1,n,l)*(1.+0.378*ew(td2n(ix,jy,1,n,l))/ &
-      psn(ix,jy,1,n,l))
-      pold=psn(ix,jy,1,n,l)
-      uvwzlev(ix,jy,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)*psn(ix,jy,1,n,l)
-        tv=tthn(ix,jy,kz,n,l)*(1.+0.608*qvhn(ix,jy,kz,n,l))
-        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
+    nxm1=nxn(l)-1 
+    nym1=nyn(l)-1 
+
+! Loop over the whole grid
+!*************************
+
+    do jy=0,nyn(l)-1
+      do ix=0,nxn(l)-1
+        tvold(ix,jy)=tt2n(ix,jy,1,n,l)*(1.+0.378*ew*(td2n(ix,jy,1,n,l))/ &
+             psn(ix,jy,1,n,l))
       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)
-
-  ! pinmconv=(h2-h1)/(p2-p1)
-
-      pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
-      ((aknew(2)+bknew(2)*psn(ix,jy,1,n,l))- &
-      (aknew(1)+bknew(1)*psn(ix,jy,1,n,l)))
-      do kz=2,nz-1
-        pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
-        ((aknew(kz+1)+bknew(kz+1)*psn(ix,jy,1,n,l))- &
-        (aknew(kz-1)+bknew(kz-1)*psn(ix,jy,1,n,l)))
-      end do
-      pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
-      ((aknew(nz)+bknew(nz)*psn(ix,jy,1,n,l))- &
-      (aknew(nz-1)+bknew(nz-1)*psn(ix,jy,1,n,l)))
-
-
-  ! Levels, where u,v,t and q are given
-  !************************************
-
-      uun(ix,jy,1,n,l)=uuhn(ix,jy,1,l)
-      vvn(ix,jy,1,n,l)=vvhn(ix,jy,1,l)
-      ttn(ix,jy,1,n,l)=tthn(ix,jy,1,n,l)
-      qvn(ix,jy,1,n,l)=qvhn(ix,jy,1,n,l)
-      pvn(ix,jy,1,n,l)=pvhn(ix,jy,1,l)
-      rhon(ix,jy,1,n,l)=rhoh(1)
-      uun(ix,jy,nz,n,l)=uuhn(ix,jy,nuvz,l)
-      vvn(ix,jy,nz,n,l)=vvhn(ix,jy,nuvz,l)
-      ttn(ix,jy,nz,n,l)=tthn(ix,jy,nuvz,n,l)
-      qvn(ix,jy,nz,n,l)=qvhn(ix,jy,nuvz,n,l)
-      pvn(ix,jy,nz,n,l)=pvhn(ix,jy,nuvz,l)
-      rhon(ix,jy,nz,n,l)=rhoh(nuvz)
-      kmin=2
-      do iz=2,nz-1
-        do kz=kmin,nuvz
-          if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
+    end do
+
+    pold(0:nxm1,0:nym1)=psn(0:nxm1,0:nym1,1,n,l)
+    uvzlev(0:nxm1,0:nym1,1)=0.
+    wzlev(0:nxm1,0:nym1,1)=0.
+    rhohn(0:nxm1,0:nym1,1)=pold(0:nxm1,0:nym1)/(r_air*tvold(0:nxm1,0:nym1))
+
+! Compute heights of eta levels
+!******************************
+
+    do kz=2,nuvz
+      pint(0:nxm1,0:nym1)=akz(kz)+bkz(kz)*psn(0:nxm1,0:nym1,1,n,l)
+      tv(0:nxm1,0:nym1)=tthn(0:nxm1,0:nym1,kz,n,l)*(1.+0.608*qvhn(0:nxm1,0:nym1,kz,n,l))
+      rhohn(0:nxm1,0:nym1,kz)=pint(0:nxm1,0:nym1)/(r_air*tv(0:nxm1,0:nym1))
+
+      where (abs(tv(0:nxm1,0:nym1)-tvold(0:nxm1,0:nym1)).gt.0.2) 
+        uvzlev(0:nxm1,0:nym1,kz)=uvzlev(0:nxm1,0:nym1,kz-1)+const*&
+             &log(pold(0:nxm1,0:nym1)/pint(0:nxm1,0:nym1))* &
+             (tv(0:nxm1,0:nym1)-tvold(0:nxm1,0:nym1))/&
+             &log(tv(0:nxm1,0:nym1)/tvold(0:nxm1,0:nym1))
+      elsewhere
+        uvzlev(0:nxm1,0:nym1,kz)=uvzlev(0:nxm1,0:nym1,kz-1)+const*&
+             &log(pold(0:nxm1,0:nym1)/pint(0:nxm1,0:nym1))*tv(0:nxm1,0:nym1)
+      endwhere
+
+      tvold(0:nxm1,0:nym1)=tv(0:nxm1,0:nym1)
+      pold(0:nxm1,0:nym1)=pint(0:nxm1,0:nym1)
+
+    end do
+
+    do kz=2,nwz-1
+      wzlev(0:nxm1,0:nym1,kz)=(uvzlev(0:nxm1,0:nym1,kz+1)+uvzlev(0:nxm1,0:nym1,kz))/2.
+    end do
+    wzlev(0:nxm1,0:nym1,nwz)=wzlev(0:nxm1,0:nym1,nwz-1)+ &
+         uvzlev(0:nxm1,0:nym1,nuvz)-uvzlev(0:nxm1,0:nym1,nuvz-1)
+
+
+    pinmconv(0:nxm1,0:nym1,1)=(uvzlev(0:nxm1,0:nym1,2))/ &
+         ((aknew(2)+bknew(2)*psn(0:nxm1,0:nym1,1,n,l))- &
+         (aknew(1)+bknew(1)*psn(0:nxm1,0:nym1,1,n,l)))
+    do kz=2,nz-1
+      pinmconv(0:nxm1,0:nym1,kz)=(uvzlev(0:nxm1,0:nym1,kz+1)-uvzlev(0:nxm1,0:nym1,kz-1))/ &
+           ((aknew(kz+1)+bknew(kz+1)*psn(0:nxm1,0:nym1,1,n,l))- &
+           (aknew(kz-1)+bknew(kz-1)*psn(0:nxm1,0:nym1,1,n,l)))
+    end do
+    pinmconv(0:nxm1,0:nym1,nz)=(uvzlev(0:nxm1,0:nym1,nz)-uvzlev(0:nxm1,0:nym1,nz-1))/ &
+         ((aknew(nz)+bknew(nz)*psn(0:nxm1,0:nym1,1,n,l))- &
+         (aknew(nz-1)+bknew(nz-1)*psn(0:nxm1,0:nym1,1,n,l)))
+
+! Levels, where u,v,t and q are given
+!************************************
+
+    uun(0:nxm1,0:nym1,1,n,l)=uuhn(0:nxm1,0:nym1,1,l)
+    vvn(0:nxm1,0:nym1,1,n,l)=vvhn(0:nxm1,0:nym1,1,l)
+    ttn(0:nxm1,0:nym1,1,n,l)=tthn(0:nxm1,0:nym1,1,n,l)
+    qvn(0:nxm1,0:nym1,1,n,l)=qvhn(0:nxm1,0:nym1,1,n,l)
+    if (readclouds_nest(l)) then
+      clwcn(0:nxm1,0:nym1,1,n,l)=clwchn(0:nxm1,0:nym1,1,n,l)
+      if (.not.sumclouds_nest(l)) ciwcn(0:nxm1,0:nym1,1,n,l)=ciwchn(0:nxm1,0:nym1,1,n,l)
+    end if
+    pvn(0:nxm1,0:nym1,1,n,l)=pvhn(0:nxm1,0:nym1,1,l)
+    rhon(0:nxm1,0:nym1,1,n,l)=rhohn(0:nxm1,0:nym1,1)
+
+    uun(0:nxm1,0:nym1,nz,n,l)=uuhn(0:nxm1,0:nym1,nuvz,l)
+    vvn(0:nxm1,0:nym1,nz,n,l)=vvhn(0:nxm1,0:nym1,nuvz,l)
+    ttn(0:nxm1,0:nym1,nz,n,l)=tthn(0:nxm1,0:nym1,nuvz,n,l)
+    qvn(0:nxm1,0:nym1,nz,n,l)=qvhn(0:nxm1,0:nym1,nuvz,n,l)
+    if (readclouds_nest(l)) then
+      clwcn(0:nxm1,0:nym1,nz,n,l)=clwchn(0:nxm1,0:nym1,nuvz,n,l)
+      if (.not.sumclouds_nest(l)) ciwcn(0:nxm1,0:nym1,nz,n,l)=ciwchn(0:nxm1,0:nym1,nuvz,n,l)
+    end if
+    pvn(0:nxm1,0:nym1,nz,n,l)=pvhn(0:nxm1,0:nym1,nuvz,l)
+    rhon(0:nxm1,0:nym1,nz,n,l)=rhohn(0:nxm1,0:nym1,nuvz)
+
+
+    kmin=2
+    idx=kmin
+    do iz=2,nz-1
+      do jy=0,nyn(l)-1
+        do ix=0,nxn(l)-1
+          if(height(iz).gt.uvzlev(ix,jy,nuvz)) then
             uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l)
             vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l)
             ttn(ix,jy,iz,n,l)=ttn(ix,jy,nz,n,l)
             qvn(ix,jy,iz,n,l)=qvn(ix,jy,nz,n,l)
+!hg adding the cloud water
+            if (readclouds_nest(l)) then
+              clwcn(ix,jy,iz,n,l)=clwcn(ix,jy,nz,n,l)
+              if (.not.sumclouds_nest(l)) ciwcn(ix,jy,iz,n,l)=ciwcn(ix,jy,nz,n,l)
+            end if
+!hg
             pvn(ix,jy,iz,n,l)=pvn(ix,jy,nz,n,l)
             rhon(ix,jy,iz,n,l)=rhon(ix,jy,nz,n,l)
-            goto 30
+          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
-          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
-          (height(iz).le.uvwzlev(ix,jy,kz))) then
-           dz1=height(iz)-uvwzlev(ix,jy,kz-1)
-           dz2=uvwzlev(ix,jy,kz)-height(iz)
-           dz=dz1+dz2
-           uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+ &
-           uuhn(ix,jy,kz,l)*dz1)/dz
-           vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+ &
-           vvhn(ix,jy,kz,l)*dz1)/dz
-           ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2+ &
-           tthn(ix,jy,kz,n,l)*dz1)/dz
-           qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2+ &
-           qvhn(ix,jy,kz,n,l)*dz1)/dz
-           pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+ &
-           pvhn(ix,jy,kz,l)*dz1)/dz
-           rhon(ix,jy,iz,n,l)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
-           kmin=kz
-           goto 30
+        enddo
+      enddo
+      do jy=0,nyn(l)-1
+        do ix=0,nxn(l)-1
+          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
+            uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+uuhn(ix,jy,kz,l)*dz1)/dz
+            vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+vvhn(ix,jy,kz,l)*dz1)/dz
+            ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2 &
+                 +tthn(ix,jy,kz,n,l)*dz1)/dz
+            qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2 &
+                 +qvhn(ix,jy,kz,n,l)*dz1)/dz
+!hg adding the cloud water
+            if (readclouds_nest(l)) then
+              clwcn(ix,jy,iz,n,l)=(clwchn(ix,jy,kz-1,n,l)*dz2+clwchn(ix,jy,kz,n,l)*dz1)/dz
+              if (.not.sumclouds_nest(l)) &
+                   &ciwcn(ix,jy,iz,n,l)=(ciwchn(ix,jy,kz-1,n,l)*dz2+ciwchn(ix,jy,kz,n,l)*dz1)/dz
+            end if
+!hg
+            pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+pvhn(ix,jy,kz,l)*dz1)/dz
+            rhon(ix,jy,iz,n,l)=(rhohn(ix,jy,kz-1)*dz2+rhohn(ix,jy,kz)*dz1)/dz
           endif
+        enddo
+      enddo
+    enddo
+
+
+
+!         do kz=kmin,nuvz
+!           if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
+!             uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l)
+!             vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l)
+!             ttn(ix,jy,iz,n,l)=ttn(ix,jy,nz,n,l)
+!             qvn(ix,jy,iz,n,l)=qvn(ix,jy,nz,n,l)
+!             pvn(ix,jy,iz,n,l)=pvn(ix,jy,nz,n,l)
+!             rhon(ix,jy,iz,n,l)=rhon(ix,jy,nz,n,l)
+!             goto 30
+!           endif
+!           if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
+!           (height(iz).le.uvwzlev(ix,jy,kz))) then
+!            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
+!            dz2=uvwzlev(ix,jy,kz)-height(iz)
+!            dz=dz1+dz2
+!            uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+ &
+!            uuhn(ix,jy,kz,l)*dz1)/dz
+!            vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+ &
+!            vvhn(ix,jy,kz,l)*dz1)/dz
+!            ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2+ &
+!            tthn(ix,jy,kz,n,l)*dz1)/dz
+!            qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2+ &
+!            qvhn(ix,jy,kz,n,l)*dz1)/dz
+!            pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+ &
+!            pvhn(ix,jy,kz,l)*dz1)/dz
+!            rhon(ix,jy,iz,n,l)=(rhohn(kz-1)*dz2+rhohn(kz)*dz1)/dz
+!            kmin=kz
+!            goto 30
+!           endif
+!         end do
+! 30      continue
+!       end do
+
+
+! Levels, where w is given
+!*************************
+
+    wwn(0:nxm1,0:nym1,1,n,l)=wwhn(0:nxm1,0:nym1,1,l)*pinmconv(0:nxm1,0:nym1,1)
+    wwn(0:nxm1,0:nym1,nz,n,l)=wwhn(0:nxm1,0:nym1,nwz,l)*pinmconv(0:nxm1,0:nym1,nz)
+    kmin=2
+    idx=kmin
+    do iz=2,nz
+      do jy=0,nym1
+        do ix=0,nxm1
+          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,nym1
+        do ix=0,nxm1
+          kz=idx(ix,jy)
+          dz1=height(iz)-wzlev(ix,jy,kz-1)
+          dz2=wzlev(ix,jy,kz)-height(iz)
+          dz=dz1+dz2
+          wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(ix,jy,kz-1)*dz2 &
+               +wwhn(ix,jy,kz,l)*pinmconv(ix,jy,kz)*dz1)/dz
+        enddo
+      enddo
+    enddo
+
+!       wwn(ix,jy,1,n,l)=wwhn(ix,jy,1,l)*pinmconv(1)
+!       wwn(ix,jy,nz,n,l)=wwhn(ix,jy,nwz,l)*pinmconv(nz)
+!       kmin=2
+!       do iz=2,nz
+!         do kz=kmin,nwz
+!           if ((height(iz).gt.wzlev(kz-1)).and. &
+!           (height(iz).le.wzlev(kz))) then
+!             dz1=height(iz)-wzlev(kz-1)
+!             dz2=wzlev(kz)-height(iz)
+!             dz=dz1+dz2
+!             wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(kz-1)*dz2 &
+!             +wwhn(ix,jy,kz,l)*pinmconv(kz)*dz1)/dz
+!             kmin=kz
+!             goto 40
+!           endif
+!         end do
+! 40      continue
+!       end do
+
+! Compute density gradients at intermediate levels
+!*************************************************
+
+    drhodzn(0:nxm1,0:nym1,1,n,l)=(rhon(0:nxm1,0:nym1,2,n,l)-rhon(0:nxm1,0:nym1,1,n,l))/ &
+         (height(2)-height(1))
+    do kz=2,nz-1
+      drhodzn(0:nxm1,0:nym1,kz,n,l)=(rhon(0:nxm1,0:nym1,kz+1,n,l)-rhon(0:nxm1,0:nym1,kz-1,n,l))/ &
+           (height(kz+1)-height(kz-1))
+    end do
+    drhodzn(0:nxm1,0:nym1,nz,n,l)=drhodzn(0:nxm1,0:nym1,nz-1,n,l)
+
+
+! drhodzn(ix,jy,1,n,l)=(rhon(ix,jy,2,n,l)-rhon(ix,jy,1,n,l))/ &
+! (height(2)-height(1))
+! do kz=2,nz-1
+!   drhodzn(ix,jy,kz,n,l)=(rhon(ix,jy,kz+1,n,l)- &
+!   rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
+! end do
+! drhodzn(ix,jy,nz,n,l)=drhodzn(ix,jy,nz-1,n,l)
+
+
+
+!   end do
+! end do
+
+
+!****************************************************************
+! Compute slope of eta levels in windward direction and resulting
+! vertical wind correction
+!****************************************************************
+
+    do jy=1,nyn(l)-2
+!    cosf=cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
+      cosf(jy)=1./cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
+    end do
+
+    kmin=2
+    idx=kmin
+    do iz=2,nz-1
+      do jy=1,nyn(l)-2
+        do ix=1,nxn(l)-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,nyn(l)-2
+        do ix=1,nxn(l)-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
+
+          wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l)+(dzdx*uun(ix,jy,iz,n,l)*dxconst*xresoln(l)*cosf(jy)+&
+               &dzdy*vvn(ix,jy,iz,n,l)*dyconst*yresoln(l))
+
         end do
-30      continue
+
       end do
-
-
-  ! Levels, where w is given
-  !*************************
-
-      wwn(ix,jy,1,n,l)=wwhn(ix,jy,1,l)*pinmconv(1)
-      wwn(ix,jy,nz,n,l)=wwhn(ix,jy,nwz,l)*pinmconv(nz)
-      kmin=2
-      do iz=2,nz
-        do kz=kmin,nwz
-          if ((height(iz).gt.wzlev(kz-1)).and. &
-          (height(iz).le.wzlev(kz))) then
-            dz1=height(iz)-wzlev(kz-1)
-            dz2=wzlev(kz)-height(iz)
-            dz=dz1+dz2
-            wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(kz-1)*dz2 &
-            +wwhn(ix,jy,kz,l)*pinmconv(kz)*dz1)/dz
-            kmin=kz
-            goto 40
-          endif
+    end do
+
+
+!   do jy=1,nyn(l)-2
+!     do ix=1,nxn(l)-2
+!       kmin=2
+!       do iz=2,nz-1
+
+!         ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/cosf(jy)
+!         vi=vvn(ix,jy,iz,n,l)*dyconst*yresoln(l)
+
+!         do kz=kmin,nz
+!           if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
+!                (height(iz).le.uvwzlev(ix,jy,kz))) then
+!             dz1=height(iz)-uvwzlev(ix,jy,kz-1)
+!             dz2=uvwzlev(ix,jy,kz)-height(iz)
+!             dz=dz1+dz2
+!             kl=kz-1
+!             klp=kz
+!             kmin=kz
+!             goto 47
+!           endif
+!         end do
+
+! 47      ix1=ix-1
+!         jy1=jy-1
+!         ixp=ix+1
+!         jyp=jy+1
+
+!         dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
+!         dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
+!         dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
+
+!         dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
+!         dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
+!         dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
+
+!         wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l)+(dzdx*ui+dzdy*vi)
+
+!       end do
+
+!     end do
+!   end do
+
+
+!write (*,*) 'initializing nested cloudsn, n:',n
+!   create a cloud and rainout/washout field, cloudsn occur where rh>80%
+
+
+!***********************************************************************************  
+    if (readclouds_nest(l)) 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(*,*) 'Nested ECMWF fields: using cloud water'
+      clwn(0:nxm1,0:nym1,:,n,l)=0.0
+!    icloud_stats(0:nxm1,0:nym1,:,n)=0.0
+      clw4n(0:nxm1,0:nym1,n,l)=0.0
+      cloudsn(0:nxm1,0:nym1,:,n,l)=0
+! If water/ice are read separately into clwc and ciwc, store sum in clwcn
+      if (.not.sumclouds_nest(l)) then 
+        clwcn(0:nxm1,0:nym1,:,n,l) = clwcn(0:nxm1,0:nym1,:,n,l) + ciwcn(0:nxm1,0:nym1,:,n,l)
+      end if
+      do jy=0,nym1
+        do ix=0,nxm1
+          lsp=lsprecn(ix,jy,1,n,l)
+          convp=convprecn(ix,jy,1,n,l)
+          prec=lsp+convp
+          tot_cloud_h=0
+! Find clouds in the vertical
+          do kz=1, nz-1 !go from top to bottom
+            if (clwcn(ix,jy,kz,n,l).gt.0) then      
+! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3 
+              clwn(ix,jy,kz,n,l)=(clwcn(ix,jy,kz,n,l)*rhon(ix,jy,kz,n,l))*(height(kz+1)-height(kz))
+              tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz)) 
+              clw4n(ix,jy,n,l) = clw4n(ix,jy,n,l)+clwn(ix,jy,kz,n,l)
+!            icloud_stats(ix,jy,4,n)= icloud_stats(ix,jy,4,n)+clw(ix,jy,kz,n)          ! Column cloud water [m3/m3]
+!           icloud_stats(ix,jy,3,n)= min(height(kz+1),height(kz))                     ! Cloud BOT height stats      [m]
+              cloudh_min=min(height(kz+1),height(kz))
+!ZHG 2015 extra for testing
+!            clh(ix,jy,kz,n)=height(kz+1)-height(kz)
+!            icloud_stats(ix,jy,1,n)=icloud_stats(ix,jy,1,n)+(height(kz+1)-height(kz)) ! Cloud total vertical extent [m]
+!            icloud_stats(ix,jy,2,n)= max(icloud_stats(ix,jy,2,n),height(kz))          ! Cloud TOP height            [m]
+!ZHG
+            endif
+          end do
+
+! If Precipitation. Define removal type in the vertical
+          if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
+
+            do kz=nz,1,-1 !go Bottom up!
+              if (clwn(ix,jy,kz,n,l).gt.0.0) then ! is in cloud
+                cloudshn(ix,jy,n,l)=cloudshn(ix,jy,n,l)+height(kz)-height(kz-1) 
+                cloudsn(ix,jy,kz,n,l)=1                               ! is a cloud
+                if (lsp.ge.convp) then
+                  cloudsn(ix,jy,kz,n,l)=3                            ! lsp in-cloud
+                else
+                  cloudsn(ix,jy,kz,n,l)=2                             ! convp in-cloud
+                endif                                              ! convective or large scale
+              else if((clwn(ix,jy,kz,n,l).le.0.0) .and. (cloudh_min.ge.height(kz))) then ! is below cloud
+                if (lsp.ge.convp) then
+                  cloudsn(ix,jy,kz,n,l)=5                             ! lsp dominated washout
+                else
+                  cloudsn(ix,jy,kz,n,l)=4                             ! convp dominated washout
+                endif                                              ! convective or large scale 
+              endif
+
+              if (height(kz).ge. 19000) then                        ! set a max height for removal
+                cloudsn(ix,jy,kz,n,l)=0
+              endif !clw>0
+            end do !nz
+          endif ! precipitation
         end do
-40      continue
       end do
-
-  ! Compute density gradients at intermediate levels
-  !*************************************************
-
-      drhodzn(ix,jy,1,n,l)=(rhon(ix,jy,2,n,l)-rhon(ix,jy,1,n,l))/ &
-      (height(2)-height(1))
-      do kz=2,nz-1
-        drhodzn(ix,jy,kz,n,l)=(rhon(ix,jy,kz+1,n,l)- &
-        rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
+!********************************************************************
+    else ! old method:
+!********************************************************************
+      write(*,*) 'Nested fields: using cloud water from Parameterization'
+      do jy=0,nyn(l)-1
+        do ix=0,nxn(l)-1
+          rain_cloud_above(ix,jy)=0
+          lsp=lsprecn(ix,jy,1,n,l)
+          convp=convprecn(ix,jy,1,n,l)
+          cloudshn(ix,jy,n,l)=0
+          do kz_inv=1,nz-1
+            kz=nz-kz_inv+1
+            pressure=rhon(ix,jy,kz,n,l)*r_air*ttn(ix,jy,kz,n,l)
+            rh=qvn(ix,jy,kz,n,l)/f_qvsat(pressure,ttn(ix,jy,kz,n,l))
+            cloudsn(ix,jy,kz,n,l)=0
+            if (rh.gt.0.8) then ! in cloud
+              if ((lsp.gt.0.01).or.(convp.gt.0.01)) then
+                rain_cloud_above(ix,jy)=1
+                cloudshn(ix,jy,n,l)=cloudshn(ix,jy,n,l)+ &
+                     height(kz)-height(kz-1)
+                if (lsp.ge.convp) then
+                  cloudsn(ix,jy,kz,n,l)=3 ! lsp dominated rainout
+                else
+                  cloudsn(ix,jy,kz,n,l)=2 ! convp dominated rainout
+                endif
+              else ! no precipitation
+                cloudsn(ix,jy,kz,n,l)=1 ! cloud
+              endif
+            else ! no cloud
+              if (rain_cloud_above(ix,jy).eq.1) then ! scavenging
+                if (lsp.ge.convp) then
+                  cloudsn(ix,jy,kz,n,l)=5 ! lsp dominated washout
+                else
+                  cloudsn(ix,jy,kz,n,l)=4 ! convp dominated washout
+                endif
+              endif
+            endif
+          end do
+        end do
       end do
-      drhodzn(ix,jy,nz,n,l)=drhodzn(ix,jy,nz-1,n,l)
-
-    end do
-  end do
-
-
-  !****************************************************************
-  ! Compute slope of eta levels in windward direction and resulting
-  ! vertical wind correction
-  !****************************************************************
-
-  do jy=1,nyn(l)-2
-    cosf=cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
-    do ix=1,nxn(l)-2
-
-      kmin=2
-      do iz=2,nz-1
-
-        ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/cosf
-        vi=vvn(ix,jy,iz,n,l)*dyconst*yresoln(l)
-
-        do kz=kmin,nz
-          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
-               (height(iz).le.uvwzlev(ix,jy,kz))) then
-            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
-            dz2=uvwzlev(ix,jy,kz)-height(iz)
-            dz=dz1+dz2
-            kl=kz-1
-            klp=kz
-            kmin=kz
-            goto 47
-          endif
-        end do
-
-47      ix1=ix-1
-        jy1=jy-1
-        ixp=ix+1
-        jyp=jy+1
-
-        dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
-        dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
-        dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
-
-        dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
-        dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
-        dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
-
-        wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l)+(dzdx*ui+dzdy*vi)
-
-      end do
-
-    end do
-  end do
-
-
-  !write (*,*) 'initializing nested cloudsn, n:',n
-  !   create a cloud and rainout/washout field, cloudsn occur where rh>80%
-  write(*,*) 'Nested fields: using cloud water from Parameterization'
-  do jy=0,nyn(l)-1
-    do ix=0,nxn(l)-1
-      rain_cloud_above=0
-      lsp=lsprecn(ix,jy,1,n,l)
-      convp=convprecn(ix,jy,1,n,l)
-      cloudsnh(ix,jy,n,l)=0
-      do kz_inv=1,nz-1
-         kz=nz-kz_inv+1
-         pressure=rhon(ix,jy,kz,n,l)*r_air*ttn(ix,jy,kz,n,l)
-         rh=qvn(ix,jy,kz,n,l)/f_qvsat(pressure,ttn(ix,jy,kz,n,l))
-         cloudsn(ix,jy,kz,n,l)=0
-         if (rh.gt.0.8) then ! in cloud
-            if ((lsp.gt.0.01).or.(convp.gt.0.01)) then
-               rain_cloud_above=1
-               cloudsnh(ix,jy,n,l)=cloudsnh(ix,jy,n,l)+ &
-               height(kz)-height(kz-1)
-               if (lsp.ge.convp) then
-                  cloudsn(ix,jy,kz,n,l)=3 ! lsp dominated rainout
-               else
-                  cloudsn(ix,jy,kz,n,l)=2 ! convp dominated rainout
-               endif
-            else ! no precipitation
-                  cloudsn(ix,jy,kz,n,l)=1 ! cloud
-            endif
-         else ! no cloud
-            if (rain_cloud_above.eq.1) then ! scavenging
-               if (lsp.ge.convp) then
-                  cloudsn(ix,jy,kz,n,l)=5 ! lsp dominated washout
-               else
-                  cloudsn(ix,jy,kz,n,l)=4 ! convp dominated washout
-               endif
-            endif
-         endif
-      end do
-    end do
-  end do
-
-  end do
+    end if
+
+  end do ! end loop over nests
 
 end subroutine verttransform_nests
Index: src/wetdepo.f90
===================================================================
--- src/wetdepo.f90	(revision 41d8574bea1b0f91e61dffd5dd1d17882407f826)
+++ src/wetdepo.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -21,46 +21,46 @@
 
 subroutine wetdepo(itime,ltsample,loutnext)
-  !                  i      i        i
-  !*****************************************************************************
-  !                                                                            *
-  ! Calculation of wet deposition using the concept of scavenging coefficients.*
-  ! For lack of detailed information, washout and rainout are jointly treated. *
-  ! It is assumed that precipitation does not occur uniformly within the whole *
-  ! grid cell, but that only a fraction of the grid cell experiences rainfall. *
-  ! This fraction is parameterized from total cloud cover and rates of large   *
-  ! scale and convective precipitation.                                        *
-  !                                                                            *
-  !    Author: A. Stohl                                                        *
-  !                                                                            *
-  !    1 December 1996                                                         *
-  !                                                                            *
-  ! Correction by Petra Seibert, Sept 2002:                                    *
-  ! use centred precipitation data for integration                             *
-  ! Code may not be correct for decay of deposition!                           *
-  !                                                                            *
-  !*****************************************************************************
-  !                                                                            *
-  ! Variables:                                                                 *
-  ! cc [0-1]           total cloud cover                                       *
-  ! convp [mm/h]       convective precipitation rate                           *
-  ! grfraction [0-1]   fraction of grid, for which precipitation occurs        *
-  ! ix,jy              indices of output grid cell for each particle           *
-  ! itime [s]          actual simulation time [s]                              *
-  ! jpart              particle index                                          *
-  ! ldeltat [s]        interval since radioactive decay was computed           *
-  ! lfr, cfr           area fraction covered by precipitation for large scale  *
-  !                    and convective precipitation (dependent on prec. rate)  *
-  ! loutnext [s]       time for which gridded deposition is next output        *
-  ! loutstep [s]       interval at which gridded deposition is output          *
-  ! lsp [mm/h]         large scale precipitation rate                          *
-  ! ltsample [s]       interval over which mass is deposited                   *
-  ! prec [mm/h]        precipitation rate in subgrid, where precipitation occurs*
-  ! wetdeposit         mass that is wet deposited                              *
-  ! wetgrid            accumulated deposited mass on output grid               *
-  ! wetscav            scavenging coefficient                                  *
-  !                                                                            *
-  ! Constants:                                                                 *
-  !                                                                            *
-  !*****************************************************************************
+!                  i      i        i
+!*****************************************************************************
+!                                                                            *
+! Calculation of wet deposition using the concept of scavenging coefficients.*
+! For lack of detailed information, washout and rainout are jointly treated. *
+! It is assumed that precipitation does not occur uniformly within the whole *
+! grid cell, but that only a fraction of the grid cell experiences rainfall. *
+! This fraction is parameterized from total cloud cover and rates of large   *
+! scale and convective precipitation.                                        *
+!                                                                            *
+!    Author: A. Stohl                                                        *
+!                                                                            *
+!    1 December 1996                                                         *
+!                                                                            *
+! Correction by Petra Seibert, Sept 2002:                                    *
+! use centred precipitation data for integration                             *
+! Code may not be correct for decay of deposition!                           *
+!                                                                            *
+!*****************************************************************************
+!                                                                            *
+! Variables:                                                                 *
+! cc [0-1]           total cloud cover                                       *
+! convp [mm/h]       convective precipitation rate                           *
+! grfraction [0-1]   fraction of grid, for which precipitation occurs        *
+! ix,jy              indices of output grid cell for each particle           *
+! itime [s]          actual simulation time [s]                              *
+! jpart              particle index                                          *
+! ldeltat [s]        interval since radioactive decay was computed           *
+! lfr, cfr           area fraction covered by precipitation for large scale  *
+!                    and convective precipitation (dependent on prec. rate)  *
+! loutnext [s]       time for which gridded deposition is next output        *
+! loutstep [s]       interval at which gridded deposition is output          *
+! lsp [mm/h]         large scale precipitation rate                          *
+! ltsample [s]       interval over which mass is deposited                   *
+! prec [mm/h]        precipitation rate in subgrid, where precipitation occurs*
+! wetdeposit         mass that is wet deposited                              *
+! wetgrid            accumulated deposited mass on output grid               *
+! wetscav            scavenging coefficient                                  *
+!                                                                            *
+! Constants:                                                                 *
+!                                                                            *
+!*****************************************************************************
 
   use point_mod
@@ -71,5 +71,6 @@
 
   integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
-  integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v
+  integer :: ngrid,itage,nage,hz,il,interp_time, n
+  integer(kind=1) :: clouds_v
   integer :: ks, kp
 !  integer :: n1,n2, icbot,ictop, indcloud !TEST
@@ -79,5 +80,5 @@
   real :: wetdeposit(maxspec),restmass
   real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
-  !save lfr,cfr
+!save lfr,cfr
 
   real, parameter :: lfr(5) = (/ 0.5,0.65,0.8,0.9,0.95/)
@@ -91,9 +92,9 @@
   integer :: blc_count, inc_count
   real    :: Si_dummy, wetscav_dummy
-
-
-
-  ! Compute interval since radioactive decay of deposited mass was computed
-  !************************************************************************
+  logical :: readclouds_this_nest
+
+
+! Compute interval since radioactive decay of deposited mass was computed
+!************************************************************************
 
   if (itime.le.loutnext) then
@@ -103,6 +104,6 @@
   endif
 
-  ! Loop over all particles
-  !************************
+! Loop over all particles
+!************************
 
   blc_count=0
@@ -118,5 +119,5 @@
     endif
 
-  ! Determine age class of the particle
+! Determine age class of the particle
     itage=abs(itra1(jpart)-itramem(jpart))
     do nage=1,nageclass
@@ -126,11 +127,11 @@
 
 
-  ! Determine which nesting level to be used
-  !*****************************************
+! Determine which nesting level to be used
+!*****************************************
 
     ngrid=0
     do j=numbnests,1,-1
       if ((xtra1(jpart).gt.xln(j)).and.(xtra1(jpart).lt.xrn(j)).and. &
-      (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
+           (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
         ngrid=j
         goto 23
@@ -140,6 +141,6 @@
 
 
-  ! Determine nested grid coordinates
-  !**********************************
+! Determine nested grid coordinates
+!**********************************
 
     if (ngrid.gt.0) then
@@ -148,4 +149,9 @@
       ix=int(xtn)
       jy=int(ytn)
+      if (readclouds_nest(ngrid)) then
+        readclouds_this_nest=.true.
+      else
+        readclouds_this_nest=.false.
+      end if
     else
       ix=int(xtra1(jpart))
@@ -154,18 +160,18 @@
 
 
-  ! Interpolate large scale precipitation, convective precipitation and
-  ! total cloud cover
-  ! Note that interpolated time refers to itime-0.5*ltsample [PS]
-  !********************************************************************
+! Interpolate large scale precipitation, convective precipitation and
+! total cloud cover
+! Note that interpolated time refers to itime-0.5*ltsample [PS]
+!********************************************************************
     interp_time=nint(itime-0.5*ltsample)
-    
+
     if (ngrid.eq.0) then
       call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
-      1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
-      memtime(1),memtime(2),interp_time,lsp,convp,cc)
+           1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
+           memtime(1),memtime(2),interp_time,lsp,convp,cc)
     else
       call interpol_rain_nests(lsprecn,convprecn,tccn, &
-      nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
-      memtime(1),memtime(2),interp_time,lsp,convp,cc)
+           nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
+           memtime(1),memtime(2),interp_time,lsp,convp,cc)
     endif
 
@@ -173,12 +179,13 @@
     if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
 
-  ! get the level were the actual particle is in
+! get the level were the actual particle is in
     do il=2,nz
       if (height(il).gt.ztra1(jpart)) then
         hz=il-1
-        goto 26
+!        goto 26
+        exit
       endif
     end do
-26  continue
+!26  continue
 
     n=memind(2)
@@ -190,20 +197,19 @@
       clouds_h=cloudsh(ix,jy,n)
     else
-      ! new removal not implemented for nests yet 
       clouds_v=cloudsn(ix,jy,hz,n,ngrid)
-      clouds_h=cloudsnh(ix,jy,n,ngrid)
-    endif
-
-  ! if there is no precipitation or the particle is above the clouds no
-  ! scavenging is done
+      clouds_h=cloudshn(ix,jy,n,ngrid)
+    endif
+
+! if there is no precipitation or the particle is above the clouds no
+! scavenging is done
 
     if (clouds_v.le.1) goto 20
-  
-  ! 1) Parameterization of the the area fraction of the grid cell where the
-  !    precipitation occurs: the absolute limit is the total cloud cover, but
-  !    for low precipitation rates, an even smaller fraction of the grid cell
-  !    is used. Large scale precipitation occurs over larger areas than
-  !    convective precipitation.
-  !**************************************************************************
+
+! 1) Parameterization of the the area fraction of the grid cell where the
+!    precipitation occurs: the absolute limit is the total cloud cover, but
+!    for low precipitation rates, an even smaller fraction of the grid cell
+!    is used. Large scale precipitation occurs over larger areas than
+!    convective precipitation.
+!**************************************************************************
 
     if (lsp.gt.20.) then
@@ -232,7 +238,7 @@
 
 
-  !ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp 
-  ! Tentatively differentiate the grfraction for lsp and convp for treating differently the two forms
-  ! for now they are treated the same
+!ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp 
+! Tentatively differentiate the grfraction for lsp and convp for treating differently the two forms
+! for now they are treated the same
     grfraction(1)=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp))
     grfraction(2)=max(0.05,cc*(lfr(i)))
@@ -240,6 +246,6 @@
 
 
-  ! 2) Computation of precipitation rate in sub-grid cell
-  !******************************************************
+! 2) Computation of precipitation rate in sub-grid cell
+!******************************************************
     prec(1)=(lsp+convp)/grfraction(1)
     prec(2)=(lsp)/grfraction(2)
@@ -247,7 +253,7 @@
 
 
-  ! 3) Computation of scavenging coefficients for all species
-  !    Computation of wet deposition
-  !**********************************************************
+! 3) Computation of scavenging coefficients for all species
+!    Computation of wet deposition
+!**********************************************************
 
     do ks=1,nspec      ! loop over species
@@ -255,5 +261,5 @@
       wetscav=0.   
 
-      !ZHG test if it nested?
+!ZHG test if it nested?
       if (ngrid.gt.0) then
         act_temp=ttn(ix,jy,hz,n,ngrid)
@@ -261,9 +267,9 @@
         act_temp=tt(ix,jy,hz,n)
       endif
-     
-
-  !***********************
-  ! BELOW CLOUD SCAVENGING
-  !***********************  
+
+
+!***********************
+! BELOW CLOUD SCAVENGING
+!***********************  
       if (clouds_v.ge.4) then !below cloud
 
@@ -274,5 +280,5 @@
           if (dquer(ks) .le. 0.) then  !is gas
             wetscav=weta(ks)*prec(1)**wetb(ks)
-            
+
 !AEROSOL
           else !is particle
@@ -281,15 +287,15 @@
 ! for particles larger than 10 um use the largest size defined in the parameterizations (10um)
             dquer_m=min(10.,dquer(ks))/1000000. !conversion from um to m
-            if (act_temp .ge. 273 .and. weta(ks).gt.0.)  then !Rain 
-              ! ZHG 2014 : Particle RAIN scavenging coefficient based on Laakso et al 2003, 
-              ! the below-cloud scavenging (rain efficienty) 
-              ! parameter A (=weta) from SPECIES file
+            if (act_temp .ge. 273. .and. weta(ks).gt.0.)  then !Rain 
+! ZHG 2014 : Particle RAIN scavenging coefficient based on Laakso et al 2003, 
+! the below-cloud scavenging (rain efficienty) 
+! parameter A (=weta) from SPECIES file
               wetscav= weta(ks)*10**(bclr(1)+ (bclr(2)*(log10(dquer_m))**(-4))+(bclr(3)*(log10(dquer_m))**(-3))+ (bclr(4)* &
                    (log10(dquer_m))**(-2))+ (bclr(5)*(log10(dquer_m))**(-1))+ bclr(6)* (prec(1))**(0.5))
 
-            elseif (act_temp .lt. 273 .and. wetb(ks).gt.0.)  then ! Snow
-              ! ZHG 2014 : Particle SNOW scavenging coefficient based on Kyro et al 2009, 
-              ! the below-cloud scavenging (Snow efficiency) 
-              ! parameter B (=wetb) from SPECIES file
+            elseif (act_temp .lt. 273. .and. wetb(ks).gt.0.)  then ! Snow
+! ZHG 2014 : Particle SNOW scavenging coefficient based on Kyro et al 2009, 
+! the below-cloud scavenging (Snow efficiency) 
+! parameter B (=wetb) from SPECIES file
               wetscav= wetb(ks)*10**(bcls(1)+ (bcls(2)*(log10(dquer_m))**(-4))+(bcls(3)*(log10(dquer_m))**(-3))+ (bcls(4)* &
                    (log10(dquer_m))**(-2))+ (bcls(5)*(log10(dquer_m))**(-1))+ bcls(6)* (prec(1))**(0.5))
@@ -304,7 +310,7 @@
 
 
-  !********************
-  ! IN CLOUD SCAVENGING
-      !******************************************************
+!********************
+! IN CLOUD SCAVENGING
+!********************
       if (clouds_v.lt.4) then ! In-cloud
 ! NIK 13 may 2015: only do incloud if positive in-cloud scavenging parameters are given in species file
@@ -315,46 +321,47 @@
           if (wetb_in(ks).lt.0.) wetb_in(ks)=0.
 
-          !ZHG 2015 Cloud liquid & ice water (CLWC+CIWC) from ECMWF
-          if (readclouds) then                  !icloud_stats(ix,jy,4,n) has units kg/m2
-!            cl =icloud_stats(ix,jy,4,n)*(grfraction(1)/cc)
-! ESO: stop using icloud_stats
-            cl =clw4(ix,jy,n)*(grfraction(1)/cc)
+!ZHG 2015 Cloud liquid & ice water (CLWC+CIWC) from ECMWF
+! nested fields
+          if (ngrid.gt.0.and.readclouds_this_nest) then
+            cl=clw4n(ix,jy,n,ngrid)*(grfraction(1)/cc)
+          else if (ngrid.eq.0.and.readclouds) then
+            cl=clw4(ix,jy,n)*(grfraction(1)/cc)
           else                                  !parameterize cloudwater m2/m3
-            !ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF
+!ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF
             cl=1.6E-6*prec(1)**0.36
           endif
 
-            !ZHG: Calculate the partition between liquid and water phase water. 
-            if (act_temp .le. 253) then
-              liq_frac=0
-            else if (act_temp .ge. 273) then
-              liq_frac=1
-            else
-              liq_frac =((act_temp-273)/(273-253))**2
-            endif
-           ! ZHG: Calculate the aerosol partition based on cloud phase and Ai and Bi
-            frac_act = liq_frac*weta_in(ks) +(1-liq_frac)*wetb_in(ks)
- 
-  !ZHG Use the activated fraction and the liqid water to calculate the washout
-
-  ! AEROSOL
-          !**************************************************
+!ZHG: Calculate the partition between liquid and water phase water. 
+          if (act_temp .le. 253.) then
+            liq_frac=0
+          else if (act_temp .ge. 273.) then
+            liq_frac=1
+          else
+            liq_frac =((act_temp-273.)/(273.-253.))**2.
+          end if
+! ZHG: Calculate the aerosol partition based on cloud phase and Ai and Bi
+          frac_act = liq_frac*weta_in(ks) +(1-liq_frac)*wetb_in(ks)
+
+!ZHG Use the activated fraction and the liqid water to calculate the washout
+
+! AEROSOL
+!**************************************************
           if (dquer(ks).gt. 0.) then ! is particle
 
             S_i= frac_act/cl
 
-          !*********************
-  ! GAS
+!*********************
+! GAS
           else ! is gas
-               
+
             cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
-            !REPLACE to switch old/ new scheme 
-            ! S_i=frac_act/cle
+!REPLACE to switch old/ new scheme 
+! S_i=frac_act/cle
             S_i=1/cle
           endif ! gas or particle
 
-  ! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
-           !OLD 
-          if (readclouds) then
+! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
+!OLD 
+          if ((readclouds.and.ngrid.eq.0).or.(readclouds_this_nest.and.ngrid.gt.0)) then
             wetscav=S_i*(prec(1)/3.6E6)
           else
@@ -362,43 +369,14 @@
           endif
 
-!ZHG 2015 TEST
-!          Si_dummy=frac_act/2E-7*prec(1)**0.36
-!           wetscav_dummy=Si_dummy*(prec(1)/3.6E6)/clouds_h
-!           if (clouds_v.lt.4) then
-!           talltest=talltest+1
-!if(talltest .eq. 1) OPEN(UNIT=199, FILE=utfil,FORM='FORMATTED',STATUS = 'UNKNOWN')
-!if(talltest .lt. 100001)  write(199,*) prec(1)/3.6E6, cl, clouds_h*2E-7*prec(1)**0.36,clouds_v,ytra1(jpart)-90
-!if(talltest .lt. 100001)  write(199,*) wetscav, wetscav_dummy,prec(1),ytra1(jpart)-90,clouds_v,cl
-!if(talltest .eq. 100001) CLOSE(199)
-!if(talltest .eq. 100001) STOP
-!
-!write(*,*)  'PREC kg/m2s CLOUD kg/m2', (prec(1)/3.6E6), cl !, '2E-7*prec(1)**0.36',  2E-7*prec(1)**0.36,'2E-7*prec(1)**0.36*clouds_h',2E-7*prec(1)**0.36*clouds_h
-!write(*,*)  'PREC kg/m2s LSP+convp kg/m2', prec(1), convp+lsp
-!write(*,*)  wetscav, wetscav_dummy
-!write(*,*) cc, grfraction(1), cc/grfraction(1)
-
-!write(*,*)  'Lmbda_old', (prec(1)/3.6E6)/(clouds_h*2E-7*prec(1)**0.36)
-
-
-!write(*,*) '**************************************************'
-!write(*,*)  'clouds_h', clouds_h, 'clouds_v',clouds_v,'abs(ltsample)', abs(ltsample)
-!write(*,*)  'readclouds', readclouds, 'wetscav',wetscav, 'wetscav_dummy', wetscav_dummy
-!write(*,*)  'S_i', S_i , 'Si_dummy', Si_dummy, 'prec(1)', prec(1) 
-
-
-!           write(*,*) 'PRECIPITATION ,cl  ECMWF , cl PARAMETIZED, clouds_v, lat' &
-!                      ,prec(1)/3.6E6, cl, clouds_h*2E-7*prec(1)**0.36,clouds_v,ytra1(jpart)-90
-
-!endif 
 
         endif ! positive in-cloud scavenging parameters given in Species file
       endif !incloud
 !END ZHG TEST
-     
-  !**************************************************
-  ! CALCULATE DEPOSITION 
-  !**************************************************
-  !     if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
-  !     +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
+
+!**************************************************
+! CALCULATE DEPOSITION 
+!**************************************************
+!     if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
+!     +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
 
       if (wetscav.gt.0.) then
@@ -421,12 +399,12 @@
       if (restmass .gt. smallnum) then
         xmass1(jpart,ks)=restmass
-  !   depostatistic
-  !   wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
-  !   depostatistic
+!   depostatistic
+!   wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
+!   depostatistic
       else
         xmass1(jpart,ks)=0.
       endif
-  !   Correct deposited mass to the last time step when radioactive decay of
-  !   gridded deposited mass was calculated
+!   Correct deposited mass to the last time step when radioactive decay of
+!   gridded deposited mass was calculated
       if (decay(ks).gt.0.) then
         wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
@@ -436,7 +414,7 @@
     end do !all species
 
-  ! Sabine Eckhardt, June 2008 create deposition runs only for forward runs
-  ! Add the wet deposition to accumulated amount on output grid and nested output grid
-  !*****************************************************************************
+! Sabine Eckhardt, June 2008 create deposition runs only for forward runs
+! Add the wet deposition to accumulated amount on output grid and nested output grid
+!*****************************************************************************
 
     if (ldirect.eq.1) then
@@ -450,5 +428,5 @@
   end do ! all particles
 
-  ! count the total number of below-cloud and in-cloud occurences:
+! count the total number of below-cloud and in-cloud occurences:
   tot_blc_count=tot_blc_count+blc_count
   tot_inc_count=tot_inc_count+inc_count
Index: src/wetdepokernel_nest.f90
===================================================================
--- src/wetdepokernel_nest.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/wetdepokernel_nest.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
@@ -20,7 +20,6 @@
 !**********************************************************************
 
-subroutine wetdepokernel_nest &
-       (nunc,deposit,x,y,nage,kp)
-  !        i      i    i i  i    i
+subroutine wetdepokernel_nest(nunc,deposit,x,y,nage,kp)
+  !                           i    i       i i i    i
   !*****************************************************************************
   !                                                                            *
@@ -55,12 +54,25 @@
   integer :: ix,jy,ixp,jyp,ks,kp,nunc,nage
 
-
+  real :: dbg_dx, dbg_dy, dbg_xoutshiftn, dbg_youtshiftn, dbg_dxoutn, dbg_dyoutn,dbg_t
 
   xl=(x*dx+xoutshiftn)/dxoutn
   yl=(y*dy+youtshiftn)/dyoutn
-  ix=int(xl)
-  jy=int(yl)
+
+  ! old:
+  ! ix=int(xl) 
+  ! jy=int(yl)
+
+  ! ESO: for xl,yl in range <-.5,-1> we get ix,jy=0 and thus negative
+  ! wx,wy as the int function rounds upwards for negative numbers.
+  ! Either use the floor function, or (perhaps more correctly?) use "(xl.gt.-0.5)" 
+  ! in place of "(ix.ge.0)" and similar for the upper boundary.
+
+  ! new:
+  ix=floor(xl)
+  jy=floor(yl)
+
   ddx=xl-real(ix)                   ! distance to left cell border
   ddy=yl-real(jy)                   ! distance to lower cell border
+
 
   if (ddx.gt.0.5) then
@@ -80,37 +92,35 @@
   endif
 
-
-  ! Determine mass fractions for four grid points
-  !**********************************************
+! Determine mass fractions for four grid points
+!**********************************************
 
   do ks=1,nspec
-
-  if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
-       (jy.le.numygridn-1)) then
-    w=wx*wy
+    if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
+         (jy.le.numygridn-1)) then
+      w=wx*wy
       wetgriduncn(ix,jy,ks,kp,nunc,nage)= &
            wetgriduncn(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w
-  endif
+    endif
 
-  if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgridn-1).and. &
-       (jyp.le.numygridn-1)) then
-    w=(1.-wx)*(1.-wy)
+    if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgridn-1).and. &
+         (jyp.le.numygridn-1)) then
+      w=(1.-wx)*(1.-wy)
       wetgriduncn(ixp,jyp,ks,kp,nunc,nage)= &
            wetgriduncn(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w
-  endif
+    endif
 
-  if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgridn-1).and. &
-       (jy.le.numygridn-1)) then
-    w=(1.-wx)*wy
+    if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgridn-1).and. &
+         (jy.le.numygridn-1)) then
+      w=(1.-wx)*wy
       wetgriduncn(ixp,jy,ks,kp,nunc,nage)= &
            wetgriduncn(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w
-  endif
+    endif
 
-  if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgridn-1).and. &
-       (jyp.le.numygridn-1)) then
-    w=wx*(1.-wy)
+    if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgridn-1).and. &
+         (jyp.le.numygridn-1)) then
+      w=wx*(1.-wy)
       wetgriduncn(ix,jyp,ks,kp,nunc,nage)= &
            wetgriduncn(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w
-  endif
+    endif
 
   end do
