MODULE init_testphys1d_mod

implicit none

contains

SUBROUTINE init_testphys1d(start1Dname,startfiname,therestart1D,therestartfi,ngrid,nlayer,odpref,               &
                           nq,q,time,psurf,u,v,temp,ndt,ptif,pks,dttestphys,zqsat,dq,dqdyn,day0,day,gru,grv,w,  &
                           play,plev,latitude,longitude,cell_area,ctrl_h2ovap,relaxtime_h2ovap,qref_h2ovap,      &
                           ctrl_h2oice,relaxtime_h2oice,qref_h2oice)

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, steps_per_sol, outputs_per_sol, 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, &
                                    tsurf, emis, qsurf, perennial_co2ice, ini_surfdat_h_slope_var, end_surfdat_h_slope_var
use infotrac,                 only: nqtot, tname, nqperes, nqfils
use read_profile_mod,         only: read_profile
use profile_temp_mod,         only: profile_temp
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, tsoil, &
                                    adsorption_soil, qsoil, igcm_h2o_ice_soil, porosity_reg,                   &
                                    ini_comsoil_h_slope_var, end_comsoil_h_slope_var, ads_massive_ice
use comvert_mod,              only: ap, bp, aps, bps, pa, preff, presnivs, pseudoalt, scaleheight
use dimradmars_mod,           only: tauvis, totcloudfrac, albedo, ini_dimradmars_mod_slope_var, end_dimradmars_mod_slope_var
use regular_lonlat_mod,       only: init_regular_lonlat
use mod_interface_dyn_phys,   only: init_interface_dyn_phys
use geometry_mod,             only: init_geometry, init_geometry_cell_area_for_outputs
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
use turb_mod,                 only: q2
use nonoro_gwd_ran_mod,       only: du_nonoro_gwd, dv_nonoro_gwd
use conf_phys_mod,            only: conf_phys
use paleoclimate_mod,         only: paleoclimate, ini_paleoclimate_h, end_paleoclimate_h, h2o_ice_depth
use comslope_mod,             only: nslope, subslope_dist, ini_comslope_h, end_comslope_h
use co2condens_mod,           only: CO2cond_ps
use callkeys_mod,             only: water, photochem, callthermos
use call_dayperi_mod,         only: call_dayperi
use initracer_mod,            only: initracer
! Mostly for XIOS outputs:
use mod_const_mpi,            only: COMM_LMDZ

implicit none

include "dimensions.h"

!=======================================================================
! Arguments
!=======================================================================
integer,      intent(in) :: ngrid, nlayer
real,         intent(in) :: odpref                     ! DOD reference pressure (Pa)
character(*), intent(in) :: start1Dname, startfiname   ! Name of starting files for 1D
logical,      intent(in) :: therestart1D, therestartfi ! Use of starting files for 1D

integer, intent(inout) :: nq

real, dimension(:,:,:), allocatable, intent(out) :: q          ! tracer mixing ratio (e.g. kg/kg)
real,                                intent(out) :: time       ! time (0<time<1; time=0.5 at noon)
real,                                intent(out) :: psurf      ! Surface pressure
real, dimension(nlayer),             intent(out) :: u, v       ! zonal, meridional wind
real, dimension(nlayer),             intent(out) :: temp       ! temperature at the middle of the layers
integer,                             intent(out) :: ndt
real,                                intent(out) :: ptif, pks
real,                                intent(out) :: dttestphys ! testphys1d timestep
real, dimension(:),     allocatable, intent(out) :: zqsat
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) :: gru, grv   ! prescribed "geostrophic" background wind
real, dimension(nlayer),             intent(out) :: w          ! "Dummy wind" in 1D
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(1),                  intent(out) :: latitude, longitude, cell_area
logical,                             intent(out) :: ctrl_h2ovap, ctrl_h2oice           ! Flags to prescribe atmospheric water vapour/ice profiles
real,                                intent(out) :: relaxtime_h2ovap, relaxtime_h2oice ! Relaxation time towards the atmospheric water vapour/ice profiles
real, dimension(nlayer),             intent(out) :: qref_h2ovap, qref_h2oice           ! User-defined atmospheric water vapour/ice profiles

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

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

! Use of starting files for 1D
logical                :: found
character(30)          :: header
real, dimension(100)   :: tab_cntrl

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

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

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

