Index: trunk/LMDZ.MARS/changelog.txt
===================================================================
--- trunk/LMDZ.MARS/changelog.txt	(revision 3055)
+++ trunk/LMDZ.MARS/changelog.txt	(revision 3056)
@@ -4215,4 +4215,5 @@
 Upgrade of "testphys1d" to Fortran 90. Cleaning of the subroutine and minor optimizations of the code.
 Correction of a bug: 'inertiedat(1,1)' was overwritten by 'inertieice'.
+From now on, in "testphys1d.F90", initialization is entirely done by a new subroutine called 'init_testphys1d' ("init_testphys1d_mod.F90").
 
 == 27/09/2023 == EM
Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90	(revision 3056)
+++ trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90	(revision 3056)
@@ -0,0 +1,727 @@
+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
+
Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90	(revision 3055)
+++ trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90	(revision 3056)
@@ -1,41 +1,22 @@
 PROGRAM testphys1d
 
-use ioipsl_getincom,          only: getin ! To use 'getin'
-use dimphy,                   only: init_dimphy
-use mod_grid_phy_lmdz,        only: regular_lonlat
-use infotrac,                 only: nqtot, tname, nqperes, nqfils
-use comsoil_h,                only: volcapa, layer, mlayer, inertiedat, inertiesoil, nsoilmx, flux_geo
-use comgeomfi_h,              only: sinlat, ini_fillgeom
-use surfdat_h,                only: albedodat, z0_default, emissiv, emisice, albedice, iceradius,     &
-                                    dtemisice, z0, zmea, zstd, zsig, zgam, zthe, phisfi, watercaptag, &
-                                    watercap, hmons, summit, base, perenial_co2ice
-use slope_mod,                only: theta_sl, psi_sl
-use comslope_mod,             only: def_slope, subslope_dist, def_slope_mean
-use phyredem,                 only: physdem0, physdem1
-use geometry_mod,             only: init_geometry
-use watersat_mod,             only: watersat
-use tracer_mod,               only: igcm_h2o_vap, igcm_h2o_ice, igcm_co2,noms
-use planete_h,                only: year_day, periheli, aphelie, peri_day, obliquit, emin_turb, lmixmin
-use comcstfi_h,               only: pi, rad, omeg, g, mugaz, rcp, r, cpp
-use time_phylmdz_mod,         only: daysec, day_step, ecritphy, iphysiq
-use dimradmars_mod,           only: tauvis, totcloudfrac
-use dust_param_mod,           only: tauscaling
-use comvert_mod,              only: ap, bp, aps, bps, pa, preff, sig, presnivs, pseudoalt, scaleheight
-use vertical_layers_mod,      only: init_vertical_layers
-use logic_mod,                only: hybrid
-use physics_distribution_mod, only: init_physics_distribution
-use regular_lonlat_mod,       only: init_regular_lonlat
-use mod_interface_dyn_phys,   only: init_interface_dyn_phys
-use phys_state_var_init_mod,  only: phys_state_var_init
-use physiq_mod,               only: physiq
-use read_profile_mod,         only: read_profile
-use inichim_newstart_mod,     only: inichim_newstart
-use phyetat0_mod,             only: phyetat0
-use netcdf,                   only: NF90_OPEN, NF90_NOERR, NF90_NOWRITE, nf90_strerror
-use iostart,                  only: open_startphy, get_var, close_startphy
-use write_output_mod,         only: write_output
+use comsoil_h,               only: inertiedat, inertiesoil, nsoilmx
+use surfdat_h,               only: albedodat, perenial_co2ice, watercap
+use comslope_mod,            only: def_slope, subslope_dist
+use phyredem,                only: physdem0, physdem1
+use watersat_mod,            only: watersat
+use tracer_mod,              only: igcm_h2o_vap, igcm_h2o_ice, noms
+use comcstfi_h,              only: pi, rad, omeg, g, mugaz, rcp, r, cpp
+use time_phylmdz_mod,        only: daysec, day_step
+use dimradmars_mod,          only: tauvis, totcloudfrac
+use dust_param_mod,          only: tauscaling
+use comvert_mod,             only: ap, bp, aps, bps, pa, preff, sig
+use physiq_mod,              only: physiq
+use phyetat0_mod,            only: phyetat0
+use write_output_mod,        only: write_output
+use init_testphys1d_mod,     only: init_testphys1d
 ! Mostly for XIOS outputs:
