module tabfi_mod
  !=========================================
  ! Author: J. Mauxion (JM in the following)
  ! Adapted from tabfi
  !=========================================

  implicit none

contains

!==============================================================================
subroutine tabfi(nid, lmodif, tab0, day_ini, lmax,                    &
                 p_rad, p_omeg, p_g, p_mugaz, p_daysec, time)
!==============================================================================

!----------------------------------------------------------------------------
! Read physical control parameters and initialize physical constants.
!----------------------------------------------------------------------------

   use comsoil_h, only : volcapa
   use surfdat_h, only : z0_default, emissiv, emisice, albedice,      &
                         iceradius, dtemisice
   use dimradmars_mod, only : tauvis
   use iostart, only : get_var
   use mod_phys_lmdz_para, only : is_parallel
   use comcstfi_h, only : g, mugaz, omeg, rad, rcp
   use time_phylmdz_mod, only : daysec, dtphys
   use planete_h, only : aphelie, emin_turb, lmixmin,                 &
                         obliquit, peri_day, periheli, year_day,       &
                         lsp2solp, iniorbit

   implicit none

!----------------------------------------------------------------------------
! Arguments
!----------------------------------------------------------------------------

   integer, intent(in) :: nid
   integer, intent(in) :: lmodif
   integer, intent(in) :: tab0

   integer, intent(out) :: day_ini
   integer, intent(out) :: lmax

   real, intent(out) :: p_rad
   real, intent(out) :: p_omeg
   real, intent(out) :: p_g
   real, intent(out) :: p_mugaz
   real, intent(out) :: p_daysec
   real, intent(out) :: time

!----------------------------------------------------------------------------
! Local variables
!----------------------------------------------------------------------------

   integer, parameter :: length = 100

   integer :: ierr

   real :: peri_ls
   real :: tab_cntrl(length)

   logical :: found

   character(len=20) :: modif
   character(len=5)  :: modname = "tabfi"

!----------------------------------------------------------------------------
! Initial message
!----------------------------------------------------------------------------

   write (*,*) "tabfi: nid=", nid, " tab0=", tab0, " Lmodif=", lmodif

!==============================================================================
! Default initialization (nid == 0)
!==============================================================================

   if (nid == 0) then

      ! Initialize control array

      tab_cntrl = 0.0

      lmax    = 0
      day_ini = 0
      time    = 0.0

!----------------------------------------------------------------------------
! Planetary constants
!----------------------------------------------------------------------------

      rad    = 3397200.0
      daysec = 88775.0
      omeg   = 4.0 * asin(1.0) / daysec
      g      = 3.72
      mugaz  = 43.49
      rcp    = 0.256793

!----------------------------------------------------------------------------
! Orbital parameters
!----------------------------------------------------------------------------

      year_day = 669.0
      periheli = 206.66
      aphelie  = 249.22
      peri_day = 485.0
      obliquit = 25.19

!----------------------------------------------------------------------------
! Boundary layer
!----------------------------------------------------------------------------

      z0_default = 1.0e-2
      emin_turb  = 1.0e-6
      lmixmin    = 30

!----------------------------------------------------------------------------
! Surface properties
!----------------------------------------------------------------------------

      emissiv = 0.95

      emisice(1) = 0.95
      emisice(2) = 0.95

      albedice(1) = 0.65
      albedice(2) = 0.65

      iceradius(1) = 100.0e-6
      iceradius(2) = 100.0e-6

      dtemisice(1) = 0.4
      dtemisice(2) = 0.4

!----------------------------------------------------------------------------
! Dust
!----------------------------------------------------------------------------

      tauvis = 0.2

!----------------------------------------------------------------------------
! Soil
!----------------------------------------------------------------------------

      volcapa = 1.0e6

!==============================================================================
! Read control array from restart file
!==============================================================================

   else

