!> \file deformation_mod_2lois.f90
!! C'est ce module qui doit etre selectionne pour faire le calcul de la deformation de la glace  
!! en utilisant une loi de deformation polynomiale
!! a deux composantes (habituellement n=3 + n=1) 
!!
!<

!> \namespace deformation_mod_2lois
!! C'est ce module qui doit etre selectionne pour faire le calcul de la
!! deformation de la glace en utilisant une loi de deformation polynomiale
!! a deux composantes (habituellement n=3 + n=1)
!! \author CatRitz
!! \date decmebre 2008
!! @note Used modules 
!! @note  - use module3D_phy
!! @note  Contient : 
!! @note  - les declarations des tableaux (etait avant dans deform_declar)
!! en fait la declaration est independante du nombre d'elements de la loi 
!! et peut etre simplement copie pour un autre nombre 
!! @note  - lecture des valeurs utilisees par namelist
!! N'est valable que pour la loi a 2 composants. Pour un nombre different de composants, 
!! le recopier et modifier 
!! @note  -  les parameters n1poly et n2poly 
!! @note  -  init_deformation en ajoutant ou supprimant des blocs de namelist
!! 
!<
module deformation_mod_2lois

use module3d_phy
implicit none


! declarations ne dependant pas du nombre de lois

real, dimension(nx,ny,nz)   :: teta           !< teta = t - tpmp
real, dimension(nx,ny,nz-1) :: ti2            !< calcule dans flow_general
real, dimension(nx,ny,nz)   :: visc           !< viscosite de la glace (pour les shelves)
real, dimension(nz)         :: edecal         !< travail (decalage de E de 1 indice)

! declaration des lois
! la valeur de n1poly et de n2poly determine le nombre de lois additionnees

integer, parameter :: n1poly=1     !< 2 lois numerotees de 1 a 2
integer, parameter :: n2poly=2


! Tous les parametres de la loi sous forme de tableau (n1poly:n2poly)

real, dimension(n1poly:n2poly) :: glen             !< l'exposant
real, dimension(n1poly:n2poly) :: ttrans           !< la temperature de transition 
real, dimension(n1poly:n2poly) :: sf               !< softening factor for flow law

! Pour chaque loi (meme exposant), il y a deux domaines avec une temperature de transition 

! 1- temperature inferieure a Ttrans
real, dimension(n1poly:n2poly) :: Q1               !< energies d'activation (temp < ttrans)  
real, dimension(n1poly:n2poly) :: Bat1             !< coefficient           (temp < ttrans)
         
! 2- temperature superieure a Ttrans
real, dimension(n1poly:n2poly) :: Q2               !< energies d'activation (temp > ttrans) 
real, dimension(n1poly:n2poly) :: Bat2             !< coefficient           (temp > ttrans) 


! declaration des tableaux utilises dans le calcul des vitesses, viscosites, ...

real, dimension(nx,ny,nz,n1poly:n2poly) :: Btt     !< coefficient au point considere
real, dimension(nx,ny,nz,n1poly:n2poly) :: Sa      !< sert dans l'integration des vitesses
real, dimension(nx,ny,nz,n1poly:n2poly) :: Sa_mx   !< Sa sur demi mailles >
real, dimension(nx,ny,nz,n1poly:n2poly) :: Sa_my   !< Sa sur demi mailles ^
real, dimension(nx,ny,nz,n1poly:n2poly) :: S2a     !< sert dans l'integration des vitesses
real, dimension(nx,ny,nz,n1poly:n2poly) :: S2a_mx  !< S2a sur demi mailles >
real, dimension(nx,ny,nz,n1poly:n2poly) :: S2a_my  !< S2a sur demi mailles ^
real, dimension(nx,ny,n1poly:n2poly) :: ddx     !< sert dans le calcul de conserv. masse
real, dimension(nx,ny,n1poly:n2poly) :: ddy     !< sert dans le calcul de conserv. masse

contains

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

!> SUBROUTINE: init_deformation
!! Routine qui lit les variables de deformation 
!>
subroutine init_deformation 

