Merge pull request #85 from pbartholomew08/fix/force-in-forces-module

Move force subroutine into forces module

Author: slaizet

src/forces.f90

@@ -203,146 +203,141 @@ subroutine restart_forces(itest1)
     endif
 
   end subroutine restart_forces
-end module forces
 
-!***********************************************************************
-subroutine force(ux1,uy1,ep1)
+  subroutine force(ux1,uy1,ep1)
 
-  !***********************************************************************
+    USE param
+    USE variables
+    USE decomp_2d
+    USE MPI
+    USE ibm_param
 
-  USE forces
-  USE param
-  USE variables
-  USE decomp_2d
-  USE MPI
-  USE ibm_param
+    use var, only : ta1, tb1, tc1, td1, di1
+    use var, only : ux2, uy2, ta2, tb2, tc2, td2, di2
 
-  use var, only : ta1, tb1, tc1, td1, di1
-  use var, only : ux2, uy2, ta2, tb2, tc2, td2, di2
+    implicit none
+    character(len=30) :: filename, filename2
+    integer :: nzmsize
+    integer                                             :: i, iv, j, k, kk, code, jj
+    integer                                             :: nvect1,nvect2,nvect3
 
-  implicit none
-  character(len=30) :: filename, filename2
-  integer :: nzmsize
-  integer                                             :: i, iv, j, k, kk, code, jj
-  integer                                             :: nvect1,nvect2,nvect3
+    real(mytype), dimension(xsize(1),xsize(2),xsize(3)),intent(in) :: ux1, uy1
+    real(mytype), dimension(xsize(1),xsize(2),xsize(3)),intent(in) :: ep1
 
-  real(mytype), dimension(xsize(1),xsize(2),xsize(3)),intent(in) :: ux1, uy1
-  real(mytype), dimension(xsize(1),xsize(2),xsize(3)),intent(in) :: ep1
+    real(mytype), dimension(ysize(1),ysize(2),ysize(3)) :: ppi2
 
-  real(mytype), dimension(ysize(1),ysize(2),ysize(3)) :: ppi2
+    real(mytype), dimension(nz) :: yLift,xDrag
+    real(mytype) :: yLift_mean,xDrag_mean
 
-  real(mytype), dimension(nz) :: yLift,xDrag
-  real(mytype) :: yLift_mean,xDrag_mean
+    real(mytype), dimension(ny-1) :: del_y
 
-  real(mytype), dimension(ny-1) :: del_y
+    real(mytype), dimension(nz) :: tunstxl, tunstyl
+    real(mytype), dimension(nz) :: tconvxl,tconvyl
+    real(mytype), dimension(nz) :: tpresxl,tpresyl
+    real(mytype), dimension(nz) :: tdiffxl,tdiffyl
 
-  real(mytype), dimension(nz) :: tunstxl, tunstyl
-  real(mytype), dimension(nz) :: tconvxl,tconvyl
-  real(mytype), dimension(nz) :: tpresxl,tpresyl
-  real(mytype), dimension(nz) :: tdiffxl,tdiffyl
+    real(mytype), dimension(nz) :: tunstx, tunsty
+    real(mytype), dimension(nz) :: tconvx,tconvy
+    real(mytype), dimension(nz) :: tpresx,tpresy
+    real(mytype), dimension(nz) :: tdiffx,tdiffy
 
-  real(mytype), dimension(nz) :: tunstx, tunsty
-  real(mytype), dimension(nz) :: tconvx,tconvy
-  real(mytype), dimension(nz) :: tpresx,tpresy
-  real(mytype), dimension(nz) :: tdiffx,tdiffy
+    real(mytype) :: uxmid,uymid,prmid
+    real(mytype) :: dudxmid,dudymid,dvdxmid,dvdymid
+    real(mytype) :: fac,tsumx,tsumy
+    real(mytype) :: fcvx,fcvy,fprx,fpry,fdix,fdiy
+    real(mytype) :: xmom,ymom
 
