Index: /branches/petra/src/FLEXPART.f90
===================================================================
--- /branches/petra/src/FLEXPART.f90	(revision 36)
+++ /branches/petra/src/FLEXPART.f90	(revision 37)
@@ -66,5 +66,5 @@
 
   ! FLEXPART version string
-  flexversion='Version 9.2_beta.r34 (2015-02-12)'
+  flexversion='Version 9.2_beta.r3? (2015-02-25)'
   verbosity=0
 
Index: /branches/petra/src/com_mod.f90
===================================================================
--- /branches/petra/src/com_mod.f90	(revision 36)
+++ /branches/petra/src/com_mod.f90	(revision 37)
@@ -32,4 +32,5 @@
 !        Modifications: 15 August 2013 IP, 
 !        2/2015 PS, add incremental deposition switch
+!        2/2015 PS, change variables for cloud diagnostics
 !                                                                              *
 !*******************************************************************************
@@ -365,8 +366,4 @@
   real :: qvh(0:nxmax-1,0:nymax-1,nuvzmax,2)
   real :: pplev(0:nxmax-1,0:nymax-1,nuvzmax,2)
-  integer(kind=1) :: clouds(0:nxmax-1,0:nymax-1,nzmax,2)
-  integer :: cloudsh(0:nxmax-1,0:nymax-1,2)
-
-
   ! uu,vv,ww [m/2]       wind components in x,y and z direction
   ! uupol,vvpol [m/s]    wind components in polar stereographic projection
@@ -377,8 +374,4 @@
   ! drhodz [kg/m2]       vertical air density gradient
   ! tth,qvh              tth,qvh on original eta levels
-  ! clouds:   no cloud, no precipitation   0
-  !      cloud, no precipitation      1
-  !      rainout  conv/lsp dominated  2/3
-  !      washout  conv/lsp dominated  4/5
 
   ! pplev for the GFS version
@@ -406,4 +399,5 @@
   real :: oli(0:nxmax-1,0:nymax-1,1,2)
   real :: diffk(0:nxmax-1,0:nymax-1,1,2)
+  integer, dimension (0:nxmax-1,0:nymax-1,2) :: icloudbot, icloudthck
 
   ! ps                   surface pressure
@@ -426,4 +420,6 @@
   ! oli [m]              inverse Obukhov length (1/L)
   ! diffk [m2/s]         diffusion coefficient at reference height
+  ! icloudbot (m)        cloud bottom height
+  ! icloudthck (m)       cloud thickness   
 
 
@@ -484,6 +480,4 @@
   real :: qvn(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
   real :: pvn(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
-  integer(kind=1) :: cloudsn(0:nxmaxn-1,0:nymaxn-1,0:nzmax,2,maxnests)
-  integer :: cloudsnh(0:nxmaxn-1,0:nymaxn-1,2,maxnests)
   real :: rhon(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
   real :: drhodzn(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
@@ -514,4 +508,5 @@
   real :: diffkn(0:nxmaxn-1,0:nymaxn-1,1,2,maxnests)
   real :: vdepn(0:nxmaxn-1,0:nymaxn-1,maxspec,2,maxnests)
+  integer, dimension (0:nxmax-1,0:nymax-1,2,maxnests) :: icloudbotn, icloudthckn
 
 
Index: /branches/petra/src/getfields.f90
===================================================================
--- /branches/petra/src/getfields.f90	(revision 36)
+++ /branches/petra/src/getfields.f90	(revision 37)
@@ -41,4 +41,7 @@
   !   Variables tth,qvh,tthn,qvhn (on eta coordinates) in common block.
   !   Function of nstop extended.
+  !
+  !  3/2015, PS: put in verbosity output
+  !
   !*****************************************************************************
   !                                                                            *
@@ -83,4 +86,5 @@
   integer :: indmin = 1
 
+  call verboutput(verbosity,1,'  getfields')
 
   ! Check, if wind fields are available for the current time step
@@ -125,9 +129,15 @@
     do indj=indmin,numbwf-1
        if (ldirect*wftime(indj+1).gt.ldirect*itime) then
+           call verboutput(verbosity,1,'    call readwind')
           call readwind(indj+1,memind(2),uuh,vvh,wwh)
+           call verboutput(verbosity,1,'    call readwind_nests')
           call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn)
+           call verboutput(verbosity,1,'    call calcpar')
           call calcpar(memind(2),uuh,vvh,pvh)
+           call verboutput(verbosity,1,'    call calcpar_nests')
           call calcpar_nests(memind(2),uuhn,vvhn,pvhn)
+           call verboutput(verbosity,1,'    call verttransform')
           call verttransform(memind(2),uuh,vvh,wwh,pvh)
+           call verboutput(verbosity,1,'    call verttransform_nests')
           call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
           memtime(2)=wftime(indj+1)
@@ -148,17 +158,29 @@
              (ldirect*wftime(indj+1).gt.ldirect*itime)) then
            memind(1)=1
+            call verboutput(verbosity,1,'    call readwind')
            call readwind(indj,memind(1),uuh,vvh,wwh)
+            call verboutput(verbosity,1,'    call readwind_nests')
            call readwind_nests(indj,memind(1),uuhn,vvhn,wwhn)
+            call verboutput(verbosity,1,'    call calcpar')
            call calcpar(memind(1),uuh,vvh,pvh)
+            call verboutput(verbosity,1,'    call calcpar_nests')
            call calcpar_nests(memind(1),uuhn,vvhn,pvhn)
+            call verboutput(verbosity,1,'    call verttransform')
            call verttransform(memind(1),uuh,vvh,wwh,pvh)
+            call verboutput(verbosity,1,'    call verttransform_nests')
            call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn)
            memtime(1)=wftime(indj)
            memind(2)=2
+            call verboutput(verbosity,1,'    call readwind')
            call readwind(indj+1,memind(2),uuh,vvh,wwh)
+            call verboutput(verbosity,1,'    call readwind_nests')
            call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn)
+            call verboutput(verbosity,1,'    call calcpar')
            call calcpar(memind(2),uuh,vvh,pvh)
+            call verboutput(verbosity,1,'    call calcpar_nests')
            call calcpar_nests(memind(2),uuhn,vvhn,pvhn)
+            call verboutput(verbosity,1,'    call verttransform')
            call verttransform(memind(2),uuh,vvh,wwh,pvh)
+            call verboutput(verbosity,1,'    call verttransform_nests')
            call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
            memtime(2)=wftime(indj+1)
@@ -175,3 +197,5 @@
   if (lwindinterv.gt.idiffmax) nstop=3
 
+  call verboutput(verbosity,1,'    leaving getfields')
+
 end subroutine getfields
Index: /branches/petra/src/gridcheck.f90
===================================================================
--- /branches/petra/src/gridcheck.f90	(revision 36)
+++ /branches/petra/src/gridcheck.f90	(revision 37)
@@ -32,4 +32,6 @@
   !             LAST UPDATE: 1997-10-10                                 *
   !                                                                     *
+  !**********************************************************************
+  !                                                                     *
   !             Update:      1999-02-08, global fields allowed, A. Stohl*
   !             CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with *
@@ -37,5 +39,8 @@
   !             CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
   !                                 ECMWF grib_api                      *
-  !                                                                     *
+  !
+  !   3/2015, PS: add GRIB2 coding for deta/dt dp/deta used by ZAMG
+  !  
+  !
   !**********************************************************************
   !                                                                     *
@@ -144,99 +149,103 @@
   if (gribVer.eq.1) then ! GRIB Edition 1
 
-  !print*,'GRiB Edition 1'
-  !read the grib2 identifiers
-  call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'level',isec1(8),iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-
-  !change code for etadot to code for omega
-  if (isec1(6).eq.77) then
-    isec1(6)=135
-  endif
-
-  !print*,isec1(6),isec1(8)
-
-  else
-
-  !print*,'GRiB Edition 2'
-  !read the grib2 identifiers
-  call grib_get_int(igrib,'discipline',discipl,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'parameterCategory',parCat,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'parameterNumber',parNum,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'level',valSurf,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'paramId',parId,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-
-  !print*,discipl,parCat,parNum,typSurf,valSurf
-
-  !convert to grib1 identifiers
-  isec1(6)=-1
-  isec1(7)=-1
-  isec1(8)=-1
-  isec1(8)=valSurf     ! level
-  if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
-    isec1(6)=130         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
-    isec1(6)=131         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
-    isec1(6)=132         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
-    isec1(6)=133         ! indicatorOfParameter
-  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
-    isec1(6)=135         ! indicatorOfParameter
-  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
-    isec1(6)=151         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
-    isec1(6)=165         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
-    isec1(6)=166         ! indicatorOfParameter
-  elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
-    isec1(6)=167         ! indicatorOfParameter
-  elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
-    isec1(6)=168         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
-    isec1(6)=141         ! indicatorOfParameter
-  elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
-    isec1(6)=164         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
-    isec1(6)=142         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
-    isec1(6)=143         ! indicatorOfParameter
-  elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
-    isec1(6)=146         ! indicatorOfParameter
-  elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
-    isec1(6)=176         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
-    isec1(6)=180         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
-    isec1(6)=181         ! indicatorOfParameter
-  elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
-    isec1(6)=129         ! indicatorOfParameter
-  elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
-    isec1(6)=160         ! indicatorOfParameter
-  elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
-       (typSurf.eq.1)) then ! LSM
-    isec1(6)=172         ! indicatorOfParameter
-  else
-    print*,'***ERROR: undefined GRiB2 message found!',discipl, &
-         parCat,parNum,typSurf
-  endif
-  if(parId .ne. isec1(6) .and. parId .ne. 77) then
-    write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
-!    stop
-  endif
-
-  endif
+    !print*,'GRiB Edition 1'
+    !read the grib2 identifiers
+    call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'level',isec1(8),iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+
+    !change code for etadot to code for omega
+    if (isec1(6).eq.77) then
+      isec1(6)=135
+    endif
+
+    !print*,isec1(6),isec1(8)
+
+  else ! GRIB 2
+
+    !print*,'GRiB Edition 2'
+    ! read the grib2 identifiers
+
+    call grib_get_int(igrib,'discipline',discipl,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'parameterCategory',parCat,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'parameterNumber',parNum,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'level',valSurf,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'paramId',parId,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+
+    !print*,discipl,parCat,parNum,typSurf,valSurf
+
+    ! convert to grib1 identifiers
+
+    isec1(6)=-1
+    isec1(7)=-1
+    isec1(8)=-1
+    isec1(8)=valSurf     ! level
+    if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
+      isec1(6)=130         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
+      isec1(6)=131         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
+      isec1(6)=132         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
+      isec1(6)=133         ! indicatorOfParameter
+    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 detadtdpdeta
+      isec1(6)=135         ! indicatorOfParameter
+    elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually detadtdpdeta
+      isec1(6)=135         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.8)) then ! Omega, actually detadtdpdeta
+      isec1(6)=135         ! indicatorOfParameter
+    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
+      isec1(6)=151         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
+      isec1(6)=165         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
+      isec1(6)=166         ! indicatorOfParameter
+    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
+      isec1(6)=167         ! indicatorOfParameter
+    elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
+      isec1(6)=168         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
+      isec1(6)=141         ! indicatorOfParameter
+    elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
+      isec1(6)=164         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
+      isec1(6)=142         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
+      isec1(6)=143         ! indicatorOfParameter
+    elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
+      isec1(6)=146         ! indicatorOfParameter
+    elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
+      isec1(6)=176         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
+      isec1(6)=180         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
+      isec1(6)=181         ! indicatorOfParameter
+    elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
+      isec1(6)=129         ! indicatorOfParameter
+    elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
+      isec1(6)=160         ! indicatorOfParameter
+    elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
+         (typSurf.eq.1)) then ! LSM
+      isec1(6)=172         ! indicatorOfParameter
+    else
+      print*,'***ERROR: undefined GRiB2 message found!',discipl, &
+           parCat,parNum,typSurf
+    endif
+    if(parId .ne. isec1(6) .and. parId .ne. 77) then
+      write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
+  !    stop
+    endif
+
+  endif ! GRIB 1 / GRIB 2 cases
 
   !get the size and data of the values array
