module phyetat0_mod

implicit none

real, save :: tab_cntrl_mod(100)

!$OMP THREADPRIVATE(tab_cntrl_mod)

!======================================================================
contains
!======================================================================

subroutine phyetat0(fichnom,tab0,Lmodif,nsoil,ngrid,nlay,nq,nqsoil,          &
                    day_ini,time0,tsurf,tsoil,albedo,emis,q2,qsurf,qsoil,    &
                    tauscaling,totcloudfrac,wstar,watercap,perennial_co2ice, &
                    def_slope,def_slope_mean,subslope_dist)

use tracer_mod,          only: noms ! tracer names
use surfdat_h,           only: phisfi, albedodat, z0, z0_default, zmea, zstd, &
                               zsig, zgam, zthe, hmons, summit, base, watercaptag
use iostart,             only: nid_start, open_startphy, close_startphy, &
                               get_field, get_var, inquire_field,        &
                               inquire_dimension, inquire_dimension_length
use nonoro_gwd_ran_mod,  only: du_nonoro_gwd, dv_nonoro_gwd
use nonoro_gwd_mix_mod,  only: du_eddymix_gwd, dv_eddymix_gwd, de_eddymix_rto, &
                               df_eddymix_flx !dr_depflux_gwd
use compute_dtau_mod,    only: dtau
use dust_rad_adjust_mod, only: dust_rad_adjust_prev,dust_rad_adjust_next
use dust_param_mod,      only: dustscaling_mode
use ioipsl_getin_p_mod,  only: getin_p
use comsoil_h,           only: flux_geo
use comslope_mod,        only: nslope, major_slope
use paleoclimate_mod,    only: paleoclimate, h2o_ice_depth, lag_co2_ice, d_coef
use comcstfi_h,          only: pi
use geometry_mod,        only: latitude

implicit none

include "callkeys.h"

!======================================================================
! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
!  Adaptation � Mars : Yann Wanherdrick
! Objet: Lecture de l etat initial pour la physique
! Modifs: Aug.2010 EM : use NetCDF90 to load variables (enables using
!                      r4 or r8 restarts independently of having compiled
!                      the GCM in r4 or r8)
!         June 2013 TN : Possibility to read files with a time axis
!         November 2013 EM : Enabeling parallel, using iostart module
!         March 2020 AD: Enabling initialization of physics without startfi
!                        flag: startphy_file
!======================================================================
integer nbsrf !Mars nbsrf a 1 au lieu de 4
parameter (nbsrf=1) ! nombre de sous-fractions pour une maille
!======================================================================
!  Arguments:
!  ---------
!  inputs:
character(*), intent(in) :: fichnom ! "startfi.nc" file
integer,      intent(in) :: tab0
integer,      intent(in) :: Lmodif
integer,      intent(in) :: nsoil   ! # of soil layers
integer,      intent(in) :: ngrid   ! # of atmospheric columns
integer,      intent(in) :: nlay    ! # of atmospheric layers
integer,      intent(in) :: nq
integer,      intent(in) :: nqsoil  ! # of tracers in the soil
integer, intent(inout) :: day_ini
real,    intent(inout) :: time0

!  outputs:
real, intent(out) :: tsurf(ngrid,nslope)              ! surface temperature
real, intent(out) :: tsoil(ngrid,nsoil,nslope)        ! soil temperature
real, intent(out) :: albedo(ngrid,2,nslope)           ! surface albedo
real, intent(out) :: emis(ngrid,nslope)               ! surface emissivity
real, intent(out) :: q2(ngrid,nlay+1)
real, intent(out) :: qsurf(ngrid,nq,nslope)           ! tracers on surface
real, intent(out) :: qsoil(ngrid,nsoil,nqsoil,nslope) ! tracers in the subsurface
real, intent(out) :: tauscaling(ngrid)                ! dust conversion factor
real, intent(out) :: totcloudfrac(ngrid)              ! total cloud fraction
real, intent(out) :: wstar(ngrid)                     ! Max vertical velocity in thermals (m/s)
real, intent(out) :: watercap(ngrid,nslope)           ! h2o_ice_cover
real, intent(out) :: perennial_co2ice(ngrid,nslope)   ! perennial co2 ice(kg/m^2)
real, intent(out) :: def_slope(nslope+1)              ! boundaries for bining of the slopes
real, intent(out) :: def_slope_mean(nslope)
real, intent(out) :: subslope_dist(ngrid,nslope)      ! undermesh statistics
!======================================================================
!  Local variables:
real :: surffield(ngrid) ! to temporarily store a surface field
real :: xmin, xmax       ! to display min and max of a field

