source: branches/GRISLIv3/SOURCES/New-remplimat/diagno-L2_mod.f90 @ 444

!module diagno_L2_mod               ! Nouvelle version, compatible remplimat 2008 Cat
module diagno_mod                   ! nom pendant les tests
 !$ USE OMP_LIB
use module3D_phy,only: nx,ny,dx
     
implicit none



real                   :: somint,test,delp,prec
real, dimension(nx,ny) :: uxb1
real, dimension(nx,ny) :: uyb1

real, dimension(nx,ny) :: uxb1ramollo
real, dimension(nx,ny) :: uyb1ramollo
real, dimension(nx,ny) :: pvi_keep

!cdc transfere dans module3d pour compatibilite avec furst_schoof_mod
!cdc  integer, dimension(nx,ny) :: imx_diag
!cdc  integer, dimension(nx,ny) :: imy_diag

integer :: nneigh=max(int(400000./dx),1) !nb. points around the grline kept for back force comp. in case we reduce the extent
integer, dimension(nx,ny) :: imx_reduce !afq -- to limit the number of points where vel is computed for backforce
integer, dimension(nx,ny) :: imy_reduce !afq -- to limit the number of points where vel is computed for backforce

integer :: nxd1,nxd2     ! domaine selon x Dans l'appel rempli_L2
integer :: nyd1,nyd2     ! domaine selon y

integer :: itour_pvi

integer :: ifail_diagno  ! pour recuperation d'erreur
integer :: ifail_diagno_ramollo  ! pour recuperation d'erreur shelf ramollo
integer :: iplus1,jplus1
integer ::  ctvisco,iumax,jumax
real  :: delumax,errmax
real  :: phiphi,bt2,d02,discr,ttau
real :: sf3,sf1,epsxxm,epsyym,epsm,sf01,sf03    ! pour le calcul de la viscosite
real :: viscm
real :: sf_shelf                                ! coef mult enhancement factor pour shelves

logical :: stopvisco,viscolin
logical :: test_visc

integer :: itour_nolin

contains    

!------------------------------------------------------------------------------------
subroutine init_diagno

use module3D_phy,only: num_rep_42,num_param,geoplace,pvimin
use deformation_mod_2lois,only: sf ! afq - 07/2023 - not clean since it is in module_choix

implicit none

integer sffunc

namelist/diagno_rheol/sf01,sf03,sffunc,pvimin

! attribution des coefficients de viscosite

! formats pour les ecritures dans 42
428 format(A)

! lecture des parametres du run                      block draghwat
!--------------------------------------------------------------------

rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
read(num_param,diagno_rheol)

write(num_rep_42,428)'!___________________________________________________________' 
write(num_rep_42,428) '&diagno_rheol              ! nom du bloc  diagno_rheol'
write(num_rep_42,*)
write(num_rep_42,*) 'sf01           = ',sf01 
write(num_rep_42,*) 'sf03           = ',sf03
write(num_rep_42,*) 'sffunc         = ',sffunc
write(num_rep_42,*) 'pvimin         = ',pvimin
write(num_rep_42,*)'/'     
write(num_rep_42,428) '! coefficients par rapport a la loi glace posee '                       
write(num_rep_42,428) '! sf01 : coefficient viscosite loi lineaire '
write(num_rep_42,428) '! sf03 : coefficient viscosite loi n=3 '
write(num_rep_42,428) '! sffunc : 1 for visco coef that depends on enh factor'
write(num_rep_42,428) '! pvimin : valeur de pvi pour les noeuds fictifs ~ 1.e3'
write(num_rep_42,428) '! tres petit par rapport aux valeurs standards ~ 1.e10'

write(num_rep_42,*)

if (sffunc.eq.1) then
     sf01 = (1 - (sf(2)-1)/10.)/sf(2)  !the sf are known via the module_choix...
     sf03 = (1 - (sf(1)-1)/10.)/sf(1)
endif

!      Precision utilisee dans de calcul
prec = 1.e-2
itour_pvi=1   ! si prend les valeurs analytiques dans le shelf

if (geoplace(1:5).eq.'mism3') then 
   sf_shelf = 1.
   itour_pvi= 0

else
   sf_shelf = 0.4
end if

return
end subroutine init_diagno

