source: branches/GRISLIv3/SOURCES/Draggings_modules/dragging_param_beta_mod.f90 @ 504

!> \file dragging_neff_contmaj_mod.f90
!! Defini les zones de stream avec 3 criteres 
!<

!> \namespace  dragging_neff_contmaj
!! Defini les zones de stream avec 3 criteres
!! \author ...
!! \date ...
!! @note Les trois criteres sont:
!! @note   * un critere sur la pression effective 
!! @note   * un critere de continuite depuis la cote
!! @note   * un critere sur la courbure du socle (si negatif, vallees) 
!! @note Used module
!! @note   - use module3D_phy
!! @note   - use param_phy_mod
!<

module dragging_param_beta

  ! Defini les zones de stream avec :
  ! * un critere sur la pression effective 
  ! * un critere de continuite depuis la cote
  ! * un critere sur la courbure du socle (si negatif, vallees) 

  use geography, only: nx,ny
  ! pour compilation de steps_time_loop :
  use fake_beta_iter_vitbil_mod, only:time_iter,time_iter_end,nb_iter_vitbil,beta_iter_vitbil

  implicit none

  real,dimension(nx,ny) :: beta_param ! local variable

  real :: betamin   ! betamin : (Pa) frottement mini sous les streams

  real :: beta_slope        ! = A in:     B = A x Neff ** K
  real :: beta_expo         ! = K in:     B = A x Neff ** K

  real,dimension(nx,ny) :: neff   ! pression effective noeuds majeurs
  real,dimension(nx,ny) :: hwatmx ! hauteur d'eau staggered grille - afq
  real,dimension(nx,ny) :: hwatmy ! hauteur d'eau staggered grille - afq

  real :: neffmin = 1.e5 ! as in neffect, but this should be the exact same variable

  real, dimension(nx,ny) :: Vcol_x           !< uniquement pour compatibilite avec spinup cat
  real, dimension(nx,ny) :: Vcol_y           !< uniquement pour compatibilite avec spinup cat
  real, dimension(nx,ny) :: Vsl_x            !< uniquement pour compatibilite avec spinup cat
  real, dimension(nx,ny) :: Vsl_y            !< uniquement pour compatibilite avec spinup cat
  logical :: corr_def = .false.              !< for deformation correction, pour compatibilite beta

contains
  !-------------------------------------------------------------------------------

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

    use module3d_phy, only:betamax,num_param,num_rep_42,inv_beta,mstream_mx,mstream_my,&
         drag_mx,drag_my,betamax_2d,niter_nolin
    use runparam, only:itracebug

    implicit none

    namelist/drag_param_beta/beta_slope,beta_expo,betamax,betamin

    if (itracebug.eq.1)  call tracebug(' dragging avec neff et slope')



    ! lecture des parametres du run                      block drag neff
    !--------------------------------------------------------------------

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

    write(num_rep_42,'(A)')'!___________________________________________________________' 
    write(num_rep_42,'(A)') '&drag_param_beta          ! nom du bloc dragging param beta'
    write(num_rep_42,*)
    write(num_rep_42,*) 'beta_slope            = ', beta_slope
    write(num_rep_42,*) 'beta_expo             = ', beta_expo
    write(num_rep_42,*) 'betamax               = ', betamax
    write(num_rep_42,*) 'betamin               = ', betamin
    write(num_rep_42,*)'/'                            
    write(num_rep_42,'(A)') '! Slope & expo of beta = - slope x Neff ** expo'
    
    inv_beta=0
    !-------------------------------------------------------------------
    ! masque stream
    mstream_mx(:,:)=1
    mstream_my(:,:)=1


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

    betamax_2d(:,:) = betamax

    niter_nolin = 1

    return
  end subroutine init_dragging
  !________________________________________________________________________________

  !> SUBROUTINE: dragging
  !! Defini la localisation des streams et le frottement basal
  !>
  !-------------------------------------------------------------------------
  subroutine dragging   ! defini le frottement basal

    use module3d_phy, only:neffmx,neffmy,hwater,betamx,betamy,betamax,flot,&
         beta_centre
    
    !$ USE OMP_LIB


    implicit none
    integer :: i,j

    !$OMP PARALLEL
    ! calcul de neff (pression effective sur noeuds majeurs)
    if (sum(neffmx(:,:)).le.0.) neffmx(:,:) =1.e8
    if (sum(neffmy(:,:)).le.0.) neffmy(:,:) =1.e8

    !$OMP DO
    do j=1,ny-1
       do i=1,nx-1
          neff(i,j)=(neffmx(i,j)+neffmx(i+1,j)+neffmy(i,j)+neffmy(i,j+1))/4
       enddo
    enddo
    !$OMP END DO
    !aurel, for the borders:
    !$OMP WORKSHARE
    neff(:,ny)=neffmin
    neff(nx,:)=neffmin
    ! calcul de hwat (staggered)
    hwatmx(:,:)=0.
    hwatmy(:,:)=0.
    !$OMP END WORKSHARE
    !$OMP DO
    do i=2,nx
       hwatmx(i,:)=(hwater(i-1,:)+hwater(i,:))/2.
    enddo
    !$OMP END DO
    !$OMP DO
    do j=2,ny
       hwatmy(:,j)=(hwater(:,j-1)+hwater(:,j))/2.
    enddo
    !$OMP END DO

    !$OMP WORKSHARE

    ! new parametrisation of beta on Neff:
    betamx(:,:)= beta_slope*(neffmx(:,:)**beta_expo)
    betamy(:,:)= beta_slope*(neffmy(:,:)**beta_expo)

    betamx(:,:)=max(betamx(:,:),betamin)
    betamy(:,:)=max(betamy(:,:),betamin)
    betamx(:,:)=min(betamx(:,:),betamax)
    betamy(:,:)=min(betamy(:,:),betamax)

    where ( hwatmx(:,:) .lt. 0.5 ) betamx(:,:) = betamax
    where ( hwatmy(:,:) .lt. 0.5 ) betamy(:,:) = betamax

    !$OMP END WORKSHARE

    ! points flottants
    !$OMP DO
    do j=2,ny
       do i=2,nx
          if (flot(i,j).and.(flot(i-1,j))) then
             betamx(i,j)=0.
          end if
          if (flot(i,j).and.(flot(i,j-1))) then
             betamy(i,j)=0.
          end if
       end do
    end do
    !$OMP END DO

    !$OMP DO
    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
    !$OMP END DO
    !$OMP END PARALLEL

    return
  end subroutine dragging
  !________________________________________________________________________________

  !> SUBROUTINE: mstream_dragging
  !! Defini la localisation des streams
  !>
  !-------------------------------------------------------------------------
  subroutine mstream_dragging   ! defini la localisation des streams

    use module3d_phy, only:fleuvemx,fleuvemy,gzmx,gzmy,flgzmx,flgzmy,flotmx,flotmy,flot,&
         betamx,betamy,betamax
    !$ USE OMP_LIB
    
    implicit none
    integer :: i,j

    !$OMP PARALLEL

    !$OMP WORKSHARE
    fleuvemx(:,:)=.false.
    fleuvemy(:,:)=.false.
    gzmx(:,:)=.true.
    gzmy(:,:)=.true.
    flgzmx(:,:)=.false.
    flgzmy(:,:)=.false.
    !$OMP END WORKSHARE


    ! calcul de gzmx 

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

    !--------- autres criteres

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

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

    fleuvemx(:,:)=gzmx(:,:)
    fleuvemy(:,:)=gzmy(:,:)
    !$OMP END WORKSHARE

    !$OMP END PARALLEL

    return
  end subroutine mstream_dragging

end module dragging_param_beta