real :: exposant_1
real :: temp_trans_1
real :: enhanc_fact_1
real :: coef_cold_1
real :: Q_cold_1
real :: coef_warm_1
real :: Q_warm_1

real :: exposant_2
real :: temp_trans_2
real :: enhanc_fact_2
real :: coef_cold_2
real :: Q_cold_2
real :: coef_warm_2
real :: Q_warm_2

namelist/loidef_1/exposant_1,temp_trans_1,enhanc_fact_1,   & 
                  coef_cold_1,Q_cold_1,coef_warm_1,Q_warm_1

namelist/loidef_2/exposant_2,temp_trans_2,enhanc_fact_2,   & 
                  coef_cold_2,Q_cold_2,coef_warm_2,Q_warm_2


! loi 1
!--------------------------------------------------------------------------------

rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
read(num_param,loidef_1)
write(num_rep_42,loidef_1)

glen(1)   = exposant_1 
ttrans(1) = temp_trans_1
sf(1)     = enhanc_fact_1  
Bat1(1)   = coef_cold_1
Q1(1)     = Q_cold_1 
Bat2(1)   = coef_warm_1
Q2(1)     = Q_warm_1 


! loi 2
!--------------------------------------------------------------------------------
  
rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
read(num_param,loidef_2)
write(num_rep_42,loidef_2)

write(num_rep_42,*)'!___________________________________________________________' 
write(num_rep_42,*)'! loi de deformation            module deformation_mod_2lois'
write(num_rep_42,*)'! exposant (glen), temperature de transition (ttrans)'
write(num_rep_42,*)'! enhancement factor (sf)'
write(num_rep_42,*)'! pour les temperatures inf. a Temp_trans :'
write(num_rep_42,*)'!                            coef_cold (Bat1) et Q_cold (Q1)'
write(num_rep_42,*)'! pour les temperatures sup. a Temp_trans :'
write(num_rep_42,*)'!                            coef_warm (Bat2) et Q_warm (Q2)'
write(num_rep_42,*)'!___________________________________________________________' 

glen(2)   = exposant_2 
ttrans(2) = temp_trans_2
sf(2)     = enhanc_fact_2  
Bat1(2)   = coef_cold_2
Q1(2)     = Q_cold_2 
Bat2(2)   = coef_warm_2
Q2(2)     = Q_warm_2 

! autre parametres ne changeant pas d'un run a l'autre
RGAS=8.314

! application des sf

do iglen= n1poly,n2poly
   bat1(iglen)=bat1(iglen)*sf(iglen)
   bat2(iglen)=bat2(iglen)*sf(iglen)
end do

ddx(:,:,:)=0.
ddy(:,:,:)=0.


end subroutine init_deformation


!--------------------------------------------------------------------
subroutine flow_general

  !$ USE OMP_LIB
  implicit none

  !$OMP PARALLEL
  !$OMP WORKSHARE
  teta(:,:,:) = t(:,:,1:nz)-tpmp(:,:,:)
  !$OMP END WORKSHARE

  !$OMP DO COLLAPSE(2)
  do k=nz-1,1,-1
     do i=2,nx-1
        do j=2,ny-1
           ti2(i,j,k) = (273.15+(tpmp(i,j,k+1)+tpmp(i,j,k))/2.)* &
                (273.15+(t(i,j,k+1)+t(i,j,k))/2.)
        end do
     end do
  end do
  !$OMP END DO
  !$OMP END PARALLEL

end subroutine flow_general
!---------------------------------------------------------------------------------------
subroutine flowlaw (iiglen)

  !$ USE OMP_LIB

  implicit none
  integer ::  iiglen   !< compteur pour la boucle flowlaw
  real :: a5 !< exposant de la loi de def +2 
  real :: a4 !< exposant de la loi de def +1 
  real :: ti1 !< utilise pour le calcul de btt a k=1
  real :: bat !< prend la valeur de bat1 ou bat2 selon les cas
  real :: q !< prend la valeur de q1 ou q2 selon les cas
  real :: aat !< variable de travail =q*tetar(i,j,k)
  real :: ssec !< variable de travail 
  real :: zetat !< variable de travail : profondeur ou t = trans(iiglen)
  real,dimension(nx,ny,nz) :: si1 !< tableau de calcul
  real,dimension(nx,ny,nz) :: si2 !< tableau de calcul
  real,dimension(nx,ny) :: tab_mx
  real,dimension(nx,ny) :: tab_my
  real,dimension(nx,ny) :: tab2d

