      SUBROUTINE inifis(ngrid,nlayer,
     $           day_ini,pdaysec,ptimestep,
     $           plat,plon,parea,
     $           prad,pg,pr,pcpp)

      use radinc_h, only : naerkind
      use datafile_mod, only: datadir
      use comdiurn_h

      !! to be conservative, but why this include here?
      use surfdat_h
      use comsaison_h

      USE comgeomfi_h

!=======================================================================
!
!   purpose:
!   -------
!
!   Initialisation for the physical parametrisations of the LMD 
!   Generic Model.
!
!   author: Frederic Hourdin 15 / 10 /93
!   -------
!   modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
!             Ehouarn Millour (oct. 2008) tracers are now identified
!              by their names and may not be contiguously
!              stored in the q(:,:,:,:) array
!             E.M. (june 2009) use getin routine to load parameters
!
!
!   arguments:
!   ----------
!
!   input:
!   ------
!
!    ngrid                 Size of the horizontal grid.
!                          All internal loops are performed on that grid.
!    nlayer                Number of vertical layers.
!    pdayref               Day of reference for the simulation
!    pday                  Number of days counted from the North. Spring
!                          equinoxe.
!
!=======================================================================
!
!-----------------------------------------------------------------------
!   declarations:
!   -------------
      use datafile_mod, only: datadir
! to use  'getin'
      USE ioipsl_getincom 
      IMPLICIT NONE
#include "dimensions.h"
#include "dimphys.h"
#include "planete.h"
#include "comcstfi.h"
#include "callkeys.h"



      REAL prad,pg,pr,pcpp,pdaysec,ptimestep
 
      INTEGER ngrid,nlayer
      REAL plat(ngrid),plon(ngrid),parea(ngrid)
      integer day_ini
      INTEGER ig,ierr
 
      EXTERNAL iniorbit,orbite
      EXTERNAL SSUM
      REAL SSUM
 
      CHARACTER ch1*12
      CHARACTER ch80*80

      logical chem, h2o
      logical :: parameter, doubleq=.false.

      real psurf,pN2 ! added by RW for Gliese 581d N2+CO2

      rad=prad
      daysec=pdaysec
      dtphys=ptimestep
      cpp=pcpp
      g=pg
      r=pr
      rcp=r/cpp

      avocado = 6.02214179e23   ! added by RW

      cursor = 1 ! added by AS in dimphys. 1 for sequential runs.

! --------------------------------------------------------
!     The usual Tests
!     --------------
      IF (nlayer.NE.nlayermx) THEN
         PRINT*,'STOP in inifis'
         PRINT*,'Probleme de dimensions :'
         PRINT*,'nlayer     = ',nlayer
         PRINT*,'nlayermx   = ',nlayermx
         STOP
      ENDIF

! --------------------------------------------------------------
!  Reading the "callphys.def" file controlling some key options
! --------------------------------------------------------------
     
      ! check that 'callphys.def' file is around
      OPEN(99,file='callphys.def',status='old',form='formatted'
     &     ,iostat=ierr)
      CLOSE(99)
      
      IF(ierr.EQ.0) THEN
         PRINT*
         PRINT*
         PRINT*,'--------------------------------------------'
         PRINT*,' inifis: Parametres pour la physique (callphys.def)'
         PRINT*,'--------------------------------------------'

         write(*,*) "Directory where external input files are:"
         ! default 'datadir' is set in "datadir_mod"
         call getin("datadir",datadir) ! default path
         write(*,*) " datadir = ",trim(datadir)

         write(*,*) "Run with or without tracer transport ?"
         tracer=.false. ! default value
         call getin("tracer",tracer)
         write(*,*) " tracer = ",tracer

         write(*,*) "Run with or without atm mass update ",
     &      " due to tracer evaporation/condensation?"
         mass_redistrib=.false. ! default value
         call getin("mass_redistrib",mass_redistrib)
         write(*,*) " mass_redistrib = ",mass_redistrib

         write(*,*) "Diurnal cycle ?"
         write(*,*) "(if diurnal=false, diurnal averaged solar heating)"
         diurnal=.true. ! default value
         call getin("diurnal",diurnal)
         write(*,*) " diurnal = ",diurnal

         write(*,*) "Seasonal cycle ?"
         write(*,*) "(if season=false, Ls stays constant, to value ",
     &   "set in 'start'"
         season=.true. ! default value
         call getin("season",season)
         write(*,*) " season = ",season

         write(*,*) "Tidally resonant rotation ?"
         tlocked=.false. ! default value
         call getin("tlocked",tlocked)
         write(*,*) "tlocked = ",tlocked

