source: LMDZ6/trunk/libf/dyn3d/check_isotopes.F90 @ 4175

Last change on this file since 4175 was 4143, checked in by dcugnet, 2 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: 7.8 KB
Line 
1SUBROUTINE check_isotopes_seq(q, ip1jmp1, err_msg)
2   USE strings_mod, ONLY: maxlen, msg, strIdx, strStack, int2str, real2str
3   USE infotrac,    ONLY: nqtot, niso, nphas, isotope, isoCheck, iqIsoPha, isoSelect, &
4                          ntiso, iH2O, nzone, tracers, isoName,  itZonIso, tnat
5   IMPLICIT NONE
6   include "dimensions.h"
7   REAL,             INTENT(INOUT) :: q(ip1jmp1,llm,nqtot)
8   INTEGER,          INTENT(IN)    :: ip1jmp1
9   CHARACTER(LEN=*), INTENT(IN)    :: err_msg    !--- Error message to display
10   CHARACTER(LEN=maxlen) :: modname, msg1, nm(2)
11   INTEGER :: ixt, ipha, k, i, iq, iiso, izon, ieau, iqeau, iqpar
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-3,  &       !--- Max. relative error
17                      qmin      = 1e-11, &
18                      deltaDmax =1000.0, &
19                      deltaDmin =-999.0, &
20                      ridicule  = 1e-12
21   INTEGER, SAVE :: iso_eau, iso_HDO, iso_O18, &
22                             iso_O17, iso_HTO
23   LOGICAL, SAVE :: first=.TRUE.
24
25   modname='check_isotopes'
26   IF(.NOT.isoCheck)    RETURN                   !--- No need to check => finished
27   IF(isoSelect('H2O')) RETURN                   !--- No H2O isotopes group found
28   IF(niso == 0)        RETURN                   !--- No isotopes => finished
29   IF(first) THEN
30      iso_eau = strIdx(isoName,'H2[16]O')
31      iso_HDO = strIdx(isoName,'H[2]HO')
32      iso_O18 = strIdx(isoName,'H2[18]O')
33      iso_O17 = strIdx(isoName,'H2[17]O')
34      iso_HTO = strIdx(isoName,'H[3]HO')
35      first = .FALSE.
36   END IF
37   CALL msg('31: err_msg='//TRIM(err_msg), modname)
38
39   !--- CHECK FOR NaNs (FOR ALL THE ISOTOPES, INCLUDING GEOGRAPHIC TAGGING TRACERS)
40   modname = 'check_isotopes:iso_verif_noNaN'
41   DO ixt = 1, ntiso
42      DO ipha = 1, nphas
43         iq = iqIsoPha(ixt,ipha)
44         DO k = 1, llm
45            DO i = 1, ip1jmp1
46               IF(ABS(q(i,k,iq)) < borne) CYCLE
47               WRITE(msg1,'(s,"(",i0,",",i0,",",i0,") = ",ES12.4)')TRIM(isoName(ixt)),i,k,iq,q(i,k,iq)
48               CALL msg(msg1, modname)
49               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
50            END DO
51         END DO
52      END DO
53   END DO
54
55   !--- CHECK CONSERVATION (MAIN ISOTOPE AND PARENT CONCENTRATIONS MUST BE EQUAL)
56   modname = 'check_isotopes:iso_verif_egalite'
57   ixt = iso_eau
58   IF(ixt /= 0) THEN
59      DO ipha = 1, nphas
60         iq = iqIsoPha(ixt,ipha)
61         iqpar = tracers(iq)%iqParent
62         DO k = 1, llm
63            DO i = 1, ip1jmp1
64               q1 = q(i,k,iqpar)
65               q2 = q(i,k,iq)
66!--- IMPROVEMENT in case at least one isotope is not negligible compared to the main isotopic form.
67!    This would be probably required to sum from smallest to highest concentrations ; the corresponding
68!    indices vector can be computed once only (in the initializations stage), using mean concentrations.
69!              q2 = SUM(q(i,k,tracers(iqPar)%iqDesc), DIM=3)
70               IF(ABS(q1-q2) <= errmax .OR. ABS(q1-q2)/MAX(MAX(ABS(q1),ABS(q2)),1e-18) <= errmaxrel) THEN
71                  q(i,k,iq) = q1                 !--- Bidouille pour convergence
72!                 q(i,k,tracers(iqPar)%iqDesc) = q(i,k,tracers(iqPar)%iqDesc) * q1 / q2
73                  CYCLE
74               END IF
75               CALL msg('ixt, iq = '//TRIM(strStack(int2str([ixt,iq]))), modname)
76               msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
77               CALL msg(TRIM(tracers(iqpar)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
78               CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
79               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
80            END DO
81         END DO
82      END DO
83   END IF
84
85   !--- CHECK DELTA ANOMALIES
86   modname = 'check_isotopes:iso_verif_aberrant'
87   ix = [ iso_HDO  ,   iso_O18 ]
88   nm = ['deltaD  ', 'deltaO18']
89   DO iiso = 1, SIZE(ix)
90      ixt = ix(iiso)
91      IF(ixt  == 0) CYCLE
92      DO ipha = 1, nphas
93         iq = iqIsoPha(ixt,ipha)
94         iqpar = tracers(iq)%iqParent
95         DO k = 1, llm
96            DO i = 1, ip1jmp1
97               q1 = q(i,k,iqpar)
98               q2 = q(i,k,iq)
99!--- IMPROVEMENT in case at least one isotope is not negligible compared to the main isotopic form.
100!    This would be probably required to sum from smallest to highest concentrations ; the corresponding
101!    indices vector can be computed once only (in the initializations stage), using mean concentrations.
102!              q2 = SUM(q(i,k,tracers(iqPar)%iqDesc), DIM=3)
103               IF(q2 <= qmin) CYCLE
104               deltaD = (q2/q1/tnat(ixt)-1.)*1000.
105               IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
106               CALL msg('ixt, iq = '//TRIM(strStack(int2str([ixt,iq]))), modname)
107               msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
108               CALL msg(TRIM(tracers(iqpar)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
109               CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
110               CALL msg(TRIM(nm(iiso))//TRIM(real2str(deltaD)), modname)
111               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
112            END DO
113         END DO
114      END DO
115   END DO
116
117   IF(nzone == 0) RETURN
118
119   !--- CHECK FOR TAGGING TRACERS DELTAD ANOMALIES
120   modname = 'check_isotopes:iso_verif_aberrant'
121   IF(iso_eau /= 0 .AND. iso_HDO /= 0) THEN
122      DO izon = 1, nzone
123         ixt  = itZonIso(izon, iso_HDO)
124         ieau = itZonIso(izon, iso_eau)
125         DO ipha = 1, nphas
126            iq    = iqIsoPha(ixt,  ipha)
127            iqeau = iqIsoPha(ieau, ipha)
128            DO k = 1, llm
129               DO i = 1, ip1jmp1
130                  q1 = q(i,k,iqeau)
131                  q2 = q(i,k,iq)
132                  IF(q2<=qmin) CYCLE
133                  deltaD = (q2/q1/tnat(iso_HDO)-1.)*1000.
134                  IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
135                  CALL msg('izon, ipha = '//TRIM(strStack(int2str([izon, ipha]))), modname)
136                  CALL msg( 'ixt, ieau = '//TRIM(strStack(int2str([ ixt, ieau]))), modname)
137                  msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
138                  CALL msg(TRIM(tracers(iqeau)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
139                  CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
140                  CALL msg('deltaD = '//TRIM(real2str(deltaD)), modname)
141                  CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
142               END DO
143            END DO
144         END DO
145      END DO
146   END IF
147
148   !--- CHECK FOR TAGGING TRACERS CONSERVATION (PARENT AND TAGGING TRACERS SUM OVER ALL REGIONS MUST BE EQUAL)
149   DO iiso = 1, niso
150      DO ipha = 1, nphas
151         iq = iqIsoPha(iiso, ipha)
152         DO k = 1, llm
153            DO i = 1, ip1jmp1
154               xiiso = q(i,k,iq)
155               xtractot = SUM(q(i, k, iqIsoPha(itZonIso(1:nzone,iiso), ipha)))
156               IF(ABS(xtractot-xiiso) > errmax .AND. ABS(xtractot-xiiso)/MAX(MAX(ABS(xtractot),ABS(xiiso)),1e-18) > errmaxrel) THEN
157                  CALL msg('Error in iso_verif_aberrant trac: '//TRIM(err_msg))
158                  CALL msg('iiso, ipha = '//TRIM(strStack(int2str([iiso, ipha]))), modname)
159                  CALL msg('q('//TRIM(strStack(int2str([i,k])))//',:) = '//TRIM(strStack(real2str(q(i,k,:)))), modname)
160                  CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
161               END IF
162               IF(ABS(xtractot) <= ridicule) CYCLE
163               DO izon = 1, nzone
164                  q(i,k,iq) = q(i,k,iq) / xtractot * xiiso !--- Bidouille pour convergence
165               END DO
166            END DO
167         END DO
168      END DO
169   END DO
170
171END SUBROUTINE check_isotopes_seq
172
Note: See TracBrowser for help on using the repository browser.