!==============================================================================
! Read control array from restart file
!==============================================================================

      call get_var("controle", tab_cntrl, found)

      if (.not. found) then

         call abort_physic(modname, &
              "tabfi: Failed reading <controle> array", 1)

      else

         write (*,*) "tabfi: tab_cntrl", tab_cntrl

      end if

!==============================================================================
! Grid / time information
!==============================================================================

      lmax    = nint(tab_cntrl(tab0 + 2))
      day_ini = tab_cntrl(tab0 + 3)
      time    = tab_cntrl(tab0 + 4)

      write (*,*) "IN tabfi day_ini =", day_ini

!==============================================================================
! Planetary constants (dynamics + physics)
!==============================================================================

      rad     = tab_cntrl(tab0 + 5)
      omeg    = tab_cntrl(tab0 + 6)
      g       = tab_cntrl(tab0 + 7)
      mugaz   = tab_cntrl(tab0 + 8)
      rcp     = tab_cntrl(tab0 + 9)
      daysec  = tab_cntrl(tab0 + 10)
      dtphys  = tab_cntrl(tab0 + 11)

!==============================================================================
! Orbital parameters
!==============================================================================

      year_day = tab_cntrl(tab0 + 14)
      periheli = tab_cntrl(tab0 + 15)
      aphelie  = tab_cntrl(tab0 + 16)
      peri_day = tab_cntrl(tab0 + 17)
      obliquit = tab_cntrl(tab0 + 18)

!==============================================================================
! Boundary layer / turbulence
!==============================================================================

      z0_default = tab_cntrl(tab0 + 19)
      lmixmin    = tab_cntrl(tab0 + 20)
      emin_turb  = tab_cntrl(tab0 + 21)

!==============================================================================
! Surface properties
!==============================================================================

      albedice(1) = tab_cntrl(tab0 + 22)
      albedice(2) = tab_cntrl(tab0 + 23)

      emisice(1)  = tab_cntrl(tab0 + 24)
      emisice(2)  = tab_cntrl(tab0 + 25)

      emissiv     = tab_cntrl(tab0 + 26)

!==============================================================================
! Dust
!==============================================================================

      tauvis = tab_cntrl(tab0 + 27)

!==============================================================================
! CO2 ice properties
!==============================================================================

      iceradius(1) = tab_cntrl(tab0 + 31)
      iceradius(2) = tab_cntrl(tab0 + 32)

      dtemisice(1) = tab_cntrl(tab0 + 33)
      dtemisice(2) = tab_cntrl(tab0 + 34)

!==============================================================================
! Soil properties
!==============================================================================

      volcapa = tab_cntrl(tab0 + 35)

!==============================================================================
! Save constants for output
!==============================================================================

      p_omeg   = omeg
      p_g      = g
      p_mugaz  = mugaz
      p_daysec = daysec
      p_rad    = rad

   end if

