source: branches/iLoveclim/SOURCES/dragging_neff_slope_mod.f90 @ 77

!> \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_neff_slope

! 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 module3d_phy
use param_phy_mod
use interface_input
use io_netcdf_grisli

implicit none

logical,dimension(nx,ny) :: fleuve
logical,dimension(nx,ny) :: cote
logical,dimension(nx,ny) :: slowssamx ! not actual stream, but ssa used as Ub
logical,dimension(nx,ny) :: slowssamy ! not actual stream, but ssa used as Ub
logical,dimension(nx,ny) :: slowssa   ! not actual stream, but ssa used as Ub

real,dimension(nx,ny) :: hires_slope ! slope comupted on a high resolution topography
character(len=100) :: slope_fich     ! fichier grille
character(len=100) :: file_ncdf      !< fichier netcdf issue des fichiers .dat

real :: valmax
integer :: imax,jmax
integer :: i_moins1,i_plus1,j_moins1,j_plus1
integer :: lmax=20
integer :: idep,jdep,iloc,jloc


real :: tostick   ! pour la glace posee
real :: betamin   ! betamin : (Pa) frottement mini sous les streams
real :: tob_ile    ! pour les iles
real :: cry_lim=50.  ! courbure limite pour le suivi des fleuves

real,dimension(nx,ny) :: neff ! pression effective noeuds majeurs
real :: seuil_neff    ! seuil sur la pression effective pour avoir glissement en zone sediment
real :: coef_gz       ! coef frottement zones stream std  (10)
real :: coef_ile      ! coef frottement zones iles   (0.1)

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.

implicit none

! local variables, defining the rugosity-enhanced flow
real :: expo_slope
real :: pente_min, pente_max

namelist/drag_neff_slope/cf,betamax,betamin,toblim,tostick,seuil_neff,coef_gz,coef_ile,slope_fich,expo_slope,pente_min,pente_max

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


! formats pour les ecritures dans 42
428 format(A)

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

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

write(num_rep_42,428)'!___________________________________________________________' 
write(num_rep_42,428) '&drag_neff_slope          ! nom du bloc dragging neff slope'
write(num_rep_42,*)
write(num_rep_42,*) 'cf               = ',cf
write(num_rep_42,*) 'betamax          = ', betamax
write(num_rep_42,*) 'betamin          = ', betamin
write(num_rep_42,*) 'toblim           = ', toblim
write(num_rep_42,*) 'tostick          = ', tostick
write(num_rep_42,*) 'seuil_neff       = ', seuil_neff
write(num_rep_42,*) 'coef_gz          = ', coef_gz
write(num_rep_42,*) 'coef_ile         = ', coef_ile
write(num_rep_42,'(A,A)') 'slope_fich = ', slope_fich
write(num_rep_42,*) 'expo_slope       = ', expo_slope
write(num_rep_42,*) 'pente_min        = ', pente_min
write(num_rep_42,*) 'pente_max        = ', pente_max
write(num_rep_42,*)'/'                            
write(num_rep_42,428) '! cf coefficient de la loi de frottement fonction Neff'
write(num_rep_42,428) '! betamax : (Pa) frottement maxi sous les streams '
write(num_rep_42,428) '! betamin : (Pa) frottement mini sous les streams '
write(num_rep_42,428) '! toblim : (Pa)  pour les iles '
write(num_rep_42,428) '! tostick : (Pa)  pour les points non flgzmx '
write(num_rep_42,428) '! seuil_neff (Pa) seuil sur la pression effective pour avoir glissement'
write(num_rep_42,428) '! coef_gz : coef frottement zones stream std'
write(num_rep_42,428) '! coef_ile : coef frottement zones iles'
write(num_rep_42,428) '! slope_fich : fichier pente bedrock'
write(num_rep_42,428) '! expo_slope : exposant fonction pente'
write(num_rep_42,428) '! pente_min : no flow enhancement below this'
write(num_rep_42,428) '! pente_max : maximum flow enhancement above this'
write(num_rep_42,*)

