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

module lect_topo_greeneem

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

real, dimension(nx,ny) :: Bsoc_bamber, S_bamber, Bsoc_new
character(len=100) :: topo_surf      !< surface file name
character(len=100) :: topo_bed       !< bedrock file name
character(len=100) :: file_ncdf      !< fichier netcdf issue des fichiers .dat

integer :: n1,n2,ndx
real :: sealevel0

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) '! 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)!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
    ! lecture adaptee aux fichiers ZBL.dat ou netcdf ou grd
    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

    S0(:,:)=max(S0(:,:),sealevel0)   ! pour etre au niveau des mers : ATTENTION si SEALEV <0
    mk0(:,:) = 1                   ! mk0=0 pour les zones interdites

    !     lecture de la topo de depart
    !     ---------------------------
    ! lecture adaptee aux fichiers ZBL.dat ou netcdf ou grd
    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
    S(:,:)=max(S(:,:),0.)   ! pour etre au niveau des mers : ATTENTION si SEALEV <0
    H(:,:)=max(H(:,:),1.)   ! pour avoir au moins 1 m

          
    ! 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