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

Last change on this file since 4608 was 4399, checked in by dcugnet, 23 months ago

Change the names of the water isotopes: no more "and?" separators for atomic numbers.

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