#include "options.inc"

c=======================================================================
c      Biharmonic friction and diffusion
c=======================================================================

      subroutine biha_hfric_u (diff_fe,diff_fn)
c---------------------------------------------------------------------------------
c      horizontal biharmonic friction for zonal momentum
c      no slip condition is possible
c---------------------------------------------------------------------------------
      use cpflame_module
      implicit none
      integer :: i,j,k,js,je
      real :: diff_fe(imt,jmt,km), diff_fn(imt,jmt,km) 
      real :: del2(imt,jmt,km),diff_tx,diff_ty,diffx
      DIFF_Tx(i,j,k) = (diff_fe(i,j,k)-diff_fe(i-1,j,k))/dx
      DIFF_Ty(i,j,k) = (diff_fn(i,j,k)-diff_fn(i,j-1,k))/dx

      js=max(2,js_pe); je = min(je_pe,jmt-1)
      diffx = sqrt(abs(ahbi))
      diff_fe(:,js_pe:je_pe,:)=0.0;diff_fn(:,js_pe:je_pe,:)=0.0
      do k=1,km
       do j=js,je
        do i=1,imt-1
         diff_fe(i,j,k)=diffx*(u(i+1,j,k,1,taum1)-u(i,j,k,1,taum1))/dx
        enddo
       enddo
      enddo
      call border_exchg3D(diff_fe,1); call setcyclic3D(diff_fe)
      do k=1,km-1
       do j=js,je
        do i=1,imt
         diff_fn(i,j,k)=diffx*(u(i,j+1,k,1,taum1)-u(i,j,k,1,taum1))/dx
     &                  *maskU(i,j+1,k)*maskU(i,j,k)
        enddo
       enddo
      enddo
      if (enable_noslip) then
       do j=js-1,je
         diff_fn(:,j,:)=diff_fn(:,j,:)-2*diffx*u(:,j,:,1,taum1)/dx
     &                   *(1-maskU(:,j+1,:))*maskU(:,j,:)
     &                                +2*diffx*u(:,j+1,:,1,taum1)/dx
     &                   *(1-maskU(:,j,:))*maskU(:,j+1,:)
       enddo
      endif
      call border_exchg3D(diff_fn,1); call setcyclic3D(diff_fn)
      del2(:,js_pe:je_pe,:)=0.0
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
	  del2(i,j,k) = (DIFF_Tx(i,j,k) + DIFF_Ty(i,j,k))*maskU(i,j,k)
        enddo
       enddo
      enddo
      call border_exchg3D(del2,1); call setcyclic3D(del2)
      diff_fe(:,js_pe:je_pe,:)=0.0;diff_fn(:,js_pe:je_pe,:)=0.0
      do k=2,km-1
       do j=js,je
        do i=1,imt-1
         diff_fe(i,j,k) =diffx*(del2(i+1,j,k)-del2(i,j,k))/dx
        enddo
       enddo
      enddo
      call border_exchg3D(diff_fe,1); call setcyclic3D(diff_fe)
      do k=1,km-1
       do j=js-1,je
        do i=2,imt-1
         diff_fn(i,j,k) = diffx*(del2(i,j+1,k) - del2(i,j,k))/dx
     &                     *maskU(i,j+1,k)*maskU(i,j,k)
        enddo
       enddo
      enddo
      if (enable_noslip) then
       do j=js-1,je
         diff_fn(:,j,:)=diff_fn(:,j,:)-2*diffx*del2(:,j,:)/dx
     &                   *(1-maskU(:,j+1,:))*maskU(:,j,:)
     &                                +2*diffx*del2(:,j+1,:)/dx
     &                   *(1-maskU(:,j,:))*maskU(:,j+1,:)
       enddo
      endif
      call border_exchg3D(diff_fn,1); call setcyclic3D(diff_fn)
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
         fu(i,j,k)=fu(i,j,k)
     &        -maskU(i,j,k)*(DIFF_Tx(i,j,k)+DIFF_Ty(i,j,k))
        enddo
       enddo
      enddo
      end subroutine biha_hfric_u




      subroutine biha_hfric_v (diff_fe,diff_fn)