!  real,dimension(nz) :: e          ! vertical coordinate in ice, scaled to h zeta
  real :: de= 1./(nz-1)
! variables openmp :  
  !$  integer :: rang
  !$ integer, dimension(:), allocatable :: tab_sync
  !$ integer :: nb_taches


  e(1)=0.
  e(nz)=1.
  
  !$OMP PARALLEL PRIVATE(bat,q,aat,ssec,zetat)
  !$ nb_taches = OMP_GET_NUM_THREADS()
  !$OMP DO 
  do k=1,nz
     if ((k.ne.1).and.(k.ne.nz)) e(k)=(k-1.)/(nz-1.)
  end do
  !$OMP END DO NOWAIT

  ! exposant de la loi de deformation utilisee : glen(iiglen)
  ! l'exposant correspondant a iiglen est defini dans deformation_mod
  a5 = glen(iiglen) + 2  
  a4 = glen(iiglen) + 1

  !    boucle pour btt a k=1
  ti1=273.15*273.15
  !$OMP DO
  do j=2,ny-1
     do i=2,nx-1
        if (t(i,j,1).le.ttrans(iiglen)) then
           btt(i,j,1,iiglen)=bat1(iiglen)*exp(q1(iiglen)/rgas/ti1*t(i,j,1))
        else
           btt(i,j,1,iiglen)=bat2(iiglen)*exp(q2(iiglen)/rgas/ti1*t(i,j,1))
        endif
     end do
  end do
  !$OMP END DO

  ! boucle pour tous les autres btt
  !$OMP DO COLLAPSE(2)
  do k=nz-1,1,-1
	do j=2,ny-1  
        do i=2,nx-1  
           if (h(i,j).gt.1.) then
              if ((teta(i,j,k).le.ttrans(iiglen)).and. &
                   (teta(i,j,k+1).le.ttrans(iiglen))) then
                 bat=bat1(iiglen)
                 q=q1(iiglen)
                 aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1)
                 aat=max(-1.8,aat)
                 aat=min(80.,aat)
                 ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
                 si1(i,j,k)=ssec/(aat+a4)*(e(k)**(aat+a4)-e(k+1)**(aat+a4))
                 si2(i,j,k)=((e(k)**(aat+a5))/(aat+a5) &
                      -(e(k+1)**(aat+a5))/(aat+a5) &
                      -(e(k+1)**(aat+a4))*e(k)+e(k+1)**(aat+a5)) &
                      * ssec/(aat+a4) 
                 btt(i,j,k+1,iiglen)=bat*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
              else if ((teta(i,j,k).ge.ttrans(iiglen)).and. &
                   (teta(i,j,k+1).ge.ttrans(iiglen))) then
                 bat=bat2(iiglen)
                 q=q2(iiglen)
                 aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1)
                 aat=max(-1.8,aat)
                 aat=min(80.,aat)
                 ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
                 si1(i,j,k)=ssec/(aat+a4)*(e(k)**(aat+a4)-e(k+1)**(aat+a4))
                 si2(i,j,k)=((e(k)**(aat+a5))/(aat+a5) &
                      -(e(k+1)**(aat+a5))/(aat+a5) &
                      -(e(k+1)**(aat+a4))*e(k)+e(k+1)**(aat+a5)) &
                      * ssec/(aat+a4) 
                 btt(i,j,k+1,iiglen)=bat*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))

                 ! cas ou la temperature de la maille est en partie au dessus et au dessous
                 ! de ttrans(iiglen)
              else if ((teta(i,j,k).lt.ttrans(iiglen)).and. &
                   (teta(i,j,k+1).gt.ttrans(iiglen))) then
                 !         calcul de la profondeur de transition
                 if (abs(teta(i,j,k)-teta(i,j,k+1)).lt.1.e-5) then
                    zetat=(e(k)+e(k+1))/2.
                 else
                    zetat=e(k+1)-(ttrans(iiglen)-teta(i,j,k+1))/ &
                         (teta(i,j,k)-teta(i,j,k+1))*de
                 endif

                 !         integration entre zeta2 et zetat, loi bat2 
                 bat=bat2(iiglen)
                 q=q2(iiglen)
                 aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1)
                 aat=max(-1.8,aat)
                 aat=min(80.,aat)
                 ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
                 si1(i,j,k)=ssec/(aat+a4)*(zetat**(aat+a4)-e(k+1)**(aat+a4))
                 si2(i,j,k)=((zetat**(aat+a5))/(aat+a5) &
                      -(e(k+1)**(aat+a5))/(aat+a5) &
                      -(e(k+1)**(aat+a4))*zetat+e(k+1)**(aat+a5)) &
                      * ssec/(aat+a4)
                 btt(i,j,k+1,iiglen)=bat*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))

                 !         integration entre zetat et zeta1, loi bat1 
                 bat=bat1(iiglen)
                 q=q1(iiglen)
                 aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1)
                 aat=max(-1.8,aat)
                 aat=min(80.,aat)
                 ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
                 si1(i,j,k)=si1(i,j,k)+ &
                      ssec/(aat+a4)*(e(k)**(aat+a4)-zetat**(aat+a4))
                 si2(i,j,k)=si2(i,j,k) &
                      +((e(k)**(aat+a5))/(aat+a5) &
                      -(zetat**(aat+a5))/(aat+a5) &
                      -(zetat**(aat+a4))*e(k)+zetat**(aat+a5)) &
                      * ssec/(aat+a4)


                 ! deuxieme cas ou la temperature de la maille est en partie au dessus et 
                 ! au dessous de ttrans(iiglen)
              else if ((teta(i,j,k).gt.ttrans(iiglen)).and. &
                   (teta(i,j,k+1).lt.ttrans(iiglen))) then
                 !         calcul de la profondeur de transition
                 if (abs(teta(i,j,k)-teta(i,j,k+1)).lt.1.e-5) then
                    zetat=(e(k)+e(k+1))/2.
                 else
                    zetat=e(k+1)-(ttrans(iiglen)-teta(i,j,k+1))/ &
                         (teta(i,j,k)-teta(i,j,k+1))*de
                 endif

                 !         integration entre zeta2 et zetat, loi bat1 
                 bat=bat1(iiglen)
                 q=q1(iiglen)
                 aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1)
                 aat=max(-1.8,aat)
                 aat=min(80.,aat)
                 ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
                 si1(i,j,k)=ssec/(aat+a4)*(zetat**(aat+a4)-e(k+1)**(aat+a4))
                 si2(i,j,k)=((zetat**(aat+a5))/(aat+a5) &
                      -(e(k+1)**(aat+a5))/(aat+a5) &
                      -(e(k+1)**(aat+a4))*zetat+e(k+1)**(aat+a5)) &
                      * ssec/(aat+a4)
                 btt(i,j,k+1,iiglen)=bat*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))

                 !         integration entre zetat et zeta1, loi bat2 
                 bat=bat2(iiglen)
                 q=q2(iiglen)
                 aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1)
                 aat=max(-1.8,aat)
                 aat=min(80.,aat)
                 ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
                 si1(i,j,k)=si1(i,j,k)+ &
                      ssec/(aat+a4)*(e(k)**(aat+a4)-zetat**(aat+a4))
                 si2(i,j,k)=si2(i,j,k) &
                      +((e(k)**(aat+a5))/(aat+a5) &
                      -(zetat**(aat+a5))/(aat+a5) &
                      -(zetat**(aat+a4))*e(k)+zetat**(aat+a5)) &
                      * ssec/(aat+a4)
              endif
           endif
        end do
     end do
  end do
  !$OMP END DO NOWAIT
  
  !     integration de sa et s2a
  !$OMP DO
  do j=1,ny
     do i=1,nx  
        sa(i,j,nz,iiglen)=0.
        s2a(i,j,nz,iiglen)=0.
        if (h(i,j).le.1.) btt(i,j,nz,iiglen)=bat2(iiglen)
     end do
  end do
  !$OMP END DO
  !$OMP END PARALLEL
  

  ! Allocation et initialisation du tableau tab_sync
  ! qui gere la synchronisation entre les differents threads
  !$ allocate(tab_sync(0:nb_taches-1))
  !$ tab_sync(0:nb_taches-1) = 1
  !$OMP PARALLEL private(i,j,k,rang) shared(tab_sync)
  !$ rang = OMP_GET_THREAD_NUM()  
  do k=nz-1,1,-1
       ! Synchronisation des threads
     !$     if (rang /= 0) then
     !$        do
     !$OMP FLUSH(tab_sync) 
     !$           if (tab_sync(rang-1)>=tab_sync(rang)+1) exit
     !$        enddo
     !$OMP FLUSH(sa)
     !$OMP FLUSH(s2a)
     !$     endif
     !$OMP DO SCHEDULE(STATIC)
     do j=1,ny
        do i=1,nx
           if (h(i,j).gt.1.) then
              sa(i,j,k,iiglen)=sa(i,j,k+1,iiglen)-si1(i,j,k)
              s2a(i,j,k,iiglen)=s2a(i,j,k+1,iiglen)+si2(i,j,k)+ &
                   sa(i,j,k+1,iiglen)*de
           else
              sa(i,j,k,iiglen)=bat2(iiglen)/a4*(1-e(k)**a4)
              s2a(i,j,k,iiglen)=bat2(iiglen)/a4*(a4/a5-e(k)+e(k)**a5/a5)
              btt(i,j,k,iiglen)=bat2(iiglen)
           endif
        end do
     end do
     !$OMP END DO NOWAIT
