[4143] | 1 | SUBROUTINE 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, & |
---|
[5214] | 5 | ntiso, iH2O, nzone, tracers, isoName, itZonIso |
---|
| 6 | USE iso_params_mod, ONLY: tnat_H216O, tnat_H217O, tnat_H218O, tnat_HDO, tnat_HTO |
---|
[5200] | 7 | #ifdef CPP_IOIPSL |
---|
| 8 | USE ioipsl, ONLY: getin |
---|
| 9 | #else |
---|
| 10 | USE ioipsl_getincom, ONLY: getin |
---|
| 11 | #endif |
---|
[4143] | 12 | IMPLICIT NONE |
---|
| 13 | include "dimensions.h" |
---|
| 14 | REAL, INTENT(INOUT) :: q(ijb_u:ije_u,llm,nqtot) |
---|
| 15 | INTEGER, INTENT(IN) :: ijb, ije !--- Can be local and different from ijb_u,ije_u, for example in qminimum |
---|
| 16 | CHARACTER(LEN=*), INTENT(IN) :: err_msg !--- Error message to display |
---|
| 17 | CHARACTER(LEN=maxlen) :: modname, msg1, nm(2) |
---|
| 18 | INTEGER :: ixt, ipha, k, i, iq, iiso, izon, ieau, iqeau, iqpar |
---|
[5214] | 19 | INTEGER, ALLOCATABLE :: ix(:) |
---|
[4367] | 20 | REAL, ALLOCATABLE, SAVE :: tnat(:) !--- OpenMP shared variable |
---|
[5214] | 21 | REAL :: xtractot, xiiso, deltaD, q1, q2 |
---|
[4143] | 22 | REAL, PARAMETER :: borne = 1e19, & |
---|
| 23 | errmax = 1e-8, & !--- Max. absolute error |
---|
| 24 | errmaxrel = 1e-3, & !--- Max. relative error |
---|
| 25 | qmin = 1e-11, & |
---|
| 26 | deltaDmax =1000.0, & |
---|
| 27 | deltaDmin =-999.0, & |
---|
| 28 | ridicule = 1e-12 |
---|
[5200] | 29 | INTEGER, SAVE :: iso_eau, iso_HDO, iso_O18, iso_O17, iso_HTO |
---|
| 30 | !$OMP THREADPRIVATE(iso_eau, iso_HDO, iso_O18, iso_O17, iso_HTO) |
---|
| 31 | LOGICAL :: ltnat1 |
---|
[4143] | 32 | LOGICAL, SAVE :: first=.TRUE. |
---|
| 33 | !$OMP THREADPRIVATE(first) |
---|
[2270] | 34 | |
---|
[4143] | 35 | modname='check_isotopes' |
---|
| 36 | IF(.NOT.isoCheck) RETURN !--- No need to check => finished |
---|
| 37 | IF(isoSelect('H2O')) RETURN !--- No H2O isotopes group found |
---|
| 38 | IF(niso == 0) RETURN !--- No isotopes => finished |
---|
| 39 | IF(first) THEN |
---|
[5190] | 40 | !$OMP MASTER |
---|
[5200] | 41 | ltnat1 = .TRUE.; CALL getin('tnateq1', ltnat1) |
---|
| 42 | ALLOCATE(tnat(niso)) |
---|
[5214] | 43 | iso_eau = strIdx(isoName,'H216O'); IF(iso_eau /= 0) tnat(iso_eau) = tnat_H216O |
---|
| 44 | iso_O17 = strIdx(isoName,'H217O'); IF(iso_O17 /= 0) tnat(iso_O17) = tnat_H217O |
---|
| 45 | iso_O18 = strIdx(isoName,'H218O'); IF(iso_O18 /= 0) tnat(iso_O18) = tnat_H218O |
---|
| 46 | iso_HDO = strIdx(isoName,'HDO'); IF(iso_HDO /= 0) tnat(iso_HDO) = tnat_HDO |
---|
| 47 | iso_HTO = strIdx(isoName,'HTO'); IF(iso_HTO /= 0) tnat(iso_HTO) = tnat_HTO |
---|
| 48 | IF(ltnat1) tnat(:) = 1.0 |
---|
[5190] | 49 | !$OMP END MASTER |
---|
| 50 | !$OMP BARRIER |
---|
[4143] | 51 | first = .FALSE. |
---|
| 52 | END IF |
---|
| 53 | CALL msg('31: err_msg='//TRIM(err_msg), modname) |
---|
[2270] | 54 | |
---|
[4143] | 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 | !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) |
---|
| 61 | DO k = 1, llm |
---|
| 62 | DO i = ijb, ije |
---|
[5214] | 63 | IF(ABS(q(i,k,iq)) < borne) CYCLE |
---|
[4143] | 64 | WRITE(msg1,'(s,"(",i0,",",i0,",",i0,") = ",ES12.4)')TRIM(isoName(ixt)),i,k,iq,q(i,k,iq) |
---|
| 65 | CALL msg(msg1, modname) |
---|
| 66 | CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1) |
---|
| 67 | END DO |
---|
| 68 | END DO |
---|
| 69 | !$OMP END DO NOWAIT |
---|
| 70 | END DO |
---|
| 71 | END DO |
---|
[2270] | 72 | |
---|
[4143] | 73 | !--- CHECK CONSERVATION (MAIN ISOTOPE AND PARENT CONCENTRATIONS MUST BE EQUAL) |
---|
| 74 | modname = 'check_isotopes:iso_verif_egalite' |
---|
| 75 | ixt = iso_eau |
---|
| 76 | IF(ixt /= 0) THEN |
---|
| 77 | DO ipha = 1, nphas |
---|
| 78 | iq = iqIsoPha(ixt,ipha) |
---|
| 79 | iqpar = tracers(iq)%iqParent |
---|
| 80 | !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) |
---|
| 81 | DO k = 1, llm |
---|
| 82 | DO i = ijb, ije |
---|
| 83 | q1 = q(i,k,iqpar) |
---|
| 84 | q2 = q(i,k,iq) |
---|
| 85 | !--- IMPROVEMENT in case at least one isotope is not negligible compared to the main isotopic form. |
---|
| 86 | ! This would be probably required to sum from smallest to highest concentrations ; the corresponding |
---|
| 87 | ! indices vector can be computed once only (in the initializations stage), using mean concentrations. |
---|
| 88 | ! q2 = SUM(q(i,k,tracers(iqPar)%iqDesc), DIM=3) |
---|
| 89 | IF(ABS(q1-q2) <= errmax .OR. ABS(q1-q2)/MAX(MAX(ABS(q1),ABS(q2)),1e-18) <= errmaxrel) THEN |
---|
| 90 | q(i,k,iq) = q1 !--- Bidouille pour convergence |
---|
| 91 | ! q(i,k,tracers(iqPar)%iqDesc) = q(i,k,tracers(iqPar)%iqDesc) * q1 / q2 |
---|
| 92 | CYCLE |
---|
| 93 | END IF |
---|
| 94 | CALL msg('ixt, iq = '//TRIM(strStack(int2str([ixt,iq]))), modname) |
---|
| 95 | msg1 = '('//TRIM(strStack(int2str([i,k])))//')' |
---|
| 96 | CALL msg(TRIM(tracers(iqpar)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname) |
---|
| 97 | CALL msg(TRIM(tracers(iq )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname) |
---|
| 98 | CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1) |
---|
| 99 | END DO |
---|
| 100 | END DO |
---|
| 101 | !$OMP END DO NOWAIT |
---|
| 102 | END DO |
---|
| 103 | END IF |
---|
[2270] | 104 | |
---|
[4143] | 105 | !--- CHECK DELTA ANOMALIES |
---|
| 106 | modname = 'check_isotopes:iso_verif_aberrant' |
---|
| 107 | ix = [ iso_HDO , iso_O18 ] |
---|
| 108 | nm = ['deltaD ', 'deltaO18'] |
---|
| 109 | DO iiso = 1, SIZE(ix) |
---|
| 110 | ixt = ix(iiso) |
---|
| 111 | IF(ixt == 0) CYCLE |
---|
| 112 | DO ipha = 1, nphas |
---|
| 113 | iq = iqIsoPha(ixt,ipha) |
---|
| 114 | iqpar = tracers(iq)%iqParent |
---|
| 115 | !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) |
---|
| 116 | DO k = 1, llm |
---|
| 117 | DO i = ijb, ije |
---|
| 118 | q1 = q(i,k,iqpar) |
---|
| 119 | q2 = q(i,k,iq) |
---|
| 120 | !--- IMPROVEMENT in case at least one isotope is not negligible compared to the main isotopic form. |
---|
| 121 | ! This would be probably required to sum from smallest to highest concentrations ; the corresponding |
---|
| 122 | ! indices vector can be computed once only (in the initializations stage), using mean concentrations. |
---|
| 123 | ! q2 = SUM(q(i,k,tracers(iqPar)%iqDesc), DIM=3) |
---|
| 124 | IF(q2 <= qmin) CYCLE |
---|
| 125 | deltaD = (q2/q1/tnat(ixt)-1.)*1000. |
---|
| 126 | IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE |
---|
| 127 | CALL msg('ixt, iq = '//TRIM(strStack(int2str([ixt,iq]))), modname) |
---|
| 128 | msg1 = '('//TRIM(strStack(int2str([i,k])))//')' |
---|
| 129 | CALL msg(TRIM(tracers(iqpar)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname) |
---|
| 130 | CALL msg(TRIM(tracers(iq )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname) |
---|
| 131 | CALL msg(TRIM(nm(iiso))//TRIM(real2str(deltaD)), modname) |
---|
| 132 | CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1) |
---|
| 133 | END DO |
---|
| 134 | END DO |
---|
| 135 | !$OMP END DO NOWAIT |
---|
| 136 | END DO |
---|
| 137 | END DO |
---|
[2270] | 138 | |
---|
[4143] | 139 | IF(nzone == 0) RETURN |
---|
[2270] | 140 | |
---|
[4143] | 141 | !--- CHECK FOR TAGGING TRACERS DELTAD ANOMALIES |
---|
| 142 | modname = 'check_isotopes:iso_verif_aberrant' |
---|
| 143 | IF(iso_eau /= 0 .AND. iso_HDO /= 0) THEN |
---|
| 144 | DO izon = 1, nzone |
---|
| 145 | ixt = itZonIso(izon, iso_HDO) |
---|
| 146 | ieau = itZonIso(izon, iso_eau) |
---|
| 147 | DO ipha = 1, nphas |
---|
| 148 | iq = iqIsoPha(ixt, ipha) |
---|
| 149 | iqeau = iqIsoPha(ieau, ipha) |
---|
| 150 | !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) |
---|
| 151 | DO k = 1, llm |
---|
| 152 | DO i = ijb, ije |
---|
| 153 | q1 = q(i,k,iqeau) |
---|
| 154 | q2 = q(i,k,iq) |
---|
| 155 | IF(q2<=qmin) CYCLE |
---|
| 156 | deltaD = (q2/q1/tnat(iso_HDO)-1.)*1000. |
---|
| 157 | IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE |
---|
| 158 | CALL msg('izon, ipha = '//TRIM(strStack(int2str([izon, ipha]))), modname) |
---|
| 159 | CALL msg( 'ixt, ieau = '//TRIM(strStack(int2str([ ixt, ieau]))), modname) |
---|
| 160 | msg1 = '('//TRIM(strStack(int2str([i,k])))//')' |
---|
| 161 | CALL msg(TRIM(tracers(iqeau)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname) |
---|
| 162 | CALL msg(TRIM(tracers(iq )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname) |
---|
| 163 | CALL msg('deltaD = '//TRIM(real2str(deltaD)), modname) |
---|
| 164 | CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1) |
---|
| 165 | END DO |
---|
| 166 | END DO |
---|
| 167 | !$OMP END DO NOWAIT |
---|
| 168 | END DO |
---|
| 169 | END DO |
---|
| 170 | END IF |
---|
[2270] | 171 | |
---|
[4143] | 172 | !--- CHECK FOR TAGGING TRACERS CONSERVATION (PARENT AND TAGGING TRACERS SUM OVER ALL REGIONS MUST BE EQUAL) |
---|
| 173 | DO iiso = 1, niso |
---|
| 174 | DO ipha = 1, nphas |
---|
| 175 | iq = iqIsoPha(iiso, ipha) |
---|
| 176 | !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) |
---|
| 177 | DO k = 1, llm |
---|
| 178 | DO i = ijb, ije |
---|
| 179 | xiiso = q(i,k,iq) |
---|
| 180 | xtractot = SUM(q(i, k, iqIsoPha(itZonIso(1:nzone,iiso), ipha))) |
---|
| 181 | IF(ABS(xtractot-xiiso) > errmax .AND. ABS(xtractot-xiiso)/MAX(MAX(ABS(xtractot),ABS(xiiso)),1e-18) > errmaxrel) THEN |
---|
| 182 | CALL msg('Error in iso_verif_aberrant trac: '//TRIM(err_msg)) |
---|
| 183 | CALL msg('iiso, ipha = '//TRIM(strStack(int2str([iiso, ipha]))), modname) |
---|
| 184 | CALL msg('q('//TRIM(strStack(int2str([i,k])))//',:) = '//TRIM(strStack(real2str(q(i,k,:)))), modname) |
---|
| 185 | CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1) |
---|
| 186 | END IF |
---|
| 187 | IF(ABS(xtractot) <= ridicule) CYCLE |
---|
| 188 | DO izon = 1, nzone |
---|
| 189 | q(i,k,iq) = q(i,k,iq) / xtractot * xiiso !--- Bidouille pour convergence |
---|
| 190 | END DO |
---|
| 191 | END DO |
---|
| 192 | END DO |
---|
| 193 | !$OMP END DO NOWAIT |
---|
| 194 | END DO |
---|
| 195 | END DO |
---|
[2270] | 196 | |
---|
[4143] | 197 | END SUBROUTINE check_isotopes |
---|
[2270] | 198 | |
---|