!tostick=1.e5   ! valeurs pour les points non flgzmx
tob_ile=betamax/2.
moteurmax=toblim

! betamax_2d depends on sub-grid slopes.
slope_fich=trim(dirnameinp)//trim(slope_fich)
call lect_input(1,'slope',1,hires_slope(:,:),slope_fich,file_ncdf)

! from slopes, we create an index between 0 & 1
! 1 is very mountainous, 0 is flat
hires_slope(:,:)=(max(min(hires_slope(:,:),pente_max),pente_min)-pente_min)/(pente_max-pente_min)
! now we compute the actual betamax used by the remplimat routines
! /!\ the slope is used to modify the beta where we have a temperate base,
! but NO ice stream... -> Slow SSA zone (SSA used to compute Ub)
betamax_2d(:,:)=max ( tostick * (1. - (1 - betamax / tostick ) * hires_slope(:,:)**(1./expo_slope)) , betamax )
!do j=1,ny
!   do i=1,nx
!      write(18745,*) betamax_2d (i,j)
!   enddo
!enddo

!-------------------------------------------------------------------
! masque stream
      mstream_mx(:,:)=1
      mstream_my(:,:)=1


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

where (socle_cry(:,:).lt.cry_lim)
   debug_3D(:,:,25)=1
elsewhere
   debug_3D(:,:,25)=0
endwhere

   
      return
      end subroutine init_dragging
!________________________________________________________________________________

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

!         les iles n'ont pas de condition neff mais ont la condition toblim
!         (ce bloc etait avant dans flottab)

!$OMP PARALLEL
!$OMP DO
do j=2,ny
   do i=2,nx
      ilemx(i,j)=ilemx(i,j).and.(abs(rog*Hmx(i,j)*sdx(i,j)).lt.toblim)
      ilemy(i,j)=ilemy(i,j).and.(abs(rog*Hmy(i,j)*sdy(i,j)).lt.toblim)
   end do
end do
!$OMP END DO

! le calcul des fleuves se fait sur les mailles majeures en partant de la cote




! le gzmx initial a ete calculé dans flottab pour les points a cheval sur la grounding line

! Partant de la cote, le fleuve est determine en remontant vers le point de hauteur d'eau max
! tant que ce point est superieur a hwatermax.

! Attention : ce systeme ne permet pas d'avoir des fleuves convergents
!             et les fleuves n'ont une épaisseur que de 1 noeud.
!$OMP WORKSHARE
fleuvemx(:,:)=.false.
fleuvemy(:,:)=.false.
fleuve(:,:)=.false.
cote(:,:)=.false.
!$OMP END WORKSHARE

! detection des cotes
!$OMP DO
do  j=2,ny-1  
   do i=2,nx-1  
      if ((.not.flot(i,j)).and. & 
      ((flot(i+1,j)).or.(flot(i,j+1)).or.(flot(i-1,j)).or.(flot(i,j-1)))) then
         cote(i,j)=.true.
      endif
   end do
end do
!$OMP END DO

! calcul de neff (pression effective sur noeuds majeurs)
!$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)=1.e5
neff(nx,:)=1.e5


