source: LMDZ6/trunk/libf/dyn3dmem/check_isotopes_loc.F90 @ 5183

Last change on this file since 5183 was 5183, checked in by dcugnet, 8 days ago
  • Remove INCA retro-compatibility with "traceur.def" (containing only water tracers but getting chemical species from an internal INCA routine).
  • The "trac_type" derived type internal to "readTracFiles_mod" is removed because a generic "keys_type" is enough: no explicit key ("%" operator) is needed, even %name.
  • The "trac_type" and "isot_type" derived types are now defined locally in "infotrac" and "infotrac_phy" (and more generally in each context: dynamic, lmdz dynamics, lmdz physics, etc.). The "readTracFiles_mod" module is now only used in these two routines:
    • few internal routines/variables (addPhase, delPhase, new2oldH2O, newHNO3, oldHNO3) are made available through "infotrac" and "infotrac_phy".
    • the "getKey" routine is only used in these two routines to define the explicit keys ("%" operator) of the local derived types "trac_type" and "isot_type". It could be in principle used outside this scope to get tracers parameters (read from "tracer.def") or isotopic parameters (read from "isotopes_params.def" - disabled for now).
  • The fortran parameters file "iso_params_mod.F90" is introduced so that "tnat" and "alpha_ideal" are defined in a single place but used in several. "ltnat1" is no longer hardcoded but defined with the *.def files parameter "tnat1"
  • Few minor changes:
    • use "infotrac_phy" instead of "infotrac" in calfis* because "tracers(:)%isAdvected" is defined in physics only.
    • "isotopes_mod" now ready for several isotopes classes (currently: only H2O)
    • isotopes class name (the name of the parent of the isotopes) is now %name and no longer %parent.
    • improvement of "getKey"