!==============================================================================
! Print control values BEFORE modifications
!==============================================================================

   write(*,*) '*****************************************************'
   write(*,*) 'Reading tab_cntrl BEFORE changes'
   write(*,*) '*****************************************************'

   write(*,'(a20, f12.2, f12.2)') '(1) ngrid?',    tab_cntrl(tab0+1), real(tab_cntrl(tab0+1))
   write(*,'(a20, f12.2, f12.2)') '(2) lmax',      tab_cntrl(tab0+2), real(lmax)
   write(*,'(a20, f12.2, f12.2)') '(3) day_ini',   tab_cntrl(tab0+3), real(day_ini)
   write(*,'(a20, f12.2, f12.2)') '(5) rad',       tab_cntrl(tab0+5), rad
   write(*,'(a20, f12.2, f12.2)') '(10) daysec',   tab_cntrl(tab0+10), daysec

   write(*,'(a20, e15.6, e15.6)') '(6) omeg',      tab_cntrl(tab0+6), omeg
   write(*,'(a20, f12.2, f12.2)') '(7) g',         tab_cntrl(tab0+7), g
   write(*,'(a20, f12.2, f12.2)') '(8) mugaz',     tab_cntrl(tab0+8), mugaz
   write(*,'(a20, f12.2, f12.2)') '(9) rcp',       tab_cntrl(tab0+9), rcp
   write(*,'(a20, e15.6, e15.6)') '(11) dtphys',   tab_cntrl(tab0+11), dtphys

   write(*,'(a20, f12.2, f12.2)') '(14) year_day', tab_cntrl(tab0+14), year_day
   write(*,'(a20, f12.2, f12.2)') '(15) periheli', tab_cntrl(tab0+15), periheli
   write(*,'(a20, f12.2, f12.2)') '(16) aphelie',  tab_cntrl(tab0+16), aphelie
   write(*,'(a20, f12.2, f12.2)') '(17) peri_day', tab_cntrl(tab0+17), peri_day
   write(*,'(a20, f12.2, f12.2)') '(18) obliquit', tab_cntrl(tab0+18), obliquit

   write(*,'(a20, e15.6, e15.6)') '(19) z0_default', tab_cntrl(tab0+19), z0_default
   write(*,'(a20, e15.6, e15.6)') '(21) emin_turb',  tab_cntrl(tab0+21), emin_turb
   write(*,'(a20, f12.2, f12.2)') '(20) lmixmin',    tab_cntrl(tab0+20), real(lmixmin)

   write(*,'(a20, f12.2, f12.2)') '(26) emissiv',    tab_cntrl(tab0+26), emissiv
   write(*,'(a20, f12.2, f12.2)') '(24) emisice(1)', tab_cntrl(tab0+24), emisice(1)
   write(*,'(a20, f12.2, f12.2)') '(25) emisice(2)', tab_cntrl(tab0+25), emisice(2)
   write(*,'(a20, f12.2, f12.2)') '(22) albedice(1)',tab_cntrl(tab0+22), albedice(1)
   write(*,'(a20, f12.2, f12.2)') '(23) albedice(2)',tab_cntrl(tab0+23), albedice(2)

   write(*,'(a20, e15.6, e15.6)') '(31) iceradius(1)', tab_cntrl(tab0+31), iceradius(1)
   write(*,'(a20, e15.6, e15.6)') '(32) iceradius(2)', tab_cntrl(tab0+32), iceradius(2)

   write(*,'(a20, f12.2, f12.2)') '(33) dtemisice(1)', tab_cntrl(tab0+33), dtemisice(1)
   write(*,'(a20, f12.2, f12.2)') '(34) dtemisice(2)', tab_cntrl(tab0+34), dtemisice(2)

   write(*,'(a20, f12.2, f12.2)') '(27) tauvis',      tab_cntrl(tab0+27), tauvis
   write(*,'(a20, f12.2, f12.2)') '(35) volcapa',      tab_cntrl(tab0+35), volcapa

   write(*,*)
   write(*,*) 'Lmodif =', lmodif

!==============================================================================
! Interactive modification section
!==============================================================================

   if (lmodif == 1) then

      write(*,*)
      write(*,*) 'Change values in tab_cntrl ?'
      write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      write(*,*) '(Current values shown above)'

      write(*,*) '(3) day_ini : initial day'
      write(*,*) '(19) z0_default : surface roughness'
      write(*,*) '(21) emin_turb : minimal PBL energy'
      write(*,*) '(20) lmixmin : mixing length'
      write(*,*) '(26) emissiv : ground emissivity'
      write(*,*) '(24-25) emisice : CO2 ice emissivity'
      write(*,*) '(22-23) albedice : CO2 ice albedo'
      write(*,*) '(31-32) iceradius : CO2 snow radius'
      write(*,*) '(33-34) dtemisice : snow metamorphism time scale'
      write(*,*) '(27) tauvis : dust optical depth'
      write(*,*) '(35) volcapa : soil heat capacity'
      write(*,*) '(18) obliquit : obliquity'
      write(*,*) '(17) peri_day : perihelion date'
      write(*,*) '(15) periheli : perihelion distance'
      write(*,*) '(16) aphelie : aphelion distance'

