!*************************************************************************************************************************************
! This script creates runoff directions from a topographic dataset based on the algorithm of Wang and Liu (2006) described in:
! An efficient method for identifying and filling surface depressions in digital elevation model for hydrologic analysis and modelling
! International Journal of Geographical Information Science, Vol. 20, No. 2, Feb. 2006, 193-213

!* It avoids endroheic basins and provide consistent runoff directions for use in the IPSL model.
!  It compares favourably with directions created with the 'Fill Sinks' algorithm implemented in 
!  GIS softwares (like ArcGIS or QGIS) which is also based on the Wang and Liu algorithm.
!  There are some differences in runoff directions related to the way I have interpreted some 
!  specific cases that occur. Basins generated by these directions are thus not identical but 
!  major basins are relatively similar.
!  The corrected topography is however rigourously similar to that obtained from GIS softwares, 
!  at the exception of the Antarctic continent where I noted some inconsistencies in these 
!  softwares due to lack of east-west periodicity and of understanding the increasing longitudinal resolution.

!* INPUTS:
!  - a topographic dataset

!* OUTPUTS:
!  - rtm: runoff directions
!  - topo: initial topographie
!  - orog_lake : corrected topographic dataset to avoid endorheic basin
!  - cmsk: coastal mask 

!*************************************************************************************************************************************

!* C Dumas - 26.11.2020

module lake_rsl

use module3D_phy,only:nx,ny,S,BSOC,H,sealevel,sealevel_2d,num_param,num_rep_42,time,dtmin,ice,mk,debug_3D,ro,row,dirnameinp,xcc,ycc,Bsoc0,S0,H0
use io_netcdf_grisli

!$ USE OMP_LIB

implicit none

!~ !* General & discrete parameters
!~ real :: undef=-9999. !undef=1.e20


!~ !* Parameters and arrays used in the Initialization step
real, dimension(nx,ny) :: orog, orog_lake, lake_level
integer, dimension(nx,ny) :: omsk, cmsk
integer, dimension(nx,ny) :: cont!, conttrack, contsave, runoff_msk, mask_search, exut
integer, dimension(nx,ny) :: ocean
integer, dimension(nx,ny) :: pot_lake
integer :: hi_res_runoff ! select 1 : interpolation on hi resolution topo or 0 : no interpolation

integer :: countcont ! nbr de continents
!~ integer :: nbocean, countocean

!~ !* Parameters and arrays of the runoff calculation
integer, dimension(nx,ny) :: rtm

!~ !* Debug mode & nc files
character(len=100) :: file_topo_hi_res, coord_topo_hi_res
!~ character(len=15) :: invar
!~ integer :: ncid, xdimid, ydimid, varid
!~ integer :: latid, lonid, rtmid, orogid, cmskid, contid, oroglakeid, lake_levelid

!* Routing code parameter
integer, dimension(8,2) :: inc

!~ !cdc variable ajoutees

!~ integer, dimension(8) :: loccont ! cont des 8 points autour du point ii jj

!~ integer :: sumcont ! somme des cont des 8 points autour du point ii jj


!
!~ integer :: maxloop ! pour debug evite boucle infinies
!~ integer,dimension(2) :: ijcoord ! coordonnees ij du point le plus bas
!~ integer,dimension(2) :: ijcoord_cmsk ! coordonnees ij du point cotier le plus bas
!~ real, dimension(8) :: locorog ! orog des 8 points voisins
integer,dimension(nx,ny) :: lake ! 1 si lac, 0 sinon
!~ integer,dimension(8) :: locexut ! 1 si exutoire => permet de gerer les lacs

!~ !integer :: contprint=0 ! pour debug et arret du print
!~ integer :: contexut ! compte les points sous le niveau du lac

!~ logical :: stop_lake
!~ integer, dimension(2) :: ij_lake
!~ integer,dimension(8) :: locomsk ! omsk des 8 voisins
!~ integer :: kloc

integer,dimension(80,2) :: inclake ! 80 voisins (+4 -4) en i,j
integer, dimension(nx,ny) :: mask_extlake ! extension du lac sous la glace
real, dimension(nx,ny) :: orog_extlake ! niveau du lac sur les zones glace a proximité (rayon de 4 pts de grille)
!~ integer, dimension(nx,ny) :: clake ! point lac avec voisin point englace
real, dimension(nx,ny) :: sealevel_2D_loc


!~ ! topo haute resolution (10km) pour calcul des lacs
!~ ! nbr de points topo hi res
integer, parameter :: nxhi = 964
integer, parameter :: nyhi = 964
integer, parameter :: ninterpo = 4 ! interpolation avec les 4 points autours

real, dimension(nxhi,nyhi) :: Bsoc0_hi, S0_hi, H0_hi ! topo haute resolution actuelle (ref)
!~ ! anomalie de topo sur la haute resolution
!~ !real, dimension(nxhi,nyhi) :: anom_topo_hires
!~ ! anomalie de topo
real,dimension(nx,ny) :: anom_Bsoc, anom_S, anom_H ! anomalie de topo vs topo ref (low res)

!~ ! matrice passage 40km vers 10km

real, dimension(nxhi,nyhi) :: xcc_hi,ycc_hi,xlong_hi,ylat_hi
real, dimension(nxhi,nyhi,4) :: poids  ! poids pour passage de low vers hi_res
integer,dimension(nxhi,nyhi,4) :: low_2_hi_i, low_2_hi_j ! correspondance grille low vers grille hi_res

