source: LMDZ6/branches/LMDZ-tracers/libf/dyn3d/check_isotopes.F90 @ 3891

Last change on this file since 3891 was 3891, checked in by dcugnet, 3 years ago
  • Bugs corrections:
    • sequential gcm fixed
    • parallel gcm compilation fixed ; to be tested
  • Some generic operations moved from infotrac to readTracFile
  • Fixed algebrical reduction routine, used in the isotopes parameters file.
  • Additional component "comp" in the tracers descriptor derived type "tra",

specifying the model component name(s) (cf. tracers sections) it belongs.

  • isotopes class selection tool fixed.
File size: 5.7 KB
Line 
1SUBROUTINE check_isotopes_seq(q, ip1jmp1, err_msg)
2  USE strings_mod, ONLY: strIdx, msg, modname, prt_level
3  USE infotrac,    ONLY: isotope, isoSelect, iH2O, isoCheck, isoName, nqtot, niso, nitr, nzon, npha, iTraPha, iZonIso, tnat
4  IMPLICIT NONE
5  include "dimensions.h"
6  REAL,             INTENT(INOUT) :: q(ip1jmp1,llm,nqtot)
7  INTEGER,          INTENT(IN)    :: ip1jmp1
8  CHARACTER(LEN=*), INTENT(IN)    :: err_msg     !--- Error message to display
9  CHARACTER(LEN=256) :: msg1
10  INTEGER :: ixt, ipha, k, i, iq, iiso, izon, ieau, iqeau
11  INTEGER, ALLOCATABLE :: ix(:)
12  REAL    :: xtractot, xiiso, deltaD, q1, q2
13  REAL, PARAMETER :: borne     = 1e19,  &
14                     errmax    = 1e-8,  &        !--- Max. absolute error
15                     errmaxrel = 1e-8,  &        !--- Max. relative error
16                     qmin      = 1e-11, &
17                     deltaDmax = 200.0, &
18                     deltaDmin =-999.9, &
19                     ridicule  = 1e-12
20  INTEGER, SAVE :: ixH2O, ixHDO, ixO18
21  LOGICAL, SAVE :: first=.TRUE.
22
23  modname = 'check_isotopes'
24  IF(first) THEN
25    iH2O = -1
26    IF(isoSelect('H2O')) RETURN
27    ixH2O = strIdx(isoName,'H2[16]O')
28    ixHDO = strIdx(isoName,'H[2]HO')
29    ixO18 = strIdx(isoName,'H2[18]O')
30    first = .FALSE.
31  ELSE
32    IF(iH2O == -1)      RETURN
33    IF(isoSelect(iH2O)) RETURN
34  END IF
35  IF(.NOT.isoCheck .OR. niso == 0) RETURN        !--- No need to check or no isotopes => finished
36
37  !--- CHECK FOR NaNs (FOR ALL THE ISOTOPES, INCLUDING GEOGRAPHIC TAGGING TRACERS)
38  DO ixt = 1, nitr
39    DO ipha = 1, npha
40      iq = iTraPha(ixt,ipha)
41      DO k = 1, llm
42        DO i = 1, ip1jmp1
43          IF(ABS(q(i,k,iq))<=borne) CYCLE
44          WRITE(msg1,'(s,"(",i0,",",i0,",",i0,") = ",ES12.4)')isoName(ixt),i,k,iq,q(i,k,iq); CALL msg(msg1)
45          CALL abort_gcm(modname, 'Error in iso_verif_noNaN: '//TRIM(err_msg), 1)
46          STOP
47        END DO
48      END DO
49    END DO
50  END DO
51
52  !--- CHECK CONSERVATION (MAIN ISOTOPE AND PARENT CONCENTRATIONS MUST BE EQUAL)
53  ixt = ixH2O
54  IF(ixt /= 0) THEN
55    DO ipha = 1, npha
56      iq = iTraPha(ixt,ipha)
57      DO k = 1, llm
58        DO i = 1, ip1jmp1
59          q1 = q(i,k,ipha); q2 = q(i,k,iq)
60          IF(ABS(q1-q2) <= errmax .OR. ABS(q1-q2)/MAX(MAX(ABS(q1),ABS(q2)),1e-18) < errmaxrel) CYCLE
61          WRITE(msg1,'("ixt = ",i0)')ixt;                                      CALL msg(msg1)
62          WRITE(msg1,'("q(",i0,",",i0,",iq=",i0,") = ",ES12.4)')i, k, iq, q2;  CALL msg(msg1)
63          WRITE(msg1,'("q(",i0,",",i0,",ipha=",i0,") = ",ES12.4)')i,k,ipha,q1; CALL msg(msg1)
64          CALL abort_gcm(modname, 'Error in iso_verif_egalite: '//TRIM(err_msg), 1)
65          q(i,k,iq) = q(i,k,ipha)                !--- Bidouille pour convergence
66        END DO
67      END DO
68    END DO
69  END IF
70
71  !--- CHECK DELTA ANOMALIES
72  ix = [ixHDO, ixO18]
73  DO iiso = 1, SIZE(ix)
74    ixt = ix(iiso)
75    IF(ixt  == 0) CYCLE
76    DO ipha = 1, npha
77      iq = iTraPha(ixt,ipha)
78      DO k = 1, llm
79        DO i = 1, ip1jmp1
80          q1 = q(i,k,ipha); q2 = q(i,k,iq)
81          IF(q2 <= qmin) CYCLE
82          deltaD = (q2/q1/tnat(ixt)-1)*1000
83          IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
84          WRITE(msg1,'("ixt = ",i0)')ixt;                                     CALL msg(msg1)
85          WRITE(msg1,'("q(",i0,",",i0,",iq=",i0,") = ",ES12.4)')i, k, iq, q2; CALL msg(msg1)
86          WRITE(msg1,'("q=",ES12.4)')q(i,k,:);                                CALL msg(msg1)
87          WRITE(msg1,'("deltaD=",ES12.4)')deltaD;                             CALL msg(msg1)
88          CALL abort_gcm(modname, 'Error in iso_verif_aberrant: '//TRIM(err_msg), 1)
89        END DO
90      END DO
91    END DO
92  END DO
93
94  !--- CHECK FOR TAGGING TRACERS DELTAD ANOMALIES
95  IF(nitr == 0) RETURN
96  IF(ixH2O /= 0 .AND. ixHDO /= 0) THEN
97    DO izon = 1, nzon
98      ixt  = iZonIso(izon, ixHDO)
99      ieau = iZonIso(izon, ixH2O)
100      DO ipha = 1, npha
101        iq    = iTraPha(ixt,  ipha)
102        iqeau = iTraPha(ieau, ipha)
103        DO k = 1, llm
104          DO i = 1, ip1jmp1
105            IF(q(i,k,iq)<=qmin) CYCLE
106            deltaD = (q(i,k,iq)/q(i,k,iqeau)/tnat(ixHDO)-1)*1000
107            IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
108            WRITE(msg1,'("izon, ipha =",2i0)')izon, ipha;                              CALL msg(msg1)
109            WRITE(msg1,'( "ixt, ieau =",2i0)') ixt, ieau;                              CALL msg(msg1)
110            WRITE(msg1,'("q(",i0,",",i0,",iq=",i0,") = ",ES12.4)')i, k, iq, q(i,k,iq); CALL msg(msg1)
111            WRITE(msg1,'("deltaD=",ES12.4)')deltaD;                                    CALL msg(msg1)
112            CALL abort_gcm(modname, 'Error in iso_verif_aberrant trac: '//TRIM(err_msg), 1)
113          END DO
114        END DO
115      END DO
116    END DO
117  END IF
118
119  !--- CHECK FOR TAGGING TRACERS CONSERVATION (PARENT AND TAGGING TRACERS SUM OVER ALL REGIONS MUST BE EQUAL)
120  DO iiso = 1, niso
121    DO ipha = 1, npha
122      iq = iTraPha(iiso, ipha)
123      DO k = 1, llm
124        DO i = 1, ip1jmp1
125          xiiso = q(i,k,iq)
126          xtractot = SUM(q(i, k, iTraPha(iZonIso(1:nzon,iiso), ipha)))
127          IF(ABS(xtractot-xiiso) > errmax .AND. ABS(xtractot-xiiso)/MAX(MAX(ABS(xtractot),ABS(xiiso)),1e-18) > errmaxrel) THEN
128            CALL msg('Error in iso_verif_aberrant trac: '//TRIM(err_msg))
129            WRITE(msg1,'("iiso, ipha =",2i0)')iiso, ipha;              CALL msg(msg1)
130            WRITE(msg1,'("i, k =",2i0)')i, k;                          CALL msg(msg1)
131            WRITE(msg1,'("q(",i0,",",i0,":) = ",ES12.4)')i,k,q(i,k,:); CALL msg(msg1)
132            STOP
133          END IF
134          IF(ABS(xtractot) <= ridicule) CYCLE
135          DO izon = 1, nzon
136            ixt = iZonIso(izon, iiso)
137            q(i,k,iq) = q(i,k,iq) / xtractot * xiiso
138          END DO
139        END DO
140      END DO
141    END DO
142  END DO
143
144END SUBROUTINE check_isotopes_seq
145
Note: See TracBrowser for help on using the repository browser.