added water threshold, write it in logfile

src/calculate_watercycle.f90

@@ -43,14 +43,16 @@ subroutine calculate_watercycle(partnumber,itime)
   real qv1(2),qvprof(2),qvi,dz,dz1,dz2
   real e_minus_p1(2), e_minus_pi
   
-  real  precip_reference(360,180),eminusp_saved
+  real  precip_reference(360,180),eminusp_saved, waterthreshold
   integer ixr,jyr,idumx,idumy,istep, watersynctime
 
 
-!    watersynctime=3600*3 ! was loutstep before, but for monthly run not possible
-   watersynctime=3600
+    watersynctime=3600*3 ! was loutstep before, but for monthly run not possible
+    waterthreshold=-.1
+!   watersynctime=3600
 !   watersynctime=loutstep
-  if ((partnumber.eq.1).and.(itime.eq.0)) write(*,*) 'Watersynctime ',watersynctime,lsynctime
+
+  if ((partnumber.eq.1).and.(itime.eq.0)) write(*,*) 'Watersynctime ',watersynctime,lsynctime,waterthreshold
 
   ! Some variables needed for temporal interpolation
   !*************************************************
@@ -229,7 +231,7 @@ subroutine calculate_watercycle(partnumber,itime)
 !     e_minus_pi=precip_reference(ix+1,jy-90+2)*(-24.)
 
 !      First timestep and criteria fullfilled, keep it or throw it?
-      if ( ((diff.ge.0).or.((e_minus_pi/24).gt.-.5)) .and. (status_q(i).eq.1) ) then 
+      if ( ((diff.ge.0).or.((e_minus_pi/24).gt.waterthreshold)) .and. (status_q(i).eq.1) ) then 
 !       if ( ((diff.ge.0).or.((e_minus_pi/24).gt.-2.0)) .and. (status_q(i).eq.1) ) then 
 !       if ( ((diff.ge.0)) .and. (status_q(i).eq.1) ) then 
 !       if  (status_q(i).eq.-100) then  ! never! all trajectories.

src/writeprecip.f90

@@ -41,7 +41,7 @@ subroutine writeprecip(itime,imem)
   character :: adate*8,atime*6
 
   integer :: ix,jy,imem
-  real :: xp1,yp1
+  real :: xp1,yp1,xp2,yp2
 
   
   if (itime.eq.0) then
@@ -60,10 +60,13 @@ subroutine writeprecip(itime,imem)
   do i=1,numpoint
     xp1=xpoint1(i)*dx+xlon0 !lat, long (real) coord
     yp1=ypoint1(i)*dy+ylat0 !lat, long (real) coord
+    xp1=xpoint2(i)*dx+xlon0 !lat, long (real) coord
+    yp1=ypoint2(i)*dy+ylat0 !lat, long (real) coord
     ix=int((xpoint1(i)+xpoint2(i))/2.)
     jy=int((ypoint1(i)+ypoint2(i))/2.)
     write(unitprecip,*)  jjjjmmdd, ihmmss, & 
-           xp1,yp1,lsprec(ix,jy,1,imem),convprec(ix,jy,1,imem) !time is the same as in the ECMWF windfield
+           (xp1+xp2)/2,(yp1+yp2)/2,lsprec(ix,jy,1,imem),convprec(ix,jy,1,imem) &
+           ,e_minus_p(ix,jy,imem) !time is the same as in the ECMWF windfield
 ! units mm/h, valid for the time given in the windfield
   end do