module init_pem1d_mod

implicit none

contains

#ifdef CPP_1D

subroutine init_pem1d(llm_in,nqtot_in,u,v,temp,q,psurf,time,ap_out,bp_out)

!-------------- init_pem1d -------------
!
! This function is a copy of the test_phys_1d.F90 initialisation part.
! It is meant to initialize all the variables in modules and for the ouput
! that can't be initialzed via the startfi file.
! 
! Basically it reads run.def and save variables values 
!
!--------------              -------------

      USE ioipsl_getincom, only: getin
      use comcstfi_h, only: pi, rad, omeg, g, mugaz, rcp, r, cpp
      use time_phylmdz_mod, only: daysec, dtphys, 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, emissiv, emisice,&
                          albedice, iceradius, dtemisice, z0,&
                          zmea, zstd, zsig, zgam, zthe, phisfi,&
                          watercaptag, watercap, hmons, summit, base,&
                          tsurf, emis,qsurf        
      use infotrac, only: nqtot, tname, nqperes,nqfils
      use regular_lonlat_mod, only: init_regular_lonlat
      use mod_grid_phy_lmdz, only : regular_lonlat
      USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff,sig,&
                            presnivs,pseudoalt,scaleheight
      use dimradmars_mod, only: tauvis,totcloudfrac
      use mod_interface_dyn_phys, only: init_interface_dyn_phys
      use geometry_mod, only: init_geometry
      use dimphy, only : init_dimphy
      USE phys_state_var_init_mod, ONLY: phys_state_var_init
      use comgeomfi_h, only: sinlat, ini_fillgeom
      use slope_mod, only: theta_sl, psi_sl
      use comslope_mod, only: def_slope,subslope_dist,def_slope_mean
      use dimradmars_mod, only: tauvis,totcloudfrac
      use dust_param_mod, only: tauscaling
      use tracer_mod, only: igcm_h2o_vap,igcm_h2o_ice,igcm_co2,noms
      USE logic_mod, ONLY: hybrid
      USE vertical_layers_mod, ONLY: init_vertical_layers
      use comsoil_h, only: volcapa, layer, mlayer, inertiedat, &
                          inertiesoil,nsoilmx,flux_geo
      USE read_profile_mod, ONLY: read_profile
      use inichim_newstart_mod, only: inichim_newstart
      use physics_distribution_mod, only: init_physics_distribution
      use iostart, only: open_startphy,get_var, close_startphy
#ifdef CPP_XIOS
      use mod_const_mpi, only: COMM_LMDZ
#endif
      include "dimensions.h"
      integer, parameter :: ngrid = 1 !(2+(jjm-1)*iim - 1/jjm)
      integer, parameter :: nlayer = llm
      include "callkeys.h"
      include "netcdf.inc"

!-------------- INPUT VARIABLES -------------

      integer, intent(in) :: llm_in,nqtot_in

!-------------- OUTPUT VARIABLES -------------

      real, intent(out) :: u(nlayer), v(nlayer) ! zonal, meridional wind
      real, intent(out) :: temp(nlayer)         ! temperature at the middle of the layers
      real, intent(out) :: q(2,llm_in,nqtot_in) ! tracer mixing ratio (e.g. kg/kg)
      real, intent(out) :: psurf(2)
      real, intent(out) :: time                 ! time (0<time<1 ; time=0.5 a midi)
      real, intent(out) :: ap_out(llm_in + 1), bp_out(llm_in + 1)


!-------------- LOCAL VARIABLES -------------

      REAL halfaxe, excentric, Lsperi
      integer :: nq=1 ! number of tracers
      real :: latitude(1), longitude(1), cell_area(1)

!   MVals: isotopes as in the dynamics (CRisi)
      INTEGER :: ifils,ipere,generation
      CHARACTER(len=30), ALLOCATABLE, DIMENSION(:) :: tnom_transp ! transporting fluid short name
      CHARACTER(len=80) :: line ! to store a line of text     
      INTEGER ierr0
      LOGICAL :: continu

      logical :: present,found
      INTEGER nlevel,nsoil,ndt
      REAL gru,grv   ! prescribed "geostrophic" background wind
      REAL pks, ptif, w(nlayer)
      REAL q2(nlayer+1)   ! Turbulent Kinetic Energy
      REAL play(nlayer)   ! Pressure at the middle of the layers (Pa)
      REAL plev(nlayer+1) ! intermediate pressure levels (pa)
      REAL zlay(nlayer)   ! altitude estimee dans les couches (km)
      REAL tmp1(0:nlayer),tmp2(0:nlayer)
      INTEGER flagthermo,flagh2o
      real atm_wat_profile, atm_wat_tau
      real, dimension(:), allocatable :: zqsat ! useful for (atm_wat_profile=2)

      INTEGER :: ierr,iq,ilayer,ilevel,isoil
      INTEGER day0,dayn          ! date initial (sol ; =0 a Ls=0) and final
      REAL day           ! date durant le run
      REAL tab_cntrl(100)
      real,parameter :: odpref=610. ! DOD reference pressure (Pa)

