!> \file dragging_calc_beta.f90
!! Module qui calcule le beta a partir de vitesses de bilan
!<

!> \namespace  dragging_calc_beta
!! Calcule le beta a partir de vitesses de bilan
!! @note Il faut partir d'un cptr.
!! \author ...
!! \date ...
!! @note Used module
!! @note   - use module3D_phy
!<
module dragging_calc_beta

  ! Calcule le beta a partir de vitesses de bilan
  ! il faut partir d'un cptr.


  use module3d_phy
  use interface_input

  implicit none

  logical, dimension(nx,ny) :: gz_centre     !< stream on major node
  real, dimension(nx,ny) :: Vcol_x           !< vertically averaged velocity x direction (balance)
  real, dimension(nx,ny) :: Vcol_y           !< vertically averaged velocity y direction (balance)
  real, dimension(nx,ny) :: Vcol2            !<  vertically averaged velocity norme² on major nodes
  real, dimension(nx,ny) :: beta_centre      !< beta on major node (average)
  real :: beta_limgz                         !< when beta gt beta_limgz -> not gzmx
  real :: beta_bord                          !< for the sides of ice streams
  real :: ubil_limgz                         !< when ubil < ubil_limgz  -> not gzmx
  real :: coefbeta                           !< coefficient to ajust beta
  real :: ecart_quad                         !< somme of quadratic difference 
  real :: umag2                              !< square of veloc. at major node
  integer :: stagger=1                       !< iteration based on centered or staggered nodes

