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.

File:
1 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
Note: See TracChangeset for help on using the changeset viewer.