source: branches/iLoveclim/SOURCES/dragging_hwat_sedim_mod.f90 @ 38

module dragging_hwat_sedim

! Defini les zones de stream avec :
! * un critere sur la hauteur d'eau 
! * un critere sur la courbure du socle (si negatif, vallees)
! * un critere sur l'epaisseur de sediments (carte de Gabi & Laske) 

use module3d_phy

implicit none
logical,dimension(nx,ny) :: fleuvemx
logical,dimension(nx,ny) :: fleuvemy
logical,dimension(nx,ny) :: fleuve
logical,dimension(nx,ny) :: fleuve_sedim
logical,dimension(nx,ny) :: fleuve_sedimx
logical,dimension(nx,ny) :: fleuve_sedimy
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) :: h_sedim    ! epaisseur de sediment
real,dimension(nx,ny) :: neff ! pression effective noeuds majeurs
real :: seuil_hwater  ! seuil sur hwater pour avoir du glissement en zone sediment
real :: seuil_sedim   ! seuil sur h_sedim pour avoir du glissement en zone sediment
real :: seuil_neff    ! seuil sur la pression effective pour avoir glissement en zone sediment
real :: coef_sedim    ! coef frottement zones sediments  (1)
real :: coef_gz       ! coef frottement zones stream std  (10)
real :: coef_ile      ! coef frottement zones iles   (0.1)
real :: toto

real :: betacotemax=100  ! valeur max pour betamx et betamy sur les cote
real :: betacotemin=4
real :: betasedmax=400    ! 400 dans c-HadCM3-evol-basal-melt-sed : cHadCM3s  et 250 dans c-HadCM3-evol-basal-melt-sed-betasedimin6
real :: betasedmin=10   ! 10 dans c-HadCM3-evol-basal-melt-sed : cHadCM3s et 6 dans c-HadCM3-evol-basal-melt-sed-betasedimin6
real :: betastreamax=500
real :: betastreamin=10
real :: betailemax=250
real :: betailemin=6

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

character(len=150) :: file_sedim   ! fichier epaisseur sediment

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

implicit none
 
namelist/drag_hwat_sedim/file_sedim,hwatstream,cf,betamax,toblim,seuil_hwater,seuil_sedim,seuil_neff,coef_sedim,coef_gz,coef_ile    

if (itracebug.eq.1) write(num_tracebug,*)' dragging avec hwatermax'

iter_beta=0

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

write(num_rep_42,428)'!___________________________________________________________' 
write(num_rep_42,428) '&drag_hwat_sedim              ! nom du bloc hwat_sedim '
write(num_rep_42,*)
write(num_rep_42,*) 'file_sedim     = ', trim(file_sedim)
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,*) 'seuil_hwater   = ', seuil_hwater
write(num_rep_42,*) 'seuil_sedim    = ', seuil_sedim
write(num_rep_42,*) 'seuil_neff     = ', seuil_neff
write(num_rep_42,*) 'coef_sedim     = ', coef_sedim
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) '! file_sedim : fichier carte epaisseur sediment'                            
write(num_rep_42,428) '! hwatstream (m) :  critere de passage en stream en partant de la cote 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,428) '! seuil_hwater (m) : seuil hwater pour avoir glissement sur zone sediment'
write(num_rep_42,428) '! seuil_sedim (m) : seuil epaisseur sediment pour avoir glissement sur zone sediment'
write(num_rep_42,428) '! seuil_neff (Pa) seuil sur la pression effective pour avoir glissement en zone sediment'
write(num_rep_42,428) '! coef_sedim : coef frottement zones stream sediments'
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

! lecture du fichier sediment :
!open(20,file=TRIM(DIRNAMEINP)//trim(file_sedim)) !   'sediment_ij_hemin40.dat')
!do j=1,ny
!   do i=1,nx
!      read(20,*),toto,toto,h_sedim(i,j)
!   enddo
!enddo
!h_sedim(:,:)=max(h_sedim(:,:),0.)  ! pas d'epaisseurs de sediments negatives
h_sedim(:,:)=10. 
  
      return
      end subroutine init_dragging
!________________________________________________________________________________


!-------------------------------------------------------------------------
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,j+1)+neffmy(i,j)+neffmy(i+1,j))/4
   enddo
enddo

! test sur hauteur d'eau et epaisseur sediment :
!where (h_sedim(:,:).GT.seuil_sedim.and.hwater(:,:).GT.seuil_hwater.and.(.not.flot(:,:))) fleuve(:,:)=.true.
where (h_sedim(:,:).GT.seuil_sedim.and.neff(:,:).LE.seuil_neff.and.hwater(:,:).GT.seuil_hwater.and.(.not.flot(:,:))) fleuve(:,:)=.true.

