!> \file icetemp_mod.f90
!! Interfaces of functions and subroutines for temperature calculation 
!<

!> \namespace icetempmod
!! Interfaces of functions and subroutines for the temperature calculation 
!<

module icetempmod

contains

  !> SUBROUTINE: init_icetemp 
  !! Initialise the temperature calculation
  !! Used modules:
  !! - Icetemp_declar
  !>

  Subroutine init_icetemp(Num_rep_42_m)
    Use Icetemp_declar, only : Nfracq,Iq,Ecart_phid

    Implicit None

    !< Arguments
    Integer,                       Intent(In)::Num_rep_42_m !< Id Of Reponse
    ! Character(Len=80) :: Filin
    !________________________________________________________________
    ! La Temperature Du Point De Congelation De L'Eau De Mer
    ! Est Tiree De Jenkins (1991), Tb Est En Degres C
    ! _______________________________________________________________

    !  Nouvelle Formulation De Fracq
    !  Nfracq=15  ! Ancienne Formulation

    Nfracq=1

    Write(Num_rep_42_m,*) 'Calcul Des Temperatures'
    Write(Num_rep_42_m,*) 'Nfracq=',Nfracq

    Iq=1
    If (Iq.Eq.1) Write(Num_rep_42_m,*) 'Chaleur Demi Maille, Iq=',Iq
    If (Iq.Eq.2) Write(Num_rep_42_m,*) 'Chaleur Centree, Iq=',Iq
    If (Iq.Eq.3) Write(Num_rep_42_m,*) 'Chaleur Pente, Iq=',Iq
    If (Iq.Eq.4) Write(Num_rep_42_m,*) 'Chaleur U_taub, Iq=',Iq
    If (Iq.Eq.5) Write(Num_rep_42_m,*) 'Chaleur U_taub_stag, Iq=',Iq
    If (Iq.Eq.6) Write(Num_rep_42_m,*) 'Chaleur Q_sia_stag,Iq=',Iq
    If (Iq.Eq.7) Write(Num_rep_42_m,*) 'Chaleur Q_all_stag,Iq=',Iq
    Write(Num_rep_42_m,*)

    Ecart_phid=0.5
    Write(Num_rep_42_m,*) ' Si Base Froide, La Chaleur De Glissement Est En Exp(Delta T * Ecart_phid)'
    Write(Num_rep_42_m,*) 'Ecart_phid=',Ecart_phid
    Write(Num_rep_42_m,*) '---------------------------------------------------------------------------'

  End Subroutine init_icetemp
  

  !-----------------------------------------------------------------------------------------------------------------------

  !> SUBROUTINE: icetemp
  !! Calculate the temperature
  !!
  !! Used modules:
  !!  -  use Icetemp_declar 
  !>

  subroutine icetemp
    !$ use OMP_LIB 
    use Icetemp_declar


    Implicit None 
    
