source: trunk/SOURCES/deformation_mod_2lois_isotherme.f90 @ 4

!> \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_isoth

use module3d_phy
implicit none


! declarations ne dependant pas du nombre de lois

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) :: sf               !< softening factor for flow law
real, dimension(n1poly:n2poly) :: B_isoth          !< coefficient




! 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

real                                    :: de= 1./(nz-1)




contains

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

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

real :: exposant_1
real :: enhanc_fact_1
real :: coef_iso_1

real :: exposant_2
real :: enhanc_fact_2
real :: coef_iso_2


namelist/loidef_1_iso/exposant_1,enhanc_fact_1,coef_iso_1
namelist/loidef_2_iso/exposant_2,enhanc_fact_2,coef_iso_2


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

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

write(num_rep_42,*)'!___________________________________________________________' 
write(num_rep_42,*)'! loi de deformation 1    module deformation_mod_2lois_isoth'
write(num_rep_42,*)
write(num_rep_42,loidef_1_iso)
write(num_rep_42,*)'! exposant (glen), coefficient et enhancement factor (sf)   '
write(num_rep_42,*)'!___________________________________________________________' 

glen(1)    = exposant_1 
sf(1)      = enhanc_fact_1  
B_isoth(1) = coef_iso_1*2.*secyear
write(6,*) 'loi de deformation B=',B_isoth(1)
write(6,*) 'gamma=', rog*(1.-ro/row)
write(6,*) 'epsilon = coef1 * H^3,   coef1 =',B_isoth(1)/2.*(rog*(1.-ro/row)/4.)**3
write(6,*) 'pvi = coef2 /H,          coef2 =',(4./(rog*(1.-ro/row)))**2 / B_isoth(1)

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

write(num_rep_42,*)'!___________________________________________________________' 
write(num_rep_42,*)'! loi de deformation 2    module deformation_mod_2lois_isoth'
write(num_rep_42,*)
write(num_rep_42,loidef_2_iso)
write(num_rep_42,*)'! exposant (glen), coefficient et enhancement factor (sf)   '
write(num_rep_42,*)'!___________________________________________________________' 

glen(2)    = exposant_2 
sf(2)      = enhanc_fact_2  
B_isoth(2) = coef_iso_2

! calcul coordonnées reduites

e(1)=0.
e(nz)=1.
do k=1,nz
   if ((k.ne.1).and.(k.ne.nz)) e(k)=(k-1.)/(nz-1.)
end do

do iglen= n1poly,n2poly
   call flowlaw(iglen)
end do

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


end subroutine init_deformation


!--------------------------------------------------------------------------------
subroutine flow_general  ! fake
end subroutine flow_general

!--------------------------------------------------------------------------------
! application des sf et calcul des btt, Sa, ...
! Tous independants de x et y

subroutine flowlaw (iiglen)

  implicit none
  integer ::  iiglen   !< compteur pour la boucle flowlaw
  real    :: n1        !< glen+1
  real    :: n2        !< glen+2



  B_isoth(iiglen)    =  B_isoth(iiglen)*sf(iiglen)
  Btt(:,:,:,iiglen)  =  B_isoth(iiglen)
  n1                 =  glen(iiglen)+1.
  n2                 =  glen(iiglen)+2.
 

  do k = 1,nz
     Sa (:,:,k,iiglen)    =  B_isoth(iiglen)/n1 * (1.-e(k)**n1 )
     S2a(:,:,k,iiglen)    =  B_isoth(iiglen)/n2 * (1+ (e(k)**n2-n2*e(k))/n1 ) 
  end do
  Sa_mx(:,:,:,iiglen)  =  Sa (:,:,:,iiglen)
  Sa_my(:,:,:,iiglen)  =  Sa (:,:,:,iiglen)
  S2a_mx(:,:,:,iiglen) =  S2a (:,:,:,iiglen)
  S2a_my(:,:,:,iiglen) =  S2a (:,:,:,iiglen)

end subroutine flowlaw

end module deformation_mod_2lois_isoth