Index: /branches/petra/src/interpol_rain.f90
===================================================================
--- /branches/petra/src/interpol_rain.f90	(revision 36)
+++ /branches/petra/src/interpol_rain.f90	(revision 37)
@@ -1,4 +1,4 @@
 !**********************************************************************
-! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
+! Copyright 1998-2015                                                 *
 ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
 ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
@@ -20,6 +20,7 @@
 !**********************************************************************
 
-subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, &
-       ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3)
+subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx,ny, &
+       memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3, &
+       icbot,ictop)
   !                          i   i   i    i    i     i   i
   !i    i    i  i    i     i      i      i     o     o     o
@@ -40,4 +41,8 @@
   !                                                                           *
   !     30 August 1996                                                        *
+  !****************************************************************************
+  ! Petra Seibert, 2011/2012-2015:
+  !  Add interpolation of cloud bottom and cloud thickness
+  !  for fix to SE's new wet scavenging scheme
   !                                                                           *
   !****************************************************************************
@@ -66,15 +71,18 @@
   !                                                                           *
   !****************************************************************************
+  
+  use com_mod, only: icloudbot, icloudthck
+  use par_mod, only: icmv
 
   implicit none
 
   integer :: nx,ny,nxmax,nymax,nzmax,memind(2),m,ix,jy,ixp,jyp
-  integer :: itime,itime1,itime2,level,indexh
+  integer :: itime,itime1,itime2,level,indexh,ip1,ip2,ip3,ip4 
+  integer :: icbot,icthck,ictop,ipsum
   real :: yy1(0:nxmax-1,0:nymax-1,nzmax,2)
   real :: yy2(0:nxmax-1,0:nymax-1,nzmax,2)
   real :: yy3(0:nxmax-1,0:nymax-1,nzmax,2)
-  real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
+  real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2),cbot(2),cthck(2)
   real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
-
 
 
@@ -112,5 +120,7 @@
   !***********************
 
+  memind_loop: &
   do m=1,2
+
     indexh=memind(m)
 
@@ -127,5 +137,45 @@
          + p3*yy3(ix ,jyp,level,indexh) &
          + p4*yy3(ixp,jyp,level,indexh)
-  end do
+
+
+! PS clouds:
+    ip1=1
+    ip2=1
+    ip3=1
+    ip4=1
+    if (icloudbot(ix ,jy ,indexh) .eq. icmv) ip1=0
+    if (icloudbot(ixp,jy ,indexh) .eq. icmv) ip2=0
+    if (icloudbot(ix ,jyp,indexh) .eq. icmv) ip3=0
+    if (icloudbot(ixp,jyp,indexh) .eq. icmv) ip4=0
+    ipsum = ip1+ip2+ip3+ip4
+    if (ipsum .eq. 0) then
+      cbot(m)=icmv
+    else
+      cbot(m)=(ip1*p1*icloudbot(ix ,jy ,indexh) &
+             + ip2*p2*icloudbot(ixp,jy ,indexh) &
+             + ip3*p3*icloudbot(ix ,jyp,indexh) &
+             + ip4*p4*icloudbot(ixp,jyp,indexh)) / ipsum
+    endif
+
+    ip1=1
+    ip2=1
+    ip3=1
+    ip4=1
+    if (icloudthck(ix ,jy ,indexh) .eq. icmv) ip1=0
+    if (icloudthck(ixp,jy ,indexh) .eq. icmv) ip2=0
+    if (icloudthck(ix ,jyp,indexh) .eq. icmv) ip3=0
+    if (icloudthck(ixp,jyp,indexh) .eq. icmv) ip4=0
+    ipsum = ip1+ip2+ip3+ip4
+    if (ipsum .eq. 0) then
+      cthck(m)=icmv
+    else
+      cthck(m)=(ip1*p1*icloudthck(ix ,jy ,indexh) &
+              + ip2*p2*icloudthck(ixp,jy ,indexh) &
+              + ip3*p3*icloudthck(ix ,jyp,indexh) &
+              + ip4*p4*icloudthck(ixp,jyp,indexh)) / ipsum
+    endif
+! PS end clouds
+
+  enddo memind_loop
 
 
@@ -142,4 +192,21 @@
   yint3=(y3(1)*dt2+y3(2)*dt1)/dt
 
+! PS clouds:
+  icbot = nint( (cbot(1)*dt2 + cbot(2)*dt1)/dt )
+  if (nint(cbot(1)) .eq. icmv) icbot=cbot(2) 
+  if (nint(cbot(2)) .eq. icmv) icbot=cbot(1) 
+
+  icthck = nint( (cthck(1)*dt2 + cthck(2)*dt1)/dt )
+  if (nint(cthck(1)) .eq. icmv) icthck=cthck(2) 
+  if (nint(cthck(2)) .eq. icmv) icthck=cthck(1) 
+
+  if (icbot .ne. icmv .and. icthck .ne. icmv) then
+    ictop = icbot + icthck ! convert cloud thickness to cloud top
+  else
+    icbot=icmv
+    ictop=icmv
+  endif
+! PS end clouds
+
 
 end subroutine interpol_rain
Index: /branches/petra/src/interpol_rain_nests.f90
===================================================================
--- /branches/petra/src/interpol_rain_nests.f90	(revision 36)
+++ /branches/petra/src/interpol_rain_nests.f90	(revision 37)
@@ -1,4 +1,4 @@
 !**********************************************************************
-! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
+! Copyright 1998-2015                                                 *
 ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
 ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
@@ -21,6 +21,7 @@
 
 subroutine interpol_rain_nests(yy1,yy2,yy3,nxmaxn,nymaxn,nzmax, &
-       maxnests,ngrid,nxn,nyn,memind,xt,yt,level,itime1,itime2,itime, &
-       yint1,yint2,yint3)
+       maxnests,ngrid,nxn,nyn, &
+       memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3, &
+       icbot,ictop)
   !                                i   i   i    i      i      i
   !   i       i    i   i    i    i  i    i     i      i      i
@@ -44,4 +45,9 @@
   !                                                                           *
   !     15 March 2000                                                         *
+  !                                                                           *
+  !****************************************************************************
+  ! Petra Seibert, 2011/2012-2015:
+  !  Add interpolation of cloud bottom and cloud thickness
+  !  for fix to SE's new wet scavenging scheme
   !                                                                           *
   !****************************************************************************
@@ -71,24 +77,25 @@
   !****************************************************************************
 
+  use com_mod, only: icloudbotn, icloudthckn
+  use par_mod, only: icmv
+
   implicit none
 
   integer :: maxnests,ngrid