!~   integer :: t1,t2,ir
!~   real                           :: temps, t_cpu_0, t_cpu_1, t_cpu, norme  
!~   
!~    ! Temps CPU de calcul initial.
!~    call cpu_time(t_cpu_0)
!~    ! Temps elapsed de reference.
!~    call system_clock(count=t1, count_rate=ir)
    
    
	 !$OMP PARALLEL
	 !$OMP WORKSHARE
    ! init
    Aa=0.    
    Bb=0.              
    Cc=1.             
    Rr=0.              
    Hh=0.              
    Tdot=0
    Chal2_x=0.
    Chal2_y=0.
    Chal2_z=0.
    Chal2_xy=0.
    Chaldef_maj=0.
    Chalglissx=0.
    Chalglissy=0.
    
    !Instruction 
    Dee=1./(Nz-1) 
    Fracq=(1.-(1.-Dee/2.)**Nfracq)/Nfracq
    Fracq=(1.-(1.-Dee/2.)**5.)/5.
    Ee(1)=0.
    Ee(Nz)=1.
    !$OMP END WORKSHARE
    
    !$OMP DO
    Do K=1,Nz
       If ((K.Ne.1).And.(K.Ne.Nz)) Ee(K)=(K-1.)/(Nz-1.)
    End Do
    !$OMP END DO
    !$OMP END PARALLEL

    !   Variables Dependant Du Pas De Temps Dtt
    Ctm=Dtt*Cm/Rom/Cpm/Dzm/Dzm
    Dttdx=Dtt/Dx
    Dx11=1./Dx

    If (Itracebug.Eq.1) Write(num_tracebug,*)' Entree Dans Routine Ice_temp Time=',Time

    !________________________________________________________________
    ! La Temperature Du Point De Congelation De L'Eau De Mer
    ! Est Tiree De Jenkins (1991), Tb Est En Degres C
    ! _______________________________________________________________
    If (Itracebug.Eq.1) Write(num_tracebug,*)' avant appel  Routine Temp_mer'
    Call Temp_mer

    ! Appel Aux Proprietes Thermiques De La Glace : Conductivite, Capacite Claorifique -Tpmp

    If (Itracebug.Eq.1) Write(num_tracebug,*)' avant appel  Routine Thermal_property'
    Call Thermal_prop_icetemp

    If (Itracebug.Eq.1) Write(num_tracebug,*)' avant appel  Routine Qprod'
    Call Qprod_ice(Iq)

  
    ! Dans Le Socle : Les elements de La Matrice Tridiag Sont Invariants Dans L'Espace
    !$OMP PARALLEL
    !$OMP DO
    Do K=Nz+1,nz+nzm-1
       Aa(K)=-Ctm
       Bb(K)=1.+2.*Ctm
       Cc(K)=-Ctm
    End Do
	 !$OMP END DO

    ! Conditions Aux Limites Pour Le Socle
    Aa(nz+nzm)=-1.
    Bb(nz+nzm)=1.
    Cc(nz+nzm)=0.

    ! Calcul D'Advection : On Traite Les Tableaux En Une Seule Ligne De Commandes
    ! les tableaux iadvec valent 1 quand la vitesse consideree est upwind et 0 sinon

    !         Advection Selon X
    !         -----------------
    !$OMP WORKSHARE
    Iadvec_w(:,:)=Nint(Max(Sign(1.0,Uxbar(:,:)),0.0))
    Iadvec_e(:,:)=Nint(Abs(Min(Sign(1.0,Uxbar(:,:)),0.0)))

    !         Advection Selon Y
    !         -----------------

    Iadvec_s(:,:)=Nint(Max(Sign(1.0,Uybar(:,:)),0.0)) 
    Iadvec_n(:,:)=Nint(Abs(Min(Sign(1.0,Uybar(:,:)),0.0)))
    !$OMP END WORKSHARE
    ! Boucle De Resolution Des Temperatures Colonne Par Colonne
	 !$OMP SINGLE
    If (Itracebug.Eq.1) Write(num_tracebug,*)' avant appel  boucle Advec_icetemp'
    !$OMP END SINGLE

	 !$OMP DO PRIVATE(Deh22,Duu,Aa,Bb,Cc,Rr,Hh,Dou)
    Do J=2,Ny-1
       Do I=2,Nx-1

 !  If (Itracebug.Eq.1) Write(num_tracebug,*)' avant appel  Advect_icetemp  i,j', i,j
          call Advec_icetemp(Nz,Nzm,Nn,Chaldef_maj(I,J,:),Cp(I,J,:),& 
               Advecx(I,J,:),Advecy(I,J,:),Advec(I,J,:),Ct(I,J,:),Iadvec_w(I,J),Iadvec_e(I+1,J),Iadvec_s(I,J),Iadvec_n(I,J+1), &
					Ux(i,j,:),Ux(i+1,j,:),Uy(i,j,:),Uy(i,j+1,:),T(i,j,:),T(i-1,j,:),T(i+1,j,:),T(i,j-1,:),T(i,j+1,:), &
					H(i,j),Ts(i,j),Aa,Bb,Cc,Rr,Dx11,Dou,DTT,Dee,Uzr(i,j,:),Dttdx)

          ! variables de calcul dans la glace
          Deh22=2.*Dee*H(I,J)**2
          Duu=Dtt/2./Dee
          Rr(nz+nzm)=-Dzm*Ghf(I,J)/Cm  ! depend de I,J Pour Un Flux Geothermique Variable

