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

Last change on this file since 4359 was 4325, checked in by dcugnet, 2 years ago
  • simplify the parser usage:
    • the getKey_init routine is now embedded in the readTracersFile routine.
    • the initIsotopes routine is now embedded in the readIsotopesFile routine.
    • the database is now unique, but can be changed using the get/setKeysDBase.
    • the derived types descriptions, originally located in trac_types_mod, are moved to readTracFiles_mod.
    • few checkings moved from infotrac to the routine testIsotopes, contained in the readIsotopesFile function from readTracFiles_mod.
    • the readTracersFiles and readIsotopesFile routines no longer use a tracers/isotopes argument.
  • remove tnat and alpha_ideal from infotrac ; use instead getKey to get them where they are used (check_isotopes, dynetat0, iniacademic)
  • the trac_type field %Childs is renamed %Children
  • move the isoSelect routine and the corresponding variables routine from infotrac and infotrac_phy to readTracFiles_mod
  • infotrac_phy routine is now fully independant of the (very similar) routine infotrac (init_infotrac_phy has no arguments left).
  • all the explicit keys of the trac_type are now included in the embedded keys database, accessible using the getKey function.
  • the getKey/addKey routines are expanded to handle vectors of integers, reals, logicals or strings.
  • few subroutines converted into functions with error return value.
  • corrections for isotopic tagging tracers mode (to be continued).