c---------------------------------------------------------------------------------
c      horizontal biharmonic friction for meridional momentum
c      no slip condition is possible
c---------------------------------------------------------------------------------
      use cpflame_module
      implicit none
      integer :: i,j,k,js,je
      real :: diff_fe(imt,jmt,km), diff_fn(imt,jmt,km) 
      real :: del2(imt,jmt,km),diff_tx,diff_ty,diffx
      DIFF_Tx(i,j,k)=(diff_fe(i,j,k)-diff_fe(i-1,j,k))/dx
      DIFF_Ty(i,j,k)=(diff_fn(i,j,k)-diff_fn(i,j-1,k))/dx

      js=max(2,js_pe); je = min(je_pe,jmt-1)
      diffx = sqrt(abs(ahbi))
      diff_fe(:,js_pe:je_pe,:)=0.0;diff_fn(:,js_pe:je_pe,:)=0.0
      do k=1,km
       do j=js,je
        do i=1,imt-1
         diff_fe(i,j,k)=diffx*(u(i+1,j,k,2,taum1)-u(i,j,k,2,taum1))/dx
     &                  *maskV(i,j,k)*maskV(i+1,j,k)
        enddo
       enddo
      enddo
      if (enable_noslip) then
       do j=js,je
       do i=1,imt-1
         diff_fe(i,j,:)=diff_fe(i,j,:)-2*diffx*u(i,j,:,2,taum1)/dx
     &                   *(1-maskV(i+1,j,:))*maskV(i,j,:)
     &                                +2*diffx*u(i+1,j,:,2,taum1)/dx
     &                   *(1-maskV(i,j,:))*maskV(i+1,j,:)
       enddo
       enddo
      endif
      call border_exchg3D(diff_fe,1); call setcyclic3D(diff_fe)
      do k=1,km-1
       do j=js,je
        do i=1,imt
         diff_fn(i,j,k)=diffx*(u(i,j+1,k,2,taum1)-u(i,j,k,2,taum1))/dx
        enddo
       enddo
      enddo
      call border_exchg3D(diff_fn,1); call setcyclic3D(diff_fn)
      del2(:,je_pe:je_pe,:)=0.0
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
	  del2(i,j,k) = (DIFF_Tx(i,j,k) + DIFF_Ty(i,j,k))*maskV(i,j,k)
        enddo
       enddo
      enddo
      call border_exchg3D(del2,1); call setcyclic3D(del2)
      diff_fe(:,js_pe:je_pe,:)=0.0;diff_fn(:,js_pe:je_pe,:)=0.0
      do k=2,km-1
       do j=js,je
        do i=1,imt-1
         diff_fe(i,j,k) =diffx*(del2(i+1,j,k)-del2(i,j,k))/dx
     &                     *maskV(i+1,j,k)*maskV(i,j,k)
        enddo
       enddo
      enddo
      if (enable_noslip) then
       do j=js,je
       do i=1,imt-1
         diff_fe(i,j,:)=diff_fe(i,j,:)-2*diffx*del2(i,j,:)/dx
     &                   *(1-maskV(i+1,j,:))*maskV(i,j,:)
     &                                +2*diffx*del2(i+1,j,:)/dx
     &                   *(1-maskV(i,j,:))*maskV(i+1,j,:)
       enddo
       enddo
      endif
      call border_exchg3D(diff_fe,1); call setcyclic3D(diff_fe)
      do k=1,km-1
       do j=js,je
        do i=2,imt-1
         diff_fn(i,j,k) = diffx*(del2(i,j+1,k) - del2(i,j,k))/dx
        enddo
       enddo
      enddo
      call border_exchg3D(diff_fn,1); call setcyclic3D(diff_fn)
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
          fv(i,j,k)= fv(i,j,k)+maskV(i,j,k)*( 
     &               -DIFF_Tx(i,j,k) - DIFF_Ty(i,j,k))
        enddo
       enddo
      enddo
      end subroutine biha_hfric_v





      subroutine biha_hfric_w (diff_fe,diff_fn)