-  real(mytype) :: uxmid,uymid,prmid
-  real(mytype) :: dudxmid,dudymid,dvdxmid,dvdymid
-  real(mytype) :: fac,tsumx,tsumy
-  real(mytype) :: fcvx,fcvy,fprx,fpry,fdix,fdiy
-  real(mytype) :: xmom,ymom
+    nvect1=xsize(1)*xsize(2)*xsize(3)
+    nvect2=ysize(1)*ysize(2)*ysize(3)
+    nvect3=zsize(1)*zsize(2)*zsize(3)
 
-  nvect1=xsize(1)*xsize(2)*xsize(3)
-  nvect2=ysize(1)*ysize(2)*ysize(3)
-  nvect3=zsize(1)*zsize(2)*zsize(3)
+    do jj = 1, ny-1
+       if (istret.eq.0) then
+          del_y(jj)=dy
+       else
+          del_y(jj)=yp(jj+1)-yp(jj) 
+       endif
+    enddo
 
-  do jj = 1, ny-1
-    if (istret.eq.0) then
-      del_y(jj)=dy
-    else
-      del_y(jj)=yp(jj+1)-yp(jj) 
+    if (itime.eq.1) then
+       do k = 1, xsize(3)
+          do j = 1, xsize(2)
+             do i = 1, xsize(1)
+                ux11(i,j,k)=ux1(i,j,k)
+                uy11(i,j,k)=uy1(i,j,k)
+             enddo
+          enddo
+       enddo
+       return
+    elseif (itime.eq.2) then
+       do k = 1, xsize(3)
+          do j = 1, xsize(2)
+             do i = 1, xsize(1)
+                ux01(i,j,k)=ux1(i,j,k)
+                uy01(i,j,k)=uy1(i,j,k)
+             enddo
+          enddo
+       enddo
+       return
     endif