! Test of incompatibility:
! if tlocked, then diurnal should be false
         if (tlocked.and.diurnal) then
           print*,'If diurnal=true, we should turn off tlocked.'
           stop
         endif

         write(*,*) "Tidal resonance ratio ?"
         nres=0          ! default value
         call getin("nres",nres)
         write(*,*) "nres = ",nres

         write(*,*) "Write some extra output to the screen ?"
         lwrite=.false. ! default value
         call getin("lwrite",lwrite)
         write(*,*) " lwrite = ",lwrite

         write(*,*) "Save statistics in file stats.nc ?"
         callstats=.true. ! default value
         call getin("callstats",callstats)
         write(*,*) " callstats = ",callstats

         write(*,*) "Test energy conservation of model physics ?"
         enertest=.false. ! default value
         call getin("enertest",enertest)
         write(*,*) " enertest = ",enertest

         write(*,*) "Check to see if cpp values used match gases.def ?"
         check_cpp_match=.true. ! default value
         call getin("check_cpp_match",check_cpp_match)
         write(*,*) " check_cpp_match = ",check_cpp_match

         write(*,*) "call radiative transfer ?"
         callrad=.true. ! default value
         call getin("callrad",callrad)
         write(*,*) " callrad = ",callrad

         write(*,*) "call correlated-k radiative transfer ?"
         corrk=.true. ! default value
         call getin("corrk",corrk)
         write(*,*) " corrk = ",corrk

         write(*,*) "call gaseous absorption in the visible bands?",
     &              "(matters only if callrad=T)"
         callgasvis=.false. ! default value
         call getin("callgasvis",callgasvis)
         write(*,*) " callgasvis = ",callgasvis
        
         write(*,*) "call continuum opacities in radiative transfer ?",
     &              "(matters only if callrad=T)"
         continuum=.true. ! default value
         call getin("continuum",continuum)
         write(*,*) " continuum = ",continuum

         write(*,*) "use analytic function for H2O continuum ?"
         H2Ocont_simple=.false. ! default value
         call getin("H2Ocont_simple",H2Ocont_simple)
         write(*,*) " H2Ocont_simple = ",H2Ocont_simple
 
         write(*,*) "call turbulent vertical diffusion ?"
         calldifv=.true. ! default value
         call getin("calldifv",calldifv)
         write(*,*) " calldifv = ",calldifv

         write(*,*) "use turbdiff instead of vdifc ?"
         UseTurbDiff=.true. ! default value
         call getin("UseTurbDiff",UseTurbDiff)
         write(*,*) " UseTurbDiff = ",UseTurbDiff

         write(*,*) "call convective adjustment ?"
         calladj=.true. ! default value
         call getin("calladj",calladj)
         write(*,*) " calladj = ",calladj

         write(*,*) "call CO2 condensation ?"
         co2cond=.false. ! default value
         call getin("co2cond",co2cond)
         write(*,*) " co2cond = ",co2cond