integer :: ig, iq, lmax, islope, nid, nvarid, ierr, i, nsrf, nqold
! integer :: isoil
! integer :: length
! parameter (length=100)
character(7) :: str7
character(2) :: str2
character(1) :: yes

real :: p_rad, p_omeg, p_g, p_mugaz, p_daysec

! flag which identifies if 'startfi.nc' file is using old names (qsurf01,...)
logical       :: oldtracernames = .false.
integer       :: count
character(30) :: txt ! to store some text

! specific for time
real, allocatable :: time(:)   ! times stored in start
integer           :: timelen   ! number of times stored in the file
integer           :: indextime ! index of selected time

integer :: edges(3),corner(3)
logical :: found
real    :: timestart  ! to pick which initial state to start from
real    :: surfemis   ! constant emissivity when no startfi
real    :: surfalbedo ! constant albedo when no startfi
real    :: watercaptag_tmp(ngrid)

! Sub-grid scale slopes
logical           :: startphy_slope ! to be retrocompatible and add the nslope dimension
real, allocatable :: default_def_slope(:)
real              :: sum_dist
real              :: current_max !var to find max distrib slope

! Variables for CO2 index
integer :: igcm_co2_tmp

character(5) :: modname="phyetat0"
!======================================================================

write(*,*) "phyetat0: startphy_file", startphy_file

if (startphy_file) then
   ! open physics initial state file:
   call open_startphy(fichnom)
   ! possibility to modify tab_cntrl in tabfi
   write(*,*)
   write(*,*) 'TABFI in phyeta0: Lmodif=',Lmodif," tab0=",tab0
   call tabfi (nid_start,Lmodif,tab0,day_ini,lmax,p_rad, &
               p_omeg,p_g,p_mugaz,p_daysec,time0)
else ! "academic" initialization of planetary parameters via tabfi
   call tabfi (0,0,0,day_ini,lmax,p_rad, &
               p_omeg,p_g,p_mugaz,p_daysec,time0)
endif ! of if (startphy_file)

if (startphy_file) then
   call get_var("def_slope",def_slope,found)
   if (.not. found) then
       startphy_slope=.false.
       write(*,*)'slope_settings: Problem while reading <def_slope>'
       write(*,*)'Checking that nslope=1'
       if (nslope == 1) then
           write(*,*)'We take default def_slope and subslope dist'
           allocate(default_def_slope(nslope + 1))
           default_def_slope(1) = -90.
           default_def_slope(2) = 90.
           subslope_dist = 1.
           def_slope = default_def_slope
       else
           write(*,*)'Variable def_slope is not present in the start and'
           write(*,*)'you are trying to run with nslope!=1'
           write(*,*)'This is not possible, you have to go through newstart'
       endif
   else
       startphy_slope=.true.
       call get_field("subslope_dist",subslope_dist,found,indextime)
       if (.not. found) then
           write(*,*)'slope_settings: Problem while reading <subslope_dist>'
           write(*,*)'We have to abort.'
           write(*,*)'Please check that nslope is coherent, you can modify this parameter with newstart'
           call abort_physic(modname,"phyetat0: Failed loading <subslope_dist>",1)
       endif
   endif
else ! (no startphy_file)
    if (nslope == 1) then
        allocate(default_def_slope(2))
        default_def_slope = 0.
        subslope_dist = 1.
        def_slope = default_def_slope
    else
        write(*,*)'Without startfi, lets first run with nslope=1'
        call abort_physic(modname,"phyetat0: No startfi and nslope!=1",1)
    endif
endif

do islope = 1,nslope
    def_slope_mean(islope) = (def_slope(islope) + def_slope(islope + 1))/2.
enddo

DO ig = 1,ngrid
    sum_dist = 0.
    DO islope = 1,nslope
        sum_dist = sum_dist + subslope_dist(ig,islope)
    ENDDO
    DO islope = 1,nslope
        subslope_dist(ig,islope) = subslope_dist(ig,islope)/sum_dist
    ENDDO
ENDDO

!Now determine the major subslope, ie. the maximal distribution

