SUBROUTINE inichim_readcallphys
!
!=======================================================================
!
!   subject:
!   --------
!
!   Initialisation for the physical parametrisations of the LMD 
!   martian atmospheric general circulation model.
!   This routine is called by inichim_newstart (which is used by programs
!   newstart and testphys1d but not the GCM)
!
!   author: Frederic Hourdin 15 / 10 /93
!   -------
!   modifications: Sebastien Lebonnois 11/06/2003 (new callphys.def)
!                  10/2008 adapted to using tracer sby name. Ehouarn Millour
!                  07/2009 use 'getin' to read callphys.def. EM
!
!   arguments:
!   ----------
!
!   input:
!   ------
!
!=======================================================================
!
!-----------------------------------------------------------------------
!   declarations:
!   -------------
! to use  'getin'
      USE ioipsl_getincom 
      IMPLICIT NONE
#include "dimensions.h"
#include "dimphys.h"
#include "planete.h"
#include "comcstfi.h"
#include "comsaison.h"
#include "comdiurn.h"
#include "comgeomfi.h"
#include "callkeys.h"
#include "surfdat.h"

      character*12 ch1
      integer ierr
      logical chem, h2o


! --------------------------------------------------------------
!  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*,' inichim_readcallphys: Parametres pour la physique',
     &          ' (callphys.def)'
         PRINT*,'--------------------------------------------'


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

         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(*,*) "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(*,*) "Save EOF profiles in file 'profiles' for ",
     &              "Climate Database?"
         calleofdump=.false. ! default value
         call getin("calleofdump",calleofdump)
         write(*,*) " calleofdump = ",calleofdump

         write(*,*) "Dust scenario:"
         iaervar=3 ! default value
         call getin("iaervar",iaervar)
         write(*,*) " iaervar = ",iaervar

         write(*,*) "Dust vertical distribution:"
         write(*,*) "(=1 Dust opt.deph read in startfi;",
     & " =2 Viking scenario; =3 MGS scenario,",
     & " =4 Mars Year 24 from TES assimilation)"
         iddist=3 ! default value
         call getin("iddist",iddist)
         write(*,*) " iddist = ",iddist

         write(*,*) "Dust top altitude (km). (Matters only if iddist=1)"
         topdustref= 90.0 ! default value
         call getin("topdustref",topdustref)
         write(*,*) " topdustref = ",topdustref


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

         write(*,*) "call NLTE radiative schemes ?",
     &              "(matters only if callrad=T)"
         callnlte=.false. ! default value
         call getin("callnlte",callnlte)
         write(*,*) " callnlte = ",callnlte
         
         write(*,*) "call CO2 NIR absorption ?",
     &              "(matters only if callrad=T)"
         callnirco2=.false. ! default value
         call getin("callnirco2",callnirco2)
         write(*,*) " callnirco2 = ",callnirco2
         
         write(*,*) "call turbulent vertical diffusion ?"
         calldifv=.true. ! default value
         call getin("calldifv",calldifv)
         write(*,*) " calldifv = ",calldifv

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

         write(*,*) "call CO2 condensation ?"
         callcond=.true. ! default value
         call getin("callcond",callcond)
         write(*,*) " callcond = ",callcond

         write(*,*)"call thermal conduction in the soil ?"
         callsoil=.true. ! default value
         call getin("callsoil",callsoil)
         write(*,*) " callsoil = ",callsoil
         

         write(*,*)"call Lott's gravity wave/subgrid topography ",
     &             "scheme ?"
         calllott=.true. ! default value
         call getin("calllott",calllott)
         write(*,*)" calllott = ",calllott


         write(*,*)"rad.transfer is computed every iradia",
     &             " physical timestep"
         iradia=1 ! default value
         call getin("iradia",iradia)
         write(*,*)" iradia = ",iradia
         

         write(*,*)"Output of the exchange coefficient mattrix ?",
     &             "(for diagnostics only)"
         callg2d=.false. ! default value
         call getin("callg2d",callg2d)
         write(*,*)" callg2d = ",callg2d

         write(*,*)"Rayleigh scattering : (should be .false. for now)"
         rayleigh=.false.
         call getin("rayleigh",rayleigh)
         write(*,*)" rayleigh = ",rayleigh


! TRACERS:

         write(*,*)"Transported dust ? (if >0, use 'dustbin' dust bins)"
         dustbin=0 ! default value
         call getin("dustbin",dustbin)
         write(*,*)" dustbin = ",dustbin

         write(*,*)"Radiatively active dust ? (matters if dustbin>0)"
         active=.false. ! default value
         call getin("active",active)
         write(*,*)" active = ",active

