source: trunk/SOURCES/New-remplimat/diagno-L2_mod.f90 @ 104

!module diagno_L2_mod               ! Nouvelle version, compatible remplimat 2008 Cat
module diagno_mod                   ! nom pendant les tests
 !$ USE OMP_LIB
use module3D_phy
use module_choix
     
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 :: 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

contains    

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

namelist/diagno_rheol/sf01,sf03,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,*) '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) '! 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,*)


!      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


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


  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



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 

!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
  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
!~                      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

  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)
       
  !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), & 
       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)
       
    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

    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(*,*)      
  endif
  

  ! 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(:,:)=uxb1(:,:)
  uybar(:,:)=uyb1(:,:)
  !$OMP END WORKSHARE


  !    calcul de tobmx et tobmy (frottement basal) apres calcul des vitesses
  !    ---------------------------------------------------------------------
  !$OMP DO
  do j=1,ny
     do i=1,nx
        tobmx(i,j)=-betamx(i,j)*uxbar(i,j)
        tobmy(i,j)=-betamy(i,j)*uybar(i,j)
     enddo
  enddo
  !$OMP END DO

  ! 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)

  return
end subroutine diagnoshelf


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

! 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 :: t1,t2,ir
  real                           :: temps, t_cpu_0, t_cpu_1, t_cpu, norme  

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

! 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
!___________________________________________________________________________
! pour imposer les conditions de mismip sur les bords du fleuve
! a appeler apres imx_imy_nx_ny

subroutine mismip_boundary_cond
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