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

Last change on this file since 4159 was 4143, checked in by dcugnet, 3 years ago
  • Some variables are renamed or replaced by direct equivalents:
    • iso_indnum -> tracers(:)%iso_iName
    • niso_possibles -> niso
    • iqiso -> iqIsoPha ; index_trac -> itZonIso
    • ok_iso_verif -> isoCheck
    • ntraceurs_zone -> nzone ; ntraciso -> ntiso
    • qperemin -> min_qparent ; masseqmin -> min_qmass ; ratiomin -> min_ratio
  • Some renamed variables are only aliased with the older name (using USE <module>, ONLY: <oldName> => <newName>) in routines where they are repeated many times.
  • Few hard-coded indexes are now computed (examples: ilic, iso, ivap, irneb, iq_vap, iq_liq, iso_H2O, iso_HDO, iso_HTO, iso_O17, iso_O18).
  • The IF(isoCheck) test is now embedded in the check_isotopes_seq and check_isotopes_loc routines (lighter calling).
File size: 8.3 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, tnat
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    :: xtractot, xiiso, deltaD, q1, q2
15   REAL, PARAMETER :: borne     = 1e19,  &
16                      errmax    = 1e-8,  &       !--- Max. absolute error
17                      errmaxrel = 1e-3,  &       !--- Max. relative error
18                      qmin      = 1e-11, &
19                      deltaDmax =1000.0, &
20                      deltaDmin =-999.0, &
21                      ridicule  = 1e-12
22   INTEGER, SAVE :: iso_eau, iso_HDO, iso_O18, & !--- OpenMP shared variables
23                             iso_O17, iso_HTO
24   LOGICAL, SAVE :: first=.TRUE.
25!$OMP THREADPRIVATE(first)
26
27   modname='check_isotopes'
28   IF(.NOT.isoCheck)    RETURN                   !--- No need to check => finished
29   IF(isoSelect('H2O')) RETURN                   !--- No H2O isotopes group found
30   IF(niso == 0)        RETURN                   !--- No isotopes => finished
31   IF(first) THEN
32!$OMP MASTER
33      iso_eau = strIdx(isoName,'H2[16]O')
34      iso_HDO = strIdx(isoName,'H[2]HO')
35      iso_O18 = strIdx(isoName,'H2[18]O')
36      iso_O17 = strIdx(isoName,'H2[17]O')
37      iso_HTO = strIdx(isoName,'H[3]HO')
38!$OMP END MASTER
39!$OMP BARRIER
40      first = .FALSE.
41   END IF
42   CALL msg('31: err_msg='//TRIM(err_msg), modname)
43
44   !--- CHECK FOR NaNs (FOR ALL THE ISOTOPES, INCLUDING GEOGRAPHIC TAGGING TRACERS)
45   modname = 'check_isotopes:iso_verif_noNaN'
46   DO ixt = 1, ntiso
47      DO ipha = 1, nphas
48         iq = iqIsoPha(ixt,ipha)
49!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
50         DO k = 1, llm
51            DO i = ijb, ije
52               IF(ABS(q(i,k,iq))<=borne) CYCLE
53               WRITE(msg1,'(s,"(",i0,",",i0,",",i0,") = ",ES12.4)')TRIM(isoName(ixt)),i,k,iq,q(i,k,iq)
54               CALL msg(msg1, modname)
55               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
56            END DO
57         END DO
58!$OMP END DO NOWAIT
59      END DO
60   END DO
61
62   !--- CHECK CONSERVATION (MAIN ISOTOPE AND PARENT CONCENTRATIONS MUST BE EQUAL)
63   modname = 'check_isotopes:iso_verif_egalite'
64   ixt = iso_eau
65   IF(ixt /= 0) THEN
66      DO ipha = 1, nphas
67         iq = iqIsoPha(ixt,ipha)
68         iqpar = tracers(iq)%iqParent
69!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
70         DO k = 1, llm
71            DO i = ijb, ije
72               q1 = q(i,k,iqpar)
73               q2 = q(i,k,iq)
74!--- IMPROVEMENT in case at least one isotope is not negligible compared to the main isotopic form.
75!    This would be probably required to sum from smallest to highest concentrations ; the corresponding
76!    indices vector can be computed once only (in the initializations stage), using mean concentrations.
77!              q2 = SUM(q(i,k,tracers(iqPar)%iqDesc), DIM=3)
78               IF(ABS(q1-q2) <= errmax .OR. ABS(q1-q2)/MAX(MAX(ABS(q1),ABS(q2)),1e-18) <= errmaxrel) THEN
79                  q(i,k,iq) = q1                 !--- Bidouille pour convergence
80!                 q(i,k,tracers(iqPar)%iqDesc) = q(i,k,tracers(iqPar)%iqDesc) * q1 / q2
81                  CYCLE
82               END IF
83               CALL msg('ixt, iq = '//TRIM(strStack(int2str([ixt,iq]))), modname)
84               msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
85               CALL msg(TRIM(tracers(iqpar)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
86               CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
87               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
88            END DO
89         END DO
90!$OMP END DO NOWAIT
91      END DO
92   END IF
93
94   !--- CHECK DELTA ANOMALIES
95   modname = 'check_isotopes:iso_verif_aberrant'
96   ix = [ iso_HDO  ,   iso_O18 ]
97   nm = ['deltaD  ', 'deltaO18']
98   DO iiso = 1, SIZE(ix)
99      ixt = ix(iiso)
100      IF(ixt  == 0) CYCLE
101      DO ipha = 1, nphas
102         iq = iqIsoPha(ixt,ipha)
103         iqpar = tracers(iq)%iqParent
104!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
105         DO k = 1, llm
106            DO i = ijb, ije
107               q1 = q(i,k,iqpar)
108               q2 = q(i,k,iq)
109!--- IMPROVEMENT in case at least one isotope is not negligible compared to the main isotopic form.
110!    This would be probably required to sum from smallest to highest concentrations ; the corresponding
111!    indices vector can be computed once only (in the initializations stage), using mean concentrations.
112!              q2 = SUM(q(i,k,tracers(iqPar)%iqDesc), DIM=3)
113               IF(q2 <= qmin) CYCLE
114               deltaD = (q2/q1/tnat(ixt)-1.)*1000.
115               IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
116               CALL msg('ixt, iq = '//TRIM(strStack(int2str([ixt,iq]))), modname)
117               msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
118               CALL msg(TRIM(tracers(iqpar)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
119               CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
120               CALL msg(TRIM(nm(iiso))//TRIM(real2str(deltaD)), modname)
121               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
122            END DO
123         END DO
124!$OMP END DO NOWAIT
125      END DO
126   END DO
127
128   IF(nzone == 0) RETURN
129
130   !--- CHECK FOR TAGGING TRACERS DELTAD ANOMALIES
131   modname = 'check_isotopes:iso_verif_aberrant'
132   IF(iso_eau /= 0 .AND. iso_HDO /= 0) THEN
133      DO izon = 1, nzone
134         ixt  = itZonIso(izon, iso_HDO)
135         ieau = itZonIso(izon, iso_eau)
136         DO ipha = 1, nphas
137            iq    = iqIsoPha(ixt,  ipha)
138            iqeau = iqIsoPha(ieau, ipha)
139!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
140            DO k = 1, llm
141               DO i = ijb, ije
142                  q1 = q(i,k,iqeau)
143                  q2 = q(i,k,iq)
144                  IF(q2<=qmin) CYCLE
145                  deltaD = (q2/q1/tnat(iso_HDO)-1.)*1000.
146                  IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
147                  CALL msg('izon, ipha = '//TRIM(strStack(int2str([izon, ipha]))), modname)
148                  CALL msg( 'ixt, ieau = '//TRIM(strStack(int2str([ ixt, ieau]))), modname)
149                  msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
150                  CALL msg(TRIM(tracers(iqeau)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
151                  CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
152                  CALL msg('deltaD = '//TRIM(real2str(deltaD)), modname)
153                  CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
154               END DO
155            END DO
156!$OMP END DO NOWAIT
157         END DO
158      END DO
159   END IF
160
161   !--- CHECK FOR TAGGING TRACERS CONSERVATION (PARENT AND TAGGING TRACERS SUM OVER ALL REGIONS MUST BE EQUAL)
162   DO iiso = 1, niso
163      DO ipha = 1, nphas
164         iq = iqIsoPha(iiso, ipha)
165!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
166         DO k = 1, llm
167            DO i = ijb, ije
168               xiiso = q(i,k,iq)
169               xtractot = SUM(q(i, k, iqIsoPha(itZonIso(1:nzone,iiso), ipha)))
170               IF(ABS(xtractot-xiiso) > errmax .AND. ABS(xtractot-xiiso)/MAX(MAX(ABS(xtractot),ABS(xiiso)),1e-18) > errmaxrel) THEN
171                  CALL msg('Error in iso_verif_aberrant trac: '//TRIM(err_msg))
172                  CALL msg('iiso, ipha = '//TRIM(strStack(int2str([iiso, ipha]))), modname)
173                  CALL msg('q('//TRIM(strStack(int2str([i,k])))//',:) = '//TRIM(strStack(real2str(q(i,k,:)))), modname)
174                  CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
175               END IF
176               IF(ABS(xtractot) <= ridicule) CYCLE
177               DO izon = 1, nzone
178                  q(i,k,iq) = q(i,k,iq) / xtractot * xiiso !--- Bidouille pour convergence
179               END DO
180            END DO
181         END DO
182!$OMP END DO NOWAIT
183      END DO
184   END DO
185
186END SUBROUTINE check_isotopes
187
Note: See TracBrowser for help on using the repository browser.