! Module Furst Schoof module furst_schoof_mod use module3d_phy use deformation_mod_2lois implicit none real :: frot_coef ! afq: fixed in Tsai but = beta in Schoof (for m_weert = 1) real :: alpha_flot integer :: gr_select real,dimension(nx,ny) :: back_force_x real,dimension(nx,ny) :: back_force_y real, parameter :: m_weert = 1.0 ! afq: is this defined elsewhere in the model??? real, parameter :: inv_mweert = 1./m_weert contains subroutine init_furst_schoof namelist/furst_schoof/frot_coef,gr_select rewind(num_param) ! pour revenir au debut du fichier param_list.dat read(num_param,furst_schoof) write(num_rep_42,*)'!___________________________________________________________' write(num_rep_42,*)'! module furst_schoof' write(num_rep_42,*) write(num_rep_42,*)'! frot_coef = ',frot_coef write(num_rep_42,*)'! gr_select = ',gr_select write(num_rep_42,*)'! frot_coef : solid friction law coef 0.035 in Ritz et al 2001' write(num_rep_42,*)'! gr_select = 1 : Tsai , 2 : Schoof' write(num_rep_42,*)'!___________________________________________________________' ! frot = 0.035 ! f in the solid friction law, 0.035 in Ritz et al. 2001 ! choix Tsai (loi de Coulomb) ou Schoof (loi de Weertman) alpha_flot= ro/row back_force_x(:,:)=0.1 back_force_y(:,:)=0.1 end subroutine init_furst_schoof subroutine interpol_glflux implicit none real :: xpos,ypos real :: Hgl ! integer :: igl real :: phi_im2, phi_im1, phi_i, phi_ip1, phi_ip2, phi_ip3 real :: phi_jm2, phi_jm1, phi_j, phi_jp1, phi_jp2, phi_jp3 real,dimension(nx,ny) :: phi_xgl, phi_ygl ! flux gr line sur maille stag integer,dimension(nx,ny) :: countx, county ! how often do we modify ux/uy real :: phi_prescr real :: toutpetit = 1e-6 real :: denom, prodscal real :: bfx, bfy !debug real,dimension(nx,ny) :: xpos_tab=9999. real,dimension(nx,ny) :: ypos_tab=9999. real,dimension(nx,ny) :: Hglx_tab=9999. real,dimension(nx,ny) :: Hgly_tab=9999. real,dimension(nx,ny) :: Hmxloc, Hmyloc ! pour eviter division par 0 real,dimension(nx,ny) :: phi_prescr_tabx, phi_prescr_taby real,dimension(nx,ny) :: uxbartemp, uybartemp real,dimension(nx,ny) :: archimtab real,dimension(nx,ny) :: imx_diag_in phi_prescr_tabx=0. phi_prescr_taby=0. Hmxloc(:,:)= max(Hmx(:,:),1.) Hmyloc(:,:)=max(Hmy(:,:),1.) uxbartemp(:,:)=0. uybartemp(:,:)=0. countx(:,:)=0. county(:,:)=0. gr_line_schoof(:,:)=0 archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot -sealevel imx_diag_in(:,:) = imx_diag(:,:) ! detection des noeuds grounding line do j= 3,ny-3 do i=3,nx-3 !archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel)) then ! grounded and with ice ! print*,'schoof gr line ij',i,j ! selon x (indice i) if (archimtab(i-1,j).LT.0..and.Uxbar(i,j).LT.0..and..not.flot_marais(i-1,j)) then !if (flot(i-1,j).and.Uxbar(i,j).LT.0.) then ! grounding line between i-1,i CAS WEST (ecoulement vers grille West) call calc_xgl4schoof(-dx,alpha_flot,H(i,j),Bsoc(i,j),H(i-1,j),Bsoc(i-1,j),Sealevel,xpos,Hgl) xpos_tab(i,j)=xpos Hglx_tab(i,j)=Hgl ! afq: the back force is on the staggered grid, the GL can be either West or East to this point. if (xpos .lt. -dx/2.) then bfx = back_force_x(i,j) else bfx = back_force_x(i+1,j) endif if (gr_select.eq.1) then ! flux de Tsai call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfx,phi_prescr) else if (gr_select.eq.2) then ! flux de Schoof ! afq: the dragging coef. is on the staggered grid, the GL can be either West or East to this point. if (xpos .lt. -dx/2.) then frot_coef = betamx(i,j) else frot_coef = betamx(i+1,j) endif call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfx,phi_prescr) else print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE' stop endif phi_prescr=-phi_prescr ! Doit être negatif dans le cas West phi_prescr_tabx(i,j)=phi_prescr if (xpos .lt. -dx/2.) then ! GL a west du i staggered, o centre, x stag ! grille ! .....x......o......x......o...G..x......O......x......o......x......o ! ! i-2 i-1 i i+1 i+2 !flux ! in out G out in !afq -- if (.not.(ilemx(i,j)).and..not.(ilemx(i-1,j))) then ! extrapolation pour avoir uxbar(i-1,j) qui sera ensuite prescrit dans call rempli_L2 phi_im2 = Uxbar(i-2,j) * Hmxloc(i-2,j) phi_xgl(i-1,j) = phi_im2 + dx * (phi_prescr-phi_im2)/(2.5 * dx + xpos) uxbartemp(i-1,j) = uxbartemp(i-1,j)+ phi_xgl(i-1,j) / Hmxloc(i-1,j) ! attention division par 0 possible ? countx (i-1,j) = countx(i-1,j)+1 imx_diag (i-1,j) = 1 gr_line_schoof(i-1,j) = 1 ! extrapolation pour avoir uxbar(i,j) phi_ip1 = Uxbar(i+1,j) * Hmxloc(i+1,j) phi_xgl(i,j) = phi_ip1 + dx * (phi_prescr - phi_ip1)/(0.5 * dx - xpos) ! print*,'uxbar(i,j)=',uxbar(i,j) uxbartemp(i,j) = uxbartemp(i,j) + phi_xgl(i,j) / Hmxloc(i,j) countx (i,j) = countx(i,j)+1 imx_diag (i,j) = 1 gr_line_schoof(i,j) = 1 ! print*,'schoof uxbar(i,j)=',uxbar(i,j) ! print*,'Hgl',Hgl ! print*,'phi_prescr',phi_prescr, phi_prescr/Hgl ! print* !afq -- end if else ! GL a Est du i staggered, o centre, x stag ! grille ! .....x......o......x......o......x...G..O......x......o......x......o ! i-2 i-1 i i+1 i+2 !flux ! in in out G out in !afq -- if (.not.(ilemx(i,j)).and..not.(ilemx(i+1,j))) then ! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2 phi_im1 = Uxbar(i-1,j) * Hmxloc(i-1,j) phi_xgl(i,j) = phi_im1 + dx * (phi_prescr-phi_im1)/(1.5 * dx + xpos) uxbartemp(i,j) = uxbartemp(i,j) + phi_xgl(i,j) / Hmxloc(i,j) ! attention division par 0 possible ? countx (i,j) = countx(i,j)+1 imx_diag (i,j) = 1 gr_line_schoof(i,j) = 1 ! extrapolation pour avoir uxbar(i+1,j) phi_ip2 = Uxbar(i+2,j) * Hmxloc(i+2,j) phi_xgl(i+1,j) = phi_ip2 + dx * (phi_prescr - phi_ip2)/(1.5 * dx - xpos) uxbartemp(i+1,j) = uxbartemp(i+1,j) + phi_xgl(i+1,j) / Hmxloc(i+1,j) countx (i+1,j) = countx(i+1,j)+1 imx_diag (i+1,j) = 1 gr_line_schoof(i+1,j) = 1 !afq -- end if end if else if (archimtab(i+1,j).LT.0..and.Uxbar(i+1,j).GT.0..and..not.flot_marais(i+1,j)) then !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) call calc_xgl4schoof(dx,alpha_flot,H(i,j),Bsoc(i,j),H(i+1,j),Bsoc(i+1,j),Sealevel,xpos,Hgl) xpos_tab(i,j)=xpos Hglx_tab(i,j)=Hgl ! afq: the back force is on the staggered grid, the GL can be either West or East to this point. if (xpos .lt. dx/2.) then bfx = back_force_x(i,j) else bfx = back_force_x(i+1,j) endif if (gr_select.eq.1) then ! flux de Tsai call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfx,phi_prescr) else if (gr_select.eq.2) then ! flux de Schoof ! afq: the dragging coef. is on the staggered grid, the GL can be either West or East to this point. if (xpos .lt. dx/2.) then frot_coef = betamx(i,j) else frot_coef = betamx(i+1,j) endif call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfx,phi_prescr) else print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE' stop endif phi_prescr_tabx(i,j)=phi_prescr if (xpos .lt. dx/2.) then ! GL a west du i staggered, o centre, x stag ! grille ! .....x......o......x......o......x......O..G...x......o......x......o ! ! i-2 i-1 i i+1 i+2 !flux ! in out G out in !afq -- if (.not.(ilemx(i,j)).and..not.(ilemx(i+1,j))) then ! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2 phi_im1 = Uxbar(i-1,j) * Hmxloc(i-1,j) phi_xgl(i,j) = phi_im1 + dx * (phi_prescr-phi_im1)/(1.5 * dx + xpos) uxbartemp(i,j) = uxbartemp(i,j) + phi_xgl(i,j) / Hmxloc(i,j) ! attention division par 0 possible ? countx (i,j) = countx(i,j)+1 imx_diag (i,j) = 1 gr_line_schoof(i,j) = 1 ! extrapolation pour avoir uxbar(i+1,j) phi_ip2 = Uxbar(i+2,j) * Hmxloc(i+2,j) phi_xgl(i+1,j) = phi_ip2 + dx * (phi_prescr - phi_ip2)/(1.5* dx - xpos) uxbartemp(i+1,j) = uxbartemp(i+1,j) + phi_xgl(i+1,j) / Hmxloc(i+1,j) countx (i+1,j) = countx(i+1,j)+1 imx_diag (i+1,j) = 1 gr_line_schoof(i+1,j) = 1 !afq -- end if else ! GL a west du i staggered, o centre, x stag ! grille ! ......x......o......x......O......x...G..o......x......o......x.....o ! ! i-1 i i+1 i+2 i+3 !flux ! in out G out in !afq -- if (.not.(ilemx(i+1,j)).and..not.(ilemx(i+2,j))) then ! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2 phi_i = Uxbar(i,j) * Hmxloc(i,j) phi_xgl(i+1,j) = phi_i + dx * (phi_prescr-phi_i)/(0.5 * dx + xpos) uxbartemp(i+1,j) = uxbartemp(i+1,j) + phi_xgl(i+1,j) / Hmxloc(i+1,j) ! attention division par 0 possible ? countx (i+1,j) = countx(i+1,j)+1 imx_diag (i+1,j) = 1 gr_line_schoof(i+1,j) = 1 ! extrapolation pour avoir uxbar(i+1,j) phi_ip3 = Uxbar(i+3,j) * Hmxloc(i+3,j) phi_xgl(i+2,j) = phi_ip3 + dx * (phi_prescr - phi_ip3)/(2.5 * dx - xpos) uxbartemp(i+2,j) = uxbartemp(i+2,j) + phi_xgl(i+2,j) / Hmxloc(i+2,j) countx (i+2,j) = countx(i+2,j)+1 imx_diag (i+2,j) = 1 gr_line_schoof(i+2,j) = 1 !afq -- end if end if end if !******************************************************************************* ! selon y (indice j) if (archimtab(i,j-1).LT.0..and.Uybar(i,j).LT.0..and..not.flot_marais(i,j-1)) then !if (flot(i,j-1).and.Uybar(i,j).LT.0.) then ! grounding line between j-1,j CAS SUD (ecoulement vers grille SUD) call calc_xgl4schoof(-dy,alpha_flot,H(i,j),Bsoc(i,j),H(i,j-1),Bsoc(i,j-1),Sealevel,ypos,Hgl) ypos_tab(i,j)=ypos Hgly_tab(i,j)=Hgl ! afq: the back force is on the staggered grid, the GL can be either South or North to this point. if (ypos .lt. -dy/2.) then bfy = back_force_y(i,j) else bfy = back_force_y(i,j+1) endif if (gr_select.eq.1) then ! flux de Tsai call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfy,phi_prescr) else if (gr_select.eq.2) then ! flux de Schoof ! afq: the dragging coef. is on the staggered grid, the GL can be either South or North to this point. if (ypos .lt. -dy/2.) then frot_coef = betamy(i,j) else frot_coef = betamy(i,j+1) endif call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfy,phi_prescr) else print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE' stop endif phi_prescr=-phi_prescr ! Doit être negatif dans le cas Sud phi_prescr_taby(i,j)=phi_prescr if (ypos .lt. -dy/2.) then ! GL au sud du j staggered, o centre, x stag ! grille ! .....x......o......x......o...G..x......O......x......o......x......o ! ! j-2 j-1 j j+1 j+2 !flux ! in out G out in !afq -- if (.not.(ilemy(i,j)).and..not.(ilemy(i,j-1))) then ! extrapolation pour avoir uybar(i,j-1) qui sera ensuite prescrit dans call rempli_L2 phi_jm2 = Uybar(i,j-2) * Hmyloc(i,j-2) phi_ygl(i,j-1) = phi_jm2 + dy * (phi_prescr-phi_jm2)/(2.5 * dy + ypos) uybartemp(i,j-1) = uybartemp(i,j-1) + phi_ygl(i,j-1) / Hmyloc(i,j-1) ! attention division par 0 possible ? county (i,j-1) = county(i,j-1)+1 imy_diag (i,j-1) = 1 gr_line_schoof(i,j-1) = 1 ! extrapolation pour avoir uybar(i,j) phi_jp1 = Uybar(i,j+1) * Hmyloc(i,j+1) phi_ygl(i,j) = phi_jp1 + dy * (phi_prescr - phi_jp1)/(0.5 * dy - ypos) uybartemp(i,j) = uybartemp(i,j) + phi_ygl(i,j) / Hmyloc(i,j) county (i,j) = county(i,j)+1 imy_diag (i,j) = 1 gr_line_schoof(i,j) = 1 !afq -- end if else ! GL au nord du j staggered, o centre, x stag ! grille ! .....x......o......x......o......x...G..O......x......o......x......o ! j-2 j-1 j j+1 j+2 !flux ! in in out G out in !afq -- if (.not.(ilemy(i,j)).and..not.(ilemy(i,j+1))) then ! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2 phi_jm1 = Uybar(i,j-1) * Hmyloc(i,j-1) phi_ygl(i,j) = phi_jm1 + dy * (phi_prescr-phi_jm1)/(1.5 * dy + ypos) uybartemp(i,j) = uybartemp(i,j) + phi_ygl(i,j) / Hmyloc(i,j) ! attention division par 0 possible ? county (i,j) = county(i,j)+1 imy_diag(i,j) = 1 gr_line_schoof(i,j) = 1 ! extrapolation pour avoir uybar(i,j+1) phi_jp2 = Uybar(i,j+2) * Hmyloc(i,j+2) phi_ygl(i,j+1) = phi_jp2 + dy * (phi_prescr - phi_jp2)/(1.5 * dy - ypos) uybartemp(i,j+1) = uybartemp(i,j+1) + phi_ygl(i,j+1) / Hmyloc(i,j+1) county (i,j+1) = county(i,j+1)+1 imy_diag (i,j+1) = 1 gr_line_schoof(i,j+1) = 1 !afq -- end if end if else if (archimtab(i,j+1).LT.0..and.Uybar(i,j+1).GT.0..and..not.flot_marais(i,j+1)) then !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) call calc_xgl4schoof(dy,alpha_flot,H(i,j),Bsoc(i,j),H(i,j+1),Bsoc(i,j+1),Sealevel,ypos,Hgl) ypos_tab(i,j)=ypos Hgly_tab(i,j)=Hgl ! afq: the back force is on the staggered grid, the GL can be either South or North to this point. if (ypos .lt. dy/2.) then bfy = back_force_y(i,j) else bfy = back_force_y(i,j+1) endif if (gr_select.eq.1) then ! flux de Tsai call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfy,phi_prescr) else if (gr_select.eq.2) then ! flux de Schoof ! afq: the dragging coef. is on the staggered grid, the GL can be either South or North to this point. if (ypos .lt. dy/2.) then frot_coef = betamy(i,j) else frot_coef = betamy(i,j+1) endif call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfy,phi_prescr) else print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE' stop endif phi_prescr_taby(i,j)=phi_prescr if (ypos .lt. dy/2.) then ! GL au sud du j staggered, o centre, x stag ! grille ! .....x......o......x......o......x......O..G...x......o......x......o ! ! j-2 j-1 j j+1 j+2 !flux ! in out G out in !afq -- if (.not.(ilemy(i,j)).and..not.(ilemy(i,j+1))) then ! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2 phi_jm1 = Uybar(i,j-1) * Hmyloc(i,j-1) phi_ygl(i,j) = phi_jm1 + dy * (phi_prescr-phi_jm1)/(1.5 * dy + ypos) uybartemp(i,j) = uybartemp(i,j) + phi_ygl(i,j) / Hmyloc(i,j) ! attention division par 0 possible ? county (i,j) = county(i,j)+1 imy_diag (i,j) = 1 gr_line_schoof(i,j) = 1 ! extrapolation pour avoir uybar(i,j+1) phi_jp2 = Uybar(i,j+2) * Hmyloc(i,j+2) phi_ygl(i,j+1) = phi_jp2 + dy * (phi_prescr - phi_jp2)/(1.5* dy - ypos) uybartemp(i,j+1) = uybartemp(i,j+1) + phi_ygl(i,j+1) / Hmyloc(i,j+1) county (i,j+1) = county(i,j+1)+1 imy_diag (i,j+1) = 1 gr_line_schoof(i,j+1) = 1 !afq -- end if else ! GL au nord du j staggered, o centre, x stag ! grille ! ......x......o......x......O......x...G..o......x......o......x.....o ! ! j-1 j j+1 j+2 j+3 !flux ! in out G out in !afq -- if (.not.(ilemy(i,j+1)).and..not.(ilemy(i,j+2))) then ! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2 phi_j = Uybar(i,j) * Hmyloc(i,j) phi_ygl(i,j+1) = phi_j + dy * (phi_prescr-phi_j)/(0.5 * dy + ypos) uybartemp(i,j+1) = uybartemp(i,j+1) + phi_ygl(i,j+1) / Hmyloc(i,j+1) ! attention division par 0 possible ? county (i,j+1) = county(i,j+1)+1 imy_diag (i,j+1) = 1 gr_line_schoof(i,j+1) = 1 ! extrapolation pour avoir uybar(i,j+1) phi_jp3 = Uybar(i,j+3) * Hmyloc(i,j+3) phi_ygl(i,j+2) = phi_jp3 + dy * (phi_prescr - phi_jp3)/(2.5 * dy - ypos) uybartemp(i,j+2) = uybartemp(i,j+2) + phi_ygl(i,j+2) / Hmyloc(i,j+2) county (i,j+2) = county(i,j+2)+1 imy_diag (i,j+2) = 1 gr_line_schoof(i,j+2) = 1 !afq -- end if end if end if end if end do end do ! afq -- if we discard the points with multiple fluxes coming, uncomment: where (countx(:,:).ge.2) uxbartemp(:,:)=uxbar(:,:)*countx(:,:) imx_diag(:,:) = imx_diag_in(:,:) end where where (county(:,:).ge.2) uybartemp(:,:)=uybar(:,:)*county(:,:) imy_diag(:,:) = imx_diag_in(:,:) end where ! afq -- where (countx(:,:).ne.0) uxbar(:,:)= max(min(uxbartemp(:,:)/countx(:,:),V_limit),-V_limit) where (county(:,:).ne.0) uybar(:,:)= max(min(uybartemp(:,:)/county(:,:),V_limit),-V_limit) do j= 3,ny-3 do i=3,nx-3 denom = sqrt((uxbartemp(i,j)**2 + uybartemp(i,j)**2))*sqrt((uxbar(i,j)**2 + uybar(i,j)**2)) if (denom.gt.toutpetit) then ! prodscal is cos theta in : U . V = u v cos theta ! this number is between -1 and 1, 1 = same direction & -1 = opposite direction ! we could use this as an artificial back force... prodscal = (uxbartemp(i,j)*uxbar(i,j) + uybartemp(i,j)*uybar(i,j))/ & (sqrt((uxbartemp(i,j)**2 + uybartemp(i,j)**2))*sqrt((uxbar(i,j)**2 + uybar(i,j)**2))) endif end do end do !!$ do j=1,ny !!$ do i=1,nx !!$ write(598,*) xpos_tab(i,j) !!$ write(599,*) ypos_tab(i,j) !!$ write(600,*) Hglx_tab(i,j) !!$ write(601,*) Hgly_tab(i,j) !!$ write(602,*) Abar(i,j) !!$ write(603,*) phi_prescr_tabx(i,j)/Hglx_tab(i,j) !!$ write(604,*) phi_prescr_taby(i,j)/Hgly_tab(i,j) !!$ enddo !!$ enddo debug_3D(:,:,66) = phi_prescr_tabx(:,:) debug_3D(:,:,67) = phi_prescr_taby(:,:) end subroutine interpol_glflux !------------------------------------------------------------------------- ! calcule la position sous maille de la ligne dechouage subroutine calc_xgl4schoof(dy,alpha,H_0,B_0,H_1,B_1,SL,xpos,Hgl) implicit none real,intent(in) :: dy !< longueur orientee de la maille real,intent(in) :: alpha !< coefficient de flottaison = rho/rhow real,intent(in) :: H_0 !< epaisseur au point pose real,intent(in) :: B_0 !< altitude socle au point pose real,intent(in) :: H_1 !< epaisseur au point flottant real,intent(in) :: B_1 !< altitude socle au point flottant real,intent(in) :: SL !< Sea level real,intent(out) :: xpos !< position de la ligne (en distance depuis le point pose) real,intent(out) :: Hgl !< Epaisseur de glace a la grounding line real :: Cgl ! variable de travail !~ Cgl = (-alpha * (H_1-H_0) - (B_1-B_0)) ! -(rho/rhow (H1-H0) + (B1-B0)) !~ if (abs(Cgl).gt.1.e-5) then !~ xpos = (B_0 + alpha * H_0 - SL) / Cgl !~ else !~ xpos = 1. - 1./abs(dy) ! verifier !~ end if !~ if (xpos.LT.0..OR.xpos.GT.1.) then !~ print*,'calc_xgl4schoof : xpos value error i,j=',i,j !~ print*, 'xpos,Cgl', xpos,Cgl !~ print*,'B_0, alpha, H_0, SL', B_0, alpha, H_0, SL !~ stop !~ endif Cgl = (B_1-B_0) + alpha * (H_1-H_0) if (abs(Cgl).gt.1.e-5) then xpos = (SL - (B_0 + alpha * H_0)) / Cgl else xpos = 0.5 ! verifier end if ! print* ! print*, 'xpos', xpos ! print*,'B_0, alpha, H_0, SL', B_0, alpha, H_0, SL ! print*,'B_1, H_1', B_1, H_1 if (xpos.LT.0..OR.xpos.GT.1.) then print*,'calc_xgl4schoof : xpos value error i,j=',i,j print*, 'xpos,Cgl', xpos,Cgl print*,'B_0, alpha, H_0, SL', B_0, alpha, H_0, SL print*,'archim=',B_1+H_1*alpha - SL !stop xpos = min(max(0.,xpos),1.) endif xpos = xpos * dy Hgl= H_0 + (H_1 - H_0) * xpos/dy end subroutine calc_xgl4schoof !-------------------------------------------------------------------- ! Schoof_flux : calcule le coefficient et le flux de Schoof !-------------------------------------------------------------------- ! Dans un premier temps, pas de terme horizontal (regarder dans toy_recul une tentative de confinement) subroutine flux_Schoof4Schoof (Hgl,A_mean,C_frot,alpha,n_glen,m_weert,back_force_coef,phi_schoof) implicit none real,intent(in) :: Hgl real,intent(in) :: A_mean real,intent(in) :: C_frot real,intent(in) :: alpha real,intent(in) :: n_glen real,intent(in) :: m_weert real,intent(in) :: back_force_coef real,intent(out) :: phi_schoof ! afq: in standard GRISLI (19/05/17), tob= -beta ub => m_weert = 1 and C_frot = beta !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)) !phi_schoof = phi_schoof * Hgl**((n_glen+3.+1/m_weert)/(1.+1./m_weert)) 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)) phi_schoof = phi_schoof * Hgl**((n_glen+3.+inv_mweert)/(1.+inv_mweert)) end subroutine flux_Schoof4Schoof !-------------------------------------------------------------------- ! Tsai_flux : calcule le coefficient et le flux de Schoof !-------------------------------------------------------------------- ! Dans un premier temps, pas de terme horizontal (regarder dans toy_recul une tentative de confinement) subroutine flux_Tsai4Schoof (Hgl,A_mean,f_frot,alpha,n_glen,back_force_coef,phi_Tsai) implicit none real,intent(in) :: Hgl real,intent(in) :: A_mean real,intent(in) :: f_frot real,intent(in) :: alpha real,intent(in) :: n_glen real,intent(in) :: back_force_coef real,intent(out) :: phi_Tsai real, parameter :: Q0 = 0.61 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.) phi_Tsai = phi_Tsai * Hgl**(n_glen+2.) end subroutine flux_Tsai4Schoof end module furst_schoof_mod