source: trunk/SOURCES/Temperature-routines/temp_col.f90 @ 24

!> \file temp_col.f90
!! Subroutines for vertical temperature calculation 
!<

!> \namespace temp_col
!!  Subroutines for vertical temperature calculation 
!<

!> SUBROUTINE: Temp_column 
!! Subroutine for vertical temperature calculation
!! Used modules:
!! - use Icetemp_declar
!! - use Tridiagmod
!>

Subroutine  Temp_col(Ii,Jj,Newtempcol_m,Ct_m,Flot_m,H_m,Tbmer_m)

  Use Icetemp_declar
  Use Tridiagmod
  
  
  Implicit None
  integer :: Ii,Jj
  Real,Dimension(Nx,Ny,Nz+Nzm),Intent(Out)  :: Newtempcol_m        !< Tableau 1d Vert. Des Temperatures En Sortie
  Real,Dimension(Nz),Intent(In)             :: Ct_m                !< Tableau 1d Vert. Conductivite Thermique
  Logical,Intent(In)                        :: Flot_m              !< Vrai Si Flot_mtant (Test D'Archimede) 'O'
  Real,Intent(In)                           :: H_m                 !< Ice Thickness  'O'
  Real,Intent(In)                           :: Tbmer_m             !< Temperature De La Mer A La Base De L'Ice Shelf


  If (H_m.Gt.10.) Then

     ! Conditions Aux Limites a La Base De La Glace
     !________________________________________________

     ! Base Froide 
     If ( ((Ibase(Ii,Jj).Eq.1).Or.(Ibase(Ii,Jj).Eq.4) &
          .Or. ((Ibase(Ii,Jj).Eq.5).And.(T(Ii,Jj,Nz).Lt.Tpmp(Ii,Jj,Nz)))) &
          .And..Not.Flot_m) Then

        Ibase(Ii,Jj)=1
        Bb(Nz)=1.

        If (Ncond.Eq.1) Then ! Avec Socle
           Dzi=H_m*Dee*Cm/Ct_m(Nz)
           Aa(Nz)=-Dzm/(Dzm+Dzi)
           Cc(Nz)=-Dzi/(Dzm+Dzi)
           Rr(Nz)=H_m*Dee*Phid(Ii,Jj)*Dzm/Ct_m(Nz)/(Dzm+Dzi)
        Else                 ! Sans Socle
           Aa(Nz)=-1.
           Cc(Nz)=0.
           Rr(Nz)=-(Ghf(Ii,Jj)-Phid(Ii,Jj))/Ct_m(Nz)*H_m*Dee
        Endif

     Else

        ! Base Temperee Ou Shelf
        ! ----------------------               
        If (Ibase(Ii,Jj).Eq.5)  Ibase(Ii,Jj)=2
        Ibase(Ii,Jj)=Max(Ibase(Ii,Jj),2)
        Aa(Nz)=0.
        Bb(Nz)=1.
        Cc(Nz)=0.

        If (.Not.Flot_m) Then
           Rr(Nz)=Tpmp(Ii,Jj,Nz)
        Else
           Rr(Nz)=Tbmer_m
        End If

     Endif


     If (Ncond.Eq.1) Then ! Avec Socle
        Do K=Nz+1,nz+nzm-1
           Rr(K)=T(Ii,Jj,K)
        End Do
        Call Tridiag (Aa,Bb,Cc,Rr,Hh,nz+nzm,Ifail1)

     Else                 ! Sans Socle

        Call Tridiag (Aa,Bb,Cc,Rr,Hh,Nz,Ifail1)
     Endif


     ! Temperature Dans Le Socle Si Ncond=0, Lineaire Avec Le Gradient Ghf
     If (Ncond.Eq.0) Then
        Do K=Nz+1,nz+nzm
           Hh(K)=Hh(Nz)-Dzm*(K-Nz)*Ghf(Ii,Jj)/Cm
        End Do
     Endif

     Do K=1,nz+nzm
        Newtempcol_m(Ii,Jj,K)=Hh(K)
     End Do
     Tbdot(Ii,Jj)=(Newtempcol_m(Ii,Jj,Nz)-T(Ii,Jj,Nz))/Dtt

     !      ------------------------------------------------- Glace Tres Fine (H<10m Ou H=0)
  Else If (H_m.Le.10.) Then

     !       Pour Eviter Des Problemes Lorsque H Est Inferieur A 10 M
     !       Profil Lineaire, Pente=Flux Geothermique
     !       Le Cas Sans Glace Est Traite De La Meme Facon  (Dou=0)

     ! ........................................ Avec Calcul Dans Le Socle
     If (Ncond.Eq.1.And..Not.Flot_m) Then   
        If (H_m.Gt.0.) Then
           ! Gradient Dans Le Socle
           Dou=(T(Ii,Jj,Nz+1)-T(Ii,Jj,Nz))/Dzm*Cm     
           Dou=Dou/Ct_m(Nz)*Dee*H_m

           Tss=Min(0.,Ts(Ii,Jj))
           Do K=1,Nz
              Newtempcol_m(Ii,Jj,K)=Tss+Dou*(K-1.)
           End Do

        Else
           ! Pas De Glace
           Tss=Ts(Ii,Jj)
           Do K=1,Nz
              Newtempcol_m(Ii,Jj,K)=Tss
           End Do
        End If

        ! Calcul Dans Le Socle Meme S'Il N'Y A Pas De Glace
        Rr(Nz)=Newtempcol_m(Ii,Jj,Nz)
        Aa(Nz)=0
        Bb(Nz)=1
        Cc(Nz)=0

        Do K=Nz+1,nz+nzm-1
           Rr(K)=T(Ii,Jj,K)
        End Do

        ! Creation De Nouveaux Tableaux Juste Pour Tridiag      
        Do K=1,Nzm+1
           Abis(K)=Aa(Nz-1+K)
           Bbis(K)=Bb(Nz-1+K)
           Cbis(K)=Cc(Nz-1+K)
           Rbis(K)=Rr(Nz-1+K)
        End Do

        Call Tridiag (Abis,Bbis,Cbis,Rbis,Hbis,Nzm+1,Ifail1)

        Do K=1,Nzm+1
           Hh(Nz-1+K)=Hbis(K)
        End Do

        ! ........................................ Sans Calcul Dans Le Socle
     Else
        ! Calotte Posee                    
        If ((H_m.Gt.0.).And..Not.Flot_m) Then
           Dou=-Ghf(Ii,Jj)/Ct_m(Nz)*Dee*H_m
           Tss=Min(0.,Ts(Ii,Jj))
           Do K=1,Nz
              Newtempcol_m(Ii,Jj,K)=Tss+Dou*(K-1.)
           End Do

           !  Shelf
        Else If ((H_m.Gt.0.).And.Flot_m) Then
           Tss=Min(0.,Ts(Ii,Jj))
           Dou=(Tbmer_m-Tss)*Dee
           Do K=1,Nz
              Newtempcol_m(Ii,Jj,K)=Tss+Dou*(K-1.)
           End Do

           ! Pas De Glace
        Else
           Tss=Ts(Ii,Jj)      
           Do K=1,Nz
              Newtempcol_m(Ii,Jj,K)=Tss
           End Do
        End If

        ! Temperature Dans Le Socle, Lineaire Avec Le Gradient Ghf
        If (.Not.Flot_m) Then
           Do K=Nz+1,nz+nzm
              Hh(K)=Newtempcol_m(Ii,Jj,Nz)-Dzm*(K-Nz)*Ghf(Ii,Jj)/Cm
           End Do
        Else
           Do K=Nz+1,nz+nzm
              Hh(K)=Tbmer_m-Dzm*(K-Nz)*Ghf(Ii,Jj)/Cm
           End Do
        End If

     Endif
     ! Fin Du Test Socle


     Do K=Nz+1,nz+nzm
        Newtempcol_m(Ii,Jj,K)=Hh(K)
     End Do

     Tbdot(Ii,Jj)=(Newtempcol_m(Ii,Jj,Nz)-T(Ii,Jj,Nz))/Dtt
     Bmelt(Ii,Jj)=0.
     Ibase(Ii,Jj)=5
     Phid(Ii,Jj)=0.


  End If

End Subroutine Temp_col