SUBROUTINE check_isotopes_seq(q, ip1jmp1, err_msg) USE strings_mod, ONLY: maxlen, msg, strIdx, strStack, num2str USE infotrac, ONLY: nqtot, niso, nphas, isotope, isoCheck, iqIsoPha, isoSelect, & ntiso, iH2O, nzone, tracers, isoName, itZonIso USE iso_params_mod, ONLY: tnat_H216O, tnat_H217O, tnat_H218O, tnat_HDO, tnat_HTO USE ioipsl, ONLY: getin USE dimensions_mod, ONLY: iim, jjm, llm, ndm IMPLICIT NONE REAL, INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) INTEGER, INTENT(IN) :: ip1jmp1 CHARACTER(LEN=*), INTENT(IN) :: err_msg !--- Error message to display CHARACTER(LEN=maxlen) :: modname, msg1, nm(2) INTEGER :: ixt, ipha, k, i, iq, iiso, izon, ieau, iqeau, iqpar INTEGER, ALLOCATABLE :: ix(:) REAL, ALLOCATABLE, SAVE :: tnat(:) REAL :: xtractot, xiiso, deltaD, q1, q2 REAL, PARAMETER :: borne = 1e19, & errmax = 1e-8, & !--- Max. absolute error errmaxrel = 1e-3, & !--- Max. relative error qmin = 1e-11, & deltaDmax =1000.0, & deltaDmin =-999.0, & ridicule = 1e-12 INTEGER, SAVE :: iso_eau, iso_O17, iso_O18, iso_HDO, iso_HTO LOGICAL :: ltnat1 LOGICAL, SAVE :: first=.TRUE. modname='check_isotopes' IF(.NOT.isoCheck) RETURN !--- No need to check => finished IF(isoSelect('H2O')) RETURN !--- No H2O isotopes group found IF(niso == 0) RETURN !--- No isotopes => finished IF(first) THEN ltnat1 = .TRUE.; CALL getin('tnateq1', ltnat1) ALLOCATE(tnat(niso)) iso_eau = strIdx(isoName,'H216O'); IF(iso_eau /= 0) tnat(iso_eau) = tnat_H216O iso_O17 = strIdx(isoName,'H217O'); IF(iso_O17 /= 0) tnat(iso_O17) = tnat_H217O iso_O18 = strIdx(isoName,'H218O'); IF(iso_O18 /= 0) tnat(iso_O18) = tnat_H218O iso_HDO = strIdx(isoName,'HDO'); IF(iso_HDO /= 0) tnat(iso_HDO) = tnat_HDO iso_HTO = strIdx(isoName,'HTO'); IF(iso_HTO /= 0) tnat(iso_HTO) = tnat_HTO IF(ltnat1) tnat(:) = 1. first = .FALSE. END IF CALL msg('31: err_msg='//TRIM(err_msg), modname) !--- CHECK FOR NaNs (FOR ALL THE ISOTOPES, INCLUDING GEOGRAPHIC TAGGING TRACERS) modname = 'check_isotopes:iso_verif_noNaN' DO ixt = 1, ntiso DO ipha = 1, nphas iq = iqIsoPha(ixt,ipha) DO k = 1, llm DO i = 1, ip1jmp1 IF(ABS(q(i,k,iq)) < borne) CYCLE WRITE(msg1,'(s,"(",i0,",",i0,",",i0,") = ",ES12.4)')TRIM(isoName(ixt)),i,k,iq,q(i,k,iq) CALL msg(msg1, modname) CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1) END DO END DO END DO END DO !--- CHECK CONSERVATION (MAIN ISOTOPE AND PARENT CONCENTRATIONS MUST BE EQUAL) modname = 'check_isotopes:iso_verif_egalite' ixt = iso_eau IF(ixt /= 0) THEN DO ipha = 1, nphas iq = iqIsoPha(ixt,ipha) iqpar = tracers(iq)%iqParent DO k = 1, llm DO i = 1, ip1jmp1 q1 = q(i,k,iqpar) q2 = q(i,k,iq) !--- IMPROVEMENT in case at least one isotope is not negligible compared to the main isotopic form. ! This would be probably required to sum from smallest to highest concentrations ; the corresponding ! indices vector can be computed once only (in the initializations stage), using mean concentrations. ! q2 = SUM(q(i,k,tracers(iqPar)%iqDesc), DIM=3) IF(ABS(q1-q2) <= errmax .OR. ABS(q1-q2)/MAX(MAX(ABS(q1),ABS(q2)),1e-18) <= errmaxrel) THEN q(i,k,iq) = q1 !--- Bidouille pour convergence ! q(i,k,tracers(iqPar)%iqDesc) = q(i,k,tracers(iqPar)%iqDesc) * q1 / q2 CYCLE END IF CALL msg('ixt, iq = '//TRIM(strStack(num2str([ixt,iq]))), modname) msg1 = '('//TRIM(strStack(num2str([i,k])))//')' CALL msg(TRIM(tracers(iqpar)%name)//TRIM(msg1)//' = '//TRIM(num2str(q1)), modname) CALL msg(TRIM(tracers(iq )%name)//TRIM(msg1)//' = '//TRIM(num2str(q2)), modname) CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1) END DO END DO END DO END IF !--- CHECK DELTA ANOMALIES modname = 'check_isotopes:iso_verif_aberrant' ix = [ iso_HDO , iso_O18 ] nm = ['deltaD ', 'deltaO18'] DO iiso = 1, SIZE(ix) ixt = ix(iiso) IF(ixt == 0) CYCLE DO ipha = 1, nphas iq = iqIsoPha(ixt,ipha) iqpar = tracers(iq)%iqParent DO k = 1, llm DO i = 1, ip1jmp1 q1 = q(i,k,iqpar) q2 = q(i,k,iq) !--- IMPROVEMENT in case at least one isotope is not negligible compared to the main isotopic form. ! This would be probably required to sum from smallest to highest concentrations ; the corresponding ! indices vector can be computed once only (in the initializations stage), using mean concentrations. ! q2 = SUM(q(i,k,tracers(iqPar)%iqDesc), DIM=3) IF(q2 <= qmin) CYCLE deltaD = (q2/q1/tnat(ixt)-1.)*1000. IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE CALL msg('ixt, iq = '//TRIM(strStack(num2str([ixt,iq]))), modname) msg1 = '('//TRIM(strStack(num2str([i,k])))//')' CALL msg(TRIM(tracers(iqpar)%name)//TRIM(msg1)//' = '//TRIM(num2str(q1)), modname) CALL msg(TRIM(tracers(iq )%name)//TRIM(msg1)//' = '//TRIM(num2str(q2)), modname) CALL msg(TRIM(nm(iiso))//TRIM(num2str(deltaD)), modname) CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1) END DO END DO END DO END DO IF(nzone == 0) RETURN !--- CHECK FOR TAGGING TRACERS DELTAD ANOMALIES modname = 'check_isotopes:iso_verif_aberrant' IF(iso_eau /= 0 .AND. iso_HDO /= 0) THEN DO izon = 1, nzone ixt = itZonIso(izon, iso_HDO) ieau = itZonIso(izon, iso_eau) DO ipha = 1, nphas iq = iqIsoPha(ixt, ipha) iqeau = iqIsoPha(ieau, ipha) DO k = 1, llm DO i = 1, ip1jmp1 q1 = q(i,k,iqeau) q2 = q(i,k,iq) IF(q2<=qmin) CYCLE deltaD = (q2/q1/tnat(iso_HDO)-1.)*1000. IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE CALL msg('izon, ipha = '//TRIM(strStack(num2str([izon, ipha]))), modname) CALL msg( 'ixt, ieau = '//TRIM(strStack(num2str([ ixt, ieau]))), modname) msg1 = '('//TRIM(strStack(num2str([i,k])))//')' CALL msg(TRIM(tracers(iqeau)%name)//TRIM(msg1)//' = '//TRIM(num2str(q1)), modname) CALL msg(TRIM(tracers(iq )%name)//TRIM(msg1)//' = '//TRIM(num2str(q2)), modname) CALL msg('deltaD = '//TRIM(num2str(deltaD)), modname) CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1) END DO END DO END DO END DO END IF !--- CHECK FOR TAGGING TRACERS CONSERVATION (PARENT AND TAGGING TRACERS SUM OVER ALL REGIONS MUST BE EQUAL) DO iiso = 1, niso DO ipha = 1, nphas iq = iqIsoPha(iiso, ipha) DO k = 1, llm DO i = 1, ip1jmp1 xiiso = q(i,k,iq) xtractot = SUM(q(i, k, iqIsoPha(itZonIso(1:nzone,iiso), ipha))) IF(ABS(xtractot-xiiso) > errmax .AND. ABS(xtractot-xiiso)/MAX(MAX(ABS(xtractot),ABS(xiiso)),1e-18) > errmaxrel) THEN CALL msg('Error in iso_verif_aberrant trac: '//TRIM(err_msg)) CALL msg('iiso, ipha = '//TRIM(strStack(num2str([iiso, ipha]))), modname) CALL msg('q('//TRIM(strStack(num2str([i,k])))//',:) = '//TRIM(strStack(num2str(q(i,k,:)))), modname) CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1) END IF IF(ABS(xtractot) <= ridicule) CYCLE DO izon = 1, nzone q(i,k,iq) = q(i,k,iq) / xtractot * xiiso !--- Bidouille pour convergence END DO END DO END DO END DO END DO END SUBROUTINE check_isotopes_seq