c---------------------------------------------------------------------------------
c      horizontal biharmonic friction for vertical momentum
c---------------------------------------------------------------------------------
      use cpflame_module
      implicit none
      integer :: i,j,k,js,je
      real :: diff_fe(imt,jmt,km), diff_fn(imt,jmt,km) 
      real :: del2(imt,jmt,km),diff_tx,diff_ty,diffx
      DIFF_Tx(i,j,k)=(diff_fe(i,j,k)-diff_fe(i-1,j,k))/dx
      DIFF_Ty(i,j,k)=(diff_fn(i,j,k)-diff_fn(i,j-1,k))/dx

      js=max(2,js_pe); je = min(je_pe,jmt-1)
      diffx = sqrt(abs(ahbi))
      diff_fe(:,js_pe:je_pe,:)=0.0;diff_fn(:,js_pe:je_pe,:)=0.0
      do k=1,km
       do j=js,je
        do i=1,imt-1
         diff_fe(i,j,k)=diffx*(u(i+1,j,k,3,taum1)-u(i,j,k,3,taum1))/dx
     &                  *maskW(i,j,k)*maskW(i+1,j,k)
        enddo
       enddo
      enddo
      call border_exchg3D(diff_fe,1); call setcyclic3D(diff_fe)
      do k=1,km-1
       do j=js,je
        do i=1,imt
         diff_fn(i,j,k)=diffx*(u(i,j+1,k,3,taum1)-u(i,j,k,3,taum1))/dx
     &                  *maskW(i,j,k)*maskW(i,j+1,k)
        enddo
       enddo
      enddo
      call border_exchg3D(diff_fn,1); call setcyclic3D(diff_fn)
      del2(:,je_pe:je_pe,:)=0.0
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
	  del2(i,j,k) = (DIFF_Tx(i,j,k) + DIFF_Ty(i,j,k))*maskW(i,j,k)
        enddo
       enddo
      enddo
      call border_exchg3D(del2,1); call setcyclic3D(del2)
      diff_fe(:,js_pe:je_pe,:)=0.0;diff_fn(:,js_pe:je_pe,:)=0.0
      do k=2,km-1
       do j=js,je
        do i=1,imt-1
         diff_fe(i,j,k) =diffx*(del2(i+1,j,k)-del2(i,j,k))/dx
     &                     *maskW(i+1,j,k)*maskW(i,j,k)
        enddo
       enddo
      enddo
      call border_exchg3D(diff_fe,1); call setcyclic3D(diff_fe)
      do k=1,km-1
       do j=js,je
        do i=2,imt-1
         diff_fn(i,j,k) = diffx*(del2(i,j+1,k) - del2(i,j,k))/dx
     &                     *maskW(i,j+1,k)*maskW(i,j,k)
        enddo
       enddo
      enddo
      call border_exchg3D(diff_fn,1); call setcyclic3D(diff_fn)
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
          fw(i,j,k)= fw(i,j,k)+maskW(i,j,k)*( 
     &               -DIFF_Tx(i,j,k) - DIFF_Ty(i,j,k))
        enddo
       enddo
      enddo
      end subroutine biha_hfric_w



      subroutine biha_vfric (diff_ft,mask,ff,n)
c---------------------------------------------------------------------------------
c       vertical biharmonic friction for  momentum component n
c       result is added to ff
c---------------------------------------------------------------------------------
      use cpflame_module
      implicit none
      integer :: i,j,k,js,je,n
      real,dimension(imt,jmt,km) :: diff_ft,mask,ff
      real :: del2(imt,jmt,km),diff_tz,diffx
      DIFF_Tz(i,j,k)=(diff_ft(i,j,k)-diff_ft(i,j,k-1))/dz

      js=max(2,js_pe); je = min(je_pe,jmt-1)
      diffx = sqrt(abs(Avbi))
      diff_ft(:,je_pe:je_pe,:)=0.0
      do k=1,km-1
       do j=js,je
        do i=2,imt-1
         diff_ft(i,j,k)=diffx*(u(i,j,k+1,n,taum1)-u(i,j,k,n,taum1))/dz
     &                  *mask(i,j,k)*mask(i,j,k+1)
        enddo
       enddo
      enddo
      del2(:,je_pe:je_pe,:)=0.0
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
	  del2(i,j,k) = DIFF_Tz(i,j,k)*mask(i,j,k)
        enddo
       enddo
      enddo
      diff_ft(:,je_pe:je_pe,:)=0.0
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
         diff_ft(i,j,k) =diffx*(del2(i,j,k+1)-del2(i,j,k))/dz
     &                     *mask(i,j,k+1)*mask(i,j,k)
        enddo
       enddo
      enddo
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
          ff(i,j,k)= ff(i,j,k)-mask(i,j,k)*DIFF_Tz(i,j,k) 
        enddo
       enddo
      enddo
      end subroutine biha_vfric



      subroutine biha_mix (diff_fe,diff_fn,diff_ft,var)