contains
  !-----------------------------------------------------------------------------------
  subroutine lect_vitbil  !  idem module spinup


    character(len=100) :: balance_vel_file        !  balance velocity file
    !  vitesses de bilan consideres comme la donnee
    !  sur laquelle calibrer

    namelist/vitbil_upwind/balance_vel_file    

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

    ! lecture des parametres du run                    
    !--------------------------------------------------------------------

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

    write(num_rep_42,428) '!___________________________________________________________' 
    write(num_rep_42,428) '!read balance velocities on staggered grid'
    write(num_rep_42,vitbil_upwind)                    
    write(num_rep_42,428) '! balance_vel_file : nom du fichier qui contient les vitesse de bilan'
    write(num_rep_42,*)

    ! read balance velocities on staggered node

    balance_vel_file  = trim(dirnameinp)//trim(balance_vel_file)

    call lect_input(2,'drag_m',1,drag_mx,balance_vel_file,trim(dirnameinp)//trim(runname)//'.nc')
    call lect_input(2,'drag_m',2,drag_my,balance_vel_file,trim(dirnameinp)//trim(runname)//'.nc')
    !call lect_datfile(nx,ny,Vcol_x,1,balance_vel_file)   ! read Vcol_x   
    !call lect_datfile(nx,ny,Vcol_y,2,balance_vel_file)   ! read Vcol_y 



    do j=2,ny-1
       do i=2,nx-1
          umag2=((Vcol_x(i,j)+Vcol_x(i+1,j))*(Vcol_x(i,j)+Vcol_x(i+1,j))) &  ! ux2
               +((Vcol_y(i,j)+Vcol_y(i,j+1))*(Vcol_y(i,j)+Vcol_y(i,j+1)))    ! uy2
          Vcol2(i,j) =umag2*0.25
       end do
    end do

    ! limit the maximum value

    Vcol_x(:,:)=max(-3900.,Vcol_x(:,:))
    Vcol_x(:,:)=min( 3900.,Vcol_x(:,:))
    Vcol_y(:,:)=max(-3900.,Vcol_y(:,:))
    Vcol_y(:,:)=min( 3900.,Vcol_y(:,:))


    debug_3D(:,:,59)=Vcol_x(:,:)
    debug_3D(:,:,60)=Vcol_y(:,:)

  end subroutine lect_vitbil
  !-----------------------------------------------------------------------------------------

  subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.

    implicit none



    if (itracebug.eq.1)  call tracebug(' Calc_beta         subroutine init_dragging')
    mstream_mx(:,:)=1
    mstream_my(:,:)=1

    beta_limgz=5.e4     
    beta_bord= 1000.        ! 
    ubil_limgz=3.

    ! coefficient permettant de modifier le basal drag.
    drag_mx(:,:)=1.
    drag_my(:,:)=1.

    call lect_vitbil

    iter_beta = 1
    iter_loop : if (iter_beta.eq.1) then  ! pour calculer initial guess de beta
       betamx(:,:)=0
       betamy(:,:)=0
       coefbeta=5.

    else if (iter_beta.eq.2) then  ! pour partir d'un beta tres fort (ne marche pas)
       betamx(:,:)=1.e5
       betamy(:,:)=1.e5
       coefbeta=10.

    else if (iter_beta.eq.3) then  ! iteration Arthern en partant de l'initial guess

       call lect_input(2,'betam',1,betamx,'beta-estime.dat',trim(dirnameinp)//trim(runname)//'.nc')
       call lect_input(2,'betam',2,betamy,'beta-estime.dat',trim(dirnameinp)//trim(runname)//'.nc')     
       !call lect_datfile(nx,ny,betamx,1,'beta-estime.dat')
       !call lect_datfile(nx,ny,betamy,2,'beta-estime.dat')

       ! average
       beta_centre(:,:)=0.
       do j=2,ny-1
          do i=2,nx-1
             beta_centre(i,j)= ((betamx(i,j)+betamx(i+1,j)) &
                  + (betamy(i,j)+betamy(i,j+1)))*0.25
          end do
       end do

       ! redistribute on staggered grid
       if (stagger.eq.0) then
          do j=2,ny
             do i=2,nx
                betamx(i,j)= min(betamx(i,j),(beta_centre(i-1,j)+beta_centre(i,j))*0.5)
                betamy(i,j)= min(betamx(i,j),(beta_centre(i,j-1)+beta_centre(i,j))*0.5)
             end do
          end do
       end if

       if (stagger.eq.2) then 
          do j=2,ny
             do i=2,nx
                if (abs(Vcol_x(i,j)).gt.50.) then     ! fleuve de glace
                   betamx(i,j+1)=min(betamx(i,j+1),beta_bord)
                   betamx(i,j-1)=min(betamx(i,j-1),beta_bord)
                   betamx(i-1,j)=min(betamx(i-1,j),beta_bord)
                   betamx(i+1,j)=min(betamx(i+1,j),beta_bord)
                end if
                if (abs(Vcol_y(i,j)).gt.50.) then     ! fleuve de glace
                   betamy(i+1,j)=min(betamy(i+1,j),beta_bord)
                   betamy(i-1,j)=min(betamy(i-1,j),beta_bord)
                   betamy(i,j+1)=min(betamy(i,j+1),beta_bord)
                   betamy(i,j-1)=min(betamy(i,j-1),beta_bord)
                end if
             end do
          end do
       end if
       coefbeta=1.e-3

    end if iter_loop
    call gzm_betacalc


    return

  end subroutine init_dragging
  !________________________________________________________________________________


  !-------------------------------------------------------------------------
  subroutine dragging   ! defini la localisation des streams et le frottement basal

    implicit none

!!$where (flgzmy(:,:))
!!$   debug_3D(:,:,84)=1
!!$elsewhere
!!$   debug_3D(:,:,84)=0
!!$endwhere


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

    call gzm_betacalc

!!$gzmx(:,:)=.true.
!!$gzmy(:,:)=.true.
!!$flgzmx(:,:)=.false.
!!$flgzmy(:,:)=.false.
!!$
!!$
!!$! points flottants : flgzmx mais pas gzmx
!!$do j=2,ny
!!$   do i=2,nx
!!$      if (flot(i,j).and.(flot(i-1,j))) then
!!$         flgzmx(i,j)=.true.
!!$         gzmx(i,j)=.false.
!!$      end if
!!$      if (flot(i,j).and.(flot(i,j-1))) then
!!$         flgzmy(:,:)=.true.
!!$         gzmy(i,j)=.false.
!!$      end if
!!$   end do
!!$end do

    ! pour źtre sūr d'avoir quelques points "Dirichlet"
!!$gzmx(:,:)=gzmx(:,:).and.(abs(vcol_x(:,:)).gt.ubil_limgz)
!!$gzmy(:,:)=gzmy(:,:).and.(abs(vcol_y(:,:)).gt.ubil_limgz )

    ! si la vitesse de deformation est deja plus grande que la vitesse de bilan
    !gzmx(:,:)=gzmx(:,:).and.(abs(uxdef).lt.abs(vcol_x(:,:)))
    !gzmy(:,:)=gzmy(:,:).and.(abs(uydef).lt.abs(vcol_y(:,:)))

    flgzmx(:,:)=flgzmx(:,:).or.gzmx(:,:)
    flgzmy(:,:)=flgzmy(:,:).or.gzmy(:,:)

    ! pour initialisation, a terme pourrait etre dans une boucle iterative


    if (iter_beta.eq.1) then 
       betamx(:,:)=0.
       betamy(:,:)=0.
    end if


    Rob:if (iter_beta.gt.1) then
       iter_beta=iter_beta+1
       if (itracebug.eq.1) then 
          call tracebug(' Dragging : iter=')
          write(num_tracebug,*) iter_beta
       endif

       ! iteration Arthern like sur beta
       !________________________________________
       ! On part d un beta issus de calc beta et eventuellement lisse par moyenne mobile
       !
       ! Vcol_x tient le role de solution Dirichlet et uxbar de solution Neumann

       ecart_quad=0.
       if (iter_beta.ge.4) then

          stag: if (stagger.eq.1) then
             do j=2,ny
                do i=2,nx
                   if (gzmx(i,j)) then
                      betamx(i,j)=betamx(i,j)-coefbeta*(Vcol_x(i,j)*Vcol_x(i,j)-uxbar(i,j)*uxbar(i,j))
                      betamx(i,j)=max(betamx(i,j),0.)
                      ecart_quad= ecart_quad  & 
                           + abs(Vcol_x(i,j)*Vcol_x(i,j)-uxbar(i,j)*uxbar(i,j))
                   end if
                   if (gzmy(i,j)) then
                      betamy(i,j)=betamy(i,j)-coefbeta*(Vcol_y(i,j)*Vcol_y(i,j)-uybar(i,j)*uybar(i,j))
                      betamy(i,j)=max(betamy(i,j),0.)
                      ecart_quad= ecart_quad  & 
                           + abs(Vcol_y(i,j)*Vcol_y(i,j)-uybar(i,j)*uybar(i,j))
                   end if
                end do
             end do
          else if (stagger.eq.0) then          ! le chemin se fait sur les noeuds centres
             do j=2,ny-1
                do i=2,nx-1
                   if (gz_centre(i,j)) then
                      umag2=((uxbar(i,j)+uxbar(i+1,j))*(uxbar(i,j)+uxbar(i+1,j))) &  ! ux2
                           +((uybar(i,j)+uybar(i,j+1))*(uybar(i,j)+uybar(i,j+1)))    ! uy2
                      umag2=umag2*0.25

                      beta_centre(i,j)=beta_centre(i,j)-coefbeta*(Vcol2(i,j)-umag2) 
                      beta_centre(i,j)=max(beta_centre(i,j),0.)            
                      ecart_quad= ecart_quad  & 
                           + abs(Vcol2(i,j)-umag2) 
                   end if
                end do
             end do
             ! redistribute on staggered grid
             do j=2,ny
                do i=2,nx
                   betamx(i,j)= (beta_centre(i-1,j)+beta_centre(i,j))*0.5
                   betamy(i,j)= (beta_centre(i,j-1)+beta_centre(i,j))*0.5
                end do
             end do

          else if (stagger.eq.2) then          ! le chemin se fait sur les noeuds centres
             beta_centre(:,:)=0.
             write(6,*) coefbeta
             do j=2,ny-1
                do i=2,nx-1
                   if (gz_centre(i,j)) then
                      umag2=((uxbar(i,j)+uxbar(i+1,j))*(uxbar(i,j)+uxbar(i+1,j))) &  ! ux2
                           +((uybar(i,j)+uybar(i,j+1))*(uybar(i,j)+uybar(i,j+1)))    ! uy2
                      umag2=umag2*0.25

                      beta_centre(i,j)=-coefbeta*(Vcol2(i,j)-umag2) 
                      beta_centre(i,j)=min(beta_centre(i,j),100.)
                      beta_centre(i,j)=max(beta_centre(i,j),-100.)
                      ecart_quad= ecart_quad  & 
                           + abs(Vcol2(i,j)-umag2) 
                   end if
                end do
             end do
             debug_3D(:,:,86)=beta_centre(:,:)

             ! redistribute on staggered grid
             do j=2,ny
                do i=2,nx
                   betamx(i,j)= betamx(i,j)+((beta_centre(i-1,j)+beta_centre(i,j)) &
                        + (beta_centre(i-1,j-1)+beta_centre(i,j-1)) &
                        + (beta_centre(i-1,j+1)+beta_centre(i,j+1)))/6. 

                   betamx(i,j)=max(betamx(i,j),0.)
                   betamx(i,j)=min(betamx(i,j),beta_limgz)

                   betamy(i,j)= betamy(i,j)+((beta_centre(i,j-1)+beta_centre(i,j)) &
                        + (beta_centre(i-1,j-1)+beta_centre(i-1,j)) &
                        + (beta_centre(i+1,j-1)+beta_centre(i+1,j)))/6. 

                   betamy(i,j)=max(betamy(i,j),0.)
                   betamy(i,j)=min(betamy(i,j),beta_limgz)
                end do
             end do

          else if (stagger.eq.-1) then  ! test dramatique !
             write(6,*) 'stagger', stagger, iter_beta
             do j=2,ny
                do i=2,nx
                   if (abs(uxbar(i,j))-abs(Vcol_x(i,j)).lt.-50.) then
                      !             write(6,*) i,j, betamx(i,j),betamx(i,j+1), betamx(i,j-1),uxbar(i,j),Vcol_x(i,j)
                      betamx(i,j)=min(betamx(i,j),100.)
                      !             write(6,*) i,j, betamx(i,j),betamx(i,j+1), betamx(i,j-1)
                   endif

                   if (abs(uybar(i,j))-abs(Vcol_y(i,j)).lt.-50.) then
                      betamy(i,j)=min(betamy(i,j),100.)
                   endif
                end do
             end do
             do j=2,ny
                do i=2,nx
                   if (betamx(i,j).eq.-1) then

                      betamx(i,j)=100.
                      betamx(i,j+1)=200
                      betamx(i,j-1)=200
                   end if
                   if (betamy(i,j).eq.-1) then
                      betamy(i,j)=100
                      betamy(i+1,j)=200
                      betamy(i-1,j)=200
                   end if
                end do
             end do
          endif stag
       else
          do j=2,ny
             do i=2,nx
                if (gzmx(i,j)) then
                   ecart_quad= ecart_quad  & 
                        + abs(Vcol_x(i,j)*Vcol_x(i,j)-uxbar(i,j)*uxbar(i,j))
                end if
                if (gzmy(i,j)) then
                   ecart_quad= ecart_quad  & 
                        + abs(Vcol_y(i,j)*Vcol_y(i,j)-uybar(i,j)*uybar(i,j))
                end if
             end do
          end do
       end if


       write(6,*) 'ecart_quad',ecart_quad
    end if Rob
    !-------------------------- fin iteration Arthern like--------------------------


    ! iteration multiplicative sur beta
    !________________________________________
    !On part dun beta tres fort et on arrete de le diminuer quand les vitesses sont ok
    !
    !   if (itracebug.eq.1) write(num_tracebug,*)'     coefbeta',coefbeta
    !
    !where (gzmx(:,:).and.(abs(uxbar(:,:)).gt.(abs(Vcol_x(:,:)))))   ! la vitesse est trop grande, augmenter le beta
    !   betamx(:,:)=betamx(:,:)*coefbeta*0.9                         ! mais moins vite pour eviter une oscillation
    ! betamx(:,:)=betamx(:,:)*coefbeta
    !elsewhere                                       !
    !   betamx(:,:)=betamx(:,:)/coefbeta
    !end where

    !where (gzmy(:,:).and.(abs(uybar(:,:)).gt.(abs(Vcol_y(:,:)))))   ! la vitesse est trop grande, augmenter le beta
    !   betamy(:,:)=betamy(:,:)*coefbeta*0.9                        ! mais moins vite pour eviter une oscillation
    !  betamy(:,:)=betamy(:,:)*coefbeta                              ! mais moins vite pour eviter une oscillation
    !elsewhere                                       !
    !   betamy(:,:)=betamy(:,:)/coefbeta
    !end where

!!$flgzmy(:,:)=.false.
!!$flgzmx(:,:)=.false.
!!$
!!$   do j=2,ny
!!$      do i=2,nx
!!$         if (flot(i,j).and.(flot(i-1,j))) then
!!$            gzmx(i,j)=.false.
!!$            flgzmx(i,j)=.true.
!!$            betamx(i,j)=0.
!!$         end if
!!$         if (flot(i,j).and.(flot(i,j-1))) then
!!$            gzmy(i,j)=.false.
!!$            flgzmy(i,j)=.true.
!!$            betamy(i,j)=0.
!!$         end if
!!$      end do
!!$   end do
!!$end if
!!$
!!$flgzmx(:,:)=flgzmx(:,:).or.gzmx(:,:)
!!$flgzmy(:,:)=flgzmy(:,:).or.gzmy(:,:)

!!$where (flgzmy(:,:))
!!$   debug_3D(:,:,84)=1
!!$elsewhere
!!$   debug_3D(:,:,84)=0
!!$endwhere



    if (itracebug.eq.1)  call tracebug(' Dragging sortie')

    return
  end subroutine dragging

  !----------------------------------------------------------------------------------
  subroutine gzm_betacalc

    ! calcul de gzmx qui ne sera plus modifie ensuite


    gzmx(:,:)=.true.
    gzmy(:,:)=.true.
    flgzmx(:,:)=.false.
    flgzmy(:,:)=.false.


    ! points flottants : flgzmx mais pas gzmx
    do j=2,ny
       do i=2,nx
          if (flot(i,j).and.(flot(i-1,j))) then
             flgzmx(i,j)=.true.
             gzmx(i,j)=.false.
             betamx(i,j)=0.
          end if
          if (flot(i,j).and.(flot(i,j-1))) then
             flgzmy(i,j)=.true.
             gzmy(i,j)=.false.
             betamy(i,j)=0.
          end if
       end do
    end do
    !--------- autres criteres

    ! points poses
    gzmx(:,:)=gzmx(:,:).and.(betamx(:,:).lt.beta_limgz)   !  Pas de calcul pour les points
    gzmy(:,:)=gzmy(:,:).and.(betamy(:,:).lt.beta_limgz)   !  au fort beta

    ! pour źtre sūr d'avoir quelques points "Dirichlet"
    gzmx(:,:)=gzmx(:,:).and.(abs(Vcol_x(:,:)).gt.ubil_limgz) ! Pas de calcul pour les pts lents
    gzmy(:,:)=gzmy(:,:).and.(abs(Vcol_y(:,:)).gt.ubil_limgz) ! Pas de calcul pour les pts lents


    flgzmx(:,:)=flgzmx(:,:).or.gzmx(:,:)
    flgzmy(:,:)=flgzmy(:,:).or.gzmy(:,:)

    gz_centre(:,:)=.false.
    do j=2,ny-1
       do i=2,nx-1
          gz_centre(i,j)=gzmx(i,j).or.gzmx(i+1,j).or.gzmy(i,j).or.gzmy(i,j+1)
       end do
    end do

    where (flgzmy(:,:))
       debug_3D(:,:,84)=1
    elsewhere
       debug_3D(:,:,84)=0
    endwhere


    if (itracebug.eq.1) write(num_tracebug,*)' Subroutine     gzm_betacalc     ', & 
         flgzmy(191,115), gzmx(218,126),gzmy(244,259),flgzmy(244,259),flotmx(3,3),H(3,3)

  end subroutine gzm_betacalc

end module dragging_calc_beta