!!$
!!$fleuve_maj: do j=2,ny-1  
!!$ifleuve:   do i=2,nx-1                     
!!$
!!$cote_detect :  if (cote(i,j)) then
!!$         idep=i
!!$         jdep=j
!!$
!!$         if (socle_cry(i,j).lt.0.) then    ! dans une vallee
!!$            fleuve(i,j)=.true. 
!!$         else
!!$            cote(i,j)=.false.
!!$            cycle ifleuve
!!$         endif
!!$
!!$suit : do l=1,lmax         ! debut de la boucle de suivi, lmax longueur maxi des fleuves
!!$           i_moins1=max(idep-1,2)
!!$           j_moins1=max(jdep-1,1)
!!$           i_plus1=min(idep+1,nx)
!!$           j_plus1=min(jdep+1,ny)
!!$           
!!$! recherche du max en suivant le socle le plus profond
!!$! * en excluant les points flottants
!!$! * et ceux qui sont deja tagges fleuves
!!$
!!$           valmax=1000.
!!$
!!$           do jloc=j_moins1,j_plus1
!!$              do iloc=i_moins1,i_plus1
!!$ 
!!$                 if ((B(iloc,jloc).lt.valmax)      & 
!!$                      .and.(.not.flot(iloc,jloc))  &
!!$                      .and.(.not.fleuve(iloc,jloc)).and.(socle_cry(iloc,jloc).lt.cry_lim)) then
!!$                    imax=iloc
!!$                    jmax=jloc
!!$                    valmax=B(iloc,jloc)
!!$                 endif
!!$              end do
!!$           end do
!!$
!!$         if ((hwater(imax,jmax).gt.hwatstream).and.(socle_cry(i,j).lt.cry_lim)) then
!!$            fleuve(imax,jmax)=.true.
!!$            idep=imax
!!$            jdep=jmax
!!$         else
!!$            fleuve(imax,jmax)=.false.
!!$            exit suit
!!$         end if
!!$ 
!!$      end do suit
!!$
!!$   end if cote_detect
!!$
!!$end do ifleuve
!!$end do fleuve_maj

! aurel, we add the neff threshold:
where ((neff(:,:).le.seuil_neff).and.(.not.flot(:,:)).and.(H(:,:).gt.1.)) fleuve(:,:)=.true.
!$OMP END WORKSHARE

!$OMP DO
do j=1,ny-1
   do i=1,nx-1
      if (fleuve(i,j)) then
         fleuvemx(i,j)=.true.
         fleuvemx(i+1,j)=.true.
         fleuvemy(i,j)=.true.
         fleuvemy(i,j+1)=.true.
      end if
   end do
end do
!$OMP END DO

! we look for the warm based points that will not be treated as stream (ub from SSA):
!$OMP WORKSHARE
slowssa(:,:)=.false.
slowssamx(:,:)=.false.
slowssamy(:,:)=.false.
!$OMP END WORKSHARE
!$OMP DO
do  j=1,ny
   do i=1,nx
      !if ((not(flot(i,j))).and.(hwater(i,j).gt.0.1)) slowssa(i,j)=.true.
      if ((not(flot(i,j))).and.(ibase(i,j).ne.1).and.(H(i,j).gt.1.)) slowssa(i,j)=.true.
   end do
end do
!$OMP END DO
!$OMP DO
do j=1,ny-1
   do i=1,nx-1
      if (slowssa(i,j)) then
         slowssamx(i,j)=.true.
         slowssamx(i+1,j)=.true.
         slowssamy(i,j)=.true.
         slowssamy(i,j+1)=.true.
      end if
   end do
end do
!$OMP END DO
!$OMP DO
do j=1,ny
   do i=1,nx

      ! the actual streams and the warm based points are gzmx now:
      if ( ((.not.ilemx(i,j)).and.(fleuvemx(i,j))) .or. ((.not.ilemx(i,j)).and.(slowssamx(i,j))) ) gzmx(i,j)=.true. 

      
! calcul du frottement basal (ce bloc etait avant dans neffect)

      if (cotemx(i,j)) then                        ! point cotier
         betamx(i,j)=cf*neffmx(i,j) 
         betamx(i,j)=min(betamx(i,j),betamax)

      else if ((gzmx(i,j)).and.(.not.cotemx(i,j))) then ! tous les points SSA

         if (fleuvemx(i,j)) then                 ! the actual streams
            betamx(i,j)=cf*neffmx(i,j)*coef_gz
            betamx(i,j)=min(betamx(i,j),betamax)
            betamx(i,j)=max(betamx(i,j),betamin)

            if (cf*neffmx(i,j).gt.betamax*2.) then  ! a stream that's becoming grounded...
               if (slowssamx(i,j)) then
                  betamx(i,j)=betamax_2d(i,j)
               else
                  gzmx(i,j)=.false.
                  betamx(i,j)=tostick
               endif
            endif

         else                    ! not an actual stream, SSA is used to compute Ub
            betamx(i,j)=betamax_2d(i,j)
         endif
         
      else if (ilemx(i,j)) then
         betamx(i,j)=cf*neffmx(i,j)*coef_ile 
         betamx(i,j)=min(betamx(i,j),tob_ile)       
      else if (flotmx(i,j)) then     ! flottant ou ile
         betamx(i,j)=0.
      else                                         ! grounded, SIA
         betamx(i,j)=tostick                       ! frottement glace posee (1 bar)
      endif

   end do
