Changeset 104

Timestamp:

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

Author:

aquiquet

Message:

Schoof: back force in x/y from ramollos.

Location:

trunk/SOURCES

Files:

• New-remplimat/diagno-L2_mod.f90 (modified) (5 diffs)
• New-remplimat/remplimat-shelves-tabTu.f90 (modified) (2 diffs)
• furst_schoof_mod.f90 (modified) (10 diffs)

trunk/SOURCES/New-remplimat/diagno-L2_mod.f90

@@ -141 +141 @@
 !~      enddo   
   
-  if (Schoof.eq.1.and.nt.GT.15000) then ! flux grounding line Schoof
-                call interpol_glflux ! calcul flux GL + interpolation sur voisins
-        endif   
-        
+  !if (Schoof.eq.1.and.nt.GT.15000) then ! flux grounding line Schoof
+  if (Schoof.eq.1) then ! flux grounding line Schoof
+     call interpol_glflux ! calcul flux GL + interpolation sur voisins
+  endif
+
 !~       do j=1,ny
 !~              do i=1,nx
@@ -153 +154 @@
 !~      print*,'ecriteure termineee !!!!!!'
 !~      read(*,*)
-        
+
   ! donnent les cas de conditions aux limites
   !
@@ -185 +186 @@
        imx_diag(nxd1:nxd2,nyd1:nyd2),imy_diag(nxd1:nxd2,nyd1:nyd2),ifail_diagno)
        
-  if (Schoof.eq.1.and.nt.GT.15000) then ! flux grounding line Schoof avec calcul de la back force par shelf ramollo
+  !if (Schoof.eq.1.and.nt.GT.15000) then ! flux grounding line Schoof avec calcul de la back force par shelf ramollo
+  if (Schoof.eq.1) then ! flux grounding line Schoof avec calcul de la back force par shelf ramollo
     pvi_keep(:,:)=pvi(:,:)       
-                where (flot(:,:).and.H(:,:).GT.2.)
-                                pvi(:,:)=1.e5    
-!                               pvi(:,:)=pvimin
-                endwhere
-                                        
-                call rempli_L2(nxd1,nxd2,nyd1,nyd2,uxbar(nxd1:nxd2,nyd1:nyd2),uybar(nxd1:nxd2,nyd1:nyd2), & 
+    where (flot(:,:).and.H(:,:).GT.2.)
+       pvi(:,:)=1.e5    
+!       pvi(:,:)=pvimin
+    endwhere
+
+    call rempli_L2(nxd1,nxd2,nyd1,nyd2,uxbar(nxd1:nxd2,nyd1:nyd2),uybar(nxd1:nxd2,nyd1:nyd2), & 
        uxb1ramollo(nxd1:nxd2,nyd1:nyd2),uyb1ramollo(nxd1:nxd2,nyd1:nyd2),  &
        imx_diag(nxd1:nxd2,nyd1:nyd2),imy_diag(nxd1:nxd2,nyd1:nyd2),ifail_diagno_ramollo)
@@ -198 +200 @@
     pvi(:,:)=pvi_keep(:,:) 
     
-    where ((uxb1ramollo(:,:)**2 + uyb1ramollo(:,:)**2).GT.0.)
-                        back_force(:,:) = sqrt(uxb1(:,:)**2+ uyb1(:,:)**2) / sqrt(uxb1ramollo(:,:)**2 + uyb1ramollo(:,:)**2)
-                elsewhere
-                        back_force(:,:)=1.
-                endwhere
-                
-                if (ifail_diagno_ramollo.gt.0) then
-                        write(6,*) ' Probleme resolution systeme L2. ramollo ifail=',ifail_diagno_ramollo
-                        STOP
-                endif
+    where (abs(uxb1ramollo(:,:)) .GT.1.e-5)
+       back_force_x(:,:) = 1.0 * abs(uxb1(:,:)) / abs(uxb1ramollo(:,:))
+    elsewhere
+       back_force_x(:,:)=1.
+    endwhere
+    where (abs(uyb1ramollo(:,:)) .GT.1.e-5)
+       back_force_y(:,:) = 1.0 * abs(uyb1(:,:)) / abs(uyb1ramollo(:,:))
+    elsewhere
+       back_force_y(:,:)=1.
+    endwhere
+
+    if (ifail_diagno_ramollo.gt.0) then
+!       write(6,*) ' Probleme resolution systeme L2. ramollo ifail=',ifail_diagno_ramollo
+!       STOP
+       write(*,*) ' Probleme resolution systeme L2. ramollo ifail=',ifail_diagno_ramollo
+       write(*,*) '          ... we go on anyway!'
+    endif
 !~   do j=1,ny
 !~              do i=1,nx
@@ -217 +226 @@
 !~              enddo           
 !~      enddo
-        
+
 !~   print*,'apres calcul rempli_L2'
 !~   read(*,*)      

trunk/SOURCES/New-remplimat/remplimat-shelves-tabTu.f90

@@ -61 +61 @@
 integer :: ifail_L2         ! pour les rapports d'erreur
   
-real, dimension(nx1:nx2,ny1:ny2),intent(inout)  :: uxprec    ! vitesse en entree routine
-real, dimension(nx1:nx2,ny1:ny2),intent(inout)  :: uyprec    ! vitesse en entree routine
+real, dimension(nx1:nx2,ny1:ny2),intent(in)  :: uxprec    ! vitesse en entree routine
+real, dimension(nx1:nx2,ny1:ny2),intent(in)  :: uyprec    ! vitesse en entree routine
 
 ! masques vitesses.
@@ -76 +76 @@
   
 
-integer, dimension(nx1:nx2,ny1:ny2),intent(inout)  :: imx    ! masque en entree routine
-integer, dimension(nx1:nx2,ny1:ny2),intent(inout)  :: imy    ! masque en entree routine
+integer, dimension(nx1:nx2,ny1:ny2),intent(in)  :: imx    ! masque en entree routine
+integer, dimension(nx1:nx2,ny1:ny2),intent(in)  :: imy    ! masque en entree routine
 
 

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