!------------------------------------------------------------------------------------
subroutine diagnoshelf !      Resolution numerique des equations diagnostiques

  use module3D_phy, only: itracebug,niter_nolin,taushelf,flot,h,pvi,coef_sflot,rog, & 
                          gr_line_schoof,geoplace,uxbar,uybar,pvimin,debug_3D,V_limit, & 
                          imx_diag,imy_diag,tobmx,tobmy,neffmx,neffmy,gr_line, schoof, & 
                          betamx,betamy,flgzmx,flgzmy,uxflgz,uyflgz,taub 
  use furst_schoof_mod, only: back_force_x,back_force_y,gr_select,frot_coef,interpol_glflux 
  use module_choix, only: btt,dragging,mstream_dragging  ! afq - 03/2023 - not clean 

  implicit none 

  integer :: diagno_grline
  integer itour_nolin
  integer i,j ! loop integers

  
  if (itracebug.eq.1)  call tracebug(' Entree dans diagnoshelf')

  do itour_nolin=1,niter_nolin

  itour_pvi=itour_pvi+1       ! boucle sur la viscosite (pour l'instant pas actif)

!  pvi(:,:)=0.
  Taushelf(:,:)=0.

  ! attention le bloc suivant est pour debug
!!$gzmx(:,:)=.false.
!!$gzmy(:,:)=.false.
!!$ilemx(:,:)=.false.
!!$ilemy(:,:)=.false.
!!$flgzmx(:,:)=flotmx(:,:)
!!$flgzmy(:,:)=flotmy(:,:)

  call dragging                    ! doit etre appele avant imx_imy
  call mstream_dragging ! afq -- mise a jour des masques


if (itour_pvi.le.1) then
  call calc_pvi                    ! calcule les viscosites integrees

!$OMP PARALLEL
!$OMP WORKSHARE
   where (flot(:,:).and.(H.gt.1))  ! valeur analytique pour les shelfs
      pvi(:,:) = (4./coef_Sflot/rog)**2/btt(:,:,1,1)/H(:,:)
   end where
!$OMP END WORKSHARE
!$OMP END PARALLEL
  ! avec couplage thermomecanique
!  write(166,*) ' apres call calc_pvi',itour_pvi

else
   call calc_pvi    
end if

  call imx_imy_nx_ny         ! pour rempli_L2 : calcule les masques imx et imy qui 

  diagno_grline=sum(gr_line(2:(nx-1),2:(ny-1)))
  if ((Schoof.ge.1).and.(diagno_grline.gt.0)) then
       call imx_imy_nx_ny_reduce(1) !afq: provides imx_reduce and imy_reduce
  endif

!cdc debug Schoof !!!!!!!!!!!!  
!~   do j=1,ny
!~              do i=1,nx
!~                      write(578,*) uxbar(i,j)
!~                      write(579,*) uybar(i,j)
!~              enddo
!~      enddo   
  
  !if (Schoof.eq.1.and.nt.GT.15000) then ! flux grounding line Schoof
! afq -- below:  if (Schoof.eq.1) then ! flux grounding line Schoof
! afq -- below:     call interpol_glflux ! calcul flux GL + interpolation sur voisins
! afq -- below:  endif

!~       do j=1,ny
!~              do i=1,nx
!~                      write(588,*) uxbar(i,j)
!~                      write(589,*) uybar(i,j)
!~              enddo
!~      enddo   
!~      print*,'ecriteure termineee !!!!!!'
!~      read(*,*)

  ! donnent les cas de conditions aux limites
  !
  ! version pour travailler sur tout le domaine nx ny

  if (geoplace(1:5).eq.'mism3') call mismip_boundary_cond


  ! appel a la routine rempl_L2 -------------------domaine nx x ny ------------
  ! 

  ! pour tout le domaine
  nxd1=1
  nxd2=nx
  nyd1=1
  nyd2=ny

  !call rempli_L2(1,nx,1,ny,uxbar,uybar,uxb1,uyb1,imx_diag,imy_diag,ifail_diagno)    
  !nxd1=15
  !nxd2=19
  !nyd1=30
  !nyd2=34

  !nxd1=35
  !nxd2=60
  !nyd1=35
  !nyd2=60

       
  !if (Schoof.eq.1.and.nt.GT.15000) then ! flux grounding line Schoof avec calcul de la back force par shelf ramollo
  if ((Schoof.ge.1).and.(diagno_grline.gt.0)) then ! flux grounding line Schoof avec calcul de la back force par shelf ramollo

    call rempli_L2(nxd1,nxd2,nyd1,nyd2,uxbar(nxd1:nxd2,nyd1:nyd2),uybar(nxd1:nxd2,nyd1:nyd2), & 
          uxb1(nxd1:nxd2,nyd1:nyd2),uyb1(nxd1:nxd2,nyd1:nyd2),  &
          imx_reduce(nxd1:nxd2,nyd1:nyd2),imy_reduce(nxd1:nxd2,nyd1:nyd2),ifail_diagno)

    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), & 
       uxb1ramollo(nxd1:nxd2,nyd1:nyd2),uyb1ramollo(nxd1:nxd2,nyd1:nyd2),  &
       imx_reduce(nxd1:nxd2,nyd1:nyd2),imy_reduce(nxd1:nxd2,nyd1:nyd2),ifail_diagno_ramollo)
       
    pvi(:,:)=pvi_keep(:,:) 
    
    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
    back_force_x(:,:) = min ( back_force_x(:,:), 1. )
    back_force_y(:,:) = min ( back_force_y(:,:), 1. )
    debug_3D(:,:,64) = back_force_x(:,:)
    debug_3D(:,:,65) = back_force_y(:,:)

    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
