Added output of dilatation rate

src/do_leaf_measurements.f90

@@ -60,7 +60,8 @@ double precision,dimension(:),allocatable :: dhdx,dhdy,dhdz
 double precision :: r,s,t,volume,w
 double precision :: exx,eyy,ezz,exy,eyz,ezx,J2d,J3d
 integer err,ierr,nproc,iproc,k,i,iint,nleaves,mpe
-double precision, dimension(:), allocatable :: e2dtemp,e3dtemp,crittemp,lodetemp,qtemp
+double precision, dimension(:), allocatable :: e2dtemp,e3dtemp,crittemp,lodetemp
+double precision, dimension(:), allocatable :: qtemp,dilatrtemp
 character*72  :: shift
 logical cutcell
 double precision :: prod
@@ -103,18 +104,21 @@ allocate (e3dtemp(nleaves),stat=err)              ; if (err.ne.0) call stop_run
 allocate (crittemp(nleaves),stat=err)             ; if (err.ne.0) call stop_run ('Error alloc crittemp in compute_e2d_crit$')
 allocate (lodetemp(nleaves),stat=err)             ; if (err.ne.0) call stop_run ('Error alloc lodetemp in compute_e2d_crit$')
 allocate (qtemp(nleaves),stat=err)                ; if (err.ne.0) call stop_run ('Error alloc qtemp in compute_e2d_crit$')
+allocate (dilatrtemp(nleaves),stat=err)             ; if (err.ne.0) call stop_run ('Error alloc dilatrtemp in compute_e2d_crit$')
 
 e2dtemp=0.d0
 e3dtemp=0.d0
 crittemp=0.d0
 lodetemp=0.d0
 qtemp=0.d0
+dilatrtemp=0.d0
 
 osolve%e2d=0.d0
 osolve%e3d=0.d0
 osolve%crit=0.d0
 osolve%lode=0.d0
 osolve%q=0.d0
+osolve%dilatr=0.d0
 
 do i=1+iproc,nleaves,nproc
 
@@ -159,6 +163,9 @@ do i=1+iproc,nleaves,nproc
       qtemp(i)=2.d0*osolve%eviscosity(i)*sqrt(second_invariant(exx,eyy,ezz,exy,eyz,ezx))
 !   end if
 
+   ! Dilatation rate
+   dilatrtemp(i)=(exx+eyy+ezz)
+
    J2d=second_invariant(exx,eyy,ezz,exy,eyz,ezx)
    J3d=third_invariant (exx,eyy,ezz,exy,eyz,ezx)
 
@@ -192,6 +199,7 @@ call mpi_reduce (e3dtemp, osolve%e3d, nleaves,mpi_double_precision,mpi_sum,0,mpi
 call mpi_reduce (crittemp,osolve%crit,nleaves,mpi_double_precision,mpi_sum,0,mpi_comm_world,ierr)
 call mpi_reduce (lodetemp,osolve%lode,nleaves,mpi_double_precision,mpi_sum,0,mpi_comm_world,ierr)
 call mpi_reduce (qtemp, osolve%q, nleaves,mpi_double_precision,mpi_sum,0,mpi_comm_world,ierr)
+call mpi_reduce (dilatrtemp, osolve%dilatr, nleaves,mpi_double_precision,mpi_sum,0,mpi_comm_world,ierr)
 
 deallocate (x,y,z)
 deallocate (vx,vy,vz)
@@ -201,6 +209,7 @@ deallocate (e2dtemp)
 deallocate (e3dtemp)
 deallocate (lodetemp) 
 deallocate (qtemp) 
+deallocate (dilatrtemp)
 
 if (iproc.eq.0 .and. params%debug>=1 .and. params%compute_qpgram) then
    write(*,'(a,f11.5)')  shift//'min(q) =', minval(osolve%q)
@@ -208,7 +217,7 @@ if (iproc.eq.0 .and. params%debug>=1 .and. params%compute_qpgram) then
 end if     
 
 call DoRuRe_vel_stats  (params%doDoRuRe,istep,iter,iter_nl,ov%nnode,ov%unode,ov%vnode,ov%wnode)
-call DoRuRe_leaf_stats (params%doDoRuRe,osolve%nleaves,osolve%e2d,osolve%e3d,osolve%pressure,osolve%q,istep,iter,iter_nl)
+call DoRuRe_leaf_stats (params%doDoRuRe,osolve%nleaves,osolve%e2d,osolve%e3d,osolve%pressure,osolve%q,osolve%dilatr,istep,iter,iter_nl)
 
 return
 end