end do
!$OMP END DO

!$OMP DO
do j=1,ny
   do i=1,nx

      ! the actual streams and the warm based points are gzmx now:
      if ( ((.not.ilemy(i,j)).and.(fleuvemy(i,j))) .or. ((.not.ilemy(i,j)).and.(slowssamy(i,j))) ) gzmy(i,j)=.true. 

      
! calcul du frottement basal (ce bloc etait avant dans neffect)

      if (cotemy(i,j)) then                        ! point cotier
         betamy(i,j)=cf*neffmy(i,j) 
         betamy(i,j)=min(betamy(i,j),betamax)

      else if ((gzmy(i,j)).and.(.not.cotemy(i,j))) then ! tous les points SSA

         if (fleuvemy(i,j)) then                 ! the actual streams
            betamy(i,j)=cf*neffmy(i,j)*coef_gz
            betamy(i,j)=min(betamy(i,j),betamax)
            betamy(i,j)=max(betamy(i,j),betamin)

            if (cf*neffmy(i,j).gt.betamax*2.) then  ! a stream that's becoming grounded...
               if (slowssamy(i,j)) then
                  betamy(i,j)=betamax_2d(i,j)
               else
                  gzmy(i,j)=.false.
                  betamy(i,j)=tostick
               endif
            endif

         else                    ! not an actual stream, SSA is used to compute Ub
            betamy(i,j)=betamax_2d(i,j)
         endif
         
      else if (ilemy(i,j)) then
         betamy(i,j)=cf*neffmy(i,j)*coef_ile 
         betamy(i,j)=min(betamy(i,j),tob_ile)       
      else if (flotmy(i,j)) then     ! flottant ou ile
         betamy(i,j)=0.
      else                                         ! grounded, SIA
         betamy(i,j)=tostick                       ! frottement glace posee (1 bar)
      endif

   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

!~ where (fleuve(:,:))
!~    debug_3D(:,:,1)=1
!~ elsewhere (flot(:,:))
!~    debug_3D(:,:,1)=2
!~ elsewhere
!~    debug_3D(:,:,1)=0
!~ endwhere
!~ 
!~ where (cote(:,:))
!~    debug_3D(:,:,2)=1
!~ elsewhere
!~    debug_3D(:,:,2)=0
!~ endwhere
!~ 
!~ where (fleuvemx(:,:))
!~    debug_3D(:,:,3)=1
!~ elsewhere
!~    debug_3D(:,:,3)=0
!~ endwhere
!~ 
!~ where (flotmx(:,:))
!~    debug_3D(:,:,3)=-1
!~ endwhere

!_____________________


!~ where (fleuvemy(:,:))
!~    debug_3D(:,:,4)=1
!~ elsewhere
!~    debug_3D(:,:,4)=0
!~ endwhere
!~ 
!~ where (flotmy(:,:))
!~    debug_3D(:,:,4)=-1
!~ end where
!~ 
!~ debug_3D(:,:,23)= abs(RO*G*HMX(:,:)*sdx(:,:)*1.e-5) 
!~ debug_3D(:,:,24)= abs(RO*G*HMY(:,:)*sdy(:,:)*1.e-5) 

return
end subroutine dragging

end module dragging_neff_slope