1 | subroutine check_isotopes_seq(q,ip1jmp1,err_msg) |
---|
2 | USE infotrac, ONLY: nqtot, nqo, niso, ntraciso, ntraceurs_zone, |
---|
3 | & ok_isotopes, ok_isotrac, use_iso, |
---|
4 | & iqiso, index_trac,indnum_fn_num, tnat |
---|
5 | implicit none |
---|
6 | |
---|
7 | #include "dimensions.h" |
---|
8 | |
---|
9 | ! inputs |
---|
10 | integer ip1jmp1 |
---|
11 | real q(ip1jmp1,llm,nqtot) |
---|
12 | character*(*) err_msg ! message d''erreur à afficher |
---|
13 | |
---|
14 | ! locals |
---|
15 | integer ixt,phase,k,i,iq,iiso,izone,ieau,iqeau |
---|
16 | real xtractot,xiiso |
---|
17 | real borne |
---|
18 | real qmin |
---|
19 | real errmax ! erreur maximale en absolu. |
---|
20 | real errmaxrel ! erreur maximale en relatif autorisée |
---|
21 | real deltaDmax,deltaDmin |
---|
22 | real ridicule |
---|
23 | parameter (borne=1e19) |
---|
24 | parameter (errmax=1e-8) |
---|
25 | parameter (errmaxrel=1e-3) |
---|
26 | parameter (qmin=1e-11) |
---|
27 | parameter (deltaDmax=1000.0,deltaDmin=-999.0) |
---|
28 | parameter (ridicule=1e-12) |
---|
29 | real deltaD |
---|
30 | |
---|
31 | if (ok_isotopes) then |
---|
32 | |
---|
33 | write(*,*) 'check_isotopes 31: err_msg=',err_msg |
---|
34 | ! verifier que rien n'est NaN |
---|
35 | do ixt=1,ntraciso |
---|
36 | do phase=1,nqo |
---|
37 | iq=iqiso(ixt,phase) |
---|
38 | do k=1,llm |
---|
39 | DO i = 1,ip1jmp1 |
---|
40 | if ((q(i,k,iq).gt.-borne).and. |
---|
41 | : (q(i,k,iq).lt.borne)) then |
---|
42 | else !if ((x(ixt,i,j).gt.-borne).and. |
---|
43 | write(*,*) 'erreur detectee par iso_verif_noNaN:' |
---|
44 | write(*,*) err_msg |
---|
45 | write(*,*) 'q,i,k,iq=',q(i,k,iq),i,k,iq |
---|
46 | write(*,*) 'borne=',borne |
---|
47 | call abort_gcm('check_isotopes_loc','plantage iso',0) |
---|
48 | endif !if ((x(ixt,i,j).gt.-borne).and. |
---|
49 | enddo !DO i = 1,ip1jmp1 |
---|
50 | enddo !do k=1,llm |
---|
51 | enddo !do phase=1,nqo |
---|
52 | enddo !do ixt=1,ntraciso |
---|
53 | |
---|
54 | !write(*,*) 'check_isotopes 52' |
---|
55 | ! verifier que l'eau normale est OK |
---|
56 | if (use_iso(1)) then |
---|
57 | ixt=indnum_fn_num(1) |
---|
58 | do phase=1,nqo |
---|
59 | iq=iqiso(ixt,phase) |
---|
60 | do k=1,llm |
---|
61 | DO i = 1,ip1jmp1 |
---|
62 | if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and. |
---|
63 | : (abs((q(i,k,phase)-q(i,k,iq))/ |
---|
64 | : max(max(abs(q(i,k,phase)),abs(q(i,k,iq))),1e-18)) |
---|
65 | : .gt.errmaxrel)) then |
---|
66 | write(*,*) 'erreur detectee par iso_verif_egalite:' |
---|
67 | write(*,*) err_msg |
---|
68 | write(*,*) 'ixt,phase=',ixt,phase |
---|
69 | write(*,*) 'q,iq,i,k=',q(i,k,iq),iq,i,k |
---|
70 | write(*,*) 'q(i,k,phase)=',q(i,k,phase) |
---|
71 | call abort_gcm('check_isotopes_loc','plantage iso',0) |
---|
72 | endif !if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and. |
---|
73 | ! bidouille pour éviter divergence: |
---|
74 | q(i,k,iq)= q(i,k,phase) |
---|
75 | enddo ! DO i = 1,ip1jmp1 |
---|
76 | enddo !do k=1,llm |
---|
77 | enddo ! do phase=1,nqo |
---|
78 | endif !if (use_iso(1)) then |
---|
79 | |
---|
80 | !write(*,*) 'check_isotopes 78' |
---|
81 | ! verifier que HDO est raisonable |
---|
82 | if (use_iso(2)) then |
---|
83 | ixt=indnum_fn_num(2) |
---|
84 | do phase=1,nqo |
---|
85 | iq=iqiso(ixt,phase) |
---|
86 | do k=1,llm |
---|
87 | DO i = 1,ip1jmp1 |
---|
88 | if (q(i,k,iq).gt.qmin) then |
---|
89 | deltaD=(q(i,k,iq)/q(i,k,phase)/tnat(2)-1)*1000 |
---|
90 | if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then |
---|
91 | write(*,*) 'erreur detectee par iso_verif_aberrant:' |
---|
92 | write(*,*) err_msg |
---|
93 | write(*,*) 'ixt,phase=',ixt,phase |
---|
94 | write(*,*) 'q,iq,i,k,=',q(i,k,iq),iq,i,k |
---|
95 | write(*,*) 'q=',q(i,k,:) |
---|
96 | write(*,*) 'deltaD=',deltaD |
---|
97 | call abort_gcm('check_isotopes_loc','plantage iso',0) |
---|
98 | endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then |
---|
99 | endif !if (q(i,k,iq).gt.qmin) then |
---|
100 | enddo !DO i = 1,ip1jmp1 |
---|
101 | enddo !do k=1,llm |
---|
102 | enddo ! do phase=1,nqo |
---|
103 | endif !if (use_iso(2)) then |
---|
104 | |
---|
105 | !write(*,*) 'check_isotopes 103' |
---|
106 | ! verifier que O18 est raisonable |
---|
107 | if (use_iso(3)) then |
---|
108 | ixt=indnum_fn_num(3) |
---|
109 | do phase=1,nqo |
---|
110 | iq=iqiso(ixt,phase) |
---|
111 | do k=1,llm |
---|
112 | DO i = 1,ip1jmp1 |
---|
113 | if (q(i,k,iq).gt.qmin) then |
---|
114 | deltaD=(q(i,k,iq)/q(i,k,phase)/tnat(3)-1)*1000 |
---|
115 | if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then |
---|
116 | write(*,*) 'erreur detectee iso_verif_aberrant O18:' |
---|
117 | write(*,*) err_msg |
---|
118 | write(*,*) 'ixt,phase=',ixt,phase |
---|
119 | write(*,*) 'q,iq,i,k,=',q(i,k,phase),iq,i,k |
---|
120 | write(*,*) 'xt=',q(i,k,:) |
---|
121 | write(*,*) 'deltaO18=',deltaD |
---|
122 | call abort_gcm('check_isotopes_loc','plantage iso',0) |
---|
123 | endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then |
---|
124 | endif !if (q(i,k,iq).gt.qmin) then |
---|
125 | enddo !DO i = 1,ip1jmp1 |
---|
126 | enddo !do k=1,llm |
---|
127 | enddo ! do phase=1,nqo |
---|
128 | endif !if (use_iso(2)) then |
---|
129 | |
---|
130 | |
---|
131 | !write(*,*) 'check_isotopes 129' |
---|
132 | if (ok_isotrac) then |
---|
133 | |
---|
134 | if (use_iso(2).and.use_iso(1)) then |
---|
135 | do izone=1,ntraceurs_zone |
---|
136 | ixt=index_trac(izone,indnum_fn_num(2)) |
---|
137 | ieau=index_trac(izone,indnum_fn_num(1)) |
---|
138 | do phase=1,nqo |
---|
139 | iq=iqiso(ixt,phase) |
---|
140 | iqeau=iqiso(ieau,phase) |
---|
141 | do k=1,llm |
---|
142 | DO i = 1,ip1jmp1 |
---|
143 | if (q(i,k,iq).gt.qmin) then |
---|
144 | deltaD=(q(i,k,iq)/q(i,k,iqeau)/tnat(2)-1)*1000 |
---|
145 | if ((deltaD.gt.deltaDmax).or. |
---|
146 | & (deltaD.lt.deltaDmin)) then |
---|
147 | write(*,*) 'erreur dans iso_verif_aberrant trac:' |
---|
148 | write(*,*) err_msg |
---|
149 | write(*,*) 'izone,phase=',izone,phase |
---|
150 | write(*,*) 'ixt,ieau=',ixt,ieau |
---|
151 | write(*,*) 'q,iq,i,k,=',q(i,k,iq),iq,i,k |
---|
152 | write(*,*) 'deltaD=',deltaD |
---|
153 | call abort_gcm('check_isotopes_loc','plantage iso',0) |
---|
154 | endif !if ((deltaD.gt.deltaDmax).or. |
---|
155 | endif !if (q(i,k,iq).gt.qmin) then |
---|
156 | enddo !DO i = 1,ip1jmp1 |
---|
157 | enddo ! do k=1,llm |
---|
158 | enddo ! do phase=1,nqo |
---|
159 | enddo !do izone=1,ntraceurs_zone |
---|
160 | endif !if (use_iso(2).and.use_iso(1)) then |
---|
161 | |
---|
162 | do iiso=1,niso |
---|
163 | do phase=1,nqo |
---|
164 | iq=iqiso(iiso,phase) |
---|
165 | do k=1,llm |
---|
166 | DO i = 1,ip1jmp1 |
---|
167 | xtractot=0.0 |
---|
168 | xiiso=q(i,k,iq) |
---|
169 | do izone=1,ntraceurs_zone |
---|
170 | iq=iqiso(index_trac(izone,iiso),phase) |
---|
171 | xtractot=xtractot+ q(i,k,iq) |
---|
172 | enddo !do izone=1,ntraceurs_zone |
---|
173 | if ((abs(xtractot-xiiso).gt.errmax).and. |
---|
174 | : (abs(xtractot-xiiso)/ |
---|
175 | : max(max(abs(xtractot),abs(xiiso)),1e-18) |
---|
176 | : .gt.errmaxrel)) then |
---|
177 | write(*,*) 'erreur detectee par iso_verif_traceurs:' |
---|
178 | write(*,*) err_msg |
---|
179 | write(*,*) 'iiso,phase=',iiso,phase |
---|
180 | write(*,*) 'i,k,=',i,k |
---|
181 | write(*,*) 'q(i,k,:)=',q(i,k,:) |
---|
182 | call abort_gcm('check_isotopes_loc','plantage iso',0) |
---|
183 | endif !if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and. |
---|
184 | |
---|
185 | ! bidouille pour éviter divergence: |
---|
186 | if (abs(xtractot).gt.ridicule) then |
---|
187 | do izone=1,ntraceurs_zone |
---|
188 | ixt=index_trac(izone,iiso) |
---|
189 | q(i,k,iq)=q(i,k,iq)/xtractot*xiiso |
---|
190 | enddo !do izone=1,ntraceurs_zone |
---|
191 | endif !if ((abs(xtractot).gt.ridicule) then |
---|
192 | enddo !DO i = 1,ip1jmp1 |
---|
193 | enddo !do k=1,llm |
---|
194 | enddo !do phase=1,nqo |
---|
195 | enddo !do iiso=1,niso |
---|
196 | |
---|
197 | endif !if (ok_isotrac) then |
---|
198 | |
---|
199 | endif ! if (ok_isotopes) |
---|
200 | !write(*,*) 'check_isotopes 198' |
---|
201 | |
---|
202 | end |
---|
203 | |
---|