!==============================================================================
! Main modification loop (replaces all GOTOs 101–119)
!==============================================================================

      do

         write(*,*)
         write(*,*) 'Enter parameter to modify (blank to exit):'
         read(*,'(a20)') modif

         if (modif(1:1) == ' ') exit

         write(*,*) trim(modif), ':'

!----------------------------------------------------------------------------
! day_ini
!----------------------------------------------------------------------------

         if (trim(modif) == 'day_ini') then

            write(*,*) 'current value:', day_ini
            write(*,*) 'enter new value:'

            do
               read(*,*, iostat=ierr) day_ini
               if (ierr == 0) exit
            end do

!----------------------------------------------------------------------------
! z0_default
!----------------------------------------------------------------------------

         else if (trim(modif) == 'z0') then

            write(*,*) 'current value:', z0_default
            write(*,*) 'enter new value:'

            do
               read(*,*, iostat=ierr) z0_default
               if (ierr == 0) exit
            end do

!----------------------------------------------------------------------------
! emin_turb
!----------------------------------------------------------------------------

         else if (trim(modif) == 'emin_turb') then

            write(*,*) 'current value:', emin_turb
            write(*,*) 'enter new value:'

            do
               read(*,*, iostat=ierr) emin_turb
               if (ierr == 0) exit
            end do

!----------------------------------------------------------------------------
! lmixmin
!----------------------------------------------------------------------------

         else if (trim(modif) == 'lmixmin') then

            write(*,*) 'current value:', lmixmin
            write(*,*) 'enter new value:'

            do
               read(*,*, iostat=ierr) lmixmin
               if (ierr == 0) exit
            end do

!----------------------------------------------------------------------------
! emissiv
!----------------------------------------------------------------------------

         else if (trim(modif) == 'emissiv') then

            write(*,*) 'current value:', emissiv
            write(*,*) 'enter new value:'

            do
               read(*,*, iostat=ierr) emissiv
               if (ierr == 0) exit
            end do

!----------------------------------------------------------------------------
! emisice (both hemispheres)
!----------------------------------------------------------------------------

         else if (trim(modif) == 'emisice') then

            write(*,*) 'current value emisice(1) North:', emisice(1)
            write(*,*) 'enter new value:'

            do
               read(*,*, iostat=ierr) emisice(1)
               if (ierr == 0) exit
            end do

            write(*,*) 'current value emisice(2) South:', emisice(2)
            write(*,*) 'enter new value:'

            do
               read(*,*, iostat=ierr) emisice(2)
               if (ierr == 0) exit
            end do

!----------------------------------------------------------------------------
! albedice
!----------------------------------------------------------------------------

         else if (trim(modif) == 'albedice') then

            write(*,*) 'current value albedice(1):', albedice(1)
            write(*,*) 'enter new value:'

            do
               read(*,*, iostat=ierr) albedice(1)
               if (ierr == 0) exit
            end do

            write(*,*) 'current value albedice(2):', albedice(2)
            write(*,*) 'enter new value:'

            do
               read(*,*, iostat=ierr) albedice(2)
               if (ierr == 0) exit
            end do

!----------------------------------------------------------------------------
! iceradius
!----------------------------------------------------------------------------

         else if (trim(modif) == 'iceradius') then

            write(*,*) 'current value iceradius(1):', iceradius(1)
            write(*,*) 'enter new value:'

            do
               read(*,*, iostat=ierr) iceradius(1)
               if (ierr == 0) exit
            end do

            write(*,*) 'current value iceradius(2):', iceradius(2)
            write(*,*) 'enter new value:'

            do
               read(*,*, iostat=ierr) iceradius(2)
               if (ierr == 0) exit
            end do

