Context Navigation

check_isotopes.F90 @ 4608

Last change on this file since 4608 was 4399, checked in by dcugnet, 23 months ago
Change the names of the water isotopes: no more "and?" separators for atomic numbers.
File size: 7.9 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format