! LL: Subsurface water ice 
real :: rho_H2O_ice           ! Density (kg/m^3) of water ice

!=======================================================================
! Code
!=======================================================================
!------------------------------------------------------
! Prescribed constants to be set here
!------------------------------------------------------
! Mars planetary constants
! ------------------------
rad = 3397200.            ! mars radius (m) ~3397200 m
daysec = 88775.           ! length of a sol (s) ~88775 s
omeg = 4.*asin(1.)/daysec ! rotation rate (rad.s-1)
g = 3.72                  ! gravity (m.s-2) ~3.72
mugaz = 43.49             ! atmosphere molar mass (g.mol-1) ~43.49
rcp = .256793             ! = r/cp ~0.256793
r = 8.314511*1000./mugaz
cpp = r/rcp
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)
Lsperi = 251.             ! Solar longitude of perihelion

! 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)

! Water ice properties
! ---------------------------------------
inertieice = 2100.
rho_H2O_ice = 920. 

! 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 ...'
    error stop
else
    write(*,*) "init_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(*,*) "init_testphys1d: error reading number of tracers"
        write(*,*) "   (first line of traceur.def) "
        error stop
    endif
    if (nq < 1) then
        write(*,*) "init_testphys1d: error number of tracers"
        write(*,*) "is nq=",nq," but must be >=1!"
        error stop
    endif
endif

