Changeset 390


Ignore:
Timestamp:
03/23/23 13:30:34 (13 months ago)
Author:
dumas
Message:

use only in module deformation_mod_2lois

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/GRISLIv3/SOURCES/deformation_mod_2lois.f90

    r224 r390  
    2727module deformation_mod_2lois 
    2828 
    29 use module3d_phy 
    30 implicit none 
    31  
    32  
    33 ! declarations ne dependant pas du nombre de lois 
    34  
    35 real, dimension(nx,ny,nz)   :: teta           !< teta = t - tpmp 
    36 real, dimension(nx,ny,nz-1) :: ti2            !< calcule dans flow_general 
    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) 
    39  
    40 ! declaration des lois 
    41 ! la valeur de n1poly et de n2poly determine le nombre de lois additionnees 
    42  
    43 integer, parameter :: n1poly=1     !< 2 lois numerotees de 1 a 2 
    44 integer, parameter :: n2poly=2 
    45  
    46  
    47 ! Tous les parametres de la loi sous forme de tableau (n1poly:n2poly) 
    48  
    49 real, dimension(n1poly:n2poly) :: glen             !< l'exposant 
    50 real, dimension(n1poly:n2poly) :: ttrans           !< la temperature de transition  
    51 real, dimension(n1poly:n2poly) :: sf               !< softening factor for flow law 
    52  
    53 ! Pour chaque loi (meme exposant), il y a deux domaines avec une temperature de transition  
    54  
    55 ! 1- temperature inferieure a Ttrans 
    56 real, dimension(n1poly:n2poly) :: Q1               !< energies d'activation (temp < ttrans)   
    57 real, dimension(n1poly:n2poly) :: Bat1             !< coefficient           (temp < ttrans) 
    58           
    59 ! 2- temperature superieure a Ttrans 
    60 real, dimension(n1poly:n2poly) :: Q2               !< energies d'activation (temp > ttrans)  
    61 real, dimension(n1poly:n2poly) :: Bat2             !< coefficient           (temp > ttrans)  
    62  
    63  
    64 ! declaration des tableaux utilises dans le calcul des vitesses, viscosites, ... 
    65  
    66 real, dimension(nx,ny,nz,n1poly:n2poly) :: Btt     !< coefficient au point considere 
    67 real, dimension(nx,ny,nz,n1poly:n2poly) :: Sa      !< sert dans l'integration des vitesses 
    68 real, dimension(nx,ny,nz,n1poly:n2poly) :: Sa_mx   !< Sa sur demi mailles > 
    69 real, dimension(nx,ny,nz,n1poly:n2poly) :: Sa_my   !< Sa sur demi mailles ^ 
    70 real, dimension(nx,ny,nz,n1poly:n2poly) :: S2a     !< sert dans l'integration des vitesses 
    71 real, dimension(nx,ny,nz,n1poly:n2poly) :: S2a_mx  !< S2a sur demi mailles > 
    72 real, dimension(nx,ny,nz,n1poly:n2poly) :: S2a_my  !< S2a sur demi mailles ^ 
    73 real, dimension(nx,ny,n1poly:n2poly) :: ddx     !< sert dans le calcul de conserv. masse 
    74 real, dimension(nx,ny,n1poly:n2poly) :: ddy     !< sert dans le calcul de conserv. masse 
     29  use module3d_phy, only:nx,ny,nz,e,num_param,num_rep_42,rgas,T,iglen,tpmp,H 
     30 
     31  implicit none 
     32 
     33 
     34  ! declarations ne dependant pas du nombre de lois 
     35 
     36  real, dimension(nx,ny,nz)   :: teta           !< teta = t - tpmp 
     37  real, dimension(nx,ny,nz-1) :: ti2            !< calcule dans flow_general 
     38  real, dimension(nx,ny,nz)   :: visc           !< viscosite de la glace (pour les shelves) 
     39  real, dimension(nz)         :: edecal         !< travail (decalage de E de 1 indice) 
     40 
     41  ! declaration des lois 
     42  ! la valeur de n1poly et de n2poly determine le nombre de lois additionnees 
     43 
     44  integer, parameter :: n1poly=1     !< 2 lois numerotees de 1 a 2 
     45  integer, parameter :: n2poly=2 
     46 
     47 
     48  ! Tous les parametres de la loi sous forme de tableau (n1poly:n2poly) 
     49 
     50  real, dimension(n1poly:n2poly) :: glen             !< l'exposant 
     51  real, dimension(n1poly:n2poly) :: ttrans           !< la temperature de transition  
     52  real, dimension(n1poly:n2poly) :: sf               !< softening factor for flow law 
     53 
     54  ! Pour chaque loi (meme exposant), il y a deux domaines avec une temperature de transition  
     55 
     56  ! 1- temperature inferieure a Ttrans 
     57  real, dimension(n1poly:n2poly) :: Q1               !< energies d'activation (temp < ttrans)   
     58  real, dimension(n1poly:n2poly) :: Bat1             !< coefficient           (temp < ttrans) 
     59 
     60  ! 2- temperature superieure a Ttrans 
     61  real, dimension(n1poly:n2poly) :: Q2               !< energies d'activation (temp > ttrans)  
     62  real, dimension(n1poly:n2poly) :: Bat2             !< coefficient           (temp > ttrans)  
     63 
     64 
     65  ! declaration des tableaux utilises dans le calcul des vitesses, viscosites, ... 
     66 
     67  real, dimension(nx,ny,nz,n1poly:n2poly) :: Btt     !< coefficient au point considere 
     68  real, dimension(nx,ny,nz,n1poly:n2poly) :: Sa      !< sert dans l'integration des vitesses 
     69  real, dimension(nx,ny,nz,n1poly:n2poly) :: Sa_mx   !< Sa sur demi mailles > 
     70  real, dimension(nx,ny,nz,n1poly:n2poly) :: Sa_my   !< Sa sur demi mailles ^ 
     71  real, dimension(nx,ny,nz,n1poly:n2poly) :: S2a     !< sert dans l'integration des vitesses 
     72  real, dimension(nx,ny,nz,n1poly:n2poly) :: S2a_mx  !< S2a sur demi mailles > 
     73  real, dimension(nx,ny,nz,n1poly:n2poly) :: S2a_my  !< S2a sur demi mailles ^ 
     74  real, dimension(nx,ny,n1poly:n2poly) :: ddx     !< sert dans le calcul de conserv. masse 
     75  real, dimension(nx,ny,n1poly:n2poly) :: ddy     !< sert dans le calcul de conserv. masse 
    7576 
    7677contains 
    7778 
    78 !------------------------------------------------------------------------------------------- 
    79  
    80 !> SUBROUTINE: init_deformation 
    81 !! Routine qui lit les variables de deformation  
    82 !> 
    83 subroutine init_deformation  
    84  
    85 real :: exposant_1 
    86 real :: temp_trans_1 
    87 real :: enhanc_fact_1 
    88 real :: coef_cold_1 
    89 real :: Q_cold_1 
    90 real :: coef_warm_1 
    91 real :: Q_warm_1 
    92  
    93 real :: exposant_2 
    94 real :: temp_trans_2 
    95 real :: enhanc_fact_2 
    96 real :: coef_cold_2 
    97 real :: Q_cold_2 
    98 real :: coef_warm_2 
    99 real :: Q_warm_2 
    100  
    101 namelist/loidef_1/exposant_1,temp_trans_1,enhanc_fact_1,   &  
    102                   coef_cold_1,Q_cold_1,coef_warm_1,Q_warm_1 
    103  
    104 namelist/loidef_2/exposant_2,temp_trans_2,enhanc_fact_2,   &  
    105                   coef_cold_2,Q_cold_2,coef_warm_2,Q_warm_2 
    106  
    107  
    108 ! loi 1 
    109 !-------------------------------------------------------------------------------- 
    110  
    111 rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
    112 read(num_param,loidef_1) 
    113 write(num_rep_42,loidef_1) 
    114  
    115 glen(1)   = exposant_1  
    116 ttrans(1) = temp_trans_1 
    117 sf(1)     = enhanc_fact_1   
    118 Bat1(1)   = coef_cold_1 
    119 Q1(1)     = Q_cold_1  
    120 Bat2(1)   = coef_warm_1 
    121 Q2(1)     = Q_warm_1  
    122  
    123  
    124 ! loi 2 
    125 !-------------------------------------------------------------------------------- 
    126    
    127 rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
    128 read(num_param,loidef_2) 
    129 write(num_rep_42,loidef_2) 
    130  
    131 write(num_rep_42,*)'!___________________________________________________________'  
    132 write(num_rep_42,*)'! loi de deformation            module deformation_mod_2lois' 
    133 write(num_rep_42,*)'! exposant (glen), temperature de transition (ttrans)' 
    134 write(num_rep_42,*)'! enhancement factor (sf)' 
    135 write(num_rep_42,*)'! pour les temperatures inf. a Temp_trans :' 
    136 write(num_rep_42,*)'!                            coef_cold (Bat1) et Q_cold (Q1)' 
    137 write(num_rep_42,*)'! pour les temperatures sup. a Temp_trans :' 
    138 write(num_rep_42,*)'!                            coef_warm (Bat2) et Q_warm (Q2)' 
    139 write(num_rep_42,*)'!___________________________________________________________'  
    140  
    141 glen(2)   = exposant_2  
    142 ttrans(2) = temp_trans_2 
    143 sf(2)     = enhanc_fact_2   
    144 Bat1(2)   = coef_cold_2 
    145 Q1(2)     = Q_cold_2  
    146 Bat2(2)   = coef_warm_2 
    147 Q2(2)     = Q_warm_2  
    148  
    149 ! autre parametres ne changeant pas d'un run a l'autre 
    150 RGAS=8.314 
    151  
    152 ! application des sf 
    153  
    154 do iglen= n1poly,n2poly 
    155    bat1(iglen)=bat1(iglen)*sf(iglen) 
    156    bat2(iglen)=bat2(iglen)*sf(iglen) 
    157 end do 
    158  
    159 ddx(:,:,:)=0. 
    160 ddy(:,:,:)=0. 
    161  
    162  
    163 end subroutine init_deformation 
    164  
    165  
    166 !-------------------------------------------------------------------- 
    167 subroutine flow_general 
    168  
    169   !$ USE OMP_LIB 
    170   implicit none 
    171  
    172   !$OMP PARALLEL 
    173   !$OMP WORKSHARE 
    174   teta(:,:,:) = t(:,:,1:nz)-tpmp(:,:,:) 
    175   !$OMP END WORKSHARE 
    176  
    177   !$OMP DO COLLAPSE(2) 
    178   do k=nz-1,1,-1 
    179      do i=2,nx-1 
    180         do j=2,ny-1 
    181            ti2(i,j,k) = (273.15+(tpmp(i,j,k+1)+tpmp(i,j,k))/2.)* & 
    182                 (273.15+(t(i,j,k+1)+t(i,j,k))/2.) 
    183         end do 
    184      end do 
    185   end do 
    186   !$OMP END DO 
    187   !$OMP END PARALLEL 
    188  
    189 end subroutine flow_general 
    190 !--------------------------------------------------------------------------------------- 
    191 subroutine flowlaw (iiglen) 
    192  
    193   !$ USE OMP_LIB 
    194  
    195   implicit none 
    196   integer ::  iiglen   !< compteur pour la boucle flowlaw 
    197   real :: a5 !< exposant de la loi de def +2  
    198   real :: a4 !< exposant de la loi de def +1  
    199   real :: ti1 !< utilise pour le calcul de btt a k=1 
    200   real :: bat !< prend la valeur de bat1 ou bat2 selon les cas 
    201   real :: q !< prend la valeur de q1 ou q2 selon les cas 
    202   real :: aat !< variable de travail =q*tetar(i,j,k) 
    203   real :: ssec !< variable de travail  
    204   real :: zetat !< variable de travail : profondeur ou t = trans(iiglen) 
    205   real,dimension(nx,ny,nz) :: si1 !< tableau de calcul 
    206   real,dimension(nx,ny,nz) :: si2 !< tableau de calcul 
    207   real,dimension(nx,ny) :: tab_mx 
    208   real,dimension(nx,ny) :: tab_my 
    209   real,dimension(nx,ny) :: tab2d 
    210  
    211 !  real,dimension(nz) :: e          ! vertical coordinate in ice, scaled to h zeta 
    212   real :: de= 1./(nz-1) 
    213 ! variables openmp :   
    214   !$  integer :: rang 
    215   !$ integer, dimension(:), allocatable :: tab_sync 
    216   !$ integer :: nb_taches 
    217  
    218  
    219   e(1)=0. 
    220   e(nz)=1. 
    221    
    222   !$OMP PARALLEL PRIVATE(bat,q,aat,ssec,zetat) 
    223   !$ nb_taches = OMP_GET_NUM_THREADS() 
    224   !$OMP DO  
    225   do k=1,nz 
    226      if ((k.ne.1).and.(k.ne.nz)) e(k)=(k-1.)/(nz-1.) 
    227   end do 
    228   !$OMP END DO NOWAIT 
    229  
    230   ! exposant de la loi de deformation utilisee : glen(iiglen) 
    231   ! l'exposant correspondant a iiglen est defini dans deformation_mod 
    232   a5 = glen(iiglen) + 2   
    233   a4 = glen(iiglen) + 1 
    234  
    235   !    boucle pour btt a k=1 
    236   ti1=273.15*273.15 
    237   !$OMP DO 
    238   do j=2,ny-1 
    239      do i=2,nx-1 
    240         if (t(i,j,1).le.ttrans(iiglen)) then 
    241            btt(i,j,1,iiglen)=bat1(iiglen)*exp(q1(iiglen)/rgas/ti1*t(i,j,1)) 
    242         else 
    243            btt(i,j,1,iiglen)=bat2(iiglen)*exp(q2(iiglen)/rgas/ti1*t(i,j,1)) 
    244         endif 
    245      end do 
    246   end do 
    247   !$OMP END DO 
    248  
    249   ! boucle pour tous les autres btt 
    250   !$OMP DO COLLAPSE(2) 
    251   do k=nz-1,1,-1 
    252         do j=2,ny-1   
    253         do i=2,nx-1   
    254            if (h(i,j).gt.1.) then 
    255               if ((teta(i,j,k).le.ttrans(iiglen)).and. & 
    256                    (teta(i,j,k+1).le.ttrans(iiglen))) then 
    257                  bat=bat1(iiglen) 
    258                  q=q1(iiglen) 
    259                  aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1) 
    260                  aat=max(-1.8,aat) 
    261                  aat=min(80.,aat) 
    262                  ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k))) 
    263                  si1(i,j,k)=ssec/(aat+a4)*(e(k)**(aat+a4)-e(k+1)**(aat+a4)) 
    264                  si2(i,j,k)=((e(k)**(aat+a5))/(aat+a5) & 
    265                       -(e(k+1)**(aat+a5))/(aat+a5) & 
    266                       -(e(k+1)**(aat+a4))*e(k)+e(k+1)**(aat+a5)) & 
    267                       * ssec/(aat+a4)  
    268                  btt(i,j,k+1,iiglen)=bat*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k))) 
    269               else if ((teta(i,j,k).ge.ttrans(iiglen)).and. & 
    270                    (teta(i,j,k+1).ge.ttrans(iiglen))) then 
    271                  bat=bat2(iiglen) 
    272                  q=q2(iiglen) 
    273                  aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1) 
    274                  aat=max(-1.8,aat) 
    275                  aat=min(80.,aat) 
    276                  ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k))) 
    277                  si1(i,j,k)=ssec/(aat+a4)*(e(k)**(aat+a4)-e(k+1)**(aat+a4)) 
    278                  si2(i,j,k)=((e(k)**(aat+a5))/(aat+a5) & 
    279                       -(e(k+1)**(aat+a5))/(aat+a5) & 
    280                       -(e(k+1)**(aat+a4))*e(k)+e(k+1)**(aat+a5)) & 
    281                       * ssec/(aat+a4)  
    282                  btt(i,j,k+1,iiglen)=bat*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k))) 
    283  
    284                  ! cas ou la temperature de la maille est en partie au dessus et au dessous 
    285                  ! de ttrans(iiglen) 
    286               else if ((teta(i,j,k).lt.ttrans(iiglen)).and. & 
    287                    (teta(i,j,k+1).gt.ttrans(iiglen))) then 
    288                  !         calcul de la profondeur de transition 
    289                  if (abs(teta(i,j,k)-teta(i,j,k+1)).lt.1.e-5) then 
    290                     zetat=(e(k)+e(k+1))/2. 
    291                  else 
    292                     zetat=e(k+1)-(ttrans(iiglen)-teta(i,j,k+1))/ & 
    293                          (teta(i,j,k)-teta(i,j,k+1))*de 
    294                  endif 
    295  
    296                  !         integration entre zeta2 et zetat, loi bat2  
    297                  bat=bat2(iiglen) 
    298                  q=q2(iiglen) 
    299                  aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1) 
    300                  aat=max(-1.8,aat) 
    301                  aat=min(80.,aat) 
    302                  ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k))) 
    303                  si1(i,j,k)=ssec/(aat+a4)*(zetat**(aat+a4)-e(k+1)**(aat+a4)) 
    304                  si2(i,j,k)=((zetat**(aat+a5))/(aat+a5) & 
    305                       -(e(k+1)**(aat+a5))/(aat+a5) & 
    306                       -(e(k+1)**(aat+a4))*zetat+e(k+1)**(aat+a5)) & 
    307                       * ssec/(aat+a4) 
    308                  btt(i,j,k+1,iiglen)=bat*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k))) 
    309  
    310                  !         integration entre zetat et zeta1, loi bat1  
    311                  bat=bat1(iiglen) 
    312                  q=q1(iiglen) 
    313                  aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1) 
    314                  aat=max(-1.8,aat) 
    315                  aat=min(80.,aat) 
    316                  ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k))) 
    317                  si1(i,j,k)=si1(i,j,k)+ & 
    318                       ssec/(aat+a4)*(e(k)**(aat+a4)-zetat**(aat+a4)) 
    319                  si2(i,j,k)=si2(i,j,k) & 
    320                       +((e(k)**(aat+a5))/(aat+a5) & 
    321                       -(zetat**(aat+a5))/(aat+a5) & 
    322                       -(zetat**(aat+a4))*e(k)+zetat**(aat+a5)) & 
    323                       * ssec/(aat+a4) 
    324  
    325  
    326                  ! deuxieme cas ou la temperature de la maille est en partie au dessus et  
    327                  ! au dessous de ttrans(iiglen) 
    328               else if ((teta(i,j,k).gt.ttrans(iiglen)).and. & 
    329                    (teta(i,j,k+1).lt.ttrans(iiglen))) then 
    330                  !         calcul de la profondeur de transition 
    331                  if (abs(teta(i,j,k)-teta(i,j,k+1)).lt.1.e-5) then 
    332                     zetat=(e(k)+e(k+1))/2. 
    333                  else 
    334                     zetat=e(k+1)-(ttrans(iiglen)-teta(i,j,k+1))/ & 
    335                          (teta(i,j,k)-teta(i,j,k+1))*de 
    336                  endif 
    337  
    338                  !         integration entre zeta2 et zetat, loi bat1  
    339                  bat=bat1(iiglen) 
    340                  q=q1(iiglen) 
    341                  aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1) 
    342                  aat=max(-1.8,aat) 
    343                  aat=min(80.,aat) 
    344                  ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k))) 
    345                  si1(i,j,k)=ssec/(aat+a4)*(zetat**(aat+a4)-e(k+1)**(aat+a4)) 
    346                  si2(i,j,k)=((zetat**(aat+a5))/(aat+a5) & 
    347                       -(e(k+1)**(aat+a5))/(aat+a5) & 
    348                       -(e(k+1)**(aat+a4))*zetat+e(k+1)**(aat+a5)) & 
    349                       * ssec/(aat+a4) 
    350                  btt(i,j,k+1,iiglen)=bat*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k))) 
    351  
    352                  !         integration entre zetat et zeta1, loi bat2  
    353                  bat=bat2(iiglen) 
    354                  q=q2(iiglen) 
    355                  aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1) 
    356                  aat=max(-1.8,aat) 
    357                  aat=min(80.,aat) 
    358                  ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k))) 
    359                  si1(i,j,k)=si1(i,j,k)+ & 
    360                       ssec/(aat+a4)*(e(k)**(aat+a4)-zetat**(aat+a4)) 
    361                  si2(i,j,k)=si2(i,j,k) & 
    362                       +((e(k)**(aat+a5))/(aat+a5) & 
    363                       -(zetat**(aat+a5))/(aat+a5) & 
    364                       -(zetat**(aat+a4))*e(k)+zetat**(aat+a5)) & 
    365                       * ssec/(aat+a4) 
    366               endif 
    367            endif 
    368         end do 
    369      end do 
    370   end do 
    371   !$OMP END DO NOWAIT 
    372    
    373   !     integration de sa et s2a 
    374   !$OMP DO 
    375   do j=1,ny 
    376      do i=1,nx   
    377         sa(i,j,nz,iiglen)=0. 
    378         s2a(i,j,nz,iiglen)=0. 
    379         if (h(i,j).le.1.) btt(i,j,nz,iiglen)=bat2(iiglen) 
    380      end do 
    381   end do 
    382   !$OMP END DO 
    383   !$OMP END PARALLEL 
    384    
    385  
    386   ! Allocation et initialisation du tableau tab_sync 
    387   ! qui gere la synchronisation entre les differents threads 
    388   !$ allocate(tab_sync(0:nb_taches-1)) 
    389   !$ tab_sync(0:nb_taches-1) = 1 
    390   !$OMP PARALLEL private(i,j,k,rang) shared(tab_sync) 
    391   !$ rang = OMP_GET_THREAD_NUM()   
    392   do k=nz-1,1,-1 
     79  !------------------------------------------------------------------------------------------- 
     80 
     81  !> SUBROUTINE: init_deformation 
     82  !! Routine qui lit les variables de deformation  
     83  !> 
     84  subroutine init_deformation  
     85 
     86    real :: exposant_1 
     87    real :: temp_trans_1 
     88    real :: enhanc_fact_1 
     89    real :: coef_cold_1 
     90    real :: Q_cold_1 
     91    real :: coef_warm_1 
     92    real :: Q_warm_1 
     93 
     94    real :: exposant_2 
     95    real :: temp_trans_2 
     96    real :: enhanc_fact_2 
     97    real :: coef_cold_2 
     98    real :: Q_cold_2 
     99    real :: coef_warm_2 
     100    real :: Q_warm_2 
     101 
     102    namelist/loidef_1/exposant_1,temp_trans_1,enhanc_fact_1,   &  
     103         coef_cold_1,Q_cold_1,coef_warm_1,Q_warm_1 
     104 
     105    namelist/loidef_2/exposant_2,temp_trans_2,enhanc_fact_2,   &  
     106         coef_cold_2,Q_cold_2,coef_warm_2,Q_warm_2 
     107 
     108 
     109    ! loi 1 
     110    !-------------------------------------------------------------------------------- 
     111 
     112    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
     113    read(num_param,loidef_1) 
     114    write(num_rep_42,loidef_1) 
     115 
     116    glen(1)   = exposant_1  
     117    ttrans(1) = temp_trans_1 
     118    sf(1)     = enhanc_fact_1   
     119    Bat1(1)   = coef_cold_1 
     120    Q1(1)     = Q_cold_1  
     121    Bat2(1)   = coef_warm_1 
     122    Q2(1)     = Q_warm_1  
     123 
     124 
     125    ! loi 2 
     126    !-------------------------------------------------------------------------------- 
     127 
     128    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
     129    read(num_param,loidef_2) 
     130    write(num_rep_42,loidef_2) 
     131 
     132    write(num_rep_42,*)'!___________________________________________________________'  
     133    write(num_rep_42,*)'! loi de deformation            module deformation_mod_2lois' 
     134    write(num_rep_42,*)'! exposant (glen), temperature de transition (ttrans)' 
     135    write(num_rep_42,*)'! enhancement factor (sf)' 
     136    write(num_rep_42,*)'! pour les temperatures inf. a Temp_trans :' 
     137    write(num_rep_42,*)'!                            coef_cold (Bat1) et Q_cold (Q1)' 
     138    write(num_rep_42,*)'! pour les temperatures sup. a Temp_trans :' 
     139    write(num_rep_42,*)'!                            coef_warm (Bat2) et Q_warm (Q2)' 
     140    write(num_rep_42,*)'!___________________________________________________________'  
     141 
     142    glen(2)   = exposant_2  
     143    ttrans(2) = temp_trans_2 
     144    sf(2)     = enhanc_fact_2   
     145    Bat1(2)   = coef_cold_2 
     146    Q1(2)     = Q_cold_2  
     147    Bat2(2)   = coef_warm_2 
     148    Q2(2)     = Q_warm_2  
     149 
     150    ! autre parametres ne changeant pas d'un run a l'autre 
     151    RGAS=8.314 
     152 
     153    ! application des sf 
     154 
     155    do iglen= n1poly,n2poly 
     156       bat1(iglen)=bat1(iglen)*sf(iglen) 
     157       bat2(iglen)=bat2(iglen)*sf(iglen) 
     158    end do 
     159 
     160    ddx(:,:,:)=0. 
     161    ddy(:,:,:)=0. 
     162 
     163 
     164  end subroutine init_deformation 
     165 
     166 
     167  !-------------------------------------------------------------------- 
     168  subroutine flow_general 
     169 
     170    !$ USE OMP_LIB 
     171    implicit none 
     172 
     173    integer :: i,j,k 
     174 
     175    !$OMP PARALLEL 
     176    !$OMP WORKSHARE 
     177    teta(:,:,:) = t(:,:,1:nz)-tpmp(:,:,:) 
     178    !$OMP END WORKSHARE 
     179 
     180    !$OMP DO COLLAPSE(2) 
     181    do k=nz-1,1,-1 
     182       do i=2,nx-1 
     183          do j=2,ny-1 
     184             ti2(i,j,k) = (273.15+(tpmp(i,j,k+1)+tpmp(i,j,k))/2.)* & 
     185                  (273.15+(t(i,j,k+1)+t(i,j,k))/2.) 
     186          end do 
     187       end do 
     188    end do 
     189    !$OMP END DO 
     190    !$OMP END PARALLEL 
     191 
     192  end subroutine flow_general 
     193  !--------------------------------------------------------------------------------------- 
     194  subroutine flowlaw (iiglen) 
     195 
     196    !$ USE OMP_LIB 
     197 
     198    implicit none 
     199 
     200    integer :: i,j,k 
     201    integer ::  iiglen   !< compteur pour la boucle flowlaw 
     202    real :: a5 !< exposant de la loi de def +2  
     203    real :: a4 !< exposant de la loi de def +1  
     204    real :: ti1 !< utilise pour le calcul de btt a k=1 
     205    real :: bat !< prend la valeur de bat1 ou bat2 selon les cas 
     206    real :: q !< prend la valeur de q1 ou q2 selon les cas 
     207    real :: aat !< variable de travail =q*tetar(i,j,k) 
     208    real :: ssec !< variable de travail  
     209    real :: zetat !< variable de travail : profondeur ou t = trans(iiglen) 
     210    real,dimension(nx,ny,nz) :: si1 !< tableau de calcul 
     211    real,dimension(nx,ny,nz) :: si2 !< tableau de calcul 
     212    real,dimension(nx,ny) :: tab_mx 
     213    real,dimension(nx,ny) :: tab_my 
     214    real,dimension(nx,ny) :: tab2d 
     215 
     216    !  real,dimension(nz) :: e          ! vertical coordinate in ice, scaled to h zeta 
     217    real :: de= 1./(nz-1) 
     218    ! variables openmp :   
     219    !$  integer :: rang 
     220    !$ integer, dimension(:), allocatable :: tab_sync 
     221    !$ integer :: nb_taches 
     222 
     223 
     224    e(1)=0. 
     225    e(nz)=1. 
     226 
     227    !$OMP PARALLEL PRIVATE(bat,q,aat,ssec,zetat) 
     228    !$ nb_taches = OMP_GET_NUM_THREADS() 
     229    !$OMP DO  
     230    do k=1,nz 
     231       if ((k.ne.1).and.(k.ne.nz)) e(k)=(k-1.)/(nz-1.) 
     232    end do 
     233    !$OMP END DO NOWAIT 
     234 
     235    ! exposant de la loi de deformation utilisee : glen(iiglen) 
     236    ! l'exposant correspondant a iiglen est defini dans deformation_mod 
     237    a5 = glen(iiglen) + 2   
     238    a4 = glen(iiglen) + 1 
     239 
     240    !    boucle pour btt a k=1 
     241    ti1=273.15*273.15 
     242    !$OMP DO 
     243    do j=2,ny-1 
     244       do i=2,nx-1 
     245          if (t(i,j,1).le.ttrans(iiglen)) then 
     246             btt(i,j,1,iiglen)=bat1(iiglen)*exp(q1(iiglen)/rgas/ti1*t(i,j,1)) 
     247          else 
     248             btt(i,j,1,iiglen)=bat2(iiglen)*exp(q2(iiglen)/rgas/ti1*t(i,j,1)) 
     249          endif 
     250       end do 
     251    end do 
     252    !$OMP END DO 
     253 
     254    ! boucle pour tous les autres btt 
     255    !$OMP DO COLLAPSE(2) 
     256    do k=nz-1,1,-1 
     257       do j=2,ny-1   
     258          do i=2,nx-1   
     259             if (h(i,j).gt.1.) then 
     260                if ((teta(i,j,k).le.ttrans(iiglen)).and. & 
     261                     (teta(i,j,k+1).le.ttrans(iiglen))) then 
     262                   bat=bat1(iiglen) 
     263                   q=q1(iiglen) 
     264                   aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1) 
     265                   aat=max(-1.8,aat) 
     266                   aat=min(80.,aat) 
     267                   ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k))) 
     268                   si1(i,j,k)=ssec/(aat+a4)*(e(k)**(aat+a4)-e(k+1)**(aat+a4)) 
     269                   si2(i,j,k)=((e(k)**(aat+a5))/(aat+a5) & 
     270                        -(e(k+1)**(aat+a5))/(aat+a5) & 
     271                        -(e(k+1)**(aat+a4))*e(k)+e(k+1)**(aat+a5)) & 
     272                        * ssec/(aat+a4)  
     273                   btt(i,j,k+1,iiglen)=bat*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k))) 
     274                else if ((teta(i,j,k).ge.ttrans(iiglen)).and. & 
     275                     (teta(i,j,k+1).ge.ttrans(iiglen))) then 
     276                   bat=bat2(iiglen) 
     277                   q=q2(iiglen) 
     278                   aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1) 
     279                   aat=max(-1.8,aat) 
     280                   aat=min(80.,aat) 
     281                   ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k))) 
     282                   si1(i,j,k)=ssec/(aat+a4)*(e(k)**(aat+a4)-e(k+1)**(aat+a4)) 
     283                   si2(i,j,k)=((e(k)**(aat+a5))/(aat+a5) & 
     284                        -(e(k+1)**(aat+a5))/(aat+a5) & 
     285                        -(e(k+1)**(aat+a4))*e(k)+e(k+1)**(aat+a5)) & 
     286                        * ssec/(aat+a4)  
     287                   btt(i,j,k+1,iiglen)=bat*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k))) 
     288 
     289                   ! cas ou la temperature de la maille est en partie au dessus et au dessous 
     290                   ! de ttrans(iiglen) 
     291                else if ((teta(i,j,k).lt.ttrans(iiglen)).and. & 
     292                     (teta(i,j,k+1).gt.ttrans(iiglen))) then 
     293                   !         calcul de la profondeur de transition 
     294                   if (abs(teta(i,j,k)-teta(i,j,k+1)).lt.1.e-5) then 
     295                      zetat=(e(k)+e(k+1))/2. 
     296                   else 
     297                      zetat=e(k+1)-(ttrans(iiglen)-teta(i,j,k+1))/ & 
     298                           (teta(i,j,k)-teta(i,j,k+1))*de 
     299                   endif 
     300 
     301                   !         integration entre zeta2 et zetat, loi bat2  
     302                   bat=bat2(iiglen) 
     303                   q=q2(iiglen) 
     304                   aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1) 
     305                   aat=max(-1.8,aat) 
     306                   aat=min(80.,aat) 
     307                   ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k))) 
     308                   si1(i,j,k)=ssec/(aat+a4)*(zetat**(aat+a4)-e(k+1)**(aat+a4)) 
     309                   si2(i,j,k)=((zetat**(aat+a5))/(aat+a5) & 
     310                        -(e(k+1)**(aat+a5))/(aat+a5) & 
     311                        -(e(k+1)**(aat+a4))*zetat+e(k+1)**(aat+a5)) & 
     312                        * ssec/(aat+a4) 
     313                   btt(i,j,k+1,iiglen)=bat*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k))) 
     314 
     315                   !         integration entre zetat et zeta1, loi bat1  
     316                   bat=bat1(iiglen) 
     317                   q=q1(iiglen) 
     318                   aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1) 
     319                   aat=max(-1.8,aat) 
     320                   aat=min(80.,aat) 
     321                   ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k))) 
     322                   si1(i,j,k)=si1(i,j,k)+ & 
     323                        ssec/(aat+a4)*(e(k)**(aat+a4)-zetat**(aat+a4)) 
     324                   si2(i,j,k)=si2(i,j,k) & 
     325                        +((e(k)**(aat+a5))/(aat+a5) & 
     326                        -(zetat**(aat+a5))/(aat+a5) & 
     327                        -(zetat**(aat+a4))*e(k)+zetat**(aat+a5)) & 
     328                        * ssec/(aat+a4) 
     329 
     330 
     331                   ! deuxieme cas ou la temperature de la maille est en partie au dessus et  
     332                   ! au dessous de ttrans(iiglen) 
     333                else if ((teta(i,j,k).gt.ttrans(iiglen)).and. & 
     334                     (teta(i,j,k+1).lt.ttrans(iiglen))) then 
     335                   !         calcul de la profondeur de transition 
     336                   if (abs(teta(i,j,k)-teta(i,j,k+1)).lt.1.e-5) then 
     337                      zetat=(e(k)+e(k+1))/2. 
     338                   else 
     339                      zetat=e(k+1)-(ttrans(iiglen)-teta(i,j,k+1))/ & 
     340                           (teta(i,j,k)-teta(i,j,k+1))*de 
     341                   endif 
     342 
     343                   !         integration entre zeta2 et zetat, loi bat1  
     344                   bat=bat1(iiglen) 
     345                   q=q1(iiglen) 
     346                   aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1) 
     347                   aat=max(-1.8,aat) 
     348                   aat=min(80.,aat) 
     349                   ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k))) 
     350                   si1(i,j,k)=ssec/(aat+a4)*(zetat**(aat+a4)-e(k+1)**(aat+a4)) 
     351                   si2(i,j,k)=((zetat**(aat+a5))/(aat+a5) & 
     352                        -(e(k+1)**(aat+a5))/(aat+a5) & 
     353                        -(e(k+1)**(aat+a4))*zetat+e(k+1)**(aat+a5)) & 
     354                        * ssec/(aat+a4) 
     355                   btt(i,j,k+1,iiglen)=bat*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k))) 
     356 
     357                   !         integration entre zetat et zeta1, loi bat2  
     358                   bat=bat2(iiglen) 
     359                   q=q2(iiglen) 
     360                   aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1) 
     361                   aat=max(-1.8,aat) 
     362                   aat=min(80.,aat) 
     363                   ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k))) 
     364                   si1(i,j,k)=si1(i,j,k)+ & 
     365                        ssec/(aat+a4)*(e(k)**(aat+a4)-zetat**(aat+a4)) 
     366                   si2(i,j,k)=si2(i,j,k) & 
     367                        +((e(k)**(aat+a5))/(aat+a5) & 
     368                        -(zetat**(aat+a5))/(aat+a5) & 
     369                        -(zetat**(aat+a4))*e(k)+zetat**(aat+a5)) & 
     370                        * ssec/(aat+a4) 
     371                endif 
     372             endif 
     373          end do 
     374       end do 
     375    end do 
     376    !$OMP END DO NOWAIT 
     377 
     378    !     integration de sa et s2a 
     379    !$OMP DO 
     380    do j=1,ny 
     381       do i=1,nx   
     382          sa(i,j,nz,iiglen)=0. 
     383          s2a(i,j,nz,iiglen)=0. 
     384          if (h(i,j).le.1.) btt(i,j,nz,iiglen)=bat2(iiglen) 
     385       end do 
     386    end do 
     387    !$OMP END DO 
     388    !$OMP END PARALLEL 
     389 
     390 
     391    ! Allocation et initialisation du tableau tab_sync 
     392    ! qui gere la synchronisation entre les differents threads 
     393    !$ allocate(tab_sync(0:nb_taches-1)) 
     394    !$ tab_sync(0:nb_taches-1) = 1 
     395    !$OMP PARALLEL private(i,j,k,rang) shared(tab_sync) 
     396    !$ rang = OMP_GET_THREAD_NUM()   
     397    do k=nz-1,1,-1 
    393398       ! Synchronisation des threads 
    394      !$     if (rang /= 0) then 
    395      !$        do 
    396      !$OMP FLUSH(tab_sync)  
    397      !$           if (tab_sync(rang-1)>=tab_sync(rang)+1) exit 
    398      !$        enddo 
    399      !$OMP FLUSH(sa) 
    400      !$OMP FLUSH(s2a) 
    401      !$     endif 
    402      !$OMP DO SCHEDULE(STATIC) 
    403      do j=1,ny 
    404         do i=1,nx 
    405            if (h(i,j).gt.1.) then 
    406               sa(i,j,k,iiglen)=sa(i,j,k+1,iiglen)-si1(i,j,k) 
    407               s2a(i,j,k,iiglen)=s2a(i,j,k+1,iiglen)+si2(i,j,k)+ & 
    408                    sa(i,j,k+1,iiglen)*de 
    409            else 
    410               sa(i,j,k,iiglen)=bat2(iiglen)/a4*(1-e(k)**a4) 
    411               s2a(i,j,k,iiglen)=bat2(iiglen)/a4*(a4/a5-e(k)+e(k)**a5/a5) 
    412               btt(i,j,k,iiglen)=bat2(iiglen) 
    413            endif 
    414         end do 
    415      end do 
    416      !$OMP END DO NOWAIT 
    417 !     ! Mis a jour du tableau tab_sync 
    418      !$     tab_sync(rang) = tab_sync(rang) + 1 
    419      !$OMP FLUSH(tab_sync,sa,s2a) 
    420   end do 
    421  
    422   !$OMP WORKSHARE 
    423   !     cas particulier des bords 
    424   sa(:,1,:,iiglen)=sa(:,2,:,iiglen) 
    425   s2a(:,1,:,iiglen)=s2a(:,2,:,iiglen) 
    426   btt(:,1,:,iiglen)=btt(:,2,:,iiglen) 
    427   sa(:,ny,:,iiglen)=sa(:,ny-1,:,iiglen) 
    428   s2a(:,ny,:,iiglen)=s2a(:,ny-1,:,iiglen) 
    429   btt(:,ny,:,iiglen)=btt(:,ny-1,:,iiglen) 
    430  
    431   sa(1,:,:,iiglen)=sa(2,:,:,iiglen) 
    432   s2a(1,:,:,iiglen)=s2a(2,:,:,iiglen) 
    433   btt(1,:,:,iiglen)=btt(2,:,:,iiglen) 
    434   sa(nx,:,:,iiglen)=sa(nx-1,:,:,iiglen) 
    435   s2a(nx,:,:,iiglen)=s2a(nx-1,:,:,iiglen) 
    436   btt(nx,:,:,iiglen)=btt(nx-1,:,:,iiglen) 
    437   !$OMP END WORKSHARE 
    438   !$OMP END PARALLEL 
    439    
    440   ! calcul des moyennes sur les mailles staggered 
    441   do k=1,nz 
    442      tab2d=sa(:,:,k,iiglen) 
    443  
    444      call moy_mxmy(nx,ny,tab2d,tab_mx,tab_my) 
    445      sa_mx(:,:,k,iglen)=tab_mx 
    446      sa_my(:,:,k,iglen)=tab_my 
    447  
    448      tab2d=s2a(:,:,k,iiglen) 
    449      call moy_mxmy(nx,ny,tab2d,tab_mx,tab_my) 
    450      s2a_mx(:,:,k,iglen)=tab_mx 
    451      s2a_my(:,:,k,iglen)=tab_my 
    452   end do 
    453   
    454 ! attribution de debug_3d pour les sorties s2a 
    455 ! loi 1 -> 190 dans description variable -> 90 dans debug_3d 
    456 !debug_3d(:,:,90+iiglen-1) = s2a(:,:,1,iiglen) 
    457  
    458 end subroutine flowlaw 
     399       !$     if (rang /= 0) then 
     400       !$        do 
     401       !$OMP FLUSH(tab_sync)  
     402       !$           if (tab_sync(rang-1)>=tab_sync(rang)+1) exit 
     403       !$        enddo 
     404       !$OMP FLUSH(sa) 
     405       !$OMP FLUSH(s2a) 
     406       !$     endif 
     407       !$OMP DO SCHEDULE(STATIC) 
     408       do j=1,ny 
     409          do i=1,nx 
     410             if (h(i,j).gt.1.) then 
     411                sa(i,j,k,iiglen)=sa(i,j,k+1,iiglen)-si1(i,j,k) 
     412                s2a(i,j,k,iiglen)=s2a(i,j,k+1,iiglen)+si2(i,j,k)+ & 
     413                     sa(i,j,k+1,iiglen)*de 
     414             else 
     415                sa(i,j,k,iiglen)=bat2(iiglen)/a4*(1-e(k)**a4) 
     416                s2a(i,j,k,iiglen)=bat2(iiglen)/a4*(a4/a5-e(k)+e(k)**a5/a5) 
     417                btt(i,j,k,iiglen)=bat2(iiglen) 
     418             endif 
     419          end do 
     420       end do 
     421       !$OMP END DO NOWAIT 
     422       !     ! Mis a jour du tableau tab_sync 
     423       !$     tab_sync(rang) = tab_sync(rang) + 1 
     424       !$OMP FLUSH(tab_sync,sa,s2a) 
     425    end do 
     426 
     427    !$OMP WORKSHARE 
     428    !     cas particulier des bords 
     429    sa(:,1,:,iiglen)=sa(:,2,:,iiglen) 
     430    s2a(:,1,:,iiglen)=s2a(:,2,:,iiglen) 
     431    btt(:,1,:,iiglen)=btt(:,2,:,iiglen) 
     432    sa(:,ny,:,iiglen)=sa(:,ny-1,:,iiglen) 
     433    s2a(:,ny,:,iiglen)=s2a(:,ny-1,:,iiglen) 
     434    btt(:,ny,:,iiglen)=btt(:,ny-1,:,iiglen) 
     435 
     436    sa(1,:,:,iiglen)=sa(2,:,:,iiglen) 
     437    s2a(1,:,:,iiglen)=s2a(2,:,:,iiglen) 
     438    btt(1,:,:,iiglen)=btt(2,:,:,iiglen) 
     439    sa(nx,:,:,iiglen)=sa(nx-1,:,:,iiglen) 
     440    s2a(nx,:,:,iiglen)=s2a(nx-1,:,:,iiglen) 
     441    btt(nx,:,:,iiglen)=btt(nx-1,:,:,iiglen) 
     442    !$OMP END WORKSHARE 
     443    !$OMP END PARALLEL 
     444 
     445    ! calcul des moyennes sur les mailles staggered 
     446    do k=1,nz 
     447       tab2d=sa(:,:,k,iiglen) 
     448 
     449       call moy_mxmy(nx,ny,tab2d,tab_mx,tab_my) 
     450       sa_mx(:,:,k,iglen)=tab_mx 
     451       sa_my(:,:,k,iglen)=tab_my 
     452 
     453       tab2d=s2a(:,:,k,iiglen) 
     454       call moy_mxmy(nx,ny,tab2d,tab_mx,tab_my) 
     455       s2a_mx(:,:,k,iglen)=tab_mx 
     456       s2a_my(:,:,k,iglen)=tab_my 
     457    end do 
     458 
     459    ! attribution de debug_3d pour les sorties s2a 
     460    ! loi 1 -> 190 dans description variable -> 90 dans debug_3d 
     461    !debug_3d(:,:,90+iiglen-1) = s2a(:,:,1,iiglen) 
     462 
     463  end subroutine flowlaw 
    459464 
    460465 
Note: See TracChangeset for help on using the changeset viewer.