!> \file lect-Ant_CISM_15_dat.f90
!!Module pour la lecture de la topographie
!<

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

module lect_topo_ant_CISM_15

  use module3D_phy
  use param_phy_mod
  use interface_input
  character(len=100) :: topo_surf      !< surface
  character(len=100) :: correc_surf    !< ice-real correction
  character(len=100) :: topo_thick     !< thickness
  character(len=100) :: topo_bed       !< bedrock
  character(len=100) :: mask_grounded  !< mask
  character(len=100) :: longitude      !< longitude
  character(len=100) :: latitude       !< latitude
  character(len=100) :: heatflux       !< geothermal heat flux

  real,dimension(nx,ny) :: work_tab 
  Character(len=100) :: file_ncdf      !< fichier netcdf issue des fichiers .dat

contains

  subroutine input_topo

    namelist/topo_CISM_15/topo_surf,correc_surf,topo_thick,topo_bed,mask_grounded,longitude,latitude,heatflux

428 format(A)

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

    write(num_rep_42,428)'!___________________________________________________________' 
    write(num_rep_42,428)'!  module  lect_topo_ant_CISM        '
    write(num_rep_42,topo_CISM_15)
    write(num_rep_42,428)'!___________________________________________________________' 


    topo_surf     = trim(dirnameinp)//trim(topo_surf)
    correc_surf   = trim(dirnameinp)//trim(correc_surf)
    topo_thick    = trim(dirnameinp)//trim(topo_thick)
    topo_bed      = trim(dirnameinp)//trim(topo_bed)
    mask_grounded = trim(dirnameinp)//trim(mask_grounded)
    longitude     = trim(dirnameinp)//trim(longitude)
    latitude      = trim(dirnameinp)//trim(latitude)
    heatflux      = trim(dirnameinp)//trim(heatflux)

    file_ncdf= trim(dirnameinp)//trim(runname)//'.nc'


    ! lecture adaptee aux fichiers ZBL.dat

    !call lect_datfile(nx,ny,Bsoc,1,topo_bed)               ! bedrock
    call lect_input(3,'Bsoc',1,Bsoc,topo_bed,file_ncdf)

    Bsoc0(:,:) = Bsoc(:,:)

    call lect_input(3,'S',1,S,topo_surf,file_ncdf)
    !call lect_datfile(nx,ny,S,1,topo_surf)                 ! surface 


    call lect_input(3,'H',1,H,topo_thick,file_ncdf)
    !call lect_datfile(nx,ny,H,1,topo_thick)                ! thickness


    !call lect_datfile(nx,ny,work_tab,1,correc_surf)        ! density correction
    call lect_input(3,'work_tab1',1,work_tab,correc_surf,file_ncdf)

    where (H(:,:).gt.0.)
       work_tab(:,:)=min(H(:,:),work_tab(:,:))   
       S(:,:) = S(:,:) - work_tab(:,:)                     ! corrected surface
       H(:,:) = H(:,:) - work_tab(:,:)                     ! corrected thickness
    end where
    H(:,:)=max(0.,H(:,:))
    H0(:,:)= H(:,:)
    S0(:,:)= S(:,:)


    call lect_input(3,'work_tab2',1,work_tab,mask_grounded,file_ncdf)
    !call lect_datfile(nx,ny,work_tab,1,mask_grounded)      ! masque

    mk0(:,:) = int(work_tab)


    ! pour l'Antarctique masque mko vrai partout (version 2006)
    MK0(:,:)=1

    ! determination des flot
    do j=1,ny
       do i=1,nx
          if ((Bsoc(i,j)+H(i,j)*ro/row -sealevel).lt.0.) then
             flot(i,j)=.true.
          else
             flot(i,j)=.false.
          endif
       enddo
    enddo

    ! le calcul des courbures du socle devrait etre dans le dragging

    call lect_input(1,'Xlong',1,Xlong,longitude,file_ncdf)
    !call lect_datfile(nx,ny,Xlong,1,longitude)      ! lecture des coordonnées geographiques

    call lect_input(1,'Ylat',1,Ylat,latitude,file_ncdf)
    !call lect_datfile(nx,ny,Ylat,1,latitude)        !

    open(22,file=longitude)
    read(22,*) i,i,xmin,xmax,ymin,xmax
    xmin=xmin/1000.
    ymin=ymin/1000.
    xmax=xmax/1000.
    ymax=ymax/1000.
    close(22)


    !     lecture du fichier de reference pour le calcul du niveau des mers : etat actuel
    !     pas de version correcte pour l'instant -> valeurs initiales
    S_sealev(:,:) = S0(:,:)
    H_sealev(:,:) = H0(:,:)
    B_sealev(:,:) = Bsoc(:,:)
    M_sealev(:,:) = Mk0(:,:)


    call lect_input(1,'ghf',1,ghf,heatflux,file_ncdf)
    !call lect_datfile(nx,ny,ghf,1,heatflux)


    ! for the purposes of SeaRISE, GTHF should be capped at a value of -0.07 w/m^2,
    ! and the field provided in the netCDF file does not include these updated values

    ghf(:,:)=min(ghf(:,:),0.07)

    ! pour passer les flux des W/m2 au J/m2/an      
    ghf(:,:)=-SECYEAR*ghf(:,:)


  end subroutine input_topo


end module lect_topo_ant_CISM_15