Changeset 4158 for LMDZ6


Ignore:
Timestamp:
May 19, 2022, 10:47:23 AM (3 years ago)
Author:
dcugnet
Message:
  • include new type "tracer.def" tracers description files in DefLists/?
  • Remove the old code (lOldCode=.TRUE.)
  • Remove the hard-coded part in isotopes_mod that was ensuring the convergence => this version will not converge with previous isotopic physics, but the

difference, only due to round off errors on "tdifrel", is acceptable.

Location:
LMDZ6/trunk
Files:
3 added
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/dyn3d_common/infotrac.F90

    r4154 r4158  
    132132   INTEGER,             SAVE, ALLOCATABLE :: conv_flg(:), &     !--- Convection     activation ; needed for INCA        (nbtr)
    133133                                              pbl_flg(:)        !--- Boundary layer activation ; needed for INCA        (nbtr)
    134    LOGICAL, PARAMETER :: lOldCode = .FALSE.
    135134
    136135CONTAINS
    … …  
    202201   
    203202   CALL msg('type_trac = "'//TRIM(type_trac)//'"', modname)
    204    IF(lOldCode) THEN
    205       str = [type_trac]; nt = 1
    206    ELSE
    207       IF(strParse(type_trac, ',', str, n=nt)) CALL abort_gcm(modname,'can''t parse "type_trac = '//TRIM(type_trac)//'"',1)
    208    END IF
     203   IF(strParse(type_trac, ',', str, n=nt)) CALL abort_gcm(modname,'can''t parse "type_trac = '//TRIM(type_trac)//'"',1)
    209204
    210205   !---------------------------------------------------------------------------------------------------------------------------
    … …  
    256251! 1) Get the numbers of: true (first order only) tracers "nqtrue", water tracers "nqo" (vapor/liquid/solid)
    257252!==============================================================================================================================
    258 
    259 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    260    IF(lOldCode) THEN         !--- "type_trac" is a single keyword => no need to loop on its parsed version "str(:)"
    261 !------------------------------------------------------------------------------------------------------------------------------
    262    !--- Determine nqtrue and (INCA only) nqo, nbtr
    263    OPEN(90, FILE='traceur.def', FORM='formatted', STATUS='old', IOSTAT=ierr)
    264    IF(ierr /= 0) CALL abort_gcm(modname, 'file "traceur.def" not found !', 1)
    265    CALL msg('File "traceur.def" successfully opened.', modname)
    266 
    267    IF(ANY(['inca','inco'] == type_trac)) THEN
    268 #ifdef INCA
    269       READ(90,*) nqo
    270       IF(nqCO2==1 .AND. nqo==4) nqo = 3
    271       CALL Init_chem_inca_trac(nqINCA)
    272       nbtr = nqINCA + nqCO2
    273       nqtrue = nbtr + nqo
    274       IF(ALL([2,3] /= nqo)) CALL abort_gcm(modname, 'Only 2 or 3 water phases allowed ; found nqo='//TRIM(int2str(nqo)), 1)
    275       ALLOCATE(hadv_inca(nqINCA), conv_flg_inca(nqINCA), solsym_inca(nqINCA))
    276       ALLOCATE(vadv_inca(nqINCA),  pbl_flg_inca(nqINCA))
    277       CALL init_transport(solsym_inca, conv_flg_inca, pbl_flg_inca, hadv_inca, vadv_inca)
    278 #endif
    279    ELSE
    280       READ(90,*) nqtrue
    281    END IF
    282 
    283    IF (planet_type=="earth" .AND. nqtrue < 2) &
    284       CALL abort_gcm('infotrac_init', 'Not enough tracers: nqtrue='//TRIM(int2str(nqtrue))//', 2 tracers is the minimum', 1)
    285 
    286    !--- Allocate variables depending on nqtrue
    287    ALLOCATE(hadv(nqtrue), vadv(nqtrue), tnom_0(nqtrue), tnom_transp(nqtrue), tracers(nqtrue))
    288 
    289    !--- Continue to read tracer.def
    290    it = 0
    291    DO iq = 1, nqtrue
    292 #ifdef INCA
    293       IF(iq > nqo+nqCO2) THEN
    294          it = it+1
    295          hadv  (iq) = hadv_inca  (it)
    296          vadv  (iq) = vadv_inca  (it)
    297          tnom_0(iq) = solsym_inca(it)
    298          tnom_transp(iq) = 'air'
    299          CYCLE
    300       END IF
    301 #endif
    302       READ(90,'(I2,X,I2,X,A)',IOSTAT=ierr) hadv(iq),vadv(iq),tchaine
    303       IF(ierr/=0) CALL abort_gcm('infotrac_init', 'Pb dans la lecture de traceur.def', 1)
    304       jq = INDEX(tchaine(1:LEN_TRIM(tchaine)),' ')
    305       CALL msg("Ancienne version de traceur.def: traceurs d'air uniquement", modname, iq==1 .AND. jq==0)
    306       CALL msg("Nouvelle version de traceur.def",                            modname, iq==1 .AND. jq/=0)
    307       CALL msg('iq, hadv, vadv, tchaine ='//TRIM(strStack(int2str([iq, hadv(iq), vadv(iq)])))//', '//TRIM(tchaine), modname)
    308       IF(jq /= 0) THEN                                               !--- Space in the string chain => new format
    309          tnom_0     (iq) = tchaine(1:jq-1)
    310          tnom_transp(iq) = tchaine(jq+1:)
    311       ELSE
    312          tnom_0     (iq) = tchaine
    313          tnom_transp(iq) = 'air'
    314       END IF
    315    END DO
    316    CLOSE(90)
    317 
    318 #ifndef INCA
    319    IF(planet_type /= 'earth') nqo = 0                                !--- Same number of tracers in dynamics and physics
    320    IF(planet_type == 'earth') nqo = COUNT(delPhase(tnom_0) == 'H2O') !--- for all planets except for Earth
    321    nbtr = nqtrue - nqo
    322 #endif
    323 
    324    CALL msg('RAW CONTENT OF "traceur.def" FILE:', modname)
    325    IF(dispTable('iiss', ['hadv  ', 'vadv  ', 'name  ', 'parent'], cat(tnom_0, tnom_transp), cat(hadv, vadv))) &
    326       CALL abort_gcm(modname, "problem with the tracers table content", 1)
    327 
    328    !--- SET FIELDS %name, %parent, %phase, %component
    329    tracers(:)%name      = old2newName(tnom_0)
    330    tracers(:)%parent    = old2newName(tnom_transp)
    331    tracers(:)%phase     = [( getPhase(tracers(iq)%name), iq=1, nqtrue )]
    332    tracers(:)%component = type_trac
    333    DO iq = 1, nqtrue
    334       IF(addPhase('H2O',tracers(iq)%phase) == tracers(iq)%name) tracers(iq)%component = 'lmdz'
    335    END DO
    336    IF(ANY(['inca','inco'] == type_trac)) tracers(1+nqo:nqCO2+nqo)%component = 'co2i'
    337    CALL setGeneration(tracers)                                       !--- SET FIELDS %iGeneration, %gen0Name
    338    WHERE(tracers(:)%iGeneration == 2) tracers(:)%type = 'tag'        !--- DEFAULT VALUE: "tracer"
    339 
    340    !--- FINALIZE
    341    DEALLOCATE(tnom_0, tnom_transp)
    342 #ifdef INCA
    343    DEALLOCATE(hadv_inca, vadv_inca, conv_flg_inca, pbl_flg_inca, solsym_inca)
    344 #endif
    345 
    346 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    347 ELSE
    348 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    349253   IF(readTracersFiles(type_trac, fType, tracers)) CALL abort_gcm(modname, 'problem with tracers file(s)',1)
    350254   IF(fType == 0) CALL abort_gcm(modname, 'Missing tracers file: "traceur.def", "tracer.def" or "tracer_<keyword>.def file.',1)
    … …  
    389293   END IF
    390294   !---------------------------------------------------------------------------------------------------------------------------
    391 END IF
    392 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    393295
    394296   CALL getKey_init(tracers)
    … …  
    495397         modname, iad == 14 .AND. ll))                 t1%iadv = 10
    496398   END DO
    497 
    498 IF(lOldCode) THEN
    499 
    500    CALL infotrac_setHeredity                !--- SET FIELDS %iqParent, %nqChilds, %iGeneration, %gen0Name, %iqDescen, %nqDescen
    501    CALL infotrac_isoinit                    !--- SET FIELDS %type, %iso_iName, %iso_iZone, %iso_iPhase
    502    CALL getKey_init(tracers, isotopes)
    503    IF(isoSelect('H2O')) RETURN              !--- Select water isotopes ; finished if no water isotopes
    504    iH2O = ixIso                             !--- Keep track of water family index
    505 
    506    !--- Remove the isotopic tracers from the tracers list passed to phytrac
    507    nbtr    = nbtr -nqo*   ntiso             !--- ISOTOPIC TAGGING TRACERS ARE NOT PASSED TO THE PHYSICS
    508    nqtottr = nqtot-nqo*(1+ntiso)            !--- NO H2O-FAMILY    TRACER  IS      PASSED TO THE PHYSICS
    509 
    510 ELSE
    511399
    512400   niso = 0; nzone=0; nphas=nqo; ntiso = 0; isoCheck=.FALSE.
    … …  
    542430      END DO
    543431   END IF
    544 
    545 END IF
    546432
    547433   !--- Convection / boundary layer activation for all tracers
    … …  
    586472
    587473END SUBROUTINE infotrac_init
    588 
    589 
    590 
    591 SUBROUTINE infotrac_setHeredity
    592    !--- Purpose: Set fields %iqParent, %nqChilds, %iGeneration, %iqDescen, %nqDescen (old method)
    593    USE strings_mod, ONLY: strIdx
    594    INTEGER               :: iq, jq, ipere, ifils
    595    INTEGER, ALLOCATABLE  :: iqfils(:,:)
    596    CHARACTER(LEN=maxlen) :: msg1, modname='infotrac_init'
    597    INCLUDE "iniprint.h"
    598 
    599    !=== SET FIELDS %iqParent, %nqChilds
    600    ALLOCATE(iqfils(nqtot,nqtot)); iqfils(:,:) = 0
    601    tracers(:)%nqChilds = 0
    602    tracers(:)%iqParent = 0
    603 
    604    DO iq = 1, nqtot
    605       msg1 = 'Tracer nr. '//TRIM(int2str(iq))//', called "'//TRIM(tracers(iq)%name)//'" is '
    606 
    607       !--- IS IT A GENERATION 0 TRACER ? IF SO, tracers(iq)%iqParent KEEPS ITS DEFAULT VALUE (0)
    608       IF(fmsg(TRIM(msg1)//' a parent', modname, tracers(iq)%parent == tran0)) CYCLE
    609 
    610       !--- TRACERS OF GENERATION > 0 ; FIND ITS PARENT INDEX
    611       ipere = strIdx(tracers(:)%name, tracers(iq)%parent)
    612       IF(ipere == 0)  CALL abort_gcm('infotrac_init', TRIM(msg1)//' an orphan', 1)
    613       IF(iq == ipere) CALL abort_gcm('infotrac_init', TRIM(msg1)//' its own parent',1)
    614 
    615       CALL msg(TRIM(msg1)//' the child of '//TRIM(tracers(ipere)%name), modname)
    616       tracers(iq)%iqParent    = ipere
    617       tracers(ipere)%nqChilds = tracers(ipere)%nqChilds+1
    618       iqfils(tracers(ipere)%nqChilds,ipere) = iq
    619    END DO
    620 
    621    CALL msg('nqGen0   = '//int2str(COUNT(tracers(:)%parent == 'air')), modname)
    622    CALL msg('nqChilds = '//strStack(int2str(tracers(:)%nqChilds)),     modname)
    623    CALL msg('iqParent = '//strStack(int2str(tracers(:)%iqParent)),     modname)
    624    CALL msg('iqChilds = '//strStack(int2str(PACK(iqfils,MASK=.TRUE.))),modname)
    625 
    626    !=== SET FIELDS %iGeneration, %iqDescen, %nqDescen
    627    tracers(:)%iGeneration = 0
    628    tracers(:)%nqDescen = 0
    629    DO iq = 1, nqtot
    630       ifils = iq
    631       DO WHILE(tracers(ifils)%iqParent > 0)
    632          ipere = tracers(ifils)%iqParent
    633          tracers(ipere)%nqDescen = tracers(ipere)%nqDescen+1   
    634          tracers(iq)%iGeneration = tracers(iq)%iGeneration+1
    635          iqfils(tracers(ipere)%nqDescen,ipere) = iq
    636          ifils = ipere
    637       END DO
    638       msg1 = 'Tracer nr. '//TRIM(int2str(iq))//', called "'//TRIM(tracers(iq)%name)//'" is '
    639       CALL msg(TRIM(msg1)//' of generation '//TRIM(int2str(tracers(iq)%iGeneration)), modname)
    640    END DO
    641    DO iq=1,nqtot
    642       tracers(iq)%iqDescen = iqfils(1:tracers(iq)%nqDescen,iq)
    643    END DO
    644 
    645    !=== SET FIELD %gen0Name
    646    DO iq = 1, nqtot
    647       jq=iq; DO WHILE(tracers(jq)%iGeneration > 0); jq=tracers(jq)%iqParent; END DO
    648       tracers(iq)%gen0Name = tracers(jq)%name
    649    END DO
    650 
    651    CALL msg('nqDescen = '//TRIM(strStack(int2str(tracers(:)%nqDescen))), modname)
    652    CALL msg('nqDescen_tot = ' //TRIM(int2str(SUM(tracers(:)%nqDescen))), modname)
    653    CALL msg('iqDescen = '//strStack(int2str(PACK(iqfils, MASK=.TRUE.))), modname)
    654 
    655 END SUBROUTINE infotrac_setHeredity
    656 
    657 
    658 
    659 SUBROUTINE infotrac_isoinit
    660 
    661 #ifdef CPP_IOIPSL
    662    USE IOIPSL
    663 #else
    664    USE ioipsl_getincom
    665 #endif
    666    USE readTracFiles_mod, ONLY: tnom_iso => newH2OIso
    667    IMPLICIT NONE
    668    INTEGER, ALLOCATABLE  :: nb_iso(:), nb_tiso(:), nb_zone(:), ix(:), iy(:)
    669    INTEGER               :: ii, ip, iq, it, iz, ixt
    670    TYPE(isot_type), POINTER :: i
    671    TYPE(trac_type), POINTER :: t(:), t1
    672    CHARACTER(LEN=maxlen)    :: tnom_trac, modname, t0
    673    CHARACTER(LEN=maxlen), ALLOCATABLE :: str(:)
    674    LOGICAL, DIMENSION(:), ALLOCATABLE :: mask
    675    REAL, ALLOCATABLE :: tnat0(:), alpha_ideal0(:)
    676    INCLUDE "iniprint.h"
    677 
    678    modname = 'infotrac_isoinit'
    679    tnat0       = [ 1.0 ,  155.76e-6, 2005.2e-6, 0.004/100., 0.0 ]    !--- Same length as tnom_iso
    680    alpha_ideal0= [ 1.0 ,  1.01,      1.006,     1.003,      1.0 ]    !--- Same length as tnom_iso
    681    ALLOCATE(isotopes(1))                                             !--- Only water
    682    nbIso = SIZE(isotopes)
    683    t => tracers
    684    i => isotopes(1)
    685    i%parent = 'H2O'
    686 
    687    !--- Effective isotopes names list (embedded in the "keys" field)
    688    mask = t%type=='tracer' .AND. t%gen0Name==addPhase('H2O', 'g') .AND. t%iGeneration==1
    689    str = PACK(delPhase(t%name), MASK=mask)
    690    CALL strReduce(str)
    691    i%niso = SIZE(str)
    692    ALLOCATE(i%keys(i%niso))
    693    i%keys(:)%name = str
    694 
    695    !--- Check whether found isotopes are known
    696    mask = [(ALL(tnom_iso /= str(ii)), ii=1, i%niso)]
    697    IF(ANY(mask)) CALL abort_gcm(modname, 'The following isotopes are unknown: '//strStack(PACK(str, MASK=mask)), 1)
    698 
    699    !--- Full isotopes list, with isotopes tagging tracers (if any) following the previous list
    700    mask = t%type=='tag'    .AND. t%gen0Name==addPhase('H2O', 'g') .AND. t%iGeneration==2
    701    str = PACK(delPhase(t%name), MASK=mask)
    702    i%ntiso = i%niso + SIZE(str)
    703    ALLOCATE(i%trac(i%ntiso))
    704    i%trac(:) = [i%keys(:)%name, str]
    705 
    706    !--- Effective tagging zones names list
    707    i%zone = strTail(str, '_', .TRUE.)
    708    CALL strReduce(i%zone)
    709    i%nzone = SIZE(i%zone)
    710    IF(i%ntiso /= i%niso*(i%nzone+1)) CALL abort_gcm(modname, 'Error in "ntiso" calculation', 1)
    711 
    712    !--- Effective phases list
    713    i%phase = ''
    714    DO ip=1,nphases; IF(strIdx(t%name, addPhase('H2O', ip))/=0) i%phase=TRIM(i%phase)//known_phases(ip:ip); END DO
    715    i%nphas = LEN_TRIM(i%phase)
    716 
    717    !--- Indexes related to isotopes
    718    DO iq = 1, nqtot
    719       t1 => tracers(iq)
    720       t0 = t1%gen0Name
    721       IF(t1%iGeneration==0 .OR. .NOT.t1%isAdvected .OR. delPhase(t0)/='H2O') CYCLE
    722       t1%iso_iGroup = 1
    723       t1%iso_iPhase = INDEX(i%phase, getPhase(t0))
    724       t1%iso_iZone = strIdx(i%zone,  strTail(t1%name, '_'))
    725       IF(t1%iso_iZone /= 0) t1%iso_iName = strIdx(i%keys(:)%name, delPhase(t1%parent))
    726       IF(t1%iso_iZone == 0) t1%iso_iName = strIdx(i%keys(:)%name, delPhase(t1%name  ))
    727    END DO
    728 
    729    !--- Get vectors, one value each "isotope%trac" element (and in the same order)
    730    ix = strIdx(tnom_iso, i%trac)
    731    iy =   PACK(ix, MASK = ix/=0)
    732    tnat        = tnat0       (iy)
    733    alpha_ideal = alpha_ideal0(iy)
    734 
    735    !--- Tests
    736    nb_iso  = [(COUNT(t%iso_iPhase == ip .AND. t%iGeneration == 1), ip=1, i%nphas)]
    737    nb_tiso = [(COUNT(t%iso_iPhase == ip .AND. t%iGeneration == 2), ip=1, i%nphas)]
    738    nb_zone = [(COUNT(t%iso_iZone  == iz), iz=1, i%nzone)]
    739    IF(ANY(nb_iso (:) /= nb_iso (1))) CALL abort_gcm(modname, 'Phases must have same number of isotopes', 1)
    740    IF(ANY(nb_tiso(:) /= nb_tiso(1))) CALL abort_gcm(modname, 'Phases must have same number of tagging tracers', 1)
    741    IF(i%nzone >= 1) THEN
    742       IF(ANY(nb_zone(:) /= nb_zone(1))) CALL abort_gcm(modname, 'Isotopes must have the same number of tagging tracers', 1)
    743    END IF
    744 
    745    !--- Isotopic checking routines activation flag
    746    i%check = .FALSE.; IF(i%niso > 0) CALL getin('ok_iso_verif', i%check)
    747 
    748    !--- Table: index in "qx(:)" of an isotope, knowing its indices "it","ip" in "isotope%iName,%iPhase"
    749    i%iqIsoPha = RESHAPE([((strIdx(t%name, TRIM(addPhase(i%trac(it),ip,i%phase))),it=1,i%ntiso),ip=1,i%nphas)],[i%ntiso,i%nphas])
    750 
    751    !--- Table: index in "isotope%tracs(:)%name" of an isotopic tagging tracer, knowing its indices "iz","ip" in "isotope%iZone,%iName"
    752    i%itZonIso = RESHAPE([((strIdx(i%trac,TRIM(i%trac(it))//'_'//TRIM(i%zone(iz))),iz=1,i%nzone),it=1,i%niso )],[i%nzone,i%niso])
    753 
    754    DO it=1,i%ntiso; CALL msg('iqIsoPha('//TRIM(int2str(it))//',:) = '//strStack(int2str(i%iqIsoPha(it,:))), modname); END DO
    755    DO iz=1,i%nzone; CALL msg('itZonIso('//TRIM(int2str(iz))//',:) = '//strStack(int2str(i%itZonIso(iz,:))), modname); END DO
    756 
    757    !--- Finalize :
    758    DEALLOCATE(nb_iso)
    759 
    760 END SUBROUTINE infotrac_isoinit
    761474
    762475

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

    r4155 r4158  
    1111  !--- Contains all isotopic variables + their initialization
    1212  !--- Isotopes-specific routines are in isotopes_routines_mod to avoid circular dependencies with isotopes_verif_mod.
    13 
    14   LOGICAL, PARAMETER :: lOldCode=.FALSE.
    1513
    1614   !--- Isotopes indices (in [1,niso] ; non-existing => 0 index)
    … …  
    165163   CALL msg('219: entree', modname)
    166164
    167    !--- Memory allocations
    168    IF(lOldCode) THEN
    169    ALLOCATE(talph1(niso), tkcin0(niso),  talps1(niso),  tnat(niso))
    170    ALLOCATE(talph2(niso), tkcin1(niso),  talps2(niso),  toce(niso))
    171    ALLOCATE(talph3(niso), tkcin2(niso), tdifrel(niso), tcorr(niso))
    172    ALLOCATE(alpha_liq_sol(niso),   Rdefault(niso),   Rmethox(niso))
    173    END IF
    174 
    175 
    176165   !--------------------------------------------------------------
    177166   ! General:
    … …  
    231220      CALL get_in('lat_max_albedo', lat_max_albedo,  100.)
    232221   END IF
    233    IF(lOldCode) &
    234    CALL get_in('deltaO18_oce',        deltaO18_oce,   0.0)
    235222   deltaO18_oce=0.0
    236223   CALL get_in('deltaP_BL',           deltaP_BL,     10.0)
    … …  
    306293   ! Parameters that depend on the nature of water isotopes:
    307294   !--------------------------------------------------------------
    308 
    309    !===========================================================================================================================
    310    IF(lOldCode) THEN
    311    !===========================================================================================================================
    312 
    313    ! Local constants
    314    fac_enrichoce18 = 0.0005            ! Then: tcorO18 = 1 + fac_enrichoce18
    315                                        !       tcorD   = 1 + fac_enrichoce18*8
    316                                        !       tcorO17 = 1 + fac_enrichoce18*0.528
    317    alpha_liq_sol_O18 = 1.00291         ! From Lehmann & Siegenthaler, 1991,
    318                                        ! Journal of Glaciology, vol 37, p 23
    319    talph1_O18 = 1137.  ; talph2_O18 = -0.4156   ; talph3_O18 = -2.0667E-3
    320    talps1_O18 = 11.839 ; talps2_O18 = -0.028244
    321    tkcin0_O18 = 0.006  ; tkcin1_O18 =  0.000285 ; tkcin2_O18 =  0.00082
    322    tdifrel_O18 = 1./0.9723
    323 
    324    ! ln(alphaeq) ratio between O18 and O17
    325    fac_coeff_eq17_liq = 0.529          ! From Amaelle
    326   !fac_coeff_eq17_ice = 0.528          ! slope MWL
    327    fac_coeff_eq17_ice = 0.529
    328 
    329    !--- Kinetic factor for surface evaporation:
    330    ! (cf: kcin = tkcin0                  if |V|<tv0cin
    331    !      kcin = tkcin1*|Vsurf| + tkcin2 if |V|>tv0cin )
    332    ! (Rq: formula discontinuous for |V|=tv0cin... )       
    333 
    334    DO ixt = 1, niso
    335       sxt=int2str(ixt)
    336       CALL msg('80: ixt='//TRIM(int2str(ixt)),modname)
    337 
    338       Rdefault(ixt) = 0.0
    339       IF(ixt == iso_HTO) THEN          !=== H[3]HO
    340          tdifrel(ixt) = 1./0.968
    341          tkcin0(ixt) = 0.01056
    342          tkcin1(ixt) = 0.0005016
    343          tkcin2(ixt) = 0.0014432
    344          tnat  (ixt) = 0.
    345          toce  (ixt) = 4.0E-19         ! Ratio T/H = 0.2 TU, Dreisigacker and Roether 1978
    346         !toce  (ixt) = 2.2222E-8       ! Corrected by Alex Cauquoin
    347         !toce  (ixt) = 1.0E-18         ! Ratio 3H/1H ocean
    348          tcorr (ixt) = 1.
    349          talph1(ixt) = 46480. ; talph2(ixt) = -103.87 ; talph3(ixt) = 0.
    350          talps1(ixt) = 46480. ; talps2(ixt) = -103.87
    351          alpha_liq_sol(ixt) = 1.
    352          Rmethox(ixt) = 0.0
    353       ELSE IF(ixt == iso_O17) THEN     !=== H2[17]O
    354          tdifrel(ixt)=1./0.98555       ! Used in 1D and in LdG's model
    355         !tdifrel(ixt)=1./0.985452      ! From Amaelle
    356         !fac_kcin=0.5145               ! From Amaelle
    357          fac_kcin= (tdifrel(ixt)-1.0)/(tdifrel_O18-1.0)
    358          tkcin0(ixt) = tkcin0_O18*fac_kcin
    359          tkcin1(ixt) = tkcin1_O18*fac_kcin
    360          tkcin2(ixt) = tkcin2_O18*fac_kcin
    361          tnat  (ixt) = 0.004/100.      ! O17 = 0.004% of oxygen
    362          pente_MWL=0.528
    363          toce  (ixt) = tnat(ixt)*(1.0+deltaO18_oce/1000.0)**pente_MWL
    364          tcorr (ixt) = 1.0+fac_enrichoce18*pente_MWL  ! From Amaelle           
    365          talph1(ixt) = talph1_O18 ; talph2(ixt) = talph2_O18 ; talph3(ixt) = talph3_O18
    366          talps1(ixt) = talps1_O18 ; talps2(ixt) = talps2_O18     
    367          alpha_liq_sol(ixt) = (alpha_liq_sol_O18)**fac_coeff_eq17_liq
    368          IF(Rdefault_smow) Rdefault(ixt) = tnat(ixt)*(-3.15/1000.0+1.0)
    369          Rmethox(ixt) = (230./1000.+1.)*tnat(ixt)     ! Zahn et al 2006
    370       ELSE IF(ixt == iso_O18) THEN     !=== H2[18]O
    371          tdifrel(ixt) = tdifrel_O18
    372          tkcin0(ixt) = tkcin0_O18
    373          tkcin1(ixt) = tkcin1_O18
    374          tkcin2(ixt) = tkcin2_O18
    375          tnat  (ixt) = 2005.2E-6
    376          toce  (ixt) = tnat(ixt)*(1.0+deltaO18_oce/1000.0)
    377          tcorr (ixt) = 1.0+fac_enrichoce18
    378          talph1(ixt) = talph1_O18 ; talph2(ixt) = talph2_O18 ; talph3(ixt) = talph3_O18
    379          talps1(ixt) = talps1_O18 ; talps2(ixt) = talps2_O18
    380          alpha_liq_sol(ixt) = alpha_liq_sol_O18   
    381          IF(Rdefault_smow) Rdefault(ixt) = tnat(ixt)*(-6.0/1000.0+1.0)
    382          Rmethox(ixt) = (130./1000.+1.)*tnat(ixt) ! Zahn et al 2006   
    383       ELSE IF(ixt == iso_HDO) THEN     !=== H[2]HO
    384          tdifrel(ixt) = 1./0.9755
    385         !fac_kcin=0.88
    386          fac_kcin= (tdifrel(ixt)-1.0)/(tdifrel_O18-1.0)
    387          tkcin0(ixt) = tkcin0_O18*fac_kcin
    388          tkcin1(ixt) = tkcin1_O18*fac_kcin
    389          tkcin2(ixt) = tkcin2_O18*fac_kcin
    390          tnat  (ixt) = 155.76E-6
    391          pente_MWL = 8.0
    392          toce  (ixt) = tnat(ixt)*(1.0+pente_MWL*deltaO18_oce/1000.0)
    393          tcorr (ixt) = 1.0+fac_enrichoce18*pente_MWL         
    394          talph1(ixt) = 24844. ; talph2(ixt) = -76.248 ; talph3(ixt) = 52.612E-3
    395          talps1(ixt) = 16288. ; talps2(ixt) = -0.0934
    396         !alpha_liq_sol(ixt)=1.0192 ZX  ! From Weston, Ralph, 1955
    397          alpha_liq_sol(ixt)=1.0212     ! From Lehmann & Siegenthaler, 1991,
    398                                        ! Journal of Glaciology, vol 37, p 23
    399          IF(Rdefault_smow) Rdefault(ixt) = tnat(ixt)*((-6.0*pente_MWL+10.0)/1000.0+1.0)
    400          Rmethox(ixt) = tnat(ixt)*(-25.0/1000.+1.) ! Zahn et al 2006
    401       ELSE IF(ixt  == iso_eau) THEN    !=== H2O[16]
    402          tkcin0(ixt) = 0.0
    403          tkcin1(ixt) = 0.0
    404          tkcin2(ixt) = 0.0
    405          tnat  (ixt) = 1.
    406          toce  (ixt)=tnat(ixt)
    407          tcorr (ixt) = 1.0
    408          tdifrel(ixt) = 1.
    409          talph1(ixt) = 0. ; talph2(ixt) = 0. ; talph3(ixt) = 0.
    410          talps1(ixt) = 0. ; talps2(ixt) = 0.
    411          alpha_liq_sol(ixt)=1.
    412          IF(Rdefault_smow) Rdefault(ixt) = tnat(ixt)*1.0
    413          Rmethox(ixt) = 1.0
    414       END IF
    415    END DO
    416    !===========================================================================================================================
    417    ELSE
    418    !===========================================================================================================================
    419 
    420295   IF(getKey('tnat',    tnat,    isoName)) CALL abort_physic(modname, 'can''t get tnat',    1)
    421296   IF(getKey('toce',    toce,    isoName)) CALL abort_physic(modname, 'can''t get toce',    1)
    … …  
    430305   IF(getKey('tkcin2',  tkcin2,  isoName)) CALL abort_physic(modname, 'can''t get tkcin2',  1)
    431306   IF(getKey('tdifrel', tdifrel, isoName)) CALL abort_physic(modname, 'can''t get tdifrel', 1)
    432    DO ixt = 1, niso
    433       IF     (ixt == iso_HTO) THEN; tdifrel(ixt) = 1./0.968
    434       ELSE IF(ixt == iso_HDO) THEN; tdifrel(ixt) = 1./0.9755
    435       ELSE IF(ixt == iso_O17) THEN; tdifrel(ixt) = 1./0.98555
    436       ELSE IF(ixt == iso_O18) THEN; tdifrel(ixt) = 1./0.9723
    437       ELSE IF(ixt == iso_eau) THEN; tdifrel(ixt) = 1.; END IF
    438    END DO
    439307   IF(getKey('alpha_liq_sol', alpha_liq_sol, isoName)) CALL abort_physic(modname, 'can''t get alpha_liq_sol',  1)
    440308   IF(getKey('Rdefault',Rdefault,isoName)) CALL abort_physic(modname, 'can''t get Rdefault',1)
    441309   IF(getKey('Rmethox', Rmethox, isoName)) CALL abort_physic(modname, 'can''t get Rmethox', 1)
    442310   IF(.NOT.Rdefault_smow) Rdefault(:) = 0.0
    443 
    444    !===========================================================================================================================
    445    END IF
    446    !===========================================================================================================================
    447311
    448312   !--- Sensitivity test: no kinetic effect in sfc evaporation
Note: See TracChangeset for help on using the changeset viewer.