Ignore:
Timestamp:
May 11, 2021, 2:10:34 PM (3 years ago)
Author:
dcugnet
Message:
  • Bugs corrections:
    • sequential gcm fixed
    • parallel gcm compilation fixed ; to be tested
  • Some generic operations moved from infotrac to readTracFile
  • Fixed algebrical reduction routine, used in the isotopes parameters file.
  • Additional component "comp" in the tracers descriptor derived type "tra",

specifying the model component name(s) (cf. tracers sections) it belongs.

  • isotopes class selection tool fixed.
File:
1 edited

Legend:

Unmodified
Added
Removed
  • LMDZ6/branches/LMDZ-tracers/libf/dyn3d_common/infotrac.F90

    r3852 r3891  
    11MODULE infotrac
    22
    3   USE       strings_mod, ONLY: msg, find, strIdx,  strFind,  strHead, dispTable, cat, get_in,  &
    4                               fmsg, test, int2str, strParse, strTail, strReduce, strStack, modname, testFile
     3  USE       strings_mod, ONLY: msg, find, strIdx,  strFind,  strHead, dispTable, testFile, cat, get_in,  &
     4                              fmsg, test, int2str, strParse, strTail, strReduce, strStack, modname
    55  USE readTracFiles_mod, ONLY: readTracersFiles, getKey_init, nphases, delPhase, aliasTracer, &
    6                         tran0, readIsotopesFile, getKey, known_phases, addPhase, indexUpdate
     6                        tran0, readIsotopesFile, getKey, known_phases, addPhase, indexUpdate, initIsotopes
    77  USE trac_types_mod,    ONLY: tra, iso, kys
    88
     
    2323  PUBLIC :: iso, isotopes, nbIso                           !--- Derived type, full isotopes families database + nb of families
    2424  PUBLIC :: isoSelect , ixIso                              !--- Isotopes family selection tool + selected family index
     25  PUBLIC :: qprntmin, massqmin, ratiomin                   !--- Min. values
    2526  !=== FOR ISOTOPES: Specific to H2O isotopes
    2627  PUBLIC :: iH2O, tnat, alpha_ideal                        !--- H2O isotopes index, natural abundance, fractionning coeff.
     
    3132  PUBLIC :: iZonIso, iTraPha                               !--- 2D index tables to get "iq" index
    3233  PUBLIC :: isoCheck                                       !--- Run isotopes checking routines
    33 
    3434  !=== FOR BOTH TRACERS AND ISOTOPES
    3535  PUBLIC :: getKey                                         !--- Get a key from "tracers" or "isotope"
     
    7575!  | type       | Type (so far: tracer or tag)                    | /           | tracer,tag             |
    7676!  | phas       | Phases list ("g"as / "l"iquid / "s"olid)        | /           | [g][l][s]              |
     77!  | comp       | Name(s) of the merged/cumulated section(s)      | /           | coma-separated names   |
    7778!  | iadv       | Advection scheme number                         | iadv        | 1-20,30 exc. 3-9,15,19 |
    7879!  | igen       | Generation (>=1)                                | /           |                        |
     
    103104
    104105
     106  REAL, PARAMETER :: qprntmin=1.E-12, massqmin=1.E-12, ratiomin=1.E-12
    105107
    106108  !=== DIMENSIONS OF THE TRACERS TABLES AND OTHER SCALAR VARIABLES
     
    113115
    114116  !=== DERIVED TYPES EMBEDDING MOST INFORMATIONS ABOUT TRACERS AND ISOTOPES FAMILIES
    115   TYPE(tra), TARGET,  SAVE, ALLOCATABLE ::  tracers(:)      !=== TRACERS DESCRIPTORS VECTOR
    116   TYPE(iso), TARGET,  SAVE, ALLOCATABLE :: isotopes(:)      !=== ISOTOPES PARAMETERS VECTOR
     117  TYPE(tra), TARGET,  SAVE, ALLOCATABLE ::  tracers(:)     !=== TRACERS DESCRIPTORS VECTOR
     118  TYPE(iso), TARGET,  SAVE, ALLOCATABLE :: isotopes(:)     !=== ISOTOPES PARAMETERS VECTOR
    117119!$OMP THREADPRIVATE(tracers, isotopes)
    118120
    119121  !=== ALIASES FOR CURRENTLY SELECTED ISOTOPES FAMILY OF VARIABLES EMBEDDED IN THE VECTOR "isotopes"
    120   TYPE(iso),          SAVE, POINTER :: isotope             !--- CURRENTLY SELECTED ISOTOPES FAMILY DESCRIPTOR
    121   INTEGER,            SAVE          :: ixIso, iH2O         !--- Index of the selected isotopes family and H2O family
    122   LOGICAL,            SAVE, POINTER :: isoCheck            !--- Flag to trigger the checking routines
    123   TYPE(kys),          SAVE, POINTER :: isoKeys(:)          !--- ONE SET OF KEYS FOR EACH ISOTOPE (LISTED IN isoName)
    124   CHARACTER(LEN=256), SAVE, POINTER :: isoName(:),       & !--- ISOTOPES NAMES FOR THE CURRENTLY SELECTED FAMILY
    125                                        isoZone(:),       & !--- TAGGING ZONES  FOR THE CURRENTLY SELECTED FAMILY
    126                                        isoPhas             !--- USED PHASES    FOR THE CURRENTLY SELECTED FAMILY
    127   INTEGER,            SAVE          :: niso, nzon, npha, & !--- NUMBER OF ISOTOPES, TAGGING ZONES AND PHASES
    128                                        nitr                !--- NUMBER OF ISOTOPES + ISOTOPIC TAGGING TRACERS
    129   INTEGER,            SAVE, POINTER :: iZonIso(:,:)        !--- INDEX IN "isoTrac" AS f(tagging zone, isotope)
    130   INTEGER,            SAVE, POINTER :: iTraPha(:,:)        !=== INDEX IN "isoTrac" AS f(isotopic tracer, phase)
     122  TYPE(iso),          SAVE, POINTER     :: isotope         !--- CURRENTLY SELECTED ISOTOPES FAMILY DESCRIPTOR
     123  INTEGER,            SAVE              :: ixIso, iH2O     !--- Index of the selected isotopes family and H2O family
     124  LOGICAL,            SAVE              :: isoCheck        !--- Flag to trigger the checking routines
     125  TYPE(kys),          SAVE, POINTER     :: isoKeys(:)      !--- ONE SET OF KEYS FOR EACH ISOTOPE (LISTED IN isoName)
     126  CHARACTER(LEN=256), SAVE, POINTER     :: isoName(:),   & !--- ISOTOPES NAMES FOR THE CURRENTLY SELECTED FAMILY
     127                                           isoZone(:),   & !--- TAGGING ZONES  FOR THE CURRENTLY SELECTED FAMILY
     128                                           isoPhas         !--- USED PHASES    FOR THE CURRENTLY SELECTED FAMILY
     129  INTEGER,            SAVE              :: niso, nzon,  & !--- NUMBER OF ISOTOPES, TAGGING ZONES AND PHASES
     130                                           npha, nitr      !--- NUMBER OF PHASES AND ISOTOPES + ISOTOPIC TAGGING TRACERS
     131  INTEGER,            SAVE, POINTER     :: iZonIso(:,:)    !--- INDEX IN "isoTrac" AS f(tagging zone, isotope)
     132  INTEGER,            SAVE, POINTER     :: iTraPha(:,:)    !--- INDEX IN "isoTrac" AS f(isotopic tracer, phase)
    131133!$OMP THREADPRIVATE(isotope, ixIso,iH2O, isoCheck, isoKeys, isoName,isoZone,isoPhas, niso,nzon,npha,nitr, iZonIso,iTraPha)
    132134
     
    137139                                            pbl_flg(:),  & !--- Boundary layer activation ; needed for INCA        (nbtr)
    138140                                         itr_indice(:),  & !--- Indexes of the tracers passed to phytrac        (nqtottr)
    139                                               niadv(:)
     141                                              niadv(:)     !--- Indexes of true tracers  (<=nqtot, such that iadv(idx)>0)
    140142  CHARACTER(LEN=8),   SAVE, ALLOCATABLE ::   solsym(:)     !--- Names from INCA                                    (nbtr)
    141 !OMP THREADPRIVATE(tnat, alpha_ideal, conv_flg, pbl_flg, itr_indice, solsym)
     143!OMP THREADPRIVATE(tnat, alpha_ideal, conv_flg, pbl_flg, itr_indice, niadv, solsym)
    142144
    143145#ifdef CPP_StratAer
     
    153155#ifdef REPROBUS
    154156  USE chem_rep,    ONLY: Init_chem_rep_trac
     157  IMPLICIT NONE
    155158#endif
    156159!==============================================================================================================================
     
    178181! Local variables
    179182  INTEGER, ALLOCATABLE :: hadv(:), hadv_inca(:), &                   !--- Horizontal/vertical transport scheme number
    180                           vadv(:), vadv_inca(:)                      !--- + specific INCA versions
    181   CHARACTER(LEN=1)   :: ph                                           !--- Phase
     183                          vadv(:), vadv_inca(:)                      !---   + specific INCA versions
    182184  CHARACTER(LEN=2)   ::   suff(9)                                    !--- Suffixes for schemes of order 3 or 4 (Prather)
    183   CHARACTER(LEN=3)   :: descrq(30)                                   !--- Advection scheme description
    184   CHARACTER(LEN=4)   :: oldH2O(3)                                    !--- Old water names
    185   CHARACTER(LEN=256) :: newH2O, iname, isoPhase                      !--- New water and isotope names, phases list
     185  CHARACTER(LEN=3)   :: descrq(30)                                   !--- Advection scheme description tags
     186  CHARACTER(LEN=4)   :: oldH2O(3)                                    !--- Old water name for the three phases
     187  CHARACTER(LEN=256) :: newH2O                                       !--- New water name
    186188  CHARACTER(LEN=256) :: msg1, msg2                                   !--- Strings for messages
    187   CHARACTER(LEN=256), ALLOCATABLE, DIMENSION(:) :: &                 !--- Temporary storage
    188              isoName, isoZone, tra0, zon0, tag0, n, p, z, str
     189  CHARACTER(LEN=256), ALLOCATABLE :: str(:)                          !--- Temporary storage
    189190  INTEGER :: fType                                                   !--- Tracers description file type ; 0: none
    190191                                                                     !--- 1: "traceur.def"  2: "tracer.def"  3: "tracer_*.def"
    191192  INTEGER :: nqtrue                                                  !--- Tracers nb from tracer.def (no higher order moments)
    192   INTEGER :: iad                                                     !--- Advection scheme
    193   INTEGER :: iH2O                                                    !--- Index in "isotopes(:)" of H2O family
    194   INTEGER :: ic,ip,iq,jq, it,nt, im,nm, ix, iz, niso, nzone, ntiso   !--- Indexes and temporary variables
    195   LOGICAL, ALLOCATABLE :: lisoGen2(:), &                             !--- Mask for second generation isotopes
    196                           lisoName(:), &                             !--- Mask for water isotopes
    197                           lisoZone(:), ll(:)                         !--- Mask for water isotopes tagging tracers
     193  INTEGER :: iad                                                     !--- Advection scheme number
     194  INTEGER :: ic, ip, iq, jq, it, nt, im, nm, ix, iz                  !--- Indexes and temporary variables
    198195  LOGICAL :: lerr
    199196  TYPE(tra), ALLOCATABLE, TARGET :: ttr(:)
    200   TYPE(tra), POINTER             :: t1, t(:)
     197  TYPE(tra), POINTER             :: t1
    201198  TYPE(iso), POINTER             :: s
    202199!------------------------------------------------------------------------------------------------------------------------------
     
    204201!------------------------------------------------------------------------------------------------------------------------------
    205202  modname = 'infotrac_init'
    206   type_trac='lmdz'!'lmdz,inca'
    207203  suff          = ['x ','y ','z ','xx','xy','xz','yy','yz','zz']
    208204  descrq( 1: 2) = ['LMV','BAK']
     
    310306        IF(nqo/=2 .AND. nqo/=3) CALL abort_gcm(modname,TRIM(msg1),1)
    311307#ifdef INCA
    312         CALL Init_chem_inca_trac(nbtr)                                   !--- Get nbtr from INCA
     308        CALL Init_chem_inca_trac(nbtr)                               !--- Get nbtr from INCA
    313309#endif
    314310        ALLOCATE(hadv_inca(nbtr), vadv_inca(nbtr), conv_flg(nbtr), pbl_flg(nbtr), solsym(nbtr))
     
    317313        CALL init_transport(hadv_inca, vadv_inca, conv_flg,   pbl_flg,   solsym)
    318314#endif
    319         nqtrue = nbtr + nqo                                              !--- Total number of tracers
     315        nqtrue = nbtr + nqo                                          !--- Total number of tracers
    320316        ALLOCATE(ttr(nqtrue)); ttr(1:nqo) = tracers(1:nqo)
    321317        DO iq = nqo+1, nqtrue
     
    405401    CALL msg('The total number of tracers needed is '//TRIM(int2str(nqtot)))
    406402  END IF
     403  CALL msg('nqtot = '//TRIM(int2str(nqtot)))
     404  CALL msg('nbtr  = '//TRIM(int2str(nbtr)))
     405  CALL msg('nqo   = '//TRIM(int2str(nqo)))
    407406  ALLOCATE(ttr(nqtot))
    408407
     
    424423    t1%lnam = t1%name; IF(iad /= 0) t1%lnam=TRIM(t1%name)//descrq(iad)
    425424
    426     !--- Defining most fields of the tracer derived type
     425    !--- Define most fields of the tracers derived type
    427426    ttr(jq)%name = t1%name
    428427    ttr(jq)%nam1 = t1%nam1
     
    478477
    479478  CALL MOVE_ALLOC(FROM=ttr, TO=tracers)
    480   t => tracers
    481 
    482   !=== VARIABLES RELATED TO GENERATIONS
    483   niadv = PACK( [(iq,iq=1,nqtot)], MASK=t(:)%iadv>=0)           !--- Indexes of "true" tracers
    484 
    485   p = PACK(delPhase(t%prnt),MASK=t%type=='tracer'.AND.t%igen==2)!--- Parents of 2nd generation isotopes
    486   CALL strReduce(p, nbIso)
    487   ALLOCATE(isotopes(nbIso))
    488 
    489   IF(nbIso==0) RETURN                                           !=== NO ISOTOPES: FINISHED
    490 
    491   CALL msg('Isotopes families required: '//strStack(p))
    492 
    493   !--- ISOTOPES RELATED VARIABLES ; NULL OR EMPTY IF NO ISOTOPES
    494   isotopes(:)%prnt = p
    495   DO ip = 1, SIZE(p)                                            !--- Loop on isotopes categories
    496     s => isotopes(ip)
    497     iname = s%prnt
    498 
    499     !=== Geographic tagging tracers descending on tracer "iname": mask, names, number
    500     lisoZone = t(:)%type=='tag'    .AND. delPhase(t(:)%nam1) == iname .AND. t(:)%igen == 3
    501     s%zone = PACK(strTail(t(:)%name,'_'), MASK = lisoZone)      !--- Tagging zones names  for isotopes category "iname"
    502     CALL strReduce(s%zone)
    503     s%nzon = SIZE(s%zone)                                       !--- Tagging zones number for isotopes category "iname"
    504 
    505     !=== Isotopes childs of tracer "iname": mask, names, number (same for each phase of "iname")
    506     lisoName = t(:)%type=='tracer' .AND. delPhase(t(:)%prnt) == iname .AND. t(:)%phas == 'g'
    507     ALLOCATE(s%keys(COUNT(lisoName)))
    508     s%keys(:)%name = PACK(delPhase(t(:)%name), MASK = lisoName)    !--- Effectively found isotopes of "iname"
    509     s%niso = SIZE(s%keys)                                       !--- Number of "effectively found isotopes of "iname"
    510     s%trac = [s%keys%name, ((TRIM(s%keys(it)%name)//'_'//TRIM(s%zone(iz)), it=1, s%niso), iz=1, s%nzon)]
    511     s%nitr = SIZE(s%trac)                                       !--- " + their geographic tracers               [ntraciso]
    512 
    513     !=== Phases for tracer "iname"
    514     s%phas = ''
    515     DO ix = 1, nphases; IF(strIdx(t%name,addPhase(iname, known_phases(ix:ix))) /= 0) s%phas = TRIM(s%phas)//ph; END DO
    516     s%npha = LEN_TRIM(s%phas)                                   !--- Equal to "nqo" for water
    517 
    518     !=== Tables giving the index in a table of effectively found items for each dynamical tracer (1<=iq<=nqtot)
    519     DO iq = 1, nqtot
    520       t1 => tracers(iq)
    521       IF(t1%nam1 /= iname) CYCLE                                 !--- Only deal with tracers descending on "iname"
    522       t1%iso_igr = ip                                            !--- Index of isotopes family in list "isotopes(:)%prnt"
    523       t1%iso_num = strIdx(s%trac, delPhase(strHead(t1%name,'_')))!--- Index of current isotope       in effective isotopes list
    524       t1%iso_zon = strIdx(s%zone,          strTail(t1%name,'_') )!--- Index of current isotope zone  in effective zones    list
    525       t1%iso_pha =  INDEX(s%phas,TRIM(t1%phas))                  !--- Index of current isotope phase in effective phases   list
    526       IF(t1%igen /= 3) t1%iso_zon = 0                            !--- Skip possible generation 2 tagging tracers
    527     END DO
    528 
    529     !=== Table used to get iq (index in dyn array, size nqtot) from the isotope and phase indexes ; the full isotopes list
    530     !    (including tagging tracers) is sorted this way:  iso1, iso2, ..., iso1_zone1, iso2_zone1, ..., iso1_zoneN, iso2_zoneN
    531     s%iTraPha = RESHAPE( [( (strIdx(t(:)%name,  addPhase(s%trac(it),s%phas(ip:ip))),     it=1, s%nitr), ip=1, s%npha)], &
    532                          [s%nitr, s%npha] )
    533 
    534     !=== Table used to get ix (index in tagging tracers isotopes list, size nitr) from the zone and isotope indexes
    535     s%iZonIso = RESHAPE( [( (strIdx(s%trac(:), TRIM(s%trac(it))//'_'//TRIM(s%zone(iz))), iz=1, s%nzon), it=1, s%niso)], &
    536                          [s%nzon, s%niso] )
    537   END DO
    538 
    539   !=== Indexes, in dynamical tracers list, of the tracers transmitted to phytrac (nqtottr non-vanishing elements)
    540   ll = delPhase(t%name)/='H2O' .AND. t%iso_num ==0              !--- Mask of tracers passed to the physics
    541   t(:)%itr = UNPACK([(iq,iq=1,COUNT(ll))], ll, [(0, iq=1, nqtot)])
    542   itr_indice = PACK(t(:)%itr, MASK = t(:)%itr/=0)               !--- Might be removed (t%itr should be enough)
    543 
    544   !=== READ PHYSICAL PARAMETERS FROM "isotopes_params.def" FILE
     479
     480  !=== Indexes of: "true" tracers, in the dynamical table of tracers transmitted to phytrac (nqtottr non-vanishing elements)
     481  niadv      = PACK([(iq,iq=1,nqtot)], MASK=tracers(:)%iadv>=0) !--- Indexes of "true" tracers
     482  itr_indice = PACK(tracers(:)%itr,    MASK=tracers(:)%itr /=0) !--- Might be removed (t%itr should be enough)
     483
     484  CALL initIsotopes(tracers, isotopes)
     485  nbIso = SIZE(isotopes); IF(nbIso==0) RETURN                   !--- No isotopes: finished.
     486
     487
     488  !=== READ PHYSICAL PARAMETERS FROM "isotopes_params.def" FILE SPECIFIC TO WATER ISOTOPES
    545489  !    DONE HERE, AND NOT ONLY IN "infotrac_phy", BECAUSE SOME PHYSICAL PARAMS ARE NEEDED FOR RESTARTS (tnat AND alpha_ideal)
    546   IF(readIsotopesFile('isotopes_params.def',isotopes)) CALL abort_gcm(modname,'Problem when reading isotopes parameters',1)
    547 print*,'coincoin'
    548 
    549   !=== Specific to water
    550490  CALL getKey_init(tracers, isotopes)
    551491  IF(isoSelect('H2O')) RETURN                                   !--- Select water isotopes ; finished if no water isotopes.
    552492  iH2O = ixIso                                                  !--- Keep track of water family index
    553   lerr = getKey('tnat' ,tnat,        isoName)
    554   lerr = getKey('alpha',alpha_ideal, isoName)
     493  IF(getKey('tnat' ,tnat,        isoName(1:niso))) CALL abort_gcm(modname,'can''t read "tnat"',1)
     494  IF(getKey('alpha',alpha_ideal, isoName(1:niso))) CALL abort_gcm(modname,'can''t read "alpha_ideal"',1)
    555495  CALL msg('end')
    556496
     
    560500!==============================================================================================================================
    561501!=== THE ROUTINE isoSelect IS USED TO SWITCH FROM AN ISOTOPE FAMILY TO ANOTHER: ISOTOPES DEPENDENT PARAMETERS ARE UPDATED
    562 !     Singe generic "isoSelect" routine, using the predefined parent index (fast version) or its name (first time).
     502!     Singe generic "isoSelect" routine, using the predefined parent index (fast version) or its name (first call).
    563503!==============================================================================================================================
    564 LOGICAL FUNCTION isoSelectByName(iName) RESULT(lerr)
    565   CHARACTER(LEN=*), INTENT(IN)  :: iName
     504LOGICAL FUNCTION isoSelectByName(iName, lVerbose) RESULT(lerr)
     505  IMPLICIT NONE
     506  CHARACTER(LEN=*),  INTENT(IN)  :: iName
     507  LOGICAL, OPTIONAL, INTENT(IN) :: lVerbose
    566508  INTEGER :: iIso
     509  LOGICAL :: lV
     510  lV = .FALSE.; IF(PRESENT(lVerbose)) lV = lVerbose
    567511  iIso = strIdx(isotopes(:)%prnt, iName)
    568   IF(test(fmsg(iIso == 0,'no isotope family named "'//TRIM(iName)//'"'),lerr)) RETURN
    569   IF(isoSelectByIndex(iIso)) RETURN
     512  lerr = iIso == 0
     513  CALL msg(lerr .AND. lV, 'no isotope family named "'//TRIM(iName)//'"')
     514  IF(lerr) RETURN
     515  lerr = isoSelectByIndex(iIso)
    570516END FUNCTION isoSelectByName
    571517!==============================================================================================================================
    572 LOGICAL FUNCTION isoSelectByIndex(iIso) RESULT(lerr)
    573   INTEGER, INTENT(IN) :: iIso
     518LOGICAL FUNCTION isoSelectByIndex(iIso, lVerbose) RESULT(lerr)
     519  IMPLICIT NONE
     520  INTEGER,           INTENT(IN) :: iIso
     521  LOGICAL, OPTIONAL, INTENT(IN) :: lVerbose
     522  LOGICAL :: lv
     523  lv = .FALSE.; IF(PRESENT(lVerbose)) lv = lVerbose
    574524  lerr = .FALSE.
    575525  IF(iIso == ixIso) RETURN                                      !--- Nothing to do if the index is already OK
    576   IF(test(fmsg(iIso<=0 .OR. iIso>=nbIso,'Inconsistent isotopes family index '//TRIM(int2str(iIso))),lerr)) RETURN
     526  lerr = iIso<=0 .OR. iIso>nbIso
     527  CALL msg(lerr .AND. lV, 'Inconsistent isotopes family index '//TRIM(int2str(iIso))//': should be > 0 and <= ' &
     528                                                               //TRIM(int2str(nbIso))//'"')
     529  IF(lerr) RETURN
    577530  ixIso = iIso                                                  !--- Update currently selected family index
    578531  isotope => isotopes(ixIso)                                    !--- Select corresponding component
    579   !--- VARIOUS ALIASES
    580   isoKeys => isotope%keys; niso = isotope%niso
    581   isoName => isotope%trac; nitr = isotope%nitr; isoCheck => isotope%check
    582   isoZone => isotope%zone; nzon = isotope%nzon; iZonIso  => isotope%iZonIso
    583   isoPhas => isotope%phas; npha = isotope%npha; iTraPha  => isotope%iTraPha
     532  isoKeys => isotope%keys;    niso     = isotope%niso
     533  isoName => isotope%trac;    nitr     = isotope%nitr
     534  isoZone => isotope%zone;    nzon     = isotope%nzon
     535  isoPhas => isotope%phas;    npha     = isotope%npha
     536  iZonIso => isotope%iZonIso; isoCheck = isotope%check
     537  iTraPha => isotope%iTraPha
    584538END FUNCTION isoSelectByIndex
    585539!==============================================================================================================================
Note: See TracChangeset for help on using the changeset viewer.