source: trunk/SOURCES/dragging_jorge_mod.f90 @ 34

module dragging_jorge

!jalv: module inspire par le dragging_hwat_contmaj mais simplifie
! Definie les zones de stream avec uniquement:
! * un critere sur la hauteur d'eau
! * presence des sediments 
 

  use module3d_phy
!  use sedim_declar

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

real :: seuil_sedim              ! seuil sur hwater pour avoir du glissement
real :: valmax
integer :: imax,jmax
integer :: i_moins1,i_plus1,j_moins1,j_plus1
integer :: lmax=20
integer :: idep,jdep,iloc,jloc

integer :: itestd

real,dimension(nx,ny) :: betamx_temp      ! pour lissage du betamx
real,dimension(nx,ny) :: betamy_temp      ! pour lissage du betamy

real,dimension(nx,ny) :: mksedim      ! jalv, hérité de HUdson

real,dimension(nx,ny) :: h_sedim      ! epaisseur de sediment


real :: toto
real :: neffstream  !jalv valeur en dessous de la quelle onpasse en stream
real :: tostick   ! pour la glace posee
real :: tob_ile    ! pour les iles
real :: cry_lim=50.  ! courbure limite pour le suivi des fleuves
contains
!-------------------------------------------------------------------------------
subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.

implicit none
 
namelist/drag_jorge/hwatstream,cf,tobmax,toblim,seuil_sedim    

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


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

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

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

write(num_rep_42,428)'!___________________________________________________________' 
write(num_rep_42,428) '&drag_jorge              ! nom du bloc drag_jorge '
write(num_rep_42,*)
write(num_rep_42,*) 'hwatstream      = ',hwatstream 
write(num_rep_42,*) 'cf             = ',cf
write(num_rep_42,*) 'tobmax         = ', tobmax
write(num_rep_42,*) 'toblim         = ', toblim
write(num_rep_42,*) 'seuil_sedim         = ', seuil_sedim
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) '! tobmax : (Pa) frottement maxi sous les streams '
write(num_rep_42,428) '! toblim : (Pa)  pour les iles '
write(num_rep_42,*)

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


neffstream=1.e6  !jalv

!-------------------------------------------------------------------
! masque stream

      mstream_mx(:,:)=1
      mstream_my(:,:)=1

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





! lecture du fichier sediment :

open(20,file=TRIM(DIRNAMEINP)//'Sediment_map_Gabi_Laske/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.1)  ! pas d'epaisseurs de sediments negatives               


where (h_sedim(:,:).gt.250) 
      mksedim(:,:)=2
endwhere

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 gzmx initial a ete calculé dans flottab pour les points a cheval sur la grounding line


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


! critere de passage en stream
do  j=2,ny-1  
   do i=2,nx-1  
      if (((hwater(i,j).gt.hwatstream).or.((mksedim(i,j).eq.2).and.(hwater(i,j).ge.seuil_sedim))).and.(.not.flot(i,j))) then

!IF ((neffmx(i,j).LT.neffstream).AND.(neffmy(i,j).LT.neffstream).and.(.not.flot(i,j))) THEN
!        IF ((((neffmx(i,j).LT.neffstream).AND.(neffmy(i,j).LT.neffstream)).or.((mksedim(i,j).eq.2.).and.(hwater(i,j).gt.seuil_sedim))).AND.(.NOT.flot(i,j))) THEN

         fleuve(i,j)=.true.
      else
         fleuve(i,j)=.false.
      endif
   enddo
enddo

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)

!.....................

!jalv
!call detect_assym(nx,ny,0,76,1,0,1,0,betamx,itestd)
!if (itestd.gt.0) then
!   write(6,*) 'avant calcul betax asymetrie sur betamx  pour time=',time
!   stop
!else
!   write(6,*) 'avant calcul betax pas d asymetrie sur betamx  pour time=',time
!end if

!call detect_assym(nx,ny,0,76,1,0,1,0,betamy,itestd)
!if (itestd.gt.0) then
!    write(6,*) 'avant calcul betay asymetrie sur betamy  pour time=',time
!    stop
!else
!    write(6,*) 'avant calcul betay pas d asymetrie sur betamy  pour time=',time
!end if



!...........................................!   
! frottement selon x (avant dans neffect)   !    ! test jalv: tous les cas ont la meme loi pour beta: beta=cf*neff
!...........................................!    ! mais le sueil maxi est different pour chaque cas
 


DO j=1,ny
  DO i=1,nx

    IF ((.NOT.ilemx(i,j)).AND.(fleuvemx(i,j))) gzmx(i,j)=.TRUE. 
    
    IF (cotemx(i,j)) THEN                        ! point cotier: frott maxi=80, mini=8

