source: LMDZ5/branches/IPSLCM5A2.1_ISO/libf/phyiso/isotopes_mod.F90 @ 4176

#ifdef ISO
! $Id: $

MODULE isotopes_mod
USE infotrac_phy, ONLY: ntraciso,niso,indnum_fn_num,ok_isotrac,use_iso, &
&       niso_possibles
IMPLICIT NONE
SAVE

! contient toutes les variables isotopiques et leur initialisation
! les routines specifiquement isotopiques sont dans
! isotopes_routines_mod pour éviter dépendance circulaire avec
! isotopes_verif_mod.


! indices des isotopes
integer, save :: iso_eau,iso_HDO,iso_O18,iso_O17,iso_HTO ! indices de 1 à niso: les isos n'existant pas sont mis à 0
!$OMP THREADPRIVATE(iso_eau,iso_HDO,iso_O18,iso_O17,iso_HTO)

integer :: iso_eau_possible,iso_HDO_possible,iso_O18_possible,iso_O17_possible,iso_HTO_possible ! indices de 1 à niso_possibles: ils correspondent aux tableaux définis dans infotrac:
! tnom_iso=(/'eau','HDO','O18','O17','HTO'/)
! ce sont ces indices qui doivent être utilisés avec use_iso, puisque use_iso est défini comme DIMENSION(niso_possibles)
parameter (iso_eau_possible=1)
parameter (iso_HDO_possible=2)
parameter (iso_O18_possible=3)
parameter (iso_O17_possible=4)
parameter (iso_HTO_possible=5)

integer, save :: ntracisoOR
!$OMP THREADPRIVATE(ntracisoOR)

! variables indépendantes des isotopes

real, save :: pxtmelt,pxtice,pxtmin,pxtmax
!$OMP THREADPRIVATE(pxtmelt,pxtice,pxtmin,pxtmax)
real, save ::  tdifexp, tv0cin, thumxt1
!$OMP THREADPRIVATE(tdifexp, tv0cin, thumxt1)
integer, save ::  ntot 
!$OMP THREADPRIVATE(ntot)
real, save ::  h_land_ice
!$OMP THREADPRIVATE(h_land_ice)
real, save ::  P_veg
!$OMP THREADPRIVATE(P_veg)
real, save ::  musi,lambda_sursat
!$OMP THREADPRIVATE(lambda_sursat)
real, save ::  Kd
!$OMP THREADPRIVATE(Kd)
real, save ::  rh_cste_surf_cond,T_cste_surf_cond
!$OMP THREADPRIVATE(rh_cste_surf_cond,T_cste_surf_cond)

logical, save ::   bidouille_anti_divergence
                ! si true, rappel régulier de xteau vers q, pour éviter dérives lentes
!$OMP THREADPRIVATE(bidouille_anti_divergence)
logical, save ::   essai_convergence 
                ! si false, on fait rigoureusement comme dans LMDZ sans isotopes, 
                ! meme si c'est génant pour les isotopes
!$OMP THREADPRIVATE(essai_convergence)
integer, save ::   initialisation_iso
                ! 0: dans fichier
                ! 1: R=0
                ! 2: R selon distill rayleigh
                ! 3: R=Rsmow
!$OMP THREADPRIVATE(initialisation_iso)
integer, save ::  modif_SST ! 0 par defaut, 1 si on veut modifier la sst
                ! 2 et 3: profils de SST
!$OMP THREADPRIVATE(modif_SST)
real, save ::  deltaTtest ! modif de la SST, uniforme.
!$OMP THREADPRIVATE(deltaTtest)
integer, save ::  modif_sic ! on met des trous dans glace de mer
!$OMP THREADPRIVATE(modif_sic)
real, save ::  deltasic ! fraction de trous minimale
!$OMP THREADPRIVATE(deltasic)
real, save ::  deltaTtestpoles
!$OMP THREADPRIVATE(deltaTtestpoles)
real, save ::  sstlatcrit
!$OMP THREADPRIVATE(sstlatcrit)
real, save ::  dsstlatcrit
!$OMP THREADPRIVATE(dsstlatcrit)
real, save ::  deltaO18_oce
!$OMP THREADPRIVATE(deltaO18_oce)
integer, save ::  albedo_prescrit ! 0 par defaut
                        ! 1 si on veut garder albedo constant