!----------------------------------------------------------------------------
! dtemisice
!----------------------------------------------------------------------------

         else if (trim(modif) == 'dtemisice') then

            write(*,*) 'current value dtemisice(1):', dtemisice(1)
            write(*,*) 'enter new value:'

            do
               read(*,*, iostat=ierr) dtemisice(1)
               if (ierr == 0) exit
            end do

            write(*,*) 'current value dtemisice(2):', dtemisice(2)
            write(*,*) 'enter new value:'

            do
               read(*,*, iostat=ierr) dtemisice(2)
               if (ierr == 0) exit
            end do

!----------------------------------------------------------------------------
! tauvis
!----------------------------------------------------------------------------

         else if (trim(modif) == 'tauvis') then

            write(*,*) 'current value:', tauvis
            write(*,*) 'enter new value:'

            do
               read(*,*, iostat=ierr) tauvis
               if (ierr == 0) exit
            end do

!----------------------------------------------------------------------------
! obliquit
!----------------------------------------------------------------------------

         else if (trim(modif) == 'obliquit') then

            write(*,*) 'current value:', obliquit
            write(*,*) 'enter new value:'

            do
               read(*,*, iostat=ierr) obliquit
               if (ierr == 0) exit
            end do

!----------------------------------------------------------------------------
! peri_day
!----------------------------------------------------------------------------

         else if (trim(modif) == 'peri_day') then

            write(*,*) 'current value:', peri_day
            write(*,*) 'enter new value:'

            do
               read(*,*, iostat=ierr) peri_day
               if (ierr == 0) exit
            end do

!----------------------------------------------------------------------------
! peri_ls (conversion)
!----------------------------------------------------------------------------

         else if (trim(modif) == 'peri_ls') then

            write(*,*) 'enter peri_ls value:'

            do
               read(*,*, iostat=ierr) peri_ls
               if (ierr == 0) exit
            end do

            write(*,*) 'peri_ls asked:', peri_ls
            write(*,*) 'conversion using orbit parameters'

            call lsp2solp(peri_ls, peri_day)

            write(*,*) 'peri_day (new value):', peri_day

!----------------------------------------------------------------------------
! periheli
!----------------------------------------------------------------------------

         else if (trim(modif) == 'periheli') then

            write(*,*) 'current value:', periheli
            write(*,*) 'enter new value:'

            do
               read(*,*, iostat=ierr) periheli
               if (ierr == 0) exit
            end do

!----------------------------------------------------------------------------
! aphelie
!----------------------------------------------------------------------------

         else if (trim(modif) == 'aphelie') then

            write(*,*) 'current value:', aphelie
            write(*,*) 'enter new value:'

            do
               read(*,*, iostat=ierr) aphelie
               if (ierr == 0) exit
            end do

!----------------------------------------------------------------------------
! volcapa
!----------------------------------------------------------------------------

         else if (trim(modif) == 'volcapa') then

            write(*,*) 'current value:', volcapa
            write(*,*) 'enter new value:'

            do
               read(*,*, iostat=ierr) volcapa
               if (ierr == 0) exit
            end do

         end if

      end do

   end if   ! end of lmodif == 1

