MODULE init_testphys1d_mod

implicit none

contains

SUBROUTINE init_testphys1d(ngrid,nq,nlayer,odpref,ndt,ptif,pks,dttestphys,startfiles_1D,q,zqsat,qsurf,dq,dqdyn, &
                           day0,day,time,psurf,tsurf,gru,grv,u,v,w,q2,play,plev,tsoil,temp,albedo,emis,         &
                           latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau)

use ioipsl_getincom,          only: getin ! To use 'getin'
use comcstfi_h,               only: pi, rad, omeg, g, mugaz, rcp, r, cpp
use time_phylmdz_mod,         only: daysec, day_step, ecritphy, iphysiq
use planete_h,                only: year_day, periheli, aphelie, peri_day, obliquit, emin_turb, lmixmin
use surfdat_h,                only: albedodat, z0_default, z0, emissiv, emisice, albedice, iceradius, dtemisice, &
                                    zmea, zstd, zsig, zgam, zthe, hmons, summit, base, phisfi, watercaptag, watercap
use infotrac,                 only: nqtot, tname, nqperes, nqfils
use read_profile_mod,         only: read_profile
use iostart,                  only: open_startphy, get_var, close_startphy
use physics_distribution_mod, only: init_physics_distribution
use comsoil_h,                only: volcapa, nsoilmx, inertiesoil, inertiedat, layer, mlayer, flux_geo
use comvert_mod,              only: ap, bp, aps, bps, pa, preff, presnivs, pseudoalt, scaleheight
use dimradmars_mod,           only: tauvis, totcloudfrac
use regular_lonlat_mod,       only: init_regular_lonlat
use mod_interface_dyn_phys,   only: init_interface_dyn_phys
use geometry_mod,             only: init_geometry
use dimphy,                   only: init_dimphy
use comgeomfi_h,              only: sinlat, ini_fillgeom
use slope_mod,                only: theta_sl, psi_sl
use comslope_mod,             only: def_slope, subslope_dist
use dust_param_mod,           only: tauscaling
use tracer_mod,               only: igcm_co2
use logic_mod,                only: hybrid
use vertical_layers_mod,      only: init_vertical_layers
use inichim_newstart_mod,     only: inichim_newstart
use mod_grid_phy_lmdz,        only: regular_lonlat
use phys_state_var_init_mod,  only: phys_state_var_init
! Mostly for XIOS outputs:
use mod_const_mpi,            only: COMM_LMDZ

implicit none

include "dimensions.h"
include "callkeys.h"

!=======================================================================
! Arguments
!=======================================================================
integer, intent(in) :: ngrid, nlayer
real,    intent(in) :: odpref ! DOD reference pressure (Pa)

integer, intent(inout) :: nq

integer,                           intent(out) :: ndt
real,                              intent(out) :: ptif, pks
real,                              intent(out) :: dttestphys    ! testphys1d timestep
logical,                           intent(out) :: startfiles_1D ! Use of "start1D.txt" and "startfi.nc" files
real, dimension(:,:), allocatable, intent(out) :: q             ! tracer mixing ratio (e.g. kg/kg)
real, dimension(:),   allocatable, intent(out) :: zqsat         ! useful for (atm_wat_profile=2)
real, dimension(:),   allocatable, intent(out) :: qsurf         ! tracer surface budget (e.g. kg.m-2)
real, dimension(:,:), allocatable, intent(out) :: dq, dqdyn     ! Physical and dynamical tandencies
integer,                           intent(out) :: day0          ! initial (sol ; =0 at Ls=0) and final date
real,                              intent(out) :: day           ! date during the run
real,                              intent(out) :: time          ! time (0<time<1 ; time=0.5 a midi)
real,                              intent(out) :: psurf         ! Surface pressure
real, dimension(1),                intent(out) :: tsurf         ! Surface temperature
real,                              intent(out) :: gru, grv      ! prescribed "geostrophic" background wind
real, dimension(nlayer),           intent(out) :: u, v, w       ! zonal, meridional wind
real, dimension(nlayer + 1),       intent(out) :: q2            ! Turbulent Kinetic Energy
real, dimension(nlayer),           intent(out) :: play          ! Pressure at the middle of the layers (Pa)
real, dimension(nlayer + 1),       intent(out) :: plev          ! intermediate pressure levels (pa)
real, dimension(nsoilmx),          intent(out) :: tsoil         ! subsurface soik temperature (K)
real, dimension(nlayer),           intent(out) :: temp          ! temperature at the middle of the layers
real, dimension(1,1),              intent(out) :: albedo        ! surface albedo
real, dimension(1),                intent(out) :: emis          ! surface layer
real, dimension(1),                intent(out) :: latitude, longitude, cell_area
real,                              intent(out) :: atm_wat_profile, atm_wat_tau

