Changeset 104 for trunk/SOURCES/furst_schoof_mod.f90

Timestamp:

05/09/17 10:20:19 (7 years ago)

Author:

aquiquet

Message:

Schoof: back force in x/y from ramollos.

File:

• trunk/SOURCES/furst_schoof_mod.f90 (modified) (10 diffs)

trunk/SOURCES/furst_schoof_mod.f90

@@ -11 +11 @@
 real :: alpha_flot
 integer :: gr_select
-real,dimension(nx,ny) :: back_force
+real,dimension(nx,ny) :: back_force_x
+real,dimension(nx,ny) :: back_force_y
 
 contains
@@ -35 +36 @@
   ! choix Tsai (loi de Coulomb) ou Schoof (loi de Weertman)
   alpha_flot= ro/row
-  back_force(:,:)=0.1
+  back_force_x(:,:)=0.1
+  back_force_y(:,:)=0.1
   
 end subroutine init_furst_schoof  
@@ -68 +70 @@
   
   real,dimension(nx,ny) :: archimtab
+  
+  real,dimension(nx,ny) :: imx_diag_in
 
   phi_prescr_tabx=0.
@@ -82 +86 @@
   archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot -sealevel
 
+  imx_diag_in(:,:) = imx_diag(:,:)
+  
   ! detection des noeuds grounding line
   do j= 3,ny-3
@@ -98 +104 @@
 
               if (gr_select.eq.1) then ! flux de Tsai
-                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force(i-1,j),phi_prescr)
+                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force_x(i-1,j),phi_prescr)
               else
                  print*,'ATTENTION FLUX AUTRE QUE TSAI NON IMPLEMENTE'
@@ -167 +173 @@
               Hglx_tab(i,j)=Hgl
               if (gr_select.eq.1) then ! flux de Tsai
-                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force(i+1,j),phi_prescr)
+                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force_x(i+1,j),phi_prescr)
               else
                  print*,'ATTENTION FLUX AUTRE QUE TSAI NON IMPLEMENTE'
@@ -232 +238 @@
               Hgly_tab(i,j)=Hgl
               if (gr_select.eq.1) then ! flux de Tsai
-                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force(i,j-1),phi_prescr)
+                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force_y(i,j-1),phi_prescr)
               else
                  print*,'ATTENTION FLUX AUTRE QUE TSAI NON IMPLEMENTE'
@@ -293 +299 @@
               Hgly_tab(i,j)=Hgl
               if (gr_select.eq.1) then ! flux de Tsai
-                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force(i,j+1),phi_prescr)
+                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force_y(i,j+1),phi_prescr)
               else
                  print*,'ATTENTION FLUX AUTRE QUE TSAI NON IMPLEMENTE'
@@ -350 +356 @@
      end do
   end do
-  where (countx(:,:).ne.0) uxbar(:,:)=uxbartemp(:,:)/countx(:,:)
-  where (county(:,:).ne.0) uybar(:,:)=uybartemp(:,:)/county(:,:)
-  
-  uxbar(:,:)= max(min(uxbar(:,:),V_limit),-V_limit)
-  uybar(:,:)= max(min(uybar(:,:),V_limit),-V_limit)
-  
+
+  ! afq -- if we discard the points with multiple fluxes coming, uncomment:
+  !where (countx(:,:).ge.2)
+  !   uxbartemp(:,:)=uxbar(i,j)*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(:,:),1000.),-1000.)
+  where (county(:,:).ne.0) uybar(:,:)= max(min(uybartemp(:,:)/county(:,:),1000.),-1000.)
+           
   do j= 3,ny-3
      do i=3,nx-3
@@ -433 +448 @@
         
   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.)
+     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