!$OMP THREADPRIVATE(albedo_prescrit)
real, save ::  lon_min_albedo,lon_max_albedo
!$OMP THREADPRIVATE(lon_min_albedo,lon_max_albedo)
real, save :: lat_min_albedo,lat_max_albedo
!$OMP THREADPRIVATE(lat_min_albedo,lat_max_albedo)
real, save ::  deltaP_BL,tdifexp_sol
!$OMP THREADPRIVATE(deltaP_BL,tdifexp_sol)
integer, save ::  ruissellement_pluie,alphak_stewart
!$OMP THREADPRIVATE(ruissellement_pluie,alphak_stewart)
integer, save ::  calendrier_guide
!$OMP THREADPRIVATE(calendrier_guide)
integer, save ::  cste_surf_cond
!$OMP THREADPRIVATE(cste_surf_cond)
real, save ::  mixlen
!$OMP THREADPRIVATE(mixlen)
integer, save ::  evap_cont_cste
!$OMP THREADPRIVATE(evap_cont_cste)
real, save ::  deltaO18_evap_cont,d_evap_cont
!$OMP THREADPRIVATE(deltaO18_evap_cont,d_evap_cont)
integer, save ::  nudge_qsol,region_nudge_qsol
!$OMP THREADPRIVATE(nudge_qsol,region_nudge_qsol)
integer, save :: nlevmaxO17
!$OMP THREADPRIVATE(nlevmaxO17)
integer, save ::  no_pce
!       real, save :: slope_limiterxy,slope_limiterz
!$OMP THREADPRIVATE(no_pce)
real, save ::  A_satlim
!$OMP THREADPRIVATE(A_satlim)
integer, save ::  ok_restrict_A_satlim,modif_ratqs
!$OMP THREADPRIVATE(ok_restrict_A_satlim,modif_ratqs)
real, save ::  Pcrit_ratqs,ratqsbasnew
!$OMP THREADPRIVATE(Pcrit_ratqs,ratqsbasnew)
real, save ::  fac_modif_evaoce
!$OMP THREADPRIVATE(fac_modif_evaoce)
integer, save ::  ok_bidouille_wake
!$OMP THREADPRIVATE(ok_bidouille_wake)
logical ::  cond_temp_env
!$OMP THREADPRIVATE(cond_temp_env)


! variables tableaux fn de niso
real, ALLOCATABLE, DIMENSION(:), save :: tnat, toce, tcorr
!$OMP THREADPRIVATE(tnat, toce, tcorr)
real, ALLOCATABLE, DIMENSION(:), save :: tdifrel
!$OMP THREADPRIVATE(tdifrel)
real, ALLOCATABLE, DIMENSION(:), save :: talph1, talph2, talph3
!$OMP THREADPRIVATE(talph1, talph2, talph3)
real, ALLOCATABLE, DIMENSION(:), save :: talps1, talps2
!$OMP THREADPRIVATE(talps1, talps2)
real, ALLOCATABLE, DIMENSION(:), save :: tkcin0, tkcin1, tkcin2
!$OMP THREADPRIVATE(tkcin0, tkcin1, tkcin2)
real, ALLOCATABLE, DIMENSION(:), save :: alpha_liq_sol
!$OMP THREADPRIVATE(alpha_liq_sol)
real, ALLOCATABLE, DIMENSION(:), save :: Rdefault, Rmethox
!$OMP THREADPRIVATE(Rdefault, Rmethox)
character*3, ALLOCATABLE, DIMENSION(:), save :: striso
!$OMP THREADPRIVATE(striso)
real, save :: fac_coeff_eq17_liq, fac_coeff_eq17_ice
!$OMP THREADPRIVATE(fac_coeff_eq17_liq, fac_coeff_eq17_ice)

      real ridicule ! valeur maximale pour qu'une variable de type
                    ! rapoport de mélange puisse être considérée comme négligeable. Si
                    ! négligeable, alors on ne verifie pas si sa compo iso esta bérrante.
      parameter (ridicule=1e-12)      