-  integer :: nxn(maxnests),nyn(maxnests),nxmaxn,nymaxn,nzmax,memind(2)
-  integer :: m,ix,jy,ixp,jyp,itime,itime1,itime2,level,indexh
+  integer :: nxn(maxnests),nyn(maxnests),nxmaxn,nymaxn,nzmax,memind(2),m
+  integer :: ix,jy,ixp,jyp
+  integer :: itime,itime1,itime2,level,indexh,ip1,ip2,ip3,ip4
+  integer :: icbot,ictop,icthck,ipsum
   real :: yy1(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
   real :: yy2(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
   real :: yy3(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
-  real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
+  real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2),cbot(2),cthck(2)
   real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
-
-
 
   ! If point at border of grid -> small displacement into grid
   !***********************************************************
 
-  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
 
 
@@ -119,5 +126,7 @@
   !***********************
 
+  memind_loop: &
   do m=1,2
+
     indexh=memind(m)
 
@@ -134,5 +143,44 @@
          + p3*yy3(ix ,jyp,level,indexh,ngrid) &
          + p4*yy3(ixp,jyp,level,indexh,ngrid)
-  end do
+
+! PS clouds:
+    ip1=1
+    ip2=1
+    ip3=1
+    ip4=1
+    if (icloudbotn(ix ,jy ,indexh,ngrid) .eq. icmv) ip1=0
+    if (icloudbotn(ixp,jy ,indexh,ngrid) .eq. icmv) ip2=0
+    if (icloudbotn(ix ,jyp,indexh,ngrid) .eq. icmv) ip3=0
+    if (icloudbotn(ixp,jyp,indexh,ngrid) .eq. icmv) ip4=0
+    ipsum = ip1+ip2+ip3+ip4
+    if (ipsum .eq. 0) then
+      cbot(m)=icmv
+    else
+      cbot(m)=(ip1*p1*icloudbotn(ix ,jy ,indexh,ngrid) &
+             + ip2*p2*icloudbotn(ixp,jy ,indexh,ngrid) &
+             + ip3*p3*icloudbotn(ix ,jyp,indexh,ngrid) &
+             + ip4*p4*icloudbotn(ixp,jyp,indexh,ngrid)) / ipsum
+    endif
+
+    ip1=1
+    ip2=1
+    ip3=1
+    ip4=1
+    if (icloudthckn(ix ,jy ,indexh,ngrid) .eq. icmv) ip1=0
+    if (icloudthckn(ixp,jy ,indexh,ngrid) .eq. icmv) ip2=0
+    if (icloudthckn(ix ,jyp,indexh,ngrid) .eq. icmv) ip3=0
+    if (icloudthckn(ixp,jyp,indexh,ngrid) .eq. icmv) ip4=0
+    ipsum = ip1+ip2+ip3+ip4
+    if (ipsum .eq. 0) then
+      cthck(m)=icmv
+    else
+      cthck(m)=(ip1*p1*icloudthckn(ix ,jy ,indexh,ngrid) &
+              + ip2*p2*icloudthckn(ixp,jy ,indexh,ngrid) &
+              + ip3*p3*icloudthckn(ix ,jyp,indexh,ngrid) &
+              + ip4*p4*icloudthckn(ixp,jyp,indexh,ngrid)) / ipsum
+    endif
+! PS end clouds
+
+  enddo memind_loop
 
 
@@ -150,3 +198,21 @@
 
 
+! PS clouds:
+  icbot = nint( (cbot(1)*dt2 + cbot(2)*dt1)/dt )
+  if (nint(cbot(1)) .eq. icmv) icbot=cbot(2) 
+  if (nint(cbot(2)) .eq. icmv) icbot=cbot(1) 
+
+  icthck = nint( (cthck(1)*dt2 + cthck(2)*dt1)/dt )
+  if (nint(cthck(1)) .eq. icmv) icthck=cthck(2) 
+  if (nint(cthck(2)) .eq. icmv) icthck=cthck(1) 
+
+  if (icbot .ne. icmv .and. icthck .ne. icmv) then
+    ictop = icbot + icthck ! convert cloud thickness to cloud top
+  else
+    icbot=icmv
+    ictop=icmv
+  endif
+! PS end clouds
+
+
 end subroutine interpol_rain_nests
Index: /branches/petra/src/par_mod.f90
===================================================================
--- /branches/petra/src/par_mod.f90	(revision 36)
+++ /branches/petra/src/par_mod.f90	(revision 37)
@@ -103,9 +103,10 @@
 
   integer,parameter :: idiffnorm=10800, idiffmax=2*idiffnorm, minstep=1
+  integer,parameter :: itagekernmin=10800
 
   ! idiffnorm [s]           normal time interval between two wind fields
   ! idiffmax [s]            maximum time interval between two wind fields
   ! minstep [s]             minimum time step to be used within FLEXPART
-
+  ! itagekernmin [s]        minimum particle age [s] for using output kernel
 
   !*****************************************************************
@@ -265,10 +266,20 @@
 
 !******************************************************
-! integer code for missing values, used in wet scavenging (PS, 2012)
+! for use in wet scavenging (PS, 2012, 2015)
 !******************************************************
 
-    !  integer icmv
-    !  parameter(icmv=-9999)
-      integer,parameter ::  icmv=-9999
+  integer,parameter :: icmv=-9999 ! integer code for missing values 
+  ! some cloud heights used when trying to construct reasonable cloud:
+  integer,parameter :: icloudtopconvmin = 6000
+  integer,parameter :: icloudtopmin = 3000
+  integer,parameter :: icloudbot1 = 500, icloudtop1 =  8000
+  integer,parameter :: icloudbot2 = 0,   icloudtop2 = 10000
+  
+  real, parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagnostics
+  real, parameter :: rhmininit = 0.90 ! standard condition for presence of clouds
+! PS       note that original by Sabine Eckhart was 80%
+! PS       however, for T<-20 C we consider saturation over ice
+! PS       so I think 90% should be enough          
+  real, parameter :: rhminx = 0.30 ! extreme condition for presence of clouds
 
 ! Parameters for testing
Index: /branches/petra/src/readcommand.f90
===================================================================
--- /branches/petra/src/readcommand.f90	(revision 36)
+++ /branches/petra/src/readcommand.f90	(revision 37)
@@ -235,13 +235,16 @@
     read(unitcommand,*,iostat=icmdstat) linit_cond
     if (icmdstat .gt. 0) &
-      print*, 'readcommand: linit_cond not read properly',icmdstat,linit_cond
-    if (old) call skplin(3,unitcommand)
-    read(unitcommand,*,iostat=icmdstat) surf_only
+      print*, 'readcommand: linit_cond not present or no read properly,&
+       & - check values: ',icmdstat,linit_cond
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*,iostat=icmdstat) icmdstat,surf_only
     if (icmdstat .gt. 0) &
-      print*, 'readcommand: linit_cond not read properly',icmdstat,surf_only
-    if (old) call skplin(3,unitcommand)
-    read(unitcommand,*,iostat=icmdstat) ldep_incr
+      print*, 'readcommand: surf_only not present or not read properly&
+       & - check values: ',icmdstat,surf_only
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*,iostat=icmdstat) icmdstat,ldep_incr
     if (icmdstat .gt. 0) &
-      print*, 'readcommand: linit_cond not read properly',icmdstat, ldep_incr
+      print*, 'readcommand: ldep_incr not present or not read properly&
+       & - check values: ',icmdstat, ldep_incr
     close(unitcommand)
 
@@ -416,17 +419,5 @@
 
 
-  ! For backward runs one releasefield for all releases makes no sense,
-  ! For quasilag and domainfill ioutputforechrelease is forbidden
-  !*****************************************************************************
-
-  if ((ldirect .lt. 0) .and. (ioutputforeachrelease .eq. 0)) then
-      write(*,*) '#### FLEXPART MODEL ERROR! FILE COMMAND:     ####'
-      write(*,*) '#### FOR BACKWARD RUNS, IOUTPUTFOREACHRLEASE ####'
-      write(*,*) '#### MUST BE SET TO ONE!                     ####'
-      stop
-  endif
-
-
-  ! For backward runs one releasefield for all releases makes no sense,
+  ! For domainfill runs one releasefield for all releases makes no sense,
   ! and is "forbidden"
   !*****************************************************************************
Index: /branches/petra/src/readspecies.f90
===================================================================
--- /branches/petra/src/readspecies.f90	(revision 36)
+++ /branches/petra/src/readspecies.f90	(revision 37)
@@ -20,10 +20,10 @@
 !**********************************************************************
 
-subroutine readspecies(id_spec,pos_spec)
+subroutine readspecies(id_spec,id_pos)
 
   !*****************************************************************************
   !                                                                            *
   !     This routine reads names and physical constants of chemical species/   *
-  !     radionuclides given in the parameter pos_spec                          *
+  !     radionuclides given in the parameter id_pos                          *
   !                                                                            *
   !   Author: A. Stohl                                                         *
@@ -32,8 +32,10 @@
   !
   !   Changes:                                                                 *
+  !                                                                            *
   !   N. Kristiansen, 31.01.2013: Including parameters for in-cloud scavenging *
   !                                                                            *
   !   HSO, 13 August 2013: added optional namelist input
   !   PS, 2/2015: access= -> position=
+  !   PS, 4/2015: quick fix for bug associated with end=22
   !                                                                            *
   !*****************************************************************************
@@ -60,5 +62,5 @@
   implicit none
 
-  integer :: i, pos_spec,j
+  integer :: i, id_pos,j
   integer :: idow,ihour,id_spec
   character(len=3) :: aspecnumb
@@ -69,5 +71,5 @@
   real :: pdsigma, pdryvel, pweightmolar, pohreact, pspec_ass, pkao
   real :: pweta_in, pwetb_in, pwetc_in, pwetd_in
-  integer :: readerror
+  integer :: ios
 
   ! declare namelist
@@ -100,122 +102,117 @@
   ! Open the SPECIES file and read species names and properties
   !************************************************************
-  specnum(pos_spec)=id_spec
-  write(aspecnumb,'(i3.3)') specnum(pos_spec)
+  specnum(id_pos)=id_spec
+  write(aspecnumb,'(i3.3)') specnum(id_pos)
   open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',form='formatted',err=998)
-  !write(*,*) 'reading SPECIES',specnum(pos_spec)
+  !write(*,*) 'reading SPECIES',specnum(id_pos)
 
   ASSSPEC=.FALSE.
 
   ! try namelist input
-  read(unitspecies,species_params,iostat=readerror)
+  read(unitspecies,species_params,iostat=ios)
   close(unitspecies)
    