!==============================================================================
! Print control values AFTER modifications
!==============================================================================

   write(*,*)
   write(*,*) '*****************************************************'
   write(*,*) 'Reading tab_cntrl AFTER changes'
   write(*,*) '*****************************************************'

   write(*,'(a20, f12.2, f12.2)') '(1) ngrid?',    tab_cntrl(tab0+1), real(tab_cntrl(tab0+1))
   write(*,'(a20, f12.2, f12.2)') '(2) lmax',      tab_cntrl(tab0+2), real(lmax)
   write(*,'(a20, f12.2, f12.2)') '(3) day_ini',   tab_cntrl(tab0+3), real(day_ini)
   write(*,'(a20, f12.2, f12.2)') '(5) rad',       tab_cntrl(tab0+5), rad
   write(*,'(a20, f12.2, f12.2)') '(10) daysec',   tab_cntrl(tab0+10), daysec

   write(*,'(a20, e15.6, e15.6)') '(6) omeg',      tab_cntrl(tab0+6), omeg
   write(*,'(a20, f12.2, f12.2)') '(7) g',         tab_cntrl(tab0+7), g
   write(*,'(a20, f12.2, f12.2)') '(8) mugaz',     tab_cntrl(tab0+8), mugaz
   write(*,'(a20, f12.2, f12.2)') '(9) rcp',       tab_cntrl(tab0+9), rcp
   write(*,'(a20, e15.6, e15.6)') '(11) dtphys',   tab_cntrl(tab0+11), dtphys

   write(*,'(a20, f12.2, f12.2)') '(14) year_day', tab_cntrl(tab0+14), year_day
   write(*,'(a20, f12.2, f12.2)') '(15) periheli', tab_cntrl(tab0+15), periheli
   write(*,'(a20, f12.2, f12.2)') '(16) aphelie',  tab_cntrl(tab0+16), aphelie
   write(*,'(a20, f12.2, f12.2)') '(17) peri_day', tab_cntrl(tab0+17), peri_day
   write(*,'(a20, f12.2, f12.2)') '(18) obliquit', tab_cntrl(tab0+18), obliquit

   write(*,'(a20, e15.6, e15.6)') '(19) z0_default', tab_cntrl(tab0+19), z0_default
   write(*,'(a20, e15.6, e15.6)') '(21) emin_turb',  tab_cntrl(tab0+21), emin_turb
   write(*,'(a20, f12.2, f12.2)') '(20) lmixmin',    tab_cntrl(tab0+20), real(lmixmin)

   write(*,'(a20, f12.2, f12.2)') '(26) emissiv',    tab_cntrl(tab0+26), emissiv
   write(*,'(a20, f12.2, f12.2)') '(24) emisice(1)', tab_cntrl(tab0+24), emisice(1)
   write(*,'(a20, f12.2, f12.2)') '(25) emisice(2)', tab_cntrl(tab0+25), emisice(2)
   write(*,'(a20, f12.2, f12.2)') '(22) albedice(1)',tab_cntrl(tab0+22), albedice(1)
   write(*,'(a20, f12.2, f12.2)') '(23) albedice(2)',tab_cntrl(tab0+23), albedice(2)

   write(*,'(a20, e15.6, e15.6)') '(31) iceradius(1)', tab_cntrl(tab0+31), iceradius(1)
   write(*,'(a20, e15.6, e15.6)') '(32) iceradius(2)', tab_cntrl(tab0+32), iceradius(2)

   write(*,'(a20, f12.2, f12.2)') '(33) dtemisice(1)', tab_cntrl(tab0+33), dtemisice(1)
   write(*,'(a20, f12.2, f12.2)') '(34) dtemisice(2)', tab_cntrl(tab0+34), dtemisice(2)

   write(*,'(a20, f12.2, f12.2)') '(27) tauvis',      tab_cntrl(tab0+27), tauvis
   write(*,'(a20, f12.2, f12.2)') '(35) volcapa',      tab_cntrl(tab0+35), volcapa

   write(*,*)
   write(*,*) 'End of tabfi modification section'

!==============================================================================
! Orbital initialization
!==============================================================================

   call iniorbit()

!==============================================================================
! Final compatibility fix (old restart files)
!==============================================================================

   if (iceradius(1) == 0.0) then

      iceradius(1) = 100.0e-6
      iceradius(2) = 100.0e-6

      dtemisice(1) = 0.4
      dtemisice(2) = 0.4

      write(*,*) 'tabfi WARNING: old initialization file detected'
      write(*,*) 'Default iceradius and dtemisice restored'

   end if

end subroutine tabfi

end module tabfi_mod