!      parameter (ridicule=1)
!
      real ridicule_rain ! valeur limite de ridicule pour les flux de pluies (rain, zrfl...)
      parameter (ridicule_rain=1e-8) ! en kg/s <-> 1e-3mm/day

      real ridicule_evap ! valeur limite de ridicule pour les evap
      parameter (ridicule_evap=ridicule_rain*1e-2) ! en kg/s <-> 1e-3mm/day

      real ridicule_qsol ! valeur limite de ridicule pour les qsol
      parameter (ridicule_qsol=ridicule_rain) ! en kg <-> 1e-8kg

      real ridicule_snow ! valeur limite de ridicule pour les snow
      parameter (ridicule_snow=ridicule_qsol) ! en kg/s <-> 1e-8kg
      
        real expb_max
        parameter (expb_max=30.0)

        ! spécifique au tritium:
        

logical, save :: ok_prod_nucl_tritium ! si oui, production de tritium par essais nucleaires
!$OMP THREADPRIVATE(ok_prod_nucl_tritium)
        integer nessai
        parameter (nessai=486)
        integer, save :: day_nucl(nessai)
!$OMP THREADPRIVATE(day_nucl)
        integer, save :: month_nucl(nessai)
!$OMP THREADPRIVATE(month_nucl)
        integer, save :: year_nucl(nessai)
!$OMP THREADPRIVATE(year_nucl)
        real, save :: lat_nucl(nessai)
!$OMP THREADPRIVATE(lat_nucl)
        real, save :: lon_nucl(nessai)
!$OMP THREADPRIVATE(lon_nucl)
        real, save :: zmin_nucl(nessai)
!$OMP THREADPRIVATE(zmin_nucl)
        real, save :: zmax_nucl(nessai)
!$OMP THREADPRIVATE(zmax_nucl)
        real, save :: HTO_nucl(nessai)
!$OMP THREADPRIVATE(HTO_nucl)

 
CONTAINS

  SUBROUTINE iso_init()
      use IOIPSL ! getin
      implicit none

! -- local variables:

      integer ixt
      ! référence O18
      real fac_enrichoce18
      real alpha_liq_sol_O18, &
     &     talph1_O18,talph2_O18,talph3_O18, &
     &     talps1_O18,talps2_O18, &
     &     tkcin0_O18,tkcin1_O18,tkcin2_O18, &
     &     tdifrel_O18  
      
      ! cas de l'O17
      real fac_kcin
      real pente_MWL
      integer ierr
      
      logical ok_nocinsat, ok_nocinsfc !sensi test
      parameter (ok_nocinsfc=.FALSE.)  ! if T: no kinetic effect in sfc evap
      parameter (ok_nocinsat=.FALSE.)  ! if T: no sursaturation effect for ice
      logical Rdefault_smow
      parameter (Rdefault_smow=.FALSE.) ! si T: Rdefault=smow; si F: nul
      ! pour le tritium
      integer iessai

! allocations mémoire

allocate (tnat(niso))
allocate (toce(niso))
allocate (tcorr(niso))
allocate (tdifrel(niso))
allocate (talph1(niso))
allocate (talph2(niso))
allocate (talph3(niso))
allocate (talps1(niso))
allocate (talps2(niso))
allocate (tkcin0(niso))
allocate (tkcin1(niso))
allocate (tkcin2(niso))
allocate (alpha_liq_sol(niso))
allocate (Rdefault(niso))
allocate (Rmethox(niso))
allocate (striso(niso))


!--------------------------------------------------------------
! General:
!--------------------------------------------------------------

! -- verif du nombre d'isotopes:
      write(*,*) 'iso_init 64: niso=',niso

! init de ntracisoOR: on écrasera en cas de ok_isotrac si complications avec
! ORCHIDEE
      ntracisoOR=ntraciso  
              
