interpolation of E-P field

src/calculate_watercycle.f90

@@ -41,6 +41,7 @@ subroutine calculate_watercycle(partnumber,itime)
   real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4
   real centerposx_real, centerposy_real, centerposx, centerposy, diff, xl(2), yl(2)
   real qv1(2),qvprof(2),qvi,dz,dz1,dz2
+  real e_minus_p1(2), e_minus_pi
 
 
   ! Some variables needed for temporal interpolation
@@ -136,6 +137,18 @@ subroutine calculate_watercycle(partnumber,itime)
       end do
       qvi=(dz1*qvprof(2)+dz2*qvprof(1))*dz
 
+      do m=1,2
+         indexh=memind(m)
+
+         e_minus_p1(m)=p1*e_minus_p(ix ,jy ,indexh) &
+                +p2*e_minus_p(ixp,jy ,indexh) &
+                +p3*e_minus_p(ix ,jyp,indexh) &
+                +p4*e_minus_p(ixp,jyp,indexh)
+
+      end do
+      e_minus_pi=(e_minus_p1(1)*dt2+dt1*e_minus_p1(2))*dtt
+
+
 !      write(*,*) 'all: ',i,diff,e_minus_p(ix,jy) &
 !              ,partnumber,val_q(i),qvi,itime,itra1(i),itage
 
@@ -160,7 +173,7 @@ subroutine calculate_watercycle(partnumber,itime)
 
 ! 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 
+      if ( ((diff.ge.0).or.((e_minus_pi/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 

src/com_mod.f90

@@ -678,7 +678,7 @@ module com_mod
   real, allocatable, dimension(:,:) :: xmass1
   real, allocatable, dimension(:,:) :: xscav_frac1
 
-  real, allocatable, dimension(:,:) :: e_minus_p ! E-P field for the watercycle
+  real, allocatable, dimension(:,:,:) :: e_minus_p ! E-P field for the watercycle
   real, allocatable, dimension(:) :: xtra1_q, ytra1_q, val_q
 
   ! eso: Moved from timemanager

src/euler_rain.f90

@@ -143,11 +143,13 @@
 
 56          continue
 
-          e_minus_p(ix,jy)=dwatercolumndt(ix,jy)+divwatercolumn(ix,jy,3)
+          n=memind(1)
+          e_minus_p(ix,jy,n)=   &
+               dwatercolumndt(ix,jy)+divwatercolumn(ix,jy,3)
 
 ! Convert to mm/day
 
-          e_minus_p(ix,jy)=e_minus_p(ix,jy)*86400.
+          e_minus_p(ix,jy,n)=e_minus_p(ix,jy,n)*86400.
 
 50        continue
 
@@ -155,23 +157,17 @@
 
 ! Set values at the poles equal to values 1 deg equatorward
 !**********************************************************
+      do n=1,2
 
       do 80 ix=1,nx-1
-        e_minus_p(ix,ny-1)=e_minus_p(ix,ny-2)
-        e_minus_p(ix,0)=e_minus_p(ix,1)
+        e_minus_p(ix,ny-1,n)=e_minus_p(ix,ny-2,n)
+        e_minus_p(ix,0,n)=e_minus_p(ix,1,n)
 80      continue
 
-      write(89,'(i8,1x,a11)') itime,' EminusP    '
+      write(89,'(i2,i8,1x,a11)') n,itime,' EminusP    '
       do 85 jy=1,ny-2
-          write(89,*) (e_minus_p(ix,jy),ix=1,nx-2)
+          write(89,*) (e_minus_p(ix,jy,n),ix=1,nx-2)
 85        continue
 
-! Use results only if we are inside the simulation period (don't use the first fields,
-! which are probably only needed for interpolation to time 0
-!*************************************************************************************
-
-      if (memtime(2).gt.0) then
-      endif
-
-
+      end do
       end

src/outgrid_init.f90

@@ -299,7 +299,7 @@ subroutine outgrid_init
           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))
+       allocate(e_minus_p(0:nxmax-1,0:nymax-1,2))
      endif
   endif