source: branches/iLoveclim/SOURCES/Old-sources/diagno-ant-0.6_mod.f90 @ 77

!> \file diagno-ant-0.6_mod.f90
!! TOOOOOOOOOOO DOOOOOOOOOOOOO
!<

!> \namespace diagno_mod 
!! TOOOOOOOOOOOOO DOOOOOOOOOOO
!! \author Cat
!! \date Avril 2007
!! @note Used module
!! @note   - use module3D_phy
!! @note   - use module_choix
!! @note   - use eq_elliptique_mod 
!<


module diagno_mod               ! mise en module Cat avril 2007

use module3D_phy
use module_choix
use eq_elliptique_mod    
implicit none



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

integer :: itour_pvi
integer :: iplus1,jplus1
integer ::  ctvisco,iumax,jumax
real  :: delumax,errmax
real  :: phiphi,bt2,d02,discr,ttau
character :: numvisco*1
real :: sf3,sf1,epsxxm,epsyym,epsm,sf01,sf03
real :: viscm


logical :: stopvisco,viscolin
logical :: test_visc

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

namelist/diagno_rheol/ sf01,sf03

! 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,*)'/'     
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,*)


!      Precision utilisee dans de calcul
prec = 1.e-2
itour_pvi=0

return
end subroutine init_diagno

!------------------------------------------------------------------------------------
subroutine diagnoshelf !(ifail)

  integer :: itest

  if (itracebug.eq.1)  call tracebug(' Entree dans diagnoshelf')
  !      Resolution numerique des equations diagnostiques

  ! initialisation epsilon
  epsilon(:,:)=1.e-5

  itour_pvi=itour_pvi+1



  !      initialisation de un et oldu
  !      --------------------------------------

  do i=1,nx-1
     do j=1,ny-1
        un(i,j) = sqrt(((uxbar(i,j)+uxbar(i+1,j))/2)**2    & 
             + ((uybar(i,j)+uybar(i,j+1))/2.)**2)     
     end do
  end do

  numvisco=' '

  ! loi polynomiale + couplage thermomécanique
  !--------------------------------------------------------------------
  !     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
  !     Attention ne marche que si la loi est la loi en n=3 + n=1
  !

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

  ! 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


  do j=1,ny
     do i=1,nx
        iplus1=min(i+1,nx)
        jplus1=min(j+1,ny)



!!$      test_visc=(flot(i,j).or.flgzmx(i,j).or.flgzmx(iplus1,j).or.     & 
!!$           flgzmy(i,j).or.flgzmy(i,jplus1) )
!!$
!!$        
!!$      if (test_visc) then
!!$         sf3=sf03      
!!$         sf1=sf01      
!!$      else
!!$         sf1=1.
!!$         sf3=1.
!!$      endif


        if (flot(i,j)) then    ! noeuds flottants
           sf3=sf03      
           sf1=sf01     
        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=3
        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
              visc(i,j,k)=1./phiphi
              ttau=2.*visc(i,j,k)*eps(i,j)
           else


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

           Taushelf(i,j)=Taushelf(i,j)+ttau

        end do



        pvi(i,j)=pvi(i,j)*H(i,j)/nz

        Taushelf(i,j)=Taushelf(i,j)/nz



     end do
  end do


  !      calcul de la viscosite integree au milieu des mailles (pvm)

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



        ! le point ij de pvm est au milieu des points majeurs
        !  (i,j),(i-1,j),(i,j-1,i-1,j-1)

        ! calcul de epsm sur les noeuds pvm

        ! ATTENTION A l'ATTRIBUTION DES SF


