module m_ECOSM_botstress
      contains
      subroutine ECOSM_botstress(n,tbs)
      !-----------------------------------------------------------------!
      ! Routine to calculate total bottom  stress, to be used in the    !
      ! sedimentation part of the NORWECOM ecosystem model.             !
      ! Last modified: 22.06.06, Cecilie Hansen                         ! 
      ! KAL -- changed to fixed-source format....
      !-----------------------------------------------------------------! 
      use mod_xc
      use mod_necessary_ecovars, only: bio_diagn2d,istress
      implicit none

      integer,intent(in) :: n

      real,intent(out),dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy):: tbs

      real    :: p_h,p_l,p_b,s_l,d_sl
      real    :: factor
      real    :: xxu,xxv,ro
      real    :: tbs_min,tbs_max,tbsx,tbsy
      integer :: i,j,k,k1

      include 'common_blocks.h'
   
      !KAL - Should use existing margin, points up to ii+margin should
      !KAL - still be valid at this point.
      !KAL - botstress uses hor. derivation so eff. margin is reduced
      !margin=0
      tbs_min=1000.
      tbs_max=-1000.
      d_sl=10. !depth boundary layer
      if(mnproc.eq.1)then
!       print*,'test inside botstress',n,margin
      endif

      do j=1-margin,jj+margin
      do i=1-margin,ii+margin
         
      if(ip(i,j)==1)then
         xxu=0.
         xxv=0.
         ro=0.
          p_b=p(i,j,kk+1)
          s_l=p_b-d_sl*onem
         do k=1,kdm
            p_h=p(i,j,k)
            p_l=p(i,j,k+1)
            if((p_h-s_l).gt.0.)then
               factor=1.*(p_l-p_h)/onem
            elseif((p_h-s_l).lt.0. .and.(p_l-s_l).gt.0.)then
               factor=(p_l-s_l)/onem
            else
               factor=0.
            endif
            xxu=xxu+(u(i,j,k,n)*real(iu(i,j))
     &     +u(i+1,j,k,n)*real(iu(i+1,j)))*factor
     &          /(max(real(iu(i,j)+iu(i+1,j)),1.))
            xxv=xxv+(v(i,j,k,n)*real(iv(i,j))
     &       +v(i,j+1,k,n)*real(iv(i,j+1)))*factor
     &          /(max(real(iv(i,j))+real(iv(i,j+1)),1.)) 
            ro=ro+(th3d(i,j,k,n)+thbase+1000)*factor
         enddo
         xxu=xxu/d_sl+(ubavg(i,j,n)*real(iu(i,j))
     &+ubavg(i+1,j,n)*real(iu(i+1,j)))
     &                 /(max(real(iu(i,j)+iu(i+1,j)),1.))

         xxv=xxv/d_sl+(vbavg(i,j,n)*real(iv(i,j))
     &+vbavg(i,j+1,n)*real(iv(i,j+1)))
     &                 /(max(real(iv(i,j)+iv(i,j+1)),1.))
            ro=ro/d_sl

         tbsx=cb*sqrt(xxu**2.+xxv**2.)*xxu              
         tbsy=cb*sqrt(xxu**2.+xxv**2.)*xxv              
         tbs(i,j)=ro*sqrt(tbsx**2.+tbsy**2.)              
         bio_diagn2d(i,j,istress)=tbs(i,j)

      endif
      enddo
      enddo

c     endif
      end subroutine ECOSM_botstress
      end module m_ECOSM_botstress