! Test of incompatibility
         if (co2cond.and.(.not.tracer)) then
            print*,'We need a CO2 ice tracer to condense CO2'
            call abort
         endif 
 
         write(*,*) "CO2 supersaturation level ?"
         co2supsat=1.0 ! default value
         call getin("co2supsat",co2supsat)
         write(*,*) " co2supsat = ",co2supsat

         write(*,*) "Radiative timescale for Newtonian cooling ?"
         tau_relax=30. ! default value
         call getin("tau_relax",tau_relax)
         write(*,*) " tau_relax = ",tau_relax
         tau_relax=tau_relax*24*3600 ! convert Earth days --> seconds

         write(*,*)"call thermal conduction in the soil ?"
         callsoil=.true. ! default value
         call getin("callsoil",callsoil)
         write(*,*) " callsoil = ",callsoil
         
         write(*,*)"Rad transfer is computed every iradia",
     &             " physical timestep"
         iradia=1 ! default value
         call getin("iradia",iradia)
         write(*,*)" iradia = ",iradia
       
         write(*,*)"Rayleigh scattering ?"
         rayleigh=.false.
         call getin("rayleigh",rayleigh)
         write(*,*)" rayleigh = ",rayleigh

         write(*,*) "Use blackbody for stellar spectrum ?"
         stelbbody=.false. ! default value
         call getin("stelbbody",stelbbody)
         write(*,*) " stelbbody = ",stelbbody

         write(*,*) "Stellar blackbody temperature ?"
         stelTbb=5800.0 ! default value
         call getin("stelTbb",stelTbb)
         write(*,*) " stelTbb = ",stelTbb

         write(*,*)"Output mean OLR in 1D?"
         meanOLR=.false.
         call getin("meanOLR",meanOLR)
         write(*,*)" meanOLR = ",meanOLR

         write(*,*)"Output spectral OLR in 3D?"
         specOLR=.false.
         call getin("specOLR",specOLR)
         write(*,*)" specOLR = ",specOLR

         write(*,*)"Operate in kastprof mode?"
         kastprof=.false.
         call getin("kastprof",kastprof)
         write(*,*)" kastprof = ",kastprof

         write(*,*)"Uniform absorption in radiative transfer?"
         graybody=.false.
         call getin("graybody",graybody)
         write(*,*)" graybody = ",graybody

! Test of incompatibility:
! if kastprof used, we must be in 1D
         if (kastprof.and.(ngrid.gt.1)) then
           print*,'kastprof can only be used in 1D!'
           call abort
         endif

         write(*,*)"Stratospheric temperature for kastprof mode?"
         Tstrat=167.0
         call getin("Tstrat",Tstrat)
         write(*,*)" Tstrat = ",Tstrat

         write(*,*)"Remove lower boundary?"
         noradsurf=.false.
         call getin("noradsurf",noradsurf)
         write(*,*)" noradsurf = ",noradsurf

! Tests of incompatibility:
         if (noradsurf.and.callsoil) then
           print*,'noradsurf not compatible with soil scheme!'
           call abort
         endif
         !if (noradsurf.and.calldifv) then
         !  print*,'noradsurf not compatible with a boundary layer!'
         !  call abort
         !endif

         write(*,*)"Use Newtonian cooling for radiative transfer?"
         newtonian=.false.
         call getin("newtonian",newtonian)
         write(*,*)" newtonian = ",newtonian

! Tests of incompatibility:
         if (newtonian.and.corrk) then
           print*,'newtonian not compatible with correlated-k!'
           call abort
         endif
         if (newtonian.and.calladj) then
           print*,'newtonian not compatible with adjustment!'
           call abort
         endif
         if (newtonian.and.calldifv) then
           print*,'newtonian not compatible with a boundary layer!'
           call abort
         endif

         write(*,*)"Test physics timescale in 1D?"
         testradtimes=.false.
         call getin("testradtimes",testradtimes)
         write(*,*)" testradtimes = ",testradtimes

! Test of incompatibility:
! if testradtimes used, we must be in 1D
         if (testradtimes.and.(ngrid.gt.1)) then
           print*,'testradtimes can only be used in 1D!'
           call abort
         endif

         write(*,*)"Default planetary temperature?"
         tplanet=215.0
         call getin("tplanet",tplanet)
         write(*,*)" tplanet = ",tplanet

         write(*,*)"Which star?"
         startype=1 ! default value = Sol
         call getin("startype",startype)
         write(*,*)" startype = ",startype

         write(*,*)"Value of stellar flux at 1 AU?"
         Fat1AU=1356.0 ! default value = Sol today
         call getin("Fat1AU",Fat1AU)
         write(*,*)" Fat1AU = ",Fat1AU