-  if ((pweightmolar.eq.-789.0).or.(readerror.ne.0)) then ! no namelist found
-
-    readerror=1
-
-    open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',err=998)
+  if ((pweightmolar.eq.-789.0).or.(ios.ne.0)) then ! no namelist found
+
+    open(unitspecies,file=path(1) (1:length(1))//'SPECIES/SPECIES_'//aspecnumb, status='old', err=998)
 
     do i=1,6
       read(unitspecies,*)
     end do
-
-    read(unitspecies,'(a10)',end=22) species(pos_spec)
-  !  write(*,*) species(pos_spec)
-    read(unitspecies,'(f18.1)',end=22) decay(pos_spec)
-  !  write(*,*) decay(pos_spec)
-    read(unitspecies,'(e18.1)',end=22) weta(pos_spec)
-  !  write(*,*) weta(pos_spec)
-    read(unitspecies,'(f18.2)',end=22) wetb(pos_spec)
-  !  write(*,*) wetb(pos_spec)
+    read(unitspecies,'(a10)',end=20) species(id_pos)
+  !  write(*,*) species(id_pos)
+    read(unitspecies,'(f18.1)',end=20) decay(id_pos)
+  !  write(*,*) decay(id_pos)
+    read(unitspecies,'(e18.1)',end=20) weta(id_pos)
+  !  write(*,*) weta(id_pos)
+    read(unitspecies,'(f18.2)',end=20) wetb(id_pos)
+  !  write(*,*) wetb(id_pos)
 
   !*** NIK 31.01.2013: including in-cloud scavening parameters
-   read(unitspecies,'(e18.1)',end=22) weta_in(pos_spec)
-  !  write(*,*) weta_in(pos_spec)
-   read(unitspecies,'(f18.2)',end=22) wetb_in(pos_spec)
-  !  write(*,*) wetb_in(pos_spec)
-   read(unitspecies,'(f18.2)',end=22) wetc_in(pos_spec)
-  !  write(*,*) wetc_in(pos_spec)
-   read(unitspecies,'(f18.2)',end=22) wetd_in(pos_spec)
-  !  write(*,*) wetd_in(pos_spec)
-
-    read(unitspecies,'(f18.1)',end=22) reldiff(pos_spec)
-  !  write(*,*) reldiff(pos_spec)
-    read(unitspecies,'(e18.1)',end=22) henry(pos_spec)
-  !  write(*,*) henry(pos_spec)
-    read(unitspecies,'(f18.1)',end=22) f0(pos_spec)
-  !  write(*,*) f0(pos_spec)
-    read(unitspecies,'(e18.1)',end=22) density(pos_spec)
-  !  write(*,*) density(pos_spec)
-    read(unitspecies,'(e18.1)',end=22) dquer(pos_spec)
-  !  write(*,*) dquer(pos_spec)
-    read(unitspecies,'(e18.1)',end=22) dsigma(pos_spec)
-  !  write(*,*) dsigma(pos_spec)
-    read(unitspecies,'(f18.2)',end=22) dryvel(pos_spec)
-  !  write(*,*) dryvel(pos_spec)
-    read(unitspecies,'(f18.2)',end=22) weightmolar(pos_spec)
-  !  write(*,*) weightmolar(pos_spec)
-    read(unitspecies,'(e18.1)',end=22) ohreact(pos_spec)
-  !  write(*,*) ohreact(pos_spec)
-    read(unitspecies,'(i18)',end=22) spec_ass(pos_spec)
-  !  write(*,*) spec_ass(pos_spec)
-    read(unitspecies,'(f18.2)',end=22) kao(pos_spec)
-  !       write(*,*) kao(pos_spec)
-
-    pspecies=species(pos_spec)
-    pdecay=decay(pos_spec)
-    pweta=weta(pos_spec)
-    pwetb=wetb(pos_spec)
-    pweta_in=weta_in(pos_spec)
-    pwetb_in=wetb_in(pos_spec)
-    pwetc_in=wetc_in(pos_spec)
-    pwetd_in=wetd_in(pos_spec)
-    preldiff=reldiff(pos_spec)
-    phenry=henry(pos_spec)
-    pf0=f0(pos_spec)
-    pdensity=density(pos_spec)
-    pdquer=dquer(pos_spec)
-    pdsigma=dsigma(pos_spec)
-    pdryvel=dryvel(pos_spec)
-    pweightmolar=weightmolar(pos_spec)
-    pohreact=ohreact(pos_spec)
-    pspec_ass=spec_ass(pos_spec)
-    pkao=kao(pos_spec)
+   read(unitspecies,'(e18.1)',end=20) weta_in(id_pos)
+  !  write(*,*) weta_in(id_pos)
+   read(unitspecies,'(f18.2)',end=20) wetb_in(id_pos)
+  !  write(*,*) wetb_in(id_pos)
+   read(unitspecies,'(f18.2)',end=20) wetc_in(id_pos)
+  !  write(*,*) wetc_in(id_pos)
+   read(unitspecies,'(f18.2)',end=20) wetd_in(id_pos)
+  !  write(*,*) wetd_in(id_pos)
+
+    read(unitspecies,'(f18.1)',end=20) reldiff(id_pos)
+  !  write(*,*) reldiff(id_pos)
+    read(unitspecies,'(e18.1)',end=20) henry(id_pos)
+  !  write(*,*) henry(id_pos)
+    read(unitspecies,'(f18.1)',end=20) f0(id_pos)
+  !  write(*,*) f0(id_pos)
+    read(unitspecies,'(e18.1)',end=20) density(id_pos)
+  !  write(*,*) density(id_pos)
+    read(unitspecies,'(e18.1)',end=20) dquer(id_pos)
+  !  write(*,*) dquer(id_pos)
+    read(unitspecies,'(e18.1)',end=20) dsigma(id_pos)
+  !  write(*,*) dsigma(id_pos)
+    read(unitspecies,'(f18.2)',end=20) dryvel(id_pos)
+  !  write(*,*) dryvel(id_pos)
+    read(unitspecies,'(f18.2)',end=20) weightmolar(id_pos)
+  !  write(*,*) weightmolar(id_pos)
+    read(unitspecies,'(e18.1)',end=20) ohreact(id_pos)
+  !  write(*,*) ohreact(id_pos)
+    read(unitspecies,'(i18)',end=20) spec_ass(id_pos)
+  !  write(*,*) spec_ass(id_pos)
+    read(unitspecies,'(f18.2)',end=20) kao(id_pos)
+  !       write(*,*) kao(id_pos)
+
+    pspecies=species(id_pos)
+    pdecay=decay(id_pos)
+    pweta=weta(id_pos)
+    pwetb=wetb(id_pos)
+    pweta_in=weta_in(id_pos)
+    pwetb_in=wetb_in(id_pos)
+    pwetc_in=wetc_in(id_pos)
+    pwetd_in=wetd_in(id_pos)
+    preldiff=reldiff(id_pos)
+    phenry=henry(id_pos)
+    pf0=f0(id_pos)
+    pdensity=density(id_pos)
+    pdquer=dquer(id_pos)
+    pdsigma=dsigma(id_pos)
+    pdryvel=dryvel(id_pos)
+    pweightmolar=weightmolar(id_pos)
+    pohreact=ohreact(id_pos)
+    pspec_ass=spec_ass(id_pos)
+    pkao=kao(id_pos)
 
   else
 
-    species(pos_spec)=pspecies
-    decay(pos_spec)=pdecay
-    weta(pos_spec)=pweta
-    wetb(pos_spec)=pwetb
-    weta_in(pos_spec)=pweta_in
-    wetb_in(pos_spec)=pwetb_in
-    wetc_in(pos_spec)=pwetc_in
-    wetd_in(pos_spec)=pwetd_in
-    reldiff(pos_spec)=preldiff
-    henry(pos_spec)=phenry
-    f0(pos_spec)=pf0
-    density(pos_spec)=pdensity
-    dquer(pos_spec)=pdquer
-    dsigma(pos_spec)=pdsigma
-    dryvel(pos_spec)=pdryvel
-    weightmolar(pos_spec)=pweightmolar
-    ohreact(pos_spec)=pohreact
-    spec_ass(pos_spec)=pspec_ass
-    kao(pos_spec)=pkao
-
-  endif
-
-  i=pos_spec
-
-  if ((weta(pos_spec).gt.0).and.(henry(pos_spec).le.0)) then
-   if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set
-  endif
-
-  if (spec_ass(pos_spec).gt.0) then
+    species(id_pos)=pspecies
+    decay(id_pos)=pdecay
+    weta(id_pos)=pweta
+    wetb(id_pos)=pwetb
+    weta_in(id_pos)=pweta_in
+    wetb_in(id_pos)=pwetb_in
+    wetc_in(id_pos)=pwetc_in
+    wetd_in(id_pos)=pwetd_in
+    reldiff(id_pos)=preldiff
+    henry(id_pos)=phenry
+    f0(id_pos)=pf0
+    density(id_pos)=pdensity
+    dquer(id_pos)=pdquer
+    dsigma(id_pos)=pdsigma
+    dryvel(id_pos)=pdryvel
+    weightmolar(id_pos)=pweightmolar
+    ohreact(id_pos)=pohreact
+    spec_ass(id_pos)=pspec_ass
+    kao(id_pos)=pkao
+
+  endif
+
+  if ((weta(id_pos).gt.0).and.(henry(id_pos).le.0)) then
+   if (dquer(id_pos).le.0) goto 996 ! no particle, no henry set
+  endif
+
+  if (spec_ass(id_pos).gt.0) then
     spec_found=.FALSE.
-    do j=1,pos_spec-1
-      if (spec_ass(pos_spec).eq.specnum(j)) then
-        spec_ass(pos_spec)=j
+    do j=1,id_pos-1
+      if (spec_ass(id_pos).eq.specnum(j)) then
+        spec_ass(id_pos)=j
         spec_found=.TRUE.
         ASSSPEC=.TRUE.
@@ -227,8 +224,8 @@
   endif
 
-  if (dsigma(i).eq.1.) dsigma(i)=1.0001   ! avoid floating exception
-  if (dsigma(i).eq.0.) dsigma(i)=1.0001   ! avoid floating exception
-
-  if ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then
+  if (dsigma(id_pos).eq.1.) dsigma(id_pos)=1.0001   ! avoid floating exception
+  if (dsigma(id_pos).eq.0.) dsigma(id_pos)=1.0001   ! avoid floating exception
+
+  if ((reldiff(id_pos).gt.0.).and.(density(id_pos).gt.0.)) then
     write(*,*) '#### FLEXPART MODEL ERROR! FILE "SPECIES"    ####'
     write(*,*) '#### IS CORRUPT. SPECIES CANNOT BE BOTH      ####'
@@ -246,21 +243,22 @@
 
   do j=1,24           ! initialize everything to no variation
-    area_hour(i,j)=1.
-    point_hour(i,j)=1.
+    area_hour(id_pos,j)=1.
+    point_hour(id_pos,j)=1.
   end do
   do j=1,7
-    area_dow(i,j)=1.
-    point_dow(i,j)=1.
+    area_dow(id_pos,j)=1.
+    point_dow(id_pos,j)=1.
   end do
 
-  if (readerror.ne.0) then ! text format input
-
-    read(unitspecies,*,end=22)
+  if (ios.ne.0) then ! text format input
+
+    read(unitspecies,*,iostat=ios)
+    if (ios .ne. 0) goto 22 ! ifort does not like err=
     do j=1,24     ! 24 hours, starting with 0-1 local time
-      read(unitspecies,*) ihour,area_hour(i,j),point_hour(i,j)
+      read(unitspecies,*) ihour,area_hour(id_pos,j),point_hour(id_pos,j)
     end do
     read(unitspecies,*)
     do j=1,7      ! 7 days of the week, starting with Monday
-      read(unitspecies,*) idow,area_dow(i,j),point_dow(i,j)
+      read(unitspecies,*) idow,area_dow(id_pos,j),point_dow(id_pos,j)
     end do
 
Index: /branches/petra/src/readwind.f90
===================================================================
--- /branches/petra/src/readwind.f90	(revision 36)
+++ /branches/petra/src/readwind.f90	(revision 37)
@@ -40,4 +40,6 @@
   !   Variables tth and qvh (on eta coordinates) in common block
   ! 2/2015, PS: add missing paramter iret in call to grib subr
+  ! 3/2015, PS: add some verbosity messages
+  ! 3/2015, PS: add GRIB2 coding for deta/dt dp/deta used by ZAMG
   !
   !**********************************************************************
@@ -113,5 +115,9 @@
   ! OPENING OF DATA FILE (GRIB CODE)
   !
-5   call grib_open_file(ifile,path(3)(1:length(3)) &
+  
+5 continue
+  call verboutput(verbosity,1,'      call grib_open_file '//path(3)(1:length(3))&
+         //trim(wfname(indj)))
+  call grib_open_file(ifile,path(3)(1:length(3)) &
          //trim(wfname(indj)),'r',iret)
   if (iret.ne.GRIB_SUCCESS) then
@@ -127,4 +133,5 @@
   ! GET NEXT FIELDS
   !
+  call verboutput(verbosity,2,'      call grib_new_from_file')
   call grib_new_from_file(ifile,igrib,iret)
   if (iret.eq.GRIB_END_OF_FILE)  then
@@ -135,107 +142,116 @@
 
   !first see if we read GRIB1 or GRIB2
+  call verboutput(verbosity,2,'      call grib_get_int')
   call grib_get_int(igrib,'editionNumber',gribVer,iret)
+  call verboutput(verbosity,2,'      call grib_check')
   call grib_check(iret,gribFunction,gribErrorMsg)
 
   if (gribVer.eq.1) then ! GRIB Edition 1
 
-  !print*,'GRiB Edition 1'
-  !read the grib2 identifiers
-  call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'level',isec1(8),iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-
-  !change code for etadot to code for omega
-  if (isec1(6).eq.77) then
-    isec1(6)=135
-  endif
-
-  conversion_factor=1.
-
-  else
-
-  !print*,'GRiB Edition 2'
-  !read the grib2 identifiers
-  call grib_get_int(igrib,'discipline',discipl,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'parameterCategory',parCat,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'parameterNumber',parNum,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'level',valSurf,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'paramId',parId,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-
-  !print*,discipl,parCat,parNum,typSurf,valSurf
-
-  !convert to grib1 identifiers
-  isec1(6)=-1
-  isec1(7)=-1
-  isec1(8)=-1
-  isec1(8)=valSurf     ! level
-  conversion_factor=1.
-  if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
-    isec1(6)=130         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
-    isec1(6)=131         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
-    isec1(6)=132         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
-    isec1(6)=133         ! indicatorOfParameter
-  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
-    isec1(6)=135         ! indicatorOfParameter
-  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
-    isec1(6)=151         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
-    isec1(6)=165         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
-    isec1(6)=166         ! indicatorOfParameter
-  elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
-    isec1(6)=167         ! indicatorOfParameter
-  elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
-    isec1(6)=168         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
-    isec1(6)=141         ! indicatorOfParameter
-    conversion_factor=1000.
-  elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
-    isec1(6)=164         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
-    isec1(6)=142         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
-    isec1(6)=143         ! indicatorOfParameter
-    conversion_factor=1000.
-  elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
-    isec1(6)=146         ! indicatorOfParameter
-  elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
-    isec1(6)=176         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
-    isec1(6)=180         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
-    isec1(6)=181         ! indicatorOfParameter
-  elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
-    isec1(6)=129         ! indicatorOfParameter
-  elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
-    isec1(6)=160         ! indicatorOfParameter
-  elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
-       (typSurf.eq.1)) then ! LSM
-    isec1(6)=172         ! indicatorOfParameter
-  else
-    print*,'***WARNING: undefined GRiB2 message found!',discipl, &
-         parCat,parNum,typSurf
-  endif
-  if(parId .ne. isec1(6) .and. parId .ne. 77) then
-    write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
-!    stop
-  endif
-
-  endif
+    !print*,'GRiB Edition 1'
+
+    ! read the grib2 identifiers
+
+    call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'level',isec1(8),iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+
+    !change code for etadot to code for omega
+    if (isec1(6).eq.77) then
+      isec1(6)=135
+    endif
+
+    conversion_factor=1.
+
+  else ! GRIB 2
+
+    !print*,'GRiB Edition 2'
+
+    ! read the grib2 identifiers
+
+    call grib_get_int(igrib,'discipline',discipl,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'parameterCategory',parCat,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'parameterNumber',parNum,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'level',valSurf,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'paramId',parId,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+
+    !print*,discipl,parCat,parNum,typSurf,valSurf
+
+    ! convert to grib1 identifiers
+
+    isec1(6)=-1
+    isec1(7)=-1
+    isec1(8)=-1
+    isec1(8)=valSurf     ! level
+    conversion_factor=1.
+    if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
+      isec1(6)=130         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
+      isec1(6)=131         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
+      isec1(6)=132         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
+      isec1(6)=133         ! indicatorOfParameter
+    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 detadtdpdeta
+      isec1(6)=135         ! indicatorOfParameter
+    elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually detadtdpdeta
+      isec1(6)=135         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.8)) then ! Omega, actually detadtdpdeta
+      isec1(6)=135         ! indicatorOfParameter
+    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
+      isec1(6)=151         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
+      isec1(6)=165         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
+      isec1(6)=166         ! indicatorOfParameter
+    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
+      isec1(6)=167         ! indicatorOfParameter
+    elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
+      isec1(6)=168         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
+      isec1(6)=141         ! indicatorOfParameter
+      conversion_factor=1000.
+    elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
+      isec1(6)=164         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
+      isec1(6)=142         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
+      isec1(6)=143         ! indicatorOfParameter
+      conversion_factor=1000.
+    elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
+      isec1(6)=146         ! indicatorOfParameter
+    elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
+      isec1(6)=176         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
+      isec1(6)=180         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
+      isec1(6)=181         ! indicatorOfParameter
+    elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
+      isec1(6)=129         ! indicatorOfParameter
+    elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
+      isec1(6)=160         ! indicatorOfParameter
+    elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
+         (typSurf.eq.1)) then ! LSM
+      isec1(6)=172         ! indicatorOfParameter
+    else
+      print*,'***WARNING: undefined GRiB2 message found!',discipl, &
+           parCat,parNum,typSurf
+    endif
+    if(parId .ne. isec1(6) .and. parId .ne. 77) then
+      write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
+  !    stop
+    endif
+
+  endif ! end GRIB 1 / GRIB 2 cases
 
   !HSO  get the size and data of the values array
@@ -351,4 +367,5 @@
   end do
 
+  call verboutput(verbosity,2,'      call grib_release')
   call grib_release(igrib)
   goto 10                      !! READ NEXT LEVEL OR PARAMETER
@@ -380,4 +397,5 @@
 
   if (xglobal) then
+    call verboutput(verbosity,2,'      call shift_field_0')
     call shift_field_0(ewss,nxfield,ny)
     call shift_field_0(nsss,nxfield,ny)
@@ -429,4 +447,5 @@
         ff10m= sqrt(u10(i,j,1,n)**2+v10(i,j,1,n)**2)
         fflev1=sqrt(uuh(i,j,3)**2+vvh(i,j,3)**2)
+        call verboutput(verbosity,2,'      call pbl_profile')
         call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, &
              tt2(i,j,1,n),tth(i,j,3,n),ff10m,fflev1, &
@@ -457,5 +476,7 @@
   if(iwmax.ne.nwz)    stop 'READWIND: NWZ NOT CONSISTENT'
 
+  call verboutput(verbosity,1,'      leaving readwind')
   return
+
 888   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
   write(*,*) ' #### ',wfname(indj),'                    #### '
Index: /branches/petra/src/readwind_nests.f90
===================================================================
--- /branches/petra/src/readwind_nests.f90	(revision 36)
+++ /branches/petra/src/readwind_nests.f90	(revision 37)
@@ -34,9 +34,13 @@
   !                                                                            *
   !*****************************************************************************
-  !  Changes, Bernd C. Krueger, Feb. 2001:                                     *
+  !  CHANGES:
+  !
+  !  2/2001, Bernc C. Krueger:                                              *
   !        Variables tthn and qvhn (on eta coordinates) in common block        *
-  !  CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with ECMWF grib_api    *
-  !  CHANGE: 03/12/2008, Harald Sodemann, update to f90 with ECMWF grib_api    *
-  ! 2/2015, PS: add missing paramter iret in call to grib subr
+  !  11/01/2008, Harald Sodemann, GRIB1/2 input with ECMWF grib_api            *
+  !  03/12/2008, Harald Sodemann, update to f90 with ECMWF grib_api            *
+  !  2/2015, PS: add missing paramter iret in call to grib subr
+  !  3/2015, PS: add verbosity output
+  !  3/2015, PS: add GRIB2 coding for deta/dt dp/deta used by ZAMG
   !
   !*****************************************************************************
@@ -84,4 +88,6 @@
   character(len=24) :: gribErrorMsg = 'Error reading grib file'
   character(len=20) :: gribFunction = 'readwind_nests'
+  
+  call verboutput(verbosity,1,'      entered readwind_nests')
 
   do l=1,numbnests
@@ -100,5 +106,8 @@
   !
 
-5   call grib_open_file(ifile,path(numpath+2*(l-1)+1) &
+5 continue
+  call verboutput(verbosity,1,'        call grib_open_file '//&
+    path(numpath+2*(l-1)+1)(1:length(numpath+2*(l-1)+1))//trim(wfnamen(l,indj)))
+  call grib_open_file(ifile,path(numpath+2*(l-1)+1) &
          (1:length(numpath+2*(l-1)+1))//trim(wfnamen(l,indj)),'r')
   if (iret.ne.GRIB_SUCCESS) then
@@ -127,101 +136,106 @@
   if (gribVer.eq.1) then ! GRIB Edition 1
 
-  !print*,'GRiB Edition 1'
-  !read the grib2 identifiers
-  call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'level',isec1(8),iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-
-  !change code for etadot to code for omega
-  if (isec1(6).eq.77) then
-    isec1(6)=135
-  endif
-
-  conversion_factor=1.
-
+    !print*,'GRiB Edition 1'
+    
+    ! read the grib2 identifiers
+    
+    call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'level',isec1(8),iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+
+    !change code for etadot to code for omega
+    if (isec1(6).eq.77) then
+      isec1(6)=135
+    endif
+
+    conversion_factor=1.
 
   else
 
-  !print*,'GRiB Edition 2'
-  !read the grib2 identifiers
-  call grib_get_int(igrib,'discipline',discipl,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'parameterCategory',parCat,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'parameterNumber',parNum,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'level',valSurf,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'paramId',parId,iret) !added by mc to make it consisitent with new readwind.f90
-  call grib_check(iret,gribFunction,gribErrorMsg) !added by mc to make it consisitent with new readwind.f90
-
-  !print*,discipl,parCat,parNum,typSurf,valSurf
-
-  !convert to grib1 identifiers
-  isec1(6)=-1
-  isec1(7)=-1
-  isec1(8)=-1
-  isec1(8)=valSurf     ! level
-   conversion_factor=1.
-  if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
-    isec1(6)=130         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
-    isec1(6)=131         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
-    isec1(6)=132         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
-    isec1(6)=133         ! indicatorOfParameter
-  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 by mc to make it consisitent with new readwind.f90
-    isec1(6)=135         ! indicatorOfParameter    !added by mc to make it consisitent with new readwind.f90
-  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
-    isec1(6)=151         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
-    isec1(6)=165         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
-    isec1(6)=166         ! indicatorOfParameter
-  elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
-    isec1(6)=167         ! indicatorOfParameter
-  elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
-    isec1(6)=168         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
-    isec1(6)=141         ! indicatorOfParameter
-    conversion_factor=1000. !added by mc to make it consisitent with new readwind.f90
-  elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC !added by mc to make it consisitent with new readwind.f90
-    isec1(6)=164         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP !added by mc to make it consisitent with new readwind.f90
-    isec1(6)=142         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
-    isec1(6)=143         ! indicatorOfParameter
-    conversion_factor=1000. !added by mc to make it consisitent with new readwind.f90
-  elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
-    isec1(6)=146         ! indicatorOfParameter
-  elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
-    isec1(6)=176         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS !added by mc to make it consisitent with new readwind.f90
-    isec1(6)=180         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS !added by mc to make it consisitent with new readwind.f90
-    isec1(6)=181         ! indicatorOfParameter
-  elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
-    isec1(6)=129         ! indicatorOfParameter
-   elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO !added by mc to make it consisitent with new readwind.f90
-    isec1(6)=160         ! indicatorOfParameter
-  elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
-       (typSurf.eq.1)) then ! LSM
-    isec1(6)=172         ! indicatorOfParameter
-  else
-    print*,'***WARNING: undefined GRiB2 message found!',discipl, &
-         parCat,parNum,typSurf
-  endif
-  if(parId .ne. isec1(6) .and. parId .ne. 77) then !added by mc to make it consisitent with new readwind.f90
-    write(*,*) 'parId',parId, 'isec1(6)',isec1(6)  !
-!    stop
-  endif
+    !print*,'GRiB Edition 2'
+
+    ! read the grib2 identifiers
+
+    call grib_get_int(igrib,'discipline',discipl,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'parameterCategory',parCat,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'parameterNumber',parNum,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'level',valSurf,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'paramId',parId,iret) !added by mc to make it consistent with new readwind.f90
+    call grib_check(iret,gribFunction,gribErrorMsg) !added by mc to make it consistent with new readwind.f90
+
+    !print*,discipl,parCat,parNum,typSurf,valSurf
+
+    !convert to grib1 identifiers
+    isec1(6)=-1
+    isec1(7)=-1
+    isec1(8)=-1
+    isec1(8)=valSurf     ! level
+     conversion_factor=1.
+    if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
+      isec1(6)=130         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
+      isec1(6)=131         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
+      isec1(6)=132         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
+      isec1(6)=133         ! indicatorOfParameter
+    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 detadtdpdeta
+      isec1(6)=135         ! indicatorOfParameter
+    elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually detadtdpdeta
+      isec1(6)=135         ! indicatorOfParameter !added by mc to make it consistent with new readwind.f90
+    elseif ((parCat.eq.2).and.(parNum.eq.8)) then ! Omega, actually detadtdpdeta
+      isec1(6)=135         ! indicatorOfParameter !added by mc to make it consistent with new readwind.f90
+    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
+      isec1(6)=151         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
+      isec1(6)=165         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
+      isec1(6)=166         ! indicatorOfParameter
+    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
+      isec1(6)=167         ! indicatorOfParameter
+    elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
+      isec1(6)=168         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
+      isec1(6)=141         ! indicatorOfParameter
+      conversion_factor=1000. !added by mc to make it consistent with new readwind.f90
+    elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC !added by mc to make it consistent with new readwind.f90
+      isec1(6)=164         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP !added by mc to make it consistent with new readwind.f90
+      isec1(6)=142         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
+      isec1(6)=143         ! indicatorOfParameter
+      conversion_factor=1000. !added by mc to make it consistent with new readwind.f90
+    elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
+      isec1(6)=146         ! indicatorOfParameter
+    elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
+      isec1(6)=176         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS !added by mc to make it consistent with new readwind.f90
+      isec1(6)=180         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS !added by mc to make it consistent with new readwind.f90
+      isec1(6)=181         ! indicatorOfParameter
+    elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
+      isec1(6)=129         ! indicatorOfParameter
+     elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO !added by mc to make it consistent with new readwind.f90
+      isec1(6)=160         ! indicatorOfParameter
+    elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
+         (typSurf.eq.1)) then ! LSM
+      isec1(6)=172         ! indicatorOfParameter
+    else
+      print*,'***WARNING: undefined GRiB2 message found!',discipl, &
+           parCat,parNum,typSurf
+    endif
+    if(parId .ne. isec1(6) .and. parId .ne. 77) then !added by mc to make it consistent with new readwind.f90
+      write(*,*) 'parId',parId, 'isec1(6)',isec1(6)  !
+  !    stop
+    endif
 
   endif
@@ -246,13 +260,13 @@
   ! CHECK GRID SPECIFICATIONS
   if(isec2(2).ne.nxn(l)) stop &
-  'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL'
+    'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL'
   if(isec2(3).ne.nyn(l)) stop &
-  'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL'
-  if(isec2(12)/2-1.ne.nlev_ec) stop 'READWIND: VERTICAL DISCRET&
-  IZATION NOT CONSISTENT FOR A NESTING LEVEL'
+    'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL'
+  if(isec2(12)/2-1.ne.nlev_ec) stop &
+    'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT FOR A NESTING LEVEL'
   endif ! ifield
 
   !HSO  get the second part of the grid dimensions only from GRiB1 messages
- if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then ! !added by mc to make it consisitent with new readwind.f90
+ if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then ! !added by mc to make it consistent with new readwind.f90
     call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
          xauxin,iret)
@@ -294,5 +308,5 @@
              zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
         if(isec1(6).eq.141) sdn(i,j,1,n,l)= &!! SNOW DEPTH
-             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consisitent with new readwind.f90!
+             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consistent with new readwind.f90!
         if(isec1(6).eq.151) msln(i,j,1,n,l)= &!! SEA LEVEL PRESS.
              zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
@@ -312,5 +326,5 @@
         endif
         if(isec1(6).eq.143) then                         !! CONVECTIVE PREC.
-          convprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consisitent with new readwind.f90
+          convprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consistent with new readwind.f90
           if (convprecn(i,j,1,n,l).lt.0.) convprecn(i,j,1,n,l)=0.
         endif
@@ -421,5 +435,8 @@
   end do
 
+  call verboutput(verbosity,1,'      leaving readwind_nests')
+
   return
+
 888   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
   write(*,*) ' #### ',wfnamen(l,indj),' FOR NESTING LEVEL  #### '
@@ -431,4 +448,5 @@
   write(*,*) ' #### ',wfnamen(l,indj),'                    #### '
   write(*,*) ' #### CANNOT BE OPENED FOR NESTING LEVEL ',l,'####'
+  call verboutput(verbosity,1,'      leaving readwind')
 
 end subroutine readwind_nests
Index: /branches/petra/src/timemanager.f90
===================================================================
--- /branches/petra/src/timemanager.f90	(revision 36)
+++ /branches/petra/src/timemanager.f90	(revision 37)
@@ -137,10 +137,8 @@
   !write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
   write(*,46) float(itime)/3600,itime,numpart
-  if (verbosity.gt.0) then
-    write (*,*) 'timemanager> starting simulation'
-    if (verbosity.gt.1) then
-      CALL SYSTEM_CLOCK(count_clock)
-      WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
-    endif     
+  call verboutput(verbosity,1,'timemanager> starting simulation')
+  if (verbosity.gt.1) then
+    CALL SYSTEM_CLOCK(count_clock)
+    WRITE(*,*) 'timemanager> SYSTEM CLOCK', (count_clock-count_clock0) / real(count_rate)
   endif
 
Index: /branches/petra/src/verboutput.f90
===================================================================
--- /branches/petra/src/verboutput.f90	(revision 37)
+++ /branches/petra/src/verboutput.f90	(revision 37)
@@ -0,0 +1,55 @@
+!**********************************************************************
+! Copyright 1998-2015                                                 *
+! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
+! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
+!                                                                     *
+! This file is part of FLEXPART.                                      *
+!                                                                     *
+! FLEXPART is free software: you can redistribute it and/or modify    *
+! it under the terms of the GNU General Public License as published by*
+! the Free Software Foundation, either version 3 of the License, or   *
+! (at your option) any later version.                                 *
+!                                                                     *
+! FLEXPART is distributed in the hope that it will be useful,         *
+! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
+! GNU General Public License for more details.                        *
+!                                                                     *
+! You should have received a copy of the GNU General Public License   *
+! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
+!**********************************************************************
+
+  subroutine verboutput (iverbosity, iverbositymin, infostring)
+  
+  !**********************************************************************
+  !                                                                     *
+  !             SUBROUTINE VERBOUTPUT                                   *
+  !                                                                     *
+  !**********************************************************************
+  !                                                                     *
+  !             AUTHOR:      PETRA sEIBERT                              *
+  !             DATE:        2015-03-07                                 *
+  !             LAST UPDATE: 2000-10-17, Andreas Stohl                  *
+  !                                                                     *
+  !**********************************************************************
+  !  Changes:
+  !
+  !**********************************************************************
+  !                                                                     *
+  ! DESCRIPTION:                                                        *
+  !                                                                     *
+  ! write debugging information to standard output
+  ! depending on the verbosity level
+  !
+  !**********************************************************************
+
+  implicit none
+  
+  integer :: iverbosity ! verbosity level of the current run
+  integer :: iverbositymin ! min. verbosity level required to produce output
+  character(*) :: infostring
+  
+  if (iverbosity .ge. iverbositymin) print*, '>'//infostring
+  
+  end subroutine verboutput
+  
Index: /branches/petra/src/verttransform.f90
===================================================================
--- /branches/petra/src/verttransform.f90	(revision 36)
+++ /branches/petra/src/verttransform.f90	(revision 37)
@@ -1,4 +1,4 @@
 !**********************************************************************
-! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
+! Copyright 1998-2015                                                 *
 ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
 ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
@@ -50,8 +50,9 @@
   ! added the variable cloud for use with scavenging - descr. in com_mod
   !*****************************************************************************
+  ! Petra Seibert, Feb 2015: Add quick fix from 2013 
+  !*****************************************************************************
   !                                                                            *
   ! Variables:                                                                 *
   ! nx,ny,nz                        field dimensions in x,y and z direction    *
-  ! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition           *
   ! uu(0:nxmax,0:nymax,nzmax,2)     wind components in x-direction [m/s]       *
   ! vv(0:nxmax,0:nymax,nzmax,2)     wind components in y-direction [m/s]       *
@@ -69,8 +70,8 @@
   implicit none
 
-  integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
-  integer :: rain_cloud_above,kz_inv
+  integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym,kz_inv
+  integer :: icloudtop
   real :: f_qvsat,pressure
-  real :: rh,lsp,convp
+  real :: rh,lsp,convp,prec,rhmin
   real :: rhoh(nuvzmax),pinmconv(nzmax)
   real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
@@ -85,5 +86,5 @@
   real,parameter :: const=r_air/ga
 
-  logical :: init = .true.
+  logical :: init = .true., lconvectprec, lsearch
 
 
@@ -477,43 +478,78 @@
 
 
-  !write (*,*) 'initializing clouds, n:',n,nymin1,nxmin1,nz
-  !   create a cloud and rainout/washout field, clouds occur where rh>80%
-  !   total cloudheight is stored at level 0
+  ! write (*,*) 'diagnosing clouds, n:',n,nymin1,nxmin1,nz
+  jy_loop: &
   do jy=0,nymin1
     do ix=0,nxmin1
-      rain_cloud_above=0
+    
       lsp=lsprec(ix,jy,1,n)
       convp=convprec(ix,jy,1,n)
-      cloudsh(ix,jy,n)=0
-      do kz_inv=1,nz-1
-         kz=nz-kz_inv+1
-         pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
-         rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
-         clouds(ix,jy,kz,n)=0
-         if (rh.gt.0.8) then ! in cloud
-            if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
-               rain_cloud_above=1
-               cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
-               height(kz)-height(kz-1)
-               if (lsp.ge.convp) then
-                  clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
-               else
-                  clouds(ix,jy,kz,n)=2 ! convp dominated rainout
-               endif
-            else ! no precipitation
-                  clouds(ix,jy,kz,n)=1 ! cloud
+      prec=lsp+convp
+      if (lsp .gt. convp) then !  prectype='lsp'
+        lconvectprec = .false.
+      else ! prectype='cp '
+        lconvectprec = .true.
+      endif
+      icloudbot(ix,jy,n)=icmv
+      icloudtop=icmv ! this is just a local variable
+      rhmin=rhmininit ! just initialise in a way that cond is true
+      lsearch=.true.
+
+      cloudsearch_loop: &
+      do while (rhmin .ge. rhminx .and. lsearch)
+        ! give up for < rhminx
+
+        kz_loop: &
+        do kz_inv=1,nz-1
+          kz=nz-kz_inv+1
+          pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
+          rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
+          if (rh .gt. rhmin) then 
+            if (icloudbot(ix,jy,n) .eq. icmv) then
+              icloudtop=nint(height(kz)) ! use int to save memory
             endif
-         else ! no cloud
-            if (rain_cloud_above.eq.1) then ! scavenging
-               if (lsp.ge.convp) then
-                  clouds(ix,jy,kz,n)=5 ! lsp dominated washout
-               else
-                  clouds(ix,jy,kz,n)=4 ! convp dominated washout
-               endif
-            endif
-         endif
-      end do
-    end do
-  end do
+            icloudbot(ix,jy,n)=nint(height(kz))
+          endif
+        end do kz_loop
+
+        ! PS try to get a cloud thicker than 50 m 
+        ! PS if there is at least precmin mm/h
+        if (prec .gt. precmin .and. &
+          ( icloudbot(ix,jy,n) .eq. icmv .or. &
+            icloudtop-icloudbot(ix,jy,n) .lt. 50)) then
+          rhmin = rhmin - 0.05
+        else
+          lsearch = .false.
+        endif
+        
+      enddo cloudsearch_loop
+
+      ! PS implement a rough fix for badly represented convection
+      ! PS is based on looking at a limited set of comparison data
+      if (lconvectprec .and. icloudtop .lt. icloudtopconvmin .or. &
+         icloudbot(ix,jy,n) .lt. icloudtopmin .and. prec .gt. precmin) then
+        if (convp .lt. 0.1) then
+          icloudbot(ix,jy,n) = icloudbot1
+          icloudtop =          icloudtop1
+        else
+          icloudbot(ix,jy,n) = icloudbot2
+          icloudtop =          icloudtop2
+        endif
+      endif
+
+      if (icloudtop .ne. icmv) then
+        icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n)
+      else
+        icloudthck(ix,jy,n) = icmv ! no cloud found whatsoever
+      endif
+
+      ! PS get rid of too thin clouds      
+      if (icloudthck(ix,jy,n) .lt. 50) then
+        icloudbot(ix,jy,n)=icmv
+        icloudthck(ix,jy,n)=icmv
+      endif
+
+    enddo
+  enddo jy_loop
 
 end subroutine verttransform
Index: /branches/petra/src/verttransform_nests.f90
===================================================================
--- /branches/petra/src/verttransform_nests.f90	(revision 36)
+++ /branches/petra/src/verttransform_nests.f90	(revision 37)
@@ -1,4 +1,4 @@
 !**********************************************************************
-! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
+! Copyright 1998-2015                                                 *
 ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
 ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
@@ -52,4 +52,6 @@
   ! add the variable cloud for use with scavenging - descr. in com_mod
   !*****************************************************************************
+  ! Petra Seibert, Feb 2015: Add quick fix from 2013 
+  !*****************************************************************************
   !                                                                            *
   ! Variables:                                                                 *
@@ -69,12 +71,13 @@
   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
+  integer :: ix,jy,kz,iz,n,l,kmin,kl,klp,ix1,jy1,ixp,jyp,kz_inv
+  integer :: icloudtop
+  real :: f_qvsat,pressure
+  real :: rh,lsp,convp,prec,rhmin
+  real :: rhoh(nuvzmax),pinmconv(nzmax)
+  real :: wzlev(nwzmax),uvwzlev(0:nxmaxn-1,0:nymaxn-1,nzmax)
+  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
   real :: dzdx,dzdy
-  real :: dzdx1,dzdx2,dzdy1,dzdy2
+  real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf
   real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
   real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
@@ -83,8 +86,10 @@
   real,parameter :: const=r_air/ga
 
+  logical :: lconvectprec, lsearch
 
   ! Loop over all nests
   !********************
 
+  numbnest_loop: &
   do l=1,numbnests
 
@@ -143,6 +148,6 @@
 
 
-  ! Levels, where u,v,t and q are given
-  !************************************
+  ! Levels where u,v,t and q are given
+  !***********************************
 
       uun(ix,jy,1,n,l)=uuhn(ix,jy,1,l)
@@ -194,6 +199,6 @@
 
 
-  ! Levels, where w is given
-  !*************************
+  ! Levels where w is given
+  !************************
 
       wwn(ix,jy,1,n,l)=wwhn(ix,jy,1,l)*pinmconv(1)
@@ -275,49 +280,81 @@
 
       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 (*,*) 'diagnosing nested clouds, n:',n
+  jy_loop: &
   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
+      prec=lsp+convp
+      if (lsp .gt. convp) then !  prectype='lsp'
+        lconvectprec = .false.
+      else ! prectype='cp '
+        lconvectprec = .true.
+      endif
+      icloudbotn(ix,jy,n,l)=icmv
+      icloudtop=icmv ! this is just a local variable
+      rhmin=rhmininit  
+      lsearch=.true.
+
+      cloudsearch_loop: &
+      do while (rhmin .ge. rhminx .and. lsearch)  
+        ! give up for < rhminx
+
+        kz_loop: &
+        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))
+           if (rh .gt. rhmin) then 
+            if (icloudbotn(ix,jy,n,l) .eq. icmv) then
+              icloudtop=nint(height(kz)) ! use int to save memory
             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
