Ignore:
Timestamp:
Sep 10, 2024, 5:14:23 PM (5 months ago)
Author:
dcugnet
Message:
  • Remove INCA retro-compatibility with "traceur.def" (containing only water tracers but getting chemical species from an internal INCA routine).
  • The "trac_type" derived type internal to "readTracFiles_mod" is removed because a generic "keys_type" is enough: no explicit key ("%" operator) is needed, even %name.
  • The "trac_type" and "isot_type" derived types are now defined locally in "infotrac" and "infotrac_phy" (and more generally in each context: dynamic, lmdz dynamics, lmdz physics, etc.). The "readTracFiles_mod" module is now only used in these two routines:
    • few internal routines/variables (addPhase, delPhase, new2oldH2O, newHNO3, oldHNO3) are made available through "infotrac" and "infotrac_phy".
    • the "getKey" routine is only used in these two routines to define the explicit keys ("%" operator) of the local derived types "trac_type" and "isot_type". It could be in principle used outside this scope to get tracers parameters (read from "tracer.def") or isotopic parameters (read from "isotopes_params.def" - disabled for now).
  • The fortran parameters file "iso_params_mod.F90" is introduced so that "tnat" and "alpha_ideal" are defined in a single place but used in several. "ltnat1" is no longer hardcoded but defined with the *.def files parameter "tnat1"
  • Few minor changes:
    • use "infotrac_phy" instead of "infotrac" in calfis* because "tracers(:)%isAdvected" is defined in physics only.
    • "isotopes_mod" now ready for several isotopes classes (currently: only H2O)
    • isotopes class name (the name of the parent of the isotopes) is now %name and no longer %parent.
    • improvement of "getKey"