! Allocate arrays
allocate(tname(nq),q(1,nlayer,nq),zqsat(nlayer))
allocate(dq(1,nlayer,nq),dqdyn(1,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
        error stop 'init_testphys1d: error reading tracer names...'
    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

#ifdef CPP_XIOS
    call init_physics_distribution(regular_lonlat,4,1,1,1,nlayer,COMM_LMDZ)
#else
    call init_physics_distribution(regular_lonlat,4,1,1,1,nlayer,1)
#endif

! Date and local time at beginning of run
! ---------------------------------------
if (.not. therestartfi) 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(trim(startfiname))
    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 !(.not. therestartfi)

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

steps_per_sol = 48 ! Default value for day_step (old name for steps_per_sol)
write(*,*)'Number of time steps per sol?'
call getin("day_step",steps_per_sol)
write(*,*) " steps_per_sol (aka day_step) = ",steps_per_sol

outputs_per_sol = steps_per_sol ! Default value for outputs_per_sol

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

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

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

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

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

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

! Coefficient to control the surface pressure change
! To be defined in "callphys.def" so that the PEM can read it asa well
! If CO2cond_ps = 0, surface pressure is kept constant; If CO2cond_ps = 1, surface pressure varies normally
! CO2cond_ps = 300*400/144.37e6 = 8.31e-4 (ratio of polar cap surface over planetary surface) is a typical value for tests
! --------------------------------------------------------------------
CO2cond_ps = 1. ! Default value
write(*,*) 'Coefficient to control the surface pressure change?'
call getin("CO2cond_ps",CO2cond_ps)
if (CO2cond_ps < 0. .or. CO2cond_ps > 1.) then
    write(*,*) 'Value for ''CO2cond_ps'' is not between 0 and 1.'
    error stop 'Please, specify a correct value!'
endif

! 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. therestartfi) then
    paleomars = .false. ! Default: no water ice reservoir
    call getin("paleomars",paleomars)
    write(*,*) "paleomars=", paleomars
    if (paleomars) then
        write(*,*) "Orbital parameters taken from callphys.def"

        write(*,*)'Obliquity?'
        call getin("obliquit",obliquit)
        write(*,*) " obliquit =",obliquit

        write(*,*) 'Eccentricity (0<e<1)?'
        call getin('eccentric',eccentric)
        write(*,*)"eccentric =",eccentric
        if (eccentric < 0. .or. eccentric > 1.) then
            write(*,*) 'Value for ''eccentric'' is not between 0 and 1.'
            error stop 'Please, specify a correct value!'
        endif

        call getin('Lsperi',Lsperi)
        write(*,*) 'Solar longitude of perihelion (0<Ls<360)?'
        write(*,*)"Lsperi =",Lsperi
        if (eccentric < 0. .or. eccentric > 360.) then
            write(*,*) 'Value for ''Lsperi'' is not between 0 and 360.'
            error stop 'Please, specify a correct value!'
        endif
        
        Lsperi = Lsperi*pi/180. ! Convert to 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 PCM"
        write(*,*) "   periheli =",periheli
        write(*,*) "   aphelie  =",aphelie
        write(*,*) "   peri_day =",peri_day
    else
        write(*,*) "Default present-day orbital parameters"
     endif
else
    write(*,*) "Orbital parameters taken from startfi.nc"
endif !(.not. therestartfi)

! 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,[1])
call init_geometry_cell_area_for_outputs(1,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)
call initracer(ngrid,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(*,*) "init_testphys1d: setting iphysiq=1"
    iphysiq = 1
endif

! Initialize local slope parameters
! ---------------------------------
nslope_tmp = nslope
write(*,*)'Number of slopes?' ! Number of subgrid scale slopes
call getin("nslope",nslope_tmp)
write(*,*) " nslope = ", nslope_tmp
if (.not. therestartfi) then
    nslope = nslope_tmp
    call end_comslope_h
    call ini_comslope_h(1,nslope)
    select case (nslope)
        case(1)
            ! Sub-slopes parameters (minimum/maximun angles)
            def_slope(1) = -90.
            def_slope(2) = 90.
            ! Fraction of subslopes in mesh
            subslope_dist(1,1) = 1.
        case(3)
            ! Sub-slopes parameters (minimum/maximun angles)
            def_slope(1) = -50.
            def_slope(2) = -3.
            def_slope(3) = 3.
            def_slope(4) = 50.
            ! Fraction of subslopes in mesh
            subslope_dist(1,1) = 0.1
            subslope_dist(1,2) = 0.8
            subslope_dist(1,3) = 0.1
        case(5)
            ! Sub-slopes parameters (minimum/maximun angles)
            def_slope(1) = -43.
            def_slope(2) = -9.
            def_slope(3) = -3.
            def_slope(4) = 3.
            def_slope(5) = 9.
            def_slope(6) = 43.
            ! Fraction of subslopes in mesh
            subslope_dist(1,1) = 0.025
            subslope_dist(1,2) = 0.075
            subslope_dist(1,3) = 0.8
            subslope_dist(1,4) = 0.075
            subslope_dist(1,5) = 0.025
        case(7)
            ! Sub-slopes parameters (minimum/maximun angles)
            def_slope(1) = -43.
            def_slope(2) = -19.
            def_slope(3) = -9.
            def_slope(4) = -3.
            def_slope(5) = 3.
            def_slope(6) = 9.
            def_slope(7) = 19.
            def_slope(8) = 43.
            ! Fraction of subslopes in mesh
            subslope_dist(1,1) = 0.01
            subslope_dist(1,2) = 0.02
            subslope_dist(1,3) = 0.06
            subslope_dist(1,4) = 0.8
            subslope_dist(1,5) = 0.06
            subslope_dist(1,6) = 0.02
            subslope_dist(1,7) = 0.01
        case default
            error stop 'Number of slopes not possible: nslope should be 1, 3, 5 or 7!'
    end select

    call end_surfdat_h_slope_var
    call ini_surfdat_h_slope_var(1,nq,nslope)
    call end_comsoil_h_slope_var
    call ini_comsoil_h_slope_var(1,nslope)
    call end_dimradmars_mod_slope_var
    call ini_dimradmars_mod_slope_var(1,nslope)
    call end_paleoclimate_h
    call ini_paleoclimate_h(1,nslope)
else
    if (nslope /= nslope_tmp) then
        error stop 'Number of slopes not compatible between "startfi.nc" and "run.def".'
    else
        write(*,*) 'Be careful: subslopes distribution and definition are taken from "startfi.nc" and not initialized by default!'
    endif
endif

! Slope inclination angle (deg) 0: horizontal, 90: vertical
theta_sl = 0. ! default: no inclination
call getin("slope_inclination",theta_sl)

! Slope orientation (deg)
! 0 == Northward, 90 == Eastward, 180 == Southward, 270 == Westward
psi_sl = 0. ! default value
call getin("slope_orientation",psi_sl)

! Initialize tracers (surface and atmosphere)
!--------------------------------------------
write(*,*) "init_testphys1d: initializing tracers"
if (.not. therestart1D) then
    call read_profile(nq,nlayer,qsurf(1,:,1),q,play)
    do iq = 1,nq
        qsurf(1,iq,:) = qsurf(1,iq,1)
    enddo
else
    do iq = 1,nq
        read(3,*,iostat = ierr) header, (qsurf(1,iq,j), j = 1,size(qsurf,3)), (q(1,ilayer,iq), ilayer = 1,nlayer)
        if (ierr /= 0) error stop 'Not enough atmospheric layers defined in the file "'//trim(start1Dname)//'" for the tracer "'//trim(header)//'"!'
        if (trim(tname(iq)) /= trim(header)) then
            write(*,*) 'Tracer names between "traceur.def" and "'//trim(start1Dname)//'" do not match!'
            write(*,*) 'Please, write the tracer names in the same order for both files.'
            error stop
        endif
    enddo
endif

! Initialize albedo / soil thermal inertia
! ----------------------------------------
if (.not. therestartfi) 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,:,:) = 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)
    write(*,*) " subsurface_ice_depth = ",ice_depth
    h2o_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. therestartfi)

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