! TRACERS:

         write(*,*)"Varying H2O cloud fraction?"
         CLFvarying=.false.     ! default value
         call getin("CLFvarying",CLFvarying)
         write(*,*)" CLFvarying = ",CLFvarying

         write(*,*)"Value of fixed H2O cloud fraction?"
         CLFfixval=1.0                ! default value
         call getin("CLFfixval",CLFfixval)
         write(*,*)" CLFfixval = ",CLFfixval

         write(*,*)"fixed radii for Cloud particles?"
         radfixed=.false. ! default value
         call getin("radfixed",radfixed)
         write(*,*)" radfixed = ",radfixed

         if(kastprof)then
            radfixed=.true.
         endif  

	 write(*,*)"Number mixing ratio of CO2 ice particles:"
         Nmix_co2=1.e6 ! default value
         call getin("Nmix_co2",Nmix_co2)
         write(*,*)" Nmix_co2 = ",Nmix_co2

!         write(*,*)"Number of radiatively active aerosols:"
!         naerkind=0. ! default value
!         call getin("naerkind",naerkind)
!         write(*,*)" naerkind = ",naerkind

         write(*,*)"Opacity of dust (if used):"
         dusttau=0. ! default value
         call getin("dusttau",dusttau)
         write(*,*)" dusttau = ",dusttau

         write(*,*)"Radiatively active CO2 aerosols?"
         aeroco2=.false.     ! default value
         call getin("aeroco2",aeroco2)
         write(*,*)" aeroco2 = ",aeroco2

         write(*,*)"Fixed CO2 aerosol distribution?"
         aerofixco2=.false.     ! default value
         call getin("aerofixco2",aerofixco2)
         write(*,*)" aerofixco2 = ",aerofixco2

         write(*,*)"Radiatively active water ice?"
         aeroh2o=.false.     ! default value
         call getin("aeroh2o",aeroh2o)
         write(*,*)" aeroh2o = ",aeroh2o

         write(*,*)"Fixed H2O aerosol distribution?"
         aerofixh2o=.false.     ! default value
         call getin("aerofixh2o",aerofixh2o)
         write(*,*)" aerofixh2o = ",aerofixh2o

         write(*,*)"Radiatively active sulfuric acid aersols?"
         aeroh2so4=.false.     ! default value
         call getin("aeroh2so4",aeroh2so4)
         write(*,*)" aeroh2so4 = ",aeroh2so4

         write(*,*)"Cloud pressure level (with kastprof only):"
         cloudlvl=0. ! default value
         call getin("cloudlvl",cloudlvl)
         write(*,*)" cloudlvl = ",cloudlvl

         write(*,*)"Is the variable gas species radiatively active?"
         Tstrat=167.0
         varactive=.false.
         call getin("varactive",varactive)
         write(*,*)" varactive = ",varactive

         write(*,*)"Is the variable gas species distribution set?"
         varfixed=.false.
         call getin("varfixed",varfixed)
         write(*,*)" varfixed = ",varfixed

         write(*,*)"What is the saturation % of the variable species?"
         satval=0.8
         call getin("satval",satval)
         write(*,*)" satval = ",satval


! Test of incompatibility:
! if varactive, then varfixed should be false
         if (varactive.and.varfixed) then
           print*,'if varactive, varfixed must be OFF!'
           stop
         endif

         write(*,*) "Gravitationnal sedimentation ?"
         sedimentation=.false. ! default value
         call getin("sedimentation",sedimentation)
         write(*,*) " sedimentation = ",sedimentation

         write(*,*) "Compute water cycle ?"
         water=.false. ! default value
         call getin("water",water)
         write(*,*) " water = ",water
         
! Test of incompatibility:
! if water is true, there should be at least a tracer
         if (water.and.(.not.tracer)) then
           print*,'if water is ON, tracer must be ON too!'
           stop
         endif

         write(*,*) "Include water condensation ?"
         watercond=.false. ! default value
         call getin("watercond",watercond)
         write(*,*) " watercond = ",watercond

