C======================================================================
      PROGRAM newstart
c=======================================================================
c
c
c   Auteur:   Christophe Hourdin/Francois Forget/Yann Wanherdrick
c   ------
c
c   Objet:  Create or modify the initial state for the LMD Titan GCM
c   -----           (fichiers NetCDF start et startfi)
c
c
c=======================================================================

      use mod_phys_lmdz_para, only: is_parallel, is_sequential,
     &                              is_mpi_root, is_omp_root,
     &                              is_master
      use infotrac, only: infotrac_init, nqtot, tname
      USE comchem_h, ONLY: nkim, nlaykim_up, ykim_up
      USE comchem_newstart_h
      USE comsoil_h, ONLY: nsoilmx, layer, mlayer, inertiedat
      USE surfdat_h, ONLY: phisfi, albedodat,
     &                     zmea, zstd, zsig, zgam, zthe
      use datafile_mod, only: datadir, surfdir
      use ioipsl_getin_p_mod, only: getin_p
      use control_mod, only: day_step, iphysiq, anneeref, planet_type
      use phyredem, only: physdem0, physdem1
      use iostart, only: open_startphy
      use filtreg_mod, only: inifilr
      USE mod_const_mpi, ONLY: COMM_LMDZ
      USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff
      USE comconst_mod, ONLY: lllm,daysec,dtvr,dtphys,cpp,kappa,
     .                        rad,omeg,g,r,pi
      USE serre_mod, ONLY: alphax
      USE temps_mod, ONLY: day_ini
      USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
      use tabfi_mod, only: tabfi
      use callkeys_mod, only: callchim
      use iniphysiq_mod, only: iniphysiq
      use phyetat0_mod, only: phyetat0
      use exner_hyb_m, only: exner_hyb
      use tracer_h
      implicit none

      include "dimensions.h"
      integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm) 
      include "paramet.h"
      include "comgeom2.h"
      include "comdissnew.h"
      include "netcdf.inc"

c=======================================================================
c   Declarations
c=======================================================================

c Variables dimension du fichier "start_archive"
c------------------------------------
      CHARACTER        relief*3


c Variables pour les lectures NetCDF des fichiers "start_archive" 
c--------------------------------------------------
      INTEGER nid_dyn, nid_fi,nid,nvarid
      INTEGER length
      parameter (length = 100)
      INTEGER tab0
      INTEGER   NB_ETATMAX
      parameter (NB_ETATMAX = 100)

      REAL  date
      REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec

c Variable histoire 
c------------------
      REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
      REAL phis(iip1,jjp1)
      REAL,ALLOCATABLE :: q(:,:,:,:)               ! champs advectes

c autre variables dynamique nouvelle grille
c------------------------------------------
      REAL pks(iip1,jjp1)
      REAL w(iip1,jjp1,llm+1)
      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
!      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
!      REAL dh(ip1jmp1,llm),dp(ip1jmp1)
      REAL phi(iip1,jjp1,llm)

      integer klatdat,klongdat
      PARAMETER (klatdat=180,klongdat=360)

c Physique sur grille scalaire 
c----------------------------
      real zmeaS(iip1,jjp1),zstdS(iip1,jjp1)
      real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1)

c variable physique
c------------------
      REAL tsurf(ngridmx)        ! surface temperature
      REAL,ALLOCATABLE :: tsoil(:,:) ! soil temperature
      REAL emis(ngridmx)        ! surface emissivity
      real emisread             ! added by RW
      REAL,ALLOCATABLE :: qsurf(:,:)
      REAL q2(ngridmx,llm+1)
!      REAL rnaturfi(ngridmx)
      real alb(iip1,jjp1),albfi(ngridmx) ! albedos
      real,ALLOCATABLE :: ith(:,:,:),ithfi(:,:) ! thermal inertia (3D)
      real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D)
      REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)

      INTEGER i,j,l,isoil,ig,idum
      real mugaz ! molar mass of the atmosphere

      integer ierr

c Variables on the new grid along scalar points 
c------------------------------------------------------
!      REAL p(iip1,jjp1)
      REAL t(iip1,jjp1,llm)
      REAL tset(iip1,jjp1,llm)
      real phisold_newgrid(iip1,jjp1)
      REAL :: teta(iip1, jjp1, llm)
      REAL :: pk(iip1,jjp1,llm)
      REAL :: pkf(iip1,jjp1,llm)
      REAL :: ps(iip1, jjp1)
      REAL :: masse(iip1,jjp1,llm)
      REAL :: xpn,xps,xppn(iim),xpps(iim)
      REAL :: p3d(iip1, jjp1, llm+1)
!      REAL dteta(ip1jmp1,llm)

c Variable de l'ancienne grille 
c------------------------------
      real time
      real tab_cntrl(100)
      real tab_cntrl_bis(100)
      
c variables diverses
c-------------------
      real choix_1,pp
      character*80      fichnom
      character*250  filestring
      integer Lmodif,iq
      character modif*20
      real z_reel(iip1,jjp1)
      real tsud,albsud,alb_bb,ith_bb,Tiso,Tabove
      real ptoto,pcap,patm,airetot,ptotn,patmn,psea
