Changeset 397


Ignore:
Timestamp:
03/24/23 16:05:58 (12 months ago)
Author:
dumas
Message:

use only in module furst_schoof_mod

Location:
branches/GRISLIv3/SOURCES
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/GRISLIv3/SOURCES/deformation_mod_2lois.f90

    r390 r397  
    2727module deformation_mod_2lois 
    2828 
    29   use module3d_phy, only:nx,ny,nz,e,num_param,num_rep_42,rgas,T,iglen,tpmp,H 
     29  use module3d_phy, only:nx,ny,nz!,e,num_param,num_rep_42,rgas,T,iglen,tpmp,H 
    3030 
    3131  implicit none 
     
    8282  !! Routine qui lit les variables de deformation  
    8383  !> 
    84   subroutine init_deformation  
    85  
     84  subroutine init_deformation 
     85     
     86    use module3d_phy, only:num_param,num_rep_42,rgas,iglen 
     87 
     88    implicit none 
     89     
    8690    real :: exposant_1 
    8791    real :: temp_trans_1 
     
    167171  !-------------------------------------------------------------------- 
    168172  subroutine flow_general 
    169  
     173     
     174    use module3d_phy, only:T,tpmp 
    170175    !$ USE OMP_LIB 
     176     
    171177    implicit none 
    172178 
     
    193199  !--------------------------------------------------------------------------------------- 
    194200  subroutine flowlaw (iiglen) 
    195  
     201     
     202    use module3d_phy, only:e,T,rgas,H,iglen 
    196203    !$ USE OMP_LIB 
    197204 
  • branches/GRISLIv3/SOURCES/furst_schoof_mod.f90

    r271 r397  
    33module furst_schoof_mod 
    44 
    5 use module3d_phy 
    6 use deformation_mod_2lois 
    7  
    8 implicit none 
    9  
    10 real :: frot_coef ! afq: fixed in Tsai but = beta in Schoof (for m_weert = 1) 
    11 real :: alpha_flot 
    12 integer :: gr_select 
    13 real,dimension(nx,ny) :: back_force_x 
    14 real,dimension(nx,ny) :: back_force_y 
    15 real, parameter :: m_weert = 1.0 ! afq: is this defined elsewhere in the model??? 
    16 real, parameter :: inv_mweert = 1./m_weert 
     5  use module3d_phy, only:nx,ny 
     6  use deformation_mod_2lois, only: 
     7 
     8  implicit none 
     9 
     10  real :: frot_coef ! afq: fixed in Tsai but = beta in Schoof (for m_weert = 1) 
     11  real :: alpha_flot 
     12  integer :: gr_select 
     13  real,dimension(nx,ny) :: back_force_x 
     14  real,dimension(nx,ny) :: back_force_y 
     15  real, parameter :: m_weert = 1.0 ! afq: is this defined elsewhere in the model??? 
     16  real, parameter :: inv_mweert = 1./m_weert 
    1717 
    1818contains 
    1919 
    20 subroutine init_furst_schoof 
    21  
    22 namelist/furst_schoof/frot_coef,gr_select 
    23  
    24 rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
    25 read(num_param,furst_schoof) 
    26 write(num_rep_42,'(A)')'!___________________________________________________________'  
    27 write(num_rep_42,'(A)') '&furst_schoof                                ! nom du bloc ' 
    28 write(num_rep_42,*) 
    29 write(num_rep_42,*)'frot_coef  = ',frot_coef 
    30 write(num_rep_42,*)'gr_select  = ',gr_select 
    31 write(num_rep_42,*)'/' 
    32 write(num_rep_42,*)'! frot_coef : solid friction law coef 0.6 in GMD 2018' 
    33 write(num_rep_42,*)'! gr_select = 1 : Tsai , 2 : Schoof' 
    34  
    35 !  frot = 0.035 ! f in the solid friction law, 0.035 in Ritz et al. 2001         
    36  
    37   ! choix Tsai (loi de Coulomb) ou Schoof (loi de Weertman) 
    38   alpha_flot= ro/row 
    39   back_force_x(:,:)=0.1 
    40   back_force_y(:,:)=0.1 
    41    
    42 end subroutine init_furst_schoof   
    43    
    44 subroutine interpol_glflux 
    45  
    46   implicit none 
    47   real :: xpos,ypos 
    48   real :: Hgl 
    49   !  integer :: igl 
    50  
    51   real :: phi_im2, phi_im1, phi_i, phi_ip1, phi_ip2, phi_ip3 
    52   real :: phi_jm2, phi_jm1, phi_j, phi_jp1, phi_jp2, phi_jp3 
    53   real,dimension(nx,ny) :: phi_xgl, phi_ygl ! flux gr line sur maille stag 
    54   integer,dimension(nx,ny) :: countx, county   ! how often do we modify ux/uy 
    55   real :: phi_prescr 
    56   real :: toutpetit = 1e-6 
    57   real :: denom, prodscal 
    58  
    59   real :: bfx, bfy 
    60  
    61   !debug 
    62   real,dimension(nx,ny) :: xpos_tab=9999. 
    63   real,dimension(nx,ny) :: ypos_tab=9999. 
    64   real,dimension(nx,ny) :: Hglx_tab=9999. 
    65   real,dimension(nx,ny) :: Hgly_tab=9999. 
    66   real,dimension(nx,ny) :: Hmxloc, Hmyloc ! pour eviter division par 0 
    67   real,dimension(nx,ny) :: phi_prescr_tabx, phi_prescr_taby 
    68  
    69   real,dimension(nx,ny) :: uxbartemp, uybartemp 
    70    
    71   real,dimension(nx,ny) :: archimtab 
    72    
    73   real,dimension(nx,ny) :: imx_diag_in 
    74  
    75   phi_prescr_tabx=0. 
    76   phi_prescr_taby=0. 
    77  
    78   Hmxloc(:,:)= max(Hmx(:,:),1.) 
    79   Hmyloc(:,:)=max(Hmy(:,:),1.) 
    80  
    81   uxbartemp(:,:)=0. 
    82   uybartemp(:,:)=0. 
    83   countx(:,:)=0. 
    84   county(:,:)=0. 
    85   gr_line_schoof(:,:)=0 
    86  
    87   archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot -sealevel_2d(:,:) 
    88  
    89   imx_diag_in(:,:) = imx_diag(:,:) 
    90    
    91   ! detection des noeuds grounding line 
    92   do j= 3,ny-3 
    93      do i=3,nx-3 
    94         !archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel 
    95         if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel_2d(i,j))) then    ! grounded and with ice 
    96            !           print*,'schoof gr line ij',i,j 
    97            ! selon x (indice i) 
    98             
    99            if (archimtab(i-1,j).LT.0..and.Uxbar(i,j).LT.0..and..not.flot_marais(i-1,j)) then 
    100            !if (flot(i-1,j).and.Uxbar(i,j).LT.0.) then    ! grounding line between i-1,i    CAS WEST (ecoulement vers grille West) 
    101  
    102               call  calc_xgl4schoof(-dx,alpha_flot,real(H(i,j)),real(Bsoc(i,j)),sealevel_2d(i,j),real(H(i-1,j)),real(Bsoc(i-1,j)),sealevel_2d(i-1,j),xpos,Hgl) 
    103               xpos_tab(i,j)=xpos 
    104               Hglx_tab(i,j)=Hgl 
    105                
    106               ! afq: the back force is on the staggered grid, the GL can be either West or East to this point. 
    107               if (xpos .lt. -dx/2.) then 
    108                  bfx = back_force_x(i,j) 
    109               else 
    110                  bfx = back_force_x(i+1,j) 
    111               endif 
    112               if (gr_select.eq.1) then ! flux de Tsai 
    113                  call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfx,phi_prescr) 
    114               else if (gr_select.eq.2) then ! flux de Schoof 
    115                  ! afq: the dragging coef. is on the staggered grid, the GL can be either West or East to this point. 
    116                  if (xpos .lt. -dx/2.) then 
    117                     frot_coef = betamx(i,j) 
    118                  else 
    119                     frot_coef = betamx(i+1,j) 
    120                  endif 
    121                  call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfx,phi_prescr) 
    122               else 
    123                  print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE' 
    124                  stop 
    125               endif 
    126               phi_prescr=-phi_prescr ! Doit être negatif dans le cas West 
    127               phi_prescr_tabx(i,j)=phi_prescr 
    128  
    129               if (xpos .lt. -dx/2.) then  !  GL a west du i staggered,                          o centre, x stag 
    130  
    131                  !  grille    ! .....x......o......x......o...G..x......O......x......o......x......o 
    132                  !            !            i-2           i-1            i            i+1           i+2 
    133                  !flux        !     in            out         G out            in 
    134  
    135  
    136 !afq --                 if (.not.(ilemx(i,j)).and..not.(ilemx(i-1,j))) then 
    137                      
    138                  ! extrapolation pour avoir uxbar(i-1,j) qui sera ensuite prescrit dans call rempli_L2 
    139  
    140                  phi_im2 = Uxbar(i-2,j) * Hmxloc(i-2,j) 
    141                  phi_xgl(i-1,j)  = phi_im2 + dx * (phi_prescr-phi_im2)/(2.5 * dx + xpos) 
    142                  uxbartemp(i-1,j)   = uxbartemp(i-1,j)+ phi_xgl(i-1,j) / Hmxloc(i-1,j) ! attention division par 0 possible ? 
    143                  countx   (i-1,j)   = countx(i-1,j)+1 
    144                  imx_diag (i-1,j)   = 1 
    145                  gr_line_schoof(i-1,j) = 1 
    146  
    147                  ! extrapolation pour avoir uxbar(i,j) 
    148  
    149                  phi_ip1 = Uxbar(i+1,j) * Hmxloc(i+1,j) 
    150                  phi_xgl(i,j) = phi_ip1 + dx * (phi_prescr - phi_ip1)/(0.5 * dx - xpos) 
    151                  !                 print*,'uxbar(i,j)=',uxbar(i,j) 
    152                  uxbartemp(i,j)   = uxbartemp(i,j) + phi_xgl(i,j) / Hmxloc(i,j) 
    153                  countx   (i,j)   = countx(i,j)+1 
    154                  imx_diag (i,j)   = 1 
    155                  gr_line_schoof(i,j) = 1 
    156                  !                 print*,'schoof uxbar(i,j)=',uxbar(i,j) 
    157                  !                 print*,'Hgl',Hgl 
    158                  !                 print*,'phi_prescr',phi_prescr, phi_prescr/Hgl 
    159                  !                 print* 
    160  
    161 !afq --                 end if 
    162                
    163               else            !  GL a Est du i staggered,                          o centre, x stag 
    164  
    165                  !  grille    ! .....x......o......x......o......x...G..O......x......o......x......o 
    166                  !            i-2           i-1            i            i+1           i+2 
    167                  !flux        !     in             in           out  G         out           in 
    168  
    169 !afq --                 if (.not.(ilemx(i,j)).and..not.(ilemx(i+1,j))) then 
    170  
    171                  ! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2 
    172  
    173                  phi_im1 = Uxbar(i-1,j) * Hmxloc(i-1,j) 
    174                  phi_xgl(i,j)  = phi_im1 + dx * (phi_prescr-phi_im1)/(1.5 * dx + xpos) 
    175                  uxbartemp(i,j)   = uxbartemp(i,j) + phi_xgl(i,j) / Hmxloc(i,j) ! attention division par 0 possible ? 
    176                  countx   (i,j)   = countx(i,j)+1 
    177                  imx_diag (i,j) = 1 
    178                  gr_line_schoof(i,j) = 1 
    179  
    180                  ! extrapolation pour avoir uxbar(i+1,j) 
    181  
    182                  phi_ip2 = Uxbar(i+2,j) * Hmxloc(i+2,j) 
    183                  phi_xgl(i+1,j) = phi_ip2 + dx * (phi_prescr - phi_ip2)/(1.5 * dx  - xpos) 
    184                  uxbartemp(i+1,j)   = uxbartemp(i+1,j) +  phi_xgl(i+1,j) / Hmxloc(i+1,j) 
    185                  countx   (i+1,j)   = countx(i+1,j)+1 
    186                  imx_diag (i+1,j)   = 1 
    187                  gr_line_schoof(i+1,j) = 1 
    188  
    189 !afq --                 end if 
    190  
    191               end if 
    192  
    193            else if (archimtab(i+1,j).LT.0..and.Uxbar(i+1,j).GT.0..and..not.flot_marais(i+1,j)) then 
    194            !else if (flot(i+1,j).and.Uxbar(i+1,j).GT.0.) then    ! grounding line between i,i+1    CAS EST (ecoulement vers grille Est) 
    195  
    196               call  calc_xgl4schoof(dx,alpha_flot,real(H(i,j)),real(Bsoc(i,j)),sealevel_2d(i,j),real(H(i+1,j)),real(Bsoc(i+1,j)),sealevel_2d(i+1,j),xpos,Hgl) 
    197               xpos_tab(i,j)=xpos 
    198               Hglx_tab(i,j)=Hgl 
    199                
    200               ! afq: the back force is on the staggered grid, the GL can be either West or East to this point. 
    201               if (xpos .lt. dx/2.) then 
    202                  bfx = back_force_x(i,j) 
    203               else 
    204                  bfx = back_force_x(i+1,j) 
    205               endif 
    206               if (gr_select.eq.1) then ! flux de Tsai 
    207                  call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfx,phi_prescr) 
    208               else if (gr_select.eq.2) then ! flux de Schoof 
    209                  ! afq: the dragging coef. is on the staggered grid, the GL can be either West or East to this point. 
    210                  if (xpos .lt. dx/2.) then 
    211                     frot_coef = betamx(i,j) 
    212                  else 
    213                     frot_coef = betamx(i+1,j) 
    214                  endif 
    215                  call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfx,phi_prescr) 
    216               else 
    217                  print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE' 
    218                  stop 
    219               endif 
    220               phi_prescr_tabx(i,j)=phi_prescr 
    221               if (xpos .lt. dx/2.) then  !  GL a west du i staggered,                          o centre, x stag 
    222  
    223                  !  grille    ! .....x......o......x......o......x......O..G...x......o......x......o 
    224                  !            !            i-2           i-1            i            i+1           i+2 
    225                  !flux        !                    in           out        G  out            in 
    226  
    227 !afq --                 if (.not.(ilemx(i,j)).and..not.(ilemx(i+1,j))) then 
    228                      
    229                  ! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2 
    230  
    231                  phi_im1 = Uxbar(i-1,j) * Hmxloc(i-1,j) 
    232                  phi_xgl(i,j)  = phi_im1 + dx * (phi_prescr-phi_im1)/(1.5 * dx + xpos) 
    233                  uxbartemp(i,j)   =  uxbartemp(i,j) + phi_xgl(i,j) / Hmxloc(i,j) ! attention division par 0 possible ? 
    234                  countx   (i,j)   = countx(i,j)+1 
    235                  imx_diag (i,j) = 1 
    236                  gr_line_schoof(i,j) = 1 
    237  
    238                  ! extrapolation pour avoir uxbar(i+1,j) 
    239  
    240                  phi_ip2 = Uxbar(i+2,j) * Hmxloc(i+2,j) 
    241                  phi_xgl(i+1,j) = phi_ip2 + dx * (phi_prescr - phi_ip2)/(1.5* dx - xpos) 
    242                  uxbartemp(i+1,j)   = uxbartemp(i+1,j) + phi_xgl(i+1,j) / Hmxloc(i+1,j) 
    243                  countx   (i+1,j)   = countx(i+1,j)+1 
    244                  imx_diag (i+1,j) = 1 
    245                  gr_line_schoof(i+1,j) = 1 
    246  
    247 !afq --                 end if 
    248                
    249               else   !  GL a west du i staggered,                          o centre, x stag 
    250  
    251                  !  grille    ! ......x......o......x......O......x...G..o......x......o......x.....o 
    252                  !            !             i-1            i            i+1           i+2          i+3 
    253                  !flux        !                     in           out  G        out            in 
    254  
    255 !afq --                 if (.not.(ilemx(i+1,j)).and..not.(ilemx(i+2,j))) then 
    256                      
    257                  ! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2 
    258  
    259                  phi_i = Uxbar(i,j) * Hmxloc(i,j) 
    260                  phi_xgl(i+1,j)  = phi_i + dx * (phi_prescr-phi_i)/(0.5 * dx + xpos) 
    261                  uxbartemp(i+1,j)   = uxbartemp(i+1,j) + phi_xgl(i+1,j) / Hmxloc(i+1,j) ! attention division par 0 possible ? 
    262                  countx   (i+1,j)   = countx(i+1,j)+1 
    263                  imx_diag (i+1,j) = 1 
    264                  gr_line_schoof(i+1,j) = 1 
    265  
    266                  ! extrapolation pour avoir uxbar(i+1,j) 
    267  
    268                  phi_ip3 = Uxbar(i+3,j) * Hmxloc(i+3,j) 
    269                  phi_xgl(i+2,j) = phi_ip3 + dx * (phi_prescr - phi_ip3)/(2.5 * dx  - xpos) 
    270                  uxbartemp(i+2,j)   = uxbartemp(i+2,j) + phi_xgl(i+2,j) / Hmxloc(i+2,j) 
    271                  countx   (i+2,j)   = countx(i+2,j)+1 
    272                  imx_diag (i+2,j) = 1 
    273                  gr_line_schoof(i+2,j) = 1 
    274  
    275 !afq --                 end if 
    276                
    277               end if 
    278  
    279            end if 
    280  
    281            !******************************************************************************* 
    282            ! selon y (indice j) 
    283            if (archimtab(i,j-1).LT.0..and.Uybar(i,j).LT.0..and..not.flot_marais(i,j-1)) then  
    284            !if (flot(i,j-1).and.Uybar(i,j).LT.0.) then    ! grounding line between j-1,j    CAS SUD (ecoulement vers grille SUD) 
    285  
    286  
    287               call  calc_xgl4schoof(-dy,alpha_flot,real(H(i,j)),real(Bsoc(i,j)),sealevel_2d(i,j),real(H(i,j-1)),real(Bsoc(i,j-1)),sealevel_2d(i,j-1),ypos,Hgl) 
    288               ypos_tab(i,j)=ypos 
    289               Hgly_tab(i,j)=Hgl 
    290                
    291               ! afq: the back force is on the staggered grid, the GL can be either South or North to this point. 
    292               if (ypos .lt. -dy/2.) then 
    293                  bfy = back_force_y(i,j) 
    294               else 
    295                  bfy = back_force_y(i,j+1) 
    296               endif 
    297               if (gr_select.eq.1) then ! flux de Tsai 
    298                  call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfy,phi_prescr) 
    299               else if (gr_select.eq.2) then ! flux de Schoof 
    300                  ! afq: the dragging coef. is on the staggered grid, the GL can be either South or North to this point. 
    301                  if (ypos .lt. -dy/2.) then 
    302                     frot_coef = betamy(i,j) 
    303                  else 
    304                     frot_coef = betamy(i,j+1) 
    305                  endif 
    306                  call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfy,phi_prescr) 
    307               else 
    308                  print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE' 
    309                  stop 
    310               endif 
    311               phi_prescr=-phi_prescr ! Doit être negatif dans le cas Sud 
    312               phi_prescr_taby(i,j)=phi_prescr 
    313               if (ypos .lt. -dy/2.) then  !  GL au sud du j staggered,                          o centre, x stag 
    314  
    315                  !  grille    ! .....x......o......x......o...G..x......O......x......o......x......o 
    316                  !            !            j-2           j-1            j            j+1           j+2 
    317                  !flux        !     in            out         G out            in 
    318  
    319 !afq --                 if (.not.(ilemy(i,j)).and..not.(ilemy(i,j-1))) then 
    320  
    321                  ! extrapolation pour avoir uybar(i,j-1) qui sera ensuite prescrit dans call rempli_L2 
    322  
    323                  phi_jm2 = Uybar(i,j-2) * Hmyloc(i,j-2) 
    324                  phi_ygl(i,j-1)  = phi_jm2 + dy * (phi_prescr-phi_jm2)/(2.5 * dy + ypos) 
    325                  uybartemp(i,j-1)   = uybartemp(i,j-1) + phi_ygl(i,j-1) / Hmyloc(i,j-1) ! attention division par 0 possible ? 
    326                  county   (i,j-1)   = county(i,j-1)+1 
    327                  imy_diag (i,j-1) = 1 
    328                  gr_line_schoof(i,j-1) = 1 
    329  
    330                  ! extrapolation pour avoir uybar(i,j) 
    331  
    332                  phi_jp1 = Uybar(i,j+1) * Hmyloc(i,j+1) 
    333                  phi_ygl(i,j) = phi_jp1 + dy * (phi_prescr - phi_jp1)/(0.5 * dy - ypos) 
    334                  uybartemp(i,j)   = uybartemp(i,j) + phi_ygl(i,j) / Hmyloc(i,j) 
    335                  county   (i,j)   = county(i,j)+1 
    336                  imy_diag (i,j) = 1 
    337                  gr_line_schoof(i,j) = 1 
    338  
    339  
    340 !afq --                 end if 
    341                
    342               else            !  GL au nord du j staggered,                          o centre, x stag 
    343  
    344                  !  grille    ! .....x......o......x......o......x...G..O......x......o......x......o 
    345                  !            j-2           j-1            j            j+1           j+2 
    346                  !flux        !     in             in           out  G         out           in 
    347  
    348 !afq --                 if (.not.(ilemy(i,j)).and..not.(ilemy(i,j+1))) then 
    349                      
    350                  ! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2 
    351  
    352                  phi_jm1 = Uybar(i,j-1) * Hmyloc(i,j-1) 
    353                  phi_ygl(i,j)  = phi_jm1 + dy * (phi_prescr-phi_jm1)/(1.5 * dy + ypos) 
    354                  uybartemp(i,j)   = uybartemp(i,j) + phi_ygl(i,j) / Hmyloc(i,j) ! attention division par 0 possible ? 
    355                  county   (i,j)   = county(i,j)+1 
    356                  imy_diag(i,j) = 1 
    357                  gr_line_schoof(i,j) = 1 
    358  
    359                  ! extrapolation pour avoir uybar(i,j+1) 
    360  
    361                  phi_jp2 = Uybar(i,j+2) * Hmyloc(i,j+2) 
    362                  phi_ygl(i,j+1) = phi_jp2 + dy * (phi_prescr - phi_jp2)/(1.5 * dy  - ypos) 
    363                  uybartemp(i,j+1)   = uybartemp(i,j+1) + phi_ygl(i,j+1) / Hmyloc(i,j+1) 
    364                  county   (i,j+1)   = county(i,j+1)+1 
    365                  imy_diag (i,j+1) = 1 
    366                  gr_line_schoof(i,j+1) = 1 
    367  
    368 !afq --                 end if 
    369                
    370               end if 
    371                
    372            else if (archimtab(i,j+1).LT.0..and.Uybar(i,j+1).GT.0..and..not.flot_marais(i,j+1)) then 
    373            !else if (flot(i,j+1).and.Uybar(i,j+1).GT.0.) then    ! grounding line between j,j+1    CAS NORD (ecoulement vers grille Nord) 
    374  
    375               call  calc_xgl4schoof(dy,alpha_flot,real(H(i,j)),real(Bsoc(i,j)),sealevel_2d(i,j),real(H(i,j+1)),real(Bsoc(i,j+1)),sealevel_2d(i,j+1),ypos,Hgl) 
    376               ypos_tab(i,j)=ypos 
    377               Hgly_tab(i,j)=Hgl 
    378                
    379               ! afq: the back force is on the staggered grid, the GL can be either South or North to this point. 
    380               if (ypos .lt. dy/2.) then 
    381                  bfy = back_force_y(i,j) 
    382               else 
    383                  bfy = back_force_y(i,j+1) 
    384               endif 
    385               if (gr_select.eq.1) then ! flux de Tsai 
    386                  call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfy,phi_prescr) 
    387               else if (gr_select.eq.2) then ! flux de Schoof 
    388                  ! afq: the dragging coef. is on the staggered grid, the GL can be either South or North to this point. 
    389                  if (ypos .lt. dy/2.) then 
    390                     frot_coef = betamy(i,j) 
    391                  else 
    392                     frot_coef = betamy(i,j+1) 
    393                  endif 
    394                  call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfy,phi_prescr) 
    395               else 
    396                  print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE' 
    397                  stop 
    398               endif 
    399               phi_prescr_taby(i,j)=phi_prescr 
    400               if (ypos .lt. dy/2.) then  !  GL au sud du j staggered,                          o centre, x stag 
    401  
    402                  !  grille    ! .....x......o......x......o......x......O..G...x......o......x......o 
    403                  !            !            j-2           j-1            j            j+1           j+2 
    404                  !flux        !                    in           out        G  out            in 
    405  
    406 !afq --                 if (.not.(ilemy(i,j)).and..not.(ilemy(i,j+1))) then 
    407  
    408                  ! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2 
    409  
    410                  phi_jm1 = Uybar(i,j-1) * Hmyloc(i,j-1) 
    411                  phi_ygl(i,j)  = phi_jm1 + dy * (phi_prescr-phi_jm1)/(1.5 * dy + ypos) 
    412                  uybartemp(i,j)   = uybartemp(i,j) + phi_ygl(i,j) / Hmyloc(i,j) ! attention division par 0 possible ? 
    413                  county   (i,j)   = county(i,j)+1 
    414                  imy_diag (i,j) = 1 
    415                  gr_line_schoof(i,j) = 1 
    416  
    417                  ! extrapolation pour avoir uybar(i,j+1) 
    418  
    419                  phi_jp2 = Uybar(i,j+2) * Hmyloc(i,j+2) 
    420                  phi_ygl(i,j+1) = phi_jp2 + dy * (phi_prescr - phi_jp2)/(1.5* dy - ypos) 
    421                  uybartemp(i,j+1)   = uybartemp(i,j+1) + phi_ygl(i,j+1) / Hmyloc(i,j+1) 
    422                  county   (i,j+1)   = county(i,j+1)+1 
    423                  imy_diag (i,j+1) = 1 
    424                  gr_line_schoof(i,j+1) = 1 
    425  
    426 !afq --                 end if 
    427                
    428               else   !  GL au nord du j staggered,                          o centre, x stag 
    429  
    430                  !  grille    ! ......x......o......x......O......x...G..o......x......o......x.....o 
    431                  !            !             j-1            j            j+1           j+2          j+3 
    432                  !flux        !                     in           out  G        out            in 
    433  
    434 !afq --                 if (.not.(ilemy(i,j+1)).and..not.(ilemy(i,j+2))) then 
    435  
    436                  ! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2 
    437  
    438                  phi_j = Uybar(i,j) * Hmyloc(i,j) 
    439                  phi_ygl(i,j+1)  = phi_j + dy * (phi_prescr-phi_j)/(0.5 * dy + ypos) 
    440                  uybartemp(i,j+1)   = uybartemp(i,j+1) +  phi_ygl(i,j+1) / Hmyloc(i,j+1) ! attention division par 0 possible ? 
    441                  county   (i,j+1)   = county(i,j+1)+1 
    442                  imy_diag (i,j+1) = 1 
    443                  gr_line_schoof(i,j+1) = 1 
    444  
    445                  ! extrapolation pour avoir uybar(i,j+1) 
    446  
    447                  phi_jp3 = Uybar(i,j+3) * Hmyloc(i,j+3) 
    448                  phi_ygl(i,j+2) = phi_jp3 + dy * (phi_prescr - phi_jp3)/(2.5 * dy  - ypos) 
    449                  uybartemp(i,j+2)   = uybartemp(i,j+2) +  phi_ygl(i,j+2) / Hmyloc(i,j+2) 
    450                  county   (i,j+2)   = county(i,j+2)+1 
    451                  imy_diag (i,j+2) = 1 
    452                  gr_line_schoof(i,j+2) = 1 
    453  
    454 !afq --                 end if 
    455  
    456               end if 
    457  
    458            end if 
    459         end if 
    460      end do 
    461   end do 
    462  
    463   ! afq -- if we discard the points with multiple fluxes coming, uncomment: 
    464   where (countx(:,:).ge.2) 
    465      uxbartemp(:,:)=uxbar(:,:)*countx(:,:) 
    466      imx_diag(:,:) = imx_diag_in(:,:) 
    467   end where 
    468   where (county(:,:).ge.2) 
    469      uybartemp(:,:)=uybar(:,:)*county(:,:) 
    470      imy_diag(:,:) = imx_diag_in(:,:) 
    471   end where 
    472   ! afq -- 
    473    
    474   where (countx(:,:).ne.0) uxbar(:,:)= max(min(uxbartemp(:,:)/countx(:,:),V_limit),-V_limit) 
    475   where (county(:,:).ne.0) uybar(:,:)= max(min(uybartemp(:,:)/county(:,:),V_limit),-V_limit) 
    476             
    477   do j= 3,ny-3 
    478      do i=3,nx-3 
    479         denom = sqrt((uxbartemp(i,j)**2 + uybartemp(i,j)**2))*sqrt((uxbar(i,j)**2 + uybar(i,j)**2)) 
    480         if (denom.gt.toutpetit) then 
    481            ! prodscal is cos theta in : U . V = u v cos theta 
    482            ! this number is between -1 and 1, 1 = same direction & -1 = opposite direction 
    483            ! we could use this as an artificial back force...  
    484            prodscal = (uxbartemp(i,j)*uxbar(i,j) + uybartemp(i,j)*uybar(i,j))/  & 
    485                 (sqrt((uxbartemp(i,j)**2 + uybartemp(i,j)**2))*sqrt((uxbar(i,j)**2 + uybar(i,j)**2))) 
    486         endif 
    487      end do 
    488   end do 
    489    
     20  subroutine init_furst_schoof 
     21     
     22    use module3d_phy, only:num_param,num_rep_42 
     23    use param_phy_mod, only:row,ro 
     24 
     25    implicit none 
     26 
     27    namelist/furst_schoof/frot_coef,gr_select 
     28 
     29    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
     30    read(num_param,furst_schoof) 
     31    write(num_rep_42,'(A)')'!___________________________________________________________'  
     32    write(num_rep_42,'(A)') '&furst_schoof                                ! nom du bloc ' 
     33    write(num_rep_42,*) 
     34    write(num_rep_42,*)'frot_coef  = ',frot_coef 
     35    write(num_rep_42,*)'gr_select  = ',gr_select 
     36    write(num_rep_42,*)'/' 
     37    write(num_rep_42,*)'! frot_coef : solid friction law coef 0.6 in GMD 2018' 
     38    write(num_rep_42,*)'! gr_select = 1 : Tsai , 2 : Schoof' 
     39 
     40    !  frot = 0.035 ! f in the solid friction law, 0.035 in Ritz et al. 2001     
     41 
     42    ! choix Tsai (loi de Coulomb) ou Schoof (loi de Weertman) 
     43    alpha_flot= ro/row 
     44    back_force_x(:,:)=0.1 
     45    back_force_y(:,:)=0.1 
     46 
     47  end subroutine init_furst_schoof 
     48 
     49  subroutine interpol_glflux 
     50     
     51    use module3d_phy, only:hmx,hmy,gr_line_schoof,Bsoc,H,sealevel_2d,imx_diag,imy_diag,& 
     52         uxbar,uybar,flot_marais,dx,dy,Abar,betamx,betamy,V_limit,debug_3D 
     53    use deformation_mod_2lois, only:glen 
     54     
     55    implicit none 
     56 
     57    integer :: i,j 
     58    real :: xpos,ypos 
     59    real :: Hgl 
     60    !  integer :: igl 
     61 
     62    real :: phi_im2, phi_im1, phi_i, phi_ip1, phi_ip2, phi_ip3 
     63    real :: phi_jm2, phi_jm1, phi_j, phi_jp1, phi_jp2, phi_jp3 
     64    real,dimension(nx,ny) :: phi_xgl, phi_ygl ! flux gr line sur maille stag 
     65    integer,dimension(nx,ny) :: countx, county   ! how often do we modify ux/uy 
     66    real :: phi_prescr 
     67    real :: toutpetit = 1e-6 
     68    real :: denom, prodscal 
     69 
     70    real :: bfx, bfy 
     71 
     72    !debug 
     73    real,dimension(nx,ny) :: xpos_tab=9999. 
     74    real,dimension(nx,ny) :: ypos_tab=9999. 
     75    real,dimension(nx,ny) :: Hglx_tab=9999. 
     76    real,dimension(nx,ny) :: Hgly_tab=9999. 
     77    real,dimension(nx,ny) :: Hmxloc, Hmyloc ! pour eviter division par 0 
     78    real,dimension(nx,ny) :: phi_prescr_tabx, phi_prescr_taby 
     79 
     80    real,dimension(nx,ny) :: uxbartemp, uybartemp 
     81 
     82    real,dimension(nx,ny) :: archimtab 
     83 
     84    real,dimension(nx,ny) :: imx_diag_in 
     85 
     86    phi_prescr_tabx=0. 
     87    phi_prescr_taby=0. 
     88 
     89    Hmxloc(:,:)= max(Hmx(:,:),1.) 
     90    Hmyloc(:,:)=max(Hmy(:,:),1.) 
     91 
     92    uxbartemp(:,:)=0. 
     93    uybartemp(:,:)=0. 
     94    countx(:,:)=0. 
     95    county(:,:)=0. 
     96    gr_line_schoof(:,:)=0 
     97 
     98    archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot -sealevel_2d(:,:) 
     99 
     100    imx_diag_in(:,:) = imx_diag(:,:) 
     101 
     102    ! detection des noeuds grounding line 
     103    do j= 3,ny-3 
     104       do i=3,nx-3 
     105          !archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel 
     106          if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel_2d(i,j))) then    ! grounded and with ice 
     107             !           print*,'schoof gr line ij',i,j 
     108             ! selon x (indice i) 
     109 
     110             if (archimtab(i-1,j).LT.0..and.Uxbar(i,j).LT.0..and..not.flot_marais(i-1,j)) then 
     111                !if (flot(i-1,j).and.Uxbar(i,j).LT.0.) then    ! grounding line between i-1,i    CAS WEST (ecoulement vers grille West) 
     112 
     113                call  calc_xgl4schoof(-dx,alpha_flot,real(H(i,j)),real(Bsoc(i,j)),sealevel_2d(i,j),real(H(i-1,j)),real(Bsoc(i-1,j)),sealevel_2d(i-1,j),xpos,Hgl) 
     114                xpos_tab(i,j)=xpos 
     115                Hglx_tab(i,j)=Hgl 
     116 
     117                ! afq: the back force is on the staggered grid, the GL can be either West or East to this point. 
     118                if (xpos .lt. -dx/2.) then 
     119                   bfx = back_force_x(i,j) 
     120                else 
     121                   bfx = back_force_x(i+1,j) 
     122                endif 
     123                if (gr_select.eq.1) then ! flux de Tsai 
     124                   call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfx,phi_prescr) 
     125                else if (gr_select.eq.2) then ! flux de Schoof 
     126                   ! afq: the dragging coef. is on the staggered grid, the GL can be either West or East to this point. 
     127                   if (xpos .lt. -dx/2.) then 
     128                      frot_coef = betamx(i,j) 
     129                   else 
     130                      frot_coef = betamx(i+1,j) 
     131                   endif 
     132                   call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfx,phi_prescr) 
     133                else 
     134                   print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE' 
     135                   stop 
     136                endif 
     137                phi_prescr=-phi_prescr ! Doit être negatif dans le cas West 
     138                phi_prescr_tabx(i,j)=phi_prescr 
     139 
     140                if (xpos .lt. -dx/2.) then  !  GL a west du i staggered,                          o centre, x stag 
     141 
     142                   !  grille    ! .....x......o......x......o...G..x......O......x......o......x......o 
     143                   !            !            i-2           i-1            i            i+1           i+2 
     144                   !flux        !     in            out         G out            in 
     145 
     146 
     147                   !afq --                 if (.not.(ilemx(i,j)).and..not.(ilemx(i-1,j))) then 
     148 
     149                   ! extrapolation pour avoir uxbar(i-1,j) qui sera ensuite prescrit dans call rempli_L2 
     150 
     151                   phi_im2 = Uxbar(i-2,j) * Hmxloc(i-2,j) 
     152                   phi_xgl(i-1,j)  = phi_im2 + dx * (phi_prescr-phi_im2)/(2.5 * dx + xpos) 
     153                   uxbartemp(i-1,j)   = uxbartemp(i-1,j)+ phi_xgl(i-1,j) / Hmxloc(i-1,j) ! attention division par 0 possible ? 
     154                   countx   (i-1,j)   = countx(i-1,j)+1 
     155                   imx_diag (i-1,j)   = 1 
     156                   gr_line_schoof(i-1,j) = 1 
     157 
     158                   ! extrapolation pour avoir uxbar(i,j) 
     159 
     160                   phi_ip1 = Uxbar(i+1,j) * Hmxloc(i+1,j) 
     161                   phi_xgl(i,j) = phi_ip1 + dx * (phi_prescr - phi_ip1)/(0.5 * dx - xpos) 
     162                   !                 print*,'uxbar(i,j)=',uxbar(i,j) 
     163                   uxbartemp(i,j)   = uxbartemp(i,j) + phi_xgl(i,j) / Hmxloc(i,j) 
     164                   countx   (i,j)   = countx(i,j)+1 
     165                   imx_diag (i,j)   = 1 
     166                   gr_line_schoof(i,j) = 1 
     167                   !                 print*,'schoof uxbar(i,j)=',uxbar(i,j) 
     168                   !                 print*,'Hgl',Hgl 
     169                   !                 print*,'phi_prescr',phi_prescr, phi_prescr/Hgl 
     170                   !                 print* 
     171 
     172                   !afq --                 end if 
     173 
     174                else            !  GL a Est du i staggered,                          o centre, x stag 
     175 
     176                   !  grille    ! .....x......o......x......o......x...G..O......x......o......x......o 
     177                   !            i-2           i-1            i            i+1           i+2 
     178                   !flux        !     in             in           out  G         out           in 
     179 
     180                   !afq --                 if (.not.(ilemx(i,j)).and..not.(ilemx(i+1,j))) then 
     181 
     182                   ! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2 
     183 
     184                   phi_im1 = Uxbar(i-1,j) * Hmxloc(i-1,j) 
     185                   phi_xgl(i,j)  = phi_im1 + dx * (phi_prescr-phi_im1)/(1.5 * dx + xpos) 
     186                   uxbartemp(i,j)   = uxbartemp(i,j) + phi_xgl(i,j) / Hmxloc(i,j) ! attention division par 0 possible ? 
     187                   countx   (i,j)   = countx(i,j)+1 
     188                   imx_diag (i,j) = 1 
     189                   gr_line_schoof(i,j) = 1 
     190 
     191                   ! extrapolation pour avoir uxbar(i+1,j) 
     192 
     193                   phi_ip2 = Uxbar(i+2,j) * Hmxloc(i+2,j) 
     194                   phi_xgl(i+1,j) = phi_ip2 + dx * (phi_prescr - phi_ip2)/(1.5 * dx  - xpos) 
     195                   uxbartemp(i+1,j)   = uxbartemp(i+1,j) +  phi_xgl(i+1,j) / Hmxloc(i+1,j) 
     196                   countx   (i+1,j)   = countx(i+1,j)+1 
     197                   imx_diag (i+1,j)   = 1 
     198                   gr_line_schoof(i+1,j) = 1 
     199 
     200                   !afq --                 end if 
     201 
     202                end if 
     203 
     204             else if (archimtab(i+1,j).LT.0..and.Uxbar(i+1,j).GT.0..and..not.flot_marais(i+1,j)) then 
     205                !else if (flot(i+1,j).and.Uxbar(i+1,j).GT.0.) then    ! grounding line between i,i+1    CAS EST (ecoulement vers grille Est) 
     206 
     207                call  calc_xgl4schoof(dx,alpha_flot,real(H(i,j)),real(Bsoc(i,j)),sealevel_2d(i,j),real(H(i+1,j)),real(Bsoc(i+1,j)),sealevel_2d(i+1,j),xpos,Hgl) 
     208                xpos_tab(i,j)=xpos 
     209                Hglx_tab(i,j)=Hgl 
     210 
     211                ! afq: the back force is on the staggered grid, the GL can be either West or East to this point. 
     212                if (xpos .lt. dx/2.) then 
     213                   bfx = back_force_x(i,j) 
     214                else 
     215                   bfx = back_force_x(i+1,j) 
     216                endif 
     217                if (gr_select.eq.1) then ! flux de Tsai 
     218                   call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfx,phi_prescr) 
     219                else if (gr_select.eq.2) then ! flux de Schoof 
     220                   ! afq: the dragging coef. is on the staggered grid, the GL can be either West or East to this point. 
     221                   if (xpos .lt. dx/2.) then 
     222                      frot_coef = betamx(i,j) 
     223                   else 
     224                      frot_coef = betamx(i+1,j) 
     225                   endif 
     226                   call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfx,phi_prescr) 
     227                else 
     228                   print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE' 
     229                   stop 
     230                endif 
     231                phi_prescr_tabx(i,j)=phi_prescr 
     232                if (xpos .lt. dx/2.) then  !  GL a west du i staggered,                          o centre, x stag 
     233 
     234                   !  grille    ! .....x......o......x......o......x......O..G...x......o......x......o 
     235                   !            !            i-2           i-1            i            i+1           i+2 
     236                   !flux        !                    in           out        G  out            in 
     237 
     238                   !afq --                 if (.not.(ilemx(i,j)).and..not.(ilemx(i+1,j))) then 
     239 
     240                   ! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2 
     241 
     242                   phi_im1 = Uxbar(i-1,j) * Hmxloc(i-1,j) 
     243                   phi_xgl(i,j)  = phi_im1 + dx * (phi_prescr-phi_im1)/(1.5 * dx + xpos) 
     244                   uxbartemp(i,j)   =  uxbartemp(i,j) + phi_xgl(i,j) / Hmxloc(i,j) ! attention division par 0 possible ? 
     245                   countx   (i,j)   = countx(i,j)+1 
     246                   imx_diag (i,j) = 1 
     247                   gr_line_schoof(i,j) = 1 
     248 
     249                   ! extrapolation pour avoir uxbar(i+1,j) 
     250 
     251                   phi_ip2 = Uxbar(i+2,j) * Hmxloc(i+2,j) 
     252                   phi_xgl(i+1,j) = phi_ip2 + dx * (phi_prescr - phi_ip2)/(1.5* dx - xpos) 
     253                   uxbartemp(i+1,j)   = uxbartemp(i+1,j) + phi_xgl(i+1,j) / Hmxloc(i+1,j) 
     254                   countx   (i+1,j)   = countx(i+1,j)+1 
     255                   imx_diag (i+1,j) = 1 
     256                   gr_line_schoof(i+1,j) = 1 
     257 
     258                   !afq --                 end if 
     259 
     260                else   !  GL a west du i staggered,                          o centre, x stag 
     261 
     262                   !  grille    ! ......x......o......x......O......x...G..o......x......o......x.....o 
     263                   !            !             i-1            i            i+1           i+2          i+3 
     264                   !flux        !                     in           out  G        out            in 
     265 
     266                   !afq --                 if (.not.(ilemx(i+1,j)).and..not.(ilemx(i+2,j))) then 
     267 
     268                   ! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2 
     269 
     270                   phi_i = Uxbar(i,j) * Hmxloc(i,j) 
     271                   phi_xgl(i+1,j)  = phi_i + dx * (phi_prescr-phi_i)/(0.5 * dx + xpos) 
     272                   uxbartemp(i+1,j)   = uxbartemp(i+1,j) + phi_xgl(i+1,j) / Hmxloc(i+1,j) ! attention division par 0 possible ? 
     273                   countx   (i+1,j)   = countx(i+1,j)+1 
     274                   imx_diag (i+1,j) = 1 
     275                   gr_line_schoof(i+1,j) = 1 
     276 
     277                   ! extrapolation pour avoir uxbar(i+1,j) 
     278 
     279                   phi_ip3 = Uxbar(i+3,j) * Hmxloc(i+3,j) 
     280                   phi_xgl(i+2,j) = phi_ip3 + dx * (phi_prescr - phi_ip3)/(2.5 * dx  - xpos) 
     281                   uxbartemp(i+2,j)   = uxbartemp(i+2,j) + phi_xgl(i+2,j) / Hmxloc(i+2,j) 
     282                   countx   (i+2,j)   = countx(i+2,j)+1 
     283                   imx_diag (i+2,j) = 1 
     284                   gr_line_schoof(i+2,j) = 1 
     285 
     286                   !afq --                 end if 
     287 
     288                end if 
     289 
     290             end if 
     291 
     292             !******************************************************************************* 
     293             ! selon y (indice j) 
     294             if (archimtab(i,j-1).LT.0..and.Uybar(i,j).LT.0..and..not.flot_marais(i,j-1)) then  
     295                !if (flot(i,j-1).and.Uybar(i,j).LT.0.) then    ! grounding line between j-1,j    CAS SUD (ecoulement vers grille SUD) 
     296 
     297 
     298                call  calc_xgl4schoof(-dy,alpha_flot,real(H(i,j)),real(Bsoc(i,j)),sealevel_2d(i,j),real(H(i,j-1)),real(Bsoc(i,j-1)),sealevel_2d(i,j-1),ypos,Hgl) 
     299                ypos_tab(i,j)=ypos 
     300                Hgly_tab(i,j)=Hgl 
     301 
     302                ! afq: the back force is on the staggered grid, the GL can be either South or North to this point. 
     303                if (ypos .lt. -dy/2.) then 
     304                   bfy = back_force_y(i,j) 
     305                else 
     306                   bfy = back_force_y(i,j+1) 
     307                endif 
     308                if (gr_select.eq.1) then ! flux de Tsai 
     309                   call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfy,phi_prescr) 
     310                else if (gr_select.eq.2) then ! flux de Schoof 
     311                   ! afq: the dragging coef. is on the staggered grid, the GL can be either South or North to this point. 
     312                   if (ypos .lt. -dy/2.) then 
     313                      frot_coef = betamy(i,j) 
     314                   else 
     315                      frot_coef = betamy(i,j+1) 
     316                   endif 
     317                   call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfy,phi_prescr) 
     318                else 
     319                   print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE' 
     320                   stop 
     321                endif 
     322                phi_prescr=-phi_prescr ! Doit être negatif dans le cas Sud 
     323                phi_prescr_taby(i,j)=phi_prescr 
     324                if (ypos .lt. -dy/2.) then  !  GL au sud du j staggered,                          o centre, x stag 
     325 
     326                   !  grille    ! .....x......o......x......o...G..x......O......x......o......x......o 
     327                   !            !            j-2           j-1            j            j+1           j+2 
     328                   !flux        !     in            out         G out            in 
     329 
     330                   !afq --                 if (.not.(ilemy(i,j)).and..not.(ilemy(i,j-1))) then 
     331 
     332                   ! extrapolation pour avoir uybar(i,j-1) qui sera ensuite prescrit dans call rempli_L2 
     333 
     334                   phi_jm2 = Uybar(i,j-2) * Hmyloc(i,j-2) 
     335                   phi_ygl(i,j-1)  = phi_jm2 + dy * (phi_prescr-phi_jm2)/(2.5 * dy + ypos) 
     336                   uybartemp(i,j-1)   = uybartemp(i,j-1) + phi_ygl(i,j-1) / Hmyloc(i,j-1) ! attention division par 0 possible ? 
     337                   county   (i,j-1)   = county(i,j-1)+1 
     338                   imy_diag (i,j-1) = 1 
     339                   gr_line_schoof(i,j-1) = 1 
     340 
     341                   ! extrapolation pour avoir uybar(i,j) 
     342 
     343                   phi_jp1 = Uybar(i,j+1) * Hmyloc(i,j+1) 
     344                   phi_ygl(i,j) = phi_jp1 + dy * (phi_prescr - phi_jp1)/(0.5 * dy - ypos) 
     345                   uybartemp(i,j)   = uybartemp(i,j) + phi_ygl(i,j) / Hmyloc(i,j) 
     346                   county   (i,j)   = county(i,j)+1 
     347                   imy_diag (i,j) = 1 
     348                   gr_line_schoof(i,j) = 1 
     349 
     350 
     351                   !afq --                 end if 
     352 
     353                else            !  GL au nord du j staggered,                          o centre, x stag 
     354 
     355                   !  grille    ! .....x......o......x......o......x...G..O......x......o......x......o 
     356                   !            j-2           j-1            j            j+1           j+2 
     357                   !flux        !     in             in           out  G         out           in 
     358 
     359                   !afq --                 if (.not.(ilemy(i,j)).and..not.(ilemy(i,j+1))) then 
     360 
     361                   ! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2 
     362 
     363                   phi_jm1 = Uybar(i,j-1) * Hmyloc(i,j-1) 
     364                   phi_ygl(i,j)  = phi_jm1 + dy * (phi_prescr-phi_jm1)/(1.5 * dy + ypos) 
     365                   uybartemp(i,j)   = uybartemp(i,j) + phi_ygl(i,j) / Hmyloc(i,j) ! attention division par 0 possible ? 
     366                   county   (i,j)   = county(i,j)+1 
     367                   imy_diag(i,j) = 1 
     368                   gr_line_schoof(i,j) = 1 
     369 
     370                   ! extrapolation pour avoir uybar(i,j+1) 
     371 
     372                   phi_jp2 = Uybar(i,j+2) * Hmyloc(i,j+2) 
     373                   phi_ygl(i,j+1) = phi_jp2 + dy * (phi_prescr - phi_jp2)/(1.5 * dy  - ypos) 
     374                   uybartemp(i,j+1)   = uybartemp(i,j+1) + phi_ygl(i,j+1) / Hmyloc(i,j+1) 
     375                   county   (i,j+1)   = county(i,j+1)+1 
     376                   imy_diag (i,j+1) = 1 
     377                   gr_line_schoof(i,j+1) = 1 
     378 
     379                   !afq --                 end if 
     380 
     381                end if 
     382 
     383             else if (archimtab(i,j+1).LT.0..and.Uybar(i,j+1).GT.0..and..not.flot_marais(i,j+1)) then 
     384                !else if (flot(i,j+1).and.Uybar(i,j+1).GT.0.) then    ! grounding line between j,j+1    CAS NORD (ecoulement vers grille Nord) 
     385 
     386                call  calc_xgl4schoof(dy,alpha_flot,real(H(i,j)),real(Bsoc(i,j)),sealevel_2d(i,j),real(H(i,j+1)),real(Bsoc(i,j+1)),sealevel_2d(i,j+1),ypos,Hgl) 
     387                ypos_tab(i,j)=ypos 
     388                Hgly_tab(i,j)=Hgl 
     389 
     390                ! afq: the back force is on the staggered grid, the GL can be either South or North to this point. 
     391                if (ypos .lt. dy/2.) then 
     392                   bfy = back_force_y(i,j) 
     393                else 
     394                   bfy = back_force_y(i,j+1) 
     395                endif 
     396                if (gr_select.eq.1) then ! flux de Tsai 
     397                   call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfy,phi_prescr) 
     398                else if (gr_select.eq.2) then ! flux de Schoof 
     399                   ! afq: the dragging coef. is on the staggered grid, the GL can be either South or North to this point. 
     400                   if (ypos .lt. dy/2.) then 
     401                      frot_coef = betamy(i,j) 
     402                   else 
     403                      frot_coef = betamy(i,j+1) 
     404                   endif 
     405                   call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfy,phi_prescr) 
     406                else 
     407                   print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE' 
     408                   stop 
     409                endif 
     410                phi_prescr_taby(i,j)=phi_prescr 
     411                if (ypos .lt. dy/2.) then  !  GL au sud du j staggered,                          o centre, x stag 
     412 
     413                   !  grille    ! .....x......o......x......o......x......O..G...x......o......x......o 
     414                   !            !            j-2           j-1            j            j+1           j+2 
     415                   !flux        !                    in           out        G  out            in 
     416 
     417                   !afq --                 if (.not.(ilemy(i,j)).and..not.(ilemy(i,j+1))) then 
     418 
     419                   ! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2 
     420 
     421                   phi_jm1 = Uybar(i,j-1) * Hmyloc(i,j-1) 
     422                   phi_ygl(i,j)  = phi_jm1 + dy * (phi_prescr-phi_jm1)/(1.5 * dy + ypos) 
     423                   uybartemp(i,j)   = uybartemp(i,j) + phi_ygl(i,j) / Hmyloc(i,j) ! attention division par 0 possible ? 
     424                   county   (i,j)   = county(i,j)+1 
     425                   imy_diag (i,j) = 1 
     426                   gr_line_schoof(i,j) = 1 
     427 
     428                   ! extrapolation pour avoir uybar(i,j+1) 
     429 
     430                   phi_jp2 = Uybar(i,j+2) * Hmyloc(i,j+2) 
     431                   phi_ygl(i,j+1) = phi_jp2 + dy * (phi_prescr - phi_jp2)/(1.5* dy - ypos) 
     432                   uybartemp(i,j+1)   = uybartemp(i,j+1) + phi_ygl(i,j+1) / Hmyloc(i,j+1) 
     433                   county   (i,j+1)   = county(i,j+1)+1 
     434                   imy_diag (i,j+1) = 1 
     435                   gr_line_schoof(i,j+1) = 1 
     436 
     437                   !afq --                 end if 
     438 
     439                else   !  GL au nord du j staggered,                          o centre, x stag 
     440 
     441                   !  grille    ! ......x......o......x......O......x...G..o......x......o......x.....o 
     442                   !            !             j-1            j            j+1           j+2          j+3 
     443                   !flux        !                     in           out  G        out            in 
     444 
     445                   !afq --                 if (.not.(ilemy(i,j+1)).and..not.(ilemy(i,j+2))) then 
     446 
     447                   ! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2 
     448 
     449                   phi_j = Uybar(i,j) * Hmyloc(i,j) 
     450                   phi_ygl(i,j+1)  = phi_j + dy * (phi_prescr-phi_j)/(0.5 * dy + ypos) 
     451                   uybartemp(i,j+1)   = uybartemp(i,j+1) +  phi_ygl(i,j+1) / Hmyloc(i,j+1) ! attention division par 0 possible ? 
     452                   county   (i,j+1)   = county(i,j+1)+1 
     453                   imy_diag (i,j+1) = 1 
     454                   gr_line_schoof(i,j+1) = 1 
     455 
     456                   ! extrapolation pour avoir uybar(i,j+1) 
     457 
     458                   phi_jp3 = Uybar(i,j+3) * Hmyloc(i,j+3) 
     459                   phi_ygl(i,j+2) = phi_jp3 + dy * (phi_prescr - phi_jp3)/(2.5 * dy  - ypos) 
     460                   uybartemp(i,j+2)   = uybartemp(i,j+2) +  phi_ygl(i,j+2) / Hmyloc(i,j+2) 
     461                   county   (i,j+2)   = county(i,j+2)+1 
     462                   imy_diag (i,j+2) = 1 
     463                   gr_line_schoof(i,j+2) = 1 
     464 
     465                   !afq --                 end if 
     466 
     467                end if 
     468 
     469             end if 
     470          end if 
     471       end do 
     472    end do 
     473 
     474    ! afq -- if we discard the points with multiple fluxes coming, uncomment: 
     475    where (countx(:,:).ge.2) 
     476       uxbartemp(:,:)=uxbar(:,:)*countx(:,:) 
     477       imx_diag(:,:) = imx_diag_in(:,:) 
     478    end where 
     479    where (county(:,:).ge.2) 
     480       uybartemp(:,:)=uybar(:,:)*county(:,:) 
     481       imy_diag(:,:) = imx_diag_in(:,:) 
     482    end where 
     483    ! afq -- 
     484 
     485    where (countx(:,:).ne.0) uxbar(:,:)= max(min(uxbartemp(:,:)/countx(:,:),V_limit),-V_limit) 
     486    where (county(:,:).ne.0) uybar(:,:)= max(min(uybartemp(:,:)/county(:,:),V_limit),-V_limit) 
     487 
     488    do j= 3,ny-3 
     489       do i=3,nx-3 
     490          denom = sqrt((uxbartemp(i,j)**2 + uybartemp(i,j)**2))*sqrt((uxbar(i,j)**2 + uybar(i,j)**2)) 
     491          if (denom.gt.toutpetit) then 
     492             ! prodscal is cos theta in : U . V = u v cos theta 
     493             ! this number is between -1 and 1, 1 = same direction & -1 = opposite direction 
     494             ! we could use this as an artificial back force...  
     495             prodscal = (uxbartemp(i,j)*uxbar(i,j) + uybartemp(i,j)*uybar(i,j))/  & 
     496                  (sqrt((uxbartemp(i,j)**2 + uybartemp(i,j)**2))*sqrt((uxbar(i,j)**2 + uybar(i,j)**2))) 
     497          endif 
     498       end do 
     499    end do 
     500 
    490501!!$  do j=1,ny 
    491502!!$     do i=1,nx 
     
    500511!!$  enddo 
    501512 
    502   debug_3D(:,:,66) = phi_prescr_tabx(:,:) 
    503   debug_3D(:,:,67) = phi_prescr_taby(:,:) 
    504    
    505 end subroutine interpol_glflux 
    506    
    507 !------------------------------------------------------------------------- 
    508 ! calcule la position sous maille de la ligne dechouage 
    509  subroutine calc_xgl4schoof(dy,alpha,H_0,B_0,SL_0,H_1,B_1,SL_1,xpos,Hgl) 
    510                  
    511  
    512   implicit none 
    513  
    514   real,intent(in)        ::    dy        !< longueur orientee de la maille  
    515   real,intent(in)        ::    alpha     !< coefficient de flottaison = rho/rhow 
    516   real,intent(in)        ::    H_0       !< epaisseur au point pose 
    517   real,intent(in)        ::    B_0       !< altitude socle au point pose 
    518   real,intent(in)        ::    SL_0      !< sea level au point pose 
    519   real,intent(in)        ::    H_1       !< epaisseur au point flottant 
    520   real,intent(in)        ::    B_1       !< altitude socle au point flottant 
    521   real,intent(in)        ::    SL_1      !< sea level au point flottant 
    522   real,intent(out)       ::    xpos      !< position de la ligne (en distance depuis le point pose) 
    523   real,intent(out)       ::    Hgl       !< epaisseur de glace a la grounding line 
    524     
    525    
    526  
    527   real                   ::    Cgl       ! variable de travail   
    528  
    529 !~   Cgl   = (-alpha * (H_1-H_0) - (B_1-B_0)) ! -(rho/rhow (H1-H0) + (B1-B0)) 
    530  
    531 !~   if (abs(Cgl).gt.1.e-5) then    
    532 !~      xpos    = (B_0 + alpha * H_0 - SL) / Cgl  
    533 !~   else 
    534 !~      xpos    = 1. - 1./abs(dy)   ! verifier 
    535 !~   end if 
    536    
    537 !~   if (xpos.LT.0..OR.xpos.GT.1.) then 
    538 !~              print*,'calc_xgl4schoof : xpos value error i,j=',i,j 
    539 !~              print*, 'xpos,Cgl', xpos,Cgl 
    540 !~              print*,'B_0, alpha, H_0, SL', B_0, alpha, H_0, SL 
    541 !~              stop 
    542 !~      endif 
    543                  
    544   Cgl   = (B_1-B_0) + alpha * (H_1-H_0) - (SL_1-SL_0) 
    545  
    546   if (abs(Cgl).gt.1.e-5) then    
    547      xpos    = (SL_0 - (B_0 + alpha * H_0)) / Cgl  
    548   else 
    549      xpos    = 0.5  ! verifier 
    550   end if 
    551    
    552 !  print* 
    553 !  print*, 'xpos', xpos 
    554 !       print*,'B_0, alpha, H_0, SL', B_0, alpha, H_0, SL 
    555 !       print*,'B_1, H_1', B_1, H_1 
    556  
    557          
    558   if (xpos.LT.0..OR.xpos.GT.1.) then 
    559      print*,'calc_xgl4schoof : xpos value error i,j=',i,j 
    560      print*, 'xpos,Cgl', xpos,Cgl 
    561      print*,'B_0, alpha, H_0, SL_0', B_0, alpha, H_0, SL_0 
    562      print*,'archim=',B_1+H_1*alpha - SL_1 
    563      !stop 
    564      xpos = min(max(0.,xpos),1.) 
    565   endif 
    566          
    567   xpos = xpos * dy 
    568    
    569   Hgl= H_0 + (H_1 - H_0) * xpos/dy 
    570  
    571 end subroutine calc_xgl4schoof 
    572  
    573  
    574 !-------------------------------------------------------------------- 
    575 ! Schoof_flux : calcule le coefficient et le flux de Schoof  
    576 !--------------------------------------------------------------------  
    577  
    578 ! Dans un premier temps, pas de terme horizontal (regarder dans toy_recul une tentative de confinement) 
    579  
    580 subroutine flux_Schoof4Schoof (Hgl,A_mean,C_frot,alpha,n_glen,m_weert,back_force_coef,phi_schoof) 
    581  
    582    
    583 implicit none 
    584  
    585 real,intent(in) :: Hgl 
    586 real,intent(in) :: A_mean 
    587 real,intent(in) :: C_frot 
    588 real,intent(in) :: alpha 
    589 real,intent(in) :: n_glen 
    590 real,intent(in) :: m_weert 
    591 real,intent(in) :: back_force_coef 
    592 real,intent(out) :: phi_schoof 
    593  
    594 ! afq: in standard GRISLI (19/05/17), tob= -beta ub => m_weert = 1 and C_frot = beta 
    595  
    596 !phi_schoof = (((A_mean * rog**(n_glen+1) * (1 - alpha)**n_glen) / (4**n_glen *C_frot)) **(1./(1.+1./m_weert)))*back_force_coef**(n_glen/(1+1./m_weert)) 
    597 !phi_schoof = phi_schoof * Hgl**((n_glen+3.+1/m_weert)/(1.+1./m_weert)) 
    598 phi_schoof = (((A_mean * rog**(n_glen+1.) * (1. - alpha)**n_glen) / (4.**n_glen*C_frot)) **(1./(1.+inv_mweert)))*back_force_coef**(n_glen/(1.+inv_mweert)) 
    599 phi_schoof = phi_schoof * Hgl**((n_glen+3.+inv_mweert)/(1.+inv_mweert)) 
    600  
    601  
    602  
    603 end subroutine flux_Schoof4Schoof 
    604  
    605 !-------------------------------------------------------------------- 
    606 ! Tsai_flux : calcule le coefficient et le flux de Schoof  
    607 !--------------------------------------------------------------------  
    608  
    609 ! Dans un premier temps, pas de terme horizontal (regarder dans toy_recul une tentative de confinement) 
    610  
    611 subroutine flux_Tsai4Schoof (Hgl,A_mean,f_frot,alpha,n_glen,back_force_coef,phi_Tsai) 
    612  
    613 implicit none 
    614  
    615 real,intent(in) :: Hgl 
    616 real,intent(in) :: A_mean 
    617 real,intent(in) :: f_frot 
    618 real,intent(in) :: alpha 
    619 real,intent(in) :: n_glen 
    620 real,intent(in) :: back_force_coef 
    621 real,intent(out) :: phi_Tsai 
    622  
    623 real, parameter :: Q0 = 0.61 
    624  
    625  
    626 phi_Tsai = Q0 * ((8. * A_mean * rog**(n_glen) * (1. - alpha)**(n_glen-1.)) / (4.**n_glen *f_frot))*back_force_coef**(n_glen-1.) 
    627 phi_Tsai = phi_Tsai * Hgl**(n_glen+2.) 
    628  
    629 end subroutine flux_Tsai4Schoof 
     513    debug_3D(:,:,66) = phi_prescr_tabx(:,:) 
     514    debug_3D(:,:,67) = phi_prescr_taby(:,:) 
     515 
     516  end subroutine interpol_glflux 
     517 
     518  !------------------------------------------------------------------------- 
     519  ! calcule la position sous maille de la ligne dechouage 
     520  subroutine calc_xgl4schoof(dy,alpha,H_0,B_0,SL_0,H_1,B_1,SL_1,xpos,Hgl) 
     521     
     522    implicit none 
     523 
     524    integer :: i,j 
     525    real,intent(in)        ::    dy        !< longueur orientee de la maille  
     526    real,intent(in)        ::    alpha     !< coefficient de flottaison = rho/rhow 
     527    real,intent(in)        ::    H_0       !< epaisseur au point pose 
     528    real,intent(in)        ::    B_0       !< altitude socle au point pose 
     529    real,intent(in)        ::    SL_0      !< sea level au point pose 
     530    real,intent(in)        ::    H_1       !< epaisseur au point flottant 
     531    real,intent(in)        ::    B_1       !< altitude socle au point flottant 
     532    real,intent(in)        ::    SL_1      !< sea level au point flottant 
     533    real,intent(out)       ::    xpos      !< position de la ligne (en distance depuis le point pose) 
     534    real,intent(out)       ::    Hgl       !< epaisseur de glace a la grounding line 
     535 
     536 
     537 
     538    real                   ::    Cgl       ! variable de travail   
     539 
     540    !~   Cgl   = (-alpha * (H_1-H_0) - (B_1-B_0)) ! -(rho/rhow (H1-H0) + (B1-B0)) 
     541 
     542    !~   if (abs(Cgl).gt.1.e-5) then    
     543    !~      xpos    = (B_0 + alpha * H_0 - SL) / Cgl  
     544    !~   else 
     545    !~      xpos    = 1. - 1./abs(dy)   ! verifier 
     546    !~   end if 
     547 
     548    !~   if (xpos.LT.0..OR.xpos.GT.1.) then 
     549    !~          print*,'calc_xgl4schoof : xpos value error i,j=',i,j 
     550    !~          print*, 'xpos,Cgl', xpos,Cgl 
     551    !~          print*,'B_0, alpha, H_0, SL', B_0, alpha, H_0, SL 
     552    !~          stop 
     553    !~  endif 
     554 
     555    Cgl   = (B_1-B_0) + alpha * (H_1-H_0) - (SL_1-SL_0) 
     556 
     557    if (abs(Cgl).gt.1.e-5) then    
     558       xpos    = (SL_0 - (B_0 + alpha * H_0)) / Cgl  
     559    else 
     560       xpos    = 0.5  ! verifier 
     561    end if 
     562 
     563    !  print* 
     564    !  print*, 'xpos', xpos 
     565    !   print*,'B_0, alpha, H_0, SL', B_0, alpha, H_0, SL 
     566    !   print*,'B_1, H_1', B_1, H_1 
     567 
     568 
     569    if (xpos.LT.0..OR.xpos.GT.1.) then 
     570       print*,'calc_xgl4schoof : xpos value error i,j=',i,j 
     571       print*, 'xpos,Cgl', xpos,Cgl 
     572       print*,'B_0, alpha, H_0, SL_0', B_0, alpha, H_0, SL_0 
     573       print*,'archim=',B_1+H_1*alpha - SL_1 
     574       !stop 
     575       xpos = min(max(0.,xpos),1.) 
     576    endif 
     577 
     578    xpos = xpos * dy 
     579 
     580    Hgl= H_0 + (H_1 - H_0) * xpos/dy 
     581 
     582  end subroutine calc_xgl4schoof 
     583 
     584 
     585  !-------------------------------------------------------------------- 
     586  ! Schoof_flux : calcule le coefficient et le flux de Schoof  
     587  !--------------------------------------------------------------------  
     588 
     589  ! Dans un premier temps, pas de terme horizontal (regarder dans toy_recul une tentative de confinement) 
     590 
     591  subroutine flux_Schoof4Schoof (Hgl,A_mean,C_frot,alpha,n_glen,m_weert,back_force_coef,phi_schoof) 
     592 
     593    use param_phy_mod, only:rog 
     594 
     595    implicit none 
     596 
     597    real,intent(in) :: Hgl 
     598    real,intent(in) :: A_mean 
     599    real,intent(in) :: C_frot 
     600    real,intent(in) :: alpha 
     601    real,intent(in) :: n_glen 
     602    real,intent(in) :: m_weert 
     603    real,intent(in) :: back_force_coef 
     604    real,intent(out) :: phi_schoof 
     605 
     606    ! afq: in standard GRISLI (19/05/17), tob= -beta ub => m_weert = 1 and C_frot = beta 
     607 
     608    !phi_schoof = (((A_mean * rog**(n_glen+1) * (1 - alpha)**n_glen) / (4**n_glen *C_frot)) **(1./(1.+1./m_weert)))*back_force_coef**(n_glen/(1+1./m_weert)) 
     609    !phi_schoof = phi_schoof * Hgl**((n_glen+3.+1/m_weert)/(1.+1./m_weert)) 
     610    phi_schoof = (((A_mean * rog**(n_glen+1.) * (1. - alpha)**n_glen) / (4.**n_glen*C_frot)) **(1./(1.+inv_mweert)))*back_force_coef**(n_glen/(1.+inv_mweert)) 
     611    phi_schoof = phi_schoof * Hgl**((n_glen+3.+inv_mweert)/(1.+inv_mweert)) 
     612 
     613 
     614 
     615  end subroutine flux_Schoof4Schoof 
     616 
     617  !-------------------------------------------------------------------- 
     618  ! Tsai_flux : calcule le coefficient et le flux de Schoof  
     619  !--------------------------------------------------------------------  
     620 
     621  ! Dans un premier temps, pas de terme horizontal (regarder dans toy_recul une tentative de confinement) 
     622 
     623  subroutine flux_Tsai4Schoof (Hgl,A_mean,f_frot,alpha,n_glen,back_force_coef,phi_Tsai) 
     624 
     625    use param_phy_mod, only:rog 
     626 
     627    implicit none 
     628 
     629    real,intent(in) :: Hgl 
     630    real,intent(in) :: A_mean 
     631    real,intent(in) :: f_frot 
     632    real,intent(in) :: alpha 
     633    real,intent(in) :: n_glen 
     634    real,intent(in) :: back_force_coef 
     635    real,intent(out) :: phi_Tsai 
     636 
     637    real, parameter :: Q0 = 0.61 
     638 
     639 
     640    phi_Tsai = Q0 * ((8. * A_mean * rog**(n_glen) * (1. - alpha)**(n_glen-1.)) / (4.**n_glen *f_frot))*back_force_coef**(n_glen-1.) 
     641    phi_Tsai = phi_Tsai * Hgl**(n_glen+2.) 
     642 
     643  end subroutine flux_Tsai4Schoof 
    630644 
    631645end module furst_schoof_mod 
Note: See TracChangeset for help on using the changeset viewer.