c---------------------------------------------------------------------------------
c      horizontal biharmonic diffusion of buoyancy
c      also vertical biha. diff. if requested
c---------------------------------------------------------------------------------
      use cpflame_module
      implicit none
      integer :: i,j,k,js,je
      real :: var(imt,jmt,km,0:2)
      real :: diff_fe(imt,jmt,km), diff_fn(imt,jmt,km) 
      real :: diff_ft(imt,jmt,km)  ,diffz
      real :: del2(imt,jmt,km),diff_tx,diff_ty,diffx,diff_tz
      DIFF_Tx(i,j,k)=(diff_fe(i,j,k)-diff_fe(i-1,j,k))/dx
      DIFF_Ty(i,j,k)=(diff_fn(i,j,k)-diff_fn(i,j-1,k))/dx
      DIFF_Tz(i,j,k)=(diff_ft(i,j,k)-diff_ft(i,j,k-1))/dz

      js=max(2,js_pe); je = min(je_pe,jmt-1)
      diffx = sqrt(abs(Khbi))
      diffz = sqrt(abs(Kvbi))
      diff_fe(:,js_pe:je_pe,:)=0.0; diff_fn(:,js_pe:je_pe,:)=0.0
      do k=1,km
       do j=js,je
        do i=1,imt-1
         diff_fe(i,j,k)=diffx*(var(i+1,j,k,taum1)-var(i,j,k,taum1))/dx
     &                  *maskU(i,j,k)
        enddo
       enddo
      enddo
      call border_exchg3D(diff_fe,1); call setcyclic3D(diff_fe)
      do k=1,km-1
       do j=js,je
        do i=1,imt
         diff_fn(i,j,k)=diffx*(var(i,j+1,k,taum1)-var(i,j,k,taum1))/dx
     &                  *maskV(i,j,k)
        enddo
       enddo
      enddo
      call border_exchg3D(diff_fn,1); call setcyclic3D(diff_fn)
      del2(:,js_pe:je_pe,:)=0.0
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
	  del2(i,j,k) = (DIFF_Tx(i,j,k) + DIFF_Ty(i,j,k))*maskT(i,j,k)
        enddo
       enddo
      enddo
      call border_exchg3D(del2,1); call setcyclic3D(del2)
      diff_fe(:,js_pe:je_pe,:)=0.0; diff_fn(:,js_pe:je_pe,:)=0.0
      do k=2,km-1
       do j=js,je
        do i=1,imt-1
         diff_fe(i,j,k) =diffx*(del2(i+1,j,k)-del2(i,j,k))/dx
     &                     *maskU(i,j,k)
        enddo
       enddo
      enddo
      call border_exchg3D(diff_fe,1); call setcyclic3D(diff_fe)
      do k=1,km-1
       do j=js,je
        do i=2,imt-1
         diff_fn(i,j,k) = diffx*(del2(i,j+1,k) - del2(i,j,k))/dx
     &                     *maskV(i,j,k)
        enddo
       enddo
      enddo
      call border_exchg3D(diff_fn,1); call setcyclic3D(diff_fn)
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
          var(i,j,k,taup1)= var(i,j,k,taup1)+c2dt*maskT(i,j,k)*( 
     &               -DIFF_Tx(i,j,k) - DIFF_Ty(i,j,k))
        enddo
       enddo
      enddo
      if (enable_vert_biha_diffusion) then
       diff_ft(:,js_pe:je_pe,:)=0.0
       do k=1,km-1
        do j=js,je
         do i=2,imt-1
          diff_ft(i,j,k)=diffz*(var(i,j,k+1,taum1)-var(i,j,k,taum1))/dz
     &                  *maskW(i,j,k)
         enddo
        enddo
       enddo
       del2(:,js_pe:je_pe,:)=0.0
       do k=2,km-1
        do j=js,je
         do i=2,imt-1
	  del2(i,j,k) = DIFF_Tz(i,j,k)*maskT(i,j,k)
         enddo
        enddo
       enddo
       diff_ft(:,js_pe:je_pe,:)=0.0
       do k=2,km-1
        do j=js,je
         do i=2,imt-1
          diff_ft(i,j,k) =diffz*(del2(i,j,k+1)-del2(i,j,k))/dz
     &                     *maskW(i,j,k)
         enddo
        enddo
       enddo
       do k=2,km-1
        do j=js,je
         do i=2,imt-1
          var(i,j,k,taup1)=var(i,j,k,taup1)
     &              -c2dt*maskT(i,j,k)*DIFF_Tz(i,j,k) 
         enddo
        enddo
       enddo
      endif
      end subroutine biha_mix