! -- Type of water isotopes:

        iso_eau=indnum_fn_num(1)
        iso_HDO=indnum_fn_num(2)
        iso_O18=indnum_fn_num(3)
        iso_O17=indnum_fn_num(4)
        iso_HTO=indnum_fn_num(5)
        write(*,*) 'iso_init 59: iso_eau=',iso_eau
        write(*,*) 'iso_HDO=',iso_HDO
        write(*,*) 'iso_O18=',iso_O18
        write(*,*) 'iso_O17=',iso_O17
        write(*,*) 'iso_HTO=',iso_HTO
        write(*,*) 'iso_init 251: use_iso=',use_iso

      ! initialisation
        lambda_sursat=0.004
        thumxt1=0.75*1.2
        ntot=20
        h_land_ice=20. ! à comparer aux 3000mm de snow_max
        P_veg=1.0
        bidouille_anti_divergence=.false.
        essai_convergence=.false.
        initialisation_iso=0
        modif_sst=0
        modif_sic=0
        deltaTtest=0.0
        deltasic=0.1
        deltaTtestpoles=0.0
        sstlatcrit=30.0
        deltaO18_oce=0.0
        albedo_prescrit=0
        lon_min_albedo=-200
        lon_max_albedo=200
        lat_min_albedo=-100
        lat_max_albedo=100
        deltaP_BL=10.0
        ruissellement_pluie=0
        alphak_stewart=1
        tdifexp_sol=0.67
        calendrier_guide=0
        cste_surf_cond=0
mixlen=35.0        
evap_cont_cste=0.0
deltaO18_evap_cont=0.0
d_evap_cont=0.0
nudge_qsol=0
region_nudge_qsol=1
nlevmaxO17=50
no_pce=0
A_satlim=1.0
ok_restrict_A_satlim=0
!        slope_limiterxy=2.0
!        slope_limiterz=2.0
modif_ratqs=0
Pcrit_ratqs=500.0
ratqsbasnew=0.05

fac_modif_evaoce=1.0
ok_bidouille_wake=0
cond_temp_env=.false.
! si oui, la temperature de cond est celle de l'environnement,
! pour eviter bugs quand temperature dans ascendances convs est
! mal calculee
ok_prod_nucl_tritium=.false.

! lecture des paramètres isotopiques:
call getin('lambda',lambda_sursat)
call getin('thumxt1',thumxt1)
call getin('ntot',ntot)
call getin('h_land_ice',h_land_ice)
call getin('P_veg',P_veg)
call getin('bidouille_anti_divergence',bidouille_anti_divergence)
call getin('essai_convergence',essai_convergence)
call getin('initialisation_iso',initialisation_iso)
!if (ok_isotrac) then      
!if (initialisation_iso.eq.0) then
!  call getin('initialisation_isotrac',initialisation_isotrac)
!endif !if (initialisation_iso.eq.0) then
!endif !if (ok_isotrac) then      
call getin('modif_sst',modif_sst)
if (modif_sst.ge.1) then
call getin('deltaTtest',deltaTtest)
if (modif_sst.ge.2) then
  call getin('deltaTtestpoles',deltaTtestpoles)
  call getin('sstlatcrit',sstlatcrit)
#ifdef ISOVERIF
      !call iso_verif_positif(sstlatcrit,'iso_init 107')
      if (sstlatcrit.lt.0.0) then
        write(*,*) 'iso_init 270: sstlatcrit=',sstlatcrit
        stop
      endif
#endif              
  if (modif_sst.ge.3) then  
      call getin('dsstlatcrit',dsstlatcrit)
#ifdef ISOVERIF
      !call iso_verif_positif(dsstlatcrit,'iso_init 110')
      if (sstlatcrit.lt.0.0) then
        write(*,*) 'iso_init 279: dsstlatcrit=',dsstlatcrit
        stop
      endif
#endif              
  endif !if (modif_sst.ge.3) then
endif !if (modif_sst.ge.2) then
endif !  if (modif_sst.ge.1) then   
call getin('modif_sic',modif_sic)
if (modif_sic.ge.1) then
call getin('deltasic',deltasic)
endif !if (modif_sic.ge.1) then