!     ! Mis a jour du tableau tab_sync
     !$     tab_sync(rang) = tab_sync(rang) + 1
     !$OMP FLUSH(tab_sync,sa,s2a)
  end do

  !$OMP WORKSHARE
  !     cas particulier des bords
  sa(:,1,:,iiglen)=sa(:,2,:,iiglen)
  s2a(:,1,:,iiglen)=s2a(:,2,:,iiglen)
  btt(:,1,:,iiglen)=btt(:,2,:,iiglen)
  sa(:,ny,:,iiglen)=sa(:,ny-1,:,iiglen)
  s2a(:,ny,:,iiglen)=s2a(:,ny-1,:,iiglen)
  btt(:,ny,:,iiglen)=btt(:,ny-1,:,iiglen)

  sa(1,:,:,iiglen)=sa(2,:,:,iiglen)
  s2a(1,:,:,iiglen)=s2a(2,:,:,iiglen)
  btt(1,:,:,iiglen)=btt(2,:,:,iiglen)
  sa(nx,:,:,iiglen)=sa(nx-1,:,:,iiglen)
  s2a(nx,:,:,iiglen)=s2a(nx-1,:,:,iiglen)
  btt(nx,:,:,iiglen)=btt(nx-1,:,:,iiglen)
  !$OMP END WORKSHARE
  !$OMP END PARALLEL
  
  ! calcul des moyennes sur les mailles staggered
  do k=1,nz
     tab2d=sa(:,:,k,iiglen)

     call moy_mxmy(nx,ny,tab2d,tab_mx,tab_my)
     sa_mx(:,:,k,iglen)=tab_mx
     sa_my(:,:,k,iglen)=tab_my

     tab2d=s2a(:,:,k,iiglen)
     call moy_mxmy(nx,ny,tab2d,tab_mx,tab_my)
     s2a_mx(:,:,k,iglen)=tab_mx
     s2a_my(:,:,k,iglen)=tab_my
  end do
 
! attribution de debug_3d pour les sorties s2a
! loi 1 -> 190 dans description variable -> 90 dans debug_3d
!debug_3d(:,:,90+iiglen-1) = s2a(:,:,1,iiglen)

end subroutine flowlaw


end module deformation_mod_2lois