! For the non-orographic gravity waves scheme
du_nonoro_gwd(1,:) = 0
dv_nonoro_gwd(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 turbulent kinetic energy
!------------------------------------
q2 = 0.

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

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

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

if (.not. therestart1D) then
    ! Altitude in km in "profile_temp": divide zlay by 1000
    zkm(0) = 0.
    zkm(1:) = zlay/1000.

    call profile_temp(nlayer + 1,zkm,teta)
    tsurf = teta(0)
    temp = teta(1:)
else
    read(3,*) header, (tsurf(1,j), j = 1,size(tsurf,2)), (temp(ilayer), ilayer = 1,nlayer)
    close(3)
endif

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

if (.not. therestartfi) qsoil = 0.

if (.not. therestartfi) then
    ! Initialize depths
    ! -----------------
    do isoil = 0,nsoil - 1
       mlayer(isoil) = 2.e-4*(1 + isoil**2.9*(1 - exp(-real(isoil)/20.))) ! mid layer depth
    enddo
    do isoil = 1,nsoil - 1
       layer(isoil) = (mlayer(isoil) + mlayer(isoil - 1))/2
    enddo
    layer(nsoil) = 2*mlayer(nsoil - 1) - mlayer(nsoil - 2)

    ! 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
            if(adsorption_soil) then
            ! add subsurface ice in qsoil if one runs with adsorption
               if(ads_massive_ice) then
                   qsoil(1,1,igcm_h2o_ice_soil,:) = (ice_depth-layer(2))/(layer(1) - layer(2))*rho_H2O_ice ! assume no porosity in the ice
                   qsoil(1,2:,igcm_h2o_ice_soil,:) = rho_H2O_ice
               else
                   qsoil(1,1,igcm_h2o_ice_soil,:) = (ice_depth-layer(2))/(layer(1) - layer(2))*porosity_reg*rho_H2O_ice
                   qsoil(1,2:,igcm_h2o_ice_soil,:) = porosity_reg*rho_H2O_ice
               endif            
            endif
        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
            
            if(adsorption_soil) then
            ! add subsurface ice in qsoil if one runs with adsorption
                qsoil(1,:iref - 1,igcm_h2o_ice_soil,:) = 0.
                if(ads_massive_ice) then
                   qsoil(1,iref+1:,igcm_h2o_ice_soil,:) = rho_H2O_ice
                   qsoil(1,iref,igcm_h2o_ice_soil,:) = (ice_depth-layer(iref))/(layer(iref-1) - layer(iref))*rho_H2O_ice
                else
                   qsoil(1,iref+1:,igcm_h2o_ice_soil,:) = porosity_reg*rho_H2O_ice
                   qsoil(1,iref,igcm_h2o_ice_soil,:) = (ice_depth-layer(iref))/(layer(iref-1) - layer(iref))*porosity_reg*rho_H2O_ice
                endif
            endif
        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
    
    do isoil = 1,nsoil
        inertiesoil(1,isoil,:) = inertiedat(1,isoil)
        tsoil(1,isoil,:) = tsurf(1,:) ! soil temperature
    enddo
endif !(.not. therestartfi)

! Initialize subsurface geothermal flux
! -------------------------------------
flux_geo = 0.
call getin("flux_geo",flux_geo(1,1))
write(*,*) " flux_geo = ",flux_geo(1,1)
flux_geo = flux_geo(1,1)

! 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 = found)
    if (found) 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,1:llm,1:nq)
    psdyn(1:2,1) = psurf
    call inichim_newstart(ngrid,nq,qdyn,qsurf(1,:,1),psdyn,flagh2o,flagthermo)
    q(1,1:llm,1:nq) = qdyn(1,1,1:llm,1:nq)
