writing out budget introducing variable for checking of 2nd timestep

src/calculate_watercycle.f90

@@ -75,6 +75,7 @@ subroutine calculate_watercycle(partnumber,itime)
        i=j ! loop over all particle
    else
        i=partnumber ! just for one particle
+       firsttimestep(i)=.true.     
    endif
 
   ! Take only valid particles, itime invalid when particle released!
@@ -183,18 +184,25 @@ subroutine calculate_watercycle(partnumber,itime)
         
 !     Do this only for the 2nd timestep i
 !     there has to be a deltaq which is negative and the column has to have precipitation
-      if (((prec_q(i)+prec).gt.0).and.(e_minus_p(ix,jy).gt..2)) then 
+!     e_minus_p units is mm/day .. convert to mm/h
+      if ( ((diff.ge.0).or.(e_minus_p(ix,jy)/24.gt.-2.0)) .and. firsttimestep(i)) then 
+!      if ( ((diff.ge.0)) .and. firsttimestep(i)) then 
          itra1(i)=-999999999 
+!        write(*,*) 'Trajectory terminated: ',i,diff,e_minus_p(ix,jy),firsttimestep(i)
+      else
+!        write(*,*) 'Trajectory accounted: ',i,diff,e_minus_p(ix,jy),firsttimestep(i),&
+!                          val_q(i),qvi,xmass1(i,1)
+         prec_q(i)=prec
+         call centerofmass(xl,yl,2,centerposx_real,centerposy_real)
+         centerposx=(centerposx_real-xlon0)/dx 
+         centerposy=(centerposy_real-ylat0)/dy
+         if (diff.gt.0) then
+            call drydepokernel(1,diff,centerposx,centerposy,nage,1)
+         else
+            call wetdepokernel(1,abs(diff),centerposx,centerposy,nage,1)
+         endif
       endif
-          prec_q(i)=prec
-          call centerofmass(xl,yl,2,centerposx_real,centerposy_real)
-          centerposx=(centerposx_real-xlon0)/dx 
-          centerposy=(centerposy_real-ylat0)/dy
-          if (diff.gt.0) then
-             call drydepokernel(1,diff,centerposx,centerposy,nage,1)
-          else
-             call wetdepokernel(1,abs(diff),centerposx,centerposy,nage,1)
-          endif
+      firsttimestep(i)=.false.
     endif
 
   endif

src/com_mod.f90

@@ -679,7 +679,8 @@ module com_mod
   real, allocatable, dimension(:,:) :: xscav_frac1
 
   real, allocatable, dimension(:,:) :: e_minus_p ! E-P field for the watercycle
-  real, allocatable, dimension(:) :: xtra1_q, ytra1_q, val_q, prec_q ! to save the previous q value and position for the WATERCYCLE mode
+  real, allocatable, dimension(:) :: xtra1_q, ytra1_q, val_q, prec_q
+  logical, 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/euler_rain.f90

-      subroutine euler_rain()
+      subroutine euler_rain(itime)
 !*******************************************************************************
 ! Compute E-P between 2 time steps from moisture budget in an atmospheric column
 ! 2 Terms: 1. Moisture column tendency between the 2 time steps,               *
@@ -14,6 +14,8 @@
       use par_mod
       use com_mod
 
+      integer itime
+
       integer ix,jy,kz,k,n,ixp,jyp,ixm,jym
       real pih,ylat,ylatp,ylatm
       parameter(pih=pi/180.)
@@ -24,21 +26,17 @@
       real watercolumn(2),watermass(nzmax,2)
       real dwatercolumndt(0:nxmax-1,0:nymax-1)
       real divwatercolumn(0:nxmax-1,0:nymax-1,3)
-      integer nmax,nprocess
+      integer nmax
       parameter(nmax=nxmax*nymax)
 
-      open(89,file='budget.dat')
       write(89,*) xlon0+dx,ylat0+dy,nx-2,ny-2,dx,dy,1,10000.,3
-      write(89,*) 1
 
-             write(*,*) nprocess
 ! Loop over the whole grid
 !*************************
 
       do 10 jy=0,ny-1
         do 10 ix=0,nx-1
 
-
 !*****************************************************************
 ! 1. Compute the precipitable water tendency in atmospheric column
 !    and the fluxes needed for the divergence term (step 2)
@@ -163,10 +161,9 @@
         e_minus_p(ix,0)=e_minus_p(ix,1)
 80      continue
 
-      write(89,'(2i8,1x,a11)') ibdate,iedate,'EminusP    '
+      write(89,'(i8,1x,a11)') itime,' EminusP    '
       do 85 jy=1,ny-2