DO ig=1,ngrid
  major_slope(ig)=1
  current_max=subslope_dist(ig,1)
  DO islope=2,nslope
    if(subslope_dist(ig,islope).GT.current_max) then
      major_slope(ig)=islope
      current_max=subslope_dist(ig,islope)
    ENDIF
  ENDDO
ENDDO

if (startphy_file) then
   ! Load surface geopotential:
   call get_field("phisfi",phisfi,found)
   if (.not.found) then
     call abort_physic(modname, &
                "phyetat0: Failed loading <phisfi>",1)
   endif
else
  phisfi(:)=0.
endif ! of if (startphy_file)
write(*,*) "phyetat0: surface geopotential <phisfi> range:", &
               minval(phisfi), maxval(phisfi)


if (startphy_file) then
   ! Load bare ground albedo:
   call get_field("albedodat",albedodat,found)
   if (.not.found) then
     call abort_physic(modname, &
                "phyetat0: Failed loading <albedodat>",1)
   endif
else ! If no startfi file, use parameter surfalbedo in def file
  surfalbedo=0.1
  call getin_p("surfalbedo_without_startfi",surfalbedo)
  print*,"surfalbedo_without_startfi",surfalbedo
  albedodat(:)=surfalbedo
endif ! of if (startphy_file)
write(*,*) "phyetat0: Bare ground albedo <albedodat> range:", &
             minval(albedodat), maxval(albedodat)

! ZMEA
if (startphy_file) then
   call get_field("ZMEA",zmea,found)
   if (.not.found) then
     call abort_physic(modname, &
                "phyetat0: Failed loading <ZMEA>",1)
   endif
else
  zmea(:)=0.
endif ! of if (startphy_file)
write(*,*) "phyetat0: <ZMEA> range:", &
            minval(zmea), maxval(zmea)

! ZSTD
if (startphy_file) then
   call get_field("ZSTD",zstd,found)
   if (.not.found) then
     call abort_physic(modname, &
                "phyetat0: Failed loading <ZSTD>",1)
   endif
else
  zstd(:)=0.
endif ! of if (startphy_file)
write(*,*) "phyetat0: <ZSTD> range:", &
            minval(zstd), maxval(zstd)

! ZSIG
if (startphy_file) then
   call get_field("ZSIG",zsig,found)
   if (.not.found) then
     call abort_physic(modname, &
                "phyetat0: Failed loading <ZSIG>",1)
   endif
else
  zsig(:)=0.
endif ! of if (startphy_file)
write(*,*) "phyetat0: <ZSIG> range:", &
            minval(zsig), maxval(zsig)

! ZGAM
if (startphy_file) then
   call get_field("ZGAM",zgam,found)
   if (.not.found) then
     call abort_physic(modname, &
                "phyetat0: Failed loading <ZGAM>",1)
   endif
else
  zgam(:)=0.
endif ! of if (startphy_file)
write(*,*) "phyetat0: <ZGAM> range:", &
            minval(zgam), maxval(zgam)

! ZTHE
if (startphy_file) then
   call get_field("ZTHE",zthe,found)
   if (.not.found) then
     call abort_physic(modname, &
                "phyetat0: Failed loading <ZTHE>",1)
   endif
else
  zthe(:)=0.
endif ! of if (startphy_file)
write(*,*) "phyetat0: <ZTHE> range:", &
             minval(zthe), maxval(zthe)

! hmons
if (startphy_file) then
   call get_field("hmons",hmons,found)
   if (.not.found) then
     write(*,*) "WARNING: phyetat0: Failed loading <hmons>"
     if (rdstorm) then
     call abort_physic(modname, &
                "phyetat0: Failed loading <hmons>",1)
     else
       write(*,*) "will continue anyway..."
       write(*,*) "because you may not need it."
       hmons(:)=0.
     end if ! if (rdstorm)
   else
     do ig=1,ngrid
       if (hmons(ig) .eq. -999999.)  hmons(ig)=0.
     enddo
   endif ! (.not.found)
else
   hmons(:)=0.
endif ! if (startphy_file)
write(*,*) "phyetat0: <hmons> range:", &
            minval(hmons), maxval(hmons)

