!> \file advec_icetemp.f90
!! Subroutines for temperature advection in the heat transfer equation
!<

!> \namespace advec_icetemp
!! Subroutines to calculate the temperature advection 
!<

!> SUBROUTINE: Advec_icetemp
!! Implementation of Routines to calculate advection in the heat transfer equation
!! for a column
!! Used modules:
!!    - use Icetemp_declar
!>

Subroutine  Advec_icetemp(Nz,Nzm,Nn,Prodq_m,Cro_m,Advecx_m,Advecy_m,Advec_m,Ct_m,Iadvec_w,Iadvec_e,Iadvec_s,Iadvec_n, &
	Uxij,Uxi_plus_un,Uyij,Uyj_plus_un,Tij,Ti_moins_un,Ti_plus_un,Tj_moins_un,Tj_plus_un,Hij,Tsij,Aa,Bb,Cc,Rr, &
	Dx11,Dou,DTT,Dee,Uzrij,Dttdx)

!  use icetemp_declar 

  Implicit None
    !< Arguments
  Integer,                     intent(in)  ::  Nz,Nzm,Nn      !< Taille des tableaux
  Real, Dimension(Nz),         intent(in)   :: Prodq_m        !< Tableau 1d Vert. heat production
  Real, Dimension(Nz),         intent(in)   :: Cro_m          !< Tableau 1d Vert. heat capcity
  Real, Dimension(Nz),         intent(inout):: Advecx_m       !< Advection selon x
  Real, Dimension(Nz),         intent(inout):: Advecy_m       !< Advection selon y
  Real, Dimension(Nz),         intent(inout):: Advec_m        !< Advection total
  Real, Dimension(Nz),         intent(in)   :: Ct_m           !< Tableau 1d Vert.  thermal cond.
  
  Integer, intent(in) :: Iadvec_w,Iadvec_e,Iadvec_s,Iadvec_n
  real, dimension(Nz), Intent(in) :: Uxij, Uxi_plus_un
  real, dimension(Nz), Intent(in) :: Uyij, Uyj_plus_un
  real, dimension(Nz+Nzm), Intent(inout) :: Tij
  real, dimension(Nz+Nzm), Intent(in) :: Ti_moins_un,Ti_plus_un,Tj_moins_un,Tj_plus_un
  real, intent(in) :: Hij
  real,intent(in) :: Tsij
  Real,Dimension(Nn),intent(inout) :: Aa    !< Work Arrays For Tridiag  !Dim Nn
  Real,Dimension(Nn),intent(inout) :: Bb    !< Work Arrays For Tridiag  !Dim Nn
  Real,Dimension(Nn),intent(inout) :: Cc    !< Work Arrays For Tridiag  !Dim Nn
  Real,Dimension(Nn),intent(inout) :: Rr    !< Work Arrays For Tridiag  !Dim Nn
  Real, intent(in) :: Dx11
  Real, intent(inout) :: Dou
  real, intent(in) :: DTT
  real, intent(in) :: Dee
  real, dimension(Nz) :: Uzrij
  real :: Dttdx
  

  
  ! variables locales
  Integer :: K
  Real :: Dzz,Dah,Ct_haut,Ct_bas

 
  Do K=2,Nz-1


     ! On Remplit Deux Tableaux Pour eviter Des Problemes D'Asymetrie
     !         Avection Selon X  (advecx_m en general > 0)
     !         ----------------
     Advecx_m(K) = Iadvec_w * Uxij(K) *          & ! ux west if upwind
          (Tij(K) - Ti_moins_un(K))  & ! west T gradient
            
          + Iadvec_e * Uxi_plus_un(K)*      & ! ux east if upwind
          (Ti_plus_un(K)-Tij(K))     ! east T gradient

     Advecx_m(K) =  Advecx_m(K) * Dx11                                        !Dx11=1/Dx


     !         Avection Selon Y
     !         ----------------
     Advecy_m(K) =   Iadvec_s * Uyij(K) *        & ! uy sud if upwind
          (Tij(K) - Tj_moins_un(K))& ! south T gradient
          
          + Iadvec_n * Uyj_plus_un(K) *   & ! uy nord is upwind
          (Tj_plus_un(K)-Tij(K))    ! north T gradient

     Advecy_m(K) =  Advecy_m(K) * Dx11

     Advec_m(K)=Advecx_m(K)+Advecy_m(K)  ! Advection Totale

  End Do

  !     Fin Boucle K, Ecriture Matrice Tridiag Dans La Glace


  !      -----------------------------Cas General (H>10m)
  thick_ice:  If (Hij.Gt.10.) Then

     !        Variables De Calcul Dans La Glace
     ! dou = dtt/dz^2
     Dou=Dtt/Dee/Dee/ Hij / Hij  

     ! Dah = dtt/dz
     Dah=Dtt /Dee /Hij

     ! thermal conductivity at mid point just below the surface                   
     Ct_bas=2.*(Ct_m(1)*Ct_m(2))/(Ct_m(1)+Ct_m(2))           

     ! surface temperature : cannot go above 0 celsius
     Tij(1)=Min(0.,Tsij)

     ! surface boundary condition
     Aa(1)=0.
     Bb(1)=1.
     Cc(1)=0.
     Rr(1)=Tij(1)

     Do K=2,Nz-1
        Dzz=Dou/Cro_m(K)

        ! Conductivite Au Milieu Des Mailles
        Ct_haut=Ct_bas
        Ct_bas=2.*(Ct_m(K)*Ct_m(K+1))/(Ct_m(K)+Ct_m(K+1))

        ! Advection Verticale Centree
        Aa(K) = -Dzz * Ct_haut - Uzrij(K) * Dah / 2.  ! lower diag
        Cc(K) = -Dzz * Ct_bas  + Uzrij(K) * Dah / 2.  ! upper diag
        Bb(K) =  1.+ Dzz * (Ct_haut+Ct_bas)                                ! diag
        
        Rr(K) =  Tij(K) + Prodq_m(K) * Dtt   ! vector
        Rr(K) =  Rr(K)- Dttdx * (Advec_m(K))

     End DO
  End If thick_ice

End Subroutine Advec_icetemp