File size: 8.8 KB
Line 
1SUBROUTINE check_isotopes(q, ijb, ije, err_msg)
2   USE parallel_lmdz
3   USE strings_mod, ONLY: maxlen, msg, strIdx, strStack, int2str, real2str
4   USE infotrac,    ONLY: nqtot, niso, nphas, isotope, isoCheck, iqIsoPha, isoSelect, &
5                          ntiso, iH2O, nzone, tracers, isoName,  itZonIso
6   USE iso_params_mod, ONLY: tnat_H216O, tnat_H217O, tnat_H218O, tnat_HDO, tnat_HTO
7#ifdef CPP_IOIPSL
8   USE ioipsl,          ONLY: getin
9#else
10   USE ioipsl_getincom, ONLY: getin
11#endif
12   IMPLICIT NONE
13   include "dimensions.h"
14   REAL,             INTENT(INOUT) :: q(ijb_u:ije_u,llm,nqtot)
15   INTEGER,          INTENT(IN)    :: ijb, ije   !--- Can be local and different from ijb_u,ije_u, for example in qminimum
16   CHARACTER(LEN=*), INTENT(IN)    :: err_msg    !--- Error message to display
17   CHARACTER(LEN=maxlen) :: modname, msg1, nm(2)
18   INTEGER :: ixt, ipha, k, i, iq, iiso, izon, ieau, iqeau, iqpar
19   INTEGER, ALLOCATABLE ::   ix(:)
20   REAL,    ALLOCATABLE, SAVE :: tnat(:)         !--- OpenMP shared variable
21   REAL    :: xtractot, xiiso, deltaD, q1, q2
22   REAL, PARAMETER :: borne     = 1e19,  &
23                      errmax    = 1e-8,  &       !--- Max. absolute error
24                      errmaxrel = 1e-3,  &       !--- Max. relative error
25                      qmin      = 1e-11, &
26                      deltaDmax =1000.0, &
27                      deltaDmin =-999.0, &
28                      ridicule  = 1e-12
29   INTEGER, SAVE :: iso_eau, iso_HDO, iso_O18, iso_O17, iso_HTO
30!$OMP THREADPRIVATE(iso_eau, iso_HDO, iso_O18, iso_O17, iso_HTO)
31   LOGICAL       :: ltnat1
32   LOGICAL, SAVE :: first=.TRUE.
33!$OMP THREADPRIVATE(first)
34
35   modname='check_isotopes'
36   IF(.NOT.isoCheck)    RETURN                   !--- No need to check => finished
37   IF(isoSelect('H2O')) RETURN                   !--- No H2O isotopes group found
38   IF(niso == 0)        RETURN                   !--- No isotopes => finished
39   IF(first) THEN
40      ltnat1 = .TRUE.; CALL getin('tnateq1', ltnat1)
41      iso_eau = strIdx(isoName,'H216O'); IF(iso_eau /= 0) tnat(iso_eau) = tnat_H216O
42      iso_O17 = strIdx(isoName,'H217O'); IF(iso_O17 /= 0) tnat(iso_O17) = tnat_H217O
43      iso_O18 = strIdx(isoName,'H218O'); IF(iso_O18 /= 0) tnat(iso_O18) = tnat_H218O
44      iso_HDO = strIdx(isoName,'HDO');   IF(iso_HDO /= 0) tnat(iso_HDO) = tnat_HDO
45      iso_HTO = strIdx(isoName,'HTO');   IF(iso_HTO /= 0) tnat(iso_HTO) = tnat_HTO
46      IF(ltnat1) tnat(:) = 1.0
47      first = .FALSE.
48   END IF
49   CALL msg('31: err_msg='//TRIM(err_msg), modname)
50
51   !--- CHECK FOR NaNs (FOR ALL THE ISOTOPES, INCLUDING GEOGRAPHIC TAGGING TRACERS)
52   modname = 'check_isotopes:iso_verif_noNaN'
53   DO ixt = 1, ntiso
54      DO ipha = 1, nphas
55         iq = iqIsoPha(ixt,ipha)
56!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
57         DO k = 1, llm
58            DO i = ijb, ije
59               IF(ABS(q(i,k,iq))<=borne) CYCLE
60               WRITE(msg1,'(s,"(",i0,",",i0,",",i0,") = ",ES12.4)')TRIM(isoName(ixt)),i,k,iq,q(i,k,iq)
61               CALL msg(msg1, modname)
62               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
63            END DO
64         END DO
65!$OMP END DO NOWAIT
66      END DO
67   END DO
68
69   !--- CHECK CONSERVATION (MAIN ISOTOPE AND PARENT CONCENTRATIONS MUST BE EQUAL)
70   modname = 'check_isotopes:iso_verif_egalite'
71   ixt = iso_eau
72   IF(ixt /= 0) THEN
73      DO ipha = 1, nphas
74         iq = iqIsoPha(ixt,ipha)
75         iqpar = tracers(iq)%iqParent
76!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
77         DO k = 1, llm
78            DO i = ijb, ije
79               q1 = q(i,k,iqpar)
80               q2 = q(i,k,iq)
81!--- IMPROVEMENT in case at least one isotope is not negligible compared to the main isotopic form.
82!    This would be probably required to sum from smallest to highest concentrations ; the corresponding
83!    indices vector can be computed once only (in the initializations stage), using mean concentrations.
84!              q2 = SUM(q(i,k,tracers(iqPar)%iqDesc), DIM=3)
85               IF(ABS(q1-q2) <= errmax .OR. ABS(q1-q2)/MAX(MAX(ABS(q1),ABS(q2)),1e-18) <= errmaxrel) THEN
86                  q(i,k,iq) = q1                 !--- Bidouille pour convergence
87!                 q(i,k,tracers(iqPar)%iqDesc) = q(i,k,tracers(iqPar)%iqDesc) * q1 / q2
88                  CYCLE
89               END IF
90               CALL msg('ixt, iq = '//TRIM(strStack(int2str([ixt,iq]))), modname)
91               msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
92               CALL msg(TRIM(tracers(iqpar)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
93               CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
94               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
95            END DO
96         END DO
97!$OMP END DO NOWAIT
98      END DO
99   END IF
100
101   !--- CHECK DELTA ANOMALIES
102   modname = 'check_isotopes:iso_verif_aberrant'
103   ix = [ iso_HDO  ,   iso_O18 ]
104   nm = ['deltaD  ', 'deltaO18']
105   DO iiso = 1, SIZE(ix)
106      ixt = ix(iiso)
107      IF(ixt  == 0) CYCLE
108      DO ipha = 1, nphas
109         iq = iqIsoPha(ixt,ipha)
110         iqpar = tracers(iq)%iqParent
111!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
112         DO k = 1, llm
113            DO i = ijb, ije
114               q1 = q(i,k,iqpar)
115               q2 = q(i,k,iq)
116!--- IMPROVEMENT in case at least one isotope is not negligible compared to the main isotopic form.
117!    This would be probably required to sum from smallest to highest concentrations ; the corresponding
118!    indices vector can be computed once only (in the initializations stage), using mean concentrations.
119!              q2 = SUM(q(i,k,tracers(iqPar)%iqDesc), DIM=3)
120               IF(q2 <= qmin) CYCLE
121               deltaD = (q2/q1/tnat(ixt)-1.)*1000.
122               IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
123               CALL msg('ixt, iq = '//TRIM(strStack(int2str([ixt,iq]))), modname)
124               msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
125               CALL msg(TRIM(tracers(iqpar)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
126               CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
127               CALL msg(TRIM(nm(iiso))//TRIM(real2str(deltaD)), modname)
128               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
129            END DO
130         END DO
131!$OMP END DO NOWAIT
132      END DO
133   END DO
134
135   IF(nzone == 0) RETURN
136
137   !--- CHECK FOR TAGGING TRACERS DELTAD ANOMALIES
138   modname = 'check_isotopes:iso_verif_aberrant'
139   IF(iso_eau /= 0 .AND. iso_HDO /= 0) THEN
140      DO izon = 1, nzone
141         ixt  = itZonIso(izon, iso_HDO)
142         ieau = itZonIso(izon, iso_eau)
143         DO ipha = 1, nphas
144            iq    = iqIsoPha(ixt,  ipha)
145            iqeau = iqIsoPha(ieau, ipha)
146!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
147            DO k = 1, llm
148               DO i = ijb, ije
149                  q1 = q(i,k,iqeau)
150                  q2 = q(i,k,iq)
151                  IF(q2<=qmin) CYCLE
152                  deltaD = (q2/q1/tnat(iso_HDO)-1.)*1000.
153                  IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
154                  CALL msg('izon, ipha = '//TRIM(strStack(int2str([izon, ipha]))), modname)
155                  CALL msg( 'ixt, ieau = '//TRIM(strStack(int2str([ ixt, ieau]))), modname)
156                  msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
157                  CALL msg(TRIM(tracers(iqeau)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
158                  CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
159                  CALL msg('deltaD = '//TRIM(real2str(deltaD)), modname)
160                  CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
161               END DO
162            END DO
163!$OMP END DO NOWAIT
164         END DO
165      END DO
166   END IF
167
168   !--- CHECK FOR TAGGING TRACERS CONSERVATION (PARENT AND TAGGING TRACERS SUM OVER ALL REGIONS MUST BE EQUAL)
169   DO iiso = 1, niso
170      DO ipha = 1, nphas
171         iq = iqIsoPha(iiso, ipha)
172!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
173         DO k = 1, llm
174            DO i = ijb, ije
175               xiiso = q(i,k,iq)
176               xtractot = SUM(q(i, k, iqIsoPha(itZonIso(1:nzone,iiso), ipha)))
177               IF(ABS(xtractot-xiiso) > errmax .AND. ABS(xtractot-xiiso)/MAX(MAX(ABS(xtractot),ABS(xiiso)),1e-18) > errmaxrel) THEN
178                  CALL msg('Error in iso_verif_aberrant trac: '//TRIM(err_msg))
179                  CALL msg('iiso, ipha = '//TRIM(strStack(int2str([iiso, ipha]))), modname)
180                  CALL msg('q('//TRIM(strStack(int2str([i,k])))//',:) = '//TRIM(strStack(real2str(q(i,k,:)))), modname)
181                  CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
182               END IF
183               IF(ABS(xtractot) <= ridicule) CYCLE
184               DO izon = 1, nzone
185                  q(i,k,iq) = q(i,k,iq) / xtractot * xiiso !--- Bidouille pour convergence
186               END DO
187            END DO
188         END DO
189!$OMP END DO NOWAIT
190      END DO
191   END DO
192
193END SUBROUTINE check_isotopes
194
Note: See TracBrowser for help on using the repository browser.