! summit
if (startphy_file) then
   call get_field("summit",summit,found)
   if (.not.found) then
     write(*,*) "WARNING: phyetat0: Failed loading <summit>"
     if (rdstorm) then
     call abort_physic(modname, &
                "phyetat0: Failed loading <summit>",1)
     else
       write(*,*) "will continue anyway..."
       write(*,*) "because you may not need it."
       summit(:)=0.
     end if
   else
     do ig=1,ngrid
       if (summit(ig) .eq. -999999.)  summit(ig)=0.
     enddo
   endif ! if (.not.found)
else
   summit(:)=0.
endif ! if (startphy_file)
write(*,*) "phyetat0: <summit> range:", &
            minval(summit), maxval(summit)

! base
if (startphy_file) then
   call get_field("base",base,found)
   if (.not.found) then
     write(*,*) "WARNING: phyetat0: Failed loading <base>"
     if (rdstorm) then
     call abort_physic(modname, &
                "phyetat0: Failed loading <base>",1)
     else
       write(*,*) "will continue anyway..."
       write(*,*) "because you may not need it."
       base(:)=0.
     end if
   else
     do ig=1,ngrid
       if (base(ig) .eq. -999999.)  base(ig)=0.
     enddo
   endif ! if(.not.found)
else
   base(:)=0.
endif ! if (startphy_file)
write(*,*) "phyetat0: <base> range:", &
            minval(base), maxval(base)

! Time axis
! obtain timestart from run.def
timestart=-9999 ! default value
call getin_p("timestart",timestart)
if (startphy_file) then
   found=inquire_dimension("Time")
   if (.not.found) then
     indextime = 1
     write(*,*) "phyetat0: No time axis found in "//trim(fichnom)
   else
     write(*,*) "phyetat0: Time axis found in "//trim(fichnom)
     timelen=inquire_dimension_length("Time")
     allocate(time(timelen))
     ! load "Time" array:
     call get_var("Time",time,found)
     if (.not.found) then
     call abort_physic(modname, &
                "phyetat0: Failed loading <Time>",1)
     endif
     ! seclect the desired time index
     IF (timestart .lt. 0) THEN  ! default: we use the last time value
       indextime = timelen
     ELSE  ! else we look for the desired value in the time axis
       indextime = 0
       DO i=1,timelen
         IF (abs(time(i) - timestart) .lt. 0.01) THEN
           indextime = i
           EXIT
         ENDIF
       ENDDO
       IF (indextime .eq. 0) THEN
         PRINT*, "Time", timestart," is not in "//trim(fichnom)//"!!"
         PRINT*, "Stored times are:"
         DO i=1,timelen
            PRINT*, time(i)
         ENDDO
         call abort_physic(modname,"phyetat0: Time error",1)
       ENDIF
     ENDIF ! of IF (timestart .lt. 0)
     ! In startfi the absolute date is day_ini + time0 + time
     ! For now on, in the GCM physics, it is day_ini + time0
     time0 = time(indextime) + time0
     day_ini = day_ini + INT(time0)
     time0 = time0 - INT(time0)
     PRINT*, "phyetat0: Selected time ",time(indextime), &
             " at index ",indextime
     DEALLOCATE(time)
   endif ! of if Time not found in file

   call ini_tab_controle_dyn_xios(day_ini)

else
  indextime = 1
endif ! if (startphy_file)

! Dust conversion factor
if (startphy_file) then
   call get_field("tauscaling",tauscaling,found,indextime)
   if (.not.found) then
     write(*,*) "phyetat0: <tauscaling> not in file"
     tauscaling(:) = 1
   endif
else
   tauscaling(:) = 1
endif ! if (startphy_file)
write(*,*) "phyetat0: dust conversion factor <tauscaling> range:", &
            minval(tauscaling), maxval(tauscaling)

! dust_rad_adjust_* for radiative rescaling of dust
if (dustscaling_mode==2) then
 if (startphy_file) then
   call get_field("dust_rad_adjust_prev",dust_rad_adjust_prev,found,indextime)
   if (.not.found) then
     write(*,*) "phyetat0: <dust_rad_adjust_prev> not in file; set to 1"
     dust_rad_adjust_prev(:) = 1
   endif
   call get_field("dust_rad_adjust_next",dust_rad_adjust_next,found,indextime)
   if (.not.found) then
     write(*,*) "phyetat0: <dust_rad_adjust_next> not in file; set to 1"
     dust_rad_adjust_next(:) = 1
   endif
 else
   dust_rad_adjust_prev(:)= 0
   dust_rad_adjust_next(:)= 0
 endif ! if (startphy_file)
 write(*,*) "phyetat0: radiative scaling coeff <dust_rad_adjust_prev> range:", &
            minval(dust_rad_adjust_prev), maxval(dust_rad_adjust_prev)
 write(*,*) "phyetat0: radiative scaling coeff <dust_rad_adjust_next> range:", &
            minval(dust_rad_adjust_next), maxval(dust_rad_adjust_next)
