source: trunk/SOURCES/Antarctique_general_files/dragging_plastic_LGM_mod.f90 @ 150

!> \file dragging-vit_bil_CISM_gen_neem_mod.f90
!! Definition des zones de streams apres lecture d'une carte evaluee offline
!<

!> \namespace dragging_neem
!! Definition zone stream. Input neem : critere sur vitesse de bilan et sur la 'courbure'
!! \author ...
!! \date ...
!! @note Defini les zones de stream avec :
!! @note   - un critere sur la hauteur d'eau 
!! @note   - un critere de fleuves actuels
!! @note Attention il faut fonctionner en additionnant SIA et Stream pour éviter les pb
!! de fleuves en pointille
!! @note Used module
!! @note   - use module3D_phy
!! @note   - use param_phy_mod
!<


module dragging_plastic_LGM

! Defini les zones de stream avec :
! * un critere sur la hauteur d'eau 
! * un critere de fleuves actuels
! attention il faut fonctionner en additionnant SIA et Stream pour éviter les pb
! de fleuves en pointille


use module3d_phy
use param_phy_mod
use interface_input


implicit none
logical,dimension(nx,ny) :: fleuvemx       !< fleuves sont les tableaux courants (dep. time)   
logical,dimension(nx,ny) :: fleuvemy
logical,dimension(nx,ny) :: fleuve      
logical,dimension(nx,ny) :: cote    
real, dimension(nx,ny) :: Vcol_x           !< uniquement pour compatibilite dans le spinup cat (aurel)
real, dimension(nx,ny) :: Vcol_y           !< uniquement pour compatibilite dans le spinup cat (aurel)
real, dimension(nx,ny) :: Vsl_x            !< pour compatibilite dragging beta (aurel)
real, dimension(nx,ny) :: Vsl_y            !< pour compatibilite dragging beta (aurel)
real, dimension(nx,ny) :: work_tab         ! pour lire les masques
real, dimension(nx,ny) :: drag_centre      ! 
real, dimension(nx,ny) :: uslid_centre     ! pour calcul beta vitesse de glissement du pas précédent
real, dimension(nx,ny) :: uslid_old        ! vitesse de glissement du pas antepenultieme pour test convergence

logical :: corr_def = .false.              !< for deformation correction, pour compatibilite beta (aurel)

character(len=100) :: masque_stream        !<  file masque streams autorises
character(len=100) :: file_ncdf            !< fichier netcdf issue des fichiers .dat

real    :: Tstream_lim                     !< limit for streaming (lower than melting point to account for spatial variability)
real    :: tob_prescrit                    !< valeur pour les ice streams
real    :: valmax
integer :: imax,jmax
integer :: i_moins1,i_plus1,j_moins1,j_plus1
integer :: lmax=20
integer :: idep,jdep,iloc,jloc

integer :: nestim                          !< pour les estimateurs de convergence
integer :: mm
real    :: uslid_diff
real    :: uslid_absdiff
real    :: uslid_moy


real :: tostick      ! pour la glace posee
real :: tob_ile      ! pour les iles

contains
!-------------------------------------------------------------------------------


!> SUBROUTINE: init_dragging
!! Cette routine fait l'initialisation du dragging.
!<   

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

implicit none

integer :: bidon

 
namelist/drag_plastic_LGM/Tstream_lim,tob_prescrit,betamax,masque_stream

if (itracebug.eq.1) write(num_tracebug,*)' init_dragging avec dragging_plastic_LGM'


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

! lecture des parametres du run                    
!--------------------------------------------------------------------

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

write(num_rep_42,428) '!___________________________________________________________' 
write(num_rep_42,428) '! module dragging_plasticLGM '
write(num_rep_42,drag_plastic_LGM)                    
write(num_rep_42,428) '! critere de temperature pour passage en stream'
write(num_rep_42,428) '! si T > Tlim'
write(num_rep_42,428) '! seulement pour les points cotiers'
write(num_rep_42,428) '! tob_prescrit : (Pa) frottement maxi sous les streams '
write(num_rep_42,428) '! masque_stream : nom du fichier qui contient le masque de streams autorises'
write(num_rep_42,*)

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

!-------------------------------------------------------------------
! masque stream : lecture du masque

masque_stream  = trim(dirnameinp)//trim(masque_stream)

!call lect_input(1,'z',1,work_tab,masque_stream,file_ncdf)

call lect_input(1,'z',1,work_tab,masque_stream,file_ncdf)

!call lect_ncfile('cry',work_tab,masque_stream)
mstream(:,:) = int(work_tab(:,:))

!i = 154
!j = 123
!if (itracebug.eq.1) write(num_tracebug,*) 'init_dragging', masque_stream, work_tab(i,j)

! elargit les streams autorisés aux noeuds mineurs adjacents
do j=1,ny
   do i=1,nx
      if (mstream(i,j).eq.1) then ! allowed ice streams
         mstream_mx(i:i+1,j) = 1
         mstream_my(i,j:j+1) = 1
      end if
   end do
end do


debug_3D(:,:,25) = mstream(:,:)


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