! Test of incompatibility:
! if active is used, then dustbin should be > 0

         if (active.and.(dustbin.lt.1)) then
           print*,'if active is used, then dustbin should > 0'
           stop
         endif

         write(*,*)"use mass and number mixing ratios to predict",
     &             " dust size ?"
         doubleq=.false. ! default value
         call getin("doubleq",doubleq)
         write(*,*)" doubleq = ",doubleq

! Test of incompatibility:
! if doubleq is used, then dustbin should be 1

         if (doubleq.and.(dustbin.ne.1)) then
           print*,'if doubleq is used, then dustbin should be 1'
           stop
         endif

         write(*,*)"dust lifted by GCM surface winds ?"
         lifting=.false. ! default value
         call getin("lifting",lifting)
         write(*,*)" lifting = ",lifting

! Test of incompatibility:
! if lifting is used, then dustbin should be > 0

         if (lifting.and.(dustbin.lt.1)) then
           print*,'if lifting is used, then dustbin should > 0'
           stop
         endif

         write(*,*)" dust lifted by dust devils ?"
         callddevil=.false. !default value
         call getin("callddevil",callddevil)
         write(*,*)" callddevil = ",callddevil
         

! Test of incompatibility:
! if dustdevil is used, then dustbin should be > 0

         if (callddevil.and.(dustbin.lt.1)) then
           print*,'if dustdevil is used, then dustbin should > 0'
           stop
         endif

         write(*,*)"Dust scavenging by CO2 snowfall ?"
         scavenging=.false. ! default value
         call getin("scavenging",scavenging)
         write(*,*)" scavenging = ",scavenging
         

! Test of incompatibility:
! if scavenging is used, then dustbin should be > 0

         if (scavenging.and.(dustbin.lt.1)) then
           print*,'if scavenging is used, then dustbin should > 0'
           stop
         endif

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

         write(*,*) "includes water ice",
     &              "(if true, 'water' must also be .true.)"

         write(*,*) "Radiatively active transported atmospheric ",
     &              "water ice ?"
         activice=.false. ! default value
         call getin("activice",activice)
         write(*,*) " activice = ",activice

         write(*,*) "Compute water cycle ?"
         water=.false. ! default value
         call getin("water",water)
         write(*,*) " water = ",water
         write(*,*) "Permanent water caps at poles ?",
     &               " .true. is RECOMMENDED"
         write(*,*) "(with .true., North cap is a source of water ",
     &   "and South pole is a cold trap)"
         caps=.true. ! default value
         call getin("caps",caps)
         write(*,*) " caps = ",caps

! Test of incompatibility:
! if activice is used, then water should be used too

         if (activice.and..not.water) then
           print*,'if activice is used, water should be used too'
           stop
         endif

         write(*,*) "photochemistry: include chemical species"
         photochem=.false. ! default value
         call getin("photochem",photochem)
         write(*,*) " photochem = ",photochem


! THERMOSPHERE

         write(*,*) "call thermosphere ?"
         callthermos=.false. ! default value
         call getin("callthermos",callthermos)
         write(*,*) " callthermos = ",callthermos
         

         write(*,*) " water included without cycle ",
     &              "(only if water=.false.)"
         thermoswater=.false. ! default value
         call getin("thermoswater",thermoswater)
         write(*,*) " thermoswater = ",thermoswater

         write(*,*) "call thermal conduction ?",
     &    " (only if callthermos=.true.)"
         callconduct=.false. ! default value
         call getin("callconduct",callconduct)
         write(*,*) " callconduct = ",callconduct

         write(*,*) "call EUV heating ?",
     &   " (only if callthermos=.true.)"
         calleuv=.false.  ! default value
         call getin("calleuv",calleuv)
         write(*,*) " calleuv = ",calleuv

         write(*,*) "call molecular viscosity ?",
     &   " (only if callthermos=.true.)"
         callmolvis=.false. ! default value
         call getin("callmolvis",callmolvis)
         write(*,*) " callmolvis = ",callmolvis

         write(*,*) "call molecular diffusion ?",
     &   " (only if callthermos=.true.)"
         callmoldiff=.false. ! default value
         call getin("callmoldiff",callmoldiff)
         write(*,*) " callmoldiff = ",callmoldiff
         

         write(*,*) "call thermospheric photochemistry ?",
     &   " (only if callthermos=.true.)"
         thermochem=.false. ! default value
         call getin("thermochem",thermochem)
         write(*,*) " thermochem = ",thermochem

         write(*,*) "date for solar flux calculation:",
     &   " (1985 < date < 2002)"
         write(*,*) "(Solar min=1996.4 ave=1993.4 max=1990.6)"
         solarcondate=1993.4 ! default value
         call getin("solarcondate",solarcondate)
         write(*,*) " solarcondate = ",solarcondate
         

         if (.not.callthermos) then
           if (thermoswater) then
             print*,'if thermoswater is set, callthermos must be true'
             stop
           endif          
           if (callconduct) then
             print*,'if callconduct is set, callthermos must be true'
             stop
           endif        
           if (calleuv) then
             print*,'if calleuv is set, callthermos must be true'
             stop
           endif         
           if (callmolvis) then
             print*,'if callmolvis is set, callthermos must be true'
             stop
           endif        
           if (callmoldiff) then
             print*,'if callmoldiff is set, callthermos must be true'
             stop
           endif          
           if (thermochem) then
             print*,'if thermochem is set, callthermos must be true'
             stop
           endif          
        endif