endif ! of if (dustscaling_mode==2)

! dtau: opacity difference between GCM and dust scenario
if (startphy_file) then
   call get_field("dtau",dtau,found,indextime)
   if (.not.found) then
     write(*,*) "phyetat0: <dtau> not in file; set to zero"
     dtau(:) = 0
   endif
else
   dtau(:)= 0
endif ! if (startphy_file)
write(*,*) "phyetat0: opacity diff wrt scenario <dtau> range:", &
            minval(dtau), maxval(dtau)

! Sub-grid cloud fraction
if (startphy_file) then
   call get_field("totcloudfrac",totcloudfrac,found,indextime)
   if (.not.found) then
     write(*,*) "phyetat0: <totcloudfrac> not in file WARNING put to 1"
     totcloudfrac(:) = 1.0 !valeur par defaut (CLFfixval par defaut)
   endif
else
   totcloudfrac(:)=1.0
endif ! if (startphy_file)
write(*,*) "phyetat0: total cloud fraction <totcloudfrac> range:", &
            minval(totcloudfrac), maxval(totcloudfrac)

! Max vertical velocity in thermals
if (startphy_file) then
   call get_field("wstar",wstar,found,indextime)
   if (.not.found) then
     write(*,*) "phyetat0: <wstar> not in file! Set to zero"
     wstar(:)=0
   endif
else
   wstar(:)=0
endif ! if (startphy_file)
write(*,*) "phyetat0: Max vertical velocity in thermals <wstar> range:", &
            minval(wstar),maxval(wstar)

! Surface temperature :
if (startphy_file) then !tsurf
   call get_field("tsurf",tsurf,found,indextime)
   if (.not.found) then
     call abort_physic(modname, &
                "phyetat0: Failed loading <tsurf>",1)
   endif
else
  tsurf(:,:)=0. ! will be updated afterwards in physiq !
endif ! of if (startphy_file)
write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
            minval(tsurf), maxval(tsurf)

! Surface albedo
if (startphy_file) then
   call get_field("albedo",albedo(:,1,:),found,indextime)
   if (.not.found) then
     write(*,*) "phyetat0: Failed loading <albedo>"
     do islope=1,nslope
       albedo(:,1,islope)=albedodat(:)
     enddo
   endif
else
   do islope=1,nslope
     albedo(:,1,islope)=albedodat(:)
   enddo
endif ! of if (startphy_file)
write(*,*) "phyetat0: Surface albedo <albedo> range:", &
            minval(albedo(:,1,:)), maxval(albedo(:,1,:))
albedo(:,2,:)=albedo(:,1,:)

! Surface emissivity
if (startphy_file) then
   call get_field("emis",emis,found,indextime)
   if (.not.found) then
     call abort_physic(modname, &
                "phyetat0: Failed loading <emis>",1)
   endif
else
  ! If no startfi file, use parameter surfemis in def file
  surfemis=0.95
  call getin_p("surfemis_without_startfi",surfemis)
  print*,"surfemis_without_startfi",surfemis
  emis(:,:)=surfemis
endif ! of if (startphy_file)
write(*,*) "phyetat0: Surface emissivity <emis> range:", &
            minval(emis), maxval(emis)


! surface roughness length (NB: z0 is a common in surfdat_h)
if (startphy_file) then
   call get_field("z0",z0,found)
   if (.not.found) then
     write(*,*) "phyetat0: Failed loading <z0>"
     write(*,*) 'will use constant value of z0_default:',z0_default
     z0(:)=z0_default
   endif
else
   z0(:)=z0_default
endif ! of if (startphy_file)
write(*,*) "phyetat0: Surface roughness <z0> range:", &
            minval(z0), maxval(z0)

