source: LMDZ6/trunk/libf/dyn3d/check_isotopes.F90 @ 5185

Last change on this file since 5185 was 5183, checked in by dcugnet, 10 months ago
  • Remove INCA retro-compatibility with "traceur.def" (containing only water tracers but getting chemical species from an internal INCA routine).
  • The "trac_type" derived type internal to "readTracFiles_mod" is removed because a generic "keys_type" is enough: no explicit key ("%" operator) is needed, even %name.
  • The "trac_type" and "isot_type" derived types are now defined locally in "infotrac" and "infotrac_phy" (and more generally in each context: dynamic, lmdz dynamics, lmdz physics, etc.). The "readTracFiles_mod" module is now only used in these two routines:
    • few internal routines/variables (addPhase, delPhase, new2oldH2O, newHNO3, oldHNO3) are made available through "infotrac" and "infotrac_phy".
    • the "getKey" routine is only used in these two routines to define the explicit keys ("%" operator) of the local derived types "trac_type" and "isot_type". It could be in principle used outside this scope to get tracers parameters (read from "tracer.def") or isotopic parameters (read from "isotopes_params.def" - disabled for now).
  • The fortran parameters file "iso_params_mod.F90" is introduced so that "tnat" and "alpha_ideal" are defined in a single place but used in several. "ltnat1" is no longer hardcoded but defined with the *.def files parameter "tnat1"
  • Few minor changes:
    • use "infotrac_phy" instead of "infotrac" in calfis* because "tracers(:)%isAdvected" is defined in physics only.
    • "isotopes_mod" now ready for several isotopes classes (currently: only H2O)
    • isotopes class name (the name of the parent of the isotopes) is now %name and no longer %parent.
    • improvement of "getKey"