+            icloudbotn(ix,jy,n,l)=nint(height(kz))
+          endif
+        end do kz_loop
+
+        ! PS try to get a cloud thicker than 50 m 
+        ! PS if there is at least precmin mm/h
+        if (prec .gt. precmin .and. &
+          ( icloudbotn(ix,jy,n,l) .eq. icmv .or. &
+            icloudtop-icloudbotn(ix,jy,n,l) .lt. 50)) then
+          rhmin = rhmin - 0.05
+        else
+          lsearch = .false.
+        endif
+        
+      enddo cloudsearch_loop
+      
+      ! PS implement a rough fix for badly represented convection
+      ! PS is based on looking at a limited set of comparison data
+      if (lconvectprec .and. icloudtop .lt. icloudtopconvmin .or. &
+         icloudbotn(ix,jy,n,l) .lt. icloudtopmin .and. prec .gt. precmin) then
+        if (convp .lt. 0.1) then
+          icloudbotn(ix,jy,n,l) = icloudbot1
+          icloudtop =             icloudtop1
+        else
+          icloudbotn(ix,jy,n,l) = icloudbot2
+          icloudtop =             icloudtop2
+        endif
+      endif
+      if (icloudtop .ne. icmv) then
+        icloudthckn(ix,jy,n,l) = icloudtop-icloudbotn(ix,jy,n,l)
+      else
+        icloudthckn(ix,jy,n,l) = icmv
+      endif
+      ! PS get rid of too thin clouds      
+      if (icloudthck(ix,jy,n) .lt. 50) then
+        icloudbotn(ix,jy,n,l)=icmv
+        icloudthckn(ix,jy,n,l)=icmv
+      endif
+
+    enddo
+  enddo jy_loop
+ 
+  enddo numbnest_loop
 
 end subroutine verttransform_nests
