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