-  enddo
-
-  if (itime.eq.1) then
-     do k = 1, xsize(3)
-        do j = 1, xsize(2)
-           do i = 1, xsize(1)
-              ux11(i,j,k)=ux1(i,j,k)
-              uy11(i,j,k)=uy1(i,j,k)
-           enddo
-        enddo
-     enddo
-     return
-  elseif (itime.eq.2) then
-     do k = 1, xsize(3)
-        do j = 1, xsize(2)
-           do i = 1, xsize(1)
-              ux01(i,j,k)=ux1(i,j,k)
-              uy01(i,j,k)=uy1(i,j,k)
-           enddo
-        enddo
-     enddo
-     return
-  endif
-
-  call derx (ta1,ux1,di1,sx,ffx,fsx,fwx,xsize(1),xsize(2),xsize(3),0,ubcx)    ! dudx
-  call derx (tb1,uy1,di1,sx,ffxp,fsxp,fwxp,xsize(1),xsize(2),xsize(3),1,ubcy) ! dvdx
-  call transpose_x_to_y(ta1,ta2) ! dudx
-  call transpose_x_to_y(tb1,tb2) ! dvdx
-
-  call transpose_x_to_y(ux1,ux2)
-  call transpose_x_to_y(uy1,uy2)
-  call transpose_x_to_y(ppi1,ppi2)
-
-  call dery (tc2,ux2,di2,sy,ffyp,fsyp,fwyp,ppy,ysize(1),ysize(2),ysize(3),1,ubcx) ! dudy
-  call dery (td2,uy2,di2,sy,ffy,fsy,fwy,ppy,ysize(1),ysize(2),ysize(3),0,ubcy)    ! dvdy
-  call transpose_y_to_x(tc2,tc1) ! dudy
-  call transpose_y_to_x(td2,td1) ! dvdy
-
-  !*****************************************************************
-  !      Drag and Lift coefficients
-  !*****************************************************************
-  do iv=1,nvol
-
-     !*****************************************************************
-     !        Calculation of the momentum terms
-     !*****************************************************************
-     !
-     !     Calculation of the momentum terms. First we integrate the 
-     !     time rate of momentum along the CV.
-     !
-     !         Excluding the body internal cells. If the centroid 
-     !         of the cell falls inside the body the cell is 
-     !         excluded.
-
-     tunstxl=zero
-     tunstyl=zero
-     do k=1,xsize(3)
-        tsumx=zero
-        tsumy=zero
-        do j=jcvlw_lx(iv),jcvup_lx(iv)
-           do i=icvlf_lx(iv),icvrt_lx(iv)
-              !     The velocity time rate has to be relative to the cell center, 
-              !     and not to the nodes, because, here, we have an integral 
-              !     relative to the volume, and, therefore, this has a sense 
-              !     of a "source".
-              !         fac   = (1.5*ux1(i,j,k)-2.0*ux01(i,j,k)+0.5*ux11(i,j,k))*epcv1(i,j,k)
-              fac   = (onepfive*ux1(i,j,k)-two*ux01(i,j,k)+half*ux11(i,j,k))*(one-ep1(i,j,k))
-              tsumx = tsumx+fac*dx*del_y(j+(xstart(2)-1))/dt    !tsumx+fac*dx*dy/dt
-              !sumx(k) = sumx(k)+dudt1*dx*dy
-
-              !         fac   = (1.5*uy1(i,j,k)-2.0*uy01(i,j,k)+0.5*uy11(i,j,k))*epcv1(i,j,k)
-              fac   = (onepfive*uy1(i,j,k)-two*uy01(i,j,k)+half*uy11(i,j,k))*(one-ep1(i,j,k))
-              tsumy = tsumy+fac*dx*del_y(j+(xstart(2)-1))/dt !tsumy+fac*dx*dy/dt
-              !sumy(k) = sumy(k)+dudt1*dx*dy
-           enddo
-        enddo
-        tunstxl(xstart(3)-1+k)=tsumx
-        tunstyl(xstart(3)-1+k)=tsumy
-     enddo
-     call MPI_ALLREDUCE(tunstxl,tunstx,nz,real_type,MPI_SUM,MPI_COMM_WORLD,code)
-     call MPI_ALLREDUCE(tunstyl,tunsty,nz,real_type,MPI_SUM,MPI_COMM_WORLD,code)
+
+    call derx (ta1,ux1,di1,sx,ffx,fsx,fwx,xsize(1),xsize(2),xsize(3),0,ubcx)    ! dudx
+    call derx (tb1,uy1,di1,sx,ffxp,fsxp,fwxp,xsize(1),xsize(2),xsize(3),1,ubcy) ! dvdx
+    call transpose_x_to_y(ta1,ta2) ! dudx
+    call transpose_x_to_y(tb1,tb2) ! dvdx
+
+    call transpose_x_to_y(ux1,ux2)
+    call transpose_x_to_y(uy1,uy2)
+    call transpose_x_to_y(ppi1,ppi2)
+
+    call dery (tc2,ux2,di2,sy,ffyp,fsyp,fwyp,ppy,ysize(1),ysize(2),ysize(3),1,ubcx) ! dudy
+    call dery (td2,uy2,di2,sy,ffy,fsy,fwy,ppy,ysize(1),ysize(2),ysize(3),0,ubcy)    ! dvdy
+    call transpose_y_to_x(tc2,tc1) ! dudy
+    call transpose_y_to_x(td2,td1) ! dvdy
+
+    !*****************************************************************
+    !      Drag and Lift coefficients
+    !*****************************************************************
+    do iv=1,nvol
+
+       !*****************************************************************
+       !        Calculation of the momentum terms
+       !*****************************************************************
+       !
+       !     Calculation of the momentum terms. First we integrate the 
+       !     time rate of momentum along the CV.
+       !
+       !         Excluding the body internal cells. If the centroid 
+       !         of the cell falls inside the body the cell is 
+       !         excluded.
+
+       tunstxl=zero
+       tunstyl=zero
+       do k=1,xsize(3)
+          tsumx=zero
+          tsumy=zero
+          do j=jcvlw_lx(iv),jcvup_lx(iv)
+             do i=icvlf_lx(iv),icvrt_lx(iv)
+                !     The velocity time rate has to be relative to the cell center, 
+                !     and not to the nodes, because, here, we have an integral 
+                !     relative to the volume, and, therefore, this has a sense 
+                !     of a "source".
+                !         fac   = (1.5*ux1(i,j,k)-2.0*ux01(i,j,k)+0.5*ux11(i,j,k))*epcv1(i,j,k)
+                fac   = (onepfive*ux1(i,j,k)-two*ux01(i,j,k)+half*ux11(i,j,k))*(one-ep1(i,j,k))
+                tsumx = tsumx+fac*dx*del_y(j+(xstart(2)-1))/dt    !tsumx+fac*dx*dy/dt
+                !sumx(k) = sumx(k)+dudt1*dx*dy
+
+                !         fac   = (1.5*uy1(i,j,k)-2.0*uy01(i,j,k)+0.5*uy11(i,j,k))*epcv1(i,j,k)
+                fac   = (onepfive*uy1(i,j,k)-two*uy01(i,j,k)+half*uy11(i,j,k))*(one-ep1(i,j,k))
+                tsumy = tsumy+fac*dx*del_y(j+(xstart(2)-1))/dt !tsumy+fac*dx*dy/dt
+                !sumy(k) = sumy(k)+dudt1*dx*dy
+             enddo
+          enddo
+          tunstxl(xstart(3)-1+k)=tsumx
+          tunstyl(xstart(3)-1+k)=tsumy
+       enddo
+       call MPI_ALLREDUCE(tunstxl,tunstx,nz,real_type,MPI_SUM,MPI_COMM_WORLD,code)
+       call MPI_ALLREDUCE(tunstyl,tunsty,nz,real_type,MPI_SUM,MPI_COMM_WORLD,code)
 
 !!$!*********************************************************************************
 !!$!     Secondly, the surface momentum fluxes