Index: /branches/petra/src/wetdepo.f90
===================================================================
--- /branches/petra/src/wetdepo.f90	(revision 36)
+++ /branches/petra/src/wetdepo.f90	(revision 37)
@@ -1,4 +1,4 @@
 !**********************************************************************
-! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
+! Copyright 1998-2015                                                 *
 ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
 ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
@@ -39,4 +39,9 @@
   ! Code may not be correct for decay of deposition!                           *
   !                                                                            *
+  ! PS, 2/2015: implement wet depo quick fix from 2011/12 in this version
+  !  it implements interpolated and improved clouds
+  !  Also, certain deficiencies for thin clouds fixed also
+  !  Pass itage to wetdepokernel
+  !                                                                            *
   !*****************************************************************************
   !                                                                            *
@@ -71,9 +76,8 @@
 
   integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
-  integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v
-  integer :: ks, kp
+  integer :: ngrid,itage,nage,kz,il,interp_time,n
+  integer :: ks, kp, n1,n2, icbot,ictop, indcloud
   real :: S_i, act_temp, cl, cle ! in cloud scavenging
-  real :: clouds_h ! cloud height for the specific grid point
-  real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav
+  real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav,wetscavold,precsub,f
   real :: wetdeposit(maxspec),restmass
   real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
@@ -97,12 +101,14 @@
   !************************
 
