source: trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90 @ 3203

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,atm_wat_profile,atm_wat_tau,CO2cond_ps)

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, &
                                    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 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, qsoil, &
                                    ini_comsoil_h_slope_var, end_comsoil_h_slope_var
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
use comslope_mod,             only: nslope, subslope_dist, ini_comslope_h, end_comslope_h
! 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)
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                        ! useful for (atm_wat_profile=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) :: 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
real,                                intent(out) :: atm_wat_profile, atm_wat_tau ! Force atmospheric water profiles
real,                                intent(out) :: CO2cond_ps                   ! Coefficient to control the surface pressure change

!=======================================================================
! Local variables
!=======================================================================
integer                   :: ierr, iq, j, 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 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 = 2100. ! ice thermal inertia
integer :: iref

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

! 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 ...'
    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

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
    open(3,file = trim(start1Dname),status = "old",action = "read")
    read(3,*) header, psurf
endif
write(*,*) " psurf = ",psurf

! Coefficient to control the surface pressure change
! To be defined in "callphys.def" so that the PEM can read it asa well
! --------------------------------------------------------------------
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

! 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. therestartfi) 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. 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)
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,qsurf)

! 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 = 1 ! default: one slope
if (.not. therestartfi) then
    ! Number of subgrid scale slopes
    write(*,*)'Number of slopes?'
    call getin("nslope",nslope)
    write(*,*) " nslope = ", nslope
    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)

    ! Fraction of subslopes in mesh
    !write(*,*)'What subslope distribution do you want to use?'
    !call getin("slope_distribution",subslope_dist)
    !write(*,*) " subslope_dist = ", subslope_dist(1,:)
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)
    do iq = 1,nq
        qsurf(1,iq,:) = qsurf(1,iq,1)
    enddo
else
    do iq = 1,nq
        read(3,*) header, (qsurf(1,iq,j), j = 1,size(qsurf,3)), (q(1,ilayer,iq), ilayer = 1,nlayer)
        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

    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)

! 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(1,:), j = 1,size(tsurf,2)), (temp(ilayer), ilayer = 1,nlayer)
    close(3)
endif

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

if (.not. therestartfi) 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

    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 soil content
! -----------------------
if (.not. therestartfi) qsoil = 0.

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

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

! 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 if (atm_wat_profile .eq. 2) then
        write(*,*) 'Prescribed atmospheric water vapor profile'
        write(*,*) 'like present day northern midlatitude'
        write(*,*) 'Unless it reaches saturation (maximal value)'
    else
        error stop 'Water vapor profile value not correct!'
    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!'
            error stop 'Please, specify atm_wat_profile'
        endif
    endif
endif

END SUBROUTINE init_testphys1d

END MODULE init_testphys1d_mod