!   If (Itracebug.Eq.1) Write(num_tracebug,*)' avant appel  temp_col  i,j', i,j
          Call temp_col(Nz,Nzm,Nn,T3d_new(i,j,:),Ct(i,j,:),&
               flot(i,j),H(i,j),Tbmer(i,j),Ibase(i,j),T(i,j,:),Tpmp(i,j,:),Aa,Bb,Cc,Rr,Hh,Ncond,Dee,Cm,&
               Dzm,Phid(i,j),Ghf(i,j),Ifail1,Tbdot(i,j),DTT,Dou,Ts(i,j))

          ! Test Permettant D'Ecrire Les Variables Quand Il Y A Une Erreur
          ! Ifail1 = 1 : Erreur Dans Tridiag
          If (Ifail1.Eq.1) Then
             Write(6,*)'Erreur Dans Tridiag'
             Write(6,*)'I = ',I,' J = ',J
             Write(6,*) 'A, B, C, R, U'
             Do K=1,nz+nzm
                Write(6,*)'K = ',K
                !Write(6,*) Aa(K),Bb(K),Cc(K),Rr(K),Hh(K)
             End Do
             Stop
          End If
       End Do
    End Do
	 !$OMP END DO

	 !$OMP DO
    Do J=1,Ny
       Do I=1,Nx
          H1(I,J)=H(I,J)
          B1(I,J)=B(I,J)
          Tpmp(I,J,1)=0.
       End Do
    End Do
    !$OMP END DO

    ! Attribution Finale De La Temperature
    ! Limitation a 3 deg De Variation Par Dtt Et Pas Plus Froid Que -70
	 !$OMP DO COLLAPSE(2) PRIVATE(Tdot)
    Do K=1,nz+nzm
       Do J=1,Ny
          Do I=1,Nx
             If (H(I,J).Gt.10.) Then

                Tdot(K)=T3d_new(I,J,K)-T(I,J,K)

                If ((K.Gt.1).And.(Tdot(K).Lt.-3.)) Tdot(K)=-3.
                If ((K.Gt.1).And.(Tdot(K).Gt.3.)) Tdot(K)=3.

                T(I,J,K)=T(I,J,K)+Tdot(K)

                If (T(I,J,K).Lt.-70.) Then
                   T(I,J,K)=-70.
                Endif
             Endif

          End Do
       End Do
    End Do
    !$OMP END DO

	 !$OMP DO COLLAPSE(2)
    Do K=2,Nz
       Do J=1,Ny
          Do I=1,Nx
             If (T(I,J,K).Gt.Tpmp(I,J,K)) Then
                T(I,J,K)=Tpmp(I,J,K)
                Ibase(I,J)=2
             Endif
          End Do
       End Do
    End Do
    !$OMP END DO

    !     Temperature Sur Les Bords :
    !     -------------------------
  
    !$OMP DO COLLAPSE(2)
    Do K=1,Nz
       Do J=1,Ny
          !        T(Nx,J,K)= T(Nx-1,J,K)
          !        T(1,J,K)= T(2,J,K) 
          T(Nx,J,K)= Ts(Nx,J)
          T(1,J,K)= Ts(1,J)
       End Do
       !
       Do I=1,Nx
          !        T(I,1,K)=T(I,2,K) 
          !        T(I,Ny,K)=T(I,Ny-1,K) 
          T(I,1,K)=Ts(I,1)
          T(I,Ny,K)=Ts(I,Ny)
       End Do
    End Do
    !$OMP END DO

    !$OMP DO COLLAPSE(2) 
    Do K=1,Nz
       Do J=1,Ny
          Do I=1,Nx
             If (H(I,J).le.1.) then 
                T(I,J,K)=Ts(I,J)
                if ((k.ge.nz).and.(flot(i,j))) T(I,J,K) =Tbmer(I,J)
             end If
                
          End Do
       End Do
    End Do
    !$OMP END DO
    !$OMP END PARALLEL
 
!~   ! Temps elapsed final
!~   call system_clock(count=t2, count_rate=ir)
!~   temps=real(t2 - t1,kind=4)/real(ir,kind=4)
!~   ! Temps CPU de calcul final
!~   call cpu_time(t_cpu_1)
!~   t_cpu = t_cpu_1 - t_cpu_0
  
!~   ! Impression du resultat.
!~   print '(//,3X,"Valeurs de nx et ny : ",I5,I5/,              &
!~            & 3X,"Temps elapsed       : ",1PE10.3," sec.",/, &
!~            & 3X,"Temps CPU           : ",1PE10.3," sec.",/, &
!~            & 3X,"Norme (PB si /= 0)  : ",1PE10.3,//)', &
!~            nx,ny,temps,t_cpu,norme


    If (Itracebug.Eq.1) Write(num_tracebug,*)' fin routine icetemp'
    return
  End Subroutine icetemp

  !-----------------------------------------------------------------------------------------------------------------------

end module icetempmod