!> \file lect-greenmint_mod.f90
!!Module pour la lecture de la topography
!<

!> \namespace lect_topo_greenmint
!! Module pour la lecture de la topography
!! \author ...
!! \date ...
!! @note Used module
!! @note   - use module3D_phy
!<

module lect_topo_greenmint

use module3D_phy
  
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) :: filin
real,dimension(nx,ny) ::  xcc , ycc      !< coordeonnes en m
real :: ghf0                             !< flux geothermique constant en mW/m2

integer :: n1,n2,ndx

contains
 
subroutine input_topo

namelist/topogreen/topo_ref,topo_dep,grid_topo,ghf0

428 format(A)

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

write(num_rep_42,428)'!___________________________________________________________' 
write(num_rep_42,428) '&topogreen                     ! module lect_topo_greenmint'
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,*)         'ghf0      =', ghf0
write(num_rep_42,*)'/'                      
write(num_rep_42,428) '! topo_ref= topo actuelle'
write(num_rep_42,428) '! ghf0 flux geothermique en mW/m2'
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)

!     lecture de la topo de référence
!     -------------------------------
! Cette topo sert a calculer le socle de reference pour l'isostasie
! voir init_iso et a avoir une surface de reference pour les temperatures

     open (20,file=topo_ref)
     read(20,*) n1,n2,ndx


     if ((n1.eq.nx).and.(n2.eq.ny).and.(ndx.eq.nint(dx/1000.))) then
                               ! le fichier correspond bien a 
                               ! la grille definie dans paradim
         read(20,*)            ! saute la ligne nom des colonnes
         do J=1,ny 
          do I=1,nx
            read (20,*) S0(I,J),H0(I,J),Bsoc0(I,J),precip(i,j)

             S0(i,j)=max(S0(i,j),0.)   ! pour être au niveau des mers

          end do
        end do
        
     else
        write(6,*) 'le fichier ne correspond pas a la grille definie'
        STOP
     end if
        
     close(20)

mk0(:,:) = 1                   ! mk0=0 pour les zones interdites



     
!     lecture de la topo de depart
!     ---------------------------
     open (20,file=topo_dep)
     read(20,*) n1,n2,ndx

     if ((n1.eq.nx).and.(n2.eq.ny).and.(ndx.eq.nint(dx/1000.))) then
                               ! le fichier correspond bien a 
                               ! la grille definie dans paradim
         read(20,*)            ! saute la ligne nom des colonnes     
         do J=1,ny 
          do I=1,nx
              read (20,*) S(I,J),H(I,J),BSOC(I,J)
              S(i,j)=max(S(i,j),0.)   ! pour etre au niveau des mers : ATTENTION si SEALEV <0
              H(i,j)=max(H(i,j),1.)   ! pour avoir au moins 1 m
          end do
        end do
  
     else
        write(6,*) 'le fichier ne correspond pas a la grille definie'
        STOP
     end if
     close(20)

! calcul de flottaison
do j=1,ny
   do i=1,nx

      archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel

      if ((ARCHIM.LT.0.).and.(H(I,J).gt.1.E-3)) then    ! le point flotte
         flot(i,j)=.true.
      else
         flot(i,j)=.false.
      end if
    
   end do
end do

! lecture des coordonnées geographiques



! les fichiers EISMINT etaient ranges i,j,lat,lon
! je n'arrive pas a retrouver la projection exacte.
! le x est donc calcule en partant de 0, 0

open(unit=20,file=grid_topo,iostat=ios)
read(20,*) n1,n2,ndx
  
     if ((n1.eq.nx).and.(n2.eq.ny).and.(ndx.eq.nint(dx/1000.))) then
                               ! le fichier correspond bien a 
                               ! la grille definie dans paradim
         read(20,*)            ! saute la ligne nom des colonnes    

         do k=1,nx*ny
            read(20,*) i,j,ylat(i,j),xlong(i,j)
         enddo
      else
        write(6,*) 'le fichier ne correspond pas a la grille definie'
        STOP
     end if
     close(20)
xmin=0.
ymin=0.

do j=1,ny
   do i=1,nx
      xcc(i,j)=xmin+dx*(i-1)
      ycc(i,j)=ymin+dx*(j-1)
   end do
end do

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


!  flux geothermique  : pour passer les flux des mW/m2 aux J/m2/an     
ghf(:,:)=-secyear/1000.*ghf0




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


end module lect_topo_greenmint