Location:
LMDZ6/trunk/libf/phylmdiso
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r4982 r5183  
    44MODULE isotopes_mod
    55   USE strings_mod,  ONLY: msg, real2str, int2str, bool2str, maxlen, strIdx, strStack
    6    USE infotrac_phy, ONLY: isoName
     6   USE infotrac_phy, ONLY: isoName, isoSelect, niso, ntiso, nbIso, isoFamilies
     7   USE iso_params_mod
    78   IMPLICIT NONE
    89   INTERFACE get_in; MODULE PROCEDURE getinp_s, getinp_i, getinp_r, getinp_l;  END INTERFACE get_in
    … …  
    1112  !--- Contains all isotopic variables + their initialization
    1213  !--- Isotopes-specific routines are in isotopes_routines_mod to avoid circular dependencies with isotopes_verif_mod.
     14
     15   !--- Negligible lower thresholds: no need to check for absurd values under these lower limits
     16   REAL, PARAMETER :: &
     17      ridicule      = 1e-12,              & ! For mixing ratios
     18      ridicule_rain = 1e-8,               & ! For rain fluxes (rain, zrfl...) in kg/s <-> 1e-3 mm/day
     19      ridicule_evap = ridicule_rain*1e-2, & ! For evaporations                in kg/s <-> 1e-3 mm/day
     20      ridicule_qsol = ridicule_rain,      & ! For qsol                        in kg <-> 1e-8 kg
     21      ridicule_snow = ridicule_qsol         ! For snow                        in kg <-> 1e-8 kg
     22   REAL, PARAMETER :: expb_max = 30.0
     23
     24   !--- Fractionation coefficients for H217O
     25   REAL, PARAMETER :: fac_coeff_eq17_liq = 0.529,    &
     26                      fac_coeff_eq17_ice = 0.529
     27
     28   !--- H218O reference
     29   REAL, PARAMETER :: fac_enrichoce18 = 0.0005,  alpha_liq_sol_O18 = 1.00291,                    &
     30                      talph1_O18 = 1137.,        talps1_O18 = 11.839,     tkcin0_O18 = 0.006,    &
     31                      talph2_O18 = -0.4156,      talps2_O18 = -0.028244,  tkcin1_O18 = 0.000285, &
     32                      talph3_O18 = -2.0667E-3,  tdifrel_O18 = 1./0.9723,  tkcin2_O18 = 0.00082
     33
     34   !--- Parameters that do not depend on the nature of water isotopes:
     35   REAL, PARAMETER :: pxtmelt = 273.15           !--- temperature at which ice formation starts
     36   REAL, PARAMETER :: pxtice  = 273.15 -  10.0   !--- temperature at which all condensate is ice:
     37   REAL, PARAMETER :: pxtmin  = 273.15 - 120.0   !--- computation done only under -120°C
     38   REAL, PARAMETER :: pxtmax  = 273.15 +  60.0   !--- computation done only  over  +60°C
     39   REAL, PARAMETER :: tdifexp = 0.58             !--- a constant for alpha_eff for equilibrium below cloud base:
     40   REAL, PARAMETER :: tv0cin  = 7.0              !--- wind threshold (m/s) for smooth/rough regime (Merlivat and Jouzel 1979)
     41   REAL, PARAMETER :: musi    = 1.0              !--- facteurs lambda et mu dans Si=musi-lambda*T
     42   REAL, PARAMETER :: Kd      = 2.5e-9    ! m2/s !--- diffusion in soil
     43   REAL, PARAMETER :: rh_cste_surf_cond = 0.6    !--- cste_surf_cond case: rhs and/or Ts set to constants
     44   REAL, PARAMETER :: T_cste_surf_cond = 288.0
    1345
    1446   !--- Isotopes indices (in [1,niso] ; non-existing => 0 index)
    … …  
    89121   !--- Vectors of length "niso"
    90122   REAL, ALLOCATABLE, DIMENSION(:), SAVE :: &
    91                     tnat, toce, tcorr, tdifrel
    92 !$OMP THREADPRIVATE(tnat, toce, tcorr, tdifrel)
     123                    alpha, tnat, toce, tcorr, tdifrel
     124!$OMP THREADPRIVATE(alpha, tnat, toce, tcorr, tdifrel)
    93125   REAL, ALLOCATABLE, DIMENSION(:), SAVE :: &
    94126                    talph1, talph2, talph3, talps1, talps2
    … …  
    100132                    alpha_liq_sol, Rdefault, Rmethox
    101133!$OMP THREADPRIVATE(alpha_liq_sol, Rdefault, Rmethox)
    102 !   REAL, SAVE ::    fac_coeff_eq17_liq, fac_coeff_eq17_ice
    103 !!$OMP THREADPRIVATE(fac_coeff_eq17_liq, fac_coeff_eq17_ice)
    104 
    105    !--- H2[18]O reference
    106    REAL, PARAMETER :: fac_enrichoce18=0.0005
    107    REAL, PARAMETER :: alpha_liq_sol_O18=1.00291
    108    REAL, PARAMETER :: talph1_O18=1137.
    109    REAL, PARAMETER :: talph2_O18=-0.4156
    110    REAL, PARAMETER :: talph3_O18=-2.0667E-3
    111    REAL, PARAMETER :: talps1_O18=11.839
    112    REAL, PARAMETER :: talps2_O18=-0.028244
    113    REAL, PARAMETER :: tdifrel_O18=1./0.9723
    114    REAL, PARAMETER :: tkcin0_O18=0.006
    115    REAL, PARAMETER :: tkcin1_O18=0.000285
    116    REAL, PARAMETER :: tkcin2_O18=0.00082
    117    REAL, PARAMETER :: fac_coeff_eq17_liq=0.529
    118    REAL, PARAMETER :: fac_coeff_eq17_ice=0.529
    119 
    120    !---- Parameters that do not depend on the nature of water isotopes:
    121    REAL, PARAMETER :: pxtmelt = 273.15 ! temperature at which ice formation starts
    122    REAL, PARAMETER :: pxtice  = 273.15-10.0 ! -- temperature at which all condensate is ice:
    123    REAL, PARAMETER :: pxtmin = 273.15 - 120.0   ! On ne calcule qu'au dessus de -120°C
    124    REAL, PARAMETER :: pxtmax = 273.15 +  60.0   ! On ne calcule qu'au dessus de +60°C
    125    REAL, PARAMETER :: tdifexp = 0.58 ! -- a constant for alpha_eff for equilibrium below cloud base:
    126    REAL, PARAMETER :: tv0cin  = 7.0 ! wind threshold (m/s) for smooth/rough regime (Merlivat and Jouzel 1979)
    127    REAL, PARAMETER :: musi=1.0  ! facteurs lambda et mu dans Si=musi-lambda*T
    128    REAL, PARAMETER :: Kd=2.5e-9 ! m2/s ! diffusion dans le sol
    129    REAL, PARAMETER :: rh_cste_surf_cond = 0.6 ! cas où cste_surf_cond: on met rhs ou/et Ts cste pour voir
    130    REAL, PARAMETER :: T_cste_surf_cond = 288.0
    131 
    132 
    133    !--- Negligible lower thresholds: no need to check for absurd values under these lower limits
    134    REAL, PARAMETER :: &
    135       ridicule      = 1e-12,              & ! For mixing ratios
    136       ridicule_rain = 1e-8,               & ! For rain fluxes (rain, zrfl...) in kg/s <-> 1e-3 mm/day
    137       ridicule_evap = ridicule_rain*1e-2, & ! For evaporations                in kg/s <-> 1e-3 mm/day
    138       ridicule_qsol = ridicule_rain,      & ! For qsol                        in kg <-> 1e-8 kg
    139       ridicule_snow = ridicule_qsol         ! For snow                        in kg <-> 1e-8 kg
    140    REAL, PARAMETER :: expb_max = 30.0
    141134
    142135   !--- Specific to HTO:
    … …  
    155148
    156149SUBROUTINE iso_init()
    157    USE infotrac_phy,       ONLY: ntiso, niso, getKey
    158     USE strings_mod,       ONLY: maxlen
    159150   IMPLICIT NONE
    160151
    161152   !=== Local variables:
    162    INTEGER :: ixt
    163 
     153   INTEGER :: ixt, ii, is
     154   LOGICAL :: ltnat1
     155   CHARACTER(LEN=maxlen) :: modname, sxt
    164156 
    165157   !--- For H2[17]O
    … …  
    170162   LOGICAL, PARAMETER ::   ok_nocinsat = .FALSE. ! if T: no sursaturation effect for ice
    171163   LOGICAL, PARAMETER :: Rdefault_smow = .FALSE. ! if T: Rdefault=smow; if F: nul
    172    LOGICAL, PARAMETER :: tnat1 = .TRUE. ! If T: all tnats are 1.
    173164
    174165   !--- For [3]H
    175166   INTEGER :: iessai
    176 
    177    CHARACTER(LEN=maxlen) :: modname, sxt
    178167
    179168   modname = 'iso_init'
    … …  
    187176   CALL msg('64: niso = '//TRIM(int2str(niso)), modname)
    188177
    189    !--- Init de ntracisoOR: on ecrasera en cas de traceurs de tagging isotopiques
    190    !                     (nzone>0) si complications avec ORCHIDEE
    191    ntracisoOR = ntiso 
    192 
    193    !--- Type of water isotopes:
    194    iso_eau = strIdx(isoName, 'H216O'); CALL msg('iso_eau='//int2str(iso_eau), modname)
    195    iso_HDO = strIdx(isoName, 'HDO');   CALL msg('iso_HDO='//int2str(iso_HDO), modname)
    196    iso_O18 = strIdx(isoName, 'H218O'); CALL msg('iso_O18='//int2str(iso_O18), modname)
    197    iso_O17 = strIdx(isoName, 'H217O'); CALL msg('iso_O17='//int2str(iso_O17), modname)
    198    iso_HTO = strIdx(isoName, 'HTO');   CALL msg('iso_HTO='//int2str(iso_HTO), modname)
    199 
    200    !--- Initialiaation: reading the isotopic parameters.
    201    CALL get_in('lambda',     lambda_sursat, 0.004)
    202    CALL get_in('thumxt1',    thumxt1,       0.75*1.2)
    203    CALL get_in('ntot',       ntot,          20,  .FALSE.)
    204    CALL get_in('h_land_ice', h_land_ice,    20., .FALSE.)
    205    CALL get_in('P_veg',      P_veg,         1.0, .FALSE.)
    206    CALL get_in('bidouille_anti_divergence', bidouille_anti_divergence, .FALSE.)
    207    CALL get_in('essai_convergence',         essai_convergence,         .FALSE.)
    208    CALL get_in('initialisation_iso',        initialisation_iso,        0)
    209 
    210 !  IF(nzone>0 .AND. initialisation_iso==0) &
    211 !      CALL get_in('initialisation_isotrac',initialisation_isotrac)
    212    CALL get_in('modif_sst',      modif_sst,         0)
    213    CALL get_in('deltaTtest',     deltaTtest,      0.0)     !--- For modif_sst>=1
    214    CALL get_in('deltaTtestpoles',deltaTtestpoles, 0.0)     !--- For modif_sst>=2
    215    CALL get_in( 'sstlatcrit',    sstlatcrit,     30.0)     !--- For modif_sst>=3
    216    CALL get_in('dsstlatcrit',   dsstlatcrit,      0.0)     !--- For modif_sst>=3
     178   DO ii = 1, nbIso
     179      CALL msg('Can''t select isotopes class "'//TRIM(isoFamilies(ii))//'"', modname, isoSelect(ii, lVerbose=.TRUE.))
     180
     181!==============================================================================================================================
     182      IF(isoFamilies(ii) == 'H2O') THEN
     183!==============================================================================================================================
     184         !--- Init de ntracisoOR: on ecrasera en cas de traceurs de tagging isotopiques
     185         !                     (nzone>0) si complications avec ORCHIDEE
     186         ntracisoOR = ntiso
     187
     188         !--- Type of water isotopes:
     189         iso_eau = strIdx(isoName, 'H216O'); CALL msg('iso_eau='//int2str(iso_eau), modname)
     190         iso_HDO = strIdx(isoName, 'HDO');   CALL msg('iso_HDO='//int2str(iso_HDO), modname)
     191         iso_O18 = strIdx(isoName, 'H218O'); CALL msg('iso_O18='//int2str(iso_O18), modname)
     192         iso_O17 = strIdx(isoName, 'H217O'); CALL msg('iso_O17='//int2str(iso_O17), modname)
     193         iso_HTO = strIdx(isoName, 'HTO');   CALL msg('iso_HTO='//int2str(iso_HTO), modname)
     194
     195         !--- Initialisation: reading the isotopic parameters.
     196         CALL get_in('lambda',     lambda_sursat, 0.004); IF(ok_nocinsat) lambda_sursat = 0.
     197         CALL get_in('thumxt1',    thumxt1,       0.75*1.2)
     198         CALL get_in('ntot',       ntot,          20,  .FALSE.)
     199         CALL get_in('h_land_ice', h_land_ice,    20., .FALSE.)
     200         CALL get_in('P_veg',      P_veg,         1.0, .FALSE.)
     201         CALL get_in('bidouille_anti_divergence', bidouille_anti_divergence, .FALSE.)
     202         CALL get_in('essai_convergence',         essai_convergence,         .FALSE.)
     203         CALL get_in('initialisation_iso',        initialisation_iso,        0)
     204
     205!        IF(nzone>0 .AND. initialisation_iso==0) &
     206!           CALL get_in('initialisation_isotrac',initialisation_isotrac)
     207         CALL get_in('modif_sst',      modif_sst,         0)
     208         CALL get_in('deltaTtest',     deltaTtest,      0.0)     !--- For modif_sst>=1
     209         CALL get_in('deltaTtestpoles',deltaTtestpoles, 0.0)     !--- For modif_sst>=2
     210         CALL get_in( 'sstlatcrit',    sstlatcrit,     30.0)     !--- For modif_sst>=3
     211         CALL get_in('dsstlatcrit',   dsstlatcrit,      0.0)     !--- For modif_sst>=3
    217212#ifdef ISOVERIF
    218    CALL msg('iso_init 270:  sstlatcrit='//real2str( sstlatcrit), modname, sstlatcrit < 0.0) !--- For modif_sst>=2
    219    CALL msg('iso_init 279: dsstlatcrit='//real2str(dsstlatcrit), modname, sstlatcrit < 0.0) !--- For modif_sst>=3
    220    IF(modif_sst >= 2 .AND. sstlatcrit < 0.0) STOP
     213         CALL msg('iso_init 270:  sstlatcrit='//real2str( sstlatcrit), modname, sstlatcrit < 0.0) !--- For modif_sst>=2
     214         CALL msg('iso_init 279: dsstlatcrit='//real2str(dsstlatcrit), modname, sstlatcrit < 0.0) !--- For modif_sst>=3
     215         IF(modif_sst >= 2 .AND. sstlatcrit < 0.0) STOP
    221216#endif             
    222217
    223    CALL get_in('modif_sic', modif_sic,  0)
    224    IF(modif_sic >= 1) &
    225    CALL get_in('deltasic',  deltasic, 0.1)
    226 
    227    CALL get_in('albedo_prescrit', albedo_prescrit, 0)
    228    IF(albedo_prescrit == 1) THEN
    229       CALL get_in('lon_min_albedo', lon_min_albedo, -200.)
    230       CALL get_in('lon_max_albedo', lon_max_albedo,  200.)
    231       CALL get_in('lat_min_albedo', lat_min_albedo, -100.)
    232       CALL get_in('lat_max_albedo', lat_max_albedo,  100.)
    233    END IF
    234    CALL get_in('deltaP_BL',           deltaP_BL,     10.0)
    235    CALL get_in('ruissellement_pluie', ruissellement_pluie, 0)
    236    CALL get_in('alphak_stewart',      alphak_stewart,      1)
    237    CALL get_in('tdifexp_sol',         tdifexp_sol,      0.67)
    238    CALL get_in('calendrier_guide',    calendrier_guide,    0)
    239    CALL get_in('cste_surf_cond',      cste_surf_cond,      0)
    240    CALL get_in('mixlen',              mixlen,           35.0)
    241    CALL get_in('evap_cont_cste',      evap_cont_cste,      0)
    242    CALL get_in('deltaO18_evap_cont',  deltaO18_evap_cont,0.0)
    243    CALL get_in('d_evap_cont',         d_evap_cont,       0.0)
    244    CALL get_in('nudge_qsol',          nudge_qsol,          0)
    245    CALL get_in('region_nudge_qsol',   region_nudge_qsol,   1)
    246    nlevmaxO17 = 50
    247    CALL msg('nlevmaxO17='//TRIM(int2str(nlevmaxO17)))
    248    CALL get_in('no_pce',   no_pce,     0)
    249    CALL get_in('A_satlim', A_satlim, 1.0)
    250    CALL get_in('ok_restrict_A_satlim', ok_restrict_A_satlim, 0)
     218         CALL get_in('modif_sic', modif_sic,  0)
     219         IF(modif_sic >= 1) &
     220         CALL get_in('deltasic',  deltasic, 0.1)
     221
     222         CALL get_in('albedo_prescrit', albedo_prescrit, 0)
     223         IF(albedo_prescrit == 1) THEN
     224            CALL get_in('lon_min_albedo', lon_min_albedo, -200.)
     225            CALL get_in('lon_max_albedo', lon_max_albedo,  200.)
     226            CALL get_in('lat_min_albedo', lat_min_albedo, -100.)
     227            CALL get_in('lat_max_albedo', lat_max_albedo,  100.)
     228         END IF
     229         CALL get_in('deltaO18_oce',        deltaO18_oce,   0.0)
     230         CALL get_in('deltaP_BL',           deltaP_BL,     10.0)
     231         CALL get_in('ruissellement_pluie', ruissellement_pluie, 0)
     232         CALL get_in('alphak_stewart',      alphak_stewart,      1)
     233         CALL get_in('tdifexp_sol',         tdifexp_sol,      0.67)
     234         CALL get_in('calendrier_guide',    calendrier_guide,    0)
     235         CALL get_in('cste_surf_cond',      cste_surf_cond,      0)
     236         CALL get_in('mixlen',              mixlen,           35.0)
     237         CALL get_in('evap_cont_cste',      evap_cont_cste,      0)
     238         CALL get_in('deltaO18_evap_cont',  deltaO18_evap_cont,0.0)
     239         CALL get_in('d_evap_cont',         d_evap_cont,       0.0)
     240         CALL get_in('nudge_qsol',          nudge_qsol,          0)
     241         CALL get_in('region_nudge_qsol',   region_nudge_qsol,   1)
     242         nlevmaxO17 = 50
     243         CALL msg('nlevmaxO17='//TRIM(int2str(nlevmaxO17)))
     244         CALL get_in('no_pce',   no_pce,     0)
     245         CALL get_in('A_satlim', A_satlim, 1.0)
     246         CALL get_in('ok_restrict_A_satlim', ok_restrict_A_satlim, 0)
    251247#ifdef ISOVERIF
    252    CALL msg(' 315: A_satlim='//real2str(A_satlim), modname, A_satlim > 1.0)
    253    IF(A_satlim > 1.0) STOP
     248         CALL msg(' 315: A_satlim='//real2str(A_satlim), modname, A_satlim > 1.0)
     249         IF(A_satlim > 1.0) STOP
    254250#endif
    255 !  CALL get_in('slope_limiterxy',   slope_limiterxy,  2.0)
    256 !  CALL get_in('slope_limiterz',    slope_limiterz,   2.0)
    257    CALL get_in('modif_ratqs',       modif_ratqs,        0)
    258    CALL get_in('Pcrit_ratqs',       Pcrit_ratqs,    500.0)
    259    CALL get_in('ratqsbasnew',       ratqsbasnew,     0.05)
    260    CALL get_in('fac_modif_evaoce',  fac_modif_evaoce, 1.0)
    261    CALL get_in('ok_bidouille_wake', ok_bidouille_wake,  0)
    262    ! si oui, la temperature de cond est celle de l'environnement, pour eviter
    263    ! bugs quand temperature dans ascendances convs est mal calculee
    264    CALL get_in('cond_temp_env',        cond_temp_env,        .FALSE.)
    265    IF(ANY(isoName == 'HTO')) &
    266    CALL get_in('ok_prod_nucl_tritium', ok_prod_nucl_tritium, .FALSE., .FALSE.)
    267 
    268    ! Ocean composition
    269    CALL get_in('deltaO18_oce',  deltaO18_oce, 0.0)
    270    
    271    CALL msg('iso_O18, iso_HDO, iso_eau = '//TRIM(strStack(int2str([iso_O18, iso_HDO, iso_eau]))), modname)
    272 
    273    !--------------------------------------------------------------
    274    ! Isotope fractionation factors and a few isotopic constants
    275    !--------------------------------------------------------------
    276    ALLOCATE(tkcin0(niso))
    277    ALLOCATE(tkcin1(niso))
    278    ALLOCATE(tkcin2(niso))
    279    ALLOCATE(tnat(niso))
    280    ALLOCATE(tdifrel(niso))
    281    ALLOCATE(toce(niso))
    282    ALLOCATE(tcorr(niso))
    283    ALLOCATE(talph1(niso))
    284    ALLOCATE(talph2(niso))
    285    ALLOCATE(talph3(niso))
    286    ALLOCATE(talps1(niso))
    287    ALLOCATE(talps2(niso))
    288    ALLOCATE(alpha_liq_sol(niso))
    289    ALLOCATE(Rdefault(niso))
    290    ALLOCATE(Rmethox(niso))
    291 
    292    do ixt=1,niso
    293      if (ixt.eq.iso_HTO) then  ! Tritium
    294        tkcin0(ixt) = 0.01056
    295        tkcin1(ixt) = 0.0005016
    296        tkcin2(ixt) = 0.0014432
    297        if (tnat1) then
    298                tnat(ixt)=1
    299        else
    300                tnat(ixt)=0.
    301        endif
    302        toce(ixt)=4.0E-19 ! rapport T/H = 0.2 TU Dreisigacker and Roether 1978
    303        tcorr(ixt)=1.
    304        tdifrel(ixt)=1./0.968
    305        talph1(ixt)=46480.
    306        talph2(ixt)=-103.87
    307        talph3(ixt)=0.
    308        talps1(ixt)=46480.
    309        talps2(ixt)=-103.87
    310        alpha_liq_sol(ixt)=1.
    311        Rmethox(ixt)=0.0
    312      endif
    313      if (ixt.eq.iso_O17) then  ! O17
    314        pente_MWL=0.528
    315        tdifrel(ixt)=1./0.98555 ! valeur utilisée en 1D et dans modèle de LdG ! tdifrel(ixt)=1./0.985452 ! donné par Amaelle
    316        fac_kcin= (tdifrel(ixt)-1.0)/(tdifrel_O18-1.0) ! fac_kcin=0.5145 ! donné par Amaelle
    317        tkcin0(ixt) = tkcin0_O18*fac_kcin
    318        tkcin1(ixt) = tkcin1_O18*fac_kcin
    319        tkcin2(ixt) = tkcin2_O18*fac_kcin
    320        if (tnat1) then
    321                tnat(ixt)=1
    322        else
    323                tnat(ixt)=0.004/100. ! O17 représente 0.004% de l'oxygène
    324        endif
    325        toce(ixt)=tnat(ixt)*(1.0+deltaO18_oce/1000.0)**pente_MWL
    326        tcorr(ixt)=1.0+fac_enrichoce18*pente_MWL ! donné par Amaelle           
    327        talph1(ixt)=talph1_O18
    328        talph2(ixt)=talph2_O18
    329        talph3(ixt)=talph3_O18
    330        talps1(ixt)=talps1_O18
    331        talps2(ixt)=talps2_O18     
    332        alpha_liq_sol(ixt)=(alpha_liq_sol_O18)**fac_coeff_eq17_liq
    333        Rdefault(ixt)=tnat(ixt)*(-3.15/1000.0+1.0)
    334        Rmethox(ixt)=(230./1000.+1.)*tnat(ixt) !Zahn et al 2006
    335      endif
    336      if (ixt.eq.iso_O18) then  ! Oxygene18
    337        tkcin0(ixt) = tkcin0_O18
    338        tkcin1(ixt) = tkcin1_O18
    339        tkcin2(ixt) = tkcin2_O18
    340        if (tnat1) then
    341                tnat(ixt)=1
    342        else
    343                tnat(ixt)=2005.2E-6
    344        endif
    345        toce(ixt)=tnat(ixt)*(1.0+deltaO18_oce/1000.0)
    346        tcorr(ixt)=1.0+fac_enrichoce18
    347        tdifrel(ixt)=tdifrel_O18
    348        talph1(ixt)=talph1_O18
    349        talph2(ixt)=talph2_O18
    350        talph3(ixt)=talph3_O18
    351        talps1(ixt)=talps1_O18
    352        talps2(ixt)=talps2_O18
    353        alpha_liq_sol(ixt)=alpha_liq_sol_O18   
    354        Rdefault(ixt)=tnat(ixt)*(-6.0/1000.0+1.0)
    355        Rmethox(ixt)=(130./1000.+1.)*tnat(ixt) !Zahn et al 2006 
    356      endif
    357      if (ixt.eq.iso_HDO) then ! Deuterium
    358        pente_MWL=8.0
    359        tdifrel(ixt)=1./0.9755 !          fac_kcin=0.88
    360        fac_kcin= (tdifrel(ixt)-1)/(tdifrel_O18-1)
    361        tkcin0(ixt) = tkcin0_O18*fac_kcin
    362        tkcin1(ixt) = tkcin1_O18*fac_kcin
    363        tkcin2(ixt) = tkcin2_O18*fac_kcin
    364        if (tnat1) then
    365                tnat(ixt)=1
    366        else
    367                tnat(ixt)=155.76E-6
    368        endif
    369        toce(ixt)=tnat(ixt)*(1.0+pente_MWL*deltaO18_oce/1000.0)
    370        tcorr(ixt)=1.0+fac_enrichoce18*pente_MWL         
    371        talph1(ixt)=24844.
    372        talph2(ixt)=-76.248
    373        talph3(ixt)=52.612E-3
    374        talps1(ixt)=16288.
    375        talps2(ixt)=-0.0934
    376        !ZXalpha_liq_sol=1.0192 ! Weston, Ralph, 1955
    377        alpha_liq_sol(ixt)=1.0212
    378        ! valeur de Lehmann & Siegenthaler, 1991, Journal of
    379        ! Glaciology, vol 37, p 23
    380        Rdefault(ixt)=tnat(ixt)*((-6.0*pente_MWL+10.0)/1000.0+1.0)
    381        Rmethox(ixt)=tnat(ixt)*(-25.0/1000.+1.) ! Zahn et al 2006
    382      endif
    383      if (ixt.eq.iso_eau) then ! Oxygene16
    384        tkcin0(ixt) = 0.0
    385        tkcin1(ixt) = 0.0
    386        tkcin2(ixt) = 0.0
    387        tnat(ixt)=1.
    388        toce(ixt)=tnat(ixt)
    389        tcorr(ixt)=1.0
    390        tdifrel(ixt)=1.
    391        talph1(ixt)=0.
    392        talph2(ixt)=0.
    393        talph3(ixt)=0.
    394        talps1(ixt)=0.
    395        talph3(ixt)=0.
    396        alpha_liq_sol(ixt)=1.
    397        Rdefault(ixt)=tnat(ixt)*1.0
    398        Rmethox(ixt)=1.0
    399      endif
    400    enddo ! ixt=1,niso
    401 
    402    IF(.NOT.Rdefault_smow) then
    403         Rdefault(:) = 0.0
    404         if (iso_eau.gt.0) Rdefault(iso_eau) = 1.0 ! correction Camille 30 mars 2023
    405    ENDIF
    406    write(*,*) 'Rdefault=',Rdefault
    407    write(*,*) 'toce=',toce
    408 
    409    !--- Sensitivity test: no kinetic effect in sfc evaporation
    410    IF(ok_nocinsfc) THEN
    411       tkcin0(1:niso) = 0.0
    412       tkcin1(1:niso) = 0.0
    413       tkcin2(1:niso) = 0.0
    414    END IF
    415 
    416    CALL msg('285: verif initialisation:', modname)
    417    DO ixt=1,niso
    418       sxt=int2str(ixt)
    419       CALL msg(' * isoName('//TRIM(sxt)//') = <'//TRIM(isoName(ixt))//'>',  modname)
    420       CALL msg(  '    tnat('//TRIM(sxt)//') = '//TRIM(real2str(tnat(ixt))), modname)
    421 !     CALL msg('    alpha_liq_sol('//TRIM(sxt)//') = '//TRIM(real2str(alpha_liq_sol(ixt))), modname)
    422 !     CALL msg(        '   tkcin0('//TRIM(sxt)//') = '//TRIM(real2str(tkcin0(ixt))),        modname)
    423 !     CALL msg(       '   tdifrel('//TRIM(sxt)//') = '//TRIM(real2str(tdifrel(ixt))),       modname)
     251!        CALL get_in('slope_limiterxy',   slope_limiterxy,  2.0)
     252!        CALL get_in('slope_limiterz',    slope_limiterz,   2.0)
     253         CALL get_in('modif_ratqs',       modif_ratqs,        0)
     254         CALL get_in('Pcrit_ratqs',       Pcrit_ratqs,    500.0)
     255         CALL get_in('ratqsbasnew',       ratqsbasnew,     0.05)
     256         CALL get_in('fac_modif_evaoce',  fac_modif_evaoce, 1.0)
     257         CALL get_in('ok_bidouille_wake', ok_bidouille_wake,  0)
     258         ! si oui, la temperature de cond est celle de l'environnement, pour eviter
     259         ! bugs quand temperature dans ascendances convs est mal calculee
     260         CALL get_in('cond_temp_env',        cond_temp_env,        .FALSE.)
     261         IF(ANY(isoName == 'HTO')) &
     262         CALL get_in('ok_prod_nucl_tritium', ok_prod_nucl_tritium, .FALSE., .FALSE.)
     263         CALL get_in('tnateq1', ltnat1, .TRUE.)
     264
     265         CALL msg('iso_O18, iso_HDO, iso_eau = '//TRIM(strStack(int2str([iso_O18, iso_HDO, iso_eau]))), modname)
     266
     267         !--------------------------------------------------------------
     268         ! Parameters that depend on the nature of water isotopes:
     269         !--------------------------------------------------------------
     270         ALLOCATE(tnat (niso), talph1(niso),  talps1(niso), tkcin0(niso), tdifrel (niso), alpha        (niso))
     271         ALLOCATE(toce (niso), talph2(niso),  talps2(niso), tkcin1(niso), Rdefault(niso), alpha_liq_sol(niso))
     272         ALLOCATE(tcorr(niso), talph3(niso),                tkcin2(niso), Rmethox (niso))
     273
     274         !=== H216O
     275         is = iso_eau
     276         IF(is /= 0) THEN
     277            tdifrel (is) = 1.0
     278            alpha   (is) = alpha_ideal_H216O
     279            tnat    (is) = tnat_H216O; IF(ltnat1) tnat(is) = 1.0
     280            toce    (is) = tnat(is)
     281            tcorr   (is) = 1.0
     282            talph1  (is) = 0.0;  talps1(is) = 0.0;  tkcin0(is)  = 0.0
     283            talph2  (is) = 0.0;  talps2(is) = 0.0;  tkcin1(is)  = 0.0
     284            talph3  (is) = 0.0;                     tkcin2(is)  = 0.0
     285            Rdefault(is) = tnat(is)*1.0
     286            Rmethox (is) = 1.0
     287            alpha_liq_sol(is) = 1.0
     288         END IF
     289
     290         !=== H217O
     291         is = iso_O17
     292         IF(is /= 0) THEN; pente_MWL = 0.528
     293            tdifrel (is) = 1./0.98555  ! used in 1D and in LdG model ; tdifrel=1./0.985452: from Amaelle
     294            alpha   (is) = alpha_ideal_H217O
     295            tnat    (is) = tnat_H217O; IF(ltnat1) tnat(is) = 1.0
     296            toce    (is) = tnat(is)*(1.0+deltaO18_oce/1000.0)**pente_MWL
     297            tcorr   (is) = 1.0+fac_enrichoce18*pente_MWL
     298            fac_kcin = (tdifrel(is)-1.0)/(tdifrel_O18-1.0)           ! fac_kcin=0.5145:     from Amaelle
     299            talph1  (is) = talph1_O18;  talps1(is) = talps1_O18;  tkcin0(is) = tkcin0_O18*fac_kcin
     300            talph2  (is) = talph2_O18;  talps2(is) = talps2_O18;  tkcin1(is) = tkcin1_O18*fac_kcin
     301            talph3  (is) = talph3_O18;                            tkcin2(is) = tkcin2_O18*fac_kcin
     302            Rdefault(is) = tnat(is)*(1.0-3.15/1000.)
     303            Rmethox (is) = tnat(is)*(1.0+230./1000.)
     304            alpha_liq_sol(is) = alpha_liq_sol_O18**fac_coeff_eq17_liq
     305         END IF
     306
     307         !=== H218O
     308         is = iso_O18
     309         IF(is /= 0) THEN
     310            tdifrel (is) = tdifrel_O18
     311            alpha   (is) = alpha_ideal_H218O
     312            tnat    (is) = tnat_H218O; IF(ltnat1) tnat(is) = 1.0
     313            toce    (is) = tnat(is)*(1.0+deltaO18_oce/1000.0)
     314            tcorr   (is) = 1.0+fac_enrichoce18
     315            talph1  (is) = talph1_O18;  talps1(is) = talps1_O18;  tkcin0(is) = tkcin0_O18
     316            talph2  (is) = talph2_O18;  talps2(is) = talps2_O18;  tkcin1(is) = tkcin1_O18
     317            talph3  (is) = talph3_O18;                            tkcin2(is) = tkcin2_O18
     318            Rdefault(is) = tnat(is)*(1.0-6.00/1000.)
     319            Rmethox (is) = tnat(is)*(1.0+130./1000.)  ! Zahn & al. 2006
     320            alpha_liq_sol(is) = alpha_liq_sol_O18
     321         END IF
     322
     323         !=== HDO
     324         is = iso_HDO
     325         IF(is /= 0) THEN; pente_MWL = 8.0
     326            tdifrel (is) = 1./0.9755                  ! fac_kcin=0.88
     327            alpha   (is) = alpha_ideal_HDO
     328            tnat    (is) = tnat_HDO; IF(ltnat1) tnat(is) = 1.0
     329            toce    (is) = tnat(is)*(1.0+deltaO18_oce/1000.0*pente_MWL)
     330            tcorr   (is) = 1.0+fac_enrichoce18*pente_MWL
     331            fac_kcin = (tdifrel(is)-1.0)/(tdifrel_O18-1.0)
     332            talph1  (is) = 24844.;      talps1(is) = 16288.;      tkcin0(is) = tkcin0_O18*fac_kcin
     333            talph2  (is) = -76.248;     talps2(is) = -0.0934;     tkcin1(is) = tkcin1_O18*fac_kcin
     334            talph3  (is) = 52.612E-3;                             tkcin2(is) = tkcin2_O18*fac_kcin
     335            Rdefault(is) = tnat(is)*(1.0+(10.0-6.0*pente_MWL)/1000.)
     336            Rmethox (is) = tnat(is)*(1.0-25.0/1000.)
     337            alpha_liq_sol(is) = 1.0212      ! Lehmann & Siegenthaler, 1991, Jo. of Glaciology, vol 37, p 23
     338                                            ! alpha_liq_sol=1.0192: Weston, Ralph, 1955
     339         END IF
     340
     341         !=== HTO
     342         is = iso_HTO
     343         IF(is /= 0) THEN
     344            tdifrel (is) = 1./0.968
     345            alpha   (is) = alpha_ideal_HTO
     346            tnat    (is) = tnat_HTO; IF(ltnat1) tnat(is) = 1.0
     347            toce    (is) = 4.0E-19          ! ratio T/H = 0.2 TU Dreisigacker & Roether 1978
     348            tcorr   (is) = 1.0
     349            talph1  (is) = 46480.;      talps1(is) = 46480.;      tkcin0(is) = 0.01056
     350            talph2  (is) = -103.87;     talps2(is) = -103.87;     tkcin1(is) = 0.0005016
     351            talph3  (is) = 0.0;                                   tkcin2(is) = 0.0014432
     352            Rdefault(is) = 0.0
     353            Rmethox (is) = 0.0
     354            alpha_liq_sol(is) = 1.0
     355         END IF
     356
     357         IF(.NOT. Rdefault_smow) THEN
     358            Rdefault(:) = 0.0; IF(iso_eau > 0) Rdefault(iso_eau) = 1.0
     359         END IF
     360         WRITE(*,*) 'Rdefault = ',Rdefault
     361         WRITE(*,*) 'toce = ', toce
     362
     363         !--- Sensitivity test: no kinetic effect in sfc evaporation
     364         IF(ok_nocinsfc) THEN
     365            tkcin0(1:niso) = 0.0
     366            tkcin1(1:niso) = 0.0
     367            tkcin2(1:niso) = 0.0
     368         END IF
     369
     370         CALL msg('285: verif initialisation:', modname)
     371         DO ixt=1,niso
     372            sxt=int2str(ixt)
     373            CALL msg(' * isoName('//TRIM(sxt)//') = <'//TRIM(isoName(ixt))//'>',  modname)
     374            CALL msg(  '    tnat('//TRIM(sxt)//') = '//TRIM(real2str(tnat(ixt))), modname)
     375!           CALL msg('    alpha_liq_sol('//TRIM(sxt)//') = '//TRIM(real2str(alpha_liq_sol(ixt))), modname)
     376!           CALL msg(        '   tkcin0('//TRIM(sxt)//') = '//TRIM(real2str(tkcin0(ixt))),        modname)
     377!           CALL msg(       '   tdifrel('//TRIM(sxt)//') = '//TRIM(real2str(tdifrel(ixt))),       modname)
     378         END DO
     379         CALL msg('69:     lambda = '//TRIM(real2str(lambda_sursat)), modname)
     380         CALL msg('69:    thumxt1 = '//TRIM(real2str(thumxt1)),       modname)
     381         CALL msg('69: h_land_ice = '//TRIM(real2str(h_land_ice)),    modname)
     382         CALL msg('69:      P_veg = '//TRIM(real2str(P_veg)),         modname)
     383!==============================================================================================================================
     384      ELSE
     385!==============================================================================================================================
     386         CALL abort_physic('"isotopes_mod" is not set up yet for isotopes family "'//TRIM(isoFamilies(ii))//'"', modname, 1)
     387!==============================================================================================================================
     388      END IF
     389!==============================================================================================================================
    424390   END DO
    425    CALL msg('69:     lambda = '//TRIM(real2str(lambda_sursat)), modname)
    426    CALL msg('69:    thumxt1 = '//TRIM(real2str(thumxt1)),       modname)
    427    CALL msg('69: h_land_ice = '//TRIM(real2str(h_land_ice)),    modname)
    428    CALL msg('69:      P_veg = '//TRIM(real2str(P_veg)),         modname)
    429391
    430392END SUBROUTINE iso_init

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

    r5084 r5183  
    1641916419   USE isotopes_mod,      ONLY: isoName,iso_HDO,iso_eau
    1642016420   USE phyetat0_get_mod,  ONLY: phyetat0_get, phyetat0_srf
    16421    USE readTracFiles_mod, ONLY: new2oldH2O
    16422    USE strings_mod,       ONLY: strIdx, strTail, maxlen, msg, int2str
     16421   USE infotrac_phy,      ONLY: new2oldH2O
     16422   USE strings_mod,       ONLY: strIdx, strHead, strTail, maxlen, msg, int2str
    1642316423#ifdef ISOVERIF
    1642416424   USE isotopes_verif_mod
    … …  
    1645916459      outiso = isoName(ixt)
    1646016460      oldIso = strTail(new2oldH2O(outiso), '_')            !--- Remove "H2O_" from "H2O_<iso>[_<tag>]"
    16461       i = INDEX(outiso, '_', .TRUE.)
    16462       oldIso2 = outiso(1:i-1)//outiso(i+1:LEN_TRIM(outiso)) ! CR 2023: on ajoute cette possibilité aussi, elle correspond au cas le plus récent.
     16461      oldIso2= TRIM(strHead(outiso,'_'))//strTail(outiso,'_') ! CR 2023: most recent possibility
    1646316462!      write(*,*) 'tmp 16541:'
    1646416463!      write(*,*) 'outiso=',outiso

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

    r4493 r5183  
    33
    44MODULE isotrac_mod
    5   USE infotrac_phy,      ONLY: niso, ntiso, nzone
    6   USE readTracFiles_mod, ONLY: delPhase
    7   USE isotopes_mod,      ONLY: ridicule, get_in
     5  USE infotrac_phy, ONLY: niso, ntiso, nzone, delPhase
     6  USE isotopes_mod, ONLY: ridicule, get_in
    87
    98  IMPLICIT NONE

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

    r5084 r5183  
    4040  USE geometry_mod,     ONLY: longitude_deg, latitude_deg
    4141  USE iostart,          ONLY: close_startphy, get_field, get_var, open_startphy
    42   USE infotrac_phy,     ONLY: nqtot, nbtr, type_trac, tracers
    43   USE readTracFiles_mod,ONLY: maxlen, new2oldH2O
     42  USE infotrac_phy,     ONLY: nqtot, nbtr, type_trac, tracers, new2oldH2O
     43  USE strings_mod,      ONLY: maxlen
    4444  USE traclmdz_mod,     ONLY: traclmdz_from_restart
    4545  USE carbon_cycle_mod, ONLY: carbon_cycle_init, carbon_cycle_cpl, carbon_cycle_tr, carbon_cycle_rad, co2_send, RCO2_glo

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

    r5169 r5183  
    3939    USE ioipsl_getin_p_mod, ONLY : getin_p
    4040    USE indice_sol_mod
    41     USE infotrac_phy, ONLY: nqtot, nbtr, nqo, tracers, type_trac,ivap,iliq,isol
    42     USE readTracFiles_mod, ONLY: addPhase
     41    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, tracers, type_trac,addPhase, ivap, iliq, isol
    4342    USE strings_mod,  ONLY: strIdx
    4443    USE iophy
Note: See TracChangeset for help on using the changeset viewer.