source: trunk/SOURCES/dragging_neff_contmaj_mod.f90 @ 56

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

! 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

implicit none
logical,dimension(nx,ny) :: fleuvemx
logical,dimension(nx,ny) :: fleuvemy
logical,dimension(nx,ny) :: fleuve
logical,dimension(nx,ny) :: cote


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
 
namelist/drag_neff_cont/hwatstream,cf,betamax,betamin,toblim,seuil_neff,coef_gz,coef_ile

if (itracebug.eq.1)  call tracebug(' dragging avec hwatermax')


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

write(num_rep_42,428)'!___________________________________________________________' 
write(num_rep_42,428) '&drag_neff_cont              ! nom du bloc dragging neff cont'
write(num_rep_42,*)
write(num_rep_42,*) 'hwatstream      = ',hwatstream 
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,*) 'seuil_neff      = ', seuil_neff
write(num_rep_42,*) 'coef_gz         = ', coef_gz
write(num_rep_42,*) 'coef_ile        = ', coef_ile
write(num_rep_42,*)'/'                            
write(num_rep_42,428) '! hwatstream (m) :  critere de passage en stream en partant de la cote'
write(num_rep_42,428) '!  si hwater > hwatstream '
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) '! 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,*)

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

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


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


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

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

fleuvemx(:,:)=.false.
fleuvemy(:,:)=.false.
fleuve(:,:)=.false.
cote(:,:)=.false.

! detection des cotes
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

! calcul de neff (pression effective sur noeuds majeurs)
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
!aurel, for the borders:
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.

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

! pas de fleuve dans les endroits trop en pente

!fleuvemx(:,:)=fleuvemx(:,:).and.(abs(rog*Hmx(:,:)*sdx(:,:)).lt.toblim)
!fleuvemy(:,:)=fleuvemy(:,:).and.(abs(rog*Hmy(:,:)*sdy(:,:)).lt.toblim)

do j=1,ny
   do i=1,nx
      if ((.not.ilemx(i,j)).and.(fleuvemx(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  ! stream
         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)

         !aurel, function of betamax and not 1500 as before:
         if (cf*neffmx(i,j).gt.betamax*2.) then
            gzmx(i,j)=.false.
            betamx(i,j)=tostick
         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


do j=1,ny                  ! le noeud 1 n'est pas attribue
   do i=1,nx   

      if ((.not.ilemy(i,j)).and.(fleuvemy(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  ! stream
         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)

         !aurel, function of betamax and not 1500 as before:
         if (cf*neffmy(i,j).gt.betamax*2.) then
            gzmy(i,j)=.false.
            betamy(i,j)=tostick
         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
      endif
      
   end do
end do


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_contmaj