+  particle_loop: &
   do jpart=1,numpart
 
     if (itra1(jpart).eq.-999999999) goto 20
-    if(ldirect.eq.1)then
+    if (ldirect.eq.1) then
       if (itra1(jpart).gt.itime) goto 20
     else
       if (itra1(jpart).lt.itime) goto 20
     endif
+    
   ! Determine age class of the particle
     itage=abs(itra1(jpart)-itramem(jpart))
@@ -123,5 +129,5 @@
         goto 23
       endif
-    end do
+    enddo
 23  continue
 
@@ -131,6 +137,6 @@
 
     if (ngrid.gt.0) then
-      xtn=(xtra1(jpart)-xln(ngrid))*xresoln(ngrid)
-      ytn=(ytra1(jpart)-yln(ngrid))*yresoln(ngrid)
+      xtn=real(xtra1(jpart)-xln(ngrid))*xresoln(ngrid)
+      ytn=real(ytra1(jpart)-yln(ngrid))*yresoln(ngrid)
       ix=int(xtn)
       jy=int(ytn)
@@ -148,41 +154,54 @@
 
     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)
+      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,icbot,ictop)
     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)
-    endif
-
-    if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
+      nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn, &
+      1,memtime(1),memtime(2),interp_time,lsp,convp,cc,icbot,ictop)
+    endif
+
+! PS 2012/2015: subtract a small value, eg 0.01 mm/h, 
+! to remove spurious precip; replaces previous code
+        prec = lsp+convp
+        precsub = 0.01
+        if (prec .lt. precsub) then
+          goto 20
+        else
+          f = (prec-precsub)/prec
+          lsp = f*lsp
+          convp = f*convp
+        endif
+
   ! get the level were the actual particle is in
     do il=2,nz
       if (height(il).gt.ztra1(jpart)) then
-        hz=il-1
+        kz=il-1
         goto 26
       endif
-    end do
+    enddo
 26  continue
 
     n=memind(2)
     if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
-    n=memind(1)
+      n=memind(1)
 
   ! if there is no precipitation or the particle is above the clouds no
   ! scavenging is done
-
-    if (ngrid.eq.0) then
-      clouds_v=clouds(ix,jy,hz,n)
-      clouds_h=cloudsh(ix,jy,n)
-    else
-      clouds_v=cloudsn(ix,jy,hz,n,ngrid)
-      clouds_h=cloudsnh(ix,jy,n,ngrid)
-    endif
-  !write(*,*) 'there is
-  !    + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz
-    if (clouds_v.le.1) goto 20
-  !write (*,*) 'there is scavenging'
+  ! PS: part of 2011/2012/2015 fix, replaces previous code
+  
+    if (ztra1(jpart) .le. float(ictop)) then
+      if (ztra1(jpart) .gt. float(icbot)) then
+        indcloud = 2 ! in-cloud
+      else
+        indcloud = 1 ! below-cloud
+      endif
+    elseif (ictop .eq. icmv) then
+      indcloud = 0 ! no cloud found, use old scheme
+    else
+      goto 20 ! above cloud
+    endif
+
 
   ! 1) Parameterization of the the area fraction of the grid cell where the
@@ -228,30 +247,26 @@
   !**********************************************************
 
+    species_loop: &
     do ks=1,nspec      ! loop over species
+    
       wetdeposit(ks)=0. 
       wetscav=0.   
 
   ! NIK 09.12.13: allowed to turn off either below cloud or in-cloud by including TWO separate blocks of if tests