File size: 8.3 KB
Line 
1SUBROUTINE check_isotopes_seq(q, ip1jmp1, err_msg)
2   USE strings_mod, ONLY: maxlen, msg, strIdx, strStack, int2str, real2str
3   USE infotrac,    ONLY: nqtot, niso, nphas, isotope, isoCheck, iqIsoPha, isoSelect, &
4                          ntiso, iH2O, nzone, tracers, isoName,  itZonIso
5   USE iso_params_mod,  ONLY: tnat_H216O, tnat_H217O, tnat_H218O, tnat_HDO, tnat_HTO
6   USE ioipsl_getincom, ONLY: getin
7   IMPLICIT NONE
8   include "dimensions.h"
9   REAL,             INTENT(INOUT) :: q(ip1jmp1,llm,nqtot)
10   INTEGER,          INTENT(IN)    :: ip1jmp1
11   CHARACTER(LEN=*), INTENT(IN)    :: err_msg    !--- Error message to display
12   CHARACTER(LEN=maxlen) :: modname, msg1, nm(2)
13   INTEGER :: ixt, ipha, k, i, iq, iiso, izon, ieau, iqeau, iqpar
14   INTEGER, ALLOCATABLE ::   ix(:)
15   REAL,    ALLOCATABLE, SAVE :: tnat(:)
16   REAL    :: xtractot, xiiso, deltaD, q1, q2
17   REAL, PARAMETER :: borne     = 1e19,  &
18                      errmax    = 1e-8,  &       !--- Max. absolute error
19                      errmaxrel = 1e-3,  &       !--- Max. relative error
20                      qmin      = 1e-11, &
21                      deltaDmax =1000.0, &
22                      deltaDmin =-999.0, &
23                      ridicule  = 1e-12
24   INTEGER, SAVE :: iso_eau, iso_O17, iso_O18, iso_HDO, iso_HTO
25   LOGICAL, SAVE :: ltnat1, first=.TRUE.
26
27   modname='check_isotopes'
28   IF(.NOT.isoCheck)    RETURN                   !--- No need to check => finished
29   IF(isoSelect('H2O')) RETURN                   !--- No H2O isotopes group found
30   IF(niso == 0)        RETURN                   !--- No isotopes => finished
31   IF(first) THEN
32      ltnat1 = .TRUE.; CALL getin('tnateq1', ltnat1)
33      ALLOCATE(tnat(niso))
34      iso_eau = strIdx(isoName,'H216O'); IF(iso_eau /= 0) tnat(iso_eau) = tnat_H216O
35      iso_O17 = strIdx(isoName,'H217O'); IF(iso_O17 /= 0) tnat(iso_O17) = tnat_H217O
36      iso_O18 = strIdx(isoName,'H218O'); IF(iso_O18 /= 0) tnat(iso_O18) = tnat_H218O
37      iso_HDO = strIdx(isoName,'HDO');   IF(iso_HDO /= 0) tnat(iso_HDO) = tnat_HDO
38      iso_HTO = strIdx(isoName,'HTO');   IF(iso_HTO /= 0) tnat(iso_HTO) = tnat_HTO
39      IF(ltnat1) tnat(:) = 1.
40      first = .FALSE.
41   END IF
42   CALL msg('31: err_msg='//TRIM(err_msg), modname)
43
44   !--- CHECK FOR NaNs (FOR ALL THE ISOTOPES, INCLUDING GEOGRAPHIC TAGGING TRACERS)
45   modname = 'check_isotopes:iso_verif_noNaN'
46   DO ixt = 1, ntiso
47      DO ipha = 1, nphas
48         iq = iqIsoPha(ixt,ipha)
49         DO k = 1, llm
50            DO i = 1, ip1jmp1
51               IF(ABS(q(i,k,iq)) < borne) CYCLE
52               WRITE(msg1,'(s,"(",i0,",",i0,",",i0,") = ",ES12.4)')TRIM(isoName(ixt)),i,k,iq,q(i,k,iq)
53               CALL msg(msg1, modname)
54               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
55            END DO
56         END DO
57      END DO
58   END DO
59
60   !--- CHECK CONSERVATION (MAIN ISOTOPE AND PARENT CONCENTRATIONS MUST BE EQUAL)
61   modname = 'check_isotopes:iso_verif_egalite'
62   ixt = iso_eau
63   IF(ixt /= 0) THEN
64      DO ipha = 1, nphas
65         iq = iqIsoPha(ixt,ipha)
66         iqpar = tracers(iq)%iqParent
67         DO k = 1, llm
68            DO i = 1, ip1jmp1
69               q1 = q(i,k,iqpar)
70               q2 = q(i,k,iq)
71!--- IMPROVEMENT in case at least one isotope is not negligible compared to the main isotopic form.
72!    This would be probably required to sum from smallest to highest concentrations ; the corresponding
73!    indices vector can be computed once only (in the initializations stage), using mean concentrations.
74!              q2 = SUM(q(i,k,tracers(iqPar)%iqDesc), DIM=3)
75               IF(ABS(q1-q2) <= errmax .OR. ABS(q1-q2)/MAX(MAX(ABS(q1),ABS(q2)),1e-18) <= errmaxrel) THEN
76                  q(i,k,iq) = q1                 !--- Bidouille pour convergence
77!                 q(i,k,tracers(iqPar)%iqDesc) = q(i,k,tracers(iqPar)%iqDesc) * q1 / q2
78                  CYCLE
79               END IF
80               CALL msg('ixt, iq = '//TRIM(strStack(int2str([ixt,iq]))), modname)
81               msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
82               CALL msg(TRIM(tracers(iqpar)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
83               CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
84               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
85            END DO
86         END DO
87      END DO
88   END IF
89
90   !--- CHECK DELTA ANOMALIES
91   modname = 'check_isotopes:iso_verif_aberrant'
92   ix = [ iso_HDO  ,   iso_O18 ]
93   nm = ['deltaD  ', 'deltaO18']
94   DO iiso = 1, SIZE(ix)
95      ixt = ix(iiso)
96      IF(ixt  == 0) CYCLE
97      DO ipha = 1, nphas
98         iq = iqIsoPha(ixt,ipha)
99         iqpar = tracers(iq)%iqParent
100         DO k = 1, llm
101            DO i = 1, ip1jmp1
102               q1 = q(i,k,iqpar)
103               q2 = q(i,k,iq)
104!--- IMPROVEMENT in case at least one isotope is not negligible compared to the main isotopic form.
105!    This would be probably required to sum from smallest to highest concentrations ; the corresponding
106!    indices vector can be computed once only (in the initializations stage), using mean concentrations.
107!              q2 = SUM(q(i,k,tracers(iqPar)%iqDesc), DIM=3)
108               IF(q2 <= qmin) CYCLE
109               deltaD = (q2/q1/tnat(ixt)-1.)*1000.
110               IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
111               CALL msg('ixt, iq = '//TRIM(strStack(int2str([ixt,iq]))), modname)
112               msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
113               CALL msg(TRIM(tracers(iqpar)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
114               CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
115               CALL msg(TRIM(nm(iiso))//TRIM(real2str(deltaD)), modname)
116               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
117            END DO
118         END DO
119      END DO
120   END DO
121
122   IF(nzone == 0) RETURN
123
124   !--- CHECK FOR TAGGING TRACERS DELTAD ANOMALIES
125   modname = 'check_isotopes:iso_verif_aberrant'
126   IF(iso_eau /= 0 .AND. iso_HDO /= 0) THEN
127      DO izon = 1, nzone
128         ixt  = itZonIso(izon, iso_HDO)
129         ieau = itZonIso(izon, iso_eau)
130         DO ipha = 1, nphas
131            iq    = iqIsoPha(ixt,  ipha)
132            iqeau = iqIsoPha(ieau, ipha)
133            DO k = 1, llm
134               DO i = 1, ip1jmp1
135                  q1 = q(i,k,iqeau)
136                  q2 = q(i,k,iq)
137                  IF(q2<=qmin) CYCLE
138                  deltaD = (q2/q1/tnat(iso_HDO)-1.)*1000.
139                  IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
140                  CALL msg('izon, ipha = '//TRIM(strStack(int2str([izon, ipha]))), modname)
141                  CALL msg( 'ixt, ieau = '//TRIM(strStack(int2str([ ixt, ieau]))), modname)
142                  msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
143                  CALL msg(TRIM(tracers(iqeau)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
144                  CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
145                  CALL msg('deltaD = '//TRIM(real2str(deltaD)), modname)
146                  CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
147               END DO
148            END DO
149         END DO
150      END DO
151   END IF
152
153   !--- CHECK FOR TAGGING TRACERS CONSERVATION (PARENT AND TAGGING TRACERS SUM OVER ALL REGIONS MUST BE EQUAL)
154   DO iiso = 1, niso
155      DO ipha = 1, nphas
156         iq = iqIsoPha(iiso, ipha)
157         DO k = 1, llm
158            DO i = 1, ip1jmp1
159               xiiso = q(i,k,iq)
160               xtractot = SUM(q(i, k, iqIsoPha(itZonIso(1:nzone,iiso), ipha)))
161               IF(ABS(xtractot-xiiso) > errmax .AND. ABS(xtractot-xiiso)/MAX(MAX(ABS(xtractot),ABS(xiiso)),1e-18) > errmaxrel) THEN
162                  CALL msg('Error in iso_verif_aberrant trac: '//TRIM(err_msg))
163                  CALL msg('iiso, ipha = '//TRIM(strStack(int2str([iiso, ipha]))), modname)
164                  CALL msg('q('//TRIM(strStack(int2str([i,k])))//',:) = '//TRIM(strStack(real2str(q(i,k,:)))), modname)
165                  CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
166               END IF
167               IF(ABS(xtractot) <= ridicule) CYCLE
168               DO izon = 1, nzone
169                  q(i,k,iq) = q(i,k,iq) / xtractot * xiiso !--- Bidouille pour convergence
170               END DO
171            END DO
172         END DO
173      END DO
174   END DO
175
176END SUBROUTINE check_isotopes_seq
177
Note: See TracBrowser for help on using the repository browser.