!        betamx(i,j)=cf*neffmx(i,j)*1.      ! clim3_39
!        betamx(i,j)=cf*neffmx(i,j)*1.      ! clim3_36
!        betamx(i,j)=cf*neffmx(i,j)*1.      ! clim3_35
!         betamx(i,j)=cf*neffmx(i,j)*1.    ! clim3_22         
!         betamx(i,j)=cf*neffmx(i,j)*0.1  !clim3_21 
        betamx(i,j)=cf*neffmx(i,j)*0.1 ! standar
!         betamx(i,j)=cf*neffmx(i,j)*0.1/h_sedim(i,j) !clim3_20
!         betamx(i,j)=min(betamx(i,j),80.)
!         betamx(i,j)=max(betamx(i,j),8.)
    ELSE IF ((gzmx(i,j)).AND.(.NOT.cotemx(i,j))) THEN  ! stream normal: frott maxi=160, mini=16                  

!        betamx(i,j)=cf*neffmx(i,j)*50.      ! clim3_39
!         betamx(i,j)=cf*neffmx(i,j)*500. ! clim3_36
!         betamx(i,j)=cf*neffmx(i,j)*100. ! clim3_22
!         betamx(i,j)=cf*neffmx(i,j)*10.   ! clim3_21
        betamx(i,j)=cf*neffmx(i,j)*10. ! standar
!         betamx(i,j)=cf*neffmx(i,j)*10./h_sedim(i,j)
!         betamx(i,j)=min(betamx(i,j),160.)           
!         betamx(i,j)=max(betamx(i,j),16.)
         
        IF ((mksedim(i,j).EQ.2).AND.(hwater(i,j).GE.seuil_sedim)) THEN  ! jalv: stream  en zone sediment

!            betamx(i,j)=cf*neffmx(i,j)*5.      ! clim3_39
!            betamx(i,j)=cf*neffmx(i,j)*50. ! clim3_36
!            betamx(i,j)=cf*neffmx(i,j)*10. ! clim3_22
!            betamx(i,j)=cf*neffmx(i,j)*1. !clim3_21
            betamx(i,j)=cf*neffmx(i,j)*1. ! standar
!            betamx(i,j)=cf*neffmx(i,j)*1./h_sedim(i,j)
!            betamx(i,j)=min(betamx(i,j),120.)          !stream sediment: frott maxi=120, mini=12 
!           betamx(i,j)=max(betamx(i,j),12.)   
        ENDIF

    ELSE IF (ilemx(i,j)) THEN                    ! ile : frott maxi=40, mini=4

!        betamx(i,j)=cf*neffmx(i,j)*0.1      ! clim3_39
!         betamx(i,j)=cf*neffmx(i,j)*0.1   ! clim3_22
        betamx(i,j)=cf*neffmx(i,j)*0.1 ! standar
!         betamx(i,j)=cf*neffmx(i,j)*0.1  !clim3_21
!         betamx(i,j)=cf*neffmx(i,j)*0.1/h_sedim(i,j)
!         betamx(i,j)=min(betamx(i,j),40.)
!         betamx(i,j)=max(betamx(i,j),4.)
    ELSE IF (flotmx(i,j)) THEN                   ! flottant : frott 0.
        betamx(i,j)=0.
    ELSE                                         ! grounded, SIA
!         betamx(i,j)=tostick                      ! frottement glace posee (1 bar) divise par mil !jalv
!         betamx(i,j)=200.                          ! jalv: le cas restant glace posee a priori: frott=150
        betamx(i,j)=cf*neffmx(i,j)*1000.
!        betamx(i,j)=cf*neffmx(i,j)*1000./h_sedim(i,j)
    ENDIF

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

   end do
end do


!................................................!
!...test..jalv..lissage..des..beta...............!
!...on..calcule..la..moyenne..ponderee..des......!
!...points..aux..alentours.......................!


!do j=1,ny
!   do i=3,nx-2

!      betamx_temp(i,j)=(betamx(i,j)+(0.3*betamx(i-1,j))+(0.2*betamx(i-2,j))+(0.3*betamx(i+1,j))+(0.2*betamx(i+2,j)))/2
     
!   end do
!end do



!betamx(:,:)=betamx_temp(:,:)

 

!...........................................!
! frottement selon y (avant dans neffect)   !   ! test jalv: tous les cas ont la meme loi pour beta: beta=cf*neff
!...........................................!   ! mais le seuil maxi est different pour chaque cas


do j=1,ny                 
   do i=1,nx   

      if ((.not.ilemy(i,j)).and.(fleuvemy(i,j))) gzmy(i,j)=.true. 

      
      IF (cotemy(i,j)) THEN                        ! point cotier: froot maxi=40, mini=4