!=======================================================================
! Local variables
!=======================================================================
integer                   :: ierr, iq, ilayer, isoil, nlevel, nsoil, flagthermo, flagh2o
integer                   :: dayn ! Final date
real, dimension(nlayer)   :: zlay ! altitude estimee dans les couches (km)
real, dimension(0:nlayer) :: tmp1, tmp2

! Dummy variables along "dynamics scalar grid"
real, dimension(:,:,:,:), allocatable :: qdyn
real, dimension(:,:),     allocatable :: psdyn

! RV & JBC: Use of "start1D.txt" and "startfi.nc" files
logical                :: found, therestart1D, therestartfi
character(len = 30)    :: header
real, dimension(100)   :: tab_cntrl
real, dimension(1,2,1) :: albedo_read ! surface albedo

! New flag to compute paleo orbital configurations + few variables JN
logical :: paleomars
real    :: halfaxe, eccentric, Lsperi

! MVals: isotopes as in the dynamics (CRisi)
integer                                        :: ifils, ipere, generation, ierr0
character(len = 30), dimension(:), allocatable :: tnom_transp ! transporting fluid short name
character(len = 80)                            :: line ! to store a line of text
logical                                        :: continu, there

! LL: Possibility to add subsurface ice
real    :: ice_depth          ! depth of the ice table, ice_depth < 0. means no ice
real    :: inertieice = 2100. ! ice thermal inertia
integer :: iref

! LL: Subsurface geothermal flux
real :: flux_geo_tmp

!=======================================================================
! Code
!=======================================================================