call getin('albedo_prescrit',albedo_prescrit)
call getin('lon_min_albedo',lon_min_albedo)
call getin('lon_max_albedo',lon_max_albedo)
call getin('lat_min_albedo',lat_min_albedo)
call getin('lat_max_albedo',lat_max_albedo)
call getin('deltaO18_oce',deltaO18_oce)
call getin('deltaP_BL',deltaP_BL)
call getin('ruissellement_pluie',ruissellement_pluie)
call getin('alphak_stewart',alphak_stewart)
call getin('tdifexp_sol',tdifexp_sol)
call getin('calendrier_guide',calendrier_guide)
call getin('cste_surf_cond',cste_surf_cond)
call getin('mixlen',mixlen)
call getin('evap_cont_cste',evap_cont_cste)
call getin('deltaO18_evap_cont',deltaO18_evap_cont)
call getin('d_evap_cont',d_evap_cont)  
call getin('nudge_qsol',nudge_qsol)
call getin('region_nudge_qsol',region_nudge_qsol)
call getin('no_pce',no_pce)
call getin('A_satlim',A_satlim)
call getin('ok_restrict_A_satlim',ok_restrict_A_satlim)
#ifdef ISOVERIF      
!call iso_verif_positif(1.0-A_satlim,'iso_init 158')
      if (A_satlim.gt.1.0) then
        write(*,*) 'iso_init 315: A_satlim=',A_satlim
        stop
      endif
#endif          
!      call getin('slope_limiterxy',slope_limiterxy) 
!      call getin('slope_limiterz',slope_limiterz) 
call getin('modif_ratqs',modif_ratqs)
call getin('Pcrit_ratqs',Pcrit_ratqs)
call getin('ratqsbasnew',ratqsbasnew)
call getin('fac_modif_evaoce',fac_modif_evaoce)
call getin('ok_bidouille_wake',ok_bidouille_wake)
call getin('cond_temp_env',cond_temp_env)
if (use_iso(iso_HTO_possible)) then
  ok_prod_nucl_tritium=.true.
  call getin('ok_prod_nucl_tritium',ok_prod_nucl_tritium)
endif

write(*,*) 'lambda,thumxt1=',lambda_sursat,thumxt1
write(*,*) 'bidouille_anti_divergence=',bidouille_anti_divergence
write(*,*) 'essai_convergence=',essai_convergence
write(*,*) 'initialisation_iso=',initialisation_iso
write(*,*) 'modif_sst=',modif_sst
if (modif_sst.ge.1) then
write(*,*) 'deltaTtest=',deltaTtest
if (modif_sst.ge.2) then  
write(*,*) 'deltaTtestpoles,sstlatcrit=', &
&           deltaTtestpoles,sstlatcrit
if (modif_sst.ge.3) then   
 write(*,*) 'dsstlatcrit=',dsstlatcrit
endif !if (modif_sst.ge.3) then
endif !if (modif_sst.ge.2) then
endif !if (modif_sst.ge.1) then
write(*,*) 'modif_sic=',modif_sic
if (modif_sic.ge.1) then  
write(*,*) 'deltasic=',deltasic
endif !if (modif_sic.ge.1) then
write(*,*) 'deltaO18_oce=',deltaO18_oce
write(*,*) 'albedo_prescrit=',albedo_prescrit
if (albedo_prescrit.eq.1) then
 write(*,*) 'lon_min_albedo,lon_max_albedo=', &
&           lon_min_albedo,lon_max_albedo
 write(*,*) 'lat_min_albedo,lat_max_albedo=', &
