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

Last change on this file since 5244 was 5214, checked in by dcugnet, 4 weeks ago

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