-        do 85 ix=1,nx-2
-          write(89,*) ix,jy,e_minus_p(ix,jy)
+          write(89,*) (e_minus_p(ix,jy),ix=1,nx-2)
 85        continue
 
 ! Use results only if we are inside the simulation period (don't use the first fields,
@@ -176,6 +173,5 @@
       if (memtime(2).gt.0) then
       endif
 
-      close(89)
 
       end

src/getfields.f90

@@ -147,7 +147,7 @@ subroutine getfields(itime,nstop,metdata_format)
         call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
         memtime(2)=wftime(indj+1)
         if (WATERCYCLE) &
-           call euler_rain() ! calculate the integrated E-P
+           call euler_rain(itime) ! calculate the integrated E-P
         nstop = 1
         goto 40
       endif

src/makefile

@@ -90,11 +90,11 @@ O_LEV_DBG = g # [0,g]
 ## LIBRARIES
 LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper -lnetcdff # -fopenmp
 
-#FFLAGS   = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(FUSER) # -Warray-bounds -fcheck=all # -march=native
-FFLAGS   = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(FUSER)  -Warray-bounds -fcheck=all # -march=native
+FFLAGS   = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(FUSER) # -Warray-bounds -fcheck=all # -march=native
+#FFLAGS   = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(FUSER)  -Warray-bounds -fcheck=all # -march=native
 
-#DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) -fbacktrace   -Wall  -fdump-core $(FUSER)   # -ffpe-trap=invalid,overflow,denormal,underflow,zero  -Warray-bounds -fcheck=all
-DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) -fbacktrace   -Wall  -fdump-core $(FUSER)  -Warray-bounds -fcheck=all
+DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) -fbacktrace   -Wall  -fdump-core $(FUSER)   # -ffpe-trap=invalid,overflow,denormal,underflow,zero  -Warray-bounds -fcheck=all
+#DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) -fbacktrace   -Wall  -fdump-core $(FUSER)  -Warray-bounds -fcheck=all
 
 LDFLAGS  = $(FFLAGS) -L$(LIBPATH1) -Wl,-rpath,$(LIBPATH1) $(LIBS) #-L$(LIBPATH2)
 LDDEBUG  = $(DBGFLAGS) -L$(LIBPATH1) $(LIBS) #-L$(LIBPATH2)

src/outgrid_init.f90

@@ -298,6 +298,9 @@ subroutine outgrid_init
      allocate(drygrid(0:max(numxgrid,numxgridn)-1, &
           0:max(numygrid,numygridn)-1),stat=stat)
      if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
+     if (WATERCYCLE) then
+       allocate(e_minus_p(0:nxmax-1,0:nymax-1))
+     endif
   endif
 
   ! Initial condition field

src/readcommand.f90

@@ -354,7 +354,7 @@ subroutine readcommand
          allocate(ytra1_q(maxpart))
          allocate(val_q(maxpart))
          allocate(prec_q(maxpart))
-         allocate(e_minus_p(0:numxgrid-1,0:numygrid-1))
+         allocate(firsttimestep(maxpart))
          WATERCYCLE=.true.
          mdomainfill=1
      end select

src/readwind_ecmwf.f90

@@ -488,7 +488,7 @@ subroutine readwind_ecmwf(indj,n,uuh,vvh,wwh)
   do i=0,nxmin1
     do j=0,nymin1
       if (WATERCYCLE) then
-        do ki=0,nuvz
+        do ki=1,nuvz
          uul(i,j,ki,n)=uuh(i,j,ki)
          vvl(i,j,ki,n)=vvh(i,j,ki)
         end do

src/timemanager.f90

@@ -126,7 +126,7 @@ subroutine timemanager(metdata_format)
   !    +                 xm_depd(maxspec,maxpointspec_act)
 
 
-  !open(88,file='TEST.dat')
+  open(89,file=path(2)(1:length(2))//'budget.dat')
 
   ! First output for time 0
   !************************
@@ -737,13 +737,15 @@ subroutine timemanager(metdata_format)
   if (ipout.eq.2) call partoutput_short(itime)     ! dump particle positions
 
       if (WATERCYCLE) then
-           call calculate_watercycle(-1,itime)   ! initilize4
+           call calculate_watercycle(-1,itime)   
           write (*,*) 'last watercycle'
       endif
      
   if (linit_cond.ge.1) call initial_cond_output(itime)   ! dump initial cond. field
 
-  !close(104)
+  if (WATERCYCLE) then
+      close(89) ! Euler rain, water budget file
+  endif
 
   ! De-allocate memory and end
   !***************************