&           lat_min_albedo,lat_max_albedo
endif !if (albedo_prescrit.eq.1) then
write(*,*) 'deltaP_BL,ruissellement_pluie,alphak_stewart=', &
&       deltaP_BL,ruissellement_pluie,alphak_stewart
write(*,*) 'cste_surf_cond=',cste_surf_cond
write(*,*) 'mixlen=',mixlen
write(*,*) 'tdifexp_sol=',tdifexp_sol
write(*,*) 'calendrier_guide=',calendrier_guide
write(*,*) 'evap_cont_cste=',evap_cont_cste
write(*,*) 'deltaO18_evap_cont,d_evap_cont=', &
&           deltaO18_evap_cont,d_evap_cont
write(*,*) 'nudge_qsol,region_nudge_qsol=', &
&  nudge_qsol,region_nudge_qsol  
write(*,*) 'nlevmaxO17=',nlevmaxO17
write(*,*) 'no_pce=',no_pce
write(*,*) 'A_satlim=',A_satlim
write(*,*) 'ok_restrict_A_satlim=',ok_restrict_A_satlim
!      write(*,*) 'slope_limiterxy=',slope_limiterxy
!      write(*,*) 'slope_limiterz=',slope_limiterz
write(*,*) 'modif_ratqs=',modif_ratqs
write(*,*) 'Pcrit_ratqs=',Pcrit_ratqs
write(*,*) 'ratqsbasnew=',ratqsbasnew
write(*,*) 'fac_modif_evaoce=',fac_modif_evaoce
write(*,*) 'ok_bidouille_wake=',ok_bidouille_wake
write(*,*) 'cond_temp_env=',cond_temp_env
write(*,*) 'ok_prod_nucl_tritium=',ok_prod_nucl_tritium
          

!--------------------------------------------------------------
! Parameters that do not depend on the nature of water isotopes:
!--------------------------------------------------------------

! -- temperature at which ice condensate starts to form (valeur ECHAM?):
pxtmelt=273.15
!      pxtmelt=273.15-10.0 ! test PHASE

! -- temperature at which all condensate is ice:
pxtice=273.15-10.0
!      pxtice=273.15-30.0 ! test PHASE

! -- minimum temperature to calculate fractionation coeff
pxtmin=273.15-120.0 ! On ne calcule qu'au dessus de -120°C
pxtmax=273.15+60.0 ! On ne calcule qu'au dessus de +60°C
! remarque: les coeffs ont été mesurés seulement jusq'à -40!

! -- a constant for alpha_eff for equilibrium below cloud base:
tdifexp=0.58
tv0cin=7.0

! facteurs lambda et mu dans Si=musi-lambda*T
musi=1.0
if (ok_nocinsat) then
lambda_sursat = 0.0 ! no sursaturation effect
endif            


! diffusion dans le sol
Kd=2.5e-9 ! m2/s   

! cas où cste_surf_cond: on met rhs ou/et Ts cste pour voir
rh_cste_surf_cond=0.6
T_cste_surf_cond=288.0

!--------------------------------------------------------------
! Parameters that depend on the nature of water isotopes:
!--------------------------------------------------------------
! ** constantes locales
fac_enrichoce18=0.0005
! on a alors tcor018=1+fac_enrichoce18
! tcorD=1+fac_enrichoce18*8
! tcorO17=1+fac_enrichoce18*0.528
alpha_liq_sol_O18=1.00291 ! valeur de Lehmann & Siegenthaler, 1991, 
  ! Journal of Glaciology, vol 37, p 23
talph1_O18=1137.
talph2_O18=-0.4156
talph3_O18=-2.0667E-3
talps1_O18=11.839
talps2_O18=-0.028244
tkcin0_O18 = 0.006
tkcin1_O18 = 0.000285
tkcin2_O18 = 0.00082
tdifrel_O18= 1./0.9723

! rapport des ln(alphaeq) entre O18 et O17
fac_coeff_eq17_liq=0.529 ! donné par Amaelle
!      fac_coeff_eq17_ice=0.528 ! slope MWL 
fac_coeff_eq17_ice=0.529


write(*,*) 'iso_O18,iso_HDO,iso_eau=',iso_O18,iso_HDO,iso_eau
do 999 ixt = 1, niso
write(*,*) 'iso_init 80: ixt=',ixt


! -- kinetic factor for surface evaporation:
! (cf: kcin = tkcin0                  if |V|<tv0cin
!      kcin = tkcin1*|Vsurf| + tkcin2 if |V|>tv0cin )
! (Rq: formula discontinuous for |V|=tv0cin... )        

