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