!!$
!!$
!!$
!!$      epsxxm=0.25*((epsxx(i,j)+epsxx(i-1,j-1))+     & 
!!$           (epsxx(i,j-1)+epsxx(i-1,j)))
!!$      epsyym=0.25*((epsyy(i,j)+epsyy(i-1,j-1))+     &      
!!$           (epsyy(i,j-1)+epsyy(i-1,j)))
!!$
!!$      epsm=epsxxm**2+epsyym**2+epsxxm*epsyym+epsxy(i,j)**2
!!$      epsm=epsm**0.5
!!$
!!$!
!!$      pvm(i,j)=0.
!!$      do k=1,nz 
!!$         BT2=((BTT(i,j,k,1)+BTT(i-1,j-1,k,1))+         &  
!!$              (BTT(i,j-1,k,1)+BTT(i-1,j,k,1)))*sf3*0.25
!!$         phiphi=((BTT(i,j,k,2)+BTT(i-1,j-1,k,2))+      &
!!$              (BTT(i,j-1,k,2)+BTT(i-1,j,k,2)))*sf1*0.25
!!$           
!!$
!!$
!!$! nouvelle version en mettant Bt2 en facteur et avec cas bt2=0
!!$
!!$         if (BT2.lt.1.e-25) then
!!$            viscm=1./phiphi
!!$         else
!!$
!!$
!!$            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.)
!!$            viscm=Bt2*ttau*ttau+phiphi
!!$            if (viscm.gt.1.e-15) then 
!!$               viscm=1./viscm
!!$            else
!!$               viscm=1.e15
!!$            endif
!!$         endif
!!$
!!$    
!!$         pvm(i,j)=pvm(i,j)+viscm
!!$      end do
!!$       
!!$      pvm(i,j)=pvm(i,j)*(Hmx(i,j)+HMX(i,j-1))/nz*0.5


        ! pour découplage des contraintes de cisaillement horizontal. Cat 30 janvier 2005
        !         pvm(i,j)=1.E6  ! donne n'importe quoi
        !         pvm(i,j)=pvm(i,j)/2.

     end do
  end do


  !      Remplissage de la matrice creuse, rearrangement et resolution:
  !      --------------------------------------------------------------

  !      plus de decoupage de l'Antarctique en 4 domaines de tailles inegales
  !      selon la geographie pour eviter les iterations    



  CALL FLOTTAB_RESCUE !test pour voir resultats avec old method
  !-----------------------------------------------------------------------

!!$! -------------------------------------------------test de symetrie
!!$call detect_assym(nx,ny,0,41,1,0,1,0,Uxbar,itest)
!!$
!!$if (itest.gt.0) then
!!$   write(6,*) 'avant remplidom asymetrie sur uxbar  pour time=',time
!!$   stop
!!$else
!!$!   write(6,*) 'apres icethick pas d asymetrie sur H  pour time=',time
!!$end if
  ! ----------------------------------------------fin test de symetrie

  ! call flottab_rescue    !Attention modif juillet 2007, pour test


  call remplidom(nx,ny,uxbar,uybar,uxb1,uyb1,ifail)

!!$! -------------------------------------------------test de symetrie
!!$call detect_assym(nx,ny,0,41,1,0,1,0,Uxb1,itest)
!!$
!!$if (itest.gt.0) then
!!$   write(6,*) 'apres remplidom asymetrie sur uxb1  pour time=',time
!!$   stop
!!$else
!!$!   write(6,*) 'apres icethick pas d asymetrie sur H  pour time=',time
!!$end if
!!$! ----------------------------------------------fin test de symetrie




  if (ifail.gt.0) then
     print*,'call flottab_rescue at time',time
     call flottab_rescue
     call remplidom(nx,ny,uxbar,uybar,uxb1,uyb1,ifail)
  endif

  !      nouvelles vitesses          attention les deux lignes suivants semblaient manquer
  uxbar(:,:)=uxb1(:,:)
  uybar(:,:)=uyb1(:,:)

  do i=2,nx-1
     do j=2,ny-1
        un(i,j) = sqrt(((uxbar(i,j)+uxbar(i+1,j))/2)**2       &  
             + ((uybar(i,j)+uybar(i,j+1))/2.)**2)
     end do
  end do



  !    calcul de tobmx et tobmy (frottement basal) apres calcul des vitesses
  !    ---------------------------------------------------------------------
  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

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

  where (flgzmy(:,:))
     uyflgz(:,:)=uybar(:,:)
  elsewhere
     uyflgz(:,:)=0.
  endwhere



  return
end subroutine diagnoshelf

end module diagno_mod