source: LMDZ6/branches/LMDZ-tracers/libf/dyn3d/check_isotopes.F90 @ 3852

Last change on this file since 3852 was 3852, checked in by dcugnet, 3 years ago

Extension of the tracers management.

The tracers files can be:

1) "traceur.def": old format, with:

  • the number of tracers on the first line
  • one line for each tracer: <tracer name> <hadv> <vadv> [<parent name>]

2) "tracer.def": new format with one section each model component.
3) "tracer_<name>.def": new format with a single section.

The formats 2 and 3 reading is driven by the "type_trac" key, which can be a

coma-separated list of components.

  • Format 2: read the sections from the "tracer.def" file.
  • format 3: read one section each "tracer_<section name>.def" file.
  • the first line of a section is "&<section name>
  • the other lines start with a tracer name followed by <key>=<val> pairs.
  • the "default" tracer name is reserved ; the other tracers of the section inherit its <key>=<val>, except for the keys that are redefined locally.

This format helps keeping the tracers files compact, thanks to the "default"
special tracer and the three levels of factorization:

  • on the tracers names: a tracer name can be a coma-separated list of tracers => all the tracers of the list have the same <key>=<val> properties
  • on the parents names: the value of the "parent" property can be a coma-separated list of tracers => only possible for geographic tagging tracers
  • on the phases: the property "phases" is [g](l][s] (gas/liquid/solid)

Read information is stored in the vector "tracers(:)", of derived type "tra".

"isotopes_params.def" is a similar file, with one section each isotopes family.
It contains a database of isotopes properties ; if there are second generation
tracers (isotopes), the corresponding sections are read.

Read information is stored in the vector "isotopes(:)", of derived type "iso".