-use mod_const_mpi,            only: init_const_mpi, COMM_LMDZ
-use parallel_lmdz,            only: init_parallel
+use mod_const_mpi,           only: init_const_mpi
+use parallel_lmdz,           only: init_parallel
 
 implicit none
@@ -56,6 +37,6 @@
 !
 !   update: 12/06/2003, including chemistry (S. Lebonnois)
-!                            and water ice (F. Montmessin)
-!           26/09/2023, conversion in F90 (JB Clément)
+!                       and water ice (F. Montmessin)
+!           27/09/2023, upgrade to F90 (JB Clément)
 !
 !=======================================================================
@@ -77,51 +58,40 @@
 ! Declarations
 !--------------------------------------------------------------
-integer, parameter                :: ngrid = 1 ! (2+(jjm-1)*iim - 1/jjm)
+integer, parameter                :: ngrid = 1     ! (2+(jjm-1)*iim - 1/jjm)
 integer, parameter                :: nlayer = llm
 real, parameter                   :: odpref = 610. ! DOD reference pressure (Pa)
-integer                           :: unitstart ! unite d'ecriture de "startfi"
-integer                           :: nlevel, nsoil, ndt
-integer                           :: ilayer, ilevel, isoil, idt, iq
+integer                           :: unitstart     ! unite d'ecriture de "startfi"
+integer                           :: ndt, ilayer, ilevel, isoil, idt, iq
 logical                           :: firstcall, lastcall
-integer                           :: day0, dayn ! initial (sol ; =0 at Ls=0) and final date
-real                              :: day        ! date during the run
-real                              :: time       ! time (0<time<1 ; time=0.5 a midi)
-real                              :: dttestphys ! testphys1d timestep
-real, dimension(nlayer)           :: play       ! Pressure at the middle of the layers (Pa)
-real, dimension(nlayer + 1)       :: plev       ! intermediate pressure levels (pa)
-real                              :: psurf      ! Surface pressure
-real, dimension(1)                :: tsurf      ! Surface temperature
-real, dimension(nlayer)           :: u, v       ! zonal, meridional wind
-real                              :: gru, grv   ! prescribed "geostrophic" background wind
-real, dimension(nlayer)           :: temp       ! temperature at the middle of the layers
-real, dimension(:,:), allocatable :: q          ! tracer mixing ratio (e.g. kg/kg)
-real, dimension(:),   allocatable :: qsurf      ! tracer surface budget (e.g. kg.m-2)
-real, dimension(nsoilmx)          :: tsoil      ! subsurface soik temperature (K)
-real, dimension(1)                :: emis       ! surface layer
-real, dimension(1,1)              :: albedo     ! surface albedo
-real, dimension(1)                :: wstar = 0. ! Thermals vertical velocity
-real, dimension(nlayer + 1)       :: q2         ! Turbulent Kinetic Energy
-real, dimension(nlayer)           :: zlay       ! altitude estimee dans les couches (km)
-
-! Physical and dynamical tandencies (e.g.  m.s-2, K/s, Pa/s)
+integer                           :: day0          ! initial (sol ; =0 at Ls=0)
+real                              :: day           ! date during the run
+real                              :: time          ! time (0<time<1 ; time=0.5 a midi)
+real                              :: dttestphys    ! testphys1d timestep
+real, dimension(nlayer)           :: play          ! Pressure at the middle of the layers (Pa)
+real, dimension(nlayer + 1)       :: plev          ! intermediate pressure levels (pa)
+real                              :: psurf         ! Surface pressure
+real, dimension(1)                :: tsurf         ! Surface temperature
+real, dimension(nlayer)           :: u, v          ! zonal, meridional wind
+real                              :: gru, grv      ! prescribed "geostrophic" background wind
+real, dimension(nlayer)           :: temp          ! temperature at the middle of the layers
+real, dimension(:,:), allocatable :: q             ! tracer mixing ratio (e.g. kg/kg)
+real, dimension(:),   allocatable :: qsurf         ! tracer surface budget (e.g. kg.m-2)
+real, dimension(nsoilmx)          :: tsoil         ! subsurface soik temperature (K)
+real, dimension(1)                :: emis          ! surface layer
+real, dimension(1,1)              :: albedo        ! surface albedo
+real, dimension(1)                :: wstar = 0.    ! Thermals vertical velocity
+real, dimension(nlayer + 1)       :: q2            ! Turbulent Kinetic Energy
+
+! Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s)
 real, dimension(nlayer)           :: du, dv, dtemp, dudyn, dvdyn, dtempdyn
 real, dimension(1)                :: dpsurf