!~                if (sqrt(uxb1(i,j)**2+ uyb1(i,j)**2).gt.0..and..not.flot(i,j)) then
!~                              write(1034,*) sqrt(uxb1(i,j)**2+ uyb1(i,j)**2) / sqrt(uxb1ramollo(i,j)**2 + uyb1ramollo(i,j)**2)
!~                else
!~                              write(1034,*) 1.
!~                      endif
!~              enddo           
!~      enddo

!~   print*,'apres calcul rempli_L2'
!~   read(*,*)
    
    call interpol_glflux ! calcul flux GL + interpolation sur voisins
    
  endif
  
  call rempli_L2(nxd1,nxd2,nyd1,nyd2,uxbar(nxd1:nxd2,nyd1:nyd2),uybar(nxd1:nxd2,nyd1:nyd2), & 
       uxb1(nxd1:nxd2,nyd1:nyd2),uyb1(nxd1:nxd2,nyd1:nyd2),  &
       imx_diag(nxd1:nxd2,nyd1:nyd2),imy_diag(nxd1:nxd2,nyd1:nyd2),ifail_diagno)
  

  ! Dans rempli_L2

  !      uxprex(n1,n2)                                                 
  !      uyprec(n1,n2) vitesses de l'iteration precedente     
  !                                                                    
  !      uxnew(n1,n2)                                                   
  !      uynew(n1,n2) uynew resultat de cette iteration 
  !
  !      imx(n1,n2) masque pour imposer les vitesses ou leur dérivee    
  !      imy(n1,n2) masque pour imposer les vitesses ou leur dérivée   

  ! eventuellement le domaine n1,n2 peut etre un sous-domaine de nx,ny
  ! attention il faudra alors appeler avec des sous-tableaux

  ! Dans l'appel uxbar -> uxprec   et Uxb1 -> Uxnew

  !---------------------------------------------------------------------------

  if (ifail_diagno.gt.0) then
     write(6,*) ' Probleme resolution systeme L2. ifail=',ifail_diagno
     STOP
  endif

  !      nouvelles vitesses         
  !$OMP PARALLEL
  !$OMP WORKSHARE
  uxbar(:,:)=max(min(uxb1(:,:),V_limit),-V_limit)
  uybar(:,:)=max(min(uyb1(:,:),V_limit),-V_limit)
  !$OMP END WORKSHARE


  !    calcul de tobmx et tobmy (frottement basal) apres calcul des vitesses
  !    ---------------------------------------------------------------------

  if (gr_select.eq.1) then !Tsai
  !$OMP WORKSHARE
     where(gr_line_schoof(:,:).gt.0)
        tobmx(:,:) = sign( neffmx(:,:) * frot_coef , -uxbar(:,:))
        tobmy(:,:) = sign( neffmy(:,:) * frot_coef , -uybar(:,:))
     elsewhere
        tobmx(:,:) = -betamx(:,:) * uxbar(:,:)
        tobmy(:,:) = -betamy(:,:) * uybar(:,:)
     end where
  !$OMP END WORKSHARE
  else
  !$OMP WORKSHARE
     tobmx(:,:) = -betamx(:,:) * uxbar(:,:)
     tobmy(:,:) = -betamy(:,:) * uybar(:,:)
  !$OMP END WORKSHARE
  end if
  !$OMP BARRIER

  ! afq -- initMIP outputs:
  !$OMP DO
  do j=2,ny-1 
     do i=2,nx-1
        taub(i,j) = sqrt( ((tobmx(i,j)+tobmx(i+1,j))*0.5)**2 &
             + ((tobmy(i,j)+tobmy(i,j+1))*0.5)**2 )
     end do
  end do
  !$OMP END DO
  !$OMP WORKSHARE
  debug_3d(:,:,117) = taub(:,:)
  !$OMP END WORKSHARE

  
  ! Mise ne réserve des vitesses flgzmx et flgzmy
  !$OMP WORKSHARE
  where (flgzmx(:,:))
     uxflgz(:,:)=uxbar(:,:)
  elsewhere
     uxflgz(:,:)=0.
  endwhere

  where (flgzmy(:,:))
     uyflgz(:,:)=uybar(:,:)
  elsewhere
     uyflgz(:,:)=0.
  endwhere
  !$OMP END WORKSHARE
  !$OMP END PARALLEL
  !i=92
  !j=152
  !write(6,*) 'time',time, uxbar(92,152),gzmx(92,152),ilemx(92,152),flotmx(92,152), flgzmx(92,152)

  end do ! afq -- sur itour_nolin

  return