File size: 8.4 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, getKey
6   IMPLICIT NONE
7   include "dimensions.h"
8   REAL,             INTENT(INOUT) :: q(ijb_u:ije_u,llm,nqtot)
9   INTEGER,          INTENT(IN)    :: ijb, ije   !--- Can be local and different from ijb_u,ije_u, for example in qminimum
10   CHARACTER(LEN=*), INTENT(IN)    :: err_msg    !--- Error message to display
11   CHARACTER(LEN=maxlen) :: modname, msg1, nm(2)
12   INTEGER :: ixt, ipha, k, i, iq, iiso, izon, ieau, iqeau, iqpar
13   INTEGER, ALLOCATABLE ::   ix(:)
14   REAL,    ALLOCATABLE :: tnat(:)               !--- OpenMP shared variable
15   REAL    :: xtractot, xiiso, deltaD, q1, q2
16   REAL, PARAMETER :: borne     = 1e19,  &
17                      errmax    = 1e-8,  &       !--- Max. absolute error
18                      errmaxrel = 1e-3,  &       !--- Max. relative error
19                      qmin      = 1e-11, &
20                      deltaDmax =1000.0, &
21                      deltaDmin =-999.0, &
22                      ridicule  = 1e-12
23   INTEGER, SAVE :: iso_eau, iso_HDO, iso_O18, & !--- OpenMP shared variables
24                             iso_O17, iso_HTO
25   LOGICAL, SAVE :: first=.TRUE.
26!$OMP THREADPRIVATE(first)
27
28   modname='check_isotopes'
29   IF(.NOT.isoCheck)    RETURN                   !--- No need to check => finished
30   IF(isoSelect('H2O')) RETURN                   !--- No H2O isotopes group found
31   IF(niso == 0)        RETURN                   !--- No isotopes => finished
32   IF(first) THEN
33!$OMP MASTER
34      iso_eau = strIdx(isoName,'H2[16]O')
35      iso_HDO = strIdx(isoName,'H[2]HO')
36      iso_O18 = strIdx(isoName,'H2[18]O')
37      iso_O17 = strIdx(isoName,'H2[17]O')
38      iso_HTO = strIdx(isoName,'H[3]HO')
39      IF(getKey('tnat', tnat)) CALL abort_gcm(modname, 'missing isotopic parameter', 1)
40!$OMP END MASTER
41!$OMP BARRIER
42      first = .FALSE.
43   END IF
44   CALL msg('31: err_msg='//TRIM(err_msg), modname)
45
46   !--- CHECK FOR NaNs (FOR ALL THE ISOTOPES, INCLUDING GEOGRAPHIC TAGGING TRACERS)
47   modname = 'check_isotopes:iso_verif_noNaN'
48   DO ixt = 1, ntiso
49      DO ipha = 1, nphas
50         iq = iqIsoPha(ixt,ipha)
51!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
52         DO k = 1, llm
53            DO i = ijb, ije
54               IF(ABS(q(i,k,iq))<=borne) CYCLE
55               WRITE(msg1,'(s,"(",i0,",",i0,",",i0,") = ",ES12.4)')TRIM(isoName(ixt)),i,k,iq,q(i,k,iq)
56               CALL msg(msg1, modname)
57               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
58            END DO
59         END DO
60!$OMP END DO NOWAIT
61      END DO
62   END DO
63
64   !--- CHECK CONSERVATION (MAIN ISOTOPE AND PARENT CONCENTRATIONS MUST BE EQUAL)
65   modname = 'check_isotopes:iso_verif_egalite'
66   ixt = iso_eau
67   IF(ixt /= 0) THEN
68      DO ipha = 1, nphas
69         iq = iqIsoPha(ixt,ipha)
70         iqpar = tracers(iq)%iqParent
71!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
72         DO k = 1, llm
73            DO i = ijb, ije
74               q1 = q(i,k,iqpar)
75               q2 = q(i,k,iq)
76!--- IMPROVEMENT in case at least one isotope is not negligible compared to the main isotopic form.
77!    This would be probably required to sum from smallest to highest concentrations ; the corresponding
78!    indices vector can be computed once only (in the initializations stage), using mean concentrations.
79!              q2 = SUM(q(i,k,tracers(iqPar)%iqDesc), DIM=3)
80               IF(ABS(q1-q2) <= errmax .OR. ABS(q1-q2)/MAX(MAX(ABS(q1),ABS(q2)),1e-18) <= errmaxrel) THEN
81                  q(i,k,iq) = q1                 !--- Bidouille pour convergence
82!                 q(i,k,tracers(iqPar)%iqDesc) = q(i,k,tracers(iqPar)%iqDesc) * q1 / q2
83                  CYCLE
84               END IF
85               CALL msg('ixt, iq = '//TRIM(strStack(int2str([ixt,iq]))), modname)
86               msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
87               CALL msg(TRIM(tracers(iqpar)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
88               CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
89               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
90            END DO
91         END DO
92!$OMP END DO NOWAIT
93      END DO
94   END IF
95
96   !--- CHECK DELTA ANOMALIES
97   modname = 'check_isotopes:iso_verif_aberrant'
98   ix = [ iso_HDO  ,   iso_O18 ]
99   nm = ['deltaD  ', 'deltaO18']
100   DO iiso = 1, SIZE(ix)
101      ixt = ix(iiso)
102      IF(ixt  == 0) CYCLE
103      DO ipha = 1, nphas
104         iq = iqIsoPha(ixt,ipha)
105         iqpar = tracers(iq)%iqParent
106!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
107         DO k = 1, llm
108            DO i = ijb, ije
109               q1 = q(i,k,iqpar)
110               q2 = q(i,k,iq)
111!--- IMPROVEMENT in case at least one isotope is not negligible compared to the main isotopic form.
112!    This would be probably required to sum from smallest to highest concentrations ; the corresponding
113!    indices vector can be computed once only (in the initializations stage), using mean concentrations.
114!              q2 = SUM(q(i,k,tracers(iqPar)%iqDesc), DIM=3)
115               IF(q2 <= qmin) CYCLE
116               deltaD = (q2/q1/tnat(ixt)-1.)*1000.
117               IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
118               CALL msg('ixt, iq = '//TRIM(strStack(int2str([ixt,iq]))), modname)
119               msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
120               CALL msg(TRIM(tracers(iqpar)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
121               CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
122               CALL msg(TRIM(nm(iiso))//TRIM(real2str(deltaD)), modname)
123               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
124            END DO
125         END DO
126!$OMP END DO NOWAIT
127      END DO
128   END DO
129
130   IF(nzone == 0) RETURN
131
132   !--- CHECK FOR TAGGING TRACERS DELTAD ANOMALIES
133   modname = 'check_isotopes:iso_verif_aberrant'
134   IF(iso_eau /= 0 .AND. iso_HDO /= 0) THEN
135      DO izon = 1, nzone
136         ixt  = itZonIso(izon, iso_HDO)
137         ieau = itZonIso(izon, iso_eau)
138         DO ipha = 1, nphas
139            iq    = iqIsoPha(ixt,  ipha)
140            iqeau = iqIsoPha(ieau, ipha)
141!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
142            DO k = 1, llm
143               DO i = ijb, ije
144                  q1 = q(i,k,iqeau)
145                  q2 = q(i,k,iq)
146                  IF(q2<=qmin) CYCLE
147                  deltaD = (q2/q1/tnat(iso_HDO)-1.)*1000.
148                  IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
149                  CALL msg('izon, ipha = '//TRIM(strStack(int2str([izon, ipha]))), modname)
150                  CALL msg( 'ixt, ieau = '//TRIM(strStack(int2str([ ixt, ieau]))), modname)
151                  msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
152                  CALL msg(TRIM(tracers(iqeau)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
153                  CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
154                  CALL msg('deltaD = '//TRIM(real2str(deltaD)), modname)
155                  CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
156               END DO
157            END DO
158!$OMP END DO NOWAIT
159         END DO
160      END DO
161   END IF
162
163   !--- CHECK FOR TAGGING TRACERS CONSERVATION (PARENT AND TAGGING TRACERS SUM OVER ALL REGIONS MUST BE EQUAL)
164   DO iiso = 1, niso
165      DO ipha = 1, nphas
166         iq = iqIsoPha(iiso, ipha)
167!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
168         DO k = 1, llm
169            DO i = ijb, ije
170               xiiso = q(i,k,iq)
171               xtractot = SUM(q(i, k, iqIsoPha(itZonIso(1:nzone,iiso), ipha)))
172               IF(ABS(xtractot-xiiso) > errmax .AND. ABS(xtractot-xiiso)/MAX(MAX(ABS(xtractot),ABS(xiiso)),1e-18) > errmaxrel) THEN
173                  CALL msg('Error in iso_verif_aberrant trac: '//TRIM(err_msg))
174                  CALL msg('iiso, ipha = '//TRIM(strStack(int2str([iiso, ipha]))), modname)
175                  CALL msg('q('//TRIM(strStack(int2str([i,k])))//',:) = '//TRIM(strStack(real2str(q(i,k,:)))), modname)
176                  CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
177               END IF
178               IF(ABS(xtractot) <= ridicule) CYCLE
179               DO izon = 1, nzone
180                  q(i,k,iq) = q(i,k,iq) / xtractot * xiiso !--- Bidouille pour convergence
181               END DO
182            END DO
183         END DO
184!$OMP END DO NOWAIT
185      END DO
186   END DO
187
188END SUBROUTINE check_isotopes
189
Note: See TracBrowser for help on using the repository browser.