!   LL: Subsurface geothermal flux
      real :: flux_geo_tmp

!   RV & JBC: Use of "start1D.txt" and "startfi.nc" files
      logical                :: therestart1D, therestartfi
      character(len = 30)    :: header

      real,allocatable :: qdyn(:,:,:,:),psdyn(:,:)


! Checking if the file "start1D.txt" exists
inquire(file ='start1D.txt',exist = therestart1D) 
if (.not. therestart1D) then
    write(*,*) 'There is no "start1D.txt" file!'
    write(*,*) 'It is required to initialize the 1D PEM.'
    stop
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

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

!     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.314511E+0 *1000.E+0/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.E+0

! 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(*,*) "pem1d: Reading file traceur.def"
        ! read number of tracers:
        read(90,*,iostat = ierr) nq
        write(*,*) "nq",nq
        nqtot = nq ! set value of nqtot (in infotrac module) as nq
        if (ierr /= 0) then
            write(*,*) "pem1d: error reading number of tracers"
            write(*,*) "   (first line of traceur.def) "
            stop
        endif
        if (nq < 1) then
            write(*,*) "pem1d: error number of tracers"
            write(*,*) "is nq=",nq," but must be >=1!"
            stop
        endif
    endif
    ! allocate arrays:
    allocate(tname(nq))
    allocate(qsurf(1,1,nq))
    allocate(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(*,*) 'init_pem1d: 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) .eq. 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.gt.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.gt.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(*,*) "init_pem1d: initializing tracers"
        do iq = 1,nq
            open(3,file = 'start1D.txt',status = "old",action = "read")
            read(3,*) header, qsurf(1,1,iq), (q(1,ilayer,iq), ilayer = 1,nlayer)
            if (tname(iq) /= trim(header)) then
                write(*,*) 'Tracer names not compatible for initialization with "start1D.txt"!'
                stop
            endif
        enddo
        q(2,:,:) = q(1,:,:)

#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

!  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     
      dtphys=daysec/day_step

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

! 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

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

!  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 dtphys 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,dtphys,&
               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(*,*) "init_pem1d: setting iphysiq=1"
        iphysiq=1
      endif

! 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.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.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.E+0
      zstd(1)=0.E+0
      zsig(1)=0.E+0
      zgam(1)=0.E+0
      zthe(1)=0.E+0
!
!  for the slope wind scheme
!  ---------------------------------
!  
      hmons(1)=0.E+0
      write(*,*)'hmons is initialized to ',hmons(1)
      summit(1)=0.E+0
      write(*,*)'summit is initialized to ',summit(1)
      base(1)=0.E+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.E+0

!   Initialization to take into account prescribed winds
!   ------------------------------------------------------
      ptif=2.E+0*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
      read(3,*) header, (u(ilayer), ilayer = 1,nlayer)
      read(3,*) header, (v(ilayer), ilayer = 1,nlayer)
      w = 0. ! default: no vertical wind

!     Initialize turbulente kinetic energy
      DO ilevel=1,nlevel
         q2(ilevel)=0.E+0
      ENDDO

!  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") then
          igcm_co2=iq
        endif
      enddo
      if (igcm_co2==0) then
        write(*,*) "init_pem1d error, missing co2 tracer!"
        stop
      endif

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

      DO ilevel=1,nlevel
        plev(ilevel)=ap(ilevel)+psurf(1)*bp(ilevel)
      ENDDO

      DO ilayer=1,nlayer
        play(ilayer)=aps(ilayer)+psurf(1)*bps(ilayer)
      ENDDO

      DO ilayer=1,nlayer
        zlay(ilayer)=-200.E+0 *r*log(play(ilayer)/plev(1))/g
      ENDDO


!  Initialize temperature profile
!  --------------------------------------
      pks=psurf(1)**rcp

! altitude in km in profile: divide zlay by 1000
      tmp1(0)=0.E+0
      DO ilayer=1,nlayer
        tmp1(ilayer)=zlay(ilayer)/1000.E+0
      ENDDO

      call profile(nlayer+1,tmp1,tmp2)

      read(3,*) header, tsurf, (temp(ilayer), ilayer = 1,nlayer)
      close(3)

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

      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
      enddo
      do isoil=1,nsoil
        layer(isoil)=2.e-4*(2.**(isoil-1)) ! 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=present)
         if (present) 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(1)
         call inichim_newstart(ngrid,nq,qdyn,qsurf,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 
! --------------------------------------------------
      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

      ap_out = ap
      bp_out = bp

end subroutine init_pem1d

#endif

end module init_pem1d_mod