source: trunk/SOURCES/Grismip6_files/lect-grismip6_mod.f90 @ 334

module lect_topo_grismip6

  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) :: ghf_fich       ! geothermal heat flux
  character(len=100) :: lon_fich       ! longitude
  character(len=100) :: lat_fich       ! latitude
  character(len=80) :: filin
  character(len=100) :: file_ncdf      !< fichier netcdf issue des fichiers .dat

  real,dimension(nx,ny) :: varloc

contains
 
subroutine input_topo

  integer :: ios

  namelist/topo_file_grismip6/topo_ref,topo_dep,ghf_fich,lon_fich,lat_fich
  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
  read(num_param,topo_file_grismip6)

  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)') 'ghf_fich = "',trim(ghf_fich),'"'
  write(num_rep_42,'(A,A,A)') 'lon_fich = "',trim(lon_fich),'"'
  write(num_rep_42,'(A,A,A)') 'lat_fich = "',trim(lat_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)') '! ghf_fich  : fichier flux geothermique'
  write(num_rep_42,'(A)') '! lon_fich  : fichier longitude'
  write(num_rep_42,'(A)') '! lat_fich  : fichier latitude'
  
  write(num_rep_42,*)

  topo_ref=trim(dirnameinp)//trim(topo_ref) 
  topo_dep=trim(dirnameinp)//trim(topo_dep) 
  ghf_fich=trim(dirnameinp)//trim(ghf_fich)
  lon_fich=trim(dirnameinp)//trim(lon_fich)
  lat_fich=trim(dirnameinp)//trim(lat_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(:,:)
  H0(:,:)=max(0.,H0(:,:))
  where(H0(:,:).le.20) H0(:,:)=0


!     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(:,:)
  H(:,:)=max(0.,H(:,:))
  where(H(:,:).le.20) H(:,:)=0

! determination des flot
  do j=1,ny
     do i=1,nx
        if ((bsoc(i,j)+h(i,j)*ro/row -sealevel_2d(i,j)).LT.0.) then
           flot(i,j)=.TRUE.
        else
           flot(i,j)=.FALSE.
        endif
     enddo
  enddo

  xmin=-720000./1000.
  ymin=-3450000./1000.
  xmax=-960000./1000.
  ymax=-570000./1000.
                       
  call lect_input(1,'ghf',1,ghf,ghf_fich,file_ncdf)

! pour passer les flux des W/m2 au J/m2/an      
  ghf(:,:)=-SECYEAR*ghf(:,:)
  
!------------------------------------------------
! pour l'Hemisphere Nord mko vrai partout (version 2006)
  MK0(:,:)=1

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

end module lect_topo_grismip6