euler_rain was not called each time, firststep replaced by age

src/boundcond_domainfill.f90

@@ -309,7 +309,6 @@ subroutine boundcond_domainfill(itime,loutend)
                 itrasplit(ipart)=itra1(ipart)+ldirect*itsplit
                 xmass1(ipart,1)=xmassperparticle
                 if (WATERCYCLE) then
-                   firsttimestep(ipart)=1
                    call calculate_watercycle(ipart,itime)
                 endif
                 if (mdomainfill.eq.2) xmass1(ipart,1)= &

src/calculate_watercycle.f90

@@ -67,7 +67,7 @@ subroutine calculate_watercycle(partnumber,itime)
       forparticle=1
   endif
 
- ! Loop about all particles
+  ! Loop about all particles
   !*************************
   do j=1,forparticle
 
@@ -79,15 +79,11 @@ subroutine calculate_watercycle(partnumber,itime)
 
   ! Take only valid particles, itime invalid when particle released!
   !*****************************************************************
-
   if (itra1(i).eq.itime) then ! trajectory itra1 is not -999999
 
-  !*****************************************************************************
-  ! Interpolate several variables (PV, specific humidity, etc.) to particle position
-  !*****************************************************************************
 
-! Determine age class of the particle - nage is used for the kernel
-!******************************************************************
+  ! Determine age class of the particle - nage is used for the kernel
+  !******************************************************************
      itage=abs(itra1(i)-itramem(i))
      do nage=1,nageclass
        if (itage.lt.lage(nage)) goto 33
@@ -108,8 +104,8 @@ subroutine calculate_watercycle(partnumber,itime)
       p3=rddx*ddy
       p4=ddx*ddy
 
-  ! Potential vorticity, specific humidity, temperature, and density
-  !*****************************************************************
+  ! Interpolate specific humidity
+  !******************************
 
       do il=2,nz
         if (height(il).gt.ztra1(i)) then
@@ -129,7 +125,6 @@ subroutine calculate_watercycle(partnumber,itime)
         do m=1,2
           indexh=memind(m)
 
-  ! Specific humidity
           qv1(m)=p1*qv(ix ,jy ,ind,indexh) &
                +p2*qv(ixp,jy ,ind,indexh) &
                +p3*qv(ix ,jyp,ind,indexh) &
@@ -141,8 +136,10 @@ subroutine calculate_watercycle(partnumber,itime)
       end do
       qvi=(dz1*qvprof(2)+dz2*qvprof(1))*dz
 
-     
-    if (firsttimestep(i).gt.1) then ! not the initialisation stage, it is assumed a value is already saved
+!      write(*,*) 'all: ',i,diff,e_minus_p(ix,jy) &
+!              ,partnumber,val_q(i),qvi,itime,itra1(i),itage
+
+    if (itage.gt.0) then ! not the initialisation stage, it is assumed a value is already saved
         
   ! Calculate humidity differences 
   !*******************************
@@ -160,14 +157,16 @@ subroutine calculate_watercycle(partnumber,itime)
 !     Do this only for the 2nd timestep 
 !     there has to be a deltaq which is negative and the column has to have precipitation
 !     e_minus_p units is mm/day .. convert to mm/h
-      if ( ((diff.ge.0).or.(e_minus_p(ix,jy).gt.-2.0)) .and. (firsttimestep(i).eq.2)) then 
+
+! we cannot be sure 3600 is the second timestep, could also be 1800!
+
+      if ( ((diff.ge.0).or.((e_minus_p(ix,jy)/24).gt.-2.0)) .and. (itage.le.3600)) then 
+!            write(*,*) 'terminated: ',i,diff,e_minus_p(ix,jy) &
+!              ,partnumber,val_q(i),qvi,itime,itra1(i),itage
          itra1(i)=-999999999 
-         firsttimestep(i)=1 
-!         write(*,*) 'terminated: ',i,diff,e_minus_p(ix,jy),firsttimestep(i) &
-!           ,partnumber,val_q(i),qvi,itime,itra1(i),itage
       else
-!         write(*,*) 'accounted: ',i,diff,e_minus_p(ix,jy),firsttimestep(i) &
-!           ,partnumber,val_q(i),qvi,itime,itra1(i),itage
+!            write(*,*) 'accounted: ',i,diff,e_minus_p(ix,jy) &
+!              ,partnumber,val_q(i),qvi,itime,itra1(i),itage
          call centerofmass(xl,yl,2,centerposx_real,centerposy_real)
          centerposx=(centerposx_real-xlon0)/dx 
          centerposy=(centerposy_real-ylat0)/dy
@@ -177,7 +176,7 @@ subroutine calculate_watercycle(partnumber,itime)
             call wetdepokernel(1,abs(diff),centerposx,centerposy,nage,1)
          endif
       endif
-    endif ! firsttimestep gt 1
+    endif 
 
     ! save the old values 
     !********************
@@ -186,10 +185,6 @@ subroutine calculate_watercycle(partnumber,itime)
      ytra1_q(i)=ytra1(i)
      val_q(i)=qvi
     
-!     write(*,*) 'saved: ',i,xtra1_q(i),ytra1_q(i),val_q(i),firsttimestep(i)
-  
-     firsttimestep(i)=firsttimestep(i)+1 ! this is a valid particle keep in the simulation
-
   endif ! for valid trajectories
   end do
 

src/com_mod.f90

@@ -680,7 +680,6 @@ module com_mod
 
   real, allocatable, dimension(:,:) :: e_minus_p ! E-P field for the watercycle
   real, allocatable, dimension(:) :: xtra1_q, ytra1_q, val_q
-  integer, allocatable, dimension(:) :: firsttimestep ! to save the previous q value and position for the WATERCYCLE mode
 
   ! eso: Moved from timemanager
   real, allocatable, dimension(:) :: uap,ucp,uzp,us,vs,ws

src/getfields.f90

@@ -184,6 +184,8 @@ subroutine getfields(itime,nstop,metdata_format)
         call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn)
         memtime(1)=wftime(indj)
         memind(2)=2
+        if (WATERCYCLE) &
+           call euler_rain(itime) ! calculate the integrated E-P
         if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
           call readwind_ecmwf(indj+1,memind(2),uuh,vvh,wwh)
         else
@@ -199,6 +201,8 @@ subroutine getfields(itime,nstop,metdata_format)
         end if
         call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
         memtime(2)=wftime(indj+1)
+        if (WATERCYCLE) &
+           call euler_rain(itime) ! calculate the integrated E-P
         nstop = 1
         goto 60
       endif

src/init_domainfill.f90

@@ -278,7 +278,6 @@ subroutine init_domainfill
                      itsplit
                 xmass1(numpart+jj,1)=colmass(ix,jy)/real(ncolumn)
                 if (WATERCYCLE) then
-                   firsttimestep(numpart+jj)=1
                    call calculate_watercycle(numpart+jj,0)
                 endif
 

src/readcommand.f90

@@ -353,7 +353,6 @@ subroutine readcommand
          allocate(xtra1_q(maxpart))
          allocate(ytra1_q(maxpart))
          allocate(val_q(maxpart))
-         allocate(firsttimestep(maxpart))
          WATERCYCLE=.true.
          mdomainfill=1
      end select