-real, dimension(:,:), allocatable :: dq
-real, dimension(:,:), allocatable :: dqdyn
+real, dimension(:,:), allocatable :: dq, dqdyn
 
 ! Various intermediate variables
-integer                   :: flagthermo, flagh2o
-real                      :: zls
-real, dimension(nlayer)   :: phi, h, s, w
-real                      :: pks, ptif
-real                      :: qtotinit,qtot
-integer                   :: ierr, aslun
-real, dimension(0:nlayer) :: tmp1, tmp2
-integer                   :: nq = 1 ! number of tracers
-real, dimension(1)        :: latitude, longitude, cell_area
-
-! Dummy variables along "dynamics scalar grid"
-real, dimension(:,:,:,:), allocatable :: qdyn
-real, dimension(:,:),     allocatable :: psdyn
+integer                 :: aslun
+real                    :: zls, pks, ptif, qtotinit, qtot
+real, dimension(nlayer) :: phi, h, s, w
+integer                 :: nq = 1 ! number of tracers
+real, dimension(1)      :: latitude, longitude, cell_area
 
 character(len = 2)  :: str2
@@ -129,31 +99,10 @@
 character(len = 44) :: txt
 
-! New flag to compute paleo orbital configurations + few variables JN
-real    :: halfaxe, excentric, Lsperi
-logical :: paleomars
+! RV & JBC: Use of "start1D.txt" and "startfi.nc" files
+logical :: startfiles_1D
 
 ! JN & JBC: Force atmospheric water profiles
 real                            :: atm_wat_profile, atm_wat_tau
 real, dimension(:), allocatable :: zqsat ! useful for (atm_wat_profile=2)
-
-! 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: Subsurface geothermal flux
-real :: flux_geo_tmp
-
-! RV & JBC: Use of "start1D.txt" and "startfi.nc" files
-logical                :: startfiles_1D, found, therestart1D, therestartfi
-character(len = 30)    :: header
-real, dimension(100)   :: tab_cntrl
-real, dimension(1,2,1) :: albedo_read ! surface albedo
-
-! 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
 !=======================================================================
 
@@ -171,608 +120,7 @@
 !call initcomgeomphy
 
-!------------------------------------------------------
-! 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
-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 mola mass (g.mol-1) ~43.49
-rcp = .256793                    ! = r/cp ~0.256793
-r = 8.314511*1000./mugaz
-cpp = r/rcp
-year_day = 669                   ! lenght 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
-excentric = 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('excentric ',excentric)
-        write(*,*)"excentric =",excentric
-        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 - excentric)
-        aphelie = halfaxe*(1 + excentric)
-        call call_dayperi(Lsperi,excentric,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
+call 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)
 
 ! Write a "startfi" file
@@ -912,2 +260,3 @@
 !***********************************************************************
 !***********************************************************************
+