! -- main:
if (ixt.eq.iso_HTO) then ! Tritium
  tkcin0(ixt) = 0.01056
  tkcin1(ixt) = 0.0005016
  tkcin2(ixt) = 0.0014432
  tnat(ixt)=0.
  !toce(ixt)=2.2222E-8 ! corrigé par Alex Cauquoin
  !toce(ixt)=1.0E-18 ! rapport 3H/1H ocean
  toce(ixt)=4.0E-19 ! rapport T/H = 0.2 TU Dreisigacker and Roether 1978
  tcorr(ixt)=1.
  tdifrel(ixt)=1./0.968
  talph1(ixt)=46480.
  talph2(ixt)=-103.87
  talph3(ixt)=0.
  talps1(ixt)=46480.
  talps2(ixt)=-103.87
  alpha_liq_sol(ixt)=1. 
  Rdefault(ixt)=0.0
  Rmethox(ixt)=0.0
  striso(ixt)='HTO'
endif     
if (ixt.eq.iso_O17) then ! Deuterium
  pente_MWL=0.528
!          tdifrel(ixt)=1./0.985452 ! donné par Amaelle 
  tdifrel(ixt)=1./0.98555 ! valeur utilisée en 1D et dans modèle de LdG
!          fac_kcin=0.5145 ! donné par Amaelle
  fac_kcin= (tdifrel(ixt)-1.0)/(tdifrel_O18-1.0)
  tkcin0(ixt) = tkcin0_O18*fac_kcin
  tkcin1(ixt) = tkcin1_O18*fac_kcin
  tkcin2(ixt) = tkcin2_O18*fac_kcin
  tnat(ixt)=0.004/100. ! O17 représente 0.004% de l'oxygène
  toce(ixt)=tnat(ixt)*(1.0+deltaO18_oce/1000.0)**pente_MWL
  tcorr(ixt)=1.0+fac_enrichoce18*pente_MWL ! donné par Amaelle            
  talph1(ixt)=talph1_O18
  talph2(ixt)=talph2_O18
  talph3(ixt)=talph3_O18
  talps1(ixt)=talps1_O18
  talps2(ixt)=talps2_O18     
  alpha_liq_sol(ixt)=(alpha_liq_sol_O18)**fac_coeff_eq17_liq 
  if (Rdefault_smow) then    
        Rdefault(ixt)=tnat(ixt)*(-3.15/1000.0+1.0)
  else
        Rdefault(ixt)=0.0
  endif
  Rmethox(ixt)=(230./1000.+1.)*tnat(ixt) !Zahn et al 2006
  striso(ixt)='O17'
endif

if (ixt.eq.iso_O18) then ! Oxygene18
  tkcin0(ixt) = tkcin0_O18
  tkcin1(ixt) = tkcin1_O18
  tkcin2(ixt) = tkcin2_O18
  tnat(ixt)=2005.2E-6
  toce(ixt)=tnat(ixt)*(1.0+deltaO18_oce/1000.0)
  tcorr(ixt)=1.0+fac_enrichoce18
  tdifrel(ixt)=tdifrel_O18
  talph1(ixt)=talph1_O18
  talph2(ixt)=talph2_O18
  talph3(ixt)=talph3_O18
  talps1(ixt)=talps1_O18
  talps2(ixt)=talps2_O18
  alpha_liq_sol(ixt)=alpha_liq_sol_O18   
  if (Rdefault_smow) then    
        Rdefault(ixt)=tnat(ixt)*(-6.0/1000.0+1.0)
  else
        Rdefault(ixt)=0.0
  endif
  Rmethox(ixt)=(130./1000.+1.)*tnat(ixt) !Zahn et al 2006   
!       write(*,*) 'iso_init 163: ZXalpha_liq_sol=',ZXalpha_liq_sol
  striso(ixt)='O18'
  write(*,*) 'isotopes_mod 519: ixt,striso(ixt)=',ixt,striso(ixt)
endif

if (ixt.eq.iso_HDO) then ! Deuterium
  pente_MWL=8.0