end subroutine diagnoshelf


!-------------------------------------------------------------------
subroutine calc_pvi

use module3d_phy, only: nx,ny,nz,itracebug,pvi,pvimin,abar,flot,gzmx,gzmy, & 
                        ilemx,ilemy,eps,taushelf,h,ramollo,debug_3d,pvm 
use module_choix, only: btt,visc  ! afq - 03/2023 - not clean 

implicit none
  
! calcule les viscosites integrees pvi et pvm
! loi polynomiale + couplage thermomécanique
!
!     Attention ne marche que si la loi est la loi en n=3 + n=1
!     y compris le pur glen (n=3) ou le pur Newtonien (n=1) 
! --------------------------------------------------------------------

!     les deformations sont supposées indépendantes de la profondeur
!     et sont calculées dans strain_rate (appelé par main)       

!     eps(i,j) -> eps
!     ttau -> tau (2eme invariant du deviateur des contraintes)
!     BT2 loi en n=3, phiphi loi en n=1 calculés dans flowlaw

!
     

! La viscosité est calculée partout y compris pour la glace posée (ou elle est moins
! précise qu'un calcul direct avec la loi de déformation.)
! Pour les noeuds posés mais ayant un voisin stream ou flottant, on calcule
! la viscosité avec stream/shelves
! le calcul se fait sur les noeuds majeurs

!$  integer :: rang ,nb_taches
!$  logical :: paral

integer i,j,k !loop integers

if (itracebug.eq.1)  call tracebug(' Calc pvi')

!$OMP PARALLEL PRIVATE(rang,iplus1,jplus1,sf3,sf1,BT2,phiphi,ttau,d02,discr)
!$ paral = OMP_IN_PARALLEL()
!$ rang=OMP_GET_THREAD_NUM()
!$ nb_taches=OMP_GET_NUM_THREADS()

!$OMP WORKSHARE
pvi(:,:)  = pvimin
Abar(:,:) = 0.
!$OMP END WORKSHARE

!$OMP DO
do j=1,ny
   do i=1,nx
      iplus1=min(i+1,nx)
      jplus1=min(j+1,ny)


      if (flot(i,j)) then    ! noeuds flottants
         sf3=sf03*sf_shelf      
         sf1=sf01*sf_shelf
     
      else if (gzmx(i,j).or.gzmx(iplus1,j).or.gzmy(i,j).or.gzmy(i,jplus1)) then
         sf1=sf01
         sf3=max(sf03,0.01)   ! pour les fleuves de glace, un peu de Glen 
     
      else if (ilemx(i,j).or.ilemx(iplus1,j).or.ilemy(i,j).or.ilemy(i,jplus1)) then
         sf1=sf01
         sf3=max(sf03,0.01)   ! pour les iles aussi

      else
!         sf1=1
!         sf3=1
         sf1=sf01             ! pour la viscosite anisotrope (ici en longitudinal)
         sf3=sf03
      endif


      do k=1,nz 

         BT2=BTT(i,j,k,1)*sf3    ! changement du sf
         phiphi=BTT(i,j,k,2)*sf1 !  changement du sf



         if (BT2.lt.1.e-25) then      ! pur newtonien
            visc(i,j,k)=1./phiphi
            ttau=2.*visc(i,j,k)*eps(i,j)
         else                         ! polynomial


!  en mettant Bt2 en facteur
            d02=eps(i,j) 
            discr=((phiphi/3.)**3.)/Bt2+d02**2
            discr=discr**0.5

            ttau=(d02+discr)**(1./3.)-(discr-d02)**(1./3.)


            ttau=ttau*Bt2**(-1./3.)


            visc(i,j,k)=Bt2*ttau*ttau+phiphi

            if (visc(i,j,k).gt.1.e-15) then 
               visc(i,j,k)=1./visc(i,j,k)
            else
               visc(i,j,k)=1.e15
            endif
         endif
         pvi(i,j)=pvi(i,j)+visc(i,j,k)
         Abar(i,j) =(Bt2/2.)**(-1./3.) + Abar(i,j)
         
         Taushelf(i,j)=Taushelf(i,j)+ttau


      end do
     

    
      pvi(i,j)  = pvi(i,j)*H(i,j)/nz
      Abar(i,j) = (Abar(i,j) /nz)**(-3.)
      

      Taushelf(i,j)=Taushelf(i,j)/nz
    
    
   end do
end do
!$OMP END DO

! cas des noeuds fictifs, si l'épaisseur est très petite 
! pvimin est très petit
!$OMP WORKSHARE
where (H(:,:).le.1.)
   pvi(:,:) = pvimin
end where

where (ramollo(:,:).ge..5)
   pvi(:,:) = pvimin
end where


debug_3D(:,:,27)=pvi(:,:)
!$OMP END WORKSHARE
! attention run 35
!-------------------- 
!!$if (time.gt.10.) then
!!$   where (flot(:,:))
!!$      pvi(:,:)=pvimin
!!$   end where
!!$end if

!  calcul de la viscosite integree au milieu des mailles (pvm)
!$OMP DO
do i=2,nx
   do j=2,ny

! les lignes suivantes pour un pvm moyenne des pvi
      pvm(i,j)=0.25*((pvi(i,j)+pvi(i-1,j-1))+    &   
           (pvi(i,j-1)+pvi(i-1,j)))

   end do
end do
!$OMP END DO
!$OMP END PARALLEL

end subroutine calc_pvi
!------------------------------------------------------------------

subroutine imx_imy_nx_ny

use module3d_phy, only: imx_diag,imy_diag,flgzmx,flgzmy 
         
implicit none

! definition des masques
! pour rempli_L2 : calcule les masques imx et imy qui 
! donnent les cas de conditions aux limites
! version pour travailler sur tout le domaine nx ny
!----------------------------------------------------



!   -34                -3   Nord              -23    
!     !----------------------------------------!
!     !                                        !
!     !              1 (prescrite)             !
!  -4 !                   ou                   ! -2
! West!              2 (L2)                    !  Est
!     !                                        !
!     !                                        !
!     !----------------------------------------!
!    -41               -1  Sud                -12 

!$OMP PARALLEL
!$OMP WORKSHARE
imx_diag(:,:)=0
imy_diag(:,:)=0

! a l'interieur du domaine
!-------------------------

where (flgzmx(:,:))
   imx_diag(:,:)=2             ! shelf ou stream        
elsewhere
   imx_diag(:,:)=1             ! vitesse imposee
end where

where (flgzmy(:,:))      
   imy_diag(:,:)=2             ! shelf ou stream
elsewhere
   imy_diag(:,:)=1             ! vitesse imposee
end where

! bord sud 
imx_diag(:,1)=-1
imy_diag(:,2)=-1

! bord nord
imx_diag(:,ny)=-3
imy_diag(:,ny)=-3

! bord Est
imx_diag(1,:)=0    ! hors domaine a cause des mailles alternees
imx_diag(2,:)=-4
imy_diag(1,:)=-4

! bord West
imx_diag(nx,:)=-2
imy_diag(nx,:)=-2

! Coins
imx_diag(2,1)=-41       ! SW
imy_diag(1,2)=-41

imx_diag(nx,1)=-12      ! SE
imy_diag(nx,2)=-12

imx_diag(nx,ny)=-23     ! NE
imy_diag(nx,ny)=-23

imx_diag(2,ny)=-34      ! NW
imy_diag(1,ny)=-34

! hors domaine
imx_diag(1,:)=0     ! hors domaine a cause des mailles alternees
imy_diag(:,1)=0     ! hors domaine a cause des mailles alternees
!$OMP END WORKSHARE
!$OMP END PARALLEL
 
end subroutine imx_imy_nx_ny


subroutine imx_imy_nx_ny_reduce(choix)

use module3d_phy, only: flot,imx_diag,imy_diag,gr_line 
         
implicit none

!afq -- For the backforce computation we do not need to compute the velocities everywhere
!       We simply compute the velocities around (nneigh) the grounding line

  integer, intent(in) :: choix

  integer i,j,nvx,nvy
  integer bordouest,bordest,bordsud,bordnord

  if (choix.eq.1) then ! we reduce the extent of velocity computation

  imx_reduce(:,:) = 1
  imy_reduce(:,:) = 1

  where (flot(:,:)) 
     imx_reduce(:,:) = imx_diag(:,:)
     imy_reduce(:,:) = imy_diag(:,:)
  endwhere

  do j=1,ny
     do i=1,nx
        if (gr_line(i,j).eq.1) then
           bordouest=max(i-nneigh,1)
           bordest=min(i+nneigh,nx)
           bordsud=max(j-nneigh,1)
           bordnord=min(j+nneigh,ny)
           do nvx=bordouest,bordest
              do nvy=bordsud,bordnord
                 imx_reduce(nvx,nvy)=imx_diag(nvx,nvy)
                 imy_reduce(nvx,nvy)=imy_diag(nvx,nvy)
              enddo
           enddo
        endif
     enddo
  enddo

! -- we need to keep imx_diag for the domain edges for land terminating geometries
  imx_reduce(:,1) = imx_diag(:,1)
  imx_reduce(:,2) = imx_diag(:,2)
  imx_reduce(:,ny) = imx_diag(:,ny)
  imx_reduce(1,:) = imx_diag(1,:)
  imx_reduce(2,:) = imx_diag(2,:)
  imx_reduce(nx,:) = imx_diag(nx,:)

  imy_reduce(:,1) = imy_diag(:,1)
  imy_reduce(:,2) = imy_diag(:,2)
  imy_reduce(:,ny) = imy_diag(:,ny)
  imy_reduce(1,:) = imy_diag(1,:)
  imy_reduce(2,:) = imy_diag(2,:)
  imy_reduce(nx,:) = imy_diag(nx,:)


  else ! we do not reduce the extent
     imx_reduce(:,:)=imx_diag(:,:)
     imy_reduce(:,:)=imy_diag(:,:)
  endif

end subroutine imx_imy_nx_ny_reduce


!___________________________________________________________________________
! pour imposer les conditions de mismip sur les bords du fleuve
! a appeler apres imx_imy_nx_ny

subroutine mismip_boundary_cond

use module3d_phy, only: itracebug,imx_diag,imy_diag,uxbar,uybar 

implicit none

if (itracebug.eq.1)  call tracebug(' Subroutine mismip_boundray_cond')

! Condition pas de flux sur les bords nord et sud

 imy_diag(:,2)      = 1
 imy_diag(:,3)      = 1
 imy_diag(:,ny-1)   = 1
 imy_diag(:,ny)     = 1


 Uybar(:,2)    = 0.
 Uybar(:,3)    = 0.
 Uybar(:,ny-1) = 0.
 Uybar(:,ny)   = 0.


! condition pas de cisaillement sur les bords nord et sud
imx_diag(:,2)    = -1
imx_diag(:,ny-1) = -3

! coins
imx_diag(2,2)     =  -41
imx_diag(nx,2)    =  -12
imx_diag(2,ny-1)  =  -34
imx_diag(nx,ny-1) =  -23
imx_diag(1,:)     =    0     ! ces points sont hors grille

end subroutine mismip_boundary_cond

!end module diagno_L2_mod
end module diagno_mod