#include "STREAMICE_OPTIONS.h"
#ifdef ALLOW_AUTODIFF
# include "AUTODIFF_OPTIONS.h"
#endif

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

CBOP
      SUBROUTINE STREAMICE_VISC_BETA_HYBRID ( myThid )

C     /============================================================\
C     | SUBROUTINE                                                 |
C     | o                                                          |
C     |============================================================|
C     |                                                            |
C     \============================================================/
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "GRID.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "STREAMICE.h"
#include "STREAMICE_CG.h"
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
#endif

C     !INPUT/OUTPUT ARGUMENTS
      INTEGER myThid
!       _RL taudx (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
!       _RL taudx (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)

#ifdef ALLOW_STREAMICE
#ifdef STREAMICE_HYBRID_STRESS

C     LOCAL VARIABLES
      INTEGER i, j, bi, bj, k, l, m
      INTEGER ikey_1
      _RL ux, uy, vx, vy, exx, eyy, exy, unorm, second_inv
      _RL ub, vb, fb, mean_u_shear, mean_v_shear, umid, vmid
      _RL omega_temp (Nr+1), u_shear(Nr+1), v_shear(Nr+1)
#ifdef STREAMICE_FLOWLINE_BUTTRESS
      _RL buttr_param, pwr
#endif

      _RL STREAMICE_BSTRESS_EXP
!       _RL total_vol_out
      external STREAMICE_BSTRESS_EXP

#ifdef STREAMICE_FLOWLINE_BUTTRESS
      buttr_param = 5**(1./3.) / (streamice_buttr_width/2.)**(4./3.)
      pwr = 1./n_glen
#endif

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1,sNy
         DO i=1,sNx
          IF (STREAMICE_hmask(i,j,bi,bj).eq.1) THEN

           umid = 0
           vmid = 0

           DO k=0,1
            DO l=0,1
             umid = umid + 0.25 *
     &        dxG(i,j+l,bi,bj)*dyG(i+k,j,bi,bj) *
     &        recip_rA(i,j,bi,bj) *
     &        U_streamice(i+k,j+l,bi,bj)
             vmid = vmid + 0.25 *
     &        dxG(i,j+l,bi,bj)*dyG(i+k,j,bi,bj) *
     &        recip_rA(i,j,bi,bj) *
     &        V_streamice(i+k,j+l,bi,bj)
            ENDDO
           ENDDO

           ux = (U_streamice(i+1,j+1,bi,bj) +
     &           U_streamice(i+1,j,bi,bj)   -
     &           U_streamice(i,j+1,bi,bj)   -
     &           U_streamice(i,j,bi,bj)) /
     &           (2. * dxF(i,j,bi,bj))
           vx = (V_streamice(i+1,j+1,bi,bj) +
     &           V_streamice(i+1,j,bi,bj)   -
     &           V_streamice(i,j+1,bi,bj)   -
     &           V_streamice(i,j,bi,bj)) /
     &           (2. * dxF(i,j,bi,bj))
           uy = (U_streamice(i+1,j+1,bi,bj) -
     &           U_streamice(i+1,j,bi,bj)   +
     &           U_streamice(i,j+1,bi,bj)   -
     &           U_streamice(i,j,bi,bj)) /
     &           (2. * dyF(i,j,bi,bj))
           vy = (V_streamice(i+1,j+1,bi,bj) -
     &           V_streamice(i+1,j,bi,bj)   +
     &           V_streamice(i,j+1,bi,bj)   -
     &           V_streamice(i,j,bi,bj)) /
     &           (2. * dyF(i,j,bi,bj))

           exx = ux + k1AtC_str(i,j,bi,bj)*vmid
           eyy = vy + k2AtC_str(i,j,bi,bj)*umid
           exy = .5*(uy+vx) +
     &      k1AtC_str(i,j,bi,bj)*umid + k2AtC_str(i,j,bi,bj)*vmid

           visc_streamice (i,j,bi,bj) = 0.0
           streamice_omega(i,j,bi,bj) = 0.0
           omega_temp (Nr+1) = 0.0
           u_shear(Nr+1) = 0.0
           v_shear(Nr+1) = 0.0

           DO m=Nr,1,-1

#ifdef ALLOW_AUTODIFF_TAMC
          act1 = bi - myBxLo(myThid)
          max1 = myBxHi(myThid) - myBxLo(myThid) + 1
          act2 = bj - myByLo(myThid)
          max2 = myByHi(myThid) - myByLo(myThid) + 1
          act3 = myThid - 1
          max3 = nTx*nTy
          act4 = ikey_dynamics - 1

          ikey_1 = m
     &         + Nr*(i-1)
     &         + Nr*sNx*(j-1)
     &         + Nr*sNx*sNy*act1
     &         + Nr*sNx*sNy*max1*act2
     &         + Nr*sNx*sNy*max1*max2*act3
     &         + Nr*sNx*sNy*max1*max2*max3*act4

CADJ STORE visc_streamice_full(i,j,m,bi,bj)
CADJ &     = comlev1_stream_hybrid, key=ikey_1
#endif

            streamice_vert_shear_uz (m) = streamice_taubx(i,j,bi,bj) /
     &       visc_streamice_full(i,j,m,bi,bj)
     &       * streamice_sigma_coord(m)

            streamice_vert_shear_vz (m) = streamice_tauby(i,j,bi,bj) /
     &       visc_streamice_full(i,j,m,bi,bj)
     &       * streamice_sigma_coord(m)

            second_inv =
     &       sqrt(exx**2+eyy**2+exx*eyy+exy**2+eps_glen_min**2+
     &            0.25 * streamice_vert_shear_uz(m)**2 +
     &            0.25 * streamice_vert_shear_vz(m)**2)

#ifdef STREAMICE_3D_GLEN_CONST
#if (defined (ALLOW_STREAMICE_OAD_FP))
            visc_full_new_si(i,j,m,bi,bj) =
#else
            visc_streamice_full(i,j,m,bi,bj) =
#endif
     &       .5 * B_glen(i,j,m,bi,bj)**2 *
     &        second_inv**((1-n_glen)/n_glen)
#else
#if (defined (ALLOW_STREAMICE_OAD_FP))
            visc_full_new_si(i,j,m,bi,bj) =
#else
            visc_streamice_full(i,j,m,bi,bj) =
#endif
     &       .5 * B_glen(i,j,bi,bj)**2 *
     &        second_inv**((1-n_glen)/n_glen)
#endif

            visc_streamice (i,j,bi,bj) = visc_streamice (i,j,bi,bj) +
     &       H_streamice(i,j,bi,bj) * streamice_delsigma (m) *
#if (defined (ALLOW_STREAMICE_OAD_FP))
     &       visc_full_new_si(i,j,m,bi,bj)
#else
     &       visc_streamice_full(i,j,m,bi,bj)
#endif

            omega_temp (m) = omega_temp(m+1) +
     &       streamice_sigma_coord(m) * streamice_delsigma(m) /
#if (defined (ALLOW_STREAMICE_OAD_FP))
     &       visc_full_new_si(i,j,m,bi,bj)
#else
     &       visc_streamice_full(i,j,m,bi,bj)
#endif

            u_shear (m) = u_shear (m+1) +
     &       streamice_vert_shear_uz (m) * streamice_delsigma (m) *
     &       H_streamice(i,j,bi,bj)

            v_shear (m) = v_shear (m+1) +
     &       streamice_vert_shear_vz (m) * streamice_delsigma (m) *
     &       H_streamice(i,j,bi,bj)

           ENDDO

           mean_u_shear = 0.0
           mean_v_shear = 0.0

           DO m=Nr,1,-1

            streamice_omega(i,j,bi,bj) = streamice_omega(i,j,bi,bj) +
     &       streamice_delsigma(m)*(omega_temp(m)+omega_temp(m+1))*.5
     &       * H_streamice(i,j,bi,bj)**2

            mean_u_shear = mean_u_shear +
     &       streamice_delsigma(m)*(u_shear(m)+u_shear(m+1))*.5

            mean_v_shear = mean_v_shear +
     &       streamice_delsigma(m)*(v_shear(m)+v_shear(m+1))*.5

           ENDDO

           streamice_u_surf(i,j,bi,bj) =
     &      u_shear(1) + umid - mean_u_shear

           streamice_v_surf(i,j,bi,bj) =
     &      v_shear(1) + vmid - mean_v_shear

           ub = umid - streamice_taubx(i,j,bi,bj) *
     &      streamice_omega(i,j,bi,bj) / H_streamice(i,j,bi,bj)

           streamice_u_bed (i,j,bi,bj) = ub

           vb = vmid - streamice_tauby(i,j,bi,bj) *
     &      streamice_omega(i,j,bi,bj) / H_streamice(i,j,bi,bj)

           streamice_v_bed (i,j,bi,bj) = vb

           unorm = sqrt(ub**2+vb**2+eps_u_min**2)

           fb = C_basal_friction(i,j,bi,bj)**2 *
     &      STREAMICE_BSTRESS_EXP (unorm,n_basal_friction) *
     &      streamice_basal_geom(i,j,bi,bj) *
     &      float_frac_streamice(i,j,bi,bj)

           tau_beta_eff_streamice(i,j,bi,bj) =
     &       fb /
     &       (1+fb*streamice_omega(i,j,bi,bj)/H_streamice(i,j,bi,bj))
#ifdef STREAMICE_FLOWLINE_BUTTRESS
           if (usestreamiceflowlineButtr) then
           tau_beta_eff_streamice(i,j,bi,bj) = 
     &      tau_beta_eff_streamice(i,j,bi,bj) + 
     &      buttr_param*B_glen(i,j,bi,bj)**2 * 
     &      STREAMICE_BSTRESS_EXP (unorm,pwr)*
     &      H_streamice(i,j,bi,bj)
          endif
#endif


          ENDIF
         ENDDO
        ENDDO
       ENDDO
      ENDDO
#ifdef ALLOW_OPENAD
#ifdef STREAMICE_HYBRID_STRESS
      tau_beta_eff_streamice(1,1,1,1) = 
     &  tau_beta_eff_streamice(1,1,1,1) +
     &  0. * streamice_u_surf(1,1,1,1) + 
     &  0. * streamice_v_surf(1,1,1,1)  
#endif
#endif

#endif
#endif
      RETURN
      END