! pbl wind variance
if (startphy_file) then
   call get_field("q2",q2,found,indextime)
   if (.not.found) then
     call abort_physic(modname, &
                "phyetat0: Failed loading <q2>",1)
   endif
else
  q2(:,:)=0.
endif ! of if (startphy_file)
write(*,*) "phyetat0: PBL wind variance <q2> range:", &
            minval(q2), maxval(q2)

! Non-orographic gravity waves
if (startphy_file) then
   call get_field("du_nonoro_gwd",du_nonoro_gwd,found,indextime)
   if (.not.found) then
      write(*,*) "phyetat0: <du_nonoro_gwd> not in file"
      du_nonoro_gwd(:,:)=0.
   endif
else
du_nonoro_gwd(:,:)=0.
endif ! if (startphy_file)
write(*,*) "phyetat0: Memory of zonal wind tendency due to non-orographic GW"
write(*,*) " <du_nonoro_gwd> range:", &
             minval(du_nonoro_gwd), maxval(du_nonoro_gwd)

if (startphy_file) then
   call get_field("dv_nonoro_gwd",dv_nonoro_gwd,found,indextime)
   if (.not.found) then
      write(*,*) "phyetat0: <dv_nonoro_gwd> not in file"
      dv_nonoro_gwd(:,:)=0.
   endif
else ! ! if (startphy_file)
dv_nonoro_gwd(:,:)=0.
endif ! if (startphy_file)
write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW"
write(*,*) " <dv_nonoro_gwd> range:", &
             minval(dv_nonoro_gwd), maxval(dv_nonoro_gwd)
if (startphy_file) then
   call get_field("du_eddymix_gwd",du_eddymix_gwd,found,indextime)
   if (.not.found) then
      write(*,*) "phyetat0: <du_eddymix_gwd> not in file"
      du_eddymix_gwd(:,:)=0.
   endif
else
du_eddymix_gwd(:,:)=0.
endif ! if (startphy_file)
write(*,*) "phyetat0: Memory of zonal wind tendency due to non-orographic GW mixing"
write(*,*) " <du_eddymix_gwd> range:", &
             minval(du_eddymix_gwd), maxval(du_eddymix_gwd)

if (startphy_file) then
   call get_field("dv_nonoro_gwd",dv_nonoro_gwd,found,indextime)
   if (.not.found) then
      write(*,*) "phyetat0: <dv_nonoro_gwd> not in file"
      dv_nonoro_gwd(:,:)=0.
   endif
else ! ! if (startphy_file)
dv_nonoro_gwd(:,:)=0.
endif ! if (startphy_file)
write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW"
write(*,*) " <dv_nonoro_gwd> range:", &
             minval(dv_nonoro_gwd), maxval(dv_nonoro_gwd)

if (startphy_file) then
   call get_field("dv_eddymix_gwd",dv_eddymix_gwd,found,indextime)
   if (.not.found) then
      write(*,*) "phyetat0: <dv_eddymix_gwd> not in file"
      dv_eddymix_gwd(:,:)=0.
   endif
else ! ! if (startphy_file)
dv_eddymix_gwd(:,:)=0.
endif ! if (startphy_file)
write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW mixing"
write(*,*) " <dv_eddymix_gwd> range:", &
             minval(dv_eddymix_gwd), maxval(dv_eddymix_gwd)

!if (startphy_file) then
!   call get_field("dr_depflux_gwd",dr_depflux_gwd,found,indextime)
!   if (.not.found) then
!      write(*,*) "phyetat0: <dr_depflux_gwd> not in file"
!      dr_depflux_gwd(:,:,:)=0.
!   endif
!else ! ! if (startphy_file)
!!dr_depflux_gwd(:,:,:)=0.
!endif ! if (startphy_file)
!write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW mixing"
!write(*,*) " <dr_depflux_gwd> range:", &
!             minval(dr_depflux_gwd), maxval(dr_depflux_gwd)

if (startphy_file) then
   call get_field("de_eddymix_rto",de_eddymix_rto,found,indextime)
   if (.not.found) then
      write(*,*) "phyetat0: <de_eddymix_rto> not in file"
      de_eddymix_rto(:,:)=0.
   endif
else ! ! if (startphy_file)
de_eddymix_rto(:,:)=0.
endif ! if (startphy_file)
write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW mixing"
write(*,*) " <de_eddymix_rto> range:", &
             minval(de_eddymix_rto), maxval(de_eddymix_rto)

