source: LMDZ6/branches/LMDZ-tracers/libf/dyn3dmem/check_isotopes_loc.F90 @ 3960

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