do j=2,ny-1
   do i=2,nx-1
      uslid_centre(i,j) = (ux(i,j,nz)+ux(i-1,j,nz))**2+(uy(i,j,nz)+uy(i,j-1,nz))**2
      uslid_centre(i,j) = uslid_centre(i,j)**0.5
      uslid_centre(i,j) = max(uslid_centre(i,j),0.1)
   end do
end do

uslid_old(:,:) = uslid_centre(:,:)

!!$call diffusiv()
!!$call  mix_SIA_L1  
!!$call strain_rate()
!!$ 
!!$do mm=1,50    
!!$   write(6,*)'time, m',time, mm
!!$   call diagnoshelf
!!$   call  mix_SIA_L1  
!!$   call strain_rate()
!!$end do


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

fleuvemx(:,:)=.false.
fleuvemy(:,:)=.false.
fleuve(:,:)=.false.
cote(:,:)=.false.
cotemx(:,:)=.false.
cotemy(:,:)=.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.
         cotemx(i:i+1,j) = .true.
         cotemy(i,j:j+1) = .true.
         gzmx(i:i+1,j)=.true.
         gzmy(i,j:j+1)=.true.
      endif
   end do
end do

! le calcul des fleuves se fait sur les mailles majeures
! normalement, les points cotiers sont deja tagges gzmx dans la routine flottab

do j=2,ny-1
   do i=2,nx-1
      uslid_centre(i,j) = (ux(i,j,nz)+ux(i-1,j,nz))**2+(uy(i,j,nz)+uy(i,j-1,nz))**2
      uslid_centre(i,j) = uslid_centre(i,j)**0.5
      uslid_centre(i,j) = max(uslid_centre(i,j),0.1)
   end do
end do

! calcul des estimateurs


nestim        = 0
uslid_diff    = 0.
uslid_absdiff = 0.
uslid_moy     = 0.


do j=2,ny-1
   do i=2,nx-1
      if ((uslid_centre(i,j).gt..2).and.(.not.flot(i,j))) then 
         nestim = nestim +1
         uslid_diff = uslid_diff + (uslid_centre(i,j)-uslid_old(i,j))
         uslid_absdiff = uslid_absdiff + abs( (uslid_centre(i,j)-uslid_old(i,j)))
         uslid_moy = uslid_moy + uslid_centre(i,j)
      end if
   end do
end do

if (nestim.gt.0) then
   write(6,*) 'conv uslid',nestim,uslid_diff/nestim,uslid_absdiff/nestim,uslid_moy/nestim 
end if

uslid_old(:,:) = uslid_centre(:,:)

do j=2,ny-1
   do i=2,nx-1
      if ( (T(i,j,nz)-Tpmp(i,j,nz).gt.Tstream_lim).and.(mstream(i,j).eq.1).and.(.not.flot(i,j))) then
         fleuve(i,j)      = .true.
         drag_centre(i,j) = tob_prescrit/uslid_centre(i,j)
         drag_centre(i,j) = max(drag_centre(i,j),10.)
         fleuvemx(i:i+1,j)=.true.
         fleuvemy(i,j:j+1)=.true.
         gzmx(i:i+1,j)=.true.
         gzmy(i,j:j+1)=.true.
      else
         drag_centre(i,j) = tostick
      endif
   end do
end do

!i = 154
!j = 123
!write(num_tracebug,*)'dans drag_plastic', mstream(i,j),drag_centre(i,j),tob_prescrit,uslid_centre(i,j)

where (flot(:,:))
   drag_centre(:,:) = 0.
end where

call  beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre)         ! redistribute on staggered grid



where (cotemx(:,:))               ! point cotier forcement beta faible
   drag_mx(:,:)=min(drag_mx(:,:),1000.)
end where

where (cotemy(:,:))               ! point cotier forcement beta faible
   drag_my(:,:)=min(drag_my(:,:),1000.)
end where

betamx(:,:)     = drag_mx(:,:)
betamy(:,:)     = drag_my(:,:)
beta_centre(:,:) = drag_centre(:,:)

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

!______________________________________________________________________________________
!>SUBROUTINE:   beta_c2stag_min  
!! redistribute central values on staggered grid
!<

  subroutine beta_c2stag(nxx,nyy,R_mx,R_my,R_centre)              ! redistribute central values on staggered grid

    implicit none
    integer                                ::  nxx,nyy                ! dimension des tableaux
    real, dimension(nxx,nyy), intent(out)  ::  R_mx                   ! valeur du beta sur maille mx
    real, dimension(nxx,nyy), intent(out)  ::  R_my                   ! valeur du beta sur maille my
    real, dimension(nxx,nyy), intent(in)   ::  R_centre               ! valeur du beta sur maille centree

    do j=2,nyy 
       do i=2,nxx
          R_mx(i,j)= (R_centre(i-1,j)+R_centre(i,j))*0.5
          R_my(i,j)= (R_centre(i,j-1)+R_centre(i,j))*0.5
       end do
    end do

  end subroutine beta_c2stag


end module dragging_plastic_LGM