if (startphy_file) then
   call get_field("df_eddymix_flx ",df_eddymix_flx ,found,indextime)
   if (.not.found) then
      write(*,*) "phyetat0: <df_eddymix_flx > not in file"
      df_eddymix_flx (:,:)=0.
   endif
else ! ! if (startphy_file)
df_eddymix_flx (:,:)=0.
endif ! if (startphy_file)
write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW mixing"
write(*,*) " <df_eddymix_flx > range:", &
             minval(df_eddymix_flx ), maxval(df_eddymix_flx )

! tracer on surface
if (nq.ge.1) then
  do iq=1,nq
    txt=noms(iq)
    if (txt.eq."h2o_vap") then
      ! There is no surface tracer for h2o_vap;
      ! "h2o_ice" should be loaded instead
      txt="h2o_ice"
      write(*,*) 'phyetat0: loading surface tracer', &
                           ' h2o_ice instead of h2o_vap'
      write(*,*) 'iq = ', iq
    endif

    if (hdo) then
      if (txt.eq."hdo_vap") then
        txt="hdo_ice"
        write(*,*) 'phyetat0: loading surface tracer', &
                           ' hdo_ice instead of hdo_vap'
      endif
    endif !hdo

    if (startphy_file) then
      if (txt.eq."co2") then
        ! We first check if co2ice exist in the startfi file (old way)
        ! CO2 ice cover
        call get_field("co2ice",qsurf(:,iq,:),found,indextime)
        ! If not, we load the corresponding tracer in qsurf
        if (.not.found) then
          call get_field(txt,qsurf(:,iq,:),found,indextime)
          if (.not.found) then
            call abort_physic(modname, &
                "phyetat0: Failed loading co2ice. There is neither the variable co2ice nor qsurf!",1)
          endif
        endif
      else ! (not the co2 tracer)
        call get_field(txt,qsurf(:,iq,:),found,indextime)
        if (.not.found) then
          write(*,*) "phyetat0: Failed loading <",trim(txt),">"
          write(*,*) "         ",trim(txt)," is set to zero"
          qsurf(:,iq,:)=0.
        endif
      endif !endif co2
    else !(not startphy_file)
      qsurf(:,iq,:)=0. ! co2ice is set to 0
    endif ! of if (startphy_file)
    write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
                 minval(qsurf(:,iq,:)), maxval(qsurf(:,iq,:))
  enddo ! of do iq=1,nq
endif ! of if (nq.ge.1)

if (startphy_file) then
   call get_field("watercap",watercap,found,indextime)
   if (.not.found) then
     write(*,*) "phyetat0: Failed loading <watercap> : ", &
                          "<watercap> is set to zero"
     watercap(:,:)=0
     write(*,*) 'Now transfer negative surface water ice to', &
                ' watercap !'
     if (nq.ge.1) then
      do iq=1,nq
       txt=noms(iq)
       if (txt.eq."h2o_ice") then
         do ig=1,ngrid
          do islope=1,nslope
            if (qsurf(ig,iq,islope).lt.0.0) then
              watercap(ig,islope) = qsurf(ig,iq,islope)
              qsurf(ig,iq,islope) = 0.0
            end if
          enddo
         end do
       endif

       if (txt.eq."hdo_ice") then
         do ig=1,ngrid
          do islope=1,nslope
            if (qsurf(ig,iq,islope).lt.0.0) then
              qsurf(ig,iq,islope) = 0.0
            end if
          enddo
         end do
       endif

      enddo
     endif ! of if (nq.ge.1)
   endif ! of if (.not.found)
else
   watercap(:,:)=0
endif ! of if (startphy_file)
write(*,*) "phyetat0: Surface water ice <watercap> range:", &
            minval(watercap), maxval(watercap)

if (startphy_file) then
  ! Call to soil_settings, in order to read soil temperatures,
  ! as well as thermal inertia and volumetric heat capacity
  call soil_settings(nid_start,ngrid,nsoil,nqsoil,tsurf,tsoil,qsoil,indextime)
else
    flux_geo(:,:) = 0.
endif ! of if (startphy_file)

