1 |
|
15273 |
SUBROUTINE 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 |
|
15273 |
modname='check_isotopes' |
27 |
✗✓ |
15273 |
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 |
|
|
|
173 |
✗✗✗✓
|
15273 |
END SUBROUTINE check_isotopes_seq |
174 |
|
|
|