real, dimension(nxhi,nyhi) :: anom_Bsoc_hi, anom_S_hi, anom_H_hi ! anomalie de topo resol std vs resol_hi
real, dimension(nxhi,nyhi) :: Bsoc_hi, S_hi, H_hi ! topo haute resolution a l'instant t
real, dimension(nxhi,nyhi) :: orog_hi ! orographie hi_res
real, dimension(nxhi,nyhi) :: S_interp_hi, H_interp_hi ! interpolation de S et H low res vers hi_res sans anomalie (pour faire topo englacee non bruitee)
integer, dimension(nxhi,nyhi) :: omsk_hi, cmsk_hi
integer, dimension(nxhi,nyhi) :: cont_hi ! masque des points continent hi_res
integer, dimension(nxhi,nyhi) :: ocean_hi ! masque des points ocean hi_res
integer, dimension(nxhi,nyhi) :: pot_lake_hi
integer, dimension(nxhi,nyhi) :: lake_hi ! 1 si lac, 0 sinon
integer, dimension(nxhi,nyhi) :: rtm_hi ! runoff direction (1 to 8)
real, dimension(nxhi,nyhi) :: lake_level_hi ! niveau du lac sur grille hi_res
real, dimension(nxhi,nyhi) :: orog_lake_hi  ! orographie a la surface des lacs

integer :: countcont_hi ! nbr de continents sur grille hi_res

contains

!> SUBROUTINE: input_rsl
!! Routine permettant d'introduire des variations locales du niveau marin
!! par la presence de lacs identifies via un algo de routage sur la topo de GRISLI
!>
subroutine input_rsl

  implicit none

  namelist/lake_rsl/hi_res_runoff, file_topo_hi_res, coord_topo_hi_res

  integer :: ii, i, j, k
  integer :: ios
  real*8, dimension(:,:), pointer   :: tab2d => null()   !< tableau 2d pointer
  real, dimension(ninterpo) :: dist_i, dist_j, dist_ij
  
  rewind(num_param)
  read(num_param,lake_rsl)

  write(num_rep_42,'(A)')'!___________________________________________________________'
  write(num_rep_42,'(A)') '&lake_rsl                                    ! nom du bloc '
  write(num_rep_42,*)
  write(num_rep_42,*) 'hi_res_runoff = ', hi_res_runoff
  write(num_rep_42,'(A,A,A)') 'file_topo_hi_res = ', trim(file_topo_hi_res)
  write(num_rep_42,'(A,A,A)') 'coord_topo_hi_res = ', trim(coord_topo_hi_res)
  write(num_rep_42,*)'/'
  write(num_rep_42,*)
  
!~   print*,'Debut input_rsl'
  
!* Continent number (outfile)
!ncfilecont = "continent_number_test-0.05.nc"
!ncfilecmsk = "coastal_mask_test-0.05.nc"
!* Runoff & corrected topo
!ncfilernff = "runoff_20Ma_Fred.nc"
!ncfiletopoout = "Topo.20Ma_Fred.OK_runoff.nc"
!ncfilernff = "runoff_NorthAm_15ka_test-0.05.nc"
!ncfiletopoout = "Topo.NorthAm_15ka_test.OK_runoff-0.05.nc"

!***************************
!*** Initialization step ***
!***************************

!~ write(*,*) "Start Initialization step"

!* The routing code (1er indice : k=1,8, 2eme indice i=1, j=2)
inc(1,1) = 0   ! indice i N
inc(1,2) = 1   ! indice j N
inc(2,1) = 1   ! indice i NE
inc(2,2) = 1   ! indice j NE
inc(3,1) = 1   ! indice i E
inc(3,2) = 0   ! indice j E
inc(4,1) = 1   ! indice i SE
inc(4,2) = -1  ! indice j SE
inc(5,1) = 0   ! indice i S
inc(5,2) = -1  ! indice j S
inc(6,1) = -1  ! indice i SW
inc(6,2) = -1  ! indice j SW
inc(7,1) = -1  ! indice i W
inc(7,2) = 0   ! indice j W
inc(8,1) = -1  ! indice i NW
inc(8,2) = 1   ! indice j NW

!~ print*
!~ print*,'rtm directions :'
!~ print*,'1 :  N'
!~ print*,'2 : NE'
!~ print*,'3 :  E'
!~ print*,'4 : SE'
!~ print*,'5 :  S'
!~ print*,'6 : SW'
!~ print*,'7 :  W'
!~ print*,'8 : NW'
!~ print*


!do ii=1,80
ii = 1
do j = -4,4
  do i = -4,4
    if (.not.((i.eq.0).and.(j.eq.0))) then
      inclake(ii,1)=i
      inclake(ii,2)=j
      ii=ii+1
    endif
  enddo
enddo