fleuve_sedim(:,:)=fleuve(:,:)  ! points uniquement fleuves sediment


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. 
!         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.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,ny
!      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
      if (fleuve_sedim(i,j)) then
         fleuve_sedimx(i,j)=.true.
         fleuve_sedimx(i+1,j)=.true.
         fleuve_sedimy(i,j)=.true.
         fleuve_sedimy(i,j+1)=.true.
      endif
   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)
!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)*0.1 
! cat         betamx(i,j)=min(betamx(i,j),tobmax)
! version Jorge
! point cotier: frott maxi=40, mini=4
!         betamx(i,j)=min(betamx(i,j),betacotemax)
!         betamx(i,j)=max(betamx(i,j),betacotemin)
! version Jorge stream  en zone sediment :
      else if (fleuve_sedimx(i,j).and.gzmx(i,j).and.(.not.cotemx(i,j))) then  ! stream  en zone sediment
! stream sediment: frott maxi=60, mini=6 
!cdc            betamx(i,j)=cf*neffmx(i,j)*3
             betamx(i,j)=cf*neffmx(i,j)*coef_sedim
!            betamx(i,j)=min(betamx(i,j),betasedmax)          
!            betamx(i,j)=max(betamx(i,j),betasedmin)
      else if ((gzmx(i,j)).and.(.not.cotemx(i,j))) then  ! stream normal
! cat         betamx(i,j)=tobmax
! cat         betamx(i,j)=cf*neffmx(i,j)*20.                  ! attention a cette modif juillet 2007
! cat         betamx(i,j)=min(betamx(i,j),tobmax)
! cat         betamx(i,j)=max(betamx(i,j),10.)
! version Jorge
! stream normal: frott maxi=80, mini=8
         betamx(i,j)=cf*neffmx(i,j)*coef_gz         
!         betamx(i,j)=min(betamx(i,j),betastreamax)           
!         betamx(i,j)=max(betamx(i,j),betastreamin)

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

      else if (ilemx(i,j)) then
         betamx(i,j)=cf*neffmx(i,j)*coef_ile
! cat         betamx(i,j)=min(betamx(i,j),tob_ile)    
! version Jorge
! ile : frott maxi=20, mini=2   
!         betamx(i,j)=min(betamx(i,j),betailemax)
!         betamx(i,j)=max(betamx(i,j),betailemin)
      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)*0.1 
! cat         betamy(i,j)=min(betamy(i,j),tobmax)
! version Jorge
! point cotier: frott maxi=40, mini=4
!         betamy(i,j)=min(betamy(i,j),betacotemax)
!         betamy(i,j)=max(betamy(i,j),betacotemin)
! version Jorge stream  en zone sediment :
      else if ((fleuve_sedimy(i,j).and.gzmx(i,j).and.(.not.cotemx(i,j)))) then  ! stream  en zone sediment
! stream sediment: frott maxi=60, mini=6 
!cdc            betamy(i,j)=cf*neffmy(i,j)*3
             betamy(i,j)=cf*neffmy(i,j)*coef_sedim
!            betamy(i,j)=min(betamy(i,j),betasedmax)          
!            betamy(i,j)=max(betamy(i,j),betasedmin)
      else if ((gzmy(i,j)).and.(.not.cotemy(i,j))) then  ! stream normal
! cat         betamy(i,j)=tobmax
! cat         betamy(i,j)=cf*neffmy(i,j)*20. 
! cat         betamy(i,j)=min(betamy(i,j),tobmax)
! cat         betamy(i,j)=max(betamy(i,j),10.)
! version Jorge
! stream normal: frott maxi=80, mini=8
         betamy(i,j)=cf*neffmy(i,j)*coef_gz                     
!         betamy(i,j)=min(betamy(i,j),betastreamax)           
!         betamy(i,j)=max(betamy(i,j),betastreamin)

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

      else if (ilemy(i,j)) then
         betamy(i,j)=cf*neffmy(i,j)*coef_ile
! cat         betamy(i,j)=min(betamy(i,j),tob_ile)     
! version Jorge
! ile : frott maxi=20, mini=2   
!         betamy(i,j)=min(betamy(i,j),betailemax)
!         betamy(i,j)=max(betamy(i,j),betailemin)
      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

!print*,'dragging debug h_sedim(78,144),hwater,seuil_sedim,seuil_hwater',h_sedim(78,144),seuil_sedim,hwater(78,144),seuil_hwater,fleuve(78,144),gzmx(78,144)
!print*,'dragging debug h_sedim(90,157),hwater,seuil_sedim,seuil_hwater',h_sedim(90,157),seuil_sedim,hwater(90,157),seuil_hwater,fleuve(90,157),gzmx(90,157),betamx(90,157),betamx(78,144)


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) 


!print*,'point (79,145)',hwater(79,145),h_sedim(79,145),flot(79,145),fleuve_sedim(79,145),fleuve(79,145),betamx(79,145),fleuve_sedimx(79,145),gzmx(79,145)

return
end subroutine dragging

end module dragging_hwat_sedim