Changeset 5214 for LMDZ6/trunk
- Timestamp:
- Sep 22, 2024, 10:07:56 PM (3 months ago)
- Location:
- LMDZ6/trunk/libf
- Files:
-
- 1 added
- 7 edited
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/dyn3d/check_isotopes.F90
r5201 r5214 2 2 USE strings_mod, ONLY: maxlen, msg, strIdx, strStack, int2str, real2str 3 3 USE infotrac, ONLY: nqtot, niso, nphas, isotope, isoCheck, iqIsoPha, isoSelect, & 4 ntiso, iH2O, nzone, tracers, isoName, itZonIso, getKey 4 ntiso, iH2O, nzone, tracers, isoName, itZonIso 5 USE iso_params_mod, ONLY: tnat_H216O, tnat_H217O, tnat_H218O, tnat_HDO, tnat_HTO 5 6 #ifdef CPP_IOIPSL 6 7 USE ioipsl, ONLY: getin … … 15 16 CHARACTER(LEN=maxlen) :: modname, msg1, nm(2) 16 17 INTEGER :: ixt, ipha, k, i, iq, iiso, izon, ieau, iqeau, iqpar 17 INTEGER, ALLOCATABLE :: ix(:)18 INTEGER, ALLOCATABLE :: ix(:) 18 19 REAL, ALLOCATABLE, SAVE :: tnat(:) 19 REAL 20 REAL :: xtractot, xiiso, deltaD, q1, q2 20 21 REAL, PARAMETER :: borne = 1e19, & 21 22 errmax = 1e-8, & !--- Max. absolute error … … 26 27 ridicule = 1e-12 27 28 INTEGER, SAVE :: iso_eau, iso_O17, iso_O18, iso_HDO, iso_HTO 28 LOGICAL, SAVE :: ltnat1, first=.TRUE. 29 LOGICAL :: ltnat1 30 LOGICAL, SAVE :: first=.TRUE. 29 31 30 32 modname='check_isotopes' … … 35 37 ltnat1 = .TRUE.; CALL getin('tnateq1', ltnat1) 36 38 ALLOCATE(tnat(niso)) 37 iso_eau = strIdx(isoName,'H216O') 38 iso_O17 = strIdx(isoName,'H217O') 39 iso_O18 = strIdx(isoName,'H218O') 40 iso_HDO = strIdx(isoName,'HDO') 41 iso_HTO = strIdx(isoName,'HTO') 42 IF(ltnat1) THEN 43 tnat(:)=1.0 44 ELSE 45 IF(getKey('tnat', tnat)) CALL abort_gcm(modname, 'missing isotopic parameter', 1) 46 END IF 39 iso_eau = strIdx(isoName,'H216O'); IF(iso_eau /= 0) tnat(iso_eau) = tnat_H216O 40 iso_O17 = strIdx(isoName,'H217O'); IF(iso_O17 /= 0) tnat(iso_O17) = tnat_H217O 41 iso_O18 = strIdx(isoName,'H218O'); IF(iso_O18 /= 0) tnat(iso_O18) = tnat_H218O 42 iso_HDO = strIdx(isoName,'HDO'); IF(iso_HDO /= 0) tnat(iso_HDO) = tnat_HDO 43 iso_HTO = strIdx(isoName,'HTO'); IF(iso_HTO /= 0) tnat(iso_HTO) = tnat_HTO 44 IF(ltnat1) tnat(:) = 1. 47 45 first = .FALSE. 48 46 END IF … … 56 54 DO k = 1, llm 57 55 DO i = 1, ip1jmp1 58 IF(ABS(q(i,k,iq)) < =borne) CYCLE56 IF(ABS(q(i,k,iq)) < borne) CYCLE 59 57 WRITE(msg1,'(s,"(",i0,",",i0,",",i0,") = ",ES12.4)')TRIM(isoName(ixt)),i,k,iq,q(i,k,iq) 60 58 CALL msg(msg1, modname) -
LMDZ6/trunk/libf/dyn3d/dynetat0.F90
r5200 r5214 7 7 !------------------------------------------------------------------------------- 8 8 USE infotrac, ONLY: nqtot, tracers, niso, iqIsoPha, iH2O, isoName, & 9 new2oldH2O, newHNO3, oldHNO3 , getKey9 new2oldH2O, newHNO3, oldHNO3 10 10 USE strings_mod, ONLY: maxlen, msg, strStack, real2str, int2str 11 11 USE netcdf, ONLY: NF90_OPEN, NF90_NOWRITE, NF90_INQ_VARID, & … … 24 24 USE ioipsl_getincom, ONLY: getin 25 25 #endif 26 USE iso_params_mod ! tnat_* and alpha_ideal_* 26 27 27 28 IMPLICIT NONE … … 166 167 CALL msg(' !!! Beware: alpha_ideal put to 1 !!!', modname) 167 168 ELSE 168 IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) & 169 CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1) 169 SELECT CASE(isoName(iName)) 170 CASE('H216O'); tnat = tnat_H216O; alpha_ideal = alpha_ideal_H216O 171 CASE('H217O'); tnat = tnat_H217O; alpha_ideal = alpha_ideal_H217O 172 CASE('H218O'); tnat = tnat_H218O; alpha_ideal = alpha_ideal_H218O 173 CASE('HDO'); tnat = tnat_HDO; alpha_ideal = alpha_ideal_HDO 174 CASE('HTO'); tnat = tnat_HTO; alpha_ideal = alpha_ideal_HTO 175 CASE DEFAULT; CALL abort_gcm(TRIM(modname),'unknown isotope "'//TRIM(isoName(iName))//'" ; check tracer.def file',1) 176 END SELECT 170 177 END IF 171 178 CALL msg('Missing tracer <'//TRIM(var)//'> => initialized with a simplified Rayleigh distillation law.', modname) -
LMDZ6/trunk/libf/dyn3d/iniacademic.F90
r5200 r5214 5 5 6 6 USE filtreg_mod, ONLY: inifilr 7 USE infotrac, ONLY: nqtot, niso, iqIsoPha, tracers, getKey,isoName, addPhase7 USE infotrac, ONLY: nqtot, niso, iqIsoPha, tracers, isoName, addPhase 8 8 USE control_mod, ONLY: day_step,planet_type 9 9 use exner_hyb_m, only: exner_hyb … … 23 23 use netcdf, only : NF90_NOWRITE,NF90_OPEN,NF90_NOERR,NF90_INQ_VARID 24 24 use netcdf, only : NF90_CLOSE, NF90_GET_VAR 25 USE iso_params_mod ! tnat_* and alpha_ideal_* 25 26 26 27 … … 327 328 WRITE(lunout, *)'In '//TRIM(modname)//': !!! Beware: alpha_ideal put to 1 !!!' 328 329 ELSE 329 IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) & 330 CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1) 330 SELECT CASE(isoName(iName)) 331 CASE('H216O'); tnat = tnat_H216O; alpha_ideal = alpha_ideal_H216O 332 CASE('H217O'); tnat = tnat_H217O; alpha_ideal = alpha_ideal_H217O 333 CASE('H218O'); tnat = tnat_H218O; alpha_ideal = alpha_ideal_H218O 334 CASE('HDO'); tnat = tnat_HDO; alpha_ideal = alpha_ideal_HDO 335 CASE('HTO'); tnat = tnat_HTO; alpha_ideal = alpha_ideal_HTO 336 CASE DEFAULT 337 CALL abort_gcm(TRIM(modname),'unknown isotope "'//TRIM(isoName(iName))//'" ; check tracer.def file',1) 338 END SELECT 331 339 END IF 332 340 q(:,:,iq) = q(:,:,iqParent)*tnat*(q(:,:,iqParent)/30.e-3)**(alpha_ideal-1.) -
LMDZ6/trunk/libf/dyn3d_common/iso_verif_dyn.F90
r5213 r5214 1 function iso_verif_noNaN_nostop(x,err_msg) 2 implicit none 3 ! si x est NaN, on affiche message 4 ! d'erreur et return 1 si erreur 5 6 ! input: 7 real x 8 character*(*) err_msg ! message d''erreur à afficher 9 10 ! output 11 real borne 12 parameter (borne=1e19) 13 integer iso_verif_noNaN_nostop 14 15 if ((x.gt.-borne).and.(x.lt.borne)) then 16 iso_verif_noNAN_nostop=0 17 else 18 write(*,*) 'erreur detectee par iso_verif_nonNaN:' 19 write(*,*) err_msg 20 write(*,*) 'x=',x 21 iso_verif_noNaN_nostop=1 22 endif 23 24 return 25 end 26 27 function iso_verif_egalite_nostop 28 : (a,b,err_msg) 29 implicit none 30 ! compare a et b. Si pas egal, on affiche message 31 ! d'erreur et stoppe 32 ! pour egalite, on verifie erreur absolue et arreur relative 33 34 ! input: 35 real a, b 36 character*(*) err_msg ! message d''erreur à afficher 37 38 ! locals 39 real errmax ! erreur maximale en absolu. 40 real errmaxrel ! erreur maximale en relatif autorisée 41 parameter (errmax=1e-8) 42 parameter (errmaxrel=1e-3) 43 44 ! output 45 integer iso_verif_egalite_nostop 46 47 iso_verif_egalite_nostop=0 48 49 if (abs(a-b).gt.errmax) then 50 if (abs((a-b)/max(max(abs(b),abs(a)),1e-18)) 51 : .gt.errmaxrel) then 52 write(*,*) 'erreur detectee par iso_verif_egalite:' 53 write(*,*) err_msg 54 write(*,*) 'a=',a 55 write(*,*) 'b=',b 56 iso_verif_egalite_nostop=1 57 endif 58 endif 59 60 return 61 end 1 LOGICAL FUNCTION iso_verif_noNaN_nostop(x, err_msg) RESULT(lerr) 2 IMPLICIT NONE 3 ! If "x" is NaN, print an error message . 4 REAL, INTENT(IN) :: x 5 CHARACTER(LEN=*), INTENT(IN) :: err_msg 6 REAL, PARAMETER :: borne = 1e19 7 lerr = x > -borne .AND. x < borne 8 IF(.NOT.lerr) RETURN 9 WRITE(*,*) 'erreur detectee par iso_verif_nonNaN:' 10 WRITE(*,*) err_msg 11 WRITE(*,*) 'x=',x 12 END FUNCTION iso_verif_noNaN_nostop 62 13 63 14 64 function iso_verif_aberrant_nostop 65 : (x,iso,q,err_msg) 15 LOGICAL FUNCTION iso_verif_egalite_nostop(a, b, err_msg) RESULT(lerr) 16 IMPLICIT NONE 17 ! Compare "a" and "b". If a/=b, print an error message. 18 ! Both absolute and relative errors are checked for equality. 19 REAL, INTENT(IN) :: a, b 20 CHARACTER(LEN=*), INTENT(IN) :: err_msg 21 REAL, PARAMETER :: errmax = 1e-8, & ! max absolute error 22 errmaxrel = 1e-3 ! max relative error 23 lerr = ABS(a-b) > errmax 24 IF(.NOT.lerr) RETURN 25 lerr = ABS( (a-b) / MAX(MAX(ABS(b), ABS(a)),1e-18) ) > errmaxrel 26 IF(.NOT.lerr) RETURN 27 WRITE(*,*) 'erreur detectee par iso_verif_egalite:' 28 WRITE(*,*) err_msg 29 WRITE(*,*) 'a=',a 30 WRITE(*,*) 'b=',b 31 END FUNCTION iso_verif_egalite_nostop 32 33 34 LOGICAL FUNCTION iso_verif_aberrant_nostop(x, iso, q, err_msg) RESULT(lerr) 66 35 #ifdef CPP_IOIPSL 67 36 USE IOIPSL, ONLY: getin 68 37 #else 69 38 USE ioipsl_getincom, ONLY: getin 70 39 #endif 71 USE infotrac, ONLY: isoName, getKey 72 implicit none 73 74 ! input: 75 real x,q 76 integer iso ! 2=HDO, 1=O18 77 character*(*) err_msg ! message d''erreur à afficher 40 USE iso_params_mod, ONLY: tnat_HDO 41 IMPLICIT NONE 42 REAL, INTENT(IN) :: x, q 43 INTEGER, INTENT(IN) :: iso ! 2=HDO, 1=O18 44 CHARACTER(LEN=*), INTENT(IN) :: err_msg 78 45 79 ! locals 80 real qmin,deltaD 81 real deltaDmax,deltaDmin,tnat 82 parameter (qmin=1e-11) 83 parameter (deltaDmax=200.0,deltaDmin=-999.9) 84 LOGICAL, SAVE :: ltnat1 85 LOGICAL, SAVE :: lFirst=.TRUE. 46 REAL, PARAMETER :: qmin = 1e-11, & 47 deltaDmax = 200.0, & 48 deltaDmin =-999.9 49 LOGICAL :: ltnat1 50 LOGICAL, SAVE :: lFirst=.TRUE. 51 REAL, SAVE :: tnat 52 REAL :: deltaD 53 IF(lFirst) THEN 54 ltnat1 = .TRUE.; CALL getin('tnateq1', ltnat1) 55 tnat = tnat_HDO; IF(ltnat1) tnat = 1.0 56 lFirst = .FALSE. 57 END IF 58 lerr = q > qmin 59 IF(.NOT.lerr) RETURN 60 deltaD = (x / q /tnat - 1.) * 1000. 61 lerr = deltaD > deltaDmax .OR. deltaD < deltaDmin 62 IF(.NOT.lerr) RETURN 63 WRITE(*,*) 'erreur detectee par iso_verif_aberrant:' 64 WRITE(*,*) err_msg 65 WRITE(*,*) 'q=',q 66 WRITE(*,*) 'deltaD=',deltaD 67 WRITE(*,*) 'iso=',iso 68 END FUNCTION iso_verif_aberrant_nostop 86 69 87 ! output88 integer iso_verif_aberrant_nostop89 90 IF(lFirst) THEN91 ltnat1 = .TRUE.; CALL getin('tnateq1', ltnat1)92 lFirst = .FALSE.93 END IF94 iso_verif_aberrant_nostop=095 96 ! verifier que HDO est raisonable97 if (q.gt.qmin) then98 IF(ltnat1) THEN99 tnat = 1.0100 ELSE IF(getKey('tnat', tnat, isoName(iso))) THEN101 err_msg = 'Missing isotopic parameter "tnat"'102 iso_verif_aberrant_nostop=1103 RETURN104 END IF105 deltaD=(x/q/tnat-1)*1000106 if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then107 write(*,*) 'erreur detectee par iso_verif_aberrant:'108 write(*,*) err_msg109 write(*,*) 'q=',q110 write(*,*) 'deltaD=',deltaD111 write(*,*) 'iso=',iso112 iso_verif_aberrant_nostop=1113 endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then114 endif !if (q(i,k,iq).gt.qmin) then115 116 117 return118 end119 -
LMDZ6/trunk/libf/dyn3dmem/check_isotopes_loc.F90
r5200 r5214 3 3 USE strings_mod, ONLY: maxlen, msg, strIdx, strStack, int2str, real2str 4 4 USE infotrac, ONLY: nqtot, niso, nphas, isotope, isoCheck, iqIsoPha, isoSelect, & 5 ntiso, iH2O, nzone, tracers, isoName, itZonIso, getKey 5 ntiso, iH2O, nzone, tracers, isoName, itZonIso 6 USE iso_params_mod, ONLY: tnat_H216O, tnat_H217O, tnat_H218O, tnat_HDO, tnat_HTO 6 7 #ifdef CPP_IOIPSL 7 8 USE ioipsl, ONLY: getin … … 16 17 CHARACTER(LEN=maxlen) :: modname, msg1, nm(2) 17 18 INTEGER :: ixt, ipha, k, i, iq, iiso, izon, ieau, iqeau, iqpar 18 INTEGER, ALLOCATABLE :: ix(:)19 INTEGER, ALLOCATABLE :: ix(:) 19 20 REAL, ALLOCATABLE, SAVE :: tnat(:) !--- OpenMP shared variable 20 REAL 21 REAL :: xtractot, xiiso, deltaD, q1, q2 21 22 REAL, PARAMETER :: borne = 1e19, & 22 23 errmax = 1e-8, & !--- Max. absolute error … … 40 41 ltnat1 = .TRUE.; CALL getin('tnateq1', ltnat1) 41 42 ALLOCATE(tnat(niso)) 42 iso_eau = strIdx(isoName,'H216O') 43 iso_O17 = strIdx(isoName,'H217O') 44 iso_O18 = strIdx(isoName,'H218O') 45 iso_HDO = strIdx(isoName,'HDO') 46 iso_HTO = strIdx(isoName,'HTO') 47 IF(ltnat1) THEN 48 tnat(:)=1.0 49 ELSE 50 IF(getKey('tnat', tnat)) CALL abort_gcm(modname, 'missing isotopic parameter', 1) 51 END IF 43 iso_eau = strIdx(isoName,'H216O'); IF(iso_eau /= 0) tnat(iso_eau) = tnat_H216O 44 iso_O17 = strIdx(isoName,'H217O'); IF(iso_O17 /= 0) tnat(iso_O17) = tnat_H217O 45 iso_O18 = strIdx(isoName,'H218O'); IF(iso_O18 /= 0) tnat(iso_O18) = tnat_H218O 46 iso_HDO = strIdx(isoName,'HDO'); IF(iso_HDO /= 0) tnat(iso_HDO) = tnat_HDO 47 iso_HTO = strIdx(isoName,'HTO'); IF(iso_HTO /= 0) tnat(iso_HTO) = tnat_HTO 48 IF(ltnat1) tnat(:) = 1.0 52 49 !$OMP END MASTER 53 50 !$OMP BARRIER … … 64 61 DO k = 1, llm 65 62 DO i = ijb, ije 66 IF(ABS(q(i,k,iq)) < =borne) CYCLE63 IF(ABS(q(i,k,iq)) < borne) CYCLE 67 64 WRITE(msg1,'(s,"(",i0,",",i0,",",i0,") = ",ES12.4)')TRIM(isoName(ixt)),i,k,iq,q(i,k,iq) 68 65 CALL msg(msg1, modname) -
LMDZ6/trunk/libf/dyn3dmem/dynetat0_loc.F90
r5200 r5214 8 8 USE parallel_lmdz 9 9 USE infotrac, ONLY: nqtot, tracers, niso, iqIsoPha, iH2O, isoName, & 10 new2oldH2O, newHNO3, oldHNO3 , getKey10 new2oldH2O, newHNO3, oldHNO3 11 11 USE strings_mod, ONLY: maxlen, msg, strStack, real2str, int2str, strIdx 12 12 USE netcdf, ONLY: NF90_OPEN, NF90_NOWRITE, NF90_INQUIRE_DIMENSION, NF90_INQ_VARID, & … … 25 25 USE ioipsl_getincom, ONLY: getin 26 26 #endif 27 USE iso_params_mod ! tnat_* and alpha_ideal_* 27 28 28 29 IMPLICIT NONE … … 191 192 CALL msg(' !!! Beware: alpha_ideal put to 1 !!!', modname) 192 193 ELSE 193 IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) & 194 CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1) 194 SELECT CASE(isoName(iName)) 195 CASE('H216O'); tnat = tnat_H216O; alpha_ideal = alpha_ideal_H216O 196 CASE('H217O'); tnat = tnat_H217O; alpha_ideal = alpha_ideal_H217O 197 CASE('H218O'); tnat = tnat_H218O; alpha_ideal = alpha_ideal_H218O 198 CASE('HDO'); tnat = tnat_HDO; alpha_ideal = alpha_ideal_HDO 199 CASE('HTO'); tnat = tnat_HTO; alpha_ideal = alpha_ideal_HTO 200 CASE DEFAULT; CALL abort_gcm(TRIM(modname),'unknown isotope "'//TRIM(isoName(iName))//'" ; check tracer.def file',1) 201 END SELECT 195 202 END IF 196 203 CALL msg('Missing tracer <'//TRIM(var)//'> => initialized with a simplified Rayleigh distillation law.', modname) -
LMDZ6/trunk/libf/dyn3dmem/iniacademic_loc.F90
r5200 r5214 5 5 6 6 USE filtreg_mod, ONLY: inifilr 7 USE infotrac, ONLY: nqtot, niso, iqIsoPha, tracers, getKey,isoName, addPhase7 USE infotrac, ONLY: nqtot, niso, iqIsoPha, tracers, isoName, addPhase 8 8 USE control_mod, ONLY: day_step,planet_type 9 9 use exner_hyb_m, only: exner_hyb … … 24 24 use netcdf, only : NF90_NOWRITE,NF90_OPEN,NF90_NOERR,NF90_INQ_VARID 25 25 use netcdf, only : NF90_CLOSE, NF90_GET_VAR 26 USE iso_params_mod ! tnat_* and alpha_ideal_* 26 27 27 28 … … 329 330 WRITE(lunout, *) 'In '//TRIM(modname)//': !!! Beware: alpha_ideal put to 1 !!!' 330 331 ELSE 331 IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) & 332 CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1) 332 SELECT CASE(isoName(iName)) 333 CASE('H216O'); tnat = tnat_H216O; alpha_ideal = alpha_ideal_H216O 334 CASE('H217O'); tnat = tnat_H217O; alpha_ideal = alpha_ideal_H217O 335 CASE('H218O'); tnat = tnat_H218O; alpha_ideal = alpha_ideal_H218O 336 CASE('HDO'); tnat = tnat_HDO; alpha_ideal = alpha_ideal_HDO 337 CASE('HTO'); tnat = tnat_HTO; alpha_ideal = alpha_ideal_HTO 338 CASE DEFAULT 339 CALL abort_gcm(TRIM(modname),'unknown isotope "'//TRIM(isoName(iName))//'" ; check tracer.def file',1) 340 END SELECT 333 341 END IF 334 342 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.) -
LMDZ6/trunk/libf/phylmdiso/isotopes_mod.F90
r5200 r5214 4 4 MODULE isotopes_mod 5 5 USE strings_mod, ONLY: msg, real2str, int2str, bool2str, maxlen, strIdx, strStack 6 USE infotrac_phy, ONLY: isoName 6 USE infotrac_phy, ONLY: isoName, niso, ntiso 7 USE iso_params_mod 7 8 IMPLICIT NONE 8 9 INTERFACE get_in; MODULE PROCEDURE getinp_s, getinp_i, getinp_r, getinp_l; END INTERFACE get_in 9 SAVE10 10 11 11 !--- Contains all isotopic variables + their initialization 12 12 !--- Isotopes-specific routines are in isotopes_routines_mod to avoid circular dependencies with isotopes_verif_mod. 13 14 !--- Negligible lower thresholds: no need to check for absurd values under these lower limits 15 REAL, PARAMETER :: & 16 ridicule = 1e-12, & ! For mixing ratios 17 ridicule_rain = 1e-8, & ! For rain fluxes (rain, zrfl...) in kg/s <-> 1e-3 mm/day 18 ridicule_evap = ridicule_rain*1e-2, & ! For evaporations in kg/s <-> 1e-3 mm/day 19 ridicule_qsol = ridicule_rain, & ! For qsol in kg <-> 1e-8 kg 20 ridicule_snow = ridicule_qsol ! For snow in kg <-> 1e-8 kg 21 REAL, PARAMETER :: expb_max = 30.0 22 23 !--- Fractionation coefficients for H217O 24 REAL, PARAMETER :: fac_coeff_eq17_liq = 0.529, & 25 fac_coeff_eq17_ice = 0.529 26 27 !--- H218O reference 28 REAL, PARAMETER :: fac_enrichoce18 = 0.0005, alpha_liq_sol_O18 = 1.00291, & 29 talph1_O18 = 1137., talps1_O18 = 11.839, tkcin0_O18 = 0.006, & 30 talph2_O18 = -0.4156, talps2_O18 = -0.028244, tkcin1_O18 = 0.000285, & 31 talph3_O18 = -2.0667E-3, tdifrel_O18 = 1./0.9723, tkcin2_O18 = 0.00082 32 33 !--- Parameters that do not depend on the nature of water isotopes: 34 REAL, PARAMETER :: pxtmelt= 273.15, & !--- temperature at which ice formation starts 35 pxtice = 273.15 - 10.0, & !--- temperature at which all condensate is ice: 36 pxtmin = 273.15 - 120.0, & !--- computation done only under -120°C 37 pxtmax = 273.15 + 60.0, & !--- computation done only over +60°C 38 tdifexp= 0.58, & !--- a constant for alpha_eff for equilibrium below cloud base: 39 tv0cin = 7.0, & !--- wind threshold (m/s) for smooth/rough regime (Merlivat and Jouzel 1979) 40 musi = 1.0, & !--- facteurs lambda et mu dans Si=musi-lambda*T 41 Kd = 2.5e-9, & !--- diffusion in soil ; m2/s 42 rh_cste_surf_cond = 0.6, & !--- cste_surf_cond case: rhs and/or Ts set to constants 43 T_cste_surf_cond = 288.0 13 44 14 45 !--- Isotopes indices (in [1,niso] ; non-existing => 0 index) … … 89 120 !--- Vectors of length "niso" 90 121 REAL, ALLOCATABLE, DIMENSION(:), SAVE :: & 91 tnat, toce, tcorr, tdifrel92 !$OMP THREADPRIVATE( tnat, toce, tcorr, tdifrel)122 alpha, tnat, toce, tcorr, tdifrel 123 !$OMP THREADPRIVATE(alpha, tnat, toce, tcorr, tdifrel) 93 124 REAL, ALLOCATABLE, DIMENSION(:), SAVE :: & 94 125 talph1, talph2, talph3, talps1, talps2 … … 100 131 alpha_liq_sol, Rdefault, Rmethox 101 132 !$OMP THREADPRIVATE(alpha_liq_sol, Rdefault, Rmethox) 102 ! REAL, SAVE :: fac_coeff_eq17_liq, fac_coeff_eq17_ice103 !!$OMP THREADPRIVATE(fac_coeff_eq17_liq, fac_coeff_eq17_ice)104 105 !--- H2[18]O reference106 REAL, PARAMETER :: fac_enrichoce18=0.0005107 REAL, PARAMETER :: alpha_liq_sol_O18=1.00291108 REAL, PARAMETER :: talph1_O18=1137.109 REAL, PARAMETER :: talph2_O18=-0.4156110 REAL, PARAMETER :: talph3_O18=-2.0667E-3111 REAL, PARAMETER :: talps1_O18=11.839112 REAL, PARAMETER :: talps2_O18=-0.028244113 REAL, PARAMETER :: tdifrel_O18=1./0.9723114 REAL, PARAMETER :: tkcin0_O18=0.006115 REAL, PARAMETER :: tkcin1_O18=0.000285116 REAL, PARAMETER :: tkcin2_O18=0.00082117 REAL, PARAMETER :: fac_coeff_eq17_liq=0.529118 REAL, PARAMETER :: fac_coeff_eq17_ice=0.529119 120 !---- Parameters that do not depend on the nature of water isotopes:121 REAL, PARAMETER :: pxtmelt = 273.15 ! temperature at which ice formation starts122 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°C124 REAL, PARAMETER :: pxtmax = 273.15 + 60.0 ! On ne calcule qu'au dessus de +60°C125 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*T128 REAL, PARAMETER :: Kd=2.5e-9 ! m2/s ! diffusion dans le sol129 REAL, PARAMETER :: rh_cste_surf_cond = 0.6 ! cas où cste_surf_cond: on met rhs ou/et Ts cste pour voir130 REAL, PARAMETER :: T_cste_surf_cond = 288.0131 132 133 !--- Negligible lower thresholds: no need to check for absurd values under these lower limits134 REAL, PARAMETER :: &135 ridicule = 1e-12, & ! For mixing ratios136 ridicule_rain = 1e-8, & ! For rain fluxes (rain, zrfl...) in kg/s <-> 1e-3 mm/day137 ridicule_evap = ridicule_rain*1e-2, & ! For evaporations in kg/s <-> 1e-3 mm/day138 ridicule_qsol = ridicule_rain, & ! For qsol in kg <-> 1e-8 kg139 ridicule_snow = ridicule_qsol ! For snow in kg <-> 1e-8 kg140 REAL, PARAMETER :: expb_max = 30.0141 133 142 134 !--- Specific to HTO: … … 155 147 156 148 SUBROUTINE iso_init() 157 USE infotrac_phy, ONLY: ntiso, niso, getKey158 USE strings_mod, ONLY: maxlen159 149 IMPLICIT NONE 160 150 161 151 !=== Local variables: 162 INTEGER :: ixt 152 INTEGER :: ixt, is 163 153 LOGICAL :: ltnat1 164 154 CHARACTER(LEN=maxlen) :: modname, sxt … … 185 175 CALL msg('64: niso = '//TRIM(int2str(niso)), modname) 186 176 187 !--- Init de ntracisoOR: on ecrasera en cas de traceurs de tagging isotopiques188 ! (nzone>0) si complications avec ORCHIDEE189 ntracisoOR = ntiso190 191 !--- Type of water isotopes:192 iso_eau = strIdx(isoName, 'H216O'); CALL msg('iso_eau='//int2str(iso_eau), modname)193 iso_HDO = strIdx(isoName, 'HDO'); CALL msg('iso_HDO='//int2str(iso_HDO), modname)194 iso_O18 = strIdx(isoName, 'H218O'); CALL msg('iso_O18='//int2str(iso_O18), modname)195 iso_O17 = strIdx(isoName, 'H217O'); CALL msg('iso_O17='//int2str(iso_O17), modname)196 iso_HTO = strIdx(isoName, 'HTO'); CALL msg('iso_HTO='//int2str(iso_HTO), modname)197 198 !--- Initialiaation: reading the isotopic parameters.199 CALL get_in('lambda', lambda_sursat, 0.004)200 CALL get_in('thumxt1', thumxt1, 0.75*1.2)201 CALL get_in('ntot', ntot, 20, .FALSE.)202 CALL get_in('h_land_ice', h_land_ice, 20., .FALSE.)203 CALL get_in('P_veg', P_veg, 1.0, .FALSE.)204 CALL get_in('bidouille_anti_divergence', bidouille_anti_divergence, .FALSE.)205 CALL get_in('essai_convergence', essai_convergence, .FALSE.)206 CALL get_in('initialisation_iso', initialisation_iso, 0)207 208 ! IF(nzone>0 .AND. initialisation_iso==0) &209 ! CALL get_in('initialisation_isotrac',initialisation_isotrac)210 CALL get_in('modif_sst', modif_sst, 0)211 CALL get_in('deltaTtest', deltaTtest, 0.0) !--- For modif_sst>=1212 CALL get_in('deltaTtestpoles',deltaTtestpoles, 0.0) !--- For modif_sst>=2213 CALL get_in( 'sstlatcrit', sstlatcrit, 30.0) !--- For modif_sst>=3214 CALL get_in('dsstlatcrit', dsstlatcrit, 0.0) !--- For modif_sst>=3177 !--- Init de ntracisoOR: on ecrasera en cas de traceurs de tagging isotopiques 178 ! (nzone>0) si complications avec ORCHIDEE 179 ntracisoOR = ntiso 180 181 !--- Type of water isotopes: 182 iso_eau = strIdx(isoName, 'H216O'); CALL msg('iso_eau='//int2str(iso_eau), modname) 183 iso_HDO = strIdx(isoName, 'HDO'); CALL msg('iso_HDO='//int2str(iso_HDO), modname) 184 iso_O18 = strIdx(isoName, 'H218O'); CALL msg('iso_O18='//int2str(iso_O18), modname) 185 iso_O17 = strIdx(isoName, 'H217O'); CALL msg('iso_O17='//int2str(iso_O17), modname) 186 iso_HTO = strIdx(isoName, 'HTO'); CALL msg('iso_HTO='//int2str(iso_HTO), modname) 187 188 !--- Initialisation: reading the isotopic parameters. 189 CALL get_in('lambda', lambda_sursat, 0.004); IF(ok_nocinsat) lambda_sursat = 0. 190 CALL get_in('thumxt1', thumxt1, 0.75*1.2) 191 CALL get_in('ntot', ntot, 20, .FALSE.) 192 CALL get_in('h_land_ice', h_land_ice, 20., .FALSE.) 193 CALL get_in('P_veg', P_veg, 1.0, .FALSE.) 194 CALL get_in('bidouille_anti_divergence', bidouille_anti_divergence, .FALSE.) 195 CALL get_in('essai_convergence', essai_convergence, .FALSE.) 196 CALL get_in('initialisation_iso', initialisation_iso, 0) 197 198 ! IF(nzone>0 .AND. initialisation_iso==0) & 199 ! CALL get_in('initialisation_isotrac',initialisation_isotrac) 200 CALL get_in('modif_sst', modif_sst, 0) 201 CALL get_in('deltaTtest', deltaTtest, 0.0) !--- For modif_sst>=1 202 CALL get_in('deltaTtestpoles',deltaTtestpoles, 0.0) !--- For modif_sst>=2 203 CALL get_in( 'sstlatcrit', sstlatcrit, 30.0) !--- For modif_sst>=3 204 CALL get_in('dsstlatcrit', dsstlatcrit, 0.0) !--- For modif_sst>=3 215 205 #ifdef ISOVERIF 216 CALL msg('iso_init 270: sstlatcrit='//real2str( sstlatcrit), modname, sstlatcrit < 0.0) !--- For modif_sst>=2217 CALL msg('iso_init 279: dsstlatcrit='//real2str(dsstlatcrit), modname, sstlatcrit < 0.0) !--- For modif_sst>=3218 IF(modif_sst >= 2 .AND. sstlatcrit < 0.0) STOP206 CALL msg('iso_init 270: sstlatcrit='//real2str( sstlatcrit), modname, sstlatcrit < 0.0) !--- For modif_sst>=2 207 CALL msg('iso_init 279: dsstlatcrit='//real2str(dsstlatcrit), modname, sstlatcrit < 0.0) !--- For modif_sst>=3 208 IF(modif_sst >= 2 .AND. sstlatcrit < 0.0) STOP 219 209 #endif 220 210 221 CALL get_in('modif_sic', modif_sic, 0) 222 IF(modif_sic >= 1) & 223 CALL get_in('deltasic', deltasic, 0.1) 224 225 CALL get_in('albedo_prescrit', albedo_prescrit, 0) 226 IF(albedo_prescrit == 1) THEN 227 CALL get_in('lon_min_albedo', lon_min_albedo, -200.) 228 CALL get_in('lon_max_albedo', lon_max_albedo, 200.) 229 CALL get_in('lat_min_albedo', lat_min_albedo, -100.) 230 CALL get_in('lat_max_albedo', lat_max_albedo, 100.) 231 END IF 232 CALL get_in('deltaP_BL', deltaP_BL, 10.0) 233 CALL get_in('ruissellement_pluie', ruissellement_pluie, 0) 234 CALL get_in('alphak_stewart', alphak_stewart, 1) 235 CALL get_in('tdifexp_sol', tdifexp_sol, 0.67) 236 CALL get_in('calendrier_guide', calendrier_guide, 0) 237 CALL get_in('cste_surf_cond', cste_surf_cond, 0) 238 CALL get_in('mixlen', mixlen, 35.0) 239 CALL get_in('evap_cont_cste', evap_cont_cste, 0) 240 CALL get_in('deltaO18_evap_cont', deltaO18_evap_cont,0.0) 241 CALL get_in('d_evap_cont', d_evap_cont, 0.0) 242 CALL get_in('nudge_qsol', nudge_qsol, 0) 243 CALL get_in('region_nudge_qsol', region_nudge_qsol, 1) 244 nlevmaxO17 = 50 245 CALL msg('nlevmaxO17='//TRIM(int2str(nlevmaxO17))) 246 CALL get_in('no_pce', no_pce, 0) 247 CALL get_in('A_satlim', A_satlim, 1.0) 248 CALL get_in('ok_restrict_A_satlim', ok_restrict_A_satlim, 0) 211 CALL get_in('modif_sic', modif_sic, 0) 212 IF(modif_sic >= 1) & 213 CALL get_in('deltasic', deltasic, 0.1) 214 215 CALL get_in('albedo_prescrit', albedo_prescrit, 0) 216 IF(albedo_prescrit == 1) THEN 217 CALL get_in('lon_min_albedo', lon_min_albedo, -200.) 218 CALL get_in('lon_max_albedo', lon_max_albedo, 200.) 219 CALL get_in('lat_min_albedo', lat_min_albedo, -100.) 220 CALL get_in('lat_max_albedo', lat_max_albedo, 100.) 221 END IF 222 CALL get_in('deltaO18_oce', deltaO18_oce, 0.0) 223 CALL get_in('deltaP_BL', deltaP_BL, 10.0) 224 CALL get_in('ruissellement_pluie', ruissellement_pluie, 0) 225 CALL get_in('alphak_stewart', alphak_stewart, 1) 226 CALL get_in('tdifexp_sol', tdifexp_sol, 0.67) 227 CALL get_in('calendrier_guide', calendrier_guide, 0) 228 CALL get_in('cste_surf_cond', cste_surf_cond, 0) 229 CALL get_in('mixlen', mixlen, 35.0) 230 CALL get_in('evap_cont_cste', evap_cont_cste, 0) 231 CALL get_in('deltaO18_evap_cont', deltaO18_evap_cont,0.0) 232 CALL get_in('d_evap_cont', d_evap_cont, 0.0) 233 CALL get_in('nudge_qsol', nudge_qsol, 0) 234 CALL get_in('region_nudge_qsol', region_nudge_qsol, 1) 235 nlevmaxO17 = 50 236 CALL msg('nlevmaxO17='//TRIM(int2str(nlevmaxO17))) 237 CALL get_in('no_pce', no_pce, 0) 238 CALL get_in('A_satlim', A_satlim, 1.0) 239 CALL get_in('ok_restrict_A_satlim', ok_restrict_A_satlim, 0) 249 240 #ifdef ISOVERIF 250 CALL msg(' 315: A_satlim='//real2str(A_satlim), modname, A_satlim > 1.0)251 IF(A_satlim > 1.0) STOP241 CALL msg(' 315: A_satlim='//real2str(A_satlim), modname, A_satlim > 1.0) 242 IF(A_satlim > 1.0) STOP 252 243 #endif 253 ! CALL get_in('slope_limiterxy', slope_limiterxy, 2.0) 254 ! CALL get_in('slope_limiterz', slope_limiterz, 2.0) 255 CALL get_in('modif_ratqs', modif_ratqs, 0) 256 CALL get_in('Pcrit_ratqs', Pcrit_ratqs, 500.0) 257 CALL get_in('ratqsbasnew', ratqsbasnew, 0.05) 258 CALL get_in('fac_modif_evaoce', fac_modif_evaoce, 1.0) 259 CALL get_in('ok_bidouille_wake', ok_bidouille_wake, 0) 260 ! si oui, la temperature de cond est celle de l'environnement, pour eviter 261 ! bugs quand temperature dans ascendances convs est mal calculee 262 CALL get_in('cond_temp_env', cond_temp_env, .FALSE.) 263 IF(ANY(isoName == 'HTO')) & 264 CALL get_in('ok_prod_nucl_tritium', ok_prod_nucl_tritium, .FALSE., .FALSE.) 265 CALL get_in('tnateq1', ltnat1, .TRUE.) 266 267 ! Ocean composition 268 CALL get_in('deltaO18_oce', deltaO18_oce, 0.0) 269 270 CALL msg('iso_O18, iso_HDO, iso_eau = '//TRIM(strStack(int2str([iso_O18, iso_HDO, iso_eau]))), modname) 271 272 !-------------------------------------------------------------- 273 ! Isotope fractionation factors and a few isotopic constants 274 !-------------------------------------------------------------- 275 ALLOCATE(tkcin0(niso)) 276 ALLOCATE(tkcin1(niso)) 277 ALLOCATE(tkcin2(niso)) 278 ALLOCATE(tnat(niso)) 279 ALLOCATE(tdifrel(niso)) 280 ALLOCATE(toce(niso)) 281 ALLOCATE(tcorr(niso)) 282 ALLOCATE(talph1(niso)) 283 ALLOCATE(talph2(niso)) 284 ALLOCATE(talph3(niso)) 285 ALLOCATE(talps1(niso)) 286 ALLOCATE(talps2(niso)) 287 ALLOCATE(alpha_liq_sol(niso)) 288 ALLOCATE(Rdefault(niso)) 289 ALLOCATE(Rmethox(niso)) 290 291 do ixt=1,niso 292 if (ixt.eq.iso_HTO) then ! Tritium 293 tkcin0(ixt) = 0.01056 294 tkcin1(ixt) = 0.0005016 295 tkcin2(ixt) = 0.0014432 296 tnat(ixt) = 0.0; IF(ltnat1) tnat(ixt)=1 297 toce(ixt)=4.0E-19 ! rapport T/H = 0.2 TU Dreisigacker and Roether 1978 298 tcorr(ixt)=1. 299 tdifrel(ixt)=1./0.968 300 talph1(ixt)=46480. 301 talph2(ixt)=-103.87 302 talph3(ixt)=0. 303 talps1(ixt)=46480. 304 talps2(ixt)=-103.87 305 alpha_liq_sol(ixt)=1. 306 Rmethox(ixt)=0.0 307 endif 308 if (ixt.eq.iso_O17) then ! O17 309 pente_MWL=0.528 310 tdifrel(ixt)=1./0.98555 ! valeur utilisée en 1D et dans modèle de LdG ! tdifrel(ixt)=1./0.985452 ! donné par Amaelle 311 fac_kcin= (tdifrel(ixt)-1.0)/(tdifrel_O18-1.0) ! fac_kcin=0.5145 ! donné par Amaelle 312 tkcin0(ixt) = tkcin0_O18*fac_kcin 313 tkcin1(ixt) = tkcin1_O18*fac_kcin 314 tkcin2(ixt) = tkcin2_O18*fac_kcin 315 tnat(ixt)=0.004/100. ! O17 représente 0.004% de l'oxygène 316 IF(ltnat1) tnat(ixt)=1 317 toce(ixt)=tnat(ixt)*(1.0+deltaO18_oce/1000.0)**pente_MWL 318 tcorr(ixt)=1.0+fac_enrichoce18*pente_MWL ! donné par Amaelle 319 talph1(ixt)=talph1_O18 320 talph2(ixt)=talph2_O18 321 talph3(ixt)=talph3_O18 322 talps1(ixt)=talps1_O18 323 talps2(ixt)=talps2_O18 324 alpha_liq_sol(ixt)=(alpha_liq_sol_O18)**fac_coeff_eq17_liq 325 Rdefault(ixt)=tnat(ixt)*(-3.15/1000.0+1.0) 326 Rmethox(ixt)=(230./1000.+1.)*tnat(ixt) !Zahn et al 2006 327 endif 328 if (ixt.eq.iso_O18) then ! Oxygene18 329 tkcin0(ixt) = tkcin0_O18 330 tkcin1(ixt) = tkcin1_O18 331 tkcin2(ixt) = tkcin2_O18 332 tnat(ixt)=2005.2E-6; IF(ltnat1) tnat(ixt)=1 333 toce(ixt)=tnat(ixt)*(1.0+deltaO18_oce/1000.0) 334 tcorr(ixt)=1.0+fac_enrichoce18 335 tdifrel(ixt)=tdifrel_O18 336 talph1(ixt)=talph1_O18 337 talph2(ixt)=talph2_O18 338 talph3(ixt)=talph3_O18 339 talps1(ixt)=talps1_O18 340 talps2(ixt)=talps2_O18 341 alpha_liq_sol(ixt)=alpha_liq_sol_O18 342 Rdefault(ixt)=tnat(ixt)*(-6.0/1000.0+1.0) 343 Rmethox(ixt)=(130./1000.+1.)*tnat(ixt) !Zahn et al 2006 344 endif 345 if (ixt.eq.iso_HDO) then ! Deuterium 346 pente_MWL=8.0 347 tdifrel(ixt)=1./0.9755 ! fac_kcin=0.88 348 fac_kcin= (tdifrel(ixt)-1)/(tdifrel_O18-1) 349 tkcin0(ixt) = tkcin0_O18*fac_kcin 350 tkcin1(ixt) = tkcin1_O18*fac_kcin 351 tkcin2(ixt) = tkcin2_O18*fac_kcin 352 tnat(ixt)=155.76E-6; IF(ltnat1) tnat(ixt)=1 353 toce(ixt)=tnat(ixt)*(1.0+pente_MWL*deltaO18_oce/1000.0) 354 tcorr(ixt)=1.0+fac_enrichoce18*pente_MWL 355 talph1(ixt)=24844. 356 talph2(ixt)=-76.248 357 talph3(ixt)=52.612E-3 358 talps1(ixt)=16288. 359 talps2(ixt)=-0.0934 360 !ZXalpha_liq_sol=1.0192 ! Weston, Ralph, 1955 361 alpha_liq_sol(ixt)=1.0212 362 ! valeur de Lehmann & Siegenthaler, 1991, Journal of 363 ! Glaciology, vol 37, p 23 364 Rdefault(ixt)=tnat(ixt)*((-6.0*pente_MWL+10.0)/1000.0+1.0) 365 Rmethox(ixt)=tnat(ixt)*(-25.0/1000.+1.) ! Zahn et al 2006 366 endif 367 if (ixt.eq.iso_eau) then ! Oxygene16 368 tkcin0(ixt) = 0.0 369 tkcin1(ixt) = 0.0 370 tkcin2(ixt) = 0.0 371 tnat(ixt)=1. 372 toce(ixt)=tnat(ixt) 373 tcorr(ixt)=1.0 374 tdifrel(ixt)=1. 375 talph1(ixt)=0. 376 talph2(ixt)=0. 377 talph3(ixt)=0. 378 talps1(ixt)=0. 379 talph3(ixt)=0. 380 alpha_liq_sol(ixt)=1. 381 Rdefault(ixt)=tnat(ixt)*1.0 382 Rmethox(ixt)=1.0 383 endif 384 enddo ! ixt=1,niso 385 386 IF(.NOT.Rdefault_smow) then 387 Rdefault(:) = 0.0 388 if (iso_eau.gt.0) Rdefault(iso_eau) = 1.0 ! correction Camille 30 mars 2023 389 ENDIF 390 write(*,*) 'Rdefault=',Rdefault 391 write(*,*) 'toce=',toce 392 393 !--- Sensitivity test: no kinetic effect in sfc evaporation 394 IF(ok_nocinsfc) THEN 395 tkcin0(1:niso) = 0.0 396 tkcin1(1:niso) = 0.0 397 tkcin2(1:niso) = 0.0 398 END IF 399 400 CALL msg('285: verif initialisation:', modname) 401 DO ixt=1,niso 402 sxt=int2str(ixt) 403 CALL msg(' * isoName('//TRIM(sxt)//') = <'//TRIM(isoName(ixt))//'>', modname) 404 CALL msg( ' tnat('//TRIM(sxt)//') = '//TRIM(real2str(tnat(ixt))), modname) 405 ! CALL msg(' alpha_liq_sol('//TRIM(sxt)//') = '//TRIM(real2str(alpha_liq_sol(ixt))), modname) 406 ! CALL msg( ' tkcin0('//TRIM(sxt)//') = '//TRIM(real2str(tkcin0(ixt))), modname) 407 ! CALL msg( ' tdifrel('//TRIM(sxt)//') = '//TRIM(real2str(tdifrel(ixt))), modname) 408 END DO 409 CALL msg('69: lambda = '//TRIM(real2str(lambda_sursat)), modname) 410 CALL msg('69: thumxt1 = '//TRIM(real2str(thumxt1)), modname) 411 CALL msg('69: h_land_ice = '//TRIM(real2str(h_land_ice)), modname) 412 CALL msg('69: P_veg = '//TRIM(real2str(P_veg)), modname) 244 ! CALL get_in('slope_limiterxy', slope_limiterxy, 2.0) 245 ! CALL get_in('slope_limiterz', slope_limiterz, 2.0) 246 CALL get_in('modif_ratqs', modif_ratqs, 0) 247 CALL get_in('Pcrit_ratqs', Pcrit_ratqs, 500.0) 248 CALL get_in('ratqsbasnew', ratqsbasnew, 0.05) 249 CALL get_in('fac_modif_evaoce', fac_modif_evaoce, 1.0) 250 CALL get_in('ok_bidouille_wake', ok_bidouille_wake, 0) 251 ! si oui, la temperature de cond est celle de l'environnement, pour eviter 252 ! bugs quand temperature dans ascendances convs est mal calculee 253 CALL get_in('cond_temp_env', cond_temp_env, .FALSE.) 254 IF(ANY(isoName == 'HTO')) & 255 CALL get_in('ok_prod_nucl_tritium', ok_prod_nucl_tritium, .FALSE., .FALSE.) 256 CALL get_in('tnateq1', ltnat1, .TRUE.) 257 258 CALL msg('iso_O18, iso_HDO, iso_eau = '//TRIM(strStack(int2str([iso_O18, iso_HDO, iso_eau]))), modname) 259 260 !-------------------------------------------------------------- 261 ! Parameters that depend on the nature of water isotopes: 262 !-------------------------------------------------------------- 263 ALLOCATE(tnat (niso), talph1(niso), talps1(niso), tkcin0(niso), tdifrel (niso), alpha (niso)) 264 ALLOCATE(toce (niso), talph2(niso), talps2(niso), tkcin1(niso), Rdefault(niso), alpha_liq_sol(niso)) 265 ALLOCATE(tcorr(niso), talph3(niso), tkcin2(niso), Rmethox (niso)) 266 267 !=== H216O 268 is = iso_eau 269 IF(is /= 0) THEN 270 tdifrel (is) = 1.0 271 alpha (is) = alpha_ideal_H216O 272 tnat (is) = tnat_H216O; IF(ltnat1) tnat(is) = 1.0 273 toce (is) = tnat(is) 274 tcorr (is) = 1.0 275 talph1 (is) = 0.0; talps1(is) = 0.0; tkcin0(is) = 0.0 276 talph2 (is) = 0.0; talps2(is) = 0.0; tkcin1(is) = 0.0 277 talph3 (is) = 0.0; tkcin2(is) = 0.0 278 Rdefault(is) = tnat(is)*1.0 279 Rmethox (is) = 1.0 280 alpha_liq_sol(is) = 1.0 281 END IF 282 283 !=== H217O 284 is = iso_O17 285 IF(is /= 0) THEN; pente_MWL = 0.528 286 tdifrel (is) = 1./0.98555 ! used in 1D and in LdG model ; tdifrel=1./0.985452: from Amaelle 287 alpha (is) = alpha_ideal_H217O 288 tnat (is) = tnat_H217O; IF(ltnat1) tnat(is) = 1.0 289 toce (is) = tnat(is)*(1.0+deltaO18_oce/1000.0)**pente_MWL 290 tcorr (is) = 1.0+fac_enrichoce18*pente_MWL 291 fac_kcin = (tdifrel(is)-1.0)/(tdifrel_O18-1.0) ! fac_kcin=0.5145: from Amaelle 292 talph1 (is) = talph1_O18; talps1(is) = talps1_O18; tkcin0(is) = tkcin0_O18*fac_kcin 293 talph2 (is) = talph2_O18; talps2(is) = talps2_O18; tkcin1(is) = tkcin1_O18*fac_kcin 294 talph3 (is) = talph3_O18; tkcin2(is) = tkcin2_O18*fac_kcin 295 Rdefault(is) = tnat(is)*(1.0-3.15/1000.) 296 Rmethox (is) = tnat(is)*(1.0+230./1000.) 297 alpha_liq_sol(is) = alpha_liq_sol_O18**fac_coeff_eq17_liq 298 END IF 299 300 !=== H218O 301 is = iso_O18 302 IF(is /= 0) THEN 303 tdifrel (is) = tdifrel_O18 304 alpha (is) = alpha_ideal_H218O 305 tnat (is) = tnat_H218O; IF(ltnat1) tnat(is) = 1.0 306 toce (is) = tnat(is)*(1.0+deltaO18_oce/1000.0) 307 tcorr (is) = 1.0+fac_enrichoce18 308 talph1 (is) = talph1_O18; talps1(is) = talps1_O18; tkcin0(is) = tkcin0_O18 309 talph2 (is) = talph2_O18; talps2(is) = talps2_O18; tkcin1(is) = tkcin1_O18 310 talph3 (is) = talph3_O18; tkcin2(is) = tkcin2_O18 311 Rdefault(is) = tnat(is)*(1.0-6.00/1000.) 312 Rmethox (is) = tnat(is)*(1.0+130./1000.) ! Zahn & al. 2006 313 alpha_liq_sol(is) = alpha_liq_sol_O18 314 END IF 315 316 !=== HDO 317 is = iso_HDO 318 IF(is /= 0) THEN; pente_MWL = 8.0 319 tdifrel (is) = 1./0.9755 ! fac_kcin=0.88 320 alpha (is) = alpha_ideal_HDO 321 tnat (is) = tnat_HDO; IF(ltnat1) tnat(is) = 1.0 322 toce (is) = tnat(is)*(1.0+deltaO18_oce/1000.0*pente_MWL) 323 tcorr (is) = 1.0+fac_enrichoce18*pente_MWL 324 fac_kcin = (tdifrel(is)-1.0)/(tdifrel_O18-1.0) 325 talph1 (is) = 24844.; talps1(is) = 16288.; tkcin0(is) = tkcin0_O18*fac_kcin 326 talph2 (is) = -76.248; talps2(is) = -0.0934; tkcin1(is) = tkcin1_O18*fac_kcin 327 talph3 (is) = 52.612E-3; tkcin2(is) = tkcin2_O18*fac_kcin 328 Rdefault(is) = tnat(is)*(1.0+(10.0-6.0*pente_MWL)/1000.) 329 Rmethox (is) = tnat(is)*(1.0-25.0/1000.) 330 alpha_liq_sol(is) = 1.0212 ! Lehmann & Siegenthaler, 1991, Jo. of Glaciology, vol 37, p 23 331 ! alpha_liq_sol=1.0192: Weston, Ralph, 1955 332 END IF 333 334 !=== HTO 335 is = iso_HTO 336 IF(is /= 0) THEN 337 tdifrel (is) = 1./0.968 338 alpha (is) = alpha_ideal_HTO 339 tnat (is) = tnat_HTO; IF(ltnat1) tnat(is) = 1.0 340 toce (is) = 4.0E-19 ! ratio T/H = 0.2 TU Dreisigacker & Roether 1978 341 tcorr (is) = 1.0 342 talph1 (is) = 46480.; talps1(is) = 46480.; tkcin0(is) = 0.01056 343 talph2 (is) = -103.87; talps2(is) = -103.87; tkcin1(is) = 0.0005016 344 talph3 (is) = 0.0; tkcin2(is) = 0.0014432 345 Rdefault(is) = 0.0 346 Rmethox (is) = 0.0 347 alpha_liq_sol(is) = 1.0 348 END IF 349 350 IF(.NOT. Rdefault_smow) THEN 351 Rdefault(:) = 0.0; IF(iso_eau > 0) Rdefault(iso_eau) = 1.0 352 END IF 353 WRITE(*,*) 'Rdefault = ',Rdefault 354 WRITE(*,*) 'toce = ', toce 355 356 !--- Sensitivity test: no kinetic effect in sfc evaporation 357 IF(ok_nocinsfc) THEN 358 tkcin0(1:niso) = 0.0 359 tkcin1(1:niso) = 0.0 360 tkcin2(1:niso) = 0.0 361 END IF 362 363 CALL msg('285: verif initialisation:', modname) 364 DO ixt=1,niso 365 sxt=int2str(ixt) 366 CALL msg(' * isoName('//TRIM(sxt)//') = <'//TRIM(isoName(ixt))//'>', modname) 367 CALL msg( ' tnat('//TRIM(sxt)//') = '//TRIM(real2str(tnat(ixt))), modname) 368 ! CALL msg(' alpha_liq_sol('//TRIM(sxt)//') = '//TRIM(real2str(alpha_liq_sol(ixt))), modname) 369 ! CALL msg( ' tkcin0('//TRIM(sxt)//') = '//TRIM(real2str(tkcin0(ixt))), modname) 370 ! CALL msg( ' tdifrel('//TRIM(sxt)//') = '//TRIM(real2str(tdifrel(ixt))), modname) 371 END DO 372 CALL msg('69: lambda = '//TRIM(real2str(lambda_sursat)), modname) 373 CALL msg('69: thumxt1 = '//TRIM(real2str(thumxt1)), modname) 374 CALL msg('69: h_land_ice = '//TRIM(real2str(h_land_ice)), modname) 375 CALL msg('69: P_veg = '//TRIM(real2str(P_veg)), modname) 413 376 414 377 END SUBROUTINE iso_init
Note: See TracChangeset
for help on using the changeset viewer.