Changeset 391


Ignore:
Timestamp:
03/23/23 13:39:51 (13 months ago)
Author:
dumas
Message:

use only in deformation_mod_2lois_isoth : compile but not tested in a simulation

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/GRISLIv3/SOURCES/deformation_mod_2lois_isotherme.f90

    r4 r391  
    2727module deformation_mod_2lois_isoth 
    2828 
    29 use module3d_phy 
    30 implicit none 
     29  use module3d_phy, only:nx,ny,nz,e,num_param,num_rep_42,secyear,iglen 
     30  use param_phy_mod, only: rog,ro,row 
     31   
     32  implicit none 
    3133 
    3234 
    33 ! declarations ne dependant pas du nombre de lois 
     35  ! declarations ne dependant pas du nombre de lois 
    3436 
    35 real, dimension(nx,ny,nz)   :: visc           !< viscosite de la glace (pour les shelves) 
    36 real, dimension(nz)         :: edecal         !< travail (decalage de E de 1 indice) 
     37  real, dimension(nx,ny,nz)   :: visc           !< viscosite de la glace (pour les shelves) 
     38  real, dimension(nz)         :: edecal         !< travail (decalage de E de 1 indice) 
    3739 
    38 ! declaration des lois 
    39 ! la valeur de n1poly et de n2poly determine le nombre de lois additionnees 
     40  ! declaration des lois 
     41  ! la valeur de n1poly et de n2poly determine le nombre de lois additionnees 
    4042 
    41 integer, parameter :: n1poly=1     !< 2 lois numerotees de 1 a 2 
    42 integer, parameter :: n2poly=2 
     43  integer, parameter :: n1poly=1     !< 2 lois numerotees de 1 a 2 
     44  integer, parameter :: n2poly=2 
    4345 
    4446 
    45 ! Tous les parametres de la loi sous forme de tableau (n1poly:n2poly) 
     47  ! Tous les parametres de la loi sous forme de tableau (n1poly:n2poly) 
    4648 
    47 real, dimension(n1poly:n2poly) :: glen             !< l'exposant 
    48 real, dimension(n1poly:n2poly) :: sf               !< softening factor for flow law 
    49 real, dimension(n1poly:n2poly) :: B_isoth          !< coefficient 
     49  real, dimension(n1poly:n2poly) :: glen             !< l'exposant 
     50  real, dimension(n1poly:n2poly) :: sf               !< softening factor for flow law 
     51  real, dimension(n1poly:n2poly) :: B_isoth          !< coefficient 
    5052 
    5153 
    5254 
    5355 
    54 ! declaration des tableaux utilises dans le calcul des vitesses, viscosites, ... 
     56  ! declaration des tableaux utilises dans le calcul des vitesses, viscosites, ... 
    5557 
    56 real, dimension(nx,ny,nz,n1poly:n2poly) :: Btt     !< coefficient au point considere 
    57 real, dimension(nx,ny,nz,n1poly:n2poly) :: Sa      !< sert dans l'integration des vitesses 
    58 real, dimension(nx,ny,nz,n1poly:n2poly) :: Sa_mx   !< Sa sur demi mailles > 
    59 real, dimension(nx,ny,nz,n1poly:n2poly) :: Sa_my   !< Sa sur demi mailles ^ 
    60 real, dimension(nx,ny,nz,n1poly:n2poly) :: S2a     !< sert dans l'integration des vitesses 
    61 real, dimension(nx,ny,nz,n1poly:n2poly) :: S2a_mx  !< S2a sur demi mailles > 
    62 real, dimension(nx,ny,nz,n1poly:n2poly) :: S2a_my  !< S2a sur demi mailles ^ 
    63 real, dimension(nx,ny,n1poly:n2poly)    :: ddx     !< sert dans le calcul de conserv. masse 
    64 real, dimension(nx,ny,n1poly:n2poly)    :: ddy     !< sert dans le calcul de conserv. masse 
     58  real, dimension(nx,ny,nz,n1poly:n2poly) :: Btt     !< coefficient au point considere 
     59  real, dimension(nx,ny,nz,n1poly:n2poly) :: Sa      !< sert dans l'integration des vitesses 
     60  real, dimension(nx,ny,nz,n1poly:n2poly) :: Sa_mx   !< Sa sur demi mailles > 
     61  real, dimension(nx,ny,nz,n1poly:n2poly) :: Sa_my   !< Sa sur demi mailles ^ 
     62  real, dimension(nx,ny,nz,n1poly:n2poly) :: S2a     !< sert dans l'integration des vitesses 
     63  real, dimension(nx,ny,nz,n1poly:n2poly) :: S2a_mx  !< S2a sur demi mailles > 
     64  real, dimension(nx,ny,nz,n1poly:n2poly) :: S2a_my  !< S2a sur demi mailles ^ 
     65  real, dimension(nx,ny,n1poly:n2poly)    :: ddx     !< sert dans le calcul de conserv. masse 
     66  real, dimension(nx,ny,n1poly:n2poly)    :: ddy     !< sert dans le calcul de conserv. masse 
    6567 
    66 real                                    :: de= 1./(nz-1) 
     68  real                                    :: de= 1./(nz-1) 
    6769 
    6870 
    … …  
    7173contains 
    7274 
    73 !------------------------------------------------------------------------------------------- 
     75  !------------------------------------------------------------------------------------------- 
    7476 
    75 !> SUBROUTINE: init_deformation 
    76 !! Routine qui lit les variables de deformation  
    77 !> 
    78 subroutine init_deformation  
     77  !> SUBROUTINE: init_deformation 
     78  !! Routine qui lit les variables de deformation  
     79  !> 
     80  subroutine init_deformation  
    7981 
    80 real :: exposant_1 
    81 real :: enhanc_fact_1 
    82 real :: coef_iso_1 
     82    integer :: k 
     83    real :: exposant_1 
     84    real :: enhanc_fact_1 
     85    real :: coef_iso_1 
    8386 
    84 real :: exposant_2 
    85 real :: enhanc_fact_2 
    86 real :: coef_iso_2 
     87    real :: exposant_2 
     88    real :: enhanc_fact_2 
     89    real :: coef_iso_2 
    8790 
    8891 
    89 namelist/loidef_1_iso/exposant_1,enhanc_fact_1,coef_iso_1 
    90 namelist/loidef_2_iso/exposant_2,enhanc_fact_2,coef_iso_2 
     92    namelist/loidef_1_iso/exposant_1,enhanc_fact_1,coef_iso_1 
     93    namelist/loidef_2_iso/exposant_2,enhanc_fact_2,coef_iso_2 
    9194 
    9295 
    93 ! loi 1 
    94 !-------------------------------------------------------------------------------- 
     96    ! loi 1 
     97    !-------------------------------------------------------------------------------- 
    9598 
    96 rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
    97 read(num_param,loidef_1_iso) 
     99    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
     100    read(num_param,loidef_1_iso) 
    98101 
    99 write(num_rep_42,*)'!___________________________________________________________'  
    100 write(num_rep_42,*)'! loi de deformation 1    module deformation_mod_2lois_isoth' 
    101 write(num_rep_42,*) 
    102 write(num_rep_42,loidef_1_iso) 
    103 write(num_rep_42,*)'! exposant (glen), coefficient et enhancement factor (sf)   ' 
    104 write(num_rep_42,*)'!___________________________________________________________'  
     102    write(num_rep_42,*)'!___________________________________________________________'  
     103    write(num_rep_42,*)'! loi de deformation 1    module deformation_mod_2lois_isoth' 
     104    write(num_rep_42,*) 
     105    write(num_rep_42,loidef_1_iso) 
     106    write(num_rep_42,*)'! exposant (glen), coefficient et enhancement factor (sf)   ' 
     107    write(num_rep_42,*)'!___________________________________________________________'  
    105108 
    106 glen(1)    = exposant_1  
    107 sf(1)      = enhanc_fact_1   
    108 B_isoth(1) = coef_iso_1*2.*secyear 
    109 write(6,*) 'loi de deformation B=',B_isoth(1) 
    110 write(6,*) 'gamma=', rog*(1.-ro/row) 
    111 write(6,*) 'epsilon = coef1 * H^3,   coef1 =',B_isoth(1)/2.*(rog*(1.-ro/row)/4.)**3 
    112 write(6,*) 'pvi = coef2 /H,          coef2 =',(4./(rog*(1.-ro/row)))**2 / B_isoth(1) 
     109    glen(1)    = exposant_1  
     110    sf(1)      = enhanc_fact_1   
     111    B_isoth(1) = coef_iso_1*2.*secyear 
     112    write(6,*) 'loi de deformation B=',B_isoth(1) 
     113    write(6,*) 'gamma=', rog*(1.-ro/row) 
     114    write(6,*) 'epsilon = coef1 * H^3,   coef1 =',B_isoth(1)/2.*(rog*(1.-ro/row)/4.)**3 
     115    write(6,*) 'pvi = coef2 /H,          coef2 =',(4./(rog*(1.-ro/row)))**2 / B_isoth(1) 
    113116 
    114 ! loi 2 
    115 !-------------------------------------------------------------------------------- 
    116    
    117 rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
    118 read(num_param,loidef_2_iso) 
     117    ! loi 2 
     118    !-------------------------------------------------------------------------------- 
    119119 
    120 write(num_rep_42,*)'!___________________________________________________________'  
    121 write(num_rep_42,*)'! loi de deformation 2    module deformation_mod_2lois_isoth' 
    122 write(num_rep_42,*) 
    123 write(num_rep_42,loidef_2_iso) 
    124 write(num_rep_42,*)'! exposant (glen), coefficient et enhancement factor (sf)   ' 
    125 write(num_rep_42,*)'!___________________________________________________________'  
     120    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
     121    read(num_param,loidef_2_iso) 
    126122 
    127 glen(2)    = exposant_2  
    128 sf(2)      = enhanc_fact_2   
    129 B_isoth(2) = coef_iso_2 
     123    write(num_rep_42,*)'!___________________________________________________________'  
     124    write(num_rep_42,*)'! loi de deformation 2    module deformation_mod_2lois_isoth' 
     125    write(num_rep_42,*) 
     126    write(num_rep_42,loidef_2_iso) 
     127    write(num_rep_42,*)'! exposant (glen), coefficient et enhancement factor (sf)   ' 
     128    write(num_rep_42,*)'!___________________________________________________________'  
    130129 
    131 ! calcul coordonnées reduites 
     130    glen(2)    = exposant_2  
     131    sf(2)      = enhanc_fact_2   
     132    B_isoth(2) = coef_iso_2 
    132133 
    133 e(1)=0. 
    134 e(nz)=1. 
    135 do k=1,nz 
    136    if ((k.ne.1).and.(k.ne.nz)) e(k)=(k-1.)/(nz-1.) 
    137 end do 
     134    ! calcul coordonnées reduites 
    138135 
    139 do iglen= n1poly,n2poly 
    140    call flowlaw(iglen) 
    141 end do 
     136    e(1)=0. 
     137    e(nz)=1. 
     138    do k=1,nz 
     139       if ((k.ne.1).and.(k.ne.nz)) e(k)=(k-1.)/(nz-1.) 
     140    end do 
    142141 
    143 ddx(:,:,:)=0. 
    144 ddy(:,:,:)=0. 
     142    do iglen= n1poly,n2poly 
     143       call flowlaw(iglen) 
     144    end do 
     145 
     146    ddx(:,:,:)=0. 
     147    ddy(:,:,:)=0. 
    145148 
    146149 
    147 end subroutine init_deformation 
     150  end subroutine init_deformation 
    148151 
    149152 
    150 !-------------------------------------------------------------------------------- 
    151 subroutine flow_general  ! fake 
    152 end subroutine flow_general 
     153  !-------------------------------------------------------------------------------- 
     154  subroutine flow_general  ! fake 
     155  end subroutine flow_general 
    153156 
    154 !-------------------------------------------------------------------------------- 
    155 ! application des sf et calcul des btt, Sa, ... 
    156 ! Tous independants de x et y 
     157  !-------------------------------------------------------------------------------- 
     158  ! application des sf et calcul des btt, Sa, ... 
     159  ! Tous independants de x et y 
    157160 
    158 subroutine flowlaw (iiglen) 
     161  subroutine flowlaw (iiglen) 
    159162 
    160   implicit none 
    161   integer ::  iiglen   !< compteur pour la boucle flowlaw 
    162   real    :: n1        !< glen+1 
    163   real    :: n2        !< glen+2 
     163    implicit none 
     164    integer :: k 
     165    integer ::  iiglen   !< compteur pour la boucle flowlaw 
     166    real    :: n1        !< glen+1 
     167    real    :: n2        !< glen+2 
    164168 
    165169 
    166170 
    167   B_isoth(iiglen)    =  B_isoth(iiglen)*sf(iiglen) 
    168   Btt(:,:,:,iiglen)  =  B_isoth(iiglen) 
    169   n1                 =  glen(iiglen)+1. 
    170   n2                 =  glen(iiglen)+2. 
    171   
     171    B_isoth(iiglen)    =  B_isoth(iiglen)*sf(iiglen) 
     172    Btt(:,:,:,iiglen)  =  B_isoth(iiglen) 
     173    n1                 =  glen(iiglen)+1. 
     174    n2                 =  glen(iiglen)+2. 
    172175 
    173   do k = 1,nz 
    174      Sa (:,:,k,iiglen)    =  B_isoth(iiglen)/n1 * (1.-e(k)**n1 ) 
    175      S2a(:,:,k,iiglen)    =  B_isoth(iiglen)/n2 * (1+ (e(k)**n2-n2*e(k))/n1 )  
    176   end do 
    177   Sa_mx(:,:,:,iiglen)  =  Sa (:,:,:,iiglen) 
    178   Sa_my(:,:,:,iiglen)  =  Sa (:,:,:,iiglen) 
    179   S2a_mx(:,:,:,iiglen) =  S2a (:,:,:,iiglen) 
    180   S2a_my(:,:,:,iiglen) =  S2a (:,:,:,iiglen) 
    181176 
    182 end subroutine flowlaw 
     177    do k = 1,nz 
     178       Sa (:,:,k,iiglen)    =  B_isoth(iiglen)/n1 * (1.-e(k)**n1 ) 
     179       S2a(:,:,k,iiglen)    =  B_isoth(iiglen)/n2 * (1+ (e(k)**n2-n2*e(k))/n1 )  
     180    end do 
     181    Sa_mx(:,:,:,iiglen)  =  Sa (:,:,:,iiglen) 
     182    Sa_my(:,:,:,iiglen)  =  Sa (:,:,:,iiglen) 
     183    S2a_mx(:,:,:,iiglen) =  S2a (:,:,:,iiglen) 
     184    S2a_my(:,:,:,iiglen) =  S2a (:,:,:,iiglen) 
     185 
     186  end subroutine flowlaw 
    183187 
    184188end module deformation_mod_2lois_isoth 
Note: See TracChangeset for help on using the changeset viewer.