source: branches/iLoveclim/SOURCES/Draggings_modules/beta_iter_vitbil_mod.f90 @ 187

!< \file beta_iter_vitbil_mod.f90
!< Module pour affiner le beta a partir d'une carte vitbil  sur les noeuds centres.
!< A appeler depuis steps_time_loop
!< a utiliser avec beta_prescr

!> \ beta_iter_vitbil
!! \author ...
!! \date ...
!! @note Used module
!! @note   - use module3D_phy
!<

module beta_iter_vitbil_mod

  use module3d_phy
  use interface_input
  use netcdf
  use io_netcdf_grisli

  implicit none

  real, dimension(nx,ny) :: Vcol_x           !< uniquement pour compatibilite
  real, dimension(nx,ny) :: Vcol_y           !<
  real, dimension(nx,ny) :: Vsl_x            !<
  real, dimension(nx,ny) :: Vsl_y            !<


  real, dimension(nx,ny) :: drag_centre      !< beta on major node (average)
  real :: beta_limgz                         !< when beta gt beta_limgz -> not gzmx
  real :: beta_min                           !< minimum value of beta
  real :: beta_mult                          !< coefficient of beta field
  integer  :: ill,jll,nmoy

  real :: maxi                               !< calcul correction dS a la grounding line
  real :: mini

  logical :: corr_def = .False.               !< for deformation correction (compatibility)


  real, dimension(nx,ny) :: Umag_vitbil             !< vitesse de bilan
  real, dimension(nx,ny) :: Uslid_vitbil            !< vitesse de bilan partie glissement
  real, dimension(nx,ny) :: Umag_direct             !< vitesse moyenne centree calculee par GRISLI
  real, dimension(nx,ny) :: Uslmag_direct           !< vitesse glissement centree calculee par GRISLI
  real, dimension(nx,ny) :: Udefmag_direct          !< vitesse deformation centree calculee par GRISLI
  real, dimension(nx,ny) :: driving_stress          !< driving stress

  character(len=100)     :: Umag_bil_file           !< fichier qui contient les donnees
  real                   :: time_iter               !< temps de debut des iterations
  real                   :: time_iter_end           !< temps de fin des iterations
  real                   :: time_reiter             !< nbr annees entre 2 iterations calcul beta
  integer                :: nb_iter_vitbil          !< nombre d'iteratiosn avant de changer de pas de temps
  real                   :: coef_iter_vitbil        !< coefficient <= 1 (rapport vitesses)
  real                   :: J_Umag                  !< estimateur ecart entre vitesses
  real                   :: J_Udef                  !< estimateur vitesse regions sans glissement
  integer                :: nb_umag                 !< nombre de points calcul J_Umag
  integer                :: nb_Udef                 !< nombre de points calcul J_Umag


contains

!-----------------------------------------------------------------------------------------
!< subroutine : init_beta_iter_vitbil
!< Cette routine fait l'initialisation de la procedure beta_iter_vitbil

subroutine init_beta_iter_vitbil

  implicit none
  
  real*8, dimension(:,:),   pointer  :: tab => null() !< tableau 2d real ecrit dans le fichier pour lect ncdf

    namelist/beta_iter_vitbil/time_iter,nb_iter_vitbil,coef_iter_vitbil,Umag_bil_file,time_reiter


    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