! Test of incompatibility:
! if watercond is used, then water should be used too
         if (watercond.and.(.not.water)) then
           print*,'if watercond is used, water should be used too'
           stop
         endif

         write(*,*) "Include water precipitation ?"
         waterrain=.false. ! default value
         call getin("waterrain",waterrain)
         write(*,*) " waterrain = ",waterrain

         write(*,*) "Include surface hydrology ?"
         hydrology=.false. ! default value
         call getin("hydrology",hydrology)
         write(*,*) " hydrology = ",hydrology

         write(*,*) "Evolve surface water sources ?"
         sourceevol=.false. ! default value
         call getin("sourceevol",sourceevol)
         write(*,*) " sourceevol = ",sourceevol

         write(*,*) "Ice evolution timestep ?"
         icetstep=100.0 ! default value
         call getin("icetstep",icetstep)
         write(*,*) " icetstep = ",icetstep

         write(*,*) "Snow albedo ?"
         albedosnow=0.5         ! default value
         call getin("albedosnow",albedosnow)
         write(*,*) " albedosnow = ",albedosnow

         write(*,*) "Maximum ice thickness ?"
         maxicethick=2.0         ! default value
         call getin("maxicethick",maxicethick)
         write(*,*) " maxicethick = ",maxicethick

         write(*,*) "Freezing point of seawater ?"
         Tsaldiff=-1.8          ! default value
         call getin("Tsaldiff",Tsaldiff)
         write(*,*) " Tsaldiff = ",Tsaldiff

         write(*,*) "Does user want to force cpp and mugaz?"
         force_cpp=.false. ! default value
         call getin("force_cpp",force_cpp)
         write(*,*) " force_cpp = ",force_cpp

         if (force_cpp) then
           mugaz = -99999.
           PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
           call getin("mugaz",mugaz)
           IF (mugaz.eq.-99999.) THEN
               PRINT *, "mugaz must be set if force_cpp = T"
               STOP
           ELSE
               write(*,*) "mugaz=",mugaz
           ENDIF
           cpp = -99999.
           PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
           call getin("cpp",cpp)
           IF (cpp.eq.-99999.) THEN
               PRINT *, "cpp must be set if force_cpp = T"
               STOP
           ELSE
               write(*,*) "cpp=",cpp
           ENDIF
         else
           mugaz=8.314*1000./pr
         endif
         call su_gases
         call calc_cpp_mugaz

         PRINT*,'--------------------------------------------'
         PRINT*
         PRINT*
      ELSE
         write(*,*)
         write(*,*) 'Cannot read file callphys.def. Is it here ?'
         stop
      ENDIF

8000  FORMAT(t5,a12,l8)
8001  FORMAT(t5,a12,i8)

      PRINT*
      PRINT*,'inifis: daysec',daysec
      PRINT*
      PRINT*,'inifis: The radiative transfer is computed:'
      PRINT*,'           each ',iradia,' physical time-step'
      PRINT*,'        or each ',iradia*dtphys,' seconds'
      PRINT*


!-----------------------------------------------------------------------
!     Some more initialization:
!     ------------------------

      ! ALLOCATE ARRAYS IN comgeomfi_h
      IF (.not. ALLOCATED(lati)) ALLOCATE(lati(ngrid))
      IF (.not. ALLOCATED(long)) ALLOCATE(long(ngrid))
      IF (.not. ALLOCATED(area)) ALLOCATE(area(ngrid))

      CALL SCOPY(ngrid,plon,1,long,1)
      CALL SCOPY(ngrid,plat,1,lati,1)
      CALL SCOPY(ngrid,parea,1,area,1)
      totarea=SSUM(ngrid,area,1)

      !! those are defined in comdiurn_h.F90
      IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
      IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
      IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
      IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))

      DO ig=1,ngrid
         sinlat(ig)=sin(plat(ig))
         coslat(ig)=cos(plat(ig))
         sinlon(ig)=sin(plon(ig))
         coslon(ig)=cos(plon(ig))
      ENDDO

      pi=2.*asin(1.) ! NB: pi is a common in comcstfi.h

      RETURN
      END