@@ -357,209 +352,210 @@ subroutine force(ux1,uy1,ep1)
 !!$!        \       CV   \
 !!$!(jcvlw) A____________D      
 
-     tconvxl=zero
-     tconvyl=zero
-     tdiffxl=zero
-     tdiffyl=zero
-     tpresxl=zero
-     tpresyl=zero
-     !BC and AD : x-pencils
-     !AD
-     if ((jcvlw(iv).ge.xstart(2)).and.(jcvlw(iv).le.xend(2))) then
-        j=jcvlw(iv)-xstart(2)+1
-        do k=1,xsize(3)
-           kk=xstart(3)-1+k
-           fcvx=zero
-           fcvy=zero
-           fpry=zero
-           fdix=zero
-           fdiy=zero
-           do i=icvlf_lx(iv),icvrt_lx(iv)-1
-              !momentum flux
-              uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k))
-              uymid = half*(uy1(i,j,k)+uy1(i+1,j,k))
-              fcvx  = fcvx -uxmid*uymid*dx
-              fcvy  = fcvy -uymid*uymid*dx
-
-              !pressure
-              prmid = half*(ppi1(i,j,k)+ppi1(i+1,j,k))
-              fpry  = fpry +prmid*dx
-
-              !viscous term
-              dudymid = half*(tc1(i,j,k)+tc1(i+1,j,k))
-              dvdxmid = half*(tb1(i,j,k)+tb1(i+1,j,k))
-              dvdymid = half*(td1(i,j,k)+td1(i+1,j,k))
-              fdix    = fdix -(xnu*(dudymid+dvdxmid)*dx)
-              fdiy    = fdiy -two*xnu*dvdymid*dx
-
-           enddo
-
-           tconvxl(kk)=tconvxl(kk)+fcvx
-           tconvyl(kk)=tconvyl(kk)+fcvy
-           tpresyl(kk)=tpresyl(kk)+fpry
-           tdiffxl(kk)=tdiffxl(kk)+fdix
-           tdiffyl(kk)=tdiffyl(kk)+fdiy
-        enddo
-     endif
-     !BC
-     if ((jcvup(iv).ge.xstart(2)).and.(jcvup(iv).le.xend(2))) then
-        j=jcvup(iv)-xstart(2)+1
-        do k=1,xsize(3)
-           kk=xstart(3)-1+k
-           fcvx=zero
-           fcvy=zero
-           fpry=zero
-           fdix=zero
-           fdiy=zero
-           do i=icvlf_lx(iv),icvrt_lx(iv)-1
-              !momentum flux
-              uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k))
-              uymid = half*(uy1(i,j,k)+uy1(i+1,j,k))
-              fcvx= fcvx +uxmid*uymid*dx
-              fcvy= fcvy +uymid*uymid*dx
-
-              !pressure
-              prmid = half*(ppi1(i,j,k)+ppi1(i+1,j,k))
-              fpry = fpry -prmid*dx
-
-              !viscous term
-              dudymid = half*(tc1(i,j,k)+tc1(i+1,j,k))
-              dvdxmid = half*(tb1(i,j,k)+tb1(i+1,j,k))
-              dvdymid = half*(td1(i,j,k)+td1(i+1,j,k))
-              fdix = fdix +(xnu*(dudymid+dvdxmid)*dx)
-              fdiy = fdiy +two*xnu*dvdymid*dx
-
-           enddo
-           tconvxl(kk)=tconvxl(kk)+fcvx
-           tconvyl(kk)=tconvyl(kk)+fcvy
-           tpresyl(kk)=tpresyl(kk)+fpry
-           tdiffxl(kk)=tdiffxl(kk)+fdix
-           tdiffyl(kk)=tdiffyl(kk)+fdiy
-        enddo
-     endif
-     !AB and DC : y-pencils
-     !AB
-     if ((icvlf(iv).ge.ystart(1)).and.(icvlf(iv).le.yend(1))) then
-        i=icvlf(iv)-ystart(1)+1
-        do k=1,ysize(3)
-           kk=ystart(3)-1+k
-           fcvx=zero
-           fcvy=zero
-           fprx=zero
-           fdix=zero
-           fdiy=zero
-           do j=jcvlw_ly(iv),jcvup_ly(iv)-1
-              !momentum flux
-              uxmid = half*(ux2(i,j,k)+ux2(i,j+1,k))
-              uymid = half*(uy2(i,j,k)+uy2(i,j+1,k))
-              fcvx= fcvx -uxmid*uxmid*del_y(j)
-              fcvy= fcvy -uxmid*uymid*del_y(j)
-
-              !pressure
-              prmid=half*(ppi2(i,j,k)+ppi2(i,j+1,k))
-              fprx = fprx +prmid*del_y(j)
-
-              !viscous term
-              dudxmid = half*(ta2(i,j,k)+ta2(i,j+1,k))
-              dudymid = half*(tc2(i,j,k)+tc2(i,j+1,k))
-              dvdxmid = half*(tb2(i,j,k)+tb2(i,j+1,k))
-              fdix = fdix -two*xnu*dudxmid*del_y(j)
-              fdiy = fdiy -xnu*(dvdxmid+dudymid)*del_y(j)
-           enddo
-           tconvxl(kk)=tconvxl(kk)+fcvx
-           tconvyl(kk)=tconvyl(kk)+fcvy
-           tpresxl(kk)=tpresxl(kk)+fprx
-           tdiffxl(kk)=tdiffxl(kk)+fdix
-           tdiffyl(kk)=tdiffyl(kk)+fdiy
-        enddo
-     endif
-     !DC
-     if ((icvrt(iv).ge.ystart(1)).and.(icvrt(iv).le.yend(1))) then
-        i=icvrt(iv)-ystart(1)+1
-        do k=1,ysize(3)
-           kk=ystart(3)-1+k
-           fcvx=zero
-           fcvy=zero
-           fprx=zero
-           fdix=zero
-           fdiy=zero
-           do j=jcvlw_ly(iv),jcvup_ly(iv)-1
-              !momentum flux
-              uxmid = half*(ux2(i,j,k)+ux2(i,j+1,k))
-              uymid = half*(uy2(i,j,k)+uy2(i,j+1,k))
-              fcvx= fcvx +uxmid*uxmid*del_y(j)
-              fcvy= fcvy +uxmid*uymid*del_y(j)
-
-              !pressure
-              prmid=half*(ppi2(i,j,k)+ppi2(i,j+1,k))
-              fprx = fprx -prmid*del_y(j)
-
-              !viscous term
-              dudxmid = half*(ta2(i,j,k)+ta2(i,j+1,k))
-              dudymid = half*(tc2(i,j,k)+tc2(i,j+1,k))
-              dvdxmid = half*(tb2(i,j,k)+tb2(i,j+1,k))
-              fdix = fdix +two*xnu*dudxmid*del_y(j)
-              fdiy = fdiy +xnu*(dvdxmid+dudymid)*del_y(j)
-           enddo
-           tconvxl(kk)=tconvxl(kk)+fcvx
-           tconvyl(kk)=tconvyl(kk)+fcvy
-           tpresxl(kk)=tpresxl(kk)+fprx
-           tdiffxl(kk)=tdiffxl(kk)+fdix
-           tdiffyl(kk)=tdiffyl(kk)+fdiy
-        enddo
-     endif
-     call MPI_ALLREDUCE(tconvxl,tconvx,nz,real_type,MPI_SUM,MPI_COMM_WORLD,code)
-     call MPI_ALLREDUCE(tconvyl,tconvy,nz,real_type,MPI_SUM,MPI_COMM_WORLD,code)
-     call MPI_ALLREDUCE(tpresxl,tpresx,nz,real_type,MPI_SUM,MPI_COMM_WORLD,code)
-     call MPI_ALLREDUCE(tpresyl,tpresy,nz,real_type,MPI_SUM,MPI_COMM_WORLD,code)
-     call MPI_ALLREDUCE(tdiffxl,tdiffx,nz,real_type,MPI_SUM,MPI_COMM_WORLD,code)
-     call MPI_ALLREDUCE(tdiffyl,tdiffy,nz,real_type,MPI_SUM,MPI_COMM_WORLD,code)
-
-     do k=1,zsize(3)
-
-        tpresx(k)=tpresx(k)/dt
-        tpresy(k)=tpresy(k)/dt
-
-        xmom    = tunstx(k)+tconvx(k)
-        ymom    = tunsty(k)+tconvy(k)
-        xDrag(k) = two*(tdiffx(k)+tpresx(k)-xmom)
-        yLift(k) = two*(tdiffy(k)+tpresy(k)-ymom)
-
-     enddo
-
-     !Edited by F. Schuch
-     xDrag_mean = sum(xDrag(:))/real(nz,mytype)
-     yLift_mean = sum(yLift(:))/real(nz,mytype)
-
-!     if ((itime==ifirst).or.(itime==0)) then
-!        if (nrank .eq. 0) then
-!        write(filename,"('aerof',I1.1)") iv
-!        open(38+(iv-1),file=filename,status='unknown',form='formatted')
-!        endif
-!     endif
-     if (nrank .eq. 0) then
-        write(38,*) t,xDrag_mean,yLift_mean
-        call flush(38)
-     endif
-     if (mod(itime, icheckpoint).eq.0) then
-        if (nrank .eq. 0) then
-           write(filename,"('forces.dat',I7.7)") itime
-           call system("cp forces.dat " //filename)
-        endif
-     endif
-  enddo
-
-  do k = 1, xsize(3)
-     do j = 1, xsize(2)
-        do i = 1, xsize(1)
-           ux11(i,j,k)=ux01(i,j,k)
-           uy11(i,j,k)=uy01(i,j,k)
-           ux01(i,j,k)=ux1(i,j,k)
-           uy01(i,j,k)=uy1(i,j,k)
-        enddo
-     enddo
-  enddo
-
-  return
-
-end subroutine force
+       tconvxl=zero
+       tconvyl=zero
+       tdiffxl=zero
+       tdiffyl=zero
+       tpresxl=zero
+       tpresyl=zero
+       !BC and AD : x-pencils
+       !AD
+       if ((jcvlw(iv).ge.xstart(2)).and.(jcvlw(iv).le.xend(2))) then
+          j=jcvlw(iv)-xstart(2)+1
+          do k=1,xsize(3)
+             kk=xstart(3)-1+k
+             fcvx=zero
+             fcvy=zero
+             fpry=zero
+             fdix=zero
+             fdiy=zero
+             do i=icvlf_lx(iv),icvrt_lx(iv)-1
+                !momentum flux
+                uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k))
+                uymid = half*(uy1(i,j,k)+uy1(i+1,j,k))
+                fcvx  = fcvx -uxmid*uymid*dx
+                fcvy  = fcvy -uymid*uymid*dx
+
+                !pressure
+                prmid = half*(ppi1(i,j,k)+ppi1(i+1,j,k))
+                fpry  = fpry +prmid*dx
+
+                !viscous term
+                dudymid = half*(tc1(i,j,k)+tc1(i+1,j,k))
+                dvdxmid = half*(tb1(i,j,k)+tb1(i+1,j,k))
+                dvdymid = half*(td1(i,j,k)+td1(i+1,j,k))
+                fdix    = fdix -(xnu*(dudymid+dvdxmid)*dx)
+                fdiy    = fdiy -two*xnu*dvdymid*dx
+
+             enddo
+
+             tconvxl(kk)=tconvxl(kk)+fcvx
+             tconvyl(kk)=tconvyl(kk)+fcvy
+             tpresyl(kk)=tpresyl(kk)+fpry
+             tdiffxl(kk)=tdiffxl(kk)+fdix
+             tdiffyl(kk)=tdiffyl(kk)+fdiy
+          enddo
+       endif
+       !BC
+       if ((jcvup(iv).ge.xstart(2)).and.(jcvup(iv).le.xend(2))) then
+          j=jcvup(iv)-xstart(2)+1
+          do k=1,xsize(3)
+             kk=xstart(3)-1+k
+             fcvx=zero
+             fcvy=zero
+             fpry=zero
+             fdix=zero
+             fdiy=zero
+             do i=icvlf_lx(iv),icvrt_lx(iv)-1
+                !momentum flux
+                uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k))
+                uymid = half*(uy1(i,j,k)+uy1(i+1,j,k))
+                fcvx= fcvx +uxmid*uymid*dx
+                fcvy= fcvy +uymid*uymid*dx
+
+                !pressure
+                prmid = half*(ppi1(i,j,k)+ppi1(i+1,j,k))
+                fpry = fpry -prmid*dx
+
+                !viscous term
+                dudymid = half*(tc1(i,j,k)+tc1(i+1,j,k))
+                dvdxmid = half*(tb1(i,j,k)+tb1(i+1,j,k))
+                dvdymid = half*(td1(i,j,k)+td1(i+1,j,k))
+                fdix = fdix +(xnu*(dudymid+dvdxmid)*dx)
+                fdiy = fdiy +two*xnu*dvdymid*dx
+
+             enddo
+             tconvxl(kk)=tconvxl(kk)+fcvx
+             tconvyl(kk)=tconvyl(kk)+fcvy
+             tpresyl(kk)=tpresyl(kk)+fpry
+             tdiffxl(kk)=tdiffxl(kk)+fdix
+             tdiffyl(kk)=tdiffyl(kk)+fdiy
+          enddo
+       endif
+       !AB and DC : y-pencils
+       !AB
+       if ((icvlf(iv).ge.ystart(1)).and.(icvlf(iv).le.yend(1))) then
+          i=icvlf(iv)-ystart(1)+1
+          do k=1,ysize(3)
+             kk=ystart(3)-1+k
+             fcvx=zero
+             fcvy=zero
+             fprx=zero
+             fdix=zero
+             fdiy=zero
+             do j=jcvlw_ly(iv),jcvup_ly(iv)-1
+                !momentum flux
+                uxmid = half*(ux2(i,j,k)+ux2(i,j+1,k))
+                uymid = half*(uy2(i,j,k)+uy2(i,j+1,k))
+                fcvx= fcvx -uxmid*uxmid*del_y(j)
+                fcvy= fcvy -uxmid*uymid*del_y(j)
+
+                !pressure
+                prmid=half*(ppi2(i,j,k)+ppi2(i,j+1,k))
+                fprx = fprx +prmid*del_y(j)
+
+                !viscous term
+                dudxmid = half*(ta2(i,j,k)+ta2(i,j+1,k))
+                dudymid = half*(tc2(i,j,k)+tc2(i,j+1,k))
+                dvdxmid = half*(tb2(i,j,k)+tb2(i,j+1,k))
+                fdix = fdix -two*xnu*dudxmid*del_y(j)
+                fdiy = fdiy -xnu*(dvdxmid+dudymid)*del_y(j)
+             enddo
+             tconvxl(kk)=tconvxl(kk)+fcvx
+             tconvyl(kk)=tconvyl(kk)+fcvy
+             tpresxl(kk)=tpresxl(kk)+fprx
+             tdiffxl(kk)=tdiffxl(kk)+fdix
+             tdiffyl(kk)=tdiffyl(kk)+fdiy
+          enddo
+       endif
+       !DC
+       if ((icvrt(iv).ge.ystart(1)).and.(icvrt(iv).le.yend(1))) then
+          i=icvrt(iv)-ystart(1)+1
+          do k=1,ysize(3)
+             kk=ystart(3)-1+k
+             fcvx=zero
+             fcvy=zero
+             fprx=zero
+             fdix=zero
+             fdiy=zero
+             do j=jcvlw_ly(iv),jcvup_ly(iv)-1
+                !momentum flux
+                uxmid = half*(ux2(i,j,k)+ux2(i,j+1,k))
+                uymid = half*(uy2(i,j,k)+uy2(i,j+1,k))
+                fcvx= fcvx +uxmid*uxmid*del_y(j)
+                fcvy= fcvy +uxmid*uymid*del_y(j)
+
+                !pressure
+                prmid=half*(ppi2(i,j,k)+ppi2(i,j+1,k))
+                fprx = fprx -prmid*del_y(j)
+
+                !viscous term
+                dudxmid = half*(ta2(i,j,k)+ta2(i,j+1,k))
+                dudymid = half*(tc2(i,j,k)+tc2(i,j+1,k))
+                dvdxmid = half*(tb2(i,j,k)+tb2(i,j+1,k))
+                fdix = fdix +two*xnu*dudxmid*del_y(j)
+                fdiy = fdiy +xnu*(dvdxmid+dudymid)*del_y(j)
+             enddo
+             tconvxl(kk)=tconvxl(kk)+fcvx
+             tconvyl(kk)=tconvyl(kk)+fcvy
+             tpresxl(kk)=tpresxl(kk)+fprx
+             tdiffxl(kk)=tdiffxl(kk)+fdix
+             tdiffyl(kk)=tdiffyl(kk)+fdiy
+          enddo
+       endif
+       call MPI_ALLREDUCE(tconvxl,tconvx,nz,real_type,MPI_SUM,MPI_COMM_WORLD,code)
+       call MPI_ALLREDUCE(tconvyl,tconvy,nz,real_type,MPI_SUM,MPI_COMM_WORLD,code)
+       call MPI_ALLREDUCE(tpresxl,tpresx,nz,real_type,MPI_SUM,MPI_COMM_WORLD,code)
+       call MPI_ALLREDUCE(tpresyl,tpresy,nz,real_type,MPI_SUM,MPI_COMM_WORLD,code)
+       call MPI_ALLREDUCE(tdiffxl,tdiffx,nz,real_type,MPI_SUM,MPI_COMM_WORLD,code)
+       call MPI_ALLREDUCE(tdiffyl,tdiffy,nz,real_type,MPI_SUM,MPI_COMM_WORLD,code)
+
+       do k=1,zsize(3)
+
+          tpresx(k)=tpresx(k)/dt
+          tpresy(k)=tpresy(k)/dt
+
+          xmom    = tunstx(k)+tconvx(k)
+          ymom    = tunsty(k)+tconvy(k)
+          xDrag(k) = two*(tdiffx(k)+tpresx(k)-xmom)
+          yLift(k) = two*(tdiffy(k)+tpresy(k)-ymom)
+
+       enddo
+
+       !Edited by F. Schuch
+       xDrag_mean = sum(xDrag(:))/real(nz,mytype)
+       yLift_mean = sum(yLift(:))/real(nz,mytype)
+
+       !     if ((itime==ifirst).or.(itime==0)) then
+       !        if (nrank .eq. 0) then
+       !        write(filename,"('aerof',I1.1)") iv
+       !        open(38+(iv-1),file=filename,status='unknown',form='formatted')
+       !        endif
+       !     endif
+       if (nrank .eq. 0) then
+          write(38,*) t,xDrag_mean,yLift_mean
+          call flush(38)
+       endif
+       if (mod(itime, icheckpoint).eq.0) then
+          if (nrank .eq. 0) then
+             write(filename,"('forces.dat',I7.7)") itime
+             call system("cp forces.dat " //filename)
+          endif
+       endif
+    enddo
+
+    do k = 1, xsize(3)
+       do j = 1, xsize(2)
+          do i = 1, xsize(1)
+             ux11(i,j,k)=ux01(i,j,k)
+             uy11(i,j,k)=uy01(i,j,k)
+             ux01(i,j,k)=ux1(i,j,k)
+             uy01(i,j,k)=uy1(i,j,k)
+          enddo
+       enddo
+    enddo
+
+    return
+
+  end subroutine force
+end module forces