if (hi_res_runoff.eq.1) then
! lecture topo hi resolution :
  call Read_Ncdf_var('Bsoc',trim(dirnameinp)//trim(file_topo_hi_res),tab2d) ! socle
  Bsoc0_hi(:,:)=tab2d(:,:)
  call Read_Ncdf_var('S',trim(dirnameinp)//trim(file_topo_hi_res),tab2d)    ! surface
  S0_hi(:,:)=tab2d(:,:)
  call Read_Ncdf_var('H',trim(dirnameinp)//trim(file_topo_hi_res),tab2d)    ! epaisseur
  H0_hi(:,:)=tab2d(:,:)
  
! lecture des coordonnées geographiques
! les coordonnees sont calculees en °dec avec GMT,
! les longitudes sont comprises entre -180 et +180 (negative a l'Ouest de
! Greenwich et positive a l'Est)
  open(unit=2004,file=trim(dirnameinp)//trim(coord_topo_hi_res),iostat=ios)
  do k=1,nxhi*nyhi
     read(2004,*) i,j,xcc_hi(i,j),ycc_hi(i,j),xlong_hi(i,j),ylat_hi(i,j)
  enddo
  close(2004)
              
!  xmin=xcc(1,1)/1000.
!  ymin=ycc(1,1)/1000.
!  xmax=xcc(nx,ny)/1000.
!  ymax=ycc(nx,ny)/1000.


! calcul de la matrice de passage de 40km => hi_res (10km)
  do j=1,nyhi
    do i=1,nxhi
      low_2_hi_i(i,j,1) = max(nint((i-1)/4.),1)
      low_2_hi_i(i,j,2) = min(low_2_hi_i(i,j,1) + 1,nx)
      low_2_hi_i(i,j,3) = max(nint((i-1)/4.),1)
      low_2_hi_i(i,j,4) = min(low_2_hi_i(i,j,1) + 1,nx)
      low_2_hi_j(i,j,1) = max(nint((j-1)/4.),1)
      low_2_hi_j(i,j,2) = min(low_2_hi_j(i,j,1) + 1,ny)
      low_2_hi_j(i,j,3) = min(low_2_hi_j(i,j,1) + 1,ny)
      low_2_hi_j(i,j,4) = max(nint((j-1)/4.),1)
! calcul de la distance entre le pt hi_res et les 4 points low_res autours
      do k=1,ninterpo
        dist_i(k)=abs(xcc(low_2_hi_i(i,j,k),low_2_hi_j(i,j,k))-xcc_hi(i,j))
        dist_j(k)=abs(ycc(low_2_hi_i(i,j,k),low_2_hi_j(i,j,k))-ycc_hi(i,j))
        dist_ij(k)=sqrt(dist_i(k)**2+dist_j(k)**2)
      enddo
      do k=1,ninterpo
        poids(i,j,k)=(sum(dist_ij(:))-dist_ij(k))/(3*sum(dist_ij(:)))
      enddo
    enddo
  enddo

! calcul de l'anomalie de topo vs la topo de reference :
  anom_Bsoc(:,:) = Bsoc(:,:) - Bsoc0(:,:)
  anom_S(:,:) = S(:,:) - S0(:,:)
  anom_H(:,:) = H(:,:) - H0(:,:)

! interpolation 40km => hi_res
  anom_Bsoc_hi(:,:)=0.
  anom_S_hi(:,:)=0.
  anom_H_hi(:,:)=0.

  do k=1,ninterpo
    do j=1,nyhi
      do i=1,nxhi
        anom_Bsoc_hi(i,j) = anom_Bsoc_hi(i,j) + poids(i,j,k)*anom_Bsoc(low_2_hi_i(i,j,k),low_2_hi_j(i,j,k))
        anom_S_hi(i,j) = anom_S_hi(i,j) + poids(i,j,k)*anom_S(low_2_hi_i(i,j,k),low_2_hi_j(i,j,k))
        anom_H_hi(i,j) = anom_H_hi(i,j) + poids(i,j,k)*anom_H(low_2_hi_i(i,j,k),low_2_hi_j(i,j,k))
      enddo
    enddo
  enddo
endif

end subroutine input_rsl


!------------------------------------------------------------------------
! subroutine rsl : caclcul du niveau marin local
!------------------------------------------------------------------------
subroutine rsl

integer :: i, j, k, l
integer :: ii, jj
! pour les lectures ncdf
real*8, dimension(:,:),   pointer    :: tab => null()       !< tableau 2d real pointer

real, dimension(nx,ny) :: orog_lake_tmp ! pour calculer orog_lake a partir de orog_lake_hi
! pour debug : ecriture resultats hi_res
character(len=100) :: ncfilecont, ncfilernff, ncfiletopoout,ncfilecmsk    ! nom fichier netcdf de sortie
character(len=1),dimension(2) :: dimname ! nom des dimensions du fichier netcdf de sortie
integer :: ncid, xdimid, ydimid, S_hivarid, H_hivarid, Bsoc_hivarid, rtm_hivarid, cont_hivarid, ocean_hivarid, cmsk_hivarid, orog_lake_hivarid, lake_hivarid
integer :: status ! erreur operations netcdf

ncfilecont = "cont-topo_hi_res.nc"
dimname(1)='x'
dimname(2)='y'

!~ print*,'Debut rsl time = ', time
! orography sur grille std :


where ((BSOC(:,:)+H(:,:)*ro/row -sealevel).LT.0.)
  orog(:,:) = BSOC(:,:) ! le point flotte et socle sous le nivx de la mer
elsewhere
  orog(:,:) = BSOC(:,:) + H(:,:)
endwhere

orog_lake(:,:)=orog(:,:) ! orographie a la surface des lacs

where (orog(:,:).GT.sealevel) ! les depressions continentales seront ajoutees a ce masque plus loin
  omsk(:,:) = 1
elsewhere
  omsk(:,:) = 0
endwhere

!~ print*,'1/ rsl'

if (hi_res_runoff.eq.1) then
! calcul de l'anomalie de topo vs la topo de reference :
  anom_Bsoc(:,:) = Bsoc(:,:) - Bsoc0(:,:)
  anom_S(:,:) = S(:,:) - S0(:,:)
  anom_H(:,:) = H(:,:) - H0(:,:)
  
  anom_Bsoc_hi(:,:) = 0.
  anom_S_hi(:,:) = 0.
  S_interp_hi(:,:) = 0.
  H_interp_hi(:,:) = 0.
  anom_H_hi(:,:) = 0.
  
  ! calcul de la topo hi_res :
  do k=1,ninterpo
    do j=1,nyhi
      do i=1,nxhi
        anom_Bsoc_hi(i,j) = anom_Bsoc_hi(i,j) + poids(i,j,k)*anom_Bsoc(low_2_hi_i(i,j,k),low_2_hi_j(i,j,k))
        anom_S_hi(i,j) = anom_S_hi(i,j) + poids(i,j,k)*anom_S(low_2_hi_i(i,j,k),low_2_hi_j(i,j,k))
        S_interp_hi(i,j) = S_interp_hi(i,j) + poids(i,j,k)*S(low_2_hi_i(i,j,k),low_2_hi_j(i,j,k))
        H_interp_hi(i,j) = H_interp_hi(i,j) + poids(i,j,k)*H(low_2_hi_i(i,j,k),low_2_hi_j(i,j,k))
!        anom_H_hi(i,j) = anom_H_hi(i,j) + poids(i,j,k)*anom_H(low_2_hi_i(i,j,k),low_2_hi_j(i,j,k))
      enddo
    enddo
  enddo
  
  Bsoc_hi(:,:) = Bsoc0_hi(:,:) + anom_Bsoc_hi(:,:)
  
  where (H_interp_hi(:,:).GT.1.) ! zones englacees
    where (Bsoc_hi(:,:).GT.S_interp_hi(:,:)) ! socle qui sort de la glace
      S_hi(:,:) = Bsoc_hi(:,:)
      H_hi(:,:) = 0.
    elsewhere ! socle sous la glace
      S_hi(:,:) = S_interp_hi(:,:)
      H_hi(:,:) = S_interp_hi(:,:) - Bsoc_hi(:,:)
    endwhere
  elsewhere ! zones non englacees
    S_hi(:,:) = max(Bsoc_hi(:,:),sealevel)
    H_hi(:,:) = 0.
  endwhere
  
! orography sur grille hi_res :
  where ((Bsoc_hi(:,:)+H_hi(:,:)*ro/row -sealevel).LT.0.)
    orog_hi(:,:) = Bsoc_hi(:,:) ! le point flotte et socle sous le nivx de la mer
  elsewhere
    orog_hi(:,:) = S_hi(:,:)    ! point terre
  endwhere
  
  orog_lake_hi(:,:)=orog_hi(:,:) ! orographie a la surface des lacs

  where (orog_hi(:,:).GT.sealevel) ! les depressions continentales seront ajoutees a ce masque plus loin
    omsk_hi(:,:) = 1
  elsewhere
    omsk_hi(:,:) = 0
  endwhere

endif

!~ print*,'2/ rsl'
! mise à jour du niveau des lacs et du routage :
if (mod(abs(TIME),500.).lt.dtmin) then
  call mask_orog_ocean(nx,ny,orog,omsk,cmsk,cont,ocean,countcont,pot_lake)

!~   print*,'3/ rsl'
  if (hi_res_runoff.eq.1) then
    call mask_orog_ocean(nxhi,nyhi,orog_hi,omsk_hi,cmsk_hi,cont_hi,ocean_hi,countcont_hi,pot_lake_hi)
    call runoff(nxhi,nyhi,countcont_hi,cont_hi,cmsk_hi,omsk_hi,pot_lake_hi,lake_hi,rtm_hi,lake_level_hi,orog_lake_hi)
! interpolation des variables sur grille std :
    orog_lake(:,:) = 0.
    orog_lake_tmp(:,:) = 0.
    lake(:,:) = 0
    do j=1,ny
      do i=1,nx
        orog_lake_tmp(i,j) = max(sum(orog_lake_hi(i*4-3:i*4,j*4-3:j*4)) / 16., orog(i,j))
!        lake(i,j) = nint(sum(lake_hi(i*4-3:i*4,j*4-3:j*4)) / 16.)  ! lac si au moins 50% des points sont lacs
        lake(i,j) = floor(sum(lake_hi(i*4-3:i*4,j*4-3:j*4)) / 16.)  ! lac si 100% des points hi_res sont lacs
      enddo
    enddo
    where (lake(:,:).EQ.1) orog_lake(:,:)=orog_lake_tmp(:,:)
  
  else
    call runoff(nx,ny,countcont,cont,cmsk,omsk,pot_lake,lake,rtm,lake_level,orog_lake)
  endif

! pour eviter la prise en compte des lacs dans les depressions a la surface de la calotte :
  where ((lake(:,:).eq.1).and.(ice.eq.1))
    lake(:,:) = 0
    orog_lake(:,:) = orog(:,:)
  endwhere


!~   print*,'3/ rsl'
  call extlake(nx,ny,lake,ice,mk,orog_lake,inclake,mask_extlake,orog_extlake)

!~   print*,'4/ rsl'

  sealevel_2D_loc(:,:) = sealevel_2D(:,:) ! pour debug : A SUPPRIMER QUAND ON REPASSERA SUR SEALEVEL_2D
  where (lake(:,:).eq.1)
    sealevel_2D_loc(:,:)=orog_lake(:,:)
  elsewhere
    sealevel_2D_loc(:,:)=sealevel
  endwhere

  where ((mask_extlake.LT.9999).AND.(ice.EQ.1).AND.(mk.EQ.0).AND.(lake.EQ.0).AND.(BSOC.LT.orog_extlake))
    sealevel_2D_loc(:,:) = orog_extlake(:,:)
  endwhere

  if (hi_res_runoff.eq.1) then

!* Save mode
    status = nf90_create(trim(ncfilecont),NF90_CLOBBER, ncid)
!  status = nf90_close(ncid)
    status = nf90_def_dim(ncid, "x", nxhi, xdimid)
    status = nf90_def_dim(ncid, "y", nyhi, ydimid)
    status = nf90_def_var(ncid, "S_hi", nf90_float, (/ xdimid, ydimid /), S_hivarid)
    status = nf90_def_var(ncid, "H_hi", nf90_float, (/ xdimid, ydimid /), H_hivarid)
    status = nf90_def_var(ncid, "Bsoc_hi", nf90_float, (/ xdimid, ydimid /), Bsoc_hivarid)
    status = nf90_def_var(ncid, "rtm_hi", nf90_float, (/ xdimid, ydimid /), rtm_hivarid)
    status = nf90_def_var(ncid, "cont_hi", nf90_float, (/ xdimid, ydimid /), cont_hivarid)
    status = nf90_def_var(ncid, "ocean_hi", nf90_float, (/ xdimid, ydimid /), ocean_hivarid)
    status = nf90_def_var(ncid, "cmsk_hi", nf90_float, (/ xdimid, ydimid /), cmsk_hivarid)
    status = nf90_def_var(ncid, "orog_lake_hi", nf90_float, (/ xdimid, ydimid /), orog_lake_hivarid)
    status = nf90_def_var(ncid, "lake_hi", nf90_float, (/ xdimid, ydimid /), lake_hivarid)
  
    status = nf90_enddef(ncid)
  
    status = nf90_put_var(ncid, S_hivarid, S_hi)
    status = nf90_put_var(ncid, H_hivarid, H_hi)
    status = nf90_put_var(ncid, Bsoc_hivarid, Bsoc_hi)
    status = nf90_put_var(ncid, rtm_hivarid, rtm_hi)
    status = nf90_put_var(ncid, cont_hivarid, cont_hi)
    status = nf90_put_var(ncid, ocean_hivarid, ocean_hi)
    status = nf90_put_var(ncid, cmsk_hivarid, cmsk_hi)
    status = nf90_put_var(ncid, orog_lake_hivarid, orog_lake_hi)
    status = nf90_put_var(ncid, lake_hivarid, lake_hi)
  
    status = nf90_close(ncid)
  endif


endif

sealevel_2D(:,:) = sealevel_2D_loc(:,:)

debug_3d(:,:,126)=sealevel_2D_loc(:,:)
debug_3d(:,:,127)=rtm(:,:)
debug_3d(:,:,128)=lake(:,:)

!~ print*,'Fin rsl'

end subroutine rsl



!----------------------------------------------------------------------------
! Subroutine de recherche des points terre, cotiers et ocean
!----------------------------------------------------------------------------
subroutine mask_orog_ocean(nx_loc,ny_loc,orog_loc,omsk_loc,cmsk_loc,cont_loc,ocean_loc,countcont_loc,pot_lake_loc)

implicit none

integer,intent(in) :: nx_loc,ny_loc
real,dimension(nx_loc,ny_loc), intent(in) :: orog_loc
integer,dimension(nx_loc,ny_loc), intent(inout) :: omsk_loc
integer,dimension(nx_loc,ny_loc), intent(out) :: cmsk_loc, cont_loc
integer, dimension(nx_loc,ny_loc), intent(out) :: ocean_loc
integer, intent(out) :: countcont_loc
integer, dimension(nx_loc,ny_loc), intent(out) :: pot_lake_loc


integer :: undefint=0
integer, dimension(nx_loc,ny_loc) :: conttrack
integer :: sumcont ! somme des cont des 8 points autour du point ii jj
integer :: maxloop ! pour debug evite boucle infinies
integer, dimension(8) :: ipinctab, jpinctab ! indice des points a scanner autour du point ii jj
integer, dimension(8) :: loccont ! cont des 8 points autour du point ii jj
integer, dimension(nx_loc,ny_loc) :: contsave

integer :: nbocean, countocean
integer :: countpot_lake ! nbr de dépressions (potential lakes)
logical :: stop_lake
integer, dimension(2) :: ij_lake

integer,dimension(8) :: locomsk ! omsk des 8 voisins
integer :: maxcont ! valeur max de cont parmi les 8 voisins
integer :: sumocean ! somme des ocean des 8 points autour du point ii jj
integer, dimension(8) :: lococean ! ocean des 8 points autour du point ii jj
integer :: maxocean ! valeur max de ocean des 8 points autour du point ii jj
integer,dimension(nx_loc,ny_loc) :: oceansave, true_ocean
integer, parameter :: oceanmax=3000 ! nbr max de zones ocean ou lac
logical, dimension(oceanmax) :: k_ocean ! permet d'identifier les ocean ouverts des lacs
integer :: nbcont
integer :: kloc
integer :: i,j,ii,jj,k,l

!* Initialize the coastal mask (cmsk_loc)
cmsk_loc = 0
!~ print*,'1/ mask_orog_ocean'

!* Discriminate continents
!~ write(*,*) "Number the different continents"
cont_loc = undefint
nbcont = 1
conttrack = 0

do j = 1, ny_loc
  do i = 1, nx_loc
!do j = 110,111
!  do i = 182,183 

    !* Locate current grid point
    ii = i; jj = j
    sumcont=0
    maxloop=0

    !* Only on continental points
    if ( omsk_loc(ii,jj) .eq. 1) then
      do k=1,8
        ipinctab(k) = max(1,min(nx_loc,ii+inc(k,1)))
! cas cyclicité selon i :
!        ipinctab(k) = ii+inc(k,1)
!        if (ipinctab(k).eq.(nx + 1)) ipinctab(k) = 1
!        if (ipinctab(k).eq.(0)) ipinctab(k) = nx
        jpinctab(k) = max(1,min(ny_loc,jj+inc(k,2)))
        loccont(k)=cont_loc(ipinctab(k),jpinctab(k))
        sumcont = sumcont + cont_loc(ipinctab(k),jpinctab(k))
      enddo
 
      if (sumcont.GT.0) then
        cont_loc(ii,jj)=minval(loccont,mask=loccont>0)
        do while (maxval(loccont,mask=loccont>0).GT.cont_loc(ii,jj).and.maxloop.lt.8) ! test si dans les 8 voisins certains ont un loccont > minval
          maxcont=maxval(loccont,mask=loccont>0)
          where (cont_loc(:,:).EQ.maxcont) cont_loc(:,:) = cont_loc(ii,jj)
          where (loccont(:).EQ.maxcont) loccont(:) = cont_loc(ii,jj)
          maxloop = maxloop + 1
        enddo
      else
        cont_loc(ii,jj) = nbcont
        nbcont = nbcont + 1
      endif
    conttrack(ii,jj) = 1
    endif
  enddo
enddo
!~ print*,'2/ mask_orog_ocean'   

!* Reordering to get an incremental list
countcont_loc = 1
contsave(:,:)=cont_loc(:,:) ! copie cont for reordonation

if (nbcont.gt.0) then ! on a trouve des continents !!!!
  do k=1,nbcont
    if (count(contsave(:,:).eq.k).gt.0) then ! there is continents associated with number k
      where (cont_loc(:,:).eq.k) cont_loc(:,:) = countcont_loc
      countcont_loc = countcont_loc + 1
    endif
  enddo
else
  print*,'CONTINENT NOT FOUND'
endif
!~ print*,'3/ mask_orog_ocean'

! recherche des lacs et oceans sous le niveau de la mer
!~ write(*,*) "Number the different lakes & oceans"
ocean_loc = undefint
nbocean = 1

do j = 1, ny_loc
  do i = 1, nx_loc
    !* Locate current grid point
    ii = i; jj = j
    sumocean=0
    maxloop=0

    !* Only on ocean points
    if ( omsk_loc(ii,jj) .eq. 0) then  ! omask=0 => orog < sealevel
      do k=1,8
        ipinctab(k) = max(1,min(nx_loc,ii+inc(k,1)))
! cas cyclicité selon i :
!        ipinctab(k) = ii+inc(k,1)
!        if (ipinctab(k).eq.(nx + 1)) ipinctab(k) = 1
!        if (ipinctab(k).eq.(0)) ipinctab(k) = nx
        jpinctab(k) = max(1,min(ny_loc,jj+inc(k,2)))
        lococean(k)=ocean_loc(ipinctab(k),jpinctab(k))
        sumocean = sumocean + ocean_loc(ipinctab(k),jpinctab(k))
      enddo

      if (sumocean.GT.0) then
        ocean_loc(ii,jj)=minval(lococean,mask=lococean>0)
        do while (maxval(lococean,mask=lococean>0).GT.ocean_loc(ii,jj).and.maxloop.lt.8) ! test si dans les 8 voisins certains ont un lococean > minval
          maxocean=maxval(lococean,mask=lococean>0)
          where (ocean_loc(:,:).EQ.maxocean) ocean_loc(:,:) = ocean_loc(ii,jj)
          where (lococean(:).EQ.maxocean) lococean(:) = ocean_loc(ii,jj)
          maxloop = maxloop + 1
        enddo
      else
        ocean_loc(ii,jj) = nbocean
        nbocean = nbocean + 1
      endif
    endif
  enddo
enddo

!~ print*,'4/ mask_orog_ocean'  

!* Reordering to get an incremental list
countocean = 1
countpot_lake = 1 ! nombre de dépression intérieur (lacs potentiels
oceansave(:,:)=ocean_loc(:,:) ! copie ocean for reordonation
true_ocean(:,:)=0 ! point vraiment ocean
pot_lake_loc(:,:)=0

k_ocean(:)=.false.
 

if (nbocean.gt.0) then ! on a trouve des oceans !!!!
! si l'ocean a un point qui touche un des bords de grille c'est bien un ocean, sinon c'est un lac
  do i=1,nx_loc
    if (oceansave(i,1).GT.0) k_ocean(oceansave(i,1))=.true.
    if (oceansave(i,ny_loc).GT.0) k_ocean(oceansave(i,ny_loc))=.true.
  enddo
  do j=1,ny_loc
    if (oceansave(1,j).GT.0) k_ocean(oceansave(1,j))=.true.
    if (oceansave(nx_loc,j).GT.0) k_ocean(oceansave(nx_loc,j))=.true.
  enddo
  
  do k=1,nbocean
    if (count(oceansave(:,:).eq.k).gt.0) then ! there is oceans points associated with number k
      if (k_ocean(k)) then  ! point is true ocean
        where (oceansave(:,:).eq.k)
          ocean_loc(:,:) = countocean
        endwhere
        countocean = countocean + 1
      else ! point is a depresion under sea level
        where (oceansave(:,:).eq.k)
          pot_lake_loc(:,:) = countpot_lake
          ocean_loc(:,:) = undefint
        endwhere
        countpot_lake = countpot_lake + 1
      endif
    endif
  enddo
else
  print*,'OCEAN NOT FOUND'
endif
!~ print*,'5/ mask_orog_ocean countpot_lake', countpot_lake
! les points pot_lake sont continent => on leur attribue le no du continent adjacent.

do l=1,countpot_lake
  stop_lake = .true.
  ij_lake=minloc(pot_lake_loc(:,:),mask=(pot_lake_loc.EQ.l))
  do while (stop_lake)
    do k=1,8
      ipinctab(k) = max(1,min(nx_loc,ij_lake(1)+inc(k,1)))
! si cyclicite selon i :
!      ipinctab(k) = ijcoord(1)+inc(k,1)
!      if (ipinctab(k).eq.(nx + 1)) ipinctab(k) = 1
!      if (ipinctab(k).eq.(0)) ipinctab(k) = nx
      jpinctab(k) = max(1,min(ny_loc,ij_lake(2)+inc(k,2)))
      locomsk(k)=omsk_loc(ipinctab(k),jpinctab(k))
      loccont(k)=cont_loc(ipinctab(k),jpinctab(k))
    enddo
    if (maxval(loccont).GT.0) then ! voisin point terre
      ! le lac appartient a ce continent
      kloc=maxloc(loccont,dim=1)
      where (pot_lake_loc(:,:).EQ.l)
        cont_loc(:,:)=cont_loc(ipinctab(kloc),jpinctab(kloc))
        omsk_loc(:,:)=1
      endwhere
      stop_lake=.false.
    else
      ij_lake(1) = ij_lake(1)+1
    endif
  enddo
enddo
!~ print*,'6/ mask_orog_ocean'

! recherche des points cotiers
where ((omsk_loc(:,:).EQ.1).AND.((eoshift(omsk_loc,shift=-1,dim=1).EQ.0) .OR. &
  (eoshift(omsk_loc,shift=1,dim=1).EQ.0) .OR. (eoshift(omsk_loc,shift=-1,dim=2).EQ.0) .OR. (eoshift(omsk_loc,shift=1,dim=2).EQ.0)))
  cmsk_loc(:,:) = 1
endwhere
!~ print*,'7/ mask_orog_ocean'

end subroutine mask_orog_ocean



!---------------------------------------------------------
! subroutine qui calcul l'écoulement (runoff) et les lacs
!---------------------------------------------------------
subroutine runoff(nx_loc,ny_loc,countcont_loc,cont_loc,cmsk_loc,omsk_loc,pot_lake_loc,lake_loc,rtm_loc,lake_level_loc,orog_lake_loc)

implicit none

integer, intent(in) :: nx_loc, ny_loc
integer, intent(in) :: countcont_loc ! nbr de continents

integer,dimension(nx_loc,ny_loc), intent(in) :: cont_loc
integer,dimension(nx_loc,ny_loc), intent(in) :: cmsk_loc
integer,dimension(nx_loc,ny_loc), intent(in) :: omsk_loc
integer,dimension(nx_loc,ny_loc), intent(in) :: pot_lake_loc

integer,dimension(nx_loc,ny_loc), intent(out) :: lake_loc ! 1 si lac, 0 sinon
integer,dimension(nx_loc,ny_loc), intent(out) :: rtm_loc ! runoff direction (1 to 8)
real,dimension(nx_loc,ny_loc), intent(out) :: lake_level_loc ! lake level
real,dimension(nx_loc,ny_loc), intent(inout) :: orog_lake_loc  ! orographie a la surface des lacs

real,dimension(nx_loc,ny_loc) :: orog_extlake ! niveau du lac sur les zones glace a proximité (rayon de 4 pts de grille)

integer,dimension(nx_loc,ny_loc) :: exut ! > 0 si exutoire, 0 sinon

integer,dimension(nx_loc,ny_loc) :: mask_search
integer,dimension(nx_loc,ny_loc) :: runoff_msk
integer :: nbcont ! no du continent étudié
integer,dimension(2) :: ijcoord ! coordonnees ij du point le plus bas

integer,dimension(2) :: ijcoord_cmsk ! coordonnees ij du point cotier le plus bas
integer :: contexut ! compte les points sous le niveau du lac
real,dimension(8) :: locorog ! orog des 8 points voisins
integer,dimension(8) :: locexut ! 1 si exutoire => permet de gerer les lacs
integer,dimension(8) :: ipinctab,jpinctab ! indice des points a scanner autour du point ii jj
integer :: countpot_lake ! nbr de dépressions (potential lakes)


integer :: k,l

!**************************
!*** Runoff calculation ***
!**************************

write(*,*) "Start runoff calculation"


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! version inspiree article
!$OMP PARALLEL
!$OMP WORKSHARE
mask_search(:,:)=0 ! masque des points a prendre en compte a l'iteration n : 0 pas traité, 1 en cours, 2 traité
runoff_msk(:,:)=0  ! initialisation du mask points non traités => 0
lake_loc(:,:)=0
exut(:,:)=0 ! point exutoire d'un lac => 1
!output(:,:)=0 ! initialisation des points de sortie vers la mer (0 pas sortie, 1 pt de sortie)
!* We work on each continent separately for convenience
!do nbcont = 1,1
!$OMP END WORKSHARE
!$OMP END PARALLEL

do nbcont = 1,countcont_loc
  write(*,*) "Working on continent nb :", nbcont
! on demarre des points cotiers en partant du plus bas puis en remontant en altitude
! mise a jour du masque des points a etudier :

! demarrage par les points cotier du continent scanné :
  where ((cont_loc.eq.nbcont).and.(cmsk_loc.eq.1)) mask_search(:,:)=1 ! masque de scan : points cote pour demarrer
  
! debut de la boucle de recherche du point a scanner 
  do while (minval(runoff_msk(:,:),mask=mask_search(:,:).EQ.1).EQ.0) ! test si il reste des points a scanner
    ijcoord(:)=minloc(orog_lake_loc,mask=(mask_search(:,:).EQ.1))   ! position du point le plus bas etudié
    
    if (count(mask=((mask_search.eq.1).and.(cmsk_loc.eq.1))).GT.0) then! il reste des points cotiers non traités
    ! position du point cotier le plus bas :
      ijcoord_cmsk(:)=minloc(orog_lake_loc,mask=((mask_search(:,:).EQ.1).and.cmsk_loc(:,:).eq.1))
      ! si altitude egale entre 2 points on traite en priorité le point cotier
      if (orog_lake_loc(ijcoord(1),ijcoord(2)).EQ.orog_lake_loc(ijcoord_cmsk(1),ijcoord_cmsk(2))) then
        ijcoord(:) = ijcoord_cmsk(:)
      endif
    endif
    contexut = 0
    do k=1,8
      ipinctab(k) = max(1,min(nx_loc,ijcoord(1)+inc(k,1)))
! si cyclicite selon i :
!      ipinctab(k) = ijcoord(1)+inc(k,1)
!      if (ipinctab(k).eq.(nx + 1)) ipinctab(k) = 1
!      if (ipinctab(k).eq.(0)) ipinctab(k) = nx
      jpinctab(k) = max(1,min(ny_loc,ijcoord(2)+inc(k,2)))
      
      ! on ajoute les points adjacents aux points scannés pour le prochain tour sauf si ce sont des points ocean ou dejà traités:
      if ((omsk_loc(ipinctab(k),jpinctab(k)).EQ.1).AND.(mask_search(ipinctab(k),jpinctab(k)).EQ.0)) then
        mask_search(ipinctab(k),jpinctab(k))=1
! cas des depressions :
! si l'un des points est plus bas que le point de base et n'est pas cotier => on remonte son orog => orog_lake_loc
        if ((orog_lake_loc(ipinctab(k),jpinctab(k))).LE.(orog_lake_loc(ijcoord(1),ijcoord(2))) &
          .AND.(cmsk_loc(ipinctab(k),jpinctab(k)).NE.1)) then
          orog_lake_loc(ipinctab(k),jpinctab(k)) = orog_lake_loc(ijcoord(1),ijcoord(2))
          lake_loc(ipinctab(k),jpinctab(k)) = 1
          contexut= contexut + 1
        endif
      endif
      locorog(k) = orog_lake_loc(ipinctab(k),jpinctab(k))
      locexut(k) = exut(ipinctab(k),jpinctab(k))
    enddo
    if ((contexut.ge.1).and.(exut(ijcoord(1),ijcoord(2)).EQ.0)) exut(ijcoord(1),ijcoord(2)) = 1
    rtm_loc(ijcoord(1),ijcoord(2))=minloc(locorog,dim=1) ! sens ecoulement donne par minloc de la topo des 8 voisins
    ! cas d'un point lac
    if (lake_loc(ijcoord(1),ijcoord(2)).eq.1) then
      exut(ijcoord(1),ijcoord(2)) = maxval(locexut) + 1
    ! ecoulement vers l'exutoire du lac
      if (maxval(locexut).GT.0) then
        rtm_loc(ijcoord(1),ijcoord(2))=minloc(locexut,dim=1,mask=(locexut.ge.1)) ! ecoulement vers l'exutoire
      endif
    endif
    ! on supprime le point attribué de mask_search :
    mask_search(ijcoord(1),ijcoord(2))=2
    ! le point a ete traité :
    runoff_msk(ijcoord(1),ijcoord(2))=1
  enddo
enddo

! pour chaque lac calcul du niveau max :
lake_level_loc(:,:)=-9999. ! initialisation
do l=1,countpot_lake
  where (pot_lake_loc(:,:).EQ.l)
    lake_level_loc(:,:)=orog_lake_loc(:,:)
  endwhere
enddo

end subroutine runoff

!--------------------------------------------------------------------
!  subroutine extlake : extension des lacs sous la glace 
!--------------------------------------------------------------------

subroutine extlake(nx_loc,ny_loc,lake_loc,ice_loc,mk_loc,orog_lake_loc,inclake_loc,mask_extlake_loc,orog_extlake_loc)

  implicit none

  integer, intent(in) :: nx_loc, ny_loc
  integer,dimension(nx_loc,ny_loc), intent(in) :: lake_loc ! 1 si lac, 0 sinon
  integer,dimension(nx_loc,ny_loc), intent(in) :: ice_loc ! 1 si glace, 0 sinon
  integer,dimension(nx_loc,ny_loc), intent(in) :: mk_loc ! 0 si terre, 1 sinon
  real,dimension(nx_loc,ny_loc), intent(in) :: orog_lake_loc ! orographie a la surface des lacs
  integer,dimension(80,2),intent(in) :: inclake_loc ! 80 voisins (+4 -4) en i,j
  integer, dimension(nx_loc,ny_loc),intent(out) :: mask_extlake_loc ! extension du lac sous la glace
  real, dimension(nx_loc,ny_loc),intent(out) :: orog_extlake_loc ! niveau du lac sur les zones glace a proximité (rayon de 4 pts de grille)

  integer, dimension(8) :: ipinctab,jpinctab ! indice des points a scanner autour du point ii jj
  integer, dimension(nx_loc,ny_loc) :: clake ! point lac avec voisin point englace
  integer :: i,j,k

! initialisation
  clake(:,:)=0
  orog_extlake_loc(:,:)=0.
  mask_extlake_loc(:,:)=9999

! recherche des points lacs en bordure de calotte => voisin pas lac englace et terre
  where ((lake_loc(:,:).EQ.1).AND.(((eoshift(lake_loc,shift=-1,dim=1).EQ.0).AND.(eoshift(ice_loc,shift=-1,dim=1).EQ.1).AND.(eoshift(mk_loc,shift=-1,dim=1).EQ.0)) .OR. &
    ((eoshift(lake_loc,shift=1,dim=1).EQ.0).AND.(eoshift(ice_loc,shift=1,dim=1).EQ.1).AND.(eoshift(mk_loc,shift=1,dim=1).EQ.0)) .OR. &
    ((eoshift(lake_loc,shift=-1,dim=2).EQ.0).AND.(eoshift(ice_loc,shift=-1,dim=2).EQ.1).AND.(eoshift(mk_loc,shift=-1,dim=2).EQ.0)) .OR. &
    ((eoshift(lake_loc,shift=1,dim=2).EQ.0).AND.(eoshift(ice_loc,shift=1,dim=2).EQ.1).AND.(eoshift(mk_loc,shift=1,dim=2).EQ.0))))
    clake(:,:) = 1
  endwhere



! point terre et englace : mk=0 et H>0
  do j=1,ny
    do i=1,nx
      if (clake(i,j).EQ.1) then
        do k=1,80
          ipinctab(k) = max(1,min(nx,i + inclake_loc(k,1)))
          jpinctab(k) = max(1,min(ny,j + inclake_loc(k,2)))
          mask_extlake_loc(ipinctab(k),jpinctab(k)) = 1
          orog_extlake_loc(ipinctab(k),jpinctab(k))=orog_lake_loc(i,j)
        enddo
      endif
    enddo
  enddo

end subroutine extlake

end module lake_rsl