Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/iniacademic.F90 =================================================================== --- LMDZ6/branches/Amaury_dev/libf/dyn3d/iniacademic.F90 (revision 5228) +++ LMDZ6/branches/Amaury_dev/libf/dyn3d/iniacademic.F90 (revision 5229) @@ -1,8 +1,6 @@ -! $Id$ - SUBROUTINE iniacademic(vcov, ucov, teta, q, masse, ps, phis, time_0) USE lmdz_filtreg, ONLY: inifilr - USE lmdz_infotrac, ONLY: nqtot, niso, iqIsoPha, tracers, getKey, isoName, addPhase + USE lmdz_infotrac, ONLY: nqtot, niso, iqIsoPha, tracers, isoName, addPhase USE control_mod, ONLY: day_step, planet_type USE exner_hyb_m, ONLY: exner_hyb @@ -20,4 +18,5 @@ USE lmdz_academic, ONLY: tetarappel, knewt_t, kfrict, knewt_g, clat4 USE lmdz_comgeom + USE iso_params_mod ! tnat_* and alpha_ideal_* ! Author: Frederic Hourdin original: 15/01/93 @@ -323,6 +322,13 @@ WRITE(lunout, *)'In '//TRIM(modname)//': !!! Beware: alpha_ideal put to 1 !!!' ELSE - IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) & - CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1) + SELECT CASE(isoName(iName)) + CASE('H216O'); tnat = tnat_H216O; alpha_ideal = alpha_ideal_H216O + CASE('H217O'); tnat = tnat_H217O; alpha_ideal = alpha_ideal_H217O + CASE('H218O'); tnat = tnat_H218O; alpha_ideal = alpha_ideal_H218O + CASE('HDO'); tnat = tnat_HDO; alpha_ideal = alpha_ideal_HDO + CASE('HTO'); tnat = tnat_HTO; alpha_ideal = alpha_ideal_HTO + CASE DEFAULT + CALL abort_gcm(TRIM(modname),'unknown isotope "'//TRIM(isoName(iName))//'" ; check tracer.def file',1) + END SELECT END IF q(:,:,iq) = q(:,:,iqParent)*tnat*(q(:,:,iqParent)/30.e-3)**(alpha_ideal-1.) Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/lmdz_check_isotopes.f90 =================================================================== --- LMDZ6/branches/Amaury_dev/libf/dyn3d/lmdz_check_isotopes.f90 (revision 5228) +++ LMDZ6/branches/Amaury_dev/libf/dyn3d/lmdz_check_isotopes.f90 (revision 5229) @@ -9,5 +9,6 @@ USE lmdz_strings, ONLY: maxlen, msg, strIdx, strStack, int2str, real2str USE lmdz_infotrac, ONLY: nqtot, niso, nphas, isotope, isoCheck, iqIsoPha, isoSelect, & - ntiso, iH2O, nzone, tracers, isoName, itZonIso, getKey + ntiso, iH2O, nzone, tracers, isoName, itZonIso + USE iso_params_mod, ONLY: tnat_H216O, tnat_H217O, tnat_H218O, tnat_HDO, tnat_HTO USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm @@ -31,5 +32,6 @@ ridicule = 1e-12 INTEGER, SAVE :: iso_eau, iso_O17, iso_O18, iso_HDO, iso_HTO - LOGICAL, SAVE :: ltnat1, first = .TRUE. + LOGICAL :: ltnat1 + LOGICAL, SAVE :: first = .TRUE. modname = 'check_isotopes' @@ -40,14 +42,10 @@ ltnat1 = .TRUE.; CALL getin('tnateq1', ltnat1) ALLOCATE(tnat(niso)) - iso_eau = strIdx(isoName, 'H216O') - iso_O17 = strIdx(isoName, 'H217O') - iso_O18 = strIdx(isoName, 'H218O') - iso_HDO = strIdx(isoName, 'HDO') - iso_HTO = strIdx(isoName, 'HTO') - IF(ltnat1) THEN - tnat(:) = 1.0 - ELSE - IF(getKey('tnat', tnat)) CALL abort_gcm(modname, 'missing isotopic parameter', 1) - END IF + iso_eau = strIdx(isoName, 'H216O'); IF(iso_eau /= 0) tnat(iso_eau) = tnat_H216O + iso_O17 = strIdx(isoName, 'H217O'); IF(iso_O17 /= 0) tnat(iso_O17) = tnat_H217O + iso_O18 = strIdx(isoName, 'H218O'); IF(iso_O18 /= 0) tnat(iso_O18) = tnat_H218O + iso_HDO = strIdx(isoName, 'HDO'); IF(iso_HDO /= 0) tnat(iso_HDO) = tnat_HDO + iso_HTO = strIdx(isoName, 'HTO'); IF(iso_HTO /= 0) tnat(iso_HTO) = tnat_HTO + IF(ltnat1) tnat(:) = 1. first = .FALSE. END IF @@ -61,5 +59,5 @@ DO k = 1, llm DO i = 1, ip1jmp1 - IF(ABS(q(i, k, iq)) <= borne) CYCLE + IF(ABS(q(i, k, iq)) < borne) CYCLE WRITE(msg1, '(s,"(",i0,",",i0,",",i0,") = ",ES12.4)')TRIM(isoName(ixt)), i, k, iq, q(i, k, iq) CALL msg(msg1, modname) Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/lmdz_dynetat0.f90 =================================================================== --- LMDZ6/branches/Amaury_dev/libf/dyn3d/lmdz_dynetat0.f90 (revision 5228) +++ LMDZ6/branches/Amaury_dev/libf/dyn3d/lmdz_dynetat0.f90 (revision 5229) @@ -12,5 +12,5 @@ !------------------------------------------------------------------------------- USE lmdz_infotrac, ONLY: nqtot, tracers, niso, iqIsoPha, iH2O, isoName, & - new2oldH2O, newHNO3, oldHNO3, getKey + new2oldH2O, newHNO3, oldHNO3 USE lmdz_strings, ONLY: maxlen, msg, strStack, real2str, int2str USE netcdf, ONLY: nf90_open, nf90_nowrite, nf90_inq_varid, nf90_close, nf90_get_var, nf90_noerr @@ -28,4 +28,5 @@ USE lmdz_strings, ONLY: strIdx USE ioipsl, ONLY: getin + USE iso_params_mod ! tnat_* and alpha_ideal_* USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm @@ -169,6 +170,12 @@ CALL msg(' !!! Beware: alpha_ideal put to 1 !!!', modname) ELSE - IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) & - CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1) + SELECT CASE(isoName(iName)) + CASE('H216O'); tnat = tnat_H216O; alpha_ideal = alpha_ideal_H216O + CASE('H217O'); tnat = tnat_H217O; alpha_ideal = alpha_ideal_H217O + CASE('H218O'); tnat = tnat_H218O; alpha_ideal = alpha_ideal_H218O + CASE('HDO'); tnat = tnat_HDO; alpha_ideal = alpha_ideal_HDO + CASE('HTO'); tnat = tnat_HTO; alpha_ideal = alpha_ideal_HTO + CASE DEFAULT; CALL abort_gcm(TRIM(modname), 'unknown isotope "' // TRIM(isoName(iName)) // '" ; check tracer.def file', 1) + END SELECT END IF CALL msg('Missing tracer <' // TRIM(var) // '> => initialized with a simplified Rayleigh distillation law.', modname) Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/iso_verif_dyn.f90 =================================================================== --- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/iso_verif_dyn.f90 (revision 5228) +++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/iso_verif_dyn.f90 (revision 5229) @@ -1,114 +1,65 @@ -function iso_verif_noNaN_nostop(x, err_msg) - IMPLICIT NONE - ! si x est NaN, on affiche message - ! d'erreur et return 1 si erreur +LOGICAL FUNCTION iso_verif_noNaN_nostop(x, err_msg) RESULT(lerr) + IMPLICIT NONE +! If "x" is NaN, print an error message . + REAL, INTENT(IN) :: x + CHARACTER(LEN=*), INTENT(IN) :: err_msg + REAL, PARAMETER :: borne = 1e19 + lerr = x > -borne .AND. x < borne + IF(.NOT.lerr) RETURN + WRITE(*,*) 'erreur detectee par iso_verif_nonNaN:' + WRITE(*,*) err_msg + WRITE(*,*) 'x=',x +END FUNCTION iso_verif_noNaN_nostop - ! input: - REAL :: x - CHARACTER(LEN = *) :: err_msg ! message d''erreur à afficher - - ! output - REAL :: borne - parameter (borne = 1e19) - INTEGER :: iso_verif_noNaN_nostop - - IF ((x>-borne).AND.(xerrmax) THEN - IF (abs((a - b) / max(max(abs(b), abs(a)), 1e-18)) & - >errmaxrel) THEN - WRITE(*, *) 'erreur detectee par iso_verif_egalite:' - WRITE(*, *) err_msg - WRITE(*, *) 'a=', a - WRITE(*, *) 'b=', b - iso_verif_egalite_nostop = 1 - endif - endif - - RETURN +LOGICAL FUNCTION iso_verif_egalite_nostop(a, b, err_msg) RESULT(lerr) + IMPLICIT NONE +! Compare "a" and "b". If a/=b, print an error message. +! Both absolute and relative errors are checked for equality. + REAL, INTENT(IN) :: a, b + CHARACTER(LEN=*), INTENT(IN) :: err_msg + REAL, PARAMETER :: errmax = 1e-8, & ! max absolute error + errmaxrel = 1e-3 ! max relative error + lerr = ABS(a-b) > errmax + IF(.NOT.lerr) RETURN + lerr = ABS( (a-b) / MAX(MAX(ABS(b), ABS(a)),1e-18) ) > errmaxrel + IF(.NOT.lerr) RETURN + WRITE(*,*) 'erreur detectee par iso_verif_egalite:' + WRITE(*,*) err_msg + WRITE(*,*) 'a=',a + WRITE(*,*) 'b=',b END FUNCTION iso_verif_egalite_nostop -function iso_verif_aberrant_nostop & - (x, iso, q, err_msg) - USE lmdz_infotrac, ONLY: isoName, getKey +LOGICAL FUNCTION iso_verif_aberrant_nostop(x, iso, q, err_msg) RESULT(lerr) USE ioipsl, ONLY: getin + USE iso_params_mod, ONLY: tnat_HDO IMPLICIT NONE - ! input: - REAL :: x, q - INTEGER :: iso ! 2=HDO, 1=O18 - CHARACTER(LEN = *) :: err_msg ! message d''erreur à afficher + REAL, INTENT(IN) :: x, q + INTEGER, INTENT(IN) :: iso ! 2=HDO, 1=O18 + CHARACTER(LEN = *), INTENT(IN) :: err_msg - ! locals - REAL :: qmin, deltaD - REAL :: deltaDmax, deltaDmin, tnat - parameter (qmin = 1e-11) - parameter (deltaDmax = 200.0, deltaDmin = -999.9) - LOGICAL, SAVE :: ltnat1 - LOGICAL, SAVE :: lFirst = .TRUE. - - ! output - INTEGER :: iso_verif_aberrant_nostop - - IF(lFirst) THEN - ltnat1 = .TRUE.; CALL getin('tnateq1', ltnat1) - lFirst = .FALSE. - END IF - iso_verif_aberrant_nostop = 0 - - ! verifier que HDO est raisonable - IF (q>qmin) THEN - IF(ltnat1) THEN - tnat = 1.0 - ELSE IF(getKey('tnat', tnat, isoName(iso))) THEN - err_msg = 'Missing isotopic parameter "tnat"' - iso_verif_aberrant_nostop = 1 - RETURN - END IF - deltaD = (x / q / tnat - 1) * 1000 - IF ((deltaD>deltaDmax).OR.(deltaD qmin + IF(.NOT.lerr) RETURN + deltaD = (x / q /tnat - 1.) * 1000. + lerr = deltaD > deltaDmax .OR. deltaD < deltaDmin + IF(.NOT.lerr) RETURN + WRITE(*,*) 'erreur detectee par iso_verif_aberrant:' + WRITE(*,*) err_msg + WRITE(*,*) 'q=',q + WRITE(*,*) 'deltaD=',deltaD + WRITE(*,*) 'iso=',iso END FUNCTION iso_verif_aberrant_nostop - Index: LMDZ6/branches/Amaury_dev/libf/dyn3dmem/check_isotopes_loc.F90 =================================================================== --- LMDZ6/branches/Amaury_dev/libf/dyn3dmem/check_isotopes_loc.F90 (revision 5228) +++ LMDZ6/branches/Amaury_dev/libf/dyn3dmem/check_isotopes_loc.F90 (revision 5229) @@ -3,5 +3,6 @@ USE lmdz_strings, ONLY: maxlen, msg, strIdx, strStack, int2str, real2str USE lmdz_infotrac, ONLY: nqtot, niso, nphas, isotope, isoCheck, iqIsoPha, isoSelect, & - ntiso, iH2O, nzone, tracers, isoName, itZonIso, getKey + ntiso, iH2O, nzone, tracers, isoName, itZonIso + USE iso_params_mod, ONLY: tnat_H216O, tnat_H217O, tnat_H218O, tnat_HDO, tnat_HTO @@ -14,7 +15,7 @@ CHARACTER(LEN=maxlen) :: modname, msg1, nm(2) INTEGER :: ixt, ipha, k, i, iq, iiso, izon, ieau, iqeau, iqpar - INTEGER, ALLOCATABLE :: ix(:) + INTEGER, ALLOCATABLE :: ix(:) REAL, ALLOCATABLE, SAVE :: tnat(:) !--- OpenMP shared variable - REAL :: xtractot, xiiso, deltaD, q1, q2 + REAL :: xtractot, xiiso, deltaD, q1, q2 REAL, PARAMETER :: borne = 1e19, & errmax = 1e-8, & !--- Max. absolute error @@ -38,14 +39,10 @@ ltnat1 = .TRUE.; CALL getin('tnateq1', ltnat1) ALLOCATE(tnat(niso)) - iso_eau = strIdx(isoName,'H216O') - iso_O17 = strIdx(isoName,'H217O') - iso_O18 = strIdx(isoName,'H218O') - iso_HDO = strIdx(isoName,'HDO') - iso_HTO = strIdx(isoName,'HTO') - IF(ltnat1) THEN - tnat(:)=1.0 - ELSE - IF(getKey('tnat', tnat)) CALL abort_gcm(modname, 'missing isotopic parameter', 1) - END IF + iso_eau = strIdx(isoName,'H216O'); IF(iso_eau /= 0) tnat(iso_eau) = tnat_H216O + iso_O17 = strIdx(isoName,'H217O'); IF(iso_O17 /= 0) tnat(iso_O17) = tnat_H217O + iso_O18 = strIdx(isoName,'H218O'); IF(iso_O18 /= 0) tnat(iso_O18) = tnat_H218O + iso_HDO = strIdx(isoName,'HDO'); IF(iso_HDO /= 0) tnat(iso_HDO) = tnat_HDO + iso_HTO = strIdx(isoName,'HTO'); IF(iso_HTO /= 0) tnat(iso_HTO) = tnat_HTO + IF(ltnat1) tnat(:) = 1.0 !$OMP END MASTER !$OMP BARRIER @@ -62,5 +59,5 @@ DO k = 1, llm DO i = ijb, ije - IF(ABS(q(i,k,iq)) <= borne) CYCLE + IF(ABS(q(i,k,iq)) < borne) CYCLE WRITE(msg1,'(s,"(",i0,",",i0,",",i0,") = ",ES12.4)')TRIM(isoName(ixt)),i,k,iq,q(i,k,iq) CALL msg(msg1, modname) Index: LMDZ6/branches/Amaury_dev/libf/dyn3dmem/dynetat0_loc.f90 =================================================================== --- LMDZ6/branches/Amaury_dev/libf/dyn3dmem/dynetat0_loc.f90 (revision 5228) +++ LMDZ6/branches/Amaury_dev/libf/dyn3dmem/dynetat0_loc.f90 (revision 5229) @@ -8,5 +8,5 @@ USE parallel_lmdz USE lmdz_infotrac, ONLY: nqtot, tracers, niso, iqIsoPha, iH2O, isoName, & - new2oldH2O, newHNO3, oldHNO3, getKey + new2oldH2O, newHNO3, oldHNO3 USE lmdz_strings, ONLY: maxlen, msg, strStack, real2str, int2str, strIdx USE netcdf, ONLY: nf90_open, nf90_nowrite, nf90_inquire_dimension, nf90_inq_varid, & @@ -23,4 +23,5 @@ USE lmdz_iniprint, ONLY: lunout, prt_level USE lmdz_comgeom + USE iso_params_mod ! tnat_* and alpha_ideal_* USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm @@ -188,6 +189,12 @@ CALL msg(' !!! Beware: alpha_ideal put to 1 !!!', modname) ELSE - IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) & - CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1) + SELECT CASE(isoName(iName)) + CASE('H216O'); tnat = tnat_H216O; alpha_ideal = alpha_ideal_H216O + CASE('H217O'); tnat = tnat_H217O; alpha_ideal = alpha_ideal_H217O + CASE('H218O'); tnat = tnat_H218O; alpha_ideal = alpha_ideal_H218O + CASE('HDO'); tnat = tnat_HDO; alpha_ideal = alpha_ideal_HDO + CASE('HTO'); tnat = tnat_HTO; alpha_ideal = alpha_ideal_HTO + CASE DEFAULT; CALL abort_gcm(TRIM(modname), 'unknown isotope "' // TRIM(isoName(iName)) // '" ; check tracer.def file', 1) + END SELECT END IF CALL msg('Missing tracer <' // TRIM(var) // '> => initialized with a simplified Rayleigh distillation law.', modname) Index: LMDZ6/branches/Amaury_dev/libf/dyn3dmem/iniacademic_loc.F90 =================================================================== --- LMDZ6/branches/Amaury_dev/libf/dyn3dmem/iniacademic_loc.F90 (revision 5228) +++ LMDZ6/branches/Amaury_dev/libf/dyn3dmem/iniacademic_loc.F90 (revision 5229) @@ -1,8 +1,6 @@ -! $Id: iniacademic.F90 1625 2012-05-09 13:14:48Z lguez $ - SUBROUTINE iniacademic_loc(vcov, ucov, teta, q, masse, ps, phis, time_0) USE lmdz_filtreg, ONLY: inifilr - USE lmdz_infotrac, ONLY: nqtot, niso, iqIsoPha, tracers, getKey, isoName, addPhase + USE lmdz_infotrac, ONLY: nqtot, niso, iqIsoPha, tracers, isoName, addPhase USE control_mod, ONLY: day_step, planet_type USE exner_hyb_m, ONLY: exner_hyb @@ -21,4 +19,5 @@ USE lmdz_academic, ONLY: tetarappel, knewt_t, kfrict, knewt_g, clat4 USE lmdz_comgeom + USE iso_params_mod ! tnat_* and alpha_ideal_* ! Author: Frederic Hourdin original: 15/01/93 @@ -322,6 +321,13 @@ WRITE(lunout, *) 'In '//TRIM(modname)//': !!! Beware: alpha_ideal put to 1 !!!' ELSE - IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) & - CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1) + SELECT CASE(isoName(iName)) + CASE('H216O'); tnat = tnat_H216O; alpha_ideal = alpha_ideal_H216O + CASE('H217O'); tnat = tnat_H217O; alpha_ideal = alpha_ideal_H217O + CASE('H218O'); tnat = tnat_H218O; alpha_ideal = alpha_ideal_H218O + CASE('HDO'); tnat = tnat_HDO; alpha_ideal = alpha_ideal_HDO + CASE('HTO'); tnat = tnat_HTO; alpha_ideal = alpha_ideal_HTO + CASE DEFAULT + CALL abort_gcm(TRIM(modname),'unknown isotope "'//TRIM(isoName(iName))//'" ; check tracer.def file',1) + END SELECT END IF q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqParent)*tnat*(q(ijb_u:ije_u,:,iqParent)/30.e-3)**(alpha_ideal-1.) Index: LMDZ6/branches/Amaury_dev/libf/misc/iso_params_mod.f90 =================================================================== --- LMDZ6/branches/Amaury_dev/libf/misc/iso_params_mod.f90 (revision 5229) +++ LMDZ6/branches/Amaury_dev/libf/misc/iso_params_mod.f90 (revision 5229) @@ -0,0 +1,10 @@ +MODULE iso_params_mod + + REAL, PARAMETER :: tnat_H216O = 1.0000, alpha_ideal_H216O = 1.000, & + tnat_H217O = 40.000e-6, alpha_ideal_H217O = 1.003, & + tnat_H218O = 2005.2e-6, alpha_ideal_H218O = 1.006, & + tnat_HDO = 155.76e-6, alpha_ideal_HDO = 1.010, & + tnat_HTO = 0.0000, alpha_ideal_HTO = 1.000 + +END MODULE iso_params_mod + Index: LMDZ6/branches/Amaury_dev/libf/phylmdiso/isotopes_mod.F90 =================================================================== --- LMDZ6/branches/Amaury_dev/libf/phylmdiso/isotopes_mod.F90 (revision 5228) +++ LMDZ6/branches/Amaury_dev/libf/phylmdiso/isotopes_mod.F90 (revision 5229) @@ -4,5 +4,6 @@ MODULE isotopes_mod USE lmdz_strings, ONLY: msg, real2str, int2str, bool2str, maxlen, strIdx, strStack - USE infotrac_phy, ONLY: isoName + USE infotrac_phy, ONLY: isoName, niso, ntiso + USE iso_params_mod IMPLICIT NONE INTERFACE get_in; MODULE PROCEDURE getinp_s, getinp_i, getinp_r, getinp_l; END INTERFACE get_in @@ -10,4 +11,35 @@ !--- Contains all isotopic variables + their initialization !--- Isotopes-specific routines are in isotopes_routines_mod to avoid circular dependencies with isotopes_verif_mod. + + !--- Negligible lower thresholds: no need to check for absurd values under these lower limits + REAL, PARAMETER :: & + ridicule = 1e-12, & ! For mixing ratios + ridicule_rain = 1e-8, & ! For rain fluxes (rain, zrfl...) in kg/s <-> 1e-3 mm/day + ridicule_evap = ridicule_rain*1e-2, & ! For evaporations in kg/s <-> 1e-3 mm/day + ridicule_qsol = ridicule_rain, & ! For qsol in kg <-> 1e-8 kg + ridicule_snow = ridicule_qsol ! For snow in kg <-> 1e-8 kg + REAL, PARAMETER :: expb_max = 30.0 + + !--- Fractionation coefficients for H217O + REAL, PARAMETER :: fac_coeff_eq17_liq = 0.529, & + fac_coeff_eq17_ice = 0.529 + + !--- H218O reference + REAL, PARAMETER :: fac_enrichoce18 = 0.0005, alpha_liq_sol_O18 = 1.00291, & + talph1_O18 = 1137., talps1_O18 = 11.839, tkcin0_O18 = 0.006, & + talph2_O18 = -0.4156, talps2_O18 = -0.028244, tkcin1_O18 = 0.000285, & + talph3_O18 = -2.0667E-3, tdifrel_O18 = 1./0.9723, tkcin2_O18 = 0.00082 + + !--- Parameters that do not depend on the nature of water isotopes: + REAL, PARAMETER :: pxtmelt= 273.15, & !--- temperature at which ice formation starts + pxtice = 273.15 - 10.0, & !--- temperature at which all condensate is ice: + pxtmin = 273.15 - 120.0, & !--- computation done only under -120°C + pxtmax = 273.15 + 60.0, & !--- computation done only over +60°C + tdifexp= 0.58, & !--- a constant for alpha_eff for equilibrium below cloud base: + tv0cin = 7.0, & !--- wind threshold (m/s) for smooth/rough regime (Merlivat and Jouzel 1979) + musi = 1.0, & !--- facteurs lambda et mu dans Si=musi-lambda*T + Kd = 2.5e-9, & !--- diffusion in soil ; m2/s + rh_cste_surf_cond = 0.6, & !--- cste_surf_cond case: rhs and/or Ts set to constants + T_cste_surf_cond = 288.0 !--- Isotopes indices (in [1,niso] ; non-existing => 0 index) @@ -88,6 +120,6 @@ !--- Vectors of length "niso" REAL, ALLOCATABLE, DIMENSION(:), SAVE :: & - tnat, toce, tcorr, tdifrel -!$OMP THREADPRIVATE(tnat, toce, tcorr, tdifrel) + alpha, tnat, toce, tcorr, tdifrel +!$OMP THREADPRIVATE(alpha, tnat, toce, tcorr, tdifrel) REAL, ALLOCATABLE, DIMENSION(:), SAVE :: & talph1, talph2, talph3, talps1, talps2 @@ -99,43 +131,4 @@ alpha_liq_sol, Rdefault, Rmethox !$OMP THREADPRIVATE(alpha_liq_sol, Rdefault, Rmethox) -! REAL, SAVE :: fac_coeff_eq17_liq, fac_coeff_eq17_ice -!!$OMP THREADPRIVATE(fac_coeff_eq17_liq, fac_coeff_eq17_ice) - - !--- H2[18]O reference - REAL, PARAMETER :: fac_enrichoce18=0.0005 - REAL, PARAMETER :: alpha_liq_sol_O18=1.00291 - REAL, PARAMETER :: talph1_O18=1137. - REAL, PARAMETER :: talph2_O18=-0.4156 - REAL, PARAMETER :: talph3_O18=-2.0667E-3 - REAL, PARAMETER :: talps1_O18=11.839 - REAL, PARAMETER :: talps2_O18=-0.028244 - REAL, PARAMETER :: tdifrel_O18=1./0.9723 - REAL, PARAMETER :: tkcin0_O18=0.006 - REAL, PARAMETER :: tkcin1_O18=0.000285 - REAL, PARAMETER :: tkcin2_O18=0.00082 - REAL, PARAMETER :: fac_coeff_eq17_liq=0.529 - REAL, PARAMETER :: fac_coeff_eq17_ice=0.529 - - !---- Parameters that do not depend on the nature of water isotopes: - REAL, PARAMETER :: pxtmelt = 273.15 ! temperature at which ice formation starts - REAL, PARAMETER :: pxtice = 273.15-10.0 ! -- temperature at which all condensate is ice: - REAL, PARAMETER :: pxtmin = 273.15 - 120.0 ! On ne calcule qu'au dessus de -120°C - REAL, PARAMETER :: pxtmax = 273.15 + 60.0 ! On ne calcule qu'au dessus de +60°C - REAL, PARAMETER :: tdifexp = 0.58 ! -- a constant for alpha_eff for equilibrium below cloud base: - REAL, PARAMETER :: tv0cin = 7.0 ! wind threshold (m/s) for smooth/rough regime (Merlivat and Jouzel 1979) - REAL, PARAMETER :: musi=1.0 ! facteurs lambda et mu dans Si=musi-lambda*T - REAL, PARAMETER :: Kd=2.5e-9 ! m2/s ! diffusion dans le sol - REAL, PARAMETER :: rh_cste_surf_cond = 0.6 ! cas où cste_surf_cond: on met rhs ou/et Ts cste pour voir - REAL, PARAMETER :: T_cste_surf_cond = 288.0 - - - !--- Negligible lower thresholds: no need to check for absurd values under these lower limits - REAL, PARAMETER :: & - ridicule = 1e-12, & ! For mixing ratios - ridicule_rain = 1e-8, & ! For rain fluxes (rain, zrfl...) in kg/s <-> 1e-3 mm/day - ridicule_evap = ridicule_rain*1e-2, & ! For evaporations in kg/s <-> 1e-3 mm/day - ridicule_qsol = ridicule_rain, & ! For qsol in kg <-> 1e-8 kg - ridicule_snow = ridicule_qsol ! For snow in kg <-> 1e-8 kg - REAL, PARAMETER :: expb_max = 30.0 !--- Specific to HTO: @@ -154,10 +147,8 @@ SUBROUTINE iso_init() - USE infotrac_phy, ONLY: ntiso, niso, getKey - USE lmdz_strings, ONLY: maxlen IMPLICIT NONE !=== Local variables: - INTEGER :: ixt + INTEGER :: ixt, is LOGICAL :: ltnat1 CHARACTER(LEN=maxlen) :: modname, sxt @@ -184,230 +175,204 @@ CALL msg('64: niso = '//TRIM(int2str(niso)), modname) - !--- Init de ntracisoOR: on ecrasera en cas de traceurs de tagging isotopiques - ! (nzone>0) si complications avec ORCHIDEE - ntracisoOR = ntiso - - !--- Type of water isotopes: - iso_eau = strIdx(isoName, 'H216O'); CALL msg('iso_eau='//int2str(iso_eau), modname) - iso_HDO = strIdx(isoName, 'HDO'); CALL msg('iso_HDO='//int2str(iso_HDO), modname) - iso_O18 = strIdx(isoName, 'H218O'); CALL msg('iso_O18='//int2str(iso_O18), modname) - iso_O17 = strIdx(isoName, 'H217O'); CALL msg('iso_O17='//int2str(iso_O17), modname) - iso_HTO = strIdx(isoName, 'HTO'); CALL msg('iso_HTO='//int2str(iso_HTO), modname) - - !--- Initialiaation: reading the isotopic parameters. - CALL get_in('lambda', lambda_sursat, 0.004) - CALL get_in('thumxt1', thumxt1, 0.75*1.2) - CALL get_in('ntot', ntot, 20, .FALSE.) - CALL get_in('h_land_ice', h_land_ice, 20., .FALSE.) - CALL get_in('P_veg', P_veg, 1.0, .FALSE.) - CALL get_in('bidouille_anti_divergence', bidouille_anti_divergence, .FALSE.) - CALL get_in('essai_convergence', essai_convergence, .FALSE.) - CALL get_in('initialisation_iso', initialisation_iso, 0) - -! IF(nzone>0 .AND. initialisation_iso==0) & -! CALL get_in('initialisation_isotrac',initialisation_isotrac) - CALL get_in('modif_sst', modif_sst, 0) - CALL get_in('deltaTtest', deltaTtest, 0.0) !--- For modif_sst>=1 - CALL get_in('deltaTtestpoles',deltaTtestpoles, 0.0) !--- For modif_sst>=2 - CALL get_in( 'sstlatcrit', sstlatcrit, 30.0) !--- For modif_sst>=3 - CALL get_in('dsstlatcrit', dsstlatcrit, 0.0) !--- For modif_sst>=3 + !--- Init de ntracisoOR: on ecrasera en cas de traceurs de tagging isotopiques + ! (nzone>0) si complications avec ORCHIDEE + ntracisoOR = ntiso + + !--- Type of water isotopes: + iso_eau = strIdx(isoName, 'H216O'); CALL msg('iso_eau='//int2str(iso_eau), modname) + iso_HDO = strIdx(isoName, 'HDO'); CALL msg('iso_HDO='//int2str(iso_HDO), modname) + iso_O18 = strIdx(isoName, 'H218O'); CALL msg('iso_O18='//int2str(iso_O18), modname) + iso_O17 = strIdx(isoName, 'H217O'); CALL msg('iso_O17='//int2str(iso_O17), modname) + iso_HTO = strIdx(isoName, 'HTO'); CALL msg('iso_HTO='//int2str(iso_HTO), modname) + + !--- Initialisation: reading the isotopic parameters. + CALL get_in('lambda', lambda_sursat, 0.004); IF(ok_nocinsat) lambda_sursat = 0. + CALL get_in('thumxt1', thumxt1, 0.75*1.2) + CALL get_in('ntot', ntot, 20, .FALSE.) + CALL get_in('h_land_ice', h_land_ice, 20., .FALSE.) + CALL get_in('P_veg', P_veg, 1.0, .FALSE.) + CALL get_in('bidouille_anti_divergence', bidouille_anti_divergence, .FALSE.) + CALL get_in('essai_convergence', essai_convergence, .FALSE.) + CALL get_in('initialisation_iso', initialisation_iso, 0) + +! IF(nzone>0 .AND. initialisation_iso==0) & +! CALL get_in('initialisation_isotrac',initialisation_isotrac) + CALL get_in('modif_sst', modif_sst, 0) + CALL get_in('deltaTtest', deltaTtest, 0.0) !--- For modif_sst>=1 + CALL get_in('deltaTtestpoles',deltaTtestpoles, 0.0) !--- For modif_sst>=2 + CALL get_in( 'sstlatcrit', sstlatcrit, 30.0) !--- For modif_sst>=3 + CALL get_in('dsstlatcrit', dsstlatcrit, 0.0) !--- For modif_sst>=3 #ifdef ISOVERIF - CALL msg('iso_init 270: sstlatcrit='//real2str( sstlatcrit), modname, sstlatcrit < 0.0) !--- For modif_sst>=2 - CALL msg('iso_init 279: dsstlatcrit='//real2str(dsstlatcrit), modname, sstlatcrit < 0.0) !--- For modif_sst>=3 - IF(modif_sst >= 2 .AND. sstlatcrit < 0.0) STOP + CALL msg('iso_init 270: sstlatcrit='//real2str( sstlatcrit), modname, sstlatcrit < 0.0) !--- For modif_sst>=2 + CALL msg('iso_init 279: dsstlatcrit='//real2str(dsstlatcrit), modname, sstlatcrit < 0.0) !--- For modif_sst>=3 + IF(modif_sst >= 2 .AND. sstlatcrit < 0.0) STOP #endif - CALL get_in('modif_sic', modif_sic, 0) - IF(modif_sic >= 1) & - CALL get_in('deltasic', deltasic, 0.1) - - CALL get_in('albedo_prescrit', albedo_prescrit, 0) - IF(albedo_prescrit == 1) THEN - CALL get_in('lon_min_albedo', lon_min_albedo, -200.) - CALL get_in('lon_max_albedo', lon_max_albedo, 200.) - CALL get_in('lat_min_albedo', lat_min_albedo, -100.) - CALL get_in('lat_max_albedo', lat_max_albedo, 100.) - END IF - CALL get_in('deltaP_BL', deltaP_BL, 10.0) - CALL get_in('ruissellement_pluie', ruissellement_pluie, 0) - CALL get_in('alphak_stewart', alphak_stewart, 1) - CALL get_in('tdifexp_sol', tdifexp_sol, 0.67) - CALL get_in('calendrier_guide', calendrier_guide, 0) - CALL get_in('cste_surf_cond', cste_surf_cond, 0) - CALL get_in('mixlen', mixlen, 35.0) - CALL get_in('evap_cont_cste', evap_cont_cste, 0) - CALL get_in('deltaO18_evap_cont', deltaO18_evap_cont,0.0) - CALL get_in('d_evap_cont', d_evap_cont, 0.0) - CALL get_in('nudge_qsol', nudge_qsol, 0) - CALL get_in('region_nudge_qsol', region_nudge_qsol, 1) - nlevmaxO17 = 50 - CALL msg('nlevmaxO17='//TRIM(int2str(nlevmaxO17))) - CALL get_in('no_pce', no_pce, 0) - CALL get_in('A_satlim', A_satlim, 1.0) - CALL get_in('ok_restrict_A_satlim', ok_restrict_A_satlim, 0) + CALL get_in('modif_sic', modif_sic, 0) + IF(modif_sic >= 1) & + CALL get_in('deltasic', deltasic, 0.1) + + CALL get_in('albedo_prescrit', albedo_prescrit, 0) + IF(albedo_prescrit == 1) THEN + CALL get_in('lon_min_albedo', lon_min_albedo, -200.) + CALL get_in('lon_max_albedo', lon_max_albedo, 200.) + CALL get_in('lat_min_albedo', lat_min_albedo, -100.) + CALL get_in('lat_max_albedo', lat_max_albedo, 100.) + END IF + CALL get_in('deltaO18_oce', deltaO18_oce, 0.0) + CALL get_in('deltaP_BL', deltaP_BL, 10.0) + CALL get_in('ruissellement_pluie', ruissellement_pluie, 0) + CALL get_in('alphak_stewart', alphak_stewart, 1) + CALL get_in('tdifexp_sol', tdifexp_sol, 0.67) + CALL get_in('calendrier_guide', calendrier_guide, 0) + CALL get_in('cste_surf_cond', cste_surf_cond, 0) + CALL get_in('mixlen', mixlen, 35.0) + CALL get_in('evap_cont_cste', evap_cont_cste, 0) + CALL get_in('deltaO18_evap_cont', deltaO18_evap_cont,0.0) + CALL get_in('d_evap_cont', d_evap_cont, 0.0) + CALL get_in('nudge_qsol', nudge_qsol, 0) + CALL get_in('region_nudge_qsol', region_nudge_qsol, 1) + nlevmaxO17 = 50 + CALL msg('nlevmaxO17='//TRIM(int2str(nlevmaxO17))) + CALL get_in('no_pce', no_pce, 0) + CALL get_in('A_satlim', A_satlim, 1.0) + CALL get_in('ok_restrict_A_satlim', ok_restrict_A_satlim, 0) #ifdef ISOVERIF - CALL msg(' 315: A_satlim='//real2str(A_satlim), modname, A_satlim > 1.0) - IF(A_satlim > 1.0) STOP + CALL msg(' 315: A_satlim='//real2str(A_satlim), modname, A_satlim > 1.0) + IF(A_satlim > 1.0) STOP #endif -! CALL get_in('slope_limiterxy', slope_limiterxy, 2.0) -! CALL get_in('slope_limiterz', slope_limiterz, 2.0) - CALL get_in('modif_ratqs', modif_ratqs, 0) - CALL get_in('Pcrit_ratqs', Pcrit_ratqs, 500.0) - CALL get_in('ratqsbasnew', ratqsbasnew, 0.05) - CALL get_in('fac_modif_evaoce', fac_modif_evaoce, 1.0) - CALL get_in('ok_bidouille_wake', ok_bidouille_wake, 0) - ! si oui, la temperature de cond est celle de l'environnement, pour eviter - ! bugs quand temperature dans ascendances convs est mal calculee - CALL get_in('cond_temp_env', cond_temp_env, .FALSE.) - IF(ANY(isoName == 'HTO')) & - CALL get_in('ok_prod_nucl_tritium', ok_prod_nucl_tritium, .FALSE., .FALSE.) - CALL get_in('tnateq1', ltnat1, .TRUE.) - - ! Ocean composition - CALL get_in('deltaO18_oce', deltaO18_oce, 0.0) - - CALL msg('iso_O18, iso_HDO, iso_eau = '//TRIM(strStack(int2str([iso_O18, iso_HDO, iso_eau]))), modname) - - !-------------------------------------------------------------- - ! Isotope fractionation factors and a few isotopic constants - !-------------------------------------------------------------- - ALLOCATE(tkcin0(niso)) - ALLOCATE(tkcin1(niso)) - ALLOCATE(tkcin2(niso)) - ALLOCATE(tnat(niso)) - ALLOCATE(tdifrel(niso)) - ALLOCATE(toce(niso)) - ALLOCATE(tcorr(niso)) - ALLOCATE(talph1(niso)) - ALLOCATE(talph2(niso)) - ALLOCATE(talph3(niso)) - ALLOCATE(talps1(niso)) - ALLOCATE(talps2(niso)) - ALLOCATE(alpha_liq_sol(niso)) - ALLOCATE(Rdefault(niso)) - ALLOCATE(Rmethox(niso)) - - DO ixt=1,niso - IF (ixt.EQ.iso_HTO) then ! Tritium - tkcin0(ixt) = 0.01056 - tkcin1(ixt) = 0.0005016 - tkcin2(ixt) = 0.0014432 - tnat(ixt) = 0.0; IF(ltnat1) tnat(ixt)=1 - 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. - Rmethox(ixt)=0.0 - endif - IF (ixt.EQ.iso_O17) then ! O17 - pente_MWL=0.528 - tdifrel(ixt)=1./0.98555 ! valeur utilisée en 1D et dans modèle de LdG ! tdifrel(ixt)=1./0.985452 ! donné par Amaelle - fac_kcin= (tdifrel(ixt)-1.0)/(tdifrel_O18-1.0) ! fac_kcin=0.5145 ! donné par Amaelle - 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 - IF(ltnat1) tnat(ixt)=1 - 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 - Rdefault(ixt)=tnat(ixt)*(-3.15/1000.0+1.0) - Rmethox(ixt)=(230./1000.+1.)*tnat(ixt) !Zahn et al 2006 - 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; IF(ltnat1) tnat(ixt)=1 - 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 - Rdefault(ixt)=tnat(ixt)*(-6.0/1000.0+1.0) - Rmethox(ixt)=(130./1000.+1.)*tnat(ixt) !Zahn et al 2006 - endif - IF (ixt.EQ.iso_HDO) then ! Deuterium - pente_MWL=8.0 - tdifrel(ixt)=1./0.9755 ! fac_kcin=0.88 - 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; IF(ltnat1) tnat(ixt)=1 - 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 - Rdefault(ixt)=tnat(ixt)*((-6.0*pente_MWL+10.0)/1000.0+1.0) - Rmethox(ixt)=tnat(ixt)*(-25.0/1000.+1.) ! Zahn et al 2006 - endif - 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. - Rdefault(ixt)=tnat(ixt)*1.0 - Rmethox(ixt)=1.0 - endif - enddo ! ixt=1,niso - - IF(.NOT.Rdefault_smow) THEN - Rdefault(:) = 0.0 - IF (iso_eau.gt.0) Rdefault(iso_eau) = 1.0 ! correction Camille 30 mars 2023 - ENDIF - WRITE(*,*) 'Rdefault=',Rdefault - WRITE(*,*) 'toce=',toce - - !--- Sensitivity test: no kinetic effect in sfc evaporation - IF(ok_nocinsfc) THEN - tkcin0(1:niso) = 0.0 - tkcin1(1:niso) = 0.0 - tkcin2(1:niso) = 0.0 - END IF - - CALL msg('285: verif initialisation:', modname) - DO ixt=1,niso - sxt=int2str(ixt) - CALL msg(' * isoName('//TRIM(sxt)//') = <'//TRIM(isoName(ixt))//'>', modname) - CALL msg( ' tnat('//TRIM(sxt)//') = '//TRIM(real2str(tnat(ixt))), modname) -! CALL msg(' alpha_liq_sol('//TRIM(sxt)//') = '//TRIM(real2str(alpha_liq_sol(ixt))), modname) -! CALL msg( ' tkcin0('//TRIM(sxt)//') = '//TRIM(real2str(tkcin0(ixt))), modname) -! CALL msg( ' tdifrel('//TRIM(sxt)//') = '//TRIM(real2str(tdifrel(ixt))), modname) - END DO - CALL msg('69: lambda = '//TRIM(real2str(lambda_sursat)), modname) - CALL msg('69: thumxt1 = '//TRIM(real2str(thumxt1)), modname) - CALL msg('69: h_land_ice = '//TRIM(real2str(h_land_ice)), modname) - CALL msg('69: P_veg = '//TRIM(real2str(P_veg)), modname) +! CALL get_in('slope_limiterxy', slope_limiterxy, 2.0) +! CALL get_in('slope_limiterz', slope_limiterz, 2.0) + CALL get_in('modif_ratqs', modif_ratqs, 0) + CALL get_in('Pcrit_ratqs', Pcrit_ratqs, 500.0) + CALL get_in('ratqsbasnew', ratqsbasnew, 0.05) + CALL get_in('fac_modif_evaoce', fac_modif_evaoce, 1.0) + CALL get_in('ok_bidouille_wake', ok_bidouille_wake, 0) + ! si oui, la temperature de cond est celle de l'environnement, pour eviter + ! bugs quand temperature dans ascendances convs est mal calculee + CALL get_in('cond_temp_env', cond_temp_env, .FALSE.) + IF(ANY(isoName == 'HTO')) & + CALL get_in('ok_prod_nucl_tritium', ok_prod_nucl_tritium, .FALSE., .FALSE.) + CALL get_in('tnateq1', ltnat1, .TRUE.) + + CALL msg('iso_O18, iso_HDO, iso_eau = '//TRIM(strStack(int2str([iso_O18, iso_HDO, iso_eau]))), modname) + + !-------------------------------------------------------------- + ! Parameters that depend on the nature of water isotopes: + !-------------------------------------------------------------- + ALLOCATE(tnat (niso), talph1(niso), talps1(niso), tkcin0(niso), tdifrel (niso), alpha (niso)) + ALLOCATE(toce (niso), talph2(niso), talps2(niso), tkcin1(niso), Rdefault(niso), alpha_liq_sol(niso)) + ALLOCATE(tcorr(niso), talph3(niso), tkcin2(niso), Rmethox (niso)) + + !=== H216O + is = iso_eau + IF(is /= 0) THEN + tdifrel (is) = 1.0 + alpha (is) = alpha_ideal_H216O + tnat (is) = tnat_H216O; IF(ltnat1) tnat(is) = 1.0 + toce (is) = tnat(is) + tcorr (is) = 1.0 + talph1 (is) = 0.0; talps1(is) = 0.0; tkcin0(is) = 0.0 + talph2 (is) = 0.0; talps2(is) = 0.0; tkcin1(is) = 0.0 + talph3 (is) = 0.0; tkcin2(is) = 0.0 + Rdefault(is) = tnat(is)*1.0 + Rmethox (is) = 1.0 + alpha_liq_sol(is) = 1.0 + END IF + + !=== H217O + is = iso_O17 + IF(is /= 0) THEN; pente_MWL = 0.528 + tdifrel (is) = 1./0.98555 ! used in 1D and in LdG model ; tdifrel=1./0.985452: from Amaelle + alpha (is) = alpha_ideal_H217O + tnat (is) = tnat_H217O; IF(ltnat1) tnat(is) = 1.0 + toce (is) = tnat(is)*(1.0+deltaO18_oce/1000.0)**pente_MWL + tcorr (is) = 1.0+fac_enrichoce18*pente_MWL + fac_kcin = (tdifrel(is)-1.0)/(tdifrel_O18-1.0) ! fac_kcin=0.5145: from Amaelle + talph1 (is) = talph1_O18; talps1(is) = talps1_O18; tkcin0(is) = tkcin0_O18*fac_kcin + talph2 (is) = talph2_O18; talps2(is) = talps2_O18; tkcin1(is) = tkcin1_O18*fac_kcin + talph3 (is) = talph3_O18; tkcin2(is) = tkcin2_O18*fac_kcin + Rdefault(is) = tnat(is)*(1.0-3.15/1000.) + Rmethox (is) = tnat(is)*(1.0+230./1000.) + alpha_liq_sol(is) = alpha_liq_sol_O18**fac_coeff_eq17_liq + END IF + + !=== H218O + is = iso_O18 + IF(is /= 0) THEN + tdifrel (is) = tdifrel_O18 + alpha (is) = alpha_ideal_H218O + tnat (is) = tnat_H218O; IF(ltnat1) tnat(is) = 1.0 + toce (is) = tnat(is)*(1.0+deltaO18_oce/1000.0) + tcorr (is) = 1.0+fac_enrichoce18 + talph1 (is) = talph1_O18; talps1(is) = talps1_O18; tkcin0(is) = tkcin0_O18 + talph2 (is) = talph2_O18; talps2(is) = talps2_O18; tkcin1(is) = tkcin1_O18 + talph3 (is) = talph3_O18; tkcin2(is) = tkcin2_O18 + Rdefault(is) = tnat(is)*(1.0-6.00/1000.) + Rmethox (is) = tnat(is)*(1.0+130./1000.) ! Zahn & al. 2006 + alpha_liq_sol(is) = alpha_liq_sol_O18 + END IF + + !=== HDO + is = iso_HDO + IF(is /= 0) THEN; pente_MWL = 8.0 + tdifrel (is) = 1./0.9755 ! fac_kcin=0.88 + alpha (is) = alpha_ideal_HDO + tnat (is) = tnat_HDO; IF(ltnat1) tnat(is) = 1.0 + toce (is) = tnat(is)*(1.0+deltaO18_oce/1000.0*pente_MWL) + tcorr (is) = 1.0+fac_enrichoce18*pente_MWL + fac_kcin = (tdifrel(is)-1.0)/(tdifrel_O18-1.0) + talph1 (is) = 24844.; talps1(is) = 16288.; tkcin0(is) = tkcin0_O18*fac_kcin + talph2 (is) = -76.248; talps2(is) = -0.0934; tkcin1(is) = tkcin1_O18*fac_kcin + talph3 (is) = 52.612E-3; tkcin2(is) = tkcin2_O18*fac_kcin + Rdefault(is) = tnat(is)*(1.0+(10.0-6.0*pente_MWL)/1000.) + Rmethox (is) = tnat(is)*(1.0-25.0/1000.) + alpha_liq_sol(is) = 1.0212 ! Lehmann & Siegenthaler, 1991, Jo. of Glaciology, vol 37, p 23 + ! alpha_liq_sol=1.0192: Weston, Ralph, 1955 + END IF + + !=== HTO + is = iso_HTO + IF(is /= 0) THEN + tdifrel (is) = 1./0.968 + alpha (is) = alpha_ideal_HTO + tnat (is) = tnat_HTO; IF(ltnat1) tnat(is) = 1.0 + toce (is) = 4.0E-19 ! ratio T/H = 0.2 TU Dreisigacker & Roether 1978 + tcorr (is) = 1.0 + talph1 (is) = 46480.; talps1(is) = 46480.; tkcin0(is) = 0.01056 + talph2 (is) = -103.87; talps2(is) = -103.87; tkcin1(is) = 0.0005016 + talph3 (is) = 0.0; tkcin2(is) = 0.0014432 + Rdefault(is) = 0.0 + Rmethox (is) = 0.0 + alpha_liq_sol(is) = 1.0 + END IF + + IF(.NOT. Rdefault_smow) THEN + Rdefault(:) = 0.0; IF(iso_eau > 0) Rdefault(iso_eau) = 1.0 + END IF + WRITE(*,*) 'Rdefault = ',Rdefault + WRITE(*,*) 'toce = ', toce + + + !--- Sensitivity test: no kinetic effect in sfc evaporation + IF(ok_nocinsfc) THEN + tkcin0(1:niso) = 0.0 + tkcin1(1:niso) = 0.0 + tkcin2(1:niso) = 0.0 + END IF + + CALL msg('285: verif initialisation:', modname) + DO ixt=1,niso + sxt=int2str(ixt) + CALL msg(' * isoName('//TRIM(sxt)//') = <'//TRIM(isoName(ixt))//'>', modname) + CALL msg( ' tnat('//TRIM(sxt)//') = '//TRIM(real2str(tnat(ixt))), modname) +! CALL msg(' alpha_liq_sol('//TRIM(sxt)//') = '//TRIM(real2str(alpha_liq_sol(ixt))), modname) +! CALL msg( ' tkcin0('//TRIM(sxt)//') = '//TRIM(real2str(tkcin0(ixt))), modname) +! CALL msg( ' tdifrel('//TRIM(sxt)//') = '//TRIM(real2str(tdifrel(ixt))), modname) + END DO + CALL msg('69: lambda = '//TRIM(real2str(lambda_sursat)), modname) + CALL msg('69: thumxt1 = '//TRIM(real2str(thumxt1)), modname) + CALL msg('69: h_land_ice = '//TRIM(real2str(h_land_ice)), modname) + CALL msg('69: P_veg = '//TRIM(real2str(P_veg)), modname) END SUBROUTINE iso_init