source: branches/iLoveclim/SOURCES/dragging_hwat_contmaj_mod.f90 @ 244

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

!> \namespace  dragging_hwat_contmaj
!! Defini les zones de stream avec 3 criteres
!! \author ...
!! \date ...
!! @note Les trois criteres sont:
!! @note   * un critere sur la hauteur d'eau 
!! @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_hwat_contmaj

! Defini les zones de stream avec :
! * un critere sur la hauteur d'eau 
! * 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 fake_beta_iter_vitbil_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 :: tob_ile    ! pour les iles
real :: cry_lim=50.  ! courbure limite pour le suivi des fleuves


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_hwat_cont/hwatstream,cf,betamax,toblim,tostick    

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


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

! lecture des parametres du run                      block draghwat
!--------------------------------------------------------------------

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

write(num_rep_42,428)'!___________________________________________________________' 
write(num_rep_42,428) '&drag_hwat_cont              ! nom du bloc hwatermax '
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,*) 'toblim         = ', toblim
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) '! seulement pour les points cotiers'
write(num_rep_42,428) '! betamax : (Pa) frottement maxi sous les streams '
write(num_rep_42,428) '! toblim : (Pa)  pour les iles '
write(num_rep_42,*)

inv_beta=0
tostick=1.e5   ! valeurs pour les points non flgzmx
tob_ile=betamax/2.
moteurmax=toblim
betamax_2d(:,:)=betamax

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

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 (hwater(i,j).gt.hwatstream) then
         if (socle_cry(i,j).lt.0.) then    ! dans une vallee
            fleuve(i,j)=.true. 
!           write(166,*)
!           write(166,*)'cote',i,j,hwater(i,j),hwatstream
         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 la hauteur d'eau maxi chez les voisins
!!$! * en excluant les points flottants
!!$! * et ceux qui sont deja tagges fleuves   !version 18 juillet 2007

!!$           valmax=0.
!!$
!!$         do jloc=j_moins1,j_plus1
!!$            do iloc=i_moins1,i_plus1
!!$ 
!!$               if ((hwater(iloc,jloc).gt.valmax) & 
!!$                    .and.(.not.flot(iloc,jloc))  &
!!$                    .and.(.not.fleuve(iloc,jloc))) then
!!$                  imax=iloc
!!$                  jmax=jloc
!!$                  valmax=hwater(iloc,jloc)
!!$               endif
!!$            end do
!!$         end do

! 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.max(hwatstream-L*40.,200.)) then  ! version 18 juil
         if ((hwater(imax,jmax).gt.hwatstream).and.(socle_cry(i,j).lt.cry_lim)) then
            fleuve(imax,jmax)=.true.
            idep=imax
            jdep=jmax
!            write(166,*) imax,jmax,hwater(imax,jmax)
         else
            fleuve(imax,jmax)=.false.
            exit suit
         end if
 
      end do suit

   end if cote_detect

end do ifleuve
end do fleuve_maj


! detection de problemes
!write(166,*)
!write(166,*) 'detection de pb dans fleuve   time=',time

do j=2,ny
   do i=2,nx
      if (fleuve(i,j).and.(socle_cry(i,j).gt.cry_lim)) then
         write(166,*) 'pb fleuve',time,i,j,hwater(i,j),socle_cry(i,j)
      endif
   end do
end 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

! 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)=betamax
         betamx(i,j)=cf*neffmx(i,j)*20.                  ! attention a cette modif juillet 2007
         betamx(i,j)=min(betamx(i,j),betamax)
         betamx(i,j)=max(betamx(i,j),10.)

         if (cf*neffmx(i,j).gt.1500.) then
            gzmx(i,j)=.false.
            betamx(i,j)=tostick
         endif


      else if (ilemx(i,j)) then
         betamx(i,j)=cf*neffmx(i,j) 
         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)=betamax
         betamy(i,j)=cf*neffmy(i,j)*20. 
         betamy(i,j)=min(betamy(i,j),betamax)
         betamy(i,j)=max(betamy(i,j),10.)

         if (cf*neffmy(i,j).gt.1500.) then
            gzmy(i,j)=.false.
            betamy(i,j)=tostick
         endif

      else if (ilemy(i,j)) then
         betamy(i,j)=cf*neffmy(i,j) 
         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_hwat_contmaj