module init_phys_1d_mod

implicit none

contains

#ifdef CPP_1D

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

!-------------- init_phys_1d -------------
!
! 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
      REAL,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

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


! check if 'run.def' file is around (otherwise reading parameters
! from callphys.def via getin() routine won't work.
      open(99,file='run.def',status='old',form='formatted',&
          iostat=ierr)
      if (ierr.ne.0) then
        write(*,*) 'Cannot find required file "run.def"'
        write(*,*) '  (which should contain some input parameters'
        write(*,*) '   along with the following line:'
        write(*,*) '   INCLUDEDEF=callphys.def'
        write(*,*) '   )'
        write(*,*) ' ... might as well stop here ...'
        stop
      else
        close(99)
      endif

! ------------------------------------------------------
!  Prescribed constants to be set here
! ------------------------------------------------------

!      if(.not. startfile_1D ) then

      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)

!      endif ! .not. startfile_1D

!     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.ne.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
          nqtot=nqtot_in ! set value of nqtot (in infotrac module) as nq
          nq=nqtot_in
        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.ne.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.ne.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(*,*) "testphys1d: initializing tracers"
        call read_profile(nq, nlayer, qsurf, q(1,:,:))
        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

      call open_startphy("startfi_evol.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


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

      day_step=48 ! default value for day_step
      PRINT *,'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
      PRINT *,'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
      PRINT *,'Surface pressure (Pa) ?'
      call getin("psurf",psurf(1))
      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

! Orbital parameters
! ------------------

!  latitude/longitude
!  ------------------
      latitude(1)=0 ! default value for latitude
      PRINT *,'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(*,*) "testphys1d: setting iphysiq=1"
        iphysiq=1
      endif 

!  Initialize albedo / soil thermal inertia
!  ----------------------------------------
!


! 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
      PRINT *,'hmons is initialized to ',hmons(1)
      summit(1)=0.E+0
      PRINT *,'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
      PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
      call getin("u",gru)
      write(*,*) " u = ",gru
      grv=0. !default value for grv
      PRINT *,'meridional northward component of the geostrophic',&
     ' wind (m/s) ?'
      call getin("v",grv)
      write(*,*) " v = ",grv

!     Initialize winds  for first time step
      DO ilayer=1,nlayer
         u(ilayer)=gru
         v(ilayer)=grv
         w(ilayer)=0 ! default: no vertical wind
      ENDDO

!     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(*,*) "testphys1d error, missing co2 tracer!"
        stop
      endif

!  Compute pressures and altitudes of atmospheric levels 
!  ----------------------------------------------------------------

!    Vertical Coordinates
!    """"""""""""""""""""
      hybrid=.true.
      PRINT *,'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)

      tsurf=tmp2(0)
      DO ilayer=1,nlayer
        temp(ilayer)=tmp2(ilayer)
      ENDDO
      

! 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
      print *,'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
      print *,'Force atmospheric water vapor profile ?'
      call getin("atm_wat_profile",atm_wat_profile)
      write(*,*) "atm_wat_profile = ", atm_wat_profile
      if (atm_wat_profile.EQ.-1) then
        write(*,*) "Free atmospheric water vapor profile"
        write(*,*) "Total water is conserved on the column"
      else if (atm_wat_profile.EQ.0) then
        write(*,*) "Dry atmospheric water vapor profile"
      else if ((atm_wat_profile.GT.0).and.(atm_wat_profile.LE.1)) then
        write(*,*) "Prescribed atmospheric water vapor MMR profile"
        write(*,*) "Unless it reaches saturation (maximal value)"
        write(*,*) "MMR chosen : ", atm_wat_profile
      endif 

      ap_out=ap
      bp_out=bp

end subroutine init_phys_1d

#endif

end module init_phys_1d_mod