!      real ssum
      character*1 yes
      logical :: flagtset=.false. ,  flagps0=.false.
      real val, val2, val3, val4 ! to store temporary variables

      INTEGER :: itau
      
      character(len=20) :: txt ! to store some text
      character(len=50) :: surfacefile ! "surface.nc" (or equivalent file)
      character(len=150) :: longline
      integer :: count
      real :: profile(llm+1) ! to store an atmospheric profile + surface value

!     added by BC for equilibrium temperature startup
      real teque

!     added by RW for nuketharsis
      real fact1
      real fact2

!     added by JVO for Titan
      REAL tankCH4(ngridmx)      

c sortie visu pour les champs dynamiques
c---------------------------------------
!      INTEGER :: visuid
!      real :: time_step,t_ops,t_wrt
!      CHARACTER*80 :: visu_file

      cpp    = 0.
      preff  = 0.
      pa     = 0. ! to ensure disaster rather than minor error if we don`t
                  ! make deliberate choice of these values elsewhere.

      planet_type="titan"

! initialize "serial/parallel" related stuff
! (required because we call tabfi() below, before calling iniphysiq)
      is_sequential=.true.
      is_parallel=.false.
      is_mpi_root=.true.
      is_omp_root=.true.
      is_master=.true.
      
! Load tracer number and names:
      call infotrac_init
! allocate arrays
      allocate(q(iip1,jjp1,llm,nqtot))
      allocate(qsurf(ngridmx,nqtot))
      
! get value of nsoilmx and allocate corresponding arrays
      ! a default value of nsoilmx is set in comsoil_h
      call getin_p("nsoilmx",nsoilmx)
      
      allocate(tsoil(ngridmx,nsoilmx))
      allocate(ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx))
      
c=======================================================================
c   Choice of the start file(s) to use
c=======================================================================
      write(*,*) 'From which kind of files do you want to create new',
     .  'start and startfi files'
      write(*,*) '    0 - from a file start_archive'
      write(*,*) '    1 - from files start and startfi'
 
c-----------------------------------------------------------------------
c   Open file(s) to modify (start or start_archive)
c-----------------------------------------------------------------------

      DO
         read(*,*,iostat=ierr) choix_1
         if ((choix_1 /= 0).OR.(choix_1 /=1)) EXIT
      ENDDO

c     Open start_archive
c     ~~~~~~~~~~~~~~~~~~~~~~~~~~
      if (choix_1.eq.0) then

        write(*,*) 'Creating start files from:'
        write(*,*) './start_archive.nc'
        write(*,*)
        fichnom = 'start_archive.nc'
        ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
        IF (ierr.NE.NF_NOERR) THEN
          write(6,*)' Problem opening file:',fichnom
          write(6,*)' ierr = ', ierr
          CALL ABORT
        ENDIF
        tab0 = 50 
        Lmodif = 1

c     OR open start and startfi files
c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      else
        write(*,*) 'Creating start files from:'
        write(*,*) './start.nc and ./startfi.nc'
        write(*,*)
        fichnom = 'start.nc'
        ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_dyn)
        IF (ierr.NE.NF_NOERR) THEN
          write(6,*)' Problem opening file:',fichnom
          write(6,*)' ierr = ', ierr
          CALL ABORT
        ENDIF
 
        fichnom = 'startfi.nc'
        ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi)
        IF (ierr.NE.NF_NOERR) THEN
          write(6,*)' Problem opening file:',fichnom
          write(6,*)' ierr = ', ierr
          CALL ABORT
        ENDIF

        tab0 = 0 
        Lmodif = 0

      endif

c=======================================================================
c  INITIALISATIONS DIVERSES
c=======================================================================

! Take care of arrays in common modules
      ! ALLOCATE ARRAYS in surfdat_h (if not already done, e.g. when using start_archive)
      IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngridmx))
      IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngridmx))
      IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngridmx))
      IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngridmx))
      IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngridmx))
      IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngridmx))
      IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngridmx))

c-----------------------------------------------------------------------
c Lecture du tableau des parametres du run (pour la dynamique)
c-----------------------------------------------------------------------
      if (choix_1.eq.0) then

        write(*,*) 'reading tab_cntrl START_ARCHIVE'
c
        ierr = NF_INQ_VARID (nid, "controle", nvarid)
#ifdef NC_DOUBLE
        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
#else
        ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
#endif
c
      else if (choix_1.eq.1) then

        write(*,*) 'reading tab_cntrl START'
c
        ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid)
#ifdef NC_DOUBLE
        ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl)
#else
        ierr = NF_GET_VAR_REAL(nid_dyn, nvarid, tab_cntrl)
#endif
c
        write(*,*) 'reading tab_cntrl STARTFI'
c
        ierr = NF_INQ_VARID (nid_fi, "controle", nvarid)
#ifdef NC_DOUBLE
        ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis)
#else
        ierr = NF_GET_VAR_REAL(nid_fi, nvarid, tab_cntrl_bis)
#endif
c
        do i=1,50
          tab_cntrl(i+50)=tab_cntrl_bis(i)
        enddo
        write(*,*) 'printing tab_cntrl', tab_cntrl
        do i=1,100
          write(*,*) i,tab_cntrl(i)
        enddo
        
        write(*,*) 'Reading file START'
        fichnom = 'start.nc'
        CALL dynetat0(fichnom,vcov,ucov,teta,q,masse,
     .       ps,phis,time)

        CALL iniconst
        CALL inigeom

! Initialize the physics
         CALL iniphysiq(iim,jjm,llm,
     &                  (jjm-1)*iim+2,comm_lmdz,
     &                  daysec,day_ini,dtphys,
     &                  rlatu,rlatv,rlonu,rlonv,
     &                  aire,cu,cv,rad,g,r,cpp,
     &                  1)

! Initialize global tracer indexes (stored in tracer.h)
! ... this has to be done before phyetat0
         CALL initracer2(nqtot,tname)

        ! Lmodif set to 0 to disable modifications possibility in phyeta0                           
        Lmodif=0
        write(*,*) 'Reading file STARTFI'
        fichnom = 'startfi.nc'
        CALL phyetat0(.true.,ngridmx,llm,fichnom,tab0,Lmodif,nsoilmx,
     .        nqtot,day_ini,time,
     .        tsurf,tsoil,emis,q2,qsurf,tankCH4)

        ! copy albedo and soil thermal inertia on (local) physics grid
        do i=1,ngridmx
          albfi(i) = albedodat(i)
          do j=1,nsoilmx
           ithfi(i,j) = inertiedat(i,j)
          enddo
        ! build a surfithfi(:) using 1st layer of ithfi(:), which might
        ! be needed later on if reinitializing soil thermal inertia
          surfithfi(i)=ithfi(i,1)
        enddo
        ! also copy albedo and soil thermal inertia on (local) dynamics grid
        ! so that options below can manipulate either (but must then ensure
        ! to correctly recast things on physics grid)
        call gr_fi_dyn(1,ngridmx,iip1,jjp1,albfi,alb)
        call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
        call gr_fi_dyn(1,ngridmx,iip1,jjp1,surfithfi,surfith)
      
      endif
c-----------------------------------------------------------------------
c                Initialisation des constantes dynamique
c-----------------------------------------------------------------------

      ! JVO 2017 : Added a +1 shift in tab_cntrl indexes to be coherent
      ! with Earth-like Titan start files ( cf dynetat0 from common )

      kappa = tab_cntrl(10) 
      etot0 = tab_cntrl(13)
      ptot0 = tab_cntrl(14)
      ztot0 = tab_cntrl(15)
      stot0 = tab_cntrl(16)
      ang0 = tab_cntrl(17)
      write(*,*) "Newstart: kappa,etot0,ptot0,ztot0,stot0,ang0"
      write(*,*) kappa,etot0,ptot0,ztot0,stot0,ang0

      ! for vertical coordinate
      preff=tab_cntrl(19) ! reference surface pressure
      pa=tab_cntrl(18)  ! reference pressure at which coord is purely pressure
      !NB: in start_archive files tab_cntrl(17)=tab_cntrl(18)=0
      write(*,*) "Newstart: preff=",preff," pa=",pa
      yes=' '
      do while ((yes.ne.'y').and.(yes.ne.'n'))
        write(*,*) "Change the values of preff and pa? (y/n)"
        read(*,fmt='(a)') yes
      end do
      if (yes.eq.'y') then
        write(*,*)"New value of reference surface pressure preff?"
        write(*,*)"   (for Titan, typically preff=140000)"
        read(*,*) preff
        write(*,*)"New value of reference pressure pa for purely"
        write(*,*)"pressure levels (for hybrid coordinates)?"
        write(*,*)"   (for Titan, typically pa=10000)"
        read(*,*) pa
      endif
c-----------------------------------------------------------------------
c   Lecture du tab_cntrl et initialisation des constantes physiques
c  - pour start:  Lmodif = 0 => pas de modifications possibles
c                  (modif dans le tabfi de readfi + loin)
c  - pour start_archive:  Lmodif = 1 => modifications possibles
c-----------------------------------------------------------------------
      if (choix_1.eq.0) then
         ! tabfi requires that input file be first opened by open_startphy(fichnom)
         fichnom = 'start_archive.nc'
         call open_startphy(fichnom)
         call tabfi (ngridmx,nid,Lmodif,tab0,day_ini,lllm,p_rad,
     .            p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
      else if (choix_1.eq.1) then
         fichnom = 'startfi.nc'
         call open_startphy(fichnom)
         Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0                           
         call tabfi (ngridmx,nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
     .            p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
      endif

      if (p_omeg .eq. -9999.) then
        p_omeg = 8.*atan(1.)/p_daysec
        print*,"new value of omega",p_omeg
      endif

      rad = p_rad
      omeg = p_omeg
      g = p_g
      cpp = p_cpp
      mugaz = p_mugaz
      daysec = p_daysec

      if (p_omeg .eq. -9999.) then
        p_omeg = 8.*atan(1.)/p_daysec
        print*,"new value of omega",p_omeg
      endif

      kappa = 8.314*1000./(p_mugaz*p_cpp) ! added by RDW

c=======================================================================
c  INITIALISATIONS DIVERSES
c=======================================================================
! Load parameters from run.def file
      CALL defrun_new( 99, .TRUE. )
! Initialize dynamics geometry and co. (which, when using start.nc may 
! have changed because of modifications (tabi, preff, pa) above)
      CALL iniconst 
      CALL inigeom
      idum=-1
      idum=0

! Initialize the physics for start_archive only
      IF (choix_1.eq.0) THEN
         CALL iniphysiq(iim,jjm,llm,
     &                  (jjm-1)*iim+2,comm_lmdz,
     &                  daysec,day_ini,dtphys,
     &                  rlatu,rlatv,rlonu,rlonv,
     &                  aire,cu,cv,rad,g,r,cpp,
     &                  1)
! Initialize global tracer indexes (stored in tracer.h)
         CALL initracer2(nqtot,tname)

      ENDIF

c=======================================================================
c   lecture topographie, albedo, inertie thermique, relief sous-maille
c=======================================================================

      if (choix_1.eq.0) then  ! for start_archive files, 
                              ! where an external "surface.nc" file is needed

c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla'))
c       write(*,*)
c       write(*,*) 'choix du relief (mola,pla)'
c       write(*,*) '(Topographie MGS MOLA, plat)'
c       read(*,fmt='(a3)') relief
        relief="mola"
c     enddo

!    First get the correct value of datadir, if not already done:
        ! default 'datadir' is set in "datafile_mod"
        call getin_p("datadir",datadir)
        write(*,*) 'Available surface data files are:'
        filestring='ls '//trim(datadir)//'/'//
     &                    trim(surfdir)//' | grep .nc'
        call system(filestring)
        ! but in ye old days these files were in datadir, so scan it as well
        ! for the sake of retro-compatibility
        filestring='ls '//trim(datadir)//'/'//' | grep .nc'
        call system(filestring)

        write(*,*) ''
        write(*,*) 'Please choose the relevant file',
     &  ' (e.g. "surface_mars.nc")'
        write(*,*) ' or "none" to not use any (and set planetary'
        write(*,*) '  albedo and surface thermal inertia)'
        read(*,fmt='(a50)') surfacefile

        if (surfacefile.ne."none") then

          CALL datareadnc(relief,surfacefile,phis,alb,surfith,
     &          zmeaS,zstdS,zsigS,zgamS,ztheS)
        else
        ! specific case when not using a "surface.nc" file
          phis(:,:)=0
          zmeaS(:,:)=0
          zstdS(:,:)=0
          zsigS(:,:)=0
          zgamS(:,:)=0
          ztheS(:,:)=0
          
          write(*,*) "Enter value of albedo of the bare ground:"
          read(*,*) alb(1,1)
          alb(:,:)=alb(1,1)
          
          write(*,*) "Enter value of thermal inertia of soil:"
          read(*,*) surfith(1,1)
          surfith(:,:)=surfith(1,1)
        
        endif ! of if (surfacefile.ne."none")

        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi)
        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea)
        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd)
        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig)
        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)

      endif ! of if (choix_1.eq.0)


c=======================================================================
c  Lecture des fichiers (start ou start_archive)
c=======================================================================

      if (choix_1.eq.0) then

        write(*,*) 'Reading file START_ARCHIVE'
        CALL lect_start_archive(ngridmx,llm,
     &   date,tsurf,tsoil,emis,q2,
     &   t,ucov,vcov,ps,teta,phisold_newgrid,
     &   q,qsurf,tankCH4,surfith,nid)
        write(*,*) "OK, read start_archive file"
        ! copy soil thermal inertia
        ithfi(:,:)=inertiedat(:,:)
        
        ierr= NF_CLOSE(nid)

      else if (choix_1.eq.1) then 
         !do nothing, start and startfi have already been read
      else 
        CALL exit(1)
      endif

      dtvr   = daysec/FLOAT(day_step)
      dtphys   = dtvr * FLOAT(iphysiq)

c=======================================================================
c 
c=======================================================================

      do ! infinite loop on list of changes

      write(*,*)
      write(*,*)
      write(*,*) 'List of possible changes :'
      write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      write(*,*)
      write(*,*) 'flat : no topography ("aquaplanet")'
      write(*,*) 'set_ps_to_preff : used if changing preff with topo'
      write(*,*) 'bilball : uniform albedo and thermal inertia'
      write(*,*) 'qname : change tracer name'
      write(*,*) 't=profile  : read temperature profile in profile.in'
      write(*,*) 'q=0 : ALL tracer =zero'
      write(*,*) 'q=x : give a specific uniform value to one tracer'
      write(*,*) 'q=profile    : specify a profile for a tracer'
      write(*,*) 'isotherm  : Isothermal Temperatures, wind set to zero'
      write(*,*) 'radequi   : Earth-like radiative equilibrium temperature
     $ profile (lat-alt) and winds set to zero'
      write(*,*) 'ptot : change total pressure'
      write(*,*) 'emis : change surface emissivity'
      write(*,*) 'therm_ini_s : Set soil thermal inertia to reference sur
     &face values'
      write(*,*) 'slab_ocean_0 : initialisation of slab 
     $ocean variables'

        write(*,*)
        write(*,*) 'Change to perform ?'
        write(*,*) '   (enter keyword or return to end)'
        write(*,*)

        read(*,fmt='(a20)') modif
        if (modif(1:1) .eq. ' ')then
         exit ! exit loop on changes
        endif

        write(*,*)
        write(*,*) trim(modif) , ' : '

c       'flat : no topography ("aquaplanet")'
c       -------------------------------------
        if (trim(modif) .eq. 'flat') then
c         set topo to zero 
          z_reel(1:iip1,1:jjp1)=0
          phis(1:iip1,1:jjp1)=g*z_reel(1:iip1,1:jjp1)
          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
          write(*,*) 'topography set to zero.'
          write(*,*) 'WARNING : the subgrid topography parameters',
     &    ' were not set to zero ! => set calllott to F'                    

c        Choice of surface pressure
         yes=' '
         do while ((yes.ne.'y').and.(yes.ne.'n'))
            write(*,*) 'Do you wish to choose homogeneous surface',
     &                 'pressure (y) or let newstart interpolate ',
     &                 ' the previous field  (n)?'
             read(*,fmt='(a)') yes
         end do
         if (yes.eq.'y') then
           flagps0=.true.
           write(*,*) 'New value for ps (Pa) ?'
 201       read(*,*,iostat=ierr) patm
            if(ierr.ne.0) goto 201
            write(*,*) patm
            if (patm.eq.-9999.) then
              patm = preff
              write(*,*) "we set patm = preff", preff
            endif
             write(*,*)
             write(*,*) ' new ps everywhere (Pa) = ', patm
             write(*,*)
             do j=1,jjp1
               do i=1,iip1
                 ps(i,j)=patm
               enddo
             enddo
         end if

c       'set_ps_to_preff' : used if changing preff with topo  
c       ----------------------------------------------------
        else if (trim(modif) .eq. 'set_ps_to_preff') then
          do j=1,jjp1
           do i=1,iip1
             ps(i,j)=preff
           enddo
          enddo


c       bilball : uniform albedo, thermal inertia
c       -----------------------------------------
        else if (trim(modif) .eq. 'bilball') then
          write(*,*) 'constante albedo and iner.therm:'
          write(*,*) 'New value for albedo (ex: 0.25) ?'
 101      read(*,*,iostat=ierr) alb_bb
          if(ierr.ne.0) goto 101
          write(*,*)
          write(*,*) ' uniform albedo (new value):',alb_bb
          write(*,*)

          write(*,*) 'New value for thermal inertia (eg: 247) ?'
 102      read(*,*,iostat=ierr) ith_bb
          if(ierr.ne.0) goto 102
          write(*,*) 'uniform thermal inertia (new value):',ith_bb
          DO j=1,jjp1
             DO i=1,iip1
                alb(i,j) = alb_bb        ! albedo
                do isoil=1,nsoilmx
                  ith(i,j,isoil) = ith_bb        ! thermal inertia
                enddo
             END DO
          END DO
!          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)

c       ptot : Modification of the total pressure: ice + current atmosphere 
c       -------------------------------------------------------------------
        else if (trim(modif).eq.'ptot') then

c         calcul de la pression totale glace + atm actuelle
          patm=0.
          airetot=0.
          pcap=0.
          DO j=1,jjp1
             DO i=1,iim
                ig=1+(j-2)*iim +i
                if(j.eq.1) ig=1
                if(j.eq.jjp1) ig=ngridmx
                patm = patm + ps(i,j)*aire(i,j)
                airetot= airetot + aire(i,j)
             ENDDO
          ENDDO
          ptoto = pcap + patm

          print*,'Current total pressure at surface (atm) ',
     &       ptoto/airetot

          print*,'new value?'
          read(*,*) ptotn
          ptotn=ptotn*airetot
          patmn=ptotn-pcap
          print*,'ptoto,patm,ptotn,patmn'
          print*,ptoto,patm,ptotn,patmn
          print*,'Mult. factor for pressure (atm only)', patmn/patm
          do j=1,jjp1
             do i=1,iip1
                ps(i,j)=ps(i,j)*patmn/patm
             enddo
          enddo



c        Correction pour la conservation des traceurs
         yes=' '
         do while ((yes.ne.'y').and.(yes.ne.'n'))
            write(*,*) 'Do you wish to conserve tracer total mass (y)',
     &              ' or tracer mixing ratio (n) ?'
             read(*,fmt='(a)') yes
         end do

         if (yes.eq.'y') then
           write(*,*) 'OK : conservation of tracer total mass'
           DO iq =1, nqtot
             DO l=1,llm
               DO j=1,jjp1
                  DO i=1,iip1
                    q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
                  ENDDO
               ENDDO
             ENDDO
           ENDDO
          else
            write(*,*) 'OK : conservation of tracer mixing ratio'
          end if

c        Correction pour la pression au niveau de la mer
         yes=' '
         do while ((yes.ne.'y').and.(yes.ne.'n'))
            write(*,*) 'Do you wish fix pressure at sea level (y)',
     &              ' or not (n) ?'
             read(*,fmt='(a)') yes
         end do

         if (yes.eq.'y') then
           write(*,*) 'Value?'
                read(*,*,iostat=ierr) psea
             DO i=1,iip1
               DO j=1,jjp1
                    ps(i,j)=psea

                  ENDDO
               ENDDO
                write(*,*) 'psea=',psea
          else
            write(*,*) 'no'
          end if


c       emis : change surface emissivity (added by RW)
c       ----------------------------------------------
        else if (trim(modif).eq.'emis') then

           print*,'new value?'
           read(*,*) emisread

           do i=1,ngridmx
              emis(i)=emisread
           enddo

c       qname : change tracer name
c       --------------------------
        else if (trim(modif).eq.'qname') then
         yes='y'
         do while (yes.eq.'y')
          write(*,*) 'Which tracer name do you want to change ?'
          do iq=1,nqtot
            write(*,'(i3,a3,a20)')iq,' : ',trim(tname(iq))
          enddo
          write(*,'(a35,i3)')
     &            '(enter tracer number; between 1 and ',nqtot
          write(*,*)' or any other value to quit this option)'
          read(*,*) iq
          if ((iq.ge.1).and.(iq.le.nqtot)) then
            write(*,*)'Change tracer name ',trim(tname(iq)),' to ?'
            read(*,*) txt
            tname(iq)=txt
            write(*,*)'Do you want to change another tracer name (y/n)?'
            read(*,'(a)') yes 
          else
! inapropiate value of iq; quit this option
            yes='n'
          endif ! of if ((iq.ge.1).and.(iq.le.nqtot))
         enddo ! of do while (yes.ne.'y')

c       q=0 : set tracers to zero
c       -------------------------
        else if (trim(modif).eq.'q=0') then
c          mise a 0 des q (traceurs)
          write(*,*) 'Tracers set to 0 (1.E-30 in fact)'
           DO iq =1, nqtot
             DO l=1,llm
               DO j=1,jjp1
                  DO i=1,iip1
                    q(i,j,l,iq)=1.e-30
                  ENDDO
               ENDDO
             ENDDO
           ENDDO

c          set surface tracers to zero
           DO iq =1, nqtot
             DO ig=1,ngridmx
                 qsurf(ig,iq)=0.
             ENDDO
           ENDDO

c       q=x : initialise tracer manually 
c       --------------------------------
        else if (trim(modif).eq.'q=x') then
             write(*,*) 'Which tracer do you want to modify ?'
             do iq=1,nqtot
               write(*,*)iq,' : ',trim(tname(iq))
             enddo
             write(*,*) '(choose between 1 and ',nqtot,')'
             read(*,*) iq 
             write(*,*)'mixing ratio of tracer ',trim(tname(iq)),
     &                 ' ? (kg/kg)'
             read(*,*) val
             DO l=1,llm
               DO j=1,jjp1
                  DO i=1,iip1
                    q(i,j,l,iq)=val
                  ENDDO
               ENDDO
             ENDDO
             write(*,*) 'SURFACE value of tracer ',trim(tname(iq)),
     &                   ' ? (kg/m2)'
             read(*,*) val
             DO ig=1,ngridmx
                 qsurf(ig,iq)=val
             ENDDO

c       t=profile : initialize temperature with a given profile
c       -------------------------------------------------------
        else if (trim(modif) .eq. 't=profile') then
             write(*,*) 'Temperature profile from ASCII file'
             write(*,*) "'profile.in' e.g. 1D output"
             write(*,*) "(one value per line in file; starting with"
             write(*,*) "surface value, the 1st atmospheric layer"
             write(*,*) "followed by 2nd, etc. up to top of atmosphere)"
             txt="profile.in"
             open(33,file=trim(txt),status='old',form='formatted',
     &            iostat=ierr)
             if (ierr.eq.0) then
               ! OK, found file 'profile_...', load the profile
               do l=1,llm+1
                 read(33,*,iostat=ierr) profile(l)
                 write(*,*) profile(l)
                 if (ierr.ne.0) then ! something went wrong
                   exit ! quit loop
                 endif
               enddo
               if (ierr.eq.0) then
                 tsurf(1:ngridmx)=profile(1)
                 tsoil(1:ngridmx,1:nsoilmx)=profile(1)
                 do l=1,llm
                   Tset(1:iip1,1:jjp1,l)=profile(l+1)
                   flagtset=.true.
                 enddo
                 ucov(1:iip1,1:jjp1,1:llm)=0.
                 vcov(1:iip1,1:jjm,1:llm)=0.
                 q2(1:ngridmx,1:llm+1)=0.
               else
                 write(*,*)'problem reading file ',trim(txt),' !'
                 write(*,*)'No modifications to temperature'
               endif
             else
               write(*,*)'Could not find file ',trim(txt),' !'
               write(*,*)'No modifications to temperature'
             endif

c       q=profile : initialize tracer with a given profile
c       --------------------------------------------------
        else if (trim(modif) .eq. 'q=profile') then
             write(*,*) 'Tracer profile will be sought in ASCII file'
             write(*,*) "'profile_tracer' where 'tracer' is tracer name"
             write(*,*) "(one value per line in file; starting with"
             write(*,*) "surface value, the 1st atmospheric layer"
             write(*,*) "followed by 2nd, etc. up to top of atmosphere)"
             write(*,*) 'Which tracer do you want to set?'
             do iq=1,nqtot
               write(*,*)iq,' : ',trim(tname(iq))
             enddo
             write(*,*) '(choose between 1 and ',nqtot,')'
             read(*,*) iq 
             if ((iq.lt.1).or.(iq.gt.nqtot)) then
               ! wrong value for iq, go back to menu
               write(*,*) "wrong input value:",iq
               cycle
             endif
             ! look for input file 'profile_tracer'
             txt="profile_"//trim(tname(iq))
             open(41,file=trim(txt),status='old',form='formatted',
     &            iostat=ierr)
             if (ierr.eq.0) then
               ! OK, found file 'profile_...', load the profile
               do l=1,llm+1
                 read(41,*,iostat=ierr) profile(l)
                 if (ierr.ne.0) then ! something went wrong
                   exit ! quit loop
                 endif
               enddo
               if (ierr.eq.0) then
                 ! initialize tracer values
                 qsurf(:,iq)=profile(1)
                 do l=1,llm
                   q(:,:,l,iq)=profile(l+1)
                 enddo
                 write(*,*)'OK, tracer ',trim(tname(iq)),' initialized '
     &                    ,'using values from file ',trim(txt)
               else
                 write(*,*)'problem reading file ',trim(txt),' !'
                 write(*,*)'No modifications to tracer ',trim(tname(iq))
               endif
             else
               write(*,*)'Could not find file ',trim(txt),' !'
               write(*,*)'No modifications to tracer ',trim(tname(iq))
             endif
             

c       isotherm : Isothermal temperatures and no winds
c       -----------------------------------------------
        else if (trim(modif) .eq. 'isotherm') then

          write(*,*)'Isothermal temperature of the atmosphere, 
     &           surface and subsurface'
          write(*,*) 'Value of this temperature ? :'
 203      read(*,*,iostat=ierr) Tiso
          if(ierr.ne.0) goto 203

          tsurf(1:ngridmx)=Tiso
          
          tsoil(1:ngridmx,1:nsoilmx)=Tiso
          
          Tset(1:iip1,1:jjp1,1:llm)=Tiso
          flagtset=.true.

          t(1:iip1,1:jjp1,1:llm)=Tiso
          !! otherwise hydrost. integrations below
          !! use the wrong temperature
          !! -- NB: Tset might be useless
        
          ucov(1:iip1,1:jjp1,1:llm)=0
          vcov(1:iip1,1:jjm,1:llm)=0
          q2(1:ngridmx,1:llm+1)=0

c       radequi : Radiative equilibrium profile of temperatures and no winds
c       --------------------------------------------------------------------
        else if (trim(modif) .eq. 'radequi') then

          write(*,*)'radiative equilibrium temperature profile'       

          do ig=1, ngridmx
             teque= 335.0-60.0*sin(latfi(ig))*sin(latfi(ig))-
     &            10.0*cos(latfi(ig))*cos(latfi(ig))
             tsurf(ig) = MAX(220.0,teque)
          end do
          do l=2,nsoilmx
             do ig=1, ngridmx
                tsoil(ig,l) = tsurf(ig)
             end do
          end do
          DO j=1,jjp1
             DO i=1,iim
                Do l=1,llm
                   teque=335.-60.0*sin(rlatu(j))*sin(rlatu(j))-
     &                  10.0*cos(rlatu(j))*cos(rlatu(j))
                   Tset(i,j,l)=MAX(220.0,teque)
                end do
             end do
          end do
          flagtset=.true.
          ucov(1:iip1,1:jjp1,1:llm)=0
          vcov(1:iip1,1:jjm,1:llm)=0
          q2(1:ngridmx,1:llm+1)=0

!       therm_ini_s: (re)-set soil thermal inertia to reference surface values
!       ----------------------------------------------------------------------

        else if (trim(modif) .eq. 'therm_ini_s') then
!          write(*,*)"surfithfi(1):",surfithfi(1)
          do isoil=1,nsoilmx
            inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
          enddo
          write(*,*)'OK: Soil thermal inertia has been reset to referenc
     &e surface values'
!          write(*,*)"inertiedat(1,1):",inertiedat(1,1)
          ithfi(:,:)=inertiedat(:,:)
         ! recast ithfi() onto ith()
         call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
! Check:
!         do i=1,iip1
!           do j=1,jjp1
!             do isoil=1,nsoilmx
!               write(77,*) i,j,isoil,"  ",ith(i,j,isoil)
!             enddo
!           enddo
!         enddo



c       slab_ocean_initialisation
c       ------------------------------------------------
        else if (trim(modif) .eq. 'slab_ocean_0') then
        
           write(*,*) "No ocean yet on Titan ! Can't use this option"
           stop

        else
          write(*,*) '       Unknown (misspelled?) option!!!'
        end if ! of if (trim(modif) .eq. '...') elseif ...



       enddo ! of do ! infinite loop on liste of changes

 999  continue


c================================================================
c   Initialisation for methane surface tank (added by JVO 2017)
c================================================================
      DO ig=1, ngridmx
         tankCH4(ig)=2.0
      ENDDO

c=======================================================================
c   Correct pressure on the new grid (menu 0)
c=======================================================================


      if ((choix_1.eq.0).and.(.not.flagps0)) then
        r = 1000.*8.31/mugaz

        do j=1,jjp1
          do i=1,iip1
             ps(i,j) = ps(i,j) *
     .            exp((phisold_newgrid(i,j)-phis(i,j)) /
     .                                  (t(i,j,1) * r))
          end do
        end do

c periodicite de ps en longitude
        do j=1,jjp1
          ps(1,j) = ps(iip1,j)
        end do
      end if

         
c=======================================================================
c=======================================================================

c=======================================================================
c    Initialisation de la physique / ecriture de newstartfi :
c=======================================================================


      CALL inifilr 
      CALL pression(ip1jmp1, ap, bp, ps, p3d)
      

c=========================================================================
c  Calcul de la dimension verticale pour la chimie  - JVO 2017
c  start_archive seulement, la grille verticale pouvant avoir ete modifiee
c==========================================================================

      IF (choix_1.eq.0) THEN

        ! Calculate the # of upper chemistry layers with the "new" pressure grid
        ! For this we use Vervack profile for upper atmosphere with dz=10km

        CALL gr_kim_vervack

        WRITE(*,*)
        WRITE(*,*) " With the compiled vertical grid we found :"
        WRITE(*,*) " Number of upper chemistry layers =", nlaykim_up
        
        IF (callchim) THEN
            
          IF (.NOT.allocated(ykim_up)) THEN
            ALLOCATE(ykim_up(nkim,ngridmx,nlaykim_up))
          ENDIF
          
          ! Regridind of upper chemistry fields and uppermost GCM
          ! chemistry tracers if needed
          
          IF (nlaykimold.ne.nlaykim_up) THEN

            WRITE(*,*) " Warning, nlaykimold=", nlaykimold
            WRITE(*,*) '  which implies that a vertical regrid on upper 
     &chemistry fields will be performed. If the ceiling is lowered I
     &will also adjust GCM uppermost layers chem. tracers.'
            WRITE(*,*)
          
            CALL vert_regrid_kim(nqtot,q)

          ELSE
            ykim_up(:,:,:) = ykim_up_oldv(:,:,:) ! No vertical recast
          ENDIF
          
        ENDIF
 
      endif ! of if (choix_1.eq.0)      
 

c-----------------------------------------------------------------------
c   Initialisation  pks:
c-----------------------------------------------------------------------

      CALL exner_hyb(ip1jmp1, ps, p3d, pks, pk, pkf)
! Calcul de la temperature potentielle teta

      if (flagtset) then
          DO l=1,llm
             DO j=1,jjp1
                DO i=1,iim
                   teta(i,j,l) = Tset(i,j,l) * cpp/pk(i,j,l)
                ENDDO
                teta (iip1,j,l)= teta (1,j,l)
             ENDDO
          ENDDO
      else if (choix_1.eq.0) then
         DO l=1,llm
            DO j=1,jjp1
               DO i=1,iim
                  teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
               ENDDO
               teta (iip1,j,l)= teta (1,j,l)
            ENDDO
         ENDDO
      endif

C Calcul intermediaire
c
      if (choix_1.eq.0) then
         CALL massdair( p3d, masse  )
c
         print *,' ALPHAX ',alphax

         DO  l = 1, llm
           DO  i    = 1, iim
             xppn(i) = aire( i, 1   ) * masse(  i     ,  1   , l )
             xpps(i) = aire( i,jjp1 ) * masse(  i     , jjp1 , l )
           ENDDO
             xpn      = SUM(xppn)/apoln
             xps      = SUM(xpps)/apols
           DO i   = 1, iip1
             masse(   i   ,   1     ,  l )   = xpn
             masse(   i   ,   jjp1  ,  l )   = xps
           ENDDO
         ENDDO
      endif
      phis(iip1,:) = phis(1,:)

      itau=0
      if (choix_1.eq.0) then
         day_ini=int(date)
      endif

c
      CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )

      CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
     *                phi,w, pbaru,pbarv,day_ini+time )

          
      CALL dynredem0("restart.nc",day_ini,phis)
      CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,masse,ps) 
C
C Ecriture etat initial physique
C

      call physdem0("restartfi.nc",lonfi,latfi,nsoilmx,ngridmx,llm,
     &              nqtot,dtphys,real(day_ini),0.0,
     &              airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
      call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
     &                dtphys,real(day_ini),
     &                tsurf,tsoil,emis,q2,qsurf,tankCH4)

c=======================================================================
c        Formats 
c=======================================================================

   1  FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema
     *rrage est differente de la valeur parametree iim =',i4//)
   2  FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema
     *rrage est differente de la valeur parametree jjm =',i4//)
   3  FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar
     *rage est differente de la valeur parametree llm =',i4//)

      write(*,*) "newstart: All is well that ends well."

      end

