module lect_topo_hemin40

  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

  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)
! formats pour les ecritures dans 42
428 format(A)
  write(num_rep_42,428)'!___________________________________________________________' 
  write(num_rep_42,428) '&topo_file                                  ! nom du bloc '
  write(num_rep_42,*)
  write(num_rep_42,'(A,A)') 'topo_ref = ', topo_ref
  write(num_rep_42,'(A,A)') 'topo_dep = ', topo_dep
  write(num_rep_42,'(A,A)') 'grid_topo =', grid_topo
  write(num_rep_42,'(A,A)') 'ghf_fich = ', ghf_fich
  write(num_rep_42,*)'/'                      
  write(num_rep_42,428) '! topo_ref= topo ref isostasie'
  write(num_rep_42,428) '! topo_dep= topo de depart'
  write(num_rep_42,428) '! grid_topo : fichier i,j,x,y,lon,lat'
  write(num_rep_42,428) '! 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)

!write(num_rep_42,*) 'file1 = ', file1
!write(num_rep_42,*) 'file2 = ', file2
!write(num_rep_42,*)'/' 
!write(num_rep_42,428) '! file1 : topo de depart'
!write(num_rep_42,428) '! file2 : topo de reference'             
!write(num_rep_42,*)

!====================================== La reponse est 42 ===========
! write(42,*)
! write(42,*)' Fichiers en entree'
! write(42,*)'----------------------'
!==================================================================== 

! dans param :
!      file1=TRIM(DIRNAMEINP)//'topo-21k.g40'     ! topo LGM ICE_5G (1=topo de depart)
!      file1=TRIM(DIRNAMEINP)//'hemin2.g40'
!      file2=TRIM(DIRNAMEINP)//'hemin2.g40'       ! topo actuelle
!      write(42,*) 'topo de depart', file1
!      write(42,*) 'topo reference', file2

      
!!$! lecture adaptee aux fichiers intercomparaison EISMINT
!!$       nxx=nx
!!$       nyy=ny

!!$!     lecture de la topo actuelle
!!$!     ---------------------------
!!$     open (20,file=TRIM(DIRNAMEINP)//file2,status='old')
!!$       
!!$     read(20,'(A80)') TITRE
!!$     read(20,*) NI,NJ,NXX,NYY,STEP
!!$     read(20,*)
!!$         do J=1,ny 
!!$          do I=1,nx
!!$             read (20,*)  S0(I,J),H0(I,J),BSOC0(I,J)
!!$             S0(i,j)=max(S0(i,j),0.)
!!$          end do
!!$        end do
!!$     close(20)
!!$
!!$     
!!$!     lecture de la topo de depart
!!$!     ---------------------------
!!$     open (20,file=TRIM(DIRNAMEINP)//file1,status='old')
!!$!         open (20,file='../INPUT-DATA/hemin.g50')
!!$     read(20,'(A80)') TITRE
!!$     read(20,*) NI,NJ,NXX,NYY,STEP
!!$     read(20,*)
!!$         do J=1,ny 
!!$          do I=1,nx
!!$             read (20,*) S(I,J),H(I,J),BSOC(I,J)
!!$          end do
!!$        end do
!!$     close(20)


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

!     lecture de la topo de départ
    call lect_input(1,'Bsoc',1,Bsoc,topo_dep,file_ncdf)    ! socle
    call lect_input(1,'S',1,S,topo_dep,file_ncdf)          ! surface
    call lect_input(1,'H',1,H,topo_dep,file_ncdf)          ! epaisseur


! calcul des courbures du socle
     call courbure(nx,ny,dx,Bsoc,bidon(:,:,1),bidon(:,:,2),bidon(:,:,3), &
          bidon(:,:,4),socle_cry,bidon(:,:,5))
     socle_cry(:,:)=socle_cry(:,:)*dx*dx

! lecture des coordonnées geographiques

!    filin=TRIM(DIRNAMEINP)//grid_topo

! 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)
    write(42,*) 'fichier grille: ', filin
              
    xmin=xcc(1,1)/1000.
    ymin=ycc(1,1)/1000.
    xmax=xcc(nx,ny)/1000.
    ymax=ycc(nx,ny)/1000.
                         
! lecture du flux geothermique de Shapiro
!    open(88,file=TRIM(DIRNAMEINP)//'ijphi_hemin40.dat')
!    open(88,file=ghf_fich)

!!$    write(42,*) 'flux geothermique Shapiro : ',TRIM(DIRNAMEINP)//'ijphi_hemin40.dat'
!!$    do k=1,nx*ny
!!$       read(88,*) i,j,ghf(i,j)
!!$    end do
!!$    close(88)

    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.*ghf(:,:)
!     write(42,*) 'flux geothermique fixe : 55 mW/m2' 
!     ghf(:,:)=-SECYEAR/1000.*55. !B6norcg2

! print*,'lect topo'
! print*,'shb',S(101,91),H(101,91),B(101,91)
! print*,'shb0',S0(101,91),H0(101,91),BSOC0(101,91)
!    Initialisation du Masque
!------------------------------------------------
! pour l'Hemisphere Nord mko vrai partout (version 2006)
    MK0(:,:)=1

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

end module lect_topo_hemin40