!          fac_kcin=0.88
  tdifrel(ixt)=1./0.9755
  fac_kcin= (tdifrel(ixt)-1)/(tdifrel_O18-1)
  tkcin0(ixt) = tkcin0_O18*fac_kcin
  tkcin1(ixt) = tkcin1_O18*fac_kcin
  tkcin2(ixt) = tkcin2_O18*fac_kcin
  tnat(ixt)=155.76E-6
  toce(ixt)=tnat(ixt)*(1.0+pente_MWL*deltaO18_oce/1000.0)
  tcorr(ixt)=1.0+fac_enrichoce18*pente_MWL          
  talph1(ixt)=24844.
  talph2(ixt)=-76.248
  talph3(ixt)=52.612E-3
  talps1(ixt)=16288.
  talps2(ixt)=-0.0934
  !ZXalpha_liq_sol=1.0192 ! Weston, Ralph, 1955
  alpha_liq_sol(ixt)=1.0212 
  ! valeur de Lehmann & Siegenthaler, 1991, Journal of
  ! Glaciology, vol 37, p 23
  if (Rdefault_smow) then    
    Rdefault(ixt)=tnat(ixt)*((-6.0*pente_MWL+10.0)/1000.0+1.0)
  else
    Rdefault(ixt)=0.0
  endif
Rmethox(ixt)=tnat(ixt)*(-25.0/1000.+1.) ! Zahn et al 2006
striso(ixt)='HDO'
  write(*,*) 'isotopes_mod 548: ixt,striso(ixt)=',ixt,striso(ixt)
endif

!       write(*,*) 'iso_init 163: ZXalpha_liq_sol=',ZXalpha_liq_sol
if (ixt.eq.iso_eau) then ! Oxygene16
  tkcin0(ixt) = 0.0
  tkcin1(ixt) = 0.0
  tkcin2(ixt) = 0.0
  tnat(ixt)=1.
  toce(ixt)=tnat(ixt)
  tcorr(ixt)=1.0
  tdifrel(ixt)=1.
  talph1(ixt)=0.
  talph2(ixt)=0.
  talph3(ixt)=0.
  talps1(ixt)=0.
  talph3(ixt)=0.
  alpha_liq_sol(ixt)=1.
  if (Rdefault_smow) then
        Rdefault(ixt)=tnat(ixt)*1.0
  else
        Rdefault(ixt)=1.0
  endif
  Rmethox(ixt)=1.0
  striso(ixt)='eau'
endif

999   continue

! test de sensibilité: 
if (ok_nocinsfc) then ! no kinetic effect in sfc evaporation
 do ixt=1,niso
  tkcin0(ixt) = 0.0 
  tkcin1(ixt) = 0.0 
  tkcin2(ixt) = 0.0
 enddo 
endif

! fermeture fichier de paramètres
close(unit=32)

! nom des isotopes

! verif
write(*,*) 'iso_init 285: verif initialisation:'

do ixt=1,niso
  write(*,*) '* striso(',ixt,')=<'//striso(ixt)//'>'
  write(*,*) 'tnat(',ixt,')=',tnat(ixt)
!          write(*,*) 'alpha_liq_sol(',ixt,')=',alpha_liq_sol(ixt)
!          write(*,*) 'tkcin0(',ixt,')=',tkcin0(ixt)
!          write(*,*) 'tdifrel(',ixt,')=',tdifrel(ixt)
enddo
write(*,*) 'iso_init 69: lambda=',lambda_sursat
write(*,*) 'iso_init 69: thumxt1=',thumxt1
write(*,*) 'iso_init 69: h_land_ice=',h_land_ice
write(*,*) 'iso_init 69: P_veg=',P_veg   
            

END SUBROUTINE iso_init

! functions basiques mises ici pour éviter dépendances circulaires
        !***********
        function double_to_real(a)
        implicit none
        double precision a
        real double_to_real

        double_to_real=a

        return
        end function double_to_real

        !***********
        function real_to_double(a)
        implicit none
        real a
        double precision real_to_double

        real_to_double=a

        return
        end function real_to_double

END MODULE isotopes_mod
#endif