!------------------------------------------------------
! Loading run parameters from "run.def" file
!------------------------------------------------------
! check if 'run.def' file is around (otherwise reading parameters
! from callphys.def via getin() routine won't work.
inquire(file = 'run.def',exist = there)
if (.not. there) then
    write(*,*) 'Cannot find required file "run.def"'
    write(*,*) '  (which should contain some input parameters along with the following line: INCLUDEDEF=callphys.def)'
    write(*,*) ' ... might as well stop here ...'
    stop
endif

write(*,*)'Do you want to use "start1D.txt" and "startfi.nc" files?'
startfiles_1D = .false.
call getin("startfiles_1D",startfiles_1D)
write(*,*) " startfiles_1D = ", startfiles_1D

if (startfiles_1D) then
    inquire(file = 'start1D.txt',exist = therestart1D)
    if (.not. therestart1D) then
        write(*,*) 'There is no "start1D.txt" file!'
        write(*,*) 'Initialization is done with default values.'
    endif
    inquire(file = 'startfi.nc',exist = therestartfi)
    if (.not. therestartfi) then
        write(*,*) 'There is no "startfi.nc" file!'
        write(*,*) 'Initialization is done with default values.'
    endif
endif

!------------------------------------------------------
! Prescribed constants to be set here
!------------------------------------------------------
pi = 2.*asin(1.)

! Mars planetary constants
! ------------------------
rad = 3397200.                   ! mars radius (m) ~3397200 m
omeg = 4.*asin(1.)/(daysec)      ! rotation rate (rad.s-1)
g = 3.72                         ! gravity (m.s-2) ~3.72
mugaz = 43.49                    ! atmosphere mola mass (g.mol-1) ~43.49
rcp = .256793                    ! = r/cp ~0.256793
r = 8.314511*1000./mugaz
cpp = r/rcp
daysec = 88775.                  ! length of a sol (s) ~88775 s
year_day = 669                   ! length of year (sols) ~668.6
periheli = 206.66                ! minimum sun-mars distance (Mkm) ~206.66
aphelie = 249.22                 ! maximum sun-mars distance (Mkm) ~249.22
halfaxe = 227.94                 ! demi-grand axe de l'ellipse
peri_day = 485.                  ! perihelion date (sols since N. Spring)
obliquit = 25.2                  ! Obliquity (deg) ~25.2
eccentric = 0.0934               ! Eccentricity (0.0934)

! Planetary Boundary Layer and Turbulence parameters
! --------------------------------------------------
z0_default = 1.e-2 ! surface roughness (m) ~0.01
emin_turb = 1.e-6  ! minimal turbulent energy ~1.e-8
lmixmin = 30       ! mixing length ~100

! cap properties and surface emissivities
! ---------------------------------------
emissiv = 0.95         ! Bare ground emissivity ~.95
emisice(1) = 0.95      ! Northern cap emissivity
emisice(2) = 0.95      ! Southern cap emisssivity
albedice(1) = 0.5      ! Northern cap albedo
albedice(2) = 0.5      ! Southern cap albedo
iceradius(1) = 100.e-6 ! mean scat radius of CO2 snow (north)
iceradius(2) = 100.e-6 ! mean scat radius of CO2 snow (south)
dtemisice(1) = 2.      ! time scale for snow metamorphism (north)
dtemisice(2) = 2.      ! time scale for snow metamorphism (south)

! mesh surface (not a very usefull quantity in 1D)
! ------------------------------------------------
cell_area(1) = 1.

! check if there is a 'traceur.def' file and process it
! load tracer names from file 'traceur.def'
open(90,file = 'traceur.def',status = 'old',form = 'formatted',iostat = ierr)
if (ierr /= 0) then
    write(*,*) 'Cannot find required file "traceur.def"'
    write(*,*) ' If you want to run with tracers, I need it'
    write(*,*) ' ... might as well stop here ...'
    stop
else
    write(*,*) "testphys1d: Reading file traceur.def"
    ! read number of tracers:
    read(90,*,iostat = ierr) nq
    nqtot = nq ! set value of nqtot (in infotrac module) as nq
    if (ierr /= 0) then
        write(*,*) "testphys1d: error reading number of tracers"
        write(*,*) "   (first line of traceur.def) "
        stop
    endif
    if (nq < 1) then
        write(*,*) "testphys1d: error number of tracers"
        write(*,*) "is nq=",nq," but must be >=1!"
        stop
    endif
endif
! allocate arrays:
allocate(tname(nq),q(nlayer,nq),zqsat(nlayer),qsurf(nq))
allocate(dq(nlayer,nq),dqdyn(nlayer,nq),tnom_transp(nq))

! read tracer names from file traceur.def
do iq = 1,nq
    read(90,'(80a)',iostat = ierr) line ! store the line from traceur.def
    if (ierr /= 0) then
        write(*,*) 'testphys1d: error reading tracer names...'
        stop
    endif
    ! if format is tnom_0, tnom_transp (isotopes)
    read(line,*,iostat = ierr0) tname(iq),tnom_transp(iq)
    if (ierr0 /= 0) then
        read(line,*) tname(iq)
        tnom_transp(iq) = 'air'
    endif
enddo
close(90)

! Isotopes: as in the 3D case we have to determine father/son relations for isotopes and carrying fluid
allocate(nqfils(nqtot))
nqperes = 0
nqfils(:) = 0
do iq = 1,nqtot
    if (tnom_transp(iq) == 'air') then
        ! ceci est un traceur père
        write(*,*) 'Le traceur',iq,', appele ',trim(tname(iq)),', est un pere'
        nqperes = nqperes + 1
    else !if (tnom_transp(iq) == 'air') then
        ! ceci est un fils. Qui est son père?
        write(*,*) 'Le traceur',iq,', appele ',trim(tname(iq)),', est un fils'
        continu = .true.
        ipere = 1
        do while (continu)
            if (tnom_transp(iq) == tname(ipere)) then
                ! Son père est ipere
                write(*,*) 'Le traceur',iq,'appele ',trim(tname(iq)),' est le fils de ',ipere,'appele ',trim(tname(ipere))
                nqfils(ipere) = nqfils(ipere) + 1
                continu = .false.
            else !if (tnom_transp(iq) == tnom_0(ipere)) then
                ipere = ipere + 1
                if (ipere > nqtot) then
                    write(*,*) 'Le traceur',iq,'appele ',trim(tname(iq)),', est orpelin.'
                    call abort_gcm('infotrac_init','Un traceur est orphelin',1)
                endif !if (ipere > nqtot) then
            endif !if (tnom_transp(iq) == tnom_0(ipere)) then
        enddo !do while (continu)
    endif !if (tnom_transp(iq) == 'air') then
enddo !do iq=1,nqtot
write(*,*) 'nqperes=',nqperes
write(*,*) 'nqfils=',nqfils

! Initialize tracers here:
write(*,*) "testphys1d: initializing tracers"
if (.not. therestart1D) then
    call read_profile(nq,nlayer,qsurf,q)
else
    do iq = 1,nq
        open(3,file = 'start1D.txt',status = "old",action = "read")
        read(3,*) header, qsurf(iq),(q(ilayer,iq), ilayer = 1,nlayer)
        if (trim(tname(iq)) /= trim(header)) then
            write(*,*) 'Tracer names not compatible for initialization with "start1D.txt"!'
            stop
        endif
    enddo
endif

call init_physics_distribution(regular_lonlat,4,1,1,1,nlayer,COMM_LMDZ)

! Date and local time at beginning of run
! ---------------------------------------
if (.not. startfiles_1D) then
    ! Date (in sols since spring solstice) at beginning of run
    day0 = 0 ! default value for day0
    write(*,*) 'Initial date (in martian sols ; =0 at Ls=0)?'
    call getin("day0",day0)
    day = float(day0)
    write(*,*) " day0 = ",day0
    ! Local time at beginning of run
    time = 0 ! default value for time
    write(*,*)'Initial local time (in hours, between 0 and 24)?'
    call getin("time",time)
    write(*,*)" time = ",time
    time = time/24. ! convert time (hours) to fraction of sol
else
    call open_startphy("startfi.nc")
    call get_var("controle",tab_cntrl,found)
    if (.not. found) then
        call abort_physic("open_startphy","tabfi: Failed reading <controle> array",1)
    else
        write(*,*)'tabfi: tab_cntrl',tab_cntrl
    endif
    day0 = tab_cntrl(3)
    day = float(day0)

    call get_var("Time",time,found)
    call close_startphy
endif !startfiles_1D

! Discretization (Definition of grid and time steps)
! --------------------------------------------------
nlevel = nlayer + 1
nsoil = nsoilmx

day_step = 48 ! default value for day_step
write(*,*)'Number of time steps per sol?'
call getin("day_step",day_step)
write(*,*) " day_step = ",day_step

ecritphy = day_step ! default value for ecritphy, output every time step

ndt = 10 ! default value for ndt
write(*,*)'Number of sols to run?'
call getin("ndt",ndt)
write(*,*) " ndt = ",ndt

dayn = day0 + ndt
ndt = ndt*day_step
dttestphys = daysec/day_step

! Imposed surface pressure
! ------------------------
psurf = 610. ! Default value for psurf
write(*,*) 'Surface pressure (Pa)?'
if (.not. therestart1D) then
    call getin("psurf",psurf)
else
    read(3,*) header, psurf
endif
write(*,*) " psurf = ",psurf

! Reference pressures
pa = 20.     ! transition pressure (for hybrid coord.)
preff = 610. ! reference surface pressure

! Aerosol properties
! ------------------
tauvis = 0.2 ! default value for tauvis (dust opacity)
write(*,'("Reference dust opacity at ",f4.0," Pa?")')odpref
call getin("tauvis",tauvis)
write(*,*) " tauvis = ",tauvis

! Orbital parameters
! ------------------
if (.not. startfiles_1D) then
    paleomars = .false. ! Default: no water ice reservoir
    call getin("paleomars",paleomars)
    if (paleomars) then
        write(*,*) "paleomars=", paleomars
        write(*,*) "Orbital parameters from callphys.def"
        write(*,*) "Enter eccentricity & Lsperi"
        write(*,*) 'Martian eccentricity (0<e<1)?'
        call getin('eccentric ',eccentric)
        write(*,*)"eccentric =",eccentric
        write(*,*) 'Solar longitude of perihelion (0<Ls<360)?'
        call getin('Lsperi',Lsperi )
        write(*,*)"Lsperi=",Lsperi
        Lsperi = Lsperi*pi/180.0 ! Put it in rad for peri_day
        periheli = halfaxe*(1 - eccentric)
        aphelie = halfaxe*(1 + eccentric)
        call call_dayperi(Lsperi,eccentric,peri_day,year_day)
        write(*,*) "Corresponding orbital params for GCM"
        write(*,*) " periheli = ",periheli
        write(*,*) " aphelie = ",aphelie
        write(*,*) "date of perihelion (sol)",peri_day
    else
        write(*,*) "paleomars=", paleomars
        write(*,*) "Default present-day orbital parameters"
        write(*,*) "Unless specified otherwise"
        write(*,*)'Min. distance Sun-Mars (Mkm)?'
        call getin("periheli",periheli)
        write(*,*) " periheli = ",periheli

        write(*,*)'Max. distance Sun-Mars (Mkm)?'
        call getin("aphelie",aphelie)
        write(*,*) " aphelie = ",aphelie

        write(*,*)'Day of perihelion?'
        call getin("periday",peri_day)
        write(*,*) " periday = ",peri_day

        write(*,*)'Obliquity?'
        call getin("obliquit",obliquit)
        write(*,*) " obliquit = ",obliquit
     endif
endif !(.not. startfiles_1D )

! Latitude/longitude
! ------------------
latitude = 0. ! default value for latitude
write(*,*)'latitude (in degrees)?'
call getin("latitude",latitude(1))
write(*,*) " latitude = ",latitude
latitude = latitude*pi/180.
longitude = 0.
longitude = longitude*pi/180.

! Some initializations (some of which have already been
! done above!) and loads parameters set in callphys.def
! and allocates some arrays
! Mars possible matter with dttestphys in input and include!!!
! Initializations below should mimick what is done in iniphysiq for 3D GCM
call init_interface_dyn_phys
call init_regular_lonlat(1,1,longitude,latitude,(/0.,0./),(/0.,0./))
call init_geometry(1,longitude,latitude,(/0.,0.,0.,0./),(/0.,0.,0.,0./),cell_area)
! Ehouarn: init_vertial_layers called later (because disvert not called yet)
!      call init_vertical_layers(nlayer,preff,scaleheight,
!     &                      ap,bp,aps,bps,presnivs,pseudoalt)
call init_dimphy(1,nlayer) ! Initialize dimphy module
call phys_state_var_init(1,llm,nq,tname,day0,dayn,time,daysec,dttestphys,rad,g,r,cpp,nqperes,nqfils)! MVals: variables isotopes
call ini_fillgeom(1,latitude,longitude,(/1.0/))
call conf_phys(1,llm,nq)

! In 1D model physics are called every time step
! ovverride iphysiq value that has been set by conf_phys
if (iphysiq /= 1) then
    write(*,*) "testphys1d: setting iphysiq=1"
    iphysiq = 1
endif

! Initialize albedo / soil thermal inertia
! ----------------------------------------
if (.not. startfiles_1D) then
    albedodat(1) = 0.2 ! default value for albedodat
    write(*,*)'Albedo of bare ground?'
    call getin("albedo",albedodat(1))
    write(*,*) " albedo = ",albedodat(1)
    albedo(1,1) = albedodat(1)

    inertiedat(1,1) = 400 ! default value for inertiedat
    write(*,*)'Soil thermal inertia (SI)?'
    call getin("inertia",inertiedat(1,1))
    write(*,*) " inertia = ",inertiedat(1,1)

    ice_depth = -1 ! default value: no ice
    call getin("subsurface_ice_depth",ice_depth)

    z0(1) = z0_default ! default value for roughness
    write(*,*) 'Surface roughness length z0 (m)?'
    call getin("z0",z0(1))
    write(*,*) " z0 = ",z0(1)
endif !(.not. startfiles_1D )

! Initialize local slope parameters (only matters if "callslope"
! is .true. in callphys.def)
! slope inclination angle (deg) 0: horizontal, 90: vertical
theta_sl(1) = 0. ! default: no inclination
call getin("slope_inclination",theta_sl(1))
! slope orientation (deg)
! 0 == Northward, 90 == Eastward, 180 == Southward, 270 == Westward
psi_sl(1) = 0. ! default value
call getin("slope_orientation",psi_sl(1))

! Sub-slopes parameters (assuming no sub-slopes distribution for now).
def_slope(1) = -90 ! minimum slope angle
def_slope(2) = 90  ! maximum slope angle
subslope_dist(1,1) = 1 ! fraction of subslopes in mesh

! For the gravity wave scheme
! ---------------------------
zmea(1) = 0.
zstd(1) = 0.
zsig(1) = 0.
zgam(1) = 0.
zthe(1) = 0.

! For the slope wind scheme
! -------------------------
hmons(1) = 0.
write(*,*)'hmons is initialized to ',hmons(1)
summit(1) = 0.
write(*,*)'summit is initialized to ',summit(1)
base(1) = 0.

! Default values initializing the coefficients calculated later
! -------------------------------------------------------------
tauscaling(1) = 1.   ! calculated in aeropacity_mod.F
totcloudfrac(1) = 1. ! calculated in watercloud_mod.F

! Specific initializations for "physiq"
! -------------------------------------
! surface geopotential is not used (or useful) since in 1D
! everything is controled by surface pressure
phisfi(1) = 0.

! Initialization to take into account prescribed winds
! ----------------------------------------------------
ptif = 2.*omeg*sinlat(1)

! Geostrophic wind
gru = 10. ! default value for gru
write(*,*)'zonal eastward component of the geostrophic wind (m/s)?'
call getin("u",gru)
write(*,*) " u = ",gru
grv = 0. !default value for grv
write(*,*)'meridional northward component of the geostrophic wind (m/s)?'
call getin("v",grv)
write(*,*) " v = ",grv

! Initialize winds for first time step
if (.not. therestart1D) then
    u(:) = gru
    v(:) = grv
else
    read(3,*) header, (u(ilayer), ilayer = 1,nlayer)
    read(3,*) header, (v(ilayer), ilayer = 1,nlayer)
endif
w = 0. ! default: no vertical wind

! Initialize turbulente kinetic energy
q2 = 0.

! CO2 ice on the surface
! ----------------------
! get the index of co2 tracer (not known at this stage)
igcm_co2 = 0
do iq = 1,nq
    if (trim(tname(iq)) == "co2") igcm_co2 = iq
enddo
if (igcm_co2 == 0) then
    write(*,*) "testphys1d error, missing co2 tracer!"
    stop
endif

if (.not. startfiles_1D) then
    qsurf(igcm_co2) = 0. ! default value for co2ice
    write(*,*)'Initial CO2 ice on the surface (kg.m-2)'
    call getin("co2ice",qsurf(igcm_co2))
    write(*,*) " co2ice = ",qsurf(igcm_co2)
endif !(.not. startfiles_1D )

! emissivity
! ----------
if (.not. startfiles_1D) then
    emis = emissiv
    if (qsurf(igcm_co2) == 1.) then
        emis = emisice(1) ! northern hemisphere
        if (latitude(1) < 0) emis = emisice(2) ! southern hemisphere
    endif
endif !(.not. startfiles_1D )

! Compute pressures and altitudes of atmospheric levels
! -----------------------------------------------------
! Vertical Coordinates
! """"""""""""""""""""
hybrid = .true.
write(*,*)'Hybrid coordinates?'
call getin("hybrid",hybrid)
write(*,*) " hybrid = ", hybrid

call disvert_noterre
! Now that disvert has been called, initialize module vertical_layers_mod
call init_vertical_layers(nlayer,preff,scaleheight,ap,bp,aps,bps,presnivs,pseudoalt)

plev(:) = ap(:) + psurf*bp(:)
play(:) = aps(:) + psurf*bps(:)
zlay(:) = -200.*r*log(play(:)/plev(1))/g


! Initialize temperature profile
! ------------------------------
pks = psurf**rcp

! Altitude in km in profile: divide zlay by 1000
tmp1(0) = 0.
tmp1(1:) = zlay(:)/1000.

call profile(nlayer + 1,tmp1,tmp2)

if (.not. therestart1D) then
    tsurf = tmp2(0)
    temp(:) = tmp2(1:)
else
    read(3,*) header, tsurf, (temp(ilayer), ilayer = 1,nlayer)
    close(3)
endif

! Initialize soil properties and temperature
! ------------------------------------------
volcapa = 1.e6 ! volumetric heat capacity

if (.not. startfiles_1D) then
    ! Initialize depths
    ! -----------------
    do isoil = 1,nsoil
        layer(isoil) = 2.e-4*(2.**(isoil - 1)) ! layer depth
    enddo

    ! Creating the new soil inertia table if there is subsurface ice:
    if (ice_depth > 0) then
        iref = 1 ! ice/regolith boundary index
        if (ice_depth < layer(1)) then
            inertiedat(1,1) = sqrt(layer(1)/((ice_depth/inertiedat(1,1)**2) + ((layer(1) - ice_depth)/inertieice**2)))
            inertiedat(1,2:) = inertieice
        else ! searching for the ice/regolith boundary:
            do isoil = 1,nsoil
                if ((ice_depth >= layer(isoil)) .and. (ice_depth < layer(isoil + 1))) then
                    iref = isoil + 1
                    exit
                endif
            enddo
            ! We then change the soil inertia table:
            inertiedat(1,:iref - 1) = inertiedat(1,1)
            ! We compute the transition in layer(iref)
            inertiedat(1,iref) = sqrt((layer(iref) - layer(iref - 1))/(((ice_depth - layer(iref - 1))/inertiedat(1,1)**2) + ((layer(iref) - ice_depth)/inertieice**2)))
            ! Finally, we compute the underlying ice:
            inertiedat(1,iref + 1:) = inertieice
        endif ! (ice_depth < layer(1))
    else ! ice_depth < 0 all is set to surface thermal inertia
        inertiedat(1,:) = inertiedat(1,1) ! soil thermal inertia
    endif ! ice_depth > 0

    inertiesoil(1,:,1) = inertiedat(1,:)

    tsoil(:) = tsurf(1) ! soil temperature
endif !(.not. startfiles_1D)

flux_geo_tmp = 0.
call getin("flux_geo",flux_geo_tmp)
flux_geo(:,:) = flux_geo_tmp

! Initialize depths
! -----------------
do isoil = 0,nsoil - 1
    mlayer(isoil) = 2.e-4*(2.**(isoil - 0.5)) ! mid-layer depth
    layer(isoil + 1) = 2.e-4*(2.**isoil)      ! layer depth
enddo

! Initialize traceurs
! -------------------
if (photochem .or. callthermos) then
    write(*,*) 'Initializing chemical species'
    ! flagthermo=0: initialize over all atmospheric layers
    flagthermo = 0
    ! check if "h2o_vap" has already been initialized
    ! (it has been if there is a "profile_h2o_vap" file around)
    inquire(file = "profile_h2o_vap",exist = there)
    if (there) then
        flagh2o = 0 ! 0: do not initialize h2o_vap
    else
        flagh2o = 1 ! 1: initialize h2o_vap in inichim_newstart
    endif

    ! hack to accomodate that inichim_newstart assumes that
    ! q & psurf arrays are on the dynamics scalar grid
    allocate(qdyn(2,1,llm,nq),psdyn(2,1))
    qdyn(1,1,1:llm,1:nq) = q(1:llm,1:nq)
    psdyn(1:2,1) = psurf
    call inichim_newstart(ngrid,nq,qdyn,qsurf,psdyn,flagh2o,flagthermo)
    q(1:llm,1:nq) = qdyn(1,1,1:llm,1:nq)
endif

! Check if the surface is a water ice reservoir
! ---------------------------------------------
if (.not. startfiles_1D) watercap(1,:) = 0 ! Initialize watercap
watercaptag(1) = .false. ! Default: no water ice reservoir
write(*,*)'Water ice cap on ground?'
call getin("watercaptag",watercaptag)
write(*,*) " watercaptag = ",watercaptag

! Check if the atmospheric water profile is specified
! ---------------------------------------------------
! Adding an option to force atmospheric water values JN
atm_wat_profile = -1. ! Default: free atm wat profile
if (water) then
    write(*,*)'Force atmospheric water vapor profile?'
    call getin('atm_wat_profile',atm_wat_profile)
    write(*,*) 'atm_wat_profile = ', atm_wat_profile
    if (abs(atm_wat_profile + 1.) < 1.e-15) then ! if == -1.
        write(*,*) 'Free atmospheric water vapor profile'
        write(*,*) 'Total water is conserved in the column'
    else if (abs(atm_wat_profile) < 1.e-15) then ! if == 0.
        write(*,*) 'Dry atmospheric water vapor profile'
    else if (0. < atm_wat_profile .and. atm_wat_profile <= 1.) then
        write(*,*) 'Prescribed atmospheric water vapor profile'
        write(*,*) 'Unless it reaches saturation (maximal value)'
    else
        write(*,*) 'Water vapor profile value not correct!'
        stop
    endif
endif

! Check if the atmospheric water profile relaxation is specified
! --------------------------------------------------------------
! Adding an option to relax atmospheric water values JBC
atm_wat_tau = -1. ! Default: no time relaxation
if (water) then
    write(*,*) 'Relax atmospheric water vapor profile?'
    call getin('atm_wat_tau',atm_wat_tau)
    write(*,*) 'atm_wat_tau = ', atm_wat_tau
    if (atm_wat_tau < 0.) then
        write(*,*) 'Atmospheric water vapor profile is not relaxed.'
    else
        if (0. <= atm_wat_profile .and. atm_wat_profile <= 1.) then
            write(*,*) 'Relaxed atmospheric water vapor profile towards ', atm_wat_profile
            write(*,*) 'Unless it reaches saturation (maximal value)'
        else
            write(*,*) 'Reference atmospheric water vapor profile not known!'
            write(*,*) 'Please, specify atm_wat_profile'
            stop
        endif
    endif
endif

END SUBROUTINE init_testphys1d

END MODULE init_testphys1d_mod