-
-
-  !****i*******************
-  ! BELOW CLOUD SCAVENGING
-  !***********************  
-      if (weta(ks).gt.0. .and. clouds_v.ge.4) then !if positive below-cloud coefficient given in SPECIES file, and if in below cloud regime
- !       for aerosols and not highliy soluble substances weta=5E-6
-        wetscav=weta(ks)*prec**wetb(ks)  ! scavenging coefficient 
-      endif !end below-cloud
-
-
-  !********************
-  ! IN CLOUD SCAVENGING
-  !********************
-      if (weta_in(ks).gt.0. .and. clouds_v.lt.4) then !if positive in-cloud coefficient given in SPECIES file, and if in in-cloud regime
-
-        if (ngrid.gt.0) then
-          act_temp=ttn(ix,jy,hz,n,ngrid)
-        else
-          act_temp=tt(ix,jy,hz,n)
-        endif
+      if (weta(ks).gt.0.) then 
+        ! positive below-cloud coefficient from SPECIES file
+        
+        if (indcloud .eq. 1) then ! below-cloud scavenging
+
+         ! for aerosols and not highly soluble substances weta=5E-6
+           wetscav=weta(ks)*prec**wetb(ks)  ! scavenging coefficient 
+
+        elseif (indcloud .eq. 2) then ! in-cloud scavenging
+
+          if (ngrid.gt.0) then
+            act_temp=ttn(ix,jy,kz,n,ngrid)
+          else
+            act_temp=tt(ix,jy,kz,n)
+          endif
 
 ! NIK 31.01.2013: SPECIES defined parameters for the in-cloud scavening 
@@ -259,34 +274,27 @@
 ! wetb_in=0.36 (default)
 ! wetc_in=0.9 (default)
-! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0-no scaling)
+! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0 - no scaling)
         cl=weta_in(ks)*prec**wetb_in(ks)
         if (dquer(ks).gt.0) then ! is particle
           S_i=wetc_in(ks)/cl
         else ! is gas
-          cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
-          S_i=1/cle
+          cle=(1.-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
+          S_i=1./cle
         endif
-        wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks)
-  !     write(*,*) 'in. wetscav:'
-  !     +  ,wetscav,cle,cl,act_temp,prec,clouds_h
-       
-      endif ! end in-cloud
-
-
+        wetscav=S_i*prec/3.6E6/wetd_in(ks)/(ictop-icbot) ! 3.6e6 converts mm/h to m/s
+! PS - prevent extrem values of wetscav:
+        wetscavold = 2.e-5*prec**0.8      
+        wetscav = min( 0.1*wetscav, wetscavold ) ! here we sacle the above wetscav. This is now twice, also in wetd_in.
+
+      else ! PS: no cloud diagnosed, old scheme, 
+
+! PS    using with fixed a,b for simplicity, one may wish to change!!
+        wetscav = 2.e-5*prec**0.8 ! was before: 1.e-4*prec**0.62
+
+      endif ! end regime cases
       
-  !**************************************************
-  ! 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
-        wetdeposit(ks)=xmass1(jpart,ks)* &
-        (1.-exp(-wetscav*abs(ltsample)))*grfraction  ! wet deposition
-      else ! if no scavenging
-        wetdeposit(ks)=0.
-      endif
-
-
+  ! calculate deposition 
+      wetdeposit(ks)=xmass1(jpart,ks)* &
+        (1.-exp(-wetscav*abs(ltsample)))*grfraction 
       restmass = xmass1(jpart,ks)-wetdeposit(ks)
       if (ioutputforeachrelease.eq.1) then
@@ -297,7 +305,6 @@
       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)
       else
         xmass1(jpart,ks)=0.
@@ -305,25 +312,28 @@
   !   Correct deposited mass to the last time step when radioactive decay of
   !   gridded deposited mass was calculated
-      if (decay(ks).gt.0.) then
+      if (decay(ks).gt.0.) &
         wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
-      endif
-
-
-
-    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
-  !*****************************************************************************
+
+    else ! weta(k) <= 0
+
+      wetdeposit(ks)=0.
+    
+    endif
+
+  enddo species_loop 
+
+ ! Sabine Eckhardt, June 2008: write deposition only in forward runs
+ ! add the wet deposition from this step to accumulated amount 
+ ! on output grid and nested output grid
 
     if (ldirect.eq.1) then
-      call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
-      real(ytra1(jpart)),nage,kp)
+      call wetdepokernel(nclass(jpart),wetdeposit, &
+        real(xtra1(jpart)),real(ytra1(jpart)),itage,nage,kp)
       if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
-        wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),nage,kp)
-    endif
-
-20  continue
-  end do ! all particles
+        wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),itage,nage,kp)
+    endif
+
+20  continue ! jump here for particles not to be treated
+  enddo particle_loop
 
 end subroutine wetdepo
Index: /branches/petra/src/wetdepokernel.f90
===================================================================
--- /branches/petra/src/wetdepokernel.f90	(revision 36)
+++ /branches/petra/src/wetdepokernel.f90	(revision 37)
@@ -1,4 +1,4 @@
 !**********************************************************************
-! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
+! Copyright 1998-2015                                                 *
 ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
 ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
@@ -20,14 +20,18 @@
 !**********************************************************************
 
-subroutine wetdepokernel(nunc,deposit,x,y,nage,kp)
-  !                          i      i    i i  i
+subroutine wetdepokernel(nunc,deposit,x,y,itage,nage,kp)
+  !                          i      i    i  i    i   i
   !*****************************************************************************
   !                                                                            *
   !     Attribution of the deposition from an individual particle to the       *
-  !     deposition fields using a uniform kernel with bandwidths dxout and dyout.*
+  !     deposition fields using a uniform kernel with bandwidths dxout         *
+  !     and dyout.                                                             *
   !                                                                            *
   !     Author: A. Stohl                                                       *
   !                                                                            *
   !     26 December 1996                                                       *
+  !
+  !  PS, 2/2015: do not use kernel for itage < 3 h
+  !   same as for concentration in conccalc.f90
   !                                                                            *
   !*****************************************************************************
@@ -48,5 +52,5 @@
 
   real :: x,y,deposit(maxspec),ddx,ddy,xl,yl,wx,wy,w
-  integer :: ix,jy,ixp,jyp,nunc,nage,ks,kp
+  integer :: ix,jy,ixp,jyp,nunc,nage,ks,kp,itage
 
   xl=(x*dx+xoutshift)/dxout
@@ -74,37 +78,50 @@
 
 
-  ! Determine mass fractions for four grid points
-  !**********************************************
+  if (itage .lt. itagekernmin) then 
+    ! no kernel, direct attribution to grid cell
+    
+    do ks=1,nspec
+      if (ix.ge.0 .and. jy.ge.0 .and. &
+          ix.le.numxgrid-1 .and. jy.le.numygrid-1) &
+        wetgridunc(ix,jy,ks,kp,nunc,nage)= &
+          wetgridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)
+    enddo
+        
+  else
+    ! Determine mass fractions for four grid points & distribute
 
-  do ks=1,nspec
+    do ks=1,nspec
 
-  if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
-       (jy.le.numygrid-1)) then
-    w=wx*wy
-      wetgridunc(ix,jy,ks,kp,nunc,nage)= &
-           wetgridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w
+      if (ix.ge.0 .and. jy.ge.0 .and. &
+          ix.le.numxgrid-1 .and. jy.le.numygrid-1) then
+        w=wx*wy
+        wetgridunc(ix,jy,ks,kp,nunc,nage)= &
+          wetgridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w
+      endif
+
+      if ((ixp.ge.0).and.(jyp.ge.0).and.&
+         (ixp.le.numxgrid-1).and. (jyp.le.numygrid-1)) then
+        w=(1.-wx)*(1.-wy)
+        wetgridunc(ixp,jyp,ks,kp,nunc,nage)= &
+          wetgridunc(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w
+      endif
+
+      if ((ixp.ge.0).and.(jy.ge.0).and. &
+          (ixp.le.numxgrid-1).and.(jy.le.numygrid-1)) then
+        w=(1.-wx)*wy
+        wetgridunc(ixp,jy,ks,kp,nunc,nage)= &
+           wetgridunc(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w
+      endif
+
+      if ((ix.ge.0).and.(jyp.ge.0).and. &
+          (ix.le.numxgrid-1).and.(jyp.le.numygrid-1)) then
+        w=wx*(1.-wy)
+        wetgridunc(ix,jyp,ks,kp,nunc,nage)= &
+          wetgridunc(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w
+      endif
+
+    end do
+
   endif
 
-  if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgrid-1).and. &
-       (jyp.le.numygrid-1)) then
-    w=(1.-wx)*(1.-wy)
-      wetgridunc(ixp,jyp,ks,kp,nunc,nage)= &
-           wetgridunc(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w
-  endif
-
-  if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgrid-1).and. &
-       (jy.le.numygrid-1)) then
-    w=(1.-wx)*wy
-      wetgridunc(ixp,jy,ks,kp,nunc,nage)= &
-           wetgridunc(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w
-  endif
-
-  if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgrid-1).and. &
-       (jyp.le.numygrid-1)) then
-    w=wx*(1.-wy)
-      wetgridunc(ix,jyp,ks,kp,nunc,nage)= &
-           wetgridunc(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w
-  endif
-  end do
-
 end subroutine wetdepokernel
Index: /branches/petra/src/wetdepokernel_nest.f90
===================================================================
--- /branches/petra/src/wetdepokernel_nest.f90	(revision 36)
+++ /branches/petra/src/wetdepokernel_nest.f90	(revision 37)
@@ -1,4 +1,4 @@
 !**********************************************************************
-! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
+! Copyright 1998-2015                                                 *
 ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
 ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
@@ -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,itage,nage,kp)
+  !                             i      i    i i   i    i   i
   !*****************************************************************************
   !                                                                            *
@@ -35,4 +34,7 @@
   !      2 September 2004: Adaptation from wetdepokernel.                      *
   !                                                                            *
+  !
+  !  PS, 2/2015: do not use kernel for itage < 3 h
+  !   same as for concentration in conccalc.f90
   !                                                                            *
   !*****************************************************************************
@@ -53,5 +55,5 @@
 
   real :: x,y,deposit(maxspec),ddx,ddy,xl,yl,wx,wy,w
-  integer :: ix,jy,ixp,jyp,ks,kp,nunc,nage
+  integer :: ix,jy,ixp,jyp,ks,kp,nunc,nage,itage
 
 
@@ -80,38 +82,49 @@
   endif
 
+  if (itage .lt. itagekernmin) then  
+   ! no kernel, direct attribution to grid cell
+    
+    do ks=1,nspec
+      if (ix.ge.0 .and. jy.ge.0 .and. &
+          ix.le.numxgridn-1 .and. jy.le.numygridn-1) &
+        wetgriduncn(ix,jy,ks,kp,nunc,nage)= &
+          wetgriduncn(ix,jy,ks,kp,nunc,nage)+deposit(ks)
+    enddo
+        
+  else
+    ! Determine mass fractions for four grid points & distribute
 
-  ! Determine mass fractions for four grid points
-  !**********************************************
+    do ks=1,nspec
 
-  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
+        wetgriduncn(ix,jy,ks,kp,nunc,nage)= &
+          wetgriduncn(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w
+      endif
 
-  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
+      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
+
+      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
+
+      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
+
+    end do
+
   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)
-      wetgriduncn(ixp,jyp,ks,kp,nunc,nage)= &
-           wetgriduncn(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w
-  endif
-
-  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
-
-  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
-
-  end do
 end subroutine wetdepokernel_nest