!          betamy(i,j)=cf*neffmy(i,j)*1.      ! clim3_39
!          betamy(i,j)=cf*neffmy(i,j)*1.      ! clim3_36        
!          betamy(i,j)=cf*neffmy(i,j)*1.      ! clim3_35
!          betamy(i,j)=cf*neffmy(i,j)*1.    ! clim3_22
!          betamy(i,j)=cf*neffmy(i,j)*0.1  !clim3_21
          betamy(i,j)=cf*neffmy(i,j)*0.1  ! standar
!         betamy(i,j)=cf*neffmy(i,j)*0.1/h_sedim(i,j) 
!         betamy(i,j)=min(betamy(i,j),80.)
!         betamy(i,j)=max(betamy(i,j),8.)
      else if ((gzmy(i,j)).and.(.not.cotemy(i,j))) then  ! stream normal: frott maxi=80, mini=8

!          betamy(i,j)=cf*neffmy(i,j)*50.      ! clim3_39
!          betamy(i,j)=cf*neffmy(i,j)*500. ! clim3_36
!          betamy(i,j)=cf*neffmy(i,j)*100. ! clim3_22
!          betamy(i,j)=cf*neffmy(i,j)*10. !clim3_21
          betamy(i,j)=cf*neffmy(i,j)*10. ! standar
!          betamy(i,j)=cf*neffmy(i,j)*10./h_sedim(i,j)
!          betamy(i,j)=min(betamy(i,j),160.)            
!          betamy(i,j)=max(betamy(i,j),16.)

         if ((mksedim(i,j).eq.2).and.(hwater(i,j).ge.seuil_sedim)) then              ! jalv: stream en zone sediment

!             betamy(i,j)=cf*neffmy(i,j)*5.      ! clim3_39
!            betamy(i,j)=cf*neffmy(i,j)*50. ! clim3_36
!            betamy(i,j)=cf*neffmy(i,j)*10. ! clim3_22
!             betamy(i,j)=cf*neffmy(i,j)*1. !clim3_21  
             betamy(i,j)=cf*neffmy(i,j)*1. ! standar 
!            betamy(i,j)=cf*neffmy(i,j)*1./h_sedim(i,j)
!            betamy(i,j)=min(betamy(i,j),120.)      !stream sediment: frott maxi=60, mini=0.
!            betamy(i,j)=max(betamy(i,j),12.)  
         endif
      
      else if (ilemy(i,j)) then                   ! ile : frott maxi=20, mini=2

!          betamy(i,j)=cf*neffmy(i,j)*0.1      ! clim3_39
!         betamy(i,j)=cf*neffmy(i,j)*0.1   ! clim3_22  
          betamy(i,j)=cf*neffmy(i,j)*0.1   ! standar
!         betamy(i,j)=cf*neffmy(i,j)*0.1  !clim3_21
!         betamy(i,j)=cf*neffmy(i,j)*0.1/h_sedim(i,j)
!         betamy(i,j)=min(betamy(i,j),40.)
!         betamy(i,j)=max(betamy(i,j),4.)
      else if (flotmy(i,j)) then                  ! flottant : frott 0
         betamy(i,j)=0.
      else                                        ! grounded, SIA
!         betamy(i,j)=tostick/1000                ! frottement glace posee divise par mil !jalv
!         betamy(i,j)=200.                         ! jalv: le cas restant glace posee a priori: frott maxi 200
         betamy(i,j)=cf*neffmy(i,j)*1000.
!         betamy(i,j)=cf*neffmy(i,j)*1000./h_sedim(i,j)
      endif

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

   end do
end do



!................................................!
!...test..jalv..lissage..des..beta...............!
!...on..calcule..la..moyenne..ponderee..des......!
!...points..aux..alentours.......................!


!do j=3,ny-2
!   do i=1,nx
      
!      betamy_temp(i,j)=(betamy(i,j)+(0.3*betamy(i,j-1))+(0.2*betamy(i,j-2))+(0.3*betamy(i,j+1))+(0.2*betamy(i,j+2)))/2
   
!   end do
!end do

!betamy(:,:)=betamy_temp(:,:)



!jalv
!call detect_assym(nx,ny,0,76,1,0,1,0,betamx,itestd)
!if (itestd.gt.0) then
!   write(6,*) 'apres dragging asymetrie sur betamx  pour time=',time
!   stop
!else
!   write(6,*) 'apres dragging pas d asymetrie sur betamx  pour time=',time
!end if

!call detect_assym(nx,ny,0,76,1,0,1,0,betamy,itestd)
!if (itestd.gt.0) then
!    write(6,*) 'apres dragging asymetrie sur betamy  pour time=',time
!    stop
!else
!    write(6,*) 'apres dragging pas d asymetrie sur betamy  pour time=',time
!end if





return
end subroutine dragging

end module dragging_jorge