endif

! Check if the surface is a water ice reservoir
! ---------------------------------------------
if (.not. therestartfi) 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

! Initialize perennial_co2ice
! ---------------------------
if (.not. therestartfi .and. paleoclimate) perennial_co2ice = 0.

! Prescription of atmospheric water vapour/ice profiles
! ----------------------------------------------
ctrl_h2ovap = .false. ; ctrl_h2oice = .false. ! Default: free atmospheric water vapour/ice profiles
relaxtime_h2ovap = -1. ; relaxtime_h2oice = -1. ! Default: no time relaxation
qref_h2ovap = 0. ; qref_h2oice = 0. ! Default
if (water) then
    write(*,*)'Force atmospheric water vapour profile (mixing ratio in kg/kg)?'
    call getin('ctrl_h2ovap',ctrl_h2ovap)
    write(*,*) 'ctrl_h2ovap = ', ctrl_h2ovap
    if (ctrl_h2ovap) then
        write(*,*) 'Atmospheric water vapour profile is prescribed'
        open(10,file = 'profile_ref_h2o_vap',status = 'old',form = 'formatted',action = 'read',iostat = ierr)
        if (ierr == 0) then
            do ilayer = 1,nlayer
                read(10,*,iostat=ierr) qref_h2ovap(ilayer)
                if (ierr /= 0) error stop 'Not enough atmospheric layers defined in the file "profile_ref_h2o_vap"!'
            enddo
        else
            error stop 'File "profile_ref_h2o_vap" was not found!'
        endif 
        close(10)

        write(*,*) 'Relax atmospheric water vapour profile with time constant (s)?'
        
        call getin('relaxtime_h2ovap',relaxtime_h2ovap)
        write(*,*) 'relaxtime_h2ovap = ', relaxtime_h2ovap
        if (relaxtime_h2ovap < 0.) then
            write(*,*) 'Atmospheric water vapour profile is not relaxed (forcing).'
        else
            write(*,*) 'Atmospheric water vapour profile is relaxed towards the profile in "profile_ref_h2o_vap"'
        endif
    else
        write(*,*) 'Free atmospheric water vapour profile'
    endif

    write(*,*)'Force atmospheric water ice profile (mixing ratio in kg/kg)?'
    call getin('ctrl_h2oice',ctrl_h2oice)
    write(*,*) 'ctrl_h2oice = ', ctrl_h2oice
    if (ctrl_h2oice) then
        write(*,*) 'Atmospheric water ice profile is prescribed'
        open(10,file = 'profile_ref_h2o_ice',status = 'old',form = 'formatted',action = 'read',iostat = ierr)
        if (ierr == 0) then
            do ilayer = 1,nlayer
                read(10,*,iostat=ierr) qref_h2oice(ilayer)
                if (ierr /= 0) error stop 'Not enough atmospheric layers defined in the file "profile_ref_h2o_ice"!'
            enddo
        else
            error stop 'File "profile_ref_h2o_ice" was not found!'
        endif 
        close(10)

        write(*,*) 'Relax atmospheric water ice profile with time constant (s)?'
        
        call getin('relaxtime_h2oice',relaxtime_h2oice)
        write(*,*) 'relaxtime_h2oice = ', relaxtime_h2oice
        if (relaxtime_h2oice < 0.) then
            write(*,*) 'Atmospheric water ice profile is not relaxed (forcing).'
        else
            write(*,*) 'Atmospheric water ice profile is relaxed towards the profile in "profile_ref_h2o_ice"'
        endif
    else
        write(*,*) 'Free atmospheric water ice profile'
    endif
endif

END SUBROUTINE init_testphys1d

END MODULE init_testphys1d_mod