Ignore:
Timestamp:
May 9, 2022, 12:35:40 PM (2 years ago)
Author:
dcugnet
Message:
  • Some variables are renamed or replaced by direct equivalents:
    • iso_indnum -> tracers(:)%iso_iName
    • niso_possibles -> niso
    • iqiso -> iqIsoPha ; index_trac -> itZonIso
    • ok_iso_verif -> isoCheck
    • ntraceurs_zone -> nzone ; ntraciso -> ntiso
    • qperemin -> min_qparent ; masseqmin -> min_qmass ; ratiomin -> min_ratio
  • Some renamed variables are only aliased with the older name (using USE <module>, ONLY: <oldName> => <newName>) in routines where they are repeated many times.
  • Few hard-coded indexes are now computed (examples: ilic, iso, ivap, irneb, iq_vap, iq_liq, iso_H2O, iso_HDO, iso_HTO, iso_O17, iso_O18).
  • The IF(isoCheck) test is now embedded in the check_isotopes_seq and check_isotopes_loc routines (lighter calling).
File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmdiso/isotopes_mod.F90

    r4124 r4143  
    33
    44MODULE isotopes_mod
    5 USE infotrac_phy, ONLY: ntraciso,niso,indnum_fn_num,use_iso, &
    6 &       niso_possibles
    7 IMPLICIT NONE
    8 SAVE
    9 
    10 ! contient toutes les variables isotopiques et leur initialisation
    11 ! les routines specifiquement isotopiques sont dans
    12 ! isotopes_routines_mod pour éviter dépendance circulaire avec
    13 ! isotopes_verif_mod.
    14 
    15 
    16 ! indices des isotopes
    17 integer, save :: iso_eau,iso_HDO,iso_O18,iso_O17,iso_HTO ! indices de 1 à niso: les isos n'existant pas sont mis à 0
    18 !$OMP THREADPRIVATE(iso_eau,iso_HDO,iso_O18,iso_O17,iso_HTO)
    19 
    20 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:
    21 ! tnom_iso=(/'eau','HDO','O18','O17','HTO'/)
    22 ! ce sont ces indices qui doivent être utilisés avec use_iso, puisque use_iso est défini comme DIMENSION(niso_possibles)
    23 parameter (iso_eau_possible=1)
    24 parameter (iso_HDO_possible=2)
    25 parameter (iso_O18_possible=3)
    26 parameter (iso_O17_possible=4)
    27 parameter (iso_HTO_possible=5)
    28 
    29 integer, save :: ntracisoOR
     5   USE strings_mod, ONLY: msg, real2str, int2str, bool2str, maxlen, strIdx, strStack
     6   IMPLICIT NONE
     7   INTERFACE get_in; MODULE PROCEDURE getinp_s, getinp_i, getinp_r, getinp_l;  END INTERFACE get_in
     8   SAVE
     9
     10  !--- Contains all isotopic variables + their initialization
     11  !--- Isotopes-specific routines are in isotopes_routines_mod to avoid circular dependencies with isotopes_verif_mod.
     12
     13   !--- Isotopes indices (in [1,niso] ; non-existing => 0 index)
     14   INTEGER, SAVE :: iso_eau, iso_HDO, iso_O18, iso_O17, iso_HTO
     15!$OMP THREADPRIVATE(iso_eau, iso_HDO, iso_O18, iso_O17, iso_HTO)
     16
     17   INTEGER, SAVE :: ntracisoOR
    3018!$OMP THREADPRIVATE(ntracisoOR)
    3119
    32 ! variables indépendantes des isotopes
    33 
    34 real, save :: pxtmelt,pxtice,pxtmin,pxtmax
    35 !$OMP THREADPRIVATE(pxtmelt,pxtice,pxtmin,pxtmax)
    36 real, save ::  tdifexp, tv0cin, thumxt1
     20   !--- Variables not depending on isotopes
     21   REAL,    SAVE :: pxtmelt, pxtice, pxtmin, pxtmax
     22!$OMP THREADPRIVATE(pxtmelt, pxtice, pxtmin, pxtmax)
     23   REAL,    SAVE :: tdifexp, tv0cin, thumxt1
    3724!$OMP THREADPRIVATE(tdifexp, tv0cin, thumxt1)
    38 integer, save :: ntot
     25   INTEGER, SAVE :: ntot
    3926!$OMP THREADPRIVATE(ntot)
    40 real, save :: h_land_ice
     27   REAL,    SAVE :: h_land_ice
    4128!$OMP THREADPRIVATE(h_land_ice)
    42 real, save :: P_veg
     29   REAL,    SAVE :: P_veg
    4330!$OMP THREADPRIVATE(P_veg)
    44 real, save ::  musi,lambda_sursat
    45 !$OMP THREADPRIVATE(lambda_sursat)
    46 real, save :: Kd
     31   REAL,    SAVE :: musi, lambda_sursat
     32!$OMP THREADPRIVATE(musi, lambda_sursat)
     33   REAL,    SAVE :: Kd
    4734!$OMP THREADPRIVATE(Kd)
    48 real, save ::  rh_cste_surf_cond,T_cste_surf_cond
    49 !$OMP THREADPRIVATE(rh_cste_surf_cond,T_cste_surf_cond)
    50 
    51 logical, save ::   bidouille_anti_divergence
    52                 ! si true, rappel régulier de xteau vers q, pour éviter dérives lentes
     35   REAL,    SAVE :: rh_cste_surf_cond, T_cste_surf_cond
     36!$OMP THREADPRIVATE(rh_cste_surf_cond, T_cste_surf_cond)
     37   LOGICAL, SAVE :: bidouille_anti_divergence    ! T: regularly, xteau <- q to avoid slow drifts
    5338!$OMP THREADPRIVATE(bidouille_anti_divergence)
    54 logical, save ::   essai_convergence
    55                 ! si false, on fait rigoureusement comme dans LMDZ sans isotopes,
    56                 ! meme si c'est génant pour les isotopes
     39   LOGICAL, SAVE :: essai_convergence            ! F: as in LMDZ without isotopes (bad for isotopes)
    5740!$OMP THREADPRIVATE(essai_convergence)
    58 integer, save ::   initialisation_iso
    59                 ! 0: dans fichier
    60                 ! 1: R=0
    61                 ! 2: R selon distill rayleigh
    62                 ! 3: R=Rsmow
     41   INTEGER, SAVE :: initialisation_iso           ! 0: file ; 1: R=0 ; 2: R=distill. Rayleigh ; 3: R=Rsmow
    6342!$OMP THREADPRIVATE(initialisation_iso)
    64 integer, save ::  modif_SST ! 0 par defaut, 1 si on veut modifier la sst
    65                 ! 2 et 3: profils de SST
     43   INTEGER, SAVE :: modif_SST                    ! 0: default ; 1: modified SST ; 2, 3: SST profiles
    6644!$OMP THREADPRIVATE(modif_SST)
    67 real, save ::  deltaTtest ! modif de la SST, uniforme.
     45   REAL,    SAVE :: deltaTtest                   ! Uniform modification of the SST
    6846!$OMP THREADPRIVATE(deltaTtest)
    69 integer, save ::  modif_sic ! on met des trous dans glace de mer
     47   INTEGER, SAVE :: modif_sic                    ! Holes in the Sea Ice
    7048!$OMP THREADPRIVATE(modif_sic)
    71 real, save ::  deltasic ! fraction de trous minimale
     49   REAL,    SAVE :: deltasic                     ! Minimal holes fraction
    7250!$OMP THREADPRIVATE(deltasic)
    73 real, save :: deltaTtestpoles
     51   REAL,    SAVE :: deltaTtestpoles
    7452!$OMP THREADPRIVATE(deltaTtestpoles)
    75 real, save ::  sstlatcrit
    76 !$OMP THREADPRIVATE(sstlatcrit)
    77 real, save ::  dsstlatcrit
    78 !$OMP THREADPRIVATE(dsstlatcrit)
    79 real, save ::  deltaO18_oce
     53   REAL,    SAVE :: sstlatcrit, dsstlatcrit
     54!$OMP THREADPRIVATE(sstlatcrit, dsstlatcrit)
     55   REAL,    SAVE :: deltaO18_oce
    8056!$OMP THREADPRIVATE(deltaO18_oce)
    81 integer, save ::  albedo_prescrit ! 0 par defaut
    82                         ! 1 si on veut garder albedo constant
     57   INTEGER, SAVE :: albedo_prescrit              ! 0: default ; 1: constant albedo
    8358!$OMP THREADPRIVATE(albedo_prescrit)
    84 real, save ::  lon_min_albedo,lon_max_albedo
    85 !$OMP THREADPRIVATE(lon_min_albedo,lon_max_albedo)
    86 real, save :: lat_min_albedo,lat_max_albedo
    87 !$OMP THREADPRIVATE(lat_min_albedo,lat_max_albedo)
    88 real, save ::  deltaP_BL,tdifexp_sol
     59   REAL,    SAVE :: lon_min_albedo, lon_max_albedo, lat_min_albedo, lat_max_albedo
     60!$OMP THREADPRIVATE(lon_min_albedo, lon_max_albedo, lat_min_albedo, lat_max_albedo)
     61   REAL,    SAVE :: deltaP_BL,tdifexp_sol
    8962!$OMP THREADPRIVATE(deltaP_BL,tdifexp_sol)
    90 integer, save ::  ruissellement_pluie,alphak_stewart
    91 !$OMP THREADPRIVATE(ruissellement_pluie,alphak_stewart)
    92 integer, save :: calendrier_guide
     63   INTEGER, SAVE :: ruissellement_pluie, alphak_stewart
     64!$OMP THREADPRIVATE(ruissellement_pluie, alphak_stewart)
     65   INTEGER, SAVE :: calendrier_guide
    9366!$OMP THREADPRIVATE(calendrier_guide)
    94 integer, save :: cste_surf_cond
     67   INTEGER, SAVE :: cste_surf_cond
    9568!$OMP THREADPRIVATE(cste_surf_cond)
    96 real, save :: mixlen
     69   REAL,    SAVE :: mixlen
    9770!$OMP THREADPRIVATE(mixlen)
    98 integer, save :: evap_cont_cste
     71   INTEGER, SAVE :: evap_cont_cste
    9972!$OMP THREADPRIVATE(evap_cont_cste)
    100 real, save ::  deltaO18_evap_cont,d_evap_cont
    101 !$OMP THREADPRIVATE(deltaO18_evap_cont,d_evap_cont)
    102 integer, save ::  nudge_qsol,region_nudge_qsol
    103 !$OMP THREADPRIVATE(nudge_qsol,region_nudge_qsol)
    104 integer, save :: nlevmaxO17
     73   REAL,    SAVE :: deltaO18_evap_cont, d_evap_cont
     74!$OMP THREADPRIVATE(deltaO18_evap_cont, d_evap_cont)
     75   INTEGER, SAVE :: nudge_qsol, region_nudge_qsol
     76!$OMP THREADPRIVATE(nudge_qsol, region_nudge_qsol)
     77   INTEGER, SAVE :: nlevmaxO17
    10578!$OMP THREADPRIVATE(nlevmaxO17)
    106 integer, save ::  no_pce
    107 !       real, save :: slope_limiterxy,slope_limiterz
     79   INTEGER, SAVE :: no_pce
    10880!$OMP THREADPRIVATE(no_pce)
    109 real, save :: A_satlim
     81   REAL,    SAVE :: A_satlim
    11082!$OMP THREADPRIVATE(A_satlim)
    111 integer, save ::  ok_restrict_A_satlim,modif_ratqs
    112 !$OMP THREADPRIVATE(ok_restrict_A_satlim,modif_ratqs)
    113 real, save ::  Pcrit_ratqs,ratqsbasnew
    114 !$OMP THREADPRIVATE(Pcrit_ratqs,ratqsbasnew)
    115 real, save :: fac_modif_evaoce
     83   INTEGER, SAVE :: ok_restrict_A_satlim, modif_ratqs
     84!$OMP THREADPRIVATE(ok_restrict_A_satlim, modif_ratqs)
     85   REAL,    SAVE :: Pcrit_ratqs, ratqsbasnew
     86!$OMP THREADPRIVATE(Pcrit_ratqs, ratqsbasnew)
     87   REAL,    SAVE :: fac_modif_evaoce
    11688!$OMP THREADPRIVATE(fac_modif_evaoce)
    117 integer, save :: ok_bidouille_wake
     89   INTEGER, SAVE :: ok_bidouille_wake
    11890!$OMP THREADPRIVATE(ok_bidouille_wake)
    119 logical :: cond_temp_env
     91   LOGICAL, SAVE :: cond_temp_env
    12092!$OMP THREADPRIVATE(cond_temp_env)
    12193
    122 
    123 ! variables tableaux fn de niso
    124 real, ALLOCATABLE, DIMENSION(:), save :: tnat, toce, tcorr
    125 !$OMP THREADPRIVATE(tnat, toce, tcorr)
    126 real, ALLOCATABLE, DIMENSION(:), save :: tdifrel
    127 !$OMP THREADPRIVATE(tdifrel)
    128 real, ALLOCATABLE, DIMENSION(:), save :: talph1, talph2, talph3
    129 !$OMP THREADPRIVATE(talph1, talph2, talph3)
    130 real, ALLOCATABLE, DIMENSION(:), save :: talps1, talps2
    131 !$OMP THREADPRIVATE(talps1, talps2)
    132 real, ALLOCATABLE, DIMENSION(:), save :: tkcin0, tkcin1, tkcin2
     94   !--- Vectors of length "niso"
     95   REAL, ALLOCATABLE, DIMENSION(:), SAVE :: &
     96                    tnat, toce, tcorr, tdifrel
     97!$OMP THREADPRIVATE(tnat, toce, tcorr, tdifrel)
     98   REAL, ALLOCATABLE, DIMENSION(:), SAVE :: &
     99                    talph1, talph2, talph3, talps1, talps2
     100!$OMP THREADPRIVATE(talph1, talph2, talph3, talps1, talps2)
     101   REAL, ALLOCATABLE, DIMENSION(:), SAVE :: &
     102                    tkcin0, tkcin1, tkcin2
    133103!$OMP THREADPRIVATE(tkcin0, tkcin1, tkcin2)
    134 real, ALLOCATABLE, DIMENSION(:), save :: alpha_liq_sol
    135 !$OMP THREADPRIVATE(alpha_liq_sol)
    136 real, ALLOCATABLE, DIMENSION(:), save :: Rdefault, Rmethox
    137 !$OMP THREADPRIVATE(Rdefault, Rmethox)
     104   REAL, ALLOCATABLE, DIMENSION(:), SAVE :: &
     105                    alpha_liq_sol, Rdefault, Rmethox
     106!$OMP THREADPRIVATE(alpha_liq_sol, Rdefault, Rmethox)
    138107character*3, ALLOCATABLE, DIMENSION(:), save :: striso
    139108!$OMP THREADPRIVATE(striso)
    140 real, save :: fac_coeff_eq17_liq, fac_coeff_eq17_ice
     109   REAL, SAVE ::    fac_coeff_eq17_liq, fac_coeff_eq17_ice
    141110!$OMP THREADPRIVATE(fac_coeff_eq17_liq, fac_coeff_eq17_ice)
    142111
    143       real ridicule ! valeur maximale pour qu'une variable de type
    144                     ! rapoport de mélange puisse être considérée comme négligeable. Si
    145                     ! négligeable, alors on ne verifie pas si sa compo iso esta bérrante.
    146       parameter (ridicule=1e-12)     
    147 !      parameter (ridicule=1)
    148 !
    149       real ridicule_rain ! valeur limite de ridicule pour les flux de pluies (rain, zrfl...)
    150       parameter (ridicule_rain=1e-8) ! en kg/s <-> 1e-3mm/day
    151 
    152       real ridicule_evap ! valeur limite de ridicule pour les evap
    153       parameter (ridicule_evap=ridicule_rain*1e-2) ! en kg/s <-> 1e-3mm/day
    154 
    155       real ridicule_qsol ! valeur limite de ridicule pour les qsol
    156       parameter (ridicule_qsol=ridicule_rain) ! en kg <-> 1e-8kg
    157 
    158       real ridicule_snow ! valeur limite de ridicule pour les snow
    159       parameter (ridicule_snow=ridicule_qsol) ! en kg/s <-> 1e-8kg
    160      
    161         real expb_max
    162         parameter (expb_max=30.0)
    163 
    164         ! spécifique au tritium:
    165        
    166 
    167 logical, save :: ok_prod_nucl_tritium ! si oui, production de tritium par essais nucleaires
     112   !--- Negligible lower thresholds: no need to check for absurd values under these lower limits
     113   REAL, PARAMETER :: &
     114      ridicule      = 1e-12,              & ! For mixing ratios
     115      ridicule_rain = 1e-8,               & ! For rain fluxes (rain, zrfl...) in kg/s <-> 1e-3 mm/day
     116      ridicule_evap = ridicule_rain*1e-2, & ! For evaporations                in kg/s <-> 1e-3 mm/day
     117      ridicule_qsol = ridicule_rain,      & ! For qsol                        in kg <-> 1e-8 kg
     118      ridicule_snow = ridicule_qsol         ! For snow                        in kg <-> 1e-8 kg
     119   REAL, PARAMETER :: expb_max = 30.0
     120!$OMP THREADPRIVATE(ridicule, ridicule_rain, ridicule_evap, ridicule_qsol, ridicule_snow, expb_max)
     121
     122   !--- Specific to HTO:
     123   LOGICAL, SAVE :: ok_prod_nucl_tritium    !--- TRUE => HTO production by nuclear tests
    168124!$OMP THREADPRIVATE(ok_prod_nucl_tritium)
    169         integer nessai
    170         parameter (nessai=486)
    171         integer, save :: day_nucl(nessai)
    172 !$OMP THREADPRIVATE(day_nucl)
    173         integer, save :: month_nucl(nessai)
    174 !$OMP THREADPRIVATE(month_nucl)
    175         integer, save :: year_nucl(nessai)
    176 !$OMP THREADPRIVATE(year_nucl)
    177         real, save :: lat_nucl(nessai)
    178 !$OMP THREADPRIVATE(lat_nucl)
    179         real, save :: lon_nucl(nessai)
    180 !$OMP THREADPRIVATE(lon_nucl)
    181         real, save :: zmin_nucl(nessai)
    182 !$OMP THREADPRIVATE(zmin_nucl)
    183         real, save :: zmax_nucl(nessai)
    184 !$OMP THREADPRIVATE(zmax_nucl)
    185         real, save :: HTO_nucl(nessai)
    186 !$OMP THREADPRIVATE(HTO_nucl)
    187 
     125   INTEGER, PARAMETER :: nessai = 486
     126   INTEGER, DIMENSION(nessai), SAVE :: &
     127                    day_nucl, month_nucl, year_nucl
     128!$OMP THREADPRIVATE(day_nucl, month_nucl, year_nucl)
     129   REAL,    DIMENSION(nessai), SAVE :: &
     130                    lat_nucl, lon_nucl, zmin_nucl, zmax_nucl, HTO_nucl
     131!$OMP THREADPRIVATE(lat_nucl, lon_nucl, zmin_nucl, zmax_nucl, HTO_nucl)
     132 
    188133 
    189134CONTAINS
    190135
    191   SUBROUTINE iso_init()
    192       use ioipsl_getin_p_mod, ONLY : getin_p
    193       implicit none
    194 
    195 ! -- local variables:
    196 
    197       integer ixt
    198       ! référence O18
    199       real fac_enrichoce18
    200       real alpha_liq_sol_O18, &
    201      &     talph1_O18,talph2_O18,talph3_O18, &
    202      &     talps1_O18,talps2_O18, &
    203      &     tkcin0_O18,tkcin1_O18,tkcin2_O18, &
    204      &     tdifrel_O18 
     136SUBROUTINE iso_init()
     137   USE ioipsl_getin_p_mod, ONLY: getin_p
     138   USE infotrac_phy,       ONLY: ntiso, niso, isoName
     139   IMPLICIT NONE
     140
     141   !=== Local variables:
     142   INTEGER :: ixt
     143
     144   !--- H2[18]O reference
     145   REAL :: fac_enrichoce18, alpha_liq_sol_O18, &
     146           talph1_O18, talph2_O18, talph3_O18, talps1_O18, talps2_O18, &
     147           tkcin0_O18, tkcin1_O18, tkcin2_O18, tdifrel_O18 
     148
     149   !--- For H2[17]O
     150   REAL    :: fac_kcin, pente_MWL
     151   INTEGER :: ierr
    205152     
    206       ! cas de l'O17
    207       real fac_kcin
    208       real pente_MWL
    209       integer ierr
    210      
    211       logical ok_nocinsat, ok_nocinsfc !sensi test
    212       parameter (ok_nocinsfc=.FALSE.)  ! if T: no kinetic effect in sfc evap
    213       parameter (ok_nocinsat=.FALSE.)  ! if T: no sursaturation effect for ice
    214       logical Rdefault_smow
    215       parameter (Rdefault_smow=.FALSE.) ! si T: Rdefault=smow; si F: nul
    216       ! pour le tritium
    217       integer iessai
    218 
    219     write(*,*) 'iso_init 219: entree'
    220 
    221 ! allocations mémoire
    222 allocate (tnat(niso))
    223 allocate (toce(niso))
    224 allocate (tcorr(niso))
    225 allocate (tdifrel(niso))
    226 allocate (talph1(niso))
    227 allocate (talph2(niso))
    228 allocate (talph3(niso))
    229 allocate (talps1(niso))
    230 allocate (talps2(niso))
    231 allocate (tkcin0(niso))
    232 allocate (tkcin1(niso))
    233 allocate (tkcin2(niso))
    234 allocate (alpha_liq_sol(niso))
    235 allocate (Rdefault(niso))
    236 allocate (Rmethox(niso))
    237 allocate (striso(niso))
    238 
    239 
    240 !--------------------------------------------------------------
    241 ! General:
    242 !--------------------------------------------------------------
    243 
    244 ! -- verif du nombre d'isotopes:
    245       write(*,*) 'iso_init 64: niso=',niso
    246 
    247 ! init de ntracisoOR: on écrasera en cas de nzone>0 si complications avec
    248 ! ORCHIDEE
    249       ntracisoOR=ntraciso 
    250              
    251 ! -- Type of water isotopes:
    252 
    253         iso_eau=indnum_fn_num(1)
    254         iso_HDO=indnum_fn_num(2)
    255         iso_O18=indnum_fn_num(3)
    256         iso_O17=indnum_fn_num(4)
    257         iso_HTO=indnum_fn_num(5)
    258         write(*,*) 'iso_init 59: iso_eau=',iso_eau
    259         write(*,*) 'iso_HDO=',iso_HDO
    260         write(*,*) 'iso_O18=',iso_O18
    261         write(*,*) 'iso_O17=',iso_O17
    262         write(*,*) 'iso_HTO=',iso_HTO
    263         write(*,*) 'iso_init 251: use_iso=',use_iso
    264 
    265       ! initialisation
    266         lambda_sursat=0.004
    267         thumxt1=0.75*1.2
    268         ntot=20
    269         h_land_ice=20. ! à comparer aux 3000mm de snow_max
    270         P_veg=1.0
    271         bidouille_anti_divergence=.false.
    272         essai_convergence=.false.
    273         initialisation_iso=0
    274         modif_sst=0
    275         modif_sic=0
    276         deltaTtest=0.0
    277         deltasic=0.1
    278         deltaTtestpoles=0.0
    279         sstlatcrit=30.0
    280         deltaO18_oce=0.0
    281         albedo_prescrit=0
    282         lon_min_albedo=-200
    283         lon_max_albedo=200
    284         lat_min_albedo=-100
    285         lat_max_albedo=100
    286         deltaP_BL=10.0
    287         ruissellement_pluie=0
    288         alphak_stewart=1
    289         tdifexp_sol=0.67
    290         calendrier_guide=0
    291         cste_surf_cond=0
    292 mixlen=35.0       
    293 evap_cont_cste=0.0
    294 deltaO18_evap_cont=0.0
    295 d_evap_cont=0.0
    296 nudge_qsol=0
    297 region_nudge_qsol=1
    298 nlevmaxO17=50
    299 no_pce=0
    300 A_satlim=1.0
    301 ok_restrict_A_satlim=0
    302 !        slope_limiterxy=2.0
    303 !        slope_limiterz=2.0
    304 modif_ratqs=0
    305 Pcrit_ratqs=500.0
    306 ratqsbasnew=0.05
    307 
    308 fac_modif_evaoce=1.0
    309 ok_bidouille_wake=0
    310 cond_temp_env=.false.
    311 ! si oui, la temperature de cond est celle de l'environnement,
    312 ! pour eviter bugs quand temperature dans ascendances convs est
    313 ! mal calculee
    314 ok_prod_nucl_tritium=.false.
    315 
    316 ! lecture des paramètres isotopiques:
    317 ! pour que ça marche en openMP, il faut utiliser getin_p. Car le getin ne peut
    318 ! être appelé que par un thread à la fois, et ça pose tout un tas de problème,
    319 ! d'où tout un tas de magouilles comme dans conf_phys_m. A terme, tout le monde
    320 ! lira par getin_p.
    321 call getin_p('lambda',lambda_sursat)
    322 call getin_p('thumxt1',thumxt1)
    323 call getin_p('ntot',ntot)
    324 call getin_p('h_land_ice',h_land_ice)
    325 call getin_p('P_veg',P_veg)
    326 call getin_p('bidouille_anti_divergence',bidouille_anti_divergence)
    327 call getin_p('essai_convergence',essai_convergence)
    328 call getin_p('initialisation_iso',initialisation_iso)
    329 !if (nzone>0) then     
    330 !if (initialisation_iso.eq.0) then
    331 !  call getin_p('initialisation_isotrac',initialisation_isotrac)
    332 !endif !if (initialisation_iso.eq.0) then
    333 !endif !if (nzone>0)
    334 call getin_p('modif_sst',modif_sst)
    335 if (modif_sst.ge.1) then
    336 call getin_p('deltaTtest',deltaTtest)
    337 if (modif_sst.ge.2) then
    338   call getin_p('deltaTtestpoles',deltaTtestpoles)
    339   call getin_p('sstlatcrit',sstlatcrit)
     153   !--- Sensitivity tests
     154   LOGICAL, PARAMETER ::   ok_nocinsfc = .FALSE. ! if T: no kinetic effect in sfc evap
     155   LOGICAL, PARAMETER ::   ok_nocinsat = .FALSE. ! if T: no sursaturation effect for ice
     156   LOGICAL, PARAMETER :: Rdefault_smow = .FALSE. ! if T: Rdefault=smow; if F: nul
     157
     158   !--- For [3]H
     159   INTEGER :: iessai
     160
     161   CHARACTER(LEN=maxlen) :: modname, sxt
     162
     163   modname = 'iso_init'
     164   CALL msg('219: entree', modname)
     165
     166   !--- Memory allocations
     167   ALLOCATE(talph1(niso), tkcin0(niso),  talps1(niso),  tnat(niso))
     168   ALLOCATE(talph2(niso), tkcin1(niso),  talps2(niso),  toce(niso))
     169   ALLOCATE(talph3(niso), tkcin2(niso), tdifrel(niso), tcorr(niso))
     170   ALLOCATE(alpha_liq_sol(niso),   Rdefault(niso),   Rmethox(niso))
     171   ALLOCATE(striso(niso))
     172
     173
     174   !--------------------------------------------------------------
     175   ! General:
     176   !--------------------------------------------------------------
     177
     178   !--- Check number of isotopes
     179   CALL msg('64: niso = '//TRIM(int2str(niso)), modname)
     180
     181   !--- Init de ntracisoOR: on ecrasera en cas de traceurs de tagging isotopiques
     182   !                     (nzone>0) si complications avec ORCHIDEE
     183   ntracisoOR = ntiso 
     184
     185   !--- Type of water isotopes:
     186   iso_eau = strIdx(isoName, 'H2[16]O'); CALL msg('59: iso_eau='//int2str(iso_eau), modname)
     187   iso_O17 = strIdx(isoName, 'H2[17]O'); CALL msg('iso_HDO='//int2str(iso_HDO), modname)
     188   iso_O18 = strIdx(isoName, 'H2[18]O'); CALL msg('iso_O18='//int2str(iso_O18), modname)
     189   iso_HDO = strIdx(isoName, 'H[2]HO');  CALL msg('iso_O17='//int2str(iso_O17), modname)
     190   iso_HTO = strIdx(isoName, 'H[3]HO');  CALL msg('iso_HTO='//int2str(iso_HTO), modname)
     191
     192   ! initialisation
     193   ! lecture des parametres isotopiques:
     194   ! pour que ca marche en openMP, il faut utiliser getin_p. Car le getin ne peut
     195   ! etre appele que par un thread a la fois, et ca pose tout un tas de problemes,
     196   ! d'ou tout un tas de magouilles comme dans conf_phys_m. A terme, tout le monde
     197   ! lira par getin_p.
     198   CALL get_in('lambda',     lambda_sursat, 0.004)
     199   CALL get_in('thumxt1',    thumxt1,       0.75*1.2)
     200   CALL get_in('ntot',       ntot,          20,  .FALSE.)
     201   CALL get_in('h_land_ice', h_land_ice,    20., .FALSE.)
     202   CALL get_in('P_veg',      P_veg,         1.0, .FALSE.)
     203   CALL get_in('bidouille_anti_divergence', bidouille_anti_divergence, .FALSE.)
     204   CALL get_in('essai_convergence',         essai_convergence,         .FALSE.)
     205   CALL get_in('initialisation_iso',        initialisation_iso,        0)
     206
     207!  IF(nzone>0 .AND. initialisation_iso==0) &
     208!      CALL get_in('initialisation_isotrac',initialisation_isotrac)
     209   CALL get_in('modif_sst',      modif_sst,         0)
     210   CALL get_in('deltaTtest',     deltaTtest,      0.0)     !--- For modif_sst>=1
     211   CALL get_in('deltaTtestpoles',deltaTtestpoles, 0.0)     !--- For modif_sst>=2
     212   CALL get_in( 'sstlatcrit',    sstlatcrit,     30.0)     !--- For modif_sst>=3
     213   CALL get_in('dsstlatcrit',   dsstlatcrit,      0.0)     !--- For modif_sst>=3
    340214#ifdef ISOVERIF
    341       !call iso_verif_positif(sstlatcrit,'iso_init 107')
    342       if (sstlatcrit.lt.0.0) then
    343         write(*,*) 'iso_init 270: sstlatcrit=',sstlatcrit
    344         stop
    345       endif
     215   CALL msg('iso_init 270:  sstlatcrit='//real2str( sstlatcrit), modname, sstlatcrit < 0.0) !--- For modif_sst>=2
     216   CALL msg('iso_init 279: dsstlatcrit='//real2str(dsstlatcrit), modname, sstlatcrit < 0.0) !--- For modif_sst>=3
     217   IF(modif_sst >= 2 .AND. sstlatcrit < 0.0) STOP
    346218#endif             
    347   if (modif_sst.ge.3) then 
    348       call getin_p('dsstlatcrit',dsstlatcrit)
     219
     220   CALL get_in('modif_sic', modif_sic,  0)
     221   IF(modif_sic >= 1) &
     222   CALL get_in('deltasic',  deltasic, 0.1)
     223
     224   CALL get_in('albedo_prescrit', albedo_prescrit, 0)
     225   IF(albedo_prescrit == 1) THEN
     226      CALL get_in('lon_min_albedo', lon_min_albedo, -200.)
     227      CALL get_in('lon_max_albedo', lon_max_albedo,  200.)
     228      CALL get_in('lat_min_albedo', lat_min_albedo, -100.)
     229      CALL get_in('lat_max_albedo', lat_max_albedo,  100.)
     230   END IF
     231   CALL get_in('deltaO18_oce',        deltaO18_oce,   0.0)
     232
     233   CALL get_in('deltaP_BL',           deltaP_BL,     10.0)
     234   CALL get_in('ruissellement_pluie', ruissellement_pluie, 0)
     235   CALL get_in('alphak_stewart',      alphak_stewart,      1)
     236   CALL get_in('tdifexp_sol',         tdifexp_sol,      0.67)
     237   CALL get_in('calendrier_guide',    calendrier_guide,    0)
     238   CALL get_in('cste_surf_cond',      cste_surf_cond,      0)
     239   CALL get_in('mixlen',              mixlen,           35.0)
     240   CALL get_in('evap_cont_cste',      evap_cont_cste,      0)
     241   CALL get_in('deltaO18_evap_cont',  deltaO18_evap_cont,0.0)
     242   CALL get_in('d_evap_cont',         d_evap_cont,       0.0)
     243   CALL get_in('nudge_qsol',          nudge_qsol,          0)
     244   CALL get_in('region_nudge_qsol',   region_nudge_qsol,   1)
     245   nlevmaxO17 = 50
     246   CALL msg('nlevmaxO17='//TRIM(int2str(nlevmaxO17)))
     247   CALL get_in('no_pce',   no_pce,     0)
     248   CALL get_in('A_satlim', A_satlim, 1.0)
     249   CALL get_in('ok_restrict_A_satlim', ok_restrict_A_satlim, 0)
    349250#ifdef ISOVERIF
    350       !call iso_verif_positif(dsstlatcrit,'iso_init 110')
    351       if (sstlatcrit.lt.0.0) then
    352         write(*,*) 'iso_init 279: dsstlatcrit=',dsstlatcrit
    353         stop
    354       endif
    355 #endif             
    356   endif !if (modif_sst.ge.3) then
    357 endif !if (modif_sst.ge.2) then
    358 endif !  if (modif_sst.ge.1) then   
    359 call getin_p('modif_sic',modif_sic)
    360 if (modif_sic.ge.1) then
    361 call getin_p('deltasic',deltasic)
    362 endif !if (modif_sic.ge.1) then
    363 
    364 call getin_p('albedo_prescrit',albedo_prescrit)
    365 call getin_p('lon_min_albedo',lon_min_albedo)
    366 call getin_p('lon_max_albedo',lon_max_albedo)
    367 call getin_p('lat_min_albedo',lat_min_albedo)
    368 call getin_p('lat_max_albedo',lat_max_albedo)
    369 call getin_p('deltaO18_oce',deltaO18_oce)
    370 call getin_p('deltaP_BL',deltaP_BL)
    371 call getin_p('ruissellement_pluie',ruissellement_pluie)
    372 call getin_p('alphak_stewart',alphak_stewart)
    373 call getin_p('tdifexp_sol',tdifexp_sol)
    374 call getin_p('calendrier_guide',calendrier_guide)
    375 call getin_p('cste_surf_cond',cste_surf_cond)
    376 call getin_p('mixlen',mixlen)
    377 call getin_p('evap_cont_cste',evap_cont_cste)
    378 call getin_p('deltaO18_evap_cont',deltaO18_evap_cont)
    379 call getin_p('d_evap_cont',d_evap_cont) 
    380 call getin_p('nudge_qsol',nudge_qsol)
    381 call getin_p('region_nudge_qsol',region_nudge_qsol)
    382 call getin_p('no_pce',no_pce)
    383 call getin_p('A_satlim',A_satlim)
    384 call getin_p('ok_restrict_A_satlim',ok_restrict_A_satlim)
    385 #ifdef ISOVERIF     
    386 !call iso_verif_positif(1.0-A_satlim,'iso_init 158')
    387       if (A_satlim.gt.1.0) then
    388         write(*,*) 'iso_init 315: A_satlim=',A_satlim
    389         stop
    390       endif
    391 #endif         
    392 !      call getin_p('slope_limiterxy',slope_limiterxy)
    393 !      call getin_p('slope_limiterz',slope_limiterz)
    394 call getin_p('modif_ratqs',modif_ratqs)
    395 call getin_p('Pcrit_ratqs',Pcrit_ratqs)
    396 call getin_p('ratqsbasnew',ratqsbasnew)
    397 call getin_p('fac_modif_evaoce',fac_modif_evaoce)
    398 call getin_p('ok_bidouille_wake',ok_bidouille_wake)
    399 call getin_p('cond_temp_env',cond_temp_env)
    400 if (use_iso(iso_HTO_possible)) then
    401   ok_prod_nucl_tritium=.true.
    402   call getin_p('ok_prod_nucl_tritium',ok_prod_nucl_tritium)
    403 endif
    404 
    405 write(*,*) 'lambda,thumxt1=',lambda_sursat,thumxt1
    406 write(*,*) 'bidouille_anti_divergence=',bidouille_anti_divergence
    407 write(*,*) 'essai_convergence=',essai_convergence
    408 write(*,*) 'initialisation_iso=',initialisation_iso
    409 write(*,*) 'modif_sst=',modif_sst
    410 if (modif_sst.ge.1) then
    411 write(*,*) 'deltaTtest=',deltaTtest
    412 if (modif_sst.ge.2) then 
    413 write(*,*) 'deltaTtestpoles,sstlatcrit=', &
    414 &           deltaTtestpoles,sstlatcrit
    415 if (modif_sst.ge.3) then   
    416  write(*,*) 'dsstlatcrit=',dsstlatcrit
    417 endif !if (modif_sst.ge.3) then
    418 endif !if (modif_sst.ge.2) then
    419 endif !if (modif_sst.ge.1) then
    420 write(*,*) 'modif_sic=',modif_sic
    421 if (modif_sic.ge.1) then 
    422 write(*,*) 'deltasic=',deltasic
    423 endif !if (modif_sic.ge.1) then
    424 write(*,*) 'deltaO18_oce=',deltaO18_oce
    425 write(*,*) 'albedo_prescrit=',albedo_prescrit
    426 if (albedo_prescrit.eq.1) then
    427  write(*,*) 'lon_min_albedo,lon_max_albedo=', &
    428 &           lon_min_albedo,lon_max_albedo
    429  write(*,*) 'lat_min_albedo,lat_max_albedo=', &
    430 &           lat_min_albedo,lat_max_albedo
    431 endif !if (albedo_prescrit.eq.1) then
    432 write(*,*) 'deltaP_BL,ruissellement_pluie,alphak_stewart=', &
    433 &       deltaP_BL,ruissellement_pluie,alphak_stewart
    434 write(*,*) 'cste_surf_cond=',cste_surf_cond
    435 write(*,*) 'mixlen=',mixlen
    436 write(*,*) 'tdifexp_sol=',tdifexp_sol
    437 write(*,*) 'calendrier_guide=',calendrier_guide
    438 write(*,*) 'evap_cont_cste=',evap_cont_cste
    439 write(*,*) 'deltaO18_evap_cont,d_evap_cont=', &
    440 &           deltaO18_evap_cont,d_evap_cont
    441 write(*,*) 'nudge_qsol,region_nudge_qsol=', &
    442 &  nudge_qsol,region_nudge_qsol 
    443 write(*,*) 'nlevmaxO17=',nlevmaxO17
    444 write(*,*) 'no_pce=',no_pce
    445 write(*,*) 'A_satlim=',A_satlim
    446 write(*,*) 'ok_restrict_A_satlim=',ok_restrict_A_satlim
    447 !      write(*,*) 'slope_limiterxy=',slope_limiterxy
    448 !      write(*,*) 'slope_limiterz=',slope_limiterz
    449 write(*,*) 'modif_ratqs=',modif_ratqs
    450 write(*,*) 'Pcrit_ratqs=',Pcrit_ratqs
    451 write(*,*) 'ratqsbasnew=',ratqsbasnew
    452 write(*,*) 'fac_modif_evaoce=',fac_modif_evaoce
    453 write(*,*) 'ok_bidouille_wake=',ok_bidouille_wake
    454 write(*,*) 'cond_temp_env=',cond_temp_env
    455 write(*,*) 'ok_prod_nucl_tritium=',ok_prod_nucl_tritium
    456          
    457 
    458 !--------------------------------------------------------------
    459 ! Parameters that do not depend on the nature of water isotopes:
    460 !--------------------------------------------------------------
    461 
    462 ! -- temperature at which ice condensate starts to form (valeur ECHAM?):
    463 pxtmelt=273.15
    464 !      pxtmelt=273.15-10.0 ! test PHASE
    465 
    466 ! -- temperature at which all condensate is ice:
    467 pxtice=273.15-10.0
    468 !      pxtice=273.15-30.0 ! test PHASE
    469 
    470 ! -- minimum temperature to calculate fractionation coeff
    471 pxtmin=273.15-120.0 ! On ne calcule qu'au dessus de -120°C
    472 pxtmax=273.15+60.0 ! On ne calcule qu'au dessus de +60°C
    473 ! remarque: les coeffs ont été mesurés seulement jusq'à -40!
    474 
    475 ! -- a constant for alpha_eff for equilibrium below cloud base:
    476 tdifexp=0.58
    477 tv0cin=7.0
    478 
    479 ! facteurs lambda et mu dans Si=musi-lambda*T
    480 musi=1.0
    481 if (ok_nocinsat) then
    482 lambda_sursat = 0.0 ! no sursaturation effect
    483 endif           
    484 
    485 
    486 ! diffusion dans le sol
    487 Kd=2.5e-9 ! m2/s   
    488 
    489 ! cas où cste_surf_cond: on met rhs ou/et Ts cste pour voir
    490 rh_cste_surf_cond=0.6
    491 T_cste_surf_cond=288.0
    492 
    493 !--------------------------------------------------------------
    494 ! Parameters that depend on the nature of water isotopes:
    495 !--------------------------------------------------------------
    496 ! ** constantes locales
    497 fac_enrichoce18=0.0005
    498 ! on a alors tcor018=1+fac_enrichoce18
    499 ! tcorD=1+fac_enrichoce18*8
    500 ! tcorO17=1+fac_enrichoce18*0.528
    501 alpha_liq_sol_O18=1.00291 ! valeur de Lehmann & Siegenthaler, 1991,
    502   ! Journal of Glaciology, vol 37, p 23
    503 talph1_O18=1137.
    504 talph2_O18=-0.4156
    505 talph3_O18=-2.0667E-3
    506 talps1_O18=11.839
    507 talps2_O18=-0.028244
    508 tkcin0_O18 = 0.006
    509 tkcin1_O18 = 0.000285
    510 tkcin2_O18 = 0.00082
    511 tdifrel_O18= 1./0.9723
    512 
    513 ! rapport des ln(alphaeq) entre O18 et O17
    514 fac_coeff_eq17_liq=0.529 ! donné par Amaelle
    515 !      fac_coeff_eq17_ice=0.528 ! slope MWL
    516 fac_coeff_eq17_ice=0.529
    517 
    518 
    519 write(*,*) 'iso_O18,iso_HDO,iso_eau=',iso_O18,iso_HDO,iso_eau
    520 do 999 ixt = 1, niso
    521 write(*,*) 'iso_init 80: ixt=',ixt
    522 
    523 
    524 ! -- kinetic factor for surface evaporation:
    525 ! (cf: kcin = tkcin0                  if |V|<tv0cin
    526 !      kcin = tkcin1*|Vsurf| + tkcin2 if |V|>tv0cin )
    527 ! (Rq: formula discontinuous for |V|=tv0cin... )       
    528 
    529 ! -- main:
    530 if (ixt.eq.iso_HTO) then ! Tritium
    531   tkcin0(ixt) = 0.01056
    532   tkcin1(ixt) = 0.0005016
    533   tkcin2(ixt) = 0.0014432
    534   tnat(ixt)=0.
    535   !toce(ixt)=2.2222E-8 ! corrigé par Alex Cauquoin
    536   !toce(ixt)=1.0E-18 ! rapport 3H/1H ocean
    537   toce(ixt)=4.0E-19 ! rapport T/H = 0.2 TU Dreisigacker and Roether 1978
    538   tcorr(ixt)=1.
    539   tdifrel(ixt)=1./0.968
    540   talph1(ixt)=46480.
    541   talph2(ixt)=-103.87
    542   talph3(ixt)=0.
    543   talps1(ixt)=46480.
    544   talps2(ixt)=-103.87
    545   alpha_liq_sol(ixt)=1.
    546   Rdefault(ixt)=0.0
    547   Rmethox(ixt)=0.0
    548   striso(ixt)='HTO'
    549 endif     
    550 if (ixt.eq.iso_O17) then ! Deuterium
    551   pente_MWL=0.528
    552 !          tdifrel(ixt)=1./0.985452 ! donné par Amaelle
    553   tdifrel(ixt)=1./0.98555 ! valeur utilisée en 1D et dans modèle de LdG
    554 !          fac_kcin=0.5145 ! donné par Amaelle
    555   fac_kcin= (tdifrel(ixt)-1.0)/(tdifrel_O18-1.0)
    556   tkcin0(ixt) = tkcin0_O18*fac_kcin
    557   tkcin1(ixt) = tkcin1_O18*fac_kcin
    558   tkcin2(ixt) = tkcin2_O18*fac_kcin
    559   tnat(ixt)=0.004/100. ! O17 représente 0.004% de l'oxygène
    560   toce(ixt)=tnat(ixt)*(1.0+deltaO18_oce/1000.0)**pente_MWL
    561   tcorr(ixt)=1.0+fac_enrichoce18*pente_MWL ! donné par Amaelle           
    562   talph1(ixt)=talph1_O18
    563   talph2(ixt)=talph2_O18
    564   talph3(ixt)=talph3_O18
    565   talps1(ixt)=talps1_O18
    566   talps2(ixt)=talps2_O18     
    567   alpha_liq_sol(ixt)=(alpha_liq_sol_O18)**fac_coeff_eq17_liq
    568   if (Rdefault_smow) then   
    569         Rdefault(ixt)=tnat(ixt)*(-3.15/1000.0+1.0)
    570   else
    571         Rdefault(ixt)=0.0
    572   endif
    573   Rmethox(ixt)=(230./1000.+1.)*tnat(ixt) !Zahn et al 2006
    574   striso(ixt)='O17'
    575 endif
    576 
    577 if (ixt.eq.iso_O18) then ! Oxygene18
    578   tkcin0(ixt) = tkcin0_O18
    579   tkcin1(ixt) = tkcin1_O18
    580   tkcin2(ixt) = tkcin2_O18
    581   tnat(ixt)=2005.2E-6
    582   toce(ixt)=tnat(ixt)*(1.0+deltaO18_oce/1000.0)
    583   tcorr(ixt)=1.0+fac_enrichoce18
    584   tdifrel(ixt)=tdifrel_O18
    585   talph1(ixt)=talph1_O18
    586   talph2(ixt)=talph2_O18
    587   talph3(ixt)=talph3_O18
    588   talps1(ixt)=talps1_O18
    589   talps2(ixt)=talps2_O18
    590   alpha_liq_sol(ixt)=alpha_liq_sol_O18   
    591   if (Rdefault_smow) then   
    592         Rdefault(ixt)=tnat(ixt)*(-6.0/1000.0+1.0)
    593   else
    594         Rdefault(ixt)=0.0
    595   endif
    596   Rmethox(ixt)=(130./1000.+1.)*tnat(ixt) !Zahn et al 2006   
    597 !       write(*,*) 'iso_init 163: ZXalpha_liq_sol=',ZXalpha_liq_sol
    598   striso(ixt)='O18'
    599   write(*,*) 'isotopes_mod 519: ixt,striso(ixt)=',ixt,striso(ixt)
    600 endif
    601 
    602 if (ixt.eq.iso_HDO) then ! Deuterium
    603   pente_MWL=8.0
    604 !          fac_kcin=0.88
    605   tdifrel(ixt)=1./0.9755
    606   fac_kcin= (tdifrel(ixt)-1)/(tdifrel_O18-1)
    607   tkcin0(ixt) = tkcin0_O18*fac_kcin
    608   tkcin1(ixt) = tkcin1_O18*fac_kcin
    609   tkcin2(ixt) = tkcin2_O18*fac_kcin
    610   tnat(ixt)=155.76E-6
    611   toce(ixt)=tnat(ixt)*(1.0+pente_MWL*deltaO18_oce/1000.0)
    612   tcorr(ixt)=1.0+fac_enrichoce18*pente_MWL         
    613   talph1(ixt)=24844.
    614   talph2(ixt)=-76.248
    615   talph3(ixt)=52.612E-3
    616   talps1(ixt)=16288.
    617   talps2(ixt)=-0.0934
    618   !ZXalpha_liq_sol=1.0192 ! Weston, Ralph, 1955
    619   alpha_liq_sol(ixt)=1.0212
    620   ! valeur de Lehmann & Siegenthaler, 1991, Journal of
    621   ! Glaciology, vol 37, p 23
    622   if (Rdefault_smow) then   
    623     Rdefault(ixt)=tnat(ixt)*((-6.0*pente_MWL+10.0)/1000.0+1.0)
    624   else
    625     Rdefault(ixt)=0.0
    626   endif
    627   Rmethox(ixt)=tnat(ixt)*(-25.0/1000.+1.) ! Zahn et al 2006
    628   striso(ixt)='HDO'
    629   write(*,*) 'isotopes_mod 548: ixt,striso(ixt)=',ixt,striso(ixt)
    630 endif
    631 
    632 !       write(*,*) 'iso_init 163: ZXalpha_liq_sol=',ZXalpha_liq_sol
    633 if (ixt.eq.iso_eau) then ! Oxygene16
    634   tkcin0(ixt) = 0.0
    635   tkcin1(ixt) = 0.0
    636   tkcin2(ixt) = 0.0
    637   tnat(ixt)=1.
    638   toce(ixt)=tnat(ixt)
    639   tcorr(ixt)=1.0
    640   tdifrel(ixt)=1.
    641   talph1(ixt)=0.
    642   talph2(ixt)=0.
    643   talph3(ixt)=0.
    644   talps1(ixt)=0.
    645   talph3(ixt)=0.
    646   alpha_liq_sol(ixt)=1.
    647   if (Rdefault_smow) then
    648         Rdefault(ixt)=tnat(ixt)*1.0
    649   else
    650         Rdefault(ixt)=1.0
    651   endif
    652   Rmethox(ixt)=1.0
    653   striso(ixt)='eau'
    654 endif
    655 
    656 999   continue
    657 
    658 ! test de sensibilité:
    659 if (ok_nocinsfc) then ! no kinetic effect in sfc evaporation
    660  do ixt=1,niso
    661   tkcin0(ixt) = 0.0
    662   tkcin1(ixt) = 0.0
    663   tkcin2(ixt) = 0.0
    664  enddo
    665 endif
    666 
    667 ! nom des isotopes
    668 
    669 ! verif
    670 write(*,*) 'iso_init 285: verif initialisation:'
    671 
    672 do ixt=1,niso
    673   write(*,*) '* striso(',ixt,')=<'//striso(ixt)//'>'
    674   write(*,*) 'tnat(',ixt,')=',tnat(ixt)
    675 !          write(*,*) 'alpha_liq_sol(',ixt,')=',alpha_liq_sol(ixt)
    676 !          write(*,*) 'tkcin0(',ixt,')=',tkcin0(ixt)
    677 !          write(*,*) 'tdifrel(',ixt,')=',tdifrel(ixt)
    678 enddo
    679 write(*,*) 'iso_init 69: lambda=',lambda_sursat
    680 write(*,*) 'iso_init 69: thumxt1=',thumxt1
    681 write(*,*) 'iso_init 69: h_land_ice=',h_land_ice
    682 write(*,*) 'iso_init 69: P_veg=',P_veg   
    683 
    684     return
     251   CALL msg(' 315: A_satlim='//real2str(A_satlim), modname, A_satlim > 1.0)
     252   IF(A_satlim > 1.0) STOP
     253#endif
     254!  CALL get_in('slope_limiterxy',   slope_limiterxy,  2.0)
     255!  CALL get_in('slope_limiterz',    slope_limiterz,   2.0)
     256   CALL get_in('modif_ratqs',       modif_ratqs,        0)
     257   CALL get_in('Pcrit_ratqs',       Pcrit_ratqs,    500.0)
     258   CALL get_in('ratqsbasnew',       ratqsbasnew,     0.05)
     259   CALL get_in('fac_modif_evaoce',  fac_modif_evaoce, 1.0)
     260   CALL get_in('ok_bidouille_wake', ok_bidouille_wake,  0)
     261   ! si oui, la temperature de cond est celle de l'environnement, pour eviter
     262   ! bugs quand temperature dans ascendances convs est mal calculee
     263   CALL get_in('cond_temp_env',        cond_temp_env,        .FALSE.)
     264   IF(ANY(isoName == 'H[3]HO')) &
     265   CALL get_in('ok_prod_nucl_tritium', ok_prod_nucl_tritium, .FALSE., .FALSE.)
     266
     267   !--------------------------------------------------------------
     268   ! Parameters that do not depend on the nature of water isotopes:
     269   !--------------------------------------------------------------
     270   ! -- temperature at which ice condensate starts to form (valeur ECHAM?):
     271   pxtmelt = 273.15
     272
     273   ! -- temperature at which all condensate is ice:
     274   pxtice  = 273.15-10.0
     275
     276   !- -- test PHASE
     277!   pxtmelt = 273.15 - 10.0
     278!   pxtice  = 273.15 - 30.0
     279
     280   ! -- minimum temperature to calculate fractionation coeff
     281   pxtmin = 273.15 - 120.0   ! On ne calcule qu'au dessus de -120°C
     282   pxtmax = 273.15 +  60.0   ! On ne calcule qu'au dessus de +60°C
     283   !    Remarque: les coeffs ont ete mesures seulement jusq'à -40!
     284
     285   ! -- a constant for alpha_eff for equilibrium below cloud base:
     286   tdifexp = 0.58
     287   tv0cin  = 7.0
     288
     289   ! facteurs lambda et mu dans Si=musi-lambda*T
     290   musi=1.0
     291   if (ok_nocinsat) lambda_sursat = 0.0          ! no sursaturation effect
     292
     293   ! diffusion dans le sol
     294   Kd=2.5e-9 ! m2/s   
     295
     296   ! cas où cste_surf_cond: on met rhs ou/et Ts cste pour voir
     297   rh_cste_surf_cond = 0.6
     298    T_cste_surf_cond = 288.0
     299   
     300   !--------------------------------------------------------------
     301   ! Parameters that depend on the nature of water isotopes:
     302   !--------------------------------------------------------------
     303   ! Local constants
     304   fac_enrichoce18 = 0.0005            ! Then: tcorO18 = 1 + fac_enrichoce18
     305                                       !       tcorD   = 1 + fac_enrichoce18*8
     306                                       !       tcorO17 = 1 + fac_enrichoce18*0.528
     307   alpha_liq_sol_O18 = 1.00291         ! From Lehmann & Siegenthaler, 1991,
     308                                       ! Journal of Glaciology, vol 37, p 23
     309   talph1_O18 = 1137.  ; talph2_O18 = -0.4156   ; talph3_O18 = -2.0667E-3
     310   talps1_O18 = 11.839 ; talps2_O18 = -0.028244
     311   tkcin0_O18 = 0.006  ; tkcin1_O18 =  0.000285 ; tkcin2_O18 =  0.00082
     312   tdifrel_O18 = 1./0.9723
     313
     314   ! ln(alphaeq) ratio between O18 and O17
     315   fac_coeff_eq17_liq = 0.529          ! From Amaelle
     316  !fac_coeff_eq17_ice = 0.528          ! slope MWL
     317   fac_coeff_eq17_ice = 0.529
     318
     319   CALL msg('iso_O18, iso_HDO, iso_eau = '//TRIM(strStack(int2str([iso_O18, iso_HDO, iso_eau]))), modname)
     320
     321   !--- Kinetic factor for surface evaporation:
     322   ! (cf: kcin = tkcin0                  if |V|<tv0cin
     323   !      kcin = tkcin1*|Vsurf| + tkcin2 if |V|>tv0cin )
     324   ! (Rq: formula discontinuous for |V|=tv0cin... )       
     325
     326   DO ixt = 1, niso
     327      sxt=int2str(ixt)
     328      WRITE(*,*) 'iso_init 80: ixt=',ixt
     329
     330      Rdefault(ixt) = 0.0
     331      IF(ixt == iso_HTO) THEN          !=== H[3]HO
     332         tdifrel(ixt) = 1./0.968
     333         tkcin0(ixt) = 0.01056
     334         tkcin1(ixt) = 0.0005016
     335         tkcin2(ixt) = 0.0014432
     336         tnat  (ixt) = 0.
     337         toce  (ixt) = 4.0E-19         ! Ratio T/H = 0.2 TU, Dreisigacker and Roether 1978
     338        !toce  (ixt) = 2.2222E-8       ! Corrected by Alex Cauquoin
     339        !toce  (ixt) = 1.0E-18         ! Ratio 3H/1H ocean
     340         tcorr (ixt) = 1.
     341         talph1(ixt) = 46480. ; talph2(ixt) = -103.87 ; talph3(ixt) = 0.
     342         talps1(ixt) = 46480. ; talps2(ixt) = -103.87
     343         alpha_liq_sol(ixt) = 1.
     344         Rmethox(ixt) = 0.0
     345         striso (ixt) = 'HTO'
     346      ELSE IF(ixt == iso_O17) THEN     !=== H2[17]O
     347         tdifrel(ixt)=1./0.98555       ! Used in 1D and in LdG's model
     348        !tdifrel(ixt)=1./0.985452      ! From Amaelle
     349        !fac_kcin=0.5145               ! From Amaelle
     350         fac_kcin= (tdifrel(ixt)-1.0)/(tdifrel_O18-1.0)
     351         tkcin0(ixt) = tkcin0_O18*fac_kcin
     352         tkcin1(ixt) = tkcin1_O18*fac_kcin
     353         tkcin2(ixt) = tkcin2_O18*fac_kcin
     354         tnat  (ixt) = 0.004/100.      ! O17 = 0.004% of oxygen
     355         pente_MWL=0.528
     356         toce  (ixt) = tnat(ixt)*(1.0+deltaO18_oce/1000.0)**pente_MWL
     357         tcorr (ixt) = 1.0+fac_enrichoce18*pente_MWL  ! From Amaelle           
     358         talph1(ixt) = talph1_O18 ; talph2(ixt) = talph2_O18 ; talph3(ixt) = talph3_O18
     359         talps1(ixt) = talps1_O18 ; talps2(ixt) = talps2_O18     
     360         alpha_liq_sol(ixt) = (alpha_liq_sol_O18)**fac_coeff_eq17_liq
     361         IF(Rdefault_smow) Rdefault(ixt) = tnat(ixt)*(-3.15/1000.0+1.0)
     362         Rmethox(ixt) = (230./1000.+1.)*tnat(ixt)     ! Zahn et al 2006
     363         striso (ixt) = 'O17'
     364      ELSE IF(ixt == iso_O18) THEN     !=== H2[18]O
     365         tdifrel(ixt) = tdifrel_O18
     366         tkcin0(ixt) = tkcin0_O18
     367         tkcin1(ixt) = tkcin1_O18
     368         tkcin2(ixt) = tkcin2_O18
     369         tnat  (ixt) = 2005.2E-6
     370         toce  (ixt) = tnat(ixt)*(1.0+deltaO18_oce/1000.0)
     371         tcorr (ixt) = 1.0+fac_enrichoce18
     372         talph1(ixt) = talph1_O18 ; talph2(ixt) = talph2_O18 ; talph3(ixt) = talph3_O18
     373         talps1(ixt) = talps1_O18 ; talps2(ixt) = talps2_O18
     374         alpha_liq_sol(ixt) = alpha_liq_sol_O18   
     375         IF(Rdefault_smow) Rdefault(ixt) = tnat(ixt)*(-6.0/1000.0+1.0)
     376         Rmethox(ixt) = (130./1000.+1.)*tnat(ixt) ! Zahn et al 2006   
     377         striso (ixt) = 'O18'
     378         CALL msg('519: ixt, striso(ixt) = '//TRIM(sxt)//', '//TRIM(striso(ixt)), modname)
     379      ELSE IF(ixt == iso_HDO) THEN     !=== H[2]HO
     380         tdifrel(ixt) = 1./0.9755
     381        !fac_kcin=0.88
     382         fac_kcin= (tdifrel(ixt)-1.0)/(tdifrel_O18-1.0)
     383         tkcin0(ixt) = tkcin0_O18*fac_kcin
     384         tkcin1(ixt) = tkcin1_O18*fac_kcin
     385         tkcin2(ixt) = tkcin2_O18*fac_kcin
     386         tnat  (ixt) = 155.76E-6
     387         pente_MWL = 8.0
     388         toce  (ixt) = tnat(ixt)*(1.0+pente_MWL*deltaO18_oce/1000.0)
     389         tcorr (ixt) = 1.0+fac_enrichoce18*pente_MWL         
     390         talph1(ixt) = 24844. ; talph2(ixt) = -76.248 ; talph3(ixt) = 52.612E-3
     391         talps1(ixt) = 16288. ; talps2(ixt) = -0.0934
     392        !alpha_liq_sol(ixt)=1.0192 ZX  ! From Weston, Ralph, 1955
     393         alpha_liq_sol(ixt)=1.0212     ! From Lehmann & Siegenthaler, 1991,
     394                                       ! Journal of Glaciology, vol 37, p 23
     395         IF(Rdefault_smow) Rdefault(ixt) = tnat(ixt)*((-6.0*pente_MWL+10.0)/1000.0+1.0)
     396         Rmethox(ixt) = tnat(ixt)*(-25.0/1000.+1.) ! Zahn et al 2006
     397         striso (ixt) = 'HDO'
     398         CALL msg('548: ixt,striso(ixt) = '//TRIM(sxt)//', '//striso(ixt), modname)
     399      ELSE IF(ixt  == iso_eau) THEN    !=== H2O[16]
     400         tkcin0(ixt) = 0.0
     401         tkcin1(ixt) = 0.0
     402         tkcin2(ixt) = 0.0
     403         tnat  (ixt) = 1.
     404         toce  (ixt)=tnat(ixt)
     405         tcorr (ixt) = 1.0
     406         tdifrel(ixt) = 1.
     407         talph1(ixt) = 0. ; talph2(ixt) = 0. ; talph3(ixt) = 0.
     408         talps1(ixt) = 0. ; talph3(ixt) = 0.
     409         alpha_liq_sol(ixt)=1.
     410         IF(Rdefault_smow) Rdefault(ixt) = tnat(ixt)*1.0
     411         Rmethox(ixt) = 1.0
     412         striso(ixt) = 'eau'
     413      END IF
     414   END DO
     415
     416   !--- Sensitivity test: no kinetic effect in sfc evaporation
     417   IF(ok_nocinsfc) THEN
     418      tkcin0(1:niso) = 0.0
     419      tkcin1(1:niso) = 0.0
     420      tkcin2(1:niso) = 0.0
     421   END IF
     422
     423   CALL msg('285: verif initialisation:', modname)
     424   DO ixt=1,niso
     425      CALL msg(' * striso('//TRIM(sxt)//') = <'//TRIM(striso(ixt))//'>',   modname)
     426      CALL msg(  '   tnat('//TRIM(sxt)//') = '//TRIM(real2str(tnat(ixt))), modname)
     427!     CALL msg('   alpha_liq_sol('//TRIM(sxt)//') = '//TRIM(real2str(alpha_liq_sol(ixt))), modname)
     428!     CALL msg(       '   tkcin0('//TRIM(sxt)//') = '//TRIM(real2str(tkcin0(ixt))),        modname)
     429!     CALL msg(      '   tdifrel('//TRIM(sxt)//') = '//TRIM(real2str(tdifrel(ixt))),       modname)
     430   END DO
     431   CALL msg('69:     lambda = '//TRIM(real2str(lambda_sursat)), modname)
     432   CALL msg('69:    thumxt1 = '//TRIM(real2str(thumxt1)),       modname)
     433   CALL msg('69: h_land_ice = '//TRIM(real2str(h_land_ice)),    modname)
     434   CALL msg('69:      P_veg = '//TRIM(real2str(P_veg)),         modname)
     435
    685436END SUBROUTINE iso_init
    686437
     438
     439SUBROUTINE getinp_s(nam, val, def, lDisp)
     440   USE ioipsl_getincom, ONLY: getin
     441   USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
     442   USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
     443   USE mod_phys_lmdz_transfert_para, ONLY : bcast
     444   CHARACTER(LEN=*),  INTENT(IN)    :: nam
     445   CHARACTER(LEN=*),  INTENT(INOUT) :: val
     446   CHARACTER(LEN=*),  INTENT(IN)    :: def
     447   LOGICAL, OPTIONAL, INTENT(IN)    :: lDisp
     448   LOGICAL :: lD
     449!$OMP BARRIER
     450   IF(.NOT.(is_mpi_root.AND.is_omp_root)) RETURN
     451   val=def; CALL getin(nam,val); CALL bcast(val)
     452   lD=.TRUE.; IF(PRESENT(lDisp)) lD=lDisp
     453   IF(lD) CALL msg(TRIM(nam)//' = '//TRIM(val))
     454END SUBROUTINE getinp_s
     455
     456SUBROUTINE getinp_i(nam, val, def, lDisp)
     457   USE ioipsl_getincom, ONLY: getin
     458   USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
     459   USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
     460   USE mod_phys_lmdz_transfert_para, ONLY : bcast
     461   CHARACTER(LEN=*),  INTENT(IN)    :: nam
     462   INTEGER,           INTENT(INOUT) :: val
     463   INTEGER,           INTENT(IN)    :: def
     464   LOGICAL, OPTIONAL, INTENT(IN)    :: lDisp
     465   LOGICAL :: lD
     466!$OMP BARRIER
     467   IF(.NOT.(is_mpi_root.AND.is_omp_root)) RETURN
     468   val=def; CALL getin(nam,val); CALL bcast(val)
     469   lD=.TRUE.; IF(PRESENT(lDisp)) lD=lDisp
     470   IF(lD) CALL msg(TRIM(nam)//' = '//TRIM(int2str(val)))
     471END SUBROUTINE getinp_i
     472
     473SUBROUTINE getinp_r(nam, val, def, lDisp)
     474   USE ioipsl_getincom, ONLY: getin
     475   USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
     476   USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
     477   USE mod_phys_lmdz_transfert_para, ONLY : bcast
     478   CHARACTER(LEN=*),  INTENT(IN)    :: nam
     479   REAL,              INTENT(INOUT) :: val
     480   REAL,              INTENT(IN)    :: def
     481   LOGICAL, OPTIONAL, INTENT(IN)    :: lDisp
     482   LOGICAL :: lD
     483!$OMP BARRIER
     484   IF(.NOT.(is_mpi_root.AND.is_omp_root)) RETURN
     485   val=def; CALL getin(nam,val); CALL bcast(val)
     486   lD=.TRUE.; IF(PRESENT(lDisp)) lD=lDisp
     487   IF(lD) CALL msg(TRIM(nam)//' = '//TRIM(real2str(val)))
     488END SUBROUTINE getinp_r
     489
     490SUBROUTINE getinp_l(nam, val, def, lDisp)
     491   USE ioipsl_getincom, ONLY: getin
     492   USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
     493   USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
     494   USE mod_phys_lmdz_transfert_para, ONLY : bcast
     495   CHARACTER(LEN=*),  INTENT(IN)    :: nam
     496   LOGICAL,           INTENT(INOUT) :: val
     497   LOGICAL,           INTENT(IN)    :: def
     498   LOGICAL, OPTIONAL, INTENT(IN)    :: lDisp
     499   LOGICAL :: lD
     500!$OMP BARRIER
     501   IF(.NOT.(is_mpi_root.AND.is_omp_root)) RETURN
     502   val=def; CALL getin(nam,val); CALL bcast(val)
     503   lD=.TRUE.; IF(PRESENT(lDisp)) lD=lDisp
     504   IF(lD) CALL msg(TRIM(nam)//' = '//TRIM(bool2str(val)))
     505END SUBROUTINE getinp_l
    687506
    688507END MODULE isotopes_mod
Note: See TracChangeset for help on using the changeset viewer.