The "getKey" function helps to get the values of the parameters stored in
"tracers" or "isotopes".
File size: 5.7 KB
Line 
1SUBROUTINE check_isotopes_seq(q, ip1jmp1, err_msg)
2  USE strings_mod, ONLY: strIdx, msg, modname, prt_level
3  USE infotrac,    ONLY: isotope, isoSelect, iH2O, isoCheck, isoName, nqtot, niso, nitr, nzon, npha, iTraPha, iZonIso, tnat
4  IMPLICIT NONE
5  include "dimensions.h"
6  REAL,             INTENT(INOUT) :: q(ip1jmp1,llm,nqtot)
7  INTEGER,          INTENT(IN)    :: ip1jmp1
8  CHARACTER(LEN=*), INTENT(IN)    :: err_msg     !--- Error message to display
9  CHARACTER(LEN=256) :: msg1
10  INTEGER :: ixt, ipha, k, i, iq, iiso, izon, ieau, iqeau
11  INTEGER, ALLOCATABLE :: ix(:)
12  REAL    :: xtractot, xiiso, deltaD, q1, q2
13  REAL, PARAMETER :: borne     = 1e19,  &
14                     errmax    = 1e-8,  &        !--- Max. absolute error
15                     errmaxrel = 1e-8,  &        !--- Max. relative error
16                     qmin      = 1e-11, &
17                     deltaDmax = 200.0, &
18                     deltaDmin =-999.9, &
19                     ridicule  = 1e-12
20  INTEGER, SAVE :: ixH2O, ixHDO, ixO18
21  LOGICAL, SAVE :: first=.TRUE.
22
23  modname = 'check_isotopes'
24  IF(first) THEN
25    IF(isoSelect('H2O')) RETURN
26    ixH2O = strIdx(isoName,'H2[16]O')
27    ixHDO = strIdx(isoName,'H[2]HO')
28    ixO18 = strIdx(isoName,'H2[18]O')
29    first = .FALSE.
30  ELSE
31    IF(isoSelect(iH2O)) RETURN
32  END IF
33  IF(.NOT.isoCheck .OR. niso == 0) RETURN        !--- No need to check or no isotopes => finished
34
35  !--- CHECK FOR NaNs (FOR ALL THE ISOTOPES, INCLUDING GEOGRAPHIC TAGGING TRACERS)
36  DO ixt = 1, nitr
37    DO ipha = 1, npha
38      iq = iTraPha(ixt,ipha)
39      DO k = 1, llm
40        DO i = 1, ip1jmp1
41          IF(ABS(q(i,k,iq))<=borne) CYCLE
42          WRITE(msg1,'(s,"(",i0,",",i0,",",i0,") = ",ES12.4)')isoName(ixt),i,k,iq,q(i,k,iq); CALL msg(msg1)
43          CALL abort_gcm(modname, 'Error in iso_verif_noNaN: '//TRIM(err_msg), 1)
44          STOP
45        END DO
46      END DO
47    END DO
48  END DO
49
50  !--- CHECK CONSERVATION (MAIN ISOTOPE AND PARENT CONCENTRATIONS MUST BE EQUAL)
51  ixt = ixH2O
52  IF(ixt /= 0) THEN
53    DO ipha = 1, npha
54      iq = iTraPha(ixt,ipha)
55      DO k = 1, llm
56        DO i = 1, ip1jmp1
57          q1 = q(i,k,ipha); q2 = q(i,k,iq)
58          IF(ABS(q1-q2) <= errmax .OR. ABS(q1-q2)/MAX(MAX(ABS(q1),ABS(q2)),1e-18) < errmaxrel) CYCLE
59          WRITE(msg1,'("ixt = ",i0)')ixt;                                      CALL msg(msg1)
60          WRITE(msg1,'("q(",i0,",",i0,",iq=",i0,") = ",ES12.4)')i, k, iq, q2;  CALL msg(msg1)
61          WRITE(msg1,'("q(",i0,",",i0,",ipha=",i0,") = ",ES12.4)')i,k,ipha,q1; CALL msg(msg1)
62          CALL abort_gcm(modname, 'Error in iso_verif_egalite: '//TRIM(err_msg), 1)
63          q(i,k,iq) = q(i,k,ipha)                !--- Bidouille pour convergence
64        END DO
65      END DO
66    END DO
67  END IF
68
69  !--- CHECK DELTA ANOMALIES
70  ix = [ixHDO, ixO18]
71  DO iiso = 1, SIZE(ix)
72    ixt = ix(iiso)
73    IF(ixt  == 0) CYCLE
74    DO ipha = 1, npha
75      iq = iTraPha(ixt,ipha)
76      DO k = 1, llm
77        DO i = 1, ip1jmp1
78          q1 = q(i,k,ipha); q2 = q(i,k,iq)
79          IF(q2 <= qmin) CYCLE
80          deltaD = (q2/q1/tnat(ixt)-1)*1000
81          IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
82          WRITE(msg1,'("ixt = ",i0)')ixt;                                     CALL msg(msg1)
83          WRITE(msg1,'("q(",i0,",",i0,",iq=",i0,") = ",ES12.4)')i, k, iq, q2; CALL msg(msg1)
84          WRITE(msg1,'("q=",ES12.4)')q(i,k,:);                                CALL msg(msg1)
85          WRITE(msg1,'("deltaD=",ES12.4)')deltaD;                             CALL msg(msg1)
86          CALL abort_gcm(modname, 'Error in iso_verif_aberrant: '//TRIM(err_msg), 1)
87        END DO
88      END DO
89    END DO
90  END DO
91
92  !--- CHECK FOR TAGGING TRACERS DELTAD ANOMALIES
93  IF(nitr == 0) RETURN
94  IF(ixH2O /= 0 .AND. ixHDO /= 0) THEN
95    DO izon = 1, nzon
96      ixt  = iZonIso(izon, ixHDO)
97      ieau = iZonIso(izon, ixH2O)
98      DO ipha = 1, npha
99        iq    = iTraPha(ixt,  ipha)
100        iqeau = iTraPha(ieau, ipha)
101        DO k = 1, llm
102          DO i = 1, ip1jmp1
103            IF(q(i,k,iq)<=qmin) CYCLE
104            deltaD = (q(i,k,iq)/q(i,k,iqeau)/tnat(ixHDO)-1)*1000
105            IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
106            WRITE(msg1,'("izon, ipha =",2i0)')izon, ipha;                              CALL msg(msg1)
107            WRITE(msg1,'( "ixt, ieau =",2i0)') ixt, ieau;                              CALL msg(msg1)
108            WRITE(msg1,'("q(",i0,",",i0,",iq=",i0,") = ",ES12.4)')i, k, iq, q(i,k,iq); CALL msg(msg1)
109            WRITE(msg1,'("deltaD=",ES12.4)')deltaD;                                    CALL msg(msg1)
110            CALL abort_gcm(modname, 'Error in iso_verif_aberrant trac: '//TRIM(err_msg), 1)
111          END DO
112        END DO
113      END DO
114    END DO
115  END IF
116
117  !--- CHECK FOR TAGGING TRACERS CONSERVATION (PARENT AND TAGGING TRACERS SUM OVER ALL REGIONS MUST BE EQUAL)
118  DO iiso = 1, niso
119    DO ipha = 1, npha
120      iq = iTraPha(iiso, ipha)
121      DO k = 1, llm
122        DO i = 1, ip1jmp1
123          xiiso = q(i,k,iq)
124          xtractot = SUM(q(i, k, iTraPha(iZonIso(1:nzon,iiso), ipha)))
125          IF(ABS(xtractot-xiiso) > errmax .AND. ABS(xtractot-xiiso)/MAX(MAX(ABS(xtractot),ABS(xiiso)),1e-18) > errmaxrel) THEN
126            CALL msg('Error in iso_verif_aberrant trac: '//TRIM(err_msg))
127            WRITE(msg1,'("iiso, ipha =",2i0)')iiso, ipha;              CALL msg(msg1)
128            WRITE(msg1,'("i, k =",2i0)')i, k;                          CALL msg(msg1)
129            WRITE(msg1,'("q(",i0,",",i0,":) = ",ES12.4)')i,k,q(i,k,:); CALL msg(msg1)
130            STOP
131          END IF
132          IF(ABS(xtractot) <= ridicule) CYCLE
133          DO izon = 1, nzon
134            ixt = iZonIso(izon, iiso)
135            q(i,k,iq) = q(i,k,iq) / xtractot * xiiso
136          END DO
137        END DO
138      END DO
139    END DO
140  END DO
141
142END SUBROUTINE check_isotopes_seq
143
Note: See TracBrowser for help on using the repository browser.