source: trunk/SOURCES/Greeneem_files/lect-greeneem_mod.f90 @ 13

module lect_topo_greeneem

use module3D_phy
use param_phy_mod  

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=100) :: filin
real, dimension(nx,ny,5) :: bidon        ! pour l'appel a courbure
!real :: ghf0                             ! flux geothermique constant en mW/m2

integer :: n1,n2,ndx

contains
 
subroutine input_topo

namelist/toponeem/topo_ref,topo_dep,grid_topo,ghf_fich

428 format(A)

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

read(num_param,toponeem)

write(num_rep_42,428)'!___________________________________________________________' 
write(num_rep_42,428) '&toponeen                     ! module lect_topo_greeneem'
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,*)         'ghf_fich  =', ghf_fich
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,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)!trim(ghf0)

sealevel0=0.  ! voir a passer dans le fichier parametre

!     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,*) 
     read(20,*) 
     read(20,*) n1,n2,ndx
     read(20,*) 

     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
         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),sealevel0)   ! pour etre au niveau des mers : ATTENTION si SEALEV <0
          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,*) 
     read(20,*) 
     read(20,*) n1,n2,ndx
     read(20,*) 


     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
  
         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,*) 
read(20,*) 
read(20,*) n1,n2,ndx
read(20,*) 

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
   
   do k=1,nx*ny
      read(20,*) i,j,xcc(i,j),ycc(i,j),xlong(i,j),ylat(i,j)
   enddo
   ! pour assurer la continuite des temperatures parametrees :
   where(xlong(:,:).lt.100.)    
      xlong(:,:)=xlong(:,:)+360.
   end where
   
else
   write(6,*) 'le fichier ne correspond pas a la grille definie'
   STOP
end if
close(20)

! xcc est en m
xmin=xcc(1,1)
ymin=ycc(1,1)
do j=1,ny
   do i=1,nx
      xcc(i,j)=xcc(i,j)*1000.
      ycc(i,j)=ycc(i,j)*1000.
   end do
end do
xmax=xcc(nx,ny)/1000
ymax=ycc(nx,ny)/1000

! lecture de la carte de flux geothermique
open (20,file=ghf_fich)
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
   do j=1,ny
      do i=1,nx
         read(20,*) ghf(i,j)
      end do
   end do
   close(20)
else
   write(6,*) 'le fichier ne correspond pas a la grille definie'
   STOP
end if


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



! 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
    
!------------------------------------------------     
end subroutine input_topo


end module lect_topo_greeneem