if (startphy_file) then
   call get_field("watercaptag",watercaptag_tmp,found,indextime)
   if (.not.found) then
     write(*,*) "phyetat0: Failed loading <watercaptag> : ", &
                          "<watercaptag> is set as defined by icelocationmode in surfini.F"
     watercaptag(:)=.false.
   else
     do ig=1,ngrid
       if(watercaptag_tmp(ig).lt.0.5) then
          watercaptag(ig)=.false.
       else
          watercaptag(ig)=.true.
       endif
     enddo
   endif
endif !startphy_file

if (paleoclimate) then
  do iq=1,nq
   txt=noms(iq)
   if (txt.eq."co2") igcm_co2_tmp = iq
  enddo

  if (startphy_file) then
! Depth of H2O lag
   call get_field("h2o_ice_depth",h2o_ice_depth,found,indextime)
   if (.not.found) then
     write(*,*) "phyetat0: Failed loading <h2o_ice_depth> : ", &
                          "<h2o_ice_depth> is set as -1 (no subsurface ice)"
     h2o_ice_depth(:,:) = -1.
   endif

! Diffusion coeficent
   call get_field("d_coef",d_coef,found,indextime)
   if (.not.found) then
     write(*,*) "phyetat0: Failed loading <d_coef> : ", &
                          "<d_coef> is set as 4e-4 (defualt value)"
     d_coef(:,:) = 4e-4
   endif

! Depth of CO2 lag
   call get_field("lag_co2_ice",lag_co2_ice,found,indextime)
   if (.not.found) then
     write(*,*) "phyetat0: Failed loading <lag_co2_ice> : ", &
                          "<lag_co2_ice> is set as -1 (no subsurface ice)"
     lag_co2_ice(:,:) = -1.
   endif

! Perennial CO2 ice
    perennial_co2ice = 0.
    call get_field("perennial_co2ice",perennial_co2ice,found,indextime)
    if (.not. found) then
        write(*,*) "phyetat0: Failed loading <perennial_co2ice> : ", &
                   "<perennial_co2ice> is set as 10m at the South Pole"
        if (abs(latitude(ngrid) - (-pi/2.)) < 1.e-5) then
            perennial_co2ice(ngrid,:) = 10*1.6e3 ! 10m which is convert to kg/m^2
            qsurf(ngrid,igcm_co2_tmp,:) = qsurf(ngrid - 1,igcm_co2_tmp,:) + perennial_co2ice(ngrid,:) ! perennial ice + frost
        endif
    endif ! not found
  else ! no startfiphyle
    h2o_ice_depth = -1.
    lag_co2_ice = -1.
    d_coef = 4.e-4
    perennial_co2ice = 0.
    if (abs(latitude(ngrid) - (-pi/2.)) < 1.e-5) then
        perennial_co2ice(ngrid,:) = 10*1.6e3 ! 10m which is convert to kg/m^2
        qsurf(ngrid,igcm_co2_tmp,:) = qsurf(ngrid - 1,igcm_co2_tmp,:) + perennial_co2ice(ngrid,:) ! perennial ice + frost
    endif
  endif !startphy_file
else
   h2o_ice_depth = -1.
   lag_co2_ice = -1.
   d_coef = 4.e-4
   perennial_co2ice = 0.
endif !paleoclimate

! close file:
if (startphy_file) call close_startphy

end subroutine phyetat0

!======================================================================
subroutine ini_tab_controle_dyn_xios(idayref)

use comcstfi_h,        only: g, mugaz, omeg, rad, rcp
use time_phylmdz_mod,  only: hour_ini, daysec, dtphys
use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat, nbp_lev

implicit none

integer*4, intent(in) :: idayref ! date (initial date for this run)

integer :: length, l
parameter (length = 100)
real    :: tab_cntrl(length) ! run parameters are stored in this array

do l = 1,length
    tab_cntrl(l) = 0.
enddo
tab_cntrl(1)  = real(nbp_lon)
tab_cntrl(2)  = real(nbp_lat-1)
tab_cntrl(3)  = real(nbp_lev)
tab_cntrl(4)  = real(idayref)
tab_cntrl(5)  = rad
tab_cntrl(6)  = omeg
tab_cntrl(7)  = g
tab_cntrl(8)  = mugaz
tab_cntrl(9)  = rcp
tab_cntrl(10) = daysec
tab_cntrl(11) = dtphys

tab_cntrl(27) = hour_ini

tab_cntrl_mod = tab_cntrl

end subroutine ini_tab_controle_dyn_xios

end module phyetat0_mod