428 format(A)
    read(num_param,beta_iter_vitbil)

    write(num_rep_42,428) '!___________________________________________________________' 
    write(num_rep_42,428) '!  module beta_iter_vitbil'
    write(num_rep_42,beta_iter_vitbil)                    
    write(num_rep_42,428) '! ' 
    write(num_rep_42,428) '!  Umag_bil_file : fichier qui contient les donnees Umag_vitbil'
    write(num_rep_42,428) '!  time_iter     : temps de debut des iterations '
    write(num_rep_42,428) '!  nb_iter_vitbil: nombre diteratation a faire avant de changer de pas de temps ' 
    write(num_rep_42,*)   '!  coef_iter_vitbil : coefficient pour le rapport vitesses  '
    write(num_rep_42,*)   '!  time_reiter : nombre d annees entre 2 iterations'
    write(num_rep_42,428) '!___________________________________________________________' 

    time_iter = time_iter + tbegin       ! time_iter est le temps de relaxation (de l'ordre de 5 ans)
                                         ! ajouter tbegin pour ne pas dependre du temps de depart.
    time_iter_end = time_iter + 1.       ! 20. ans en version std Cat et Seb
    print*,'debug init_beta_iter_vitbil time_iter=', time_iter, time_iter_end, time

    Umag_bil_file = trim(dirnameinp)//trim(Umag_bil_file)

    write(6,*) Umag_bil_file
!    call lect_input(1,'Umag_vitbil',1,Umag_vitbil,Umag_bil_file,'')
    call Read_Ncdf_var('Umag_vitbil',Umag_bil_file,tab)
    Umag_vitbil(:,:) = tab(:,:)
    
    debug_3D(:,:,58) = Umag_vitbil(:,:)

  end subroutine init_beta_iter_vitbil




!______________________________________________________________________________________
!>SUBROUTINE:   beta_min_value 
!! valeur mini que peut avoir beta (en Pa an m-1) 
!<

  subroutine beta_min_value(val_betamin)

  real :: val_betamin     ! valeur mini que peut avoir beta (en Pa an m-1)

    drag_centre(:,:) = max(drag_centre(:,:),val_betamin)
    drag_mx(:,:) = max(drag_mx(:,:),val_betamin)
    drag_my(:,:) = max(drag_my(:,:),val_betamin)
    where(flot(:,:))
       drag_centre(:,:) = 0.
    end where

  end subroutine beta_min_value

!>SUBROUTINE:   beta_c2stag  
!! redistribute central values on staggered grid
!<

  subroutine beta_c2stag(nxx,nyy,R_mx,R_my,R_centre)              ! redistribute central values on staggered grid

    implicit none
    integer                                ::  nxx,nyy                ! dimension des tableaux
    real, dimension(nxx,nyy), intent(out)  ::  R_mx                   ! valeur du beta sur maille mx
    real, dimension(nxx,nyy), intent(out)  ::  R_my                   ! valeur du beta sur maille my
    real, dimension(nxx,nyy), intent(in)   ::  R_centre               ! valeur du beta sur maille centree

    do j=2,nyy 
       do i=2,nxx
          R_mx(i,j)= (R_centre(i-1,j)+R_centre(i,j))*0.5
          R_my(i,j)= (R_centre(i,j-1)+R_centre(i,j))*0.5
       end do
    end do

  end subroutine beta_c2stag


!-----------------------------------------------------------------------------------------
!< subroutine : beta_iter_vitbil
!< fait des iterations pour affiner la valeur de beta_centre pour s'ajuster sur vitbil
!-----------------------------------------------------------------------------------------
subroutine beta_iter_vitbil(m_iter)

implicit none
integer :: m_iter ! indice bloucle iteration


! calcul des vitesses cibles :
!if ((time.eq.time_iter).and.(time.gt.time_reiter).and.(m_iter.eq.1)) then
if ((abs(time-time_iter).lt.dtmin).and.(time.gt.time_reiter).and.(m_iter.eq.1)) then
  Umag_direct(:,:)    = (((uxbar(:,:)+eoshift(uxbar(:,:),shift=1,boundary=0.,dim=1))**2+ &        ! moyenne
                          (uybar(:,:)+eoshift(uybar(:,:),shift=1,boundary=0.,dim=2))**2)**0.5)*0.5
  Umag_vitbil(:,:)= min(Umag_direct(:,:) * min(1.e3,H(:,:)/H0(:,:)), 5000.)
  print*,'*****debug beta_iter_vitbil calcul de Umag_vitbil',time,m_iter
  do j=1,ny
    do i=1,nx
      write(266,*) Umag_vitbil(i,j)
    enddo
  enddo
endif

! calcule la vitesse centree venant du calcul direct

  Umag_direct(:,:)    = (((uxbar(:,:)+eoshift(uxbar(:,:),shift=1,boundary=0.,dim=1))**2+ &        ! moyenne
                          (uybar(:,:)+eoshift(uybar(:,:),shift=1,boundary=0.,dim=2))**2)**0.5)*0.5

  Uslmag_direct(:,:)  = (((ux(:,:,nz)+eoshift(ux(:,:,nz),shift=1,boundary=0.,dim=1))**2+ &        ! glissement 
                          (uy(:,:,nz)+eoshift(uy(:,:,nz),shift=1,boundary=0.,dim=2))**2)**0.5)*0.5

  Udefmag_direct(:,:) = max( (Umag_direct(:,:)- Uslmag_direct(:,:)) , 0.)                          ! deformation


! guess du beta dans les zones où GRISLI a une vitesse de glissement nulle

  call slope_surf
  driving_stress(:,:) = rog * H(:,:) * slope(:,:)   !  driving stress


  Uslid_vitbil(:,:)   = Umag_vitbil(:,:) - Udefmag_direct(:,:)  ! calcule la vitesse de glissement "bilan"


! Pas de glissement dans GRISLI

where ((Uslmag_direct(:,:).le.0.).and.(Uslid_vitbil(:,:).gt.1.e-3))     ! pas de glissement dans GRISLI
  beta_centre(:,:) = driving_stress(:,:) / Uslid_vitbil(:,:)            ! mais du glissement "bilan"
end where


where ((Uslmag_direct(:,:).le.0.).and.(Uslid_vitbil(:,:).le.1.e-3))     ! ni l'un ni l'autre
   Uslid_vitbil(:,:) = 0.
   beta_centre(:,:)  = beta_limgz
end where


! Glissement dans GRISLI

  where (Uslid_vitbil(:,:).le.1.e-3)    ! la vitesse de bilan est deja trop rapide
     Uslid_vitbil(:,:) = 0.
     beta_centre(:,:)  = beta_limgz

  elsewhere (Uslmag_direct(:,:).gt.0.)
     beta_centre(:,:) = beta_centre(:,:) * coef_iter_vitbil * Uslmag_direct(:,:) / Uslid_vitbil(:,:)

  end where

  where(flot(:,:))
     beta_centre(:,:) = 0.
  end where



  beta_centre(:,:) = min(beta_centre(:,:),beta_limgz)
  
    call  beta_c2stag(nx,ny,betamx,betamy,beta_centre)  ! redistribue sur les mailles alternees
    call  beta_min_value(beta_min)                      !  valeur mini que peut avoir beta (en Pa an m-1)
 

 debug_3D(:,:,59) = betamx(:,:)
 debug_3D(:,:,60) = beta_centre(:,:)

! estimateurs
J_Umag  = 0.
J_Udef   = 0.
nb_umag = 0
nb_udef = 0

do j=2,ny-1
   do i=2,nx-1

      if (.not.flot(i,j)) then                       ! calcul sur les points poses
         nb_umag = nb_umag + 1
         J_Umag = (Umag_direct(i,j) - Umag_vitbil(i,j))**2. + J_Umag

         if (Uslid_vitbil(i,j).le.1.e-3) then       ! calcul sur la partie deformation suelement
            nb_udef = nb_udef + 1
            J_Udef = (Umag_direct(i,j) - Umag_vitbil(i,j))**2. + J_Udef
         end if
      end if

   end do
end do
      J_Umag = (J_Umag**0.5) / nb_umag
      J_Udef = (J_Udef**0.5) / nb_udef


drag_centre(:,:) = beta_centre(:,:)
drag_mx(:,:)     = betamx(:,:)
drag_my(:,:)     = betamy(:,:)


end subroutine beta_iter_vitbil
!-----------------------------------------------------------------------------------------
end module beta_iter_vitbil_mod