module lect_topo_alps

  use module3D_phy
  use interface_input
  use io_netcdf_grisli

  implicit none
  
  character(len=100) :: topo_dep       ! Topo de départ
  character(len=100) :: topo_ref       ! Topo de référence
  character(len=100) :: grid_topo      ! fichier grille
  character(len=100) :: ghf_fich       ! fichier grille
  character(len=80) :: filin
  real, dimension(nx,ny,5) :: bidon          ! pour l'appel a courbure
  character(len=100) :: file_ncdf      !< fichier netcdf issue des fichiers .dat

contains
 
subroutine input_topo

  integer :: ios
  real, dimension(nx,ny) :: varloc
  
  namelist/topo_file/topo_ref,topo_dep,grid_topo !,ghf_fich
  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
  read(num_param,topo_file)

  write(num_rep_42,'(A)')'!___________________________________________________________' 
  write(num_rep_42,'(A)') '&topo_file                                  ! nom du bloc '
  write(num_rep_42,*)
  write(num_rep_42,'(A,A,A)') 'topo_ref = "',trim(topo_ref),'"'
  write(num_rep_42,'(A,A,A)') 'topo_dep = "',trim(topo_dep),'"'
  write(num_rep_42,'(A,A,A)') 'grid_topo = "',trim(grid_topo),'"'
  !write(num_rep_42,'(A,A,A)') 'ghf_fich = "',trim(ghf_fich),'"'
  write(num_rep_42,*)'/'                      
  write(num_rep_42,'(A)') '! topo_ref= topo ref isostasie'
  write(num_rep_42,'(A)') '! topo_dep= topo de depart'
  write(num_rep_42,'(A)') '! grid_topo : fichier i,j,x,y,lon,lat'
  write(num_rep_42,'(A)') '! ghf_fich  : fichier flux geothermique'  
  write(num_rep_42,*)

  topo_ref=trim(dirnameinp)//trim(topo_ref) 
  topo_dep=trim(dirnameinp)//trim(topo_dep) 
  grid_topo=trim(dirnameinp)//trim(grid_topo)
  !ghf_fich=trim(dirnameinp)//trim(ghf_fich)



!     lecture de la topo de référence
  call lect_input(1,'Bsoc',1,varloc,topo_ref,file_ncdf)    ! socle
  Bsoc0(:,:)=varloc(:,:)
  call lect_input(1,'S',1,varloc,topo_ref,file_ncdf)          ! surface
  S0(:,:)=varloc(:,:)
  call lect_input(1,'H',1,varloc,topo_ref,file_ncdf)          ! epaisseur
  H0(:,:)=varloc(:,:)
  
!     lecture de la topo de départ
  call lect_input(1,'Bsoc',1,varloc,topo_dep,file_ncdf)    ! socle
  Bsoc(:,:)=varloc(:,:)
  call lect_input(1,'S',1,varloc,topo_dep,file_ncdf)          ! surface
  S(:,:)=varloc(:,:)
  call lect_input(1,'H',1,varloc,topo_dep,file_ncdf)          ! epaisseur
  H(:,:)=varloc(:,:)


! 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=grid_topo,iostat=ios)
  do k=1,nx*ny
     read(2004,*) i,j,XCC(i,j),YCC(i,j),XLONG(i,j),YLAT(i,j)
  enddo
  close(2004)

  where(xlong(:,:).lt.0.) xlong(:,:)=xlong(:,:)+360

  xmin=xcc(1,1)/1000.
  ymin=ycc(1,1)/1000.
  xmax=xcc(nx,ny)/1000.
  ymax=ycc(nx,ny)/1000.

  !!call lect_input(1,'ghf',1,ghf,ghf_fich,file_ncdf)

! pour passer les flux des mW/m2 au J/m2/an      
  ghf(:,:)=-SECYEAR/1000.* 55 !ghf(:,:)
!     write(42,*) 'flux geothermique fixe : 55 mW/m2' 
!     ghf(:,:)=-SECYEAR/1000.*55. !B6norcg2




!    Initialisation du Masque
!------------------------------------------------
! pour l'Hemisphere Nord mko vrai partout (version 2006)
  MK0(:,:)=1

!------------------------------------------------     
end subroutine input_topo

end module lect_topo_alps