source: branches/GRISLIv3/SOURCES/deformation_mod_2lois_isotherme.f90 @ 425

!> \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, only:nx,ny,nz,e,num_param,num_rep_42,secyear,iglen
  use param_phy_mod, only: rog,ro,row
  
  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 

    integer :: k
    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 :: k
    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