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