! Test of incompatibility:
! if photochem is used, then water should be used too

         if (photochem.and..not.water) then
           print*,'if photochem is used, water should be used too'
           stop
         endif

! if callthermos is used, then thermoswater should be used too 
! (if water not used already)

         if (callthermos .and. .not.water) then
           if (callthermos .and. .not.thermoswater) then
             print*,'if callthermos is used, water or thermoswater 
     &               should be used too'
             stop
           endif
         endif

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

      pi=2.*asin(1.)

!     managing the tracers, and tests:
!     -------------------------------

      if(tracer) then

!          when photochem is used, nqchem_min is the rank
!          of the first chemical species

! Ehouarn: nqchem_min is now meaningless and no longer used
!       nqchem_min = 1
       if (photochem .or. callthermos) then
         chem = .true.
!        if (doubleq) then
!          nqchem_min = 3
!        else
!          nqchem_min = dustbin+1
!        end if
       end if

       if (water .or. thermoswater) h2o = .true.

c          TESTS

       print*,'inichim_readcallphys: TRACERS:'

       if ((doubleq).and.(h2o).and.
     $     (chem)) then
!         print*,' 1: dust ; 2: dust (doubleq)'
!         print*,' 3 to ',nqmx-2,': chemistry'
!         print*,nqmx-1,': water ice ; ',nqmx,': water vapor'
         print*,' 2 dust tracers (doubleq)'
         print*,' 1 water vapour tracer'
         print*,' 1 water ice tracer'
         print*,nqmx-4,' chemistry tracers'
       endif

       if ((doubleq).and.(h2o).and.
     $     .not.(chem)) then
!         print*,' 1: dust ; 2: dust (doubleq)'
!         print*,nqmx-1,': water ice ; ',nqmx,': water vapor'
         print*,' 2 dust tracers (doubleq)'
         print*,' 1 water vapour tracer'
         print*,' 1 water ice tracer'
         if (nqmx.ne.4) then
           print*,'nqmx should be 4 with these options.'
		   print*,'(or check callphys.def)'
           stop
         endif
       endif

       if ((doubleq).and..not.(h2o)) then
!         print*,' 1: dust ; 2: dust (doubleq)'
         print*,' 2 dust tracers (doubleq)'
         if (nqmx.ne.2) then
           print*,'nqmx should be 2 with these options...'
		   print*,'(or check callphys.def)'
           stop
         endif
       endif

       if (.not.(doubleq).and.(h2o).and.
     $     (chem)) then
         if (dustbin.gt.0) then
!           print*,' 1 to ',dustbin,': dust bins'
           print*,dustbin,' dust bins'
         endif
!         print*,nqchem_min,' to ',nqmx-2,': chemistry'
!         print*,nqmx-1,': water ice ; ',nqmx,': water vapor'
         print*,nqmx-2-dustbin,' chemistry tracers'
         print*,' 1 water vapour tracer'
         print*,' 1 water ice tracer'
       endif

       if (.not.(doubleq).and.(h2o).and.
     $     .not.(chem)) then
         if (dustbin.gt.0) then
!           print*,' 1 to ',dustbin,': dust bins'
           print*,dustbin,' dust bins'
         endif
!         print*,nqmx-1,': water ice ; ',nqmx,': water vapor'
         print*,' 1 water vapour tracer'
         print*,' 1 water ice tracer'
         if (nqmx.ne.(dustbin+2)) then
           print*,'nqmx should be ',(dustbin+2),
     $            ' with these options...'
		   print*,'(or check callphys.def)'
           stop
         endif
       endif

       if (.not.(doubleq).and..not.(h2o)) then
         if (dustbin.gt.0) then
!           print*,' 1 to ',dustbin,': dust bins'
           print*,dustbin,' dust bins'
           if (nqmx.ne.dustbin) then
             print*,'nqmx should be ',dustbin,
     $              ' with these options...'
             print*,'(or check callphys.def)'
             stop
           endif
         else
           print*,'dustbin=',dustbin,
     $            ': tracer should be F with these options...'
     $           ,'UNLESS you just want to move tracers around '
         endif
       endif

      endif ! of if (tracer)

      RETURN
      END