[2270] | 1 | subroutine check_isotopes(q,ijb,ije,err_msg) |
---|
| 2 | USE infotrac |
---|
| 3 | USE parallel_lmdz |
---|
| 4 | implicit none |
---|
| 5 | |
---|
| 6 | #include "dimensions.h" |
---|
| 7 | |
---|
| 8 | ! inputs |
---|
| 9 | integer ijb,ije ! peut être local et différent de ijb_u,ije_u, ex: dans qminimum |
---|
| 10 | real q(ijb_u:ije_u,llm,nqtot) |
---|
| 11 | character*(*) err_msg ! message d''erreur à afficher |
---|
| 12 | |
---|
| 13 | ! locals |
---|
[3331] | 14 | integer ixt,phase,k,i,iq,iiso,izone,ieau,iqeau,ixt2,iq2 |
---|
[2270] | 15 | real xtractot,xiiso |
---|
| 16 | real borne |
---|
| 17 | real qmin |
---|
| 18 | real errmax ! erreur maximale en absolu. |
---|
| 19 | real errmaxrel ! erreur maximale en relatif autorisée |
---|
| 20 | real deltaDmax,deltaDmin |
---|
| 21 | real ridicule |
---|
| 22 | parameter (borne=1e19) |
---|
| 23 | parameter (errmax=1e-8) |
---|
| 24 | parameter (errmaxrel=1e-3) |
---|
| 25 | parameter (qmin=1e-11) |
---|
[3331] | 26 | parameter (deltaDmax=1000.0,deltaDmin=-999.0) |
---|
[2270] | 27 | parameter (ridicule=1e-12) |
---|
| 28 | real deltaD |
---|
[3331] | 29 | real dexcessmax,dexcessmin |
---|
| 30 | parameter (dexcessmax=6000.0,dexcessmin=-100.0) |
---|
[2270] | 31 | |
---|
| 32 | if (ok_isotopes) then |
---|
| 33 | |
---|
[2281] | 34 | !write(*,*) 'check_isotopes 31: err_msg=',err_msg |
---|
[2270] | 35 | ! verifier que rien n'est NaN |
---|
| 36 | do ixt=1,ntraciso |
---|
| 37 | do phase=1,nqo |
---|
| 38 | iq=iqiso(ixt,phase) |
---|
[2281] | 39 | c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) |
---|
[2270] | 40 | do k=1,llm |
---|
| 41 | DO i = ijb,ije |
---|
| 42 | if ((q(i,k,iq).gt.-borne).and. |
---|
| 43 | : (q(i,k,iq).lt.borne)) then |
---|
| 44 | else !if ((x(ixt,i,j).gt.-borne).and. |
---|
| 45 | write(*,*) 'erreur detectee par iso_verif_noNaN:' |
---|
| 46 | write(*,*) err_msg |
---|
| 47 | write(*,*) 'q,i,k,iq=',q(i,k,iq),i,k,iq |
---|
| 48 | write(*,*) 'borne=',borne |
---|
| 49 | stop |
---|
| 50 | endif !if ((x(ixt,i,j).gt.-borne).and. |
---|
| 51 | enddo !DO i = ijb,ije |
---|
| 52 | enddo !do k=1,llm |
---|
[2281] | 53 | c$OMP END DO NOWAIT |
---|
[2270] | 54 | enddo !do phase=1,nqo |
---|
| 55 | enddo !do ixt=1,ntraciso |
---|
| 56 | |
---|
[2281] | 57 | !write(*,*) 'check_isotopes 52' |
---|
[2270] | 58 | ! verifier que l'eau normale est OK |
---|
| 59 | if (use_iso(1)) then |
---|
| 60 | ixt=indnum_fn_num(1) |
---|
| 61 | do phase=1,nqo |
---|
| 62 | iq=iqiso(ixt,phase) |
---|
[2281] | 63 | c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) |
---|
[2270] | 64 | do k=1,llm |
---|
| 65 | DO i = ijb,ije |
---|
| 66 | if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and. |
---|
| 67 | : (abs((q(i,k,phase)-q(i,k,iq))/ |
---|
| 68 | : max(max(abs(q(i,k,phase)),abs(q(i,k,iq))),1e-18)) |
---|
| 69 | : .gt.errmaxrel)) then |
---|
| 70 | write(*,*) 'erreur detectee par iso_verif_egalite:' |
---|
| 71 | write(*,*) err_msg |
---|
[2281] | 72 | write(*,*) 'ixt,phase,ijb=',ixt,phase,ijb |
---|
[2270] | 73 | write(*,*) 'q,iq,i,k=',q(i,k,iq),iq,i,k |
---|
| 74 | write(*,*) 'q(i,k,phase)=',q(i,k,phase) |
---|
| 75 | stop |
---|
| 76 | endif !if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and. |
---|
| 77 | ! bidouille pour éviter divergence: |
---|
| 78 | q(i,k,iq)= q(i,k,phase) |
---|
| 79 | enddo ! DO i = ijb,ije |
---|
| 80 | enddo !do k=1,llm |
---|
[2281] | 81 | c$OMP END DO NOWAIT |
---|
[2270] | 82 | enddo ! do phase=1,nqo |
---|
| 83 | endif !if (use_iso(1)) then |
---|
| 84 | |
---|
[2281] | 85 | !write(*,*) 'check_isotopes 78' |
---|
[2270] | 86 | ! verifier que HDO est raisonable |
---|
| 87 | if (use_iso(2)) then |
---|
| 88 | ixt=indnum_fn_num(2) |
---|
| 89 | do phase=1,nqo |
---|
| 90 | iq=iqiso(ixt,phase) |
---|
[2281] | 91 | c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) |
---|
[2270] | 92 | do k=1,llm |
---|
| 93 | DO i = ijb,ije |
---|
| 94 | if (q(i,k,iq).gt.qmin) then |
---|
| 95 | deltaD=(q(i,k,iq)/q(i,k,phase)/tnat(2)-1)*1000 |
---|
| 96 | if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then |
---|
| 97 | write(*,*) 'erreur detectee par iso_verif_aberrant:' |
---|
| 98 | write(*,*) err_msg |
---|
| 99 | write(*,*) 'ixt,phase=',ixt,phase |
---|
| 100 | write(*,*) 'q,iq,i,k,=',q(i,k,iq),iq,i,k |
---|
| 101 | write(*,*) 'q=',q(i,k,:) |
---|
| 102 | write(*,*) 'deltaD=',deltaD |
---|
| 103 | stop |
---|
| 104 | endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then |
---|
| 105 | endif !if (q(i,k,iq).gt.qmin) then |
---|
| 106 | enddo !DO i = ijb,ije |
---|
| 107 | enddo !do k=1,llm |
---|
[2281] | 108 | c$OMP END DO NOWAIT |
---|
[2270] | 109 | enddo ! do phase=1,nqo |
---|
| 110 | endif !if (use_iso(2)) then |
---|
| 111 | |
---|
[2281] | 112 | !write(*,*) 'check_isotopes 103' |
---|
[2270] | 113 | ! verifier que O18 est raisonable |
---|
| 114 | if (use_iso(3)) then |
---|
| 115 | ixt=indnum_fn_num(3) |
---|
| 116 | do phase=1,nqo |
---|
| 117 | iq=iqiso(ixt,phase) |
---|
[2281] | 118 | c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) |
---|
[2270] | 119 | do k=1,llm |
---|
| 120 | DO i = ijb,ije |
---|
| 121 | if (q(i,k,iq).gt.qmin) then |
---|
| 122 | deltaD=(q(i,k,iq)/q(i,k,phase)/tnat(3)-1)*1000 |
---|
| 123 | if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then |
---|
| 124 | write(*,*) 'erreur detectee iso_verif_aberrant O18:' |
---|
| 125 | write(*,*) err_msg |
---|
| 126 | write(*,*) 'ixt,phase=',ixt,phase |
---|
| 127 | write(*,*) 'q,iq,i,k,=',q(i,k,phase),iq,i,k |
---|
| 128 | write(*,*) 'xt=',q(i,k,:) |
---|
| 129 | write(*,*) 'deltaO18=',deltaD |
---|
| 130 | stop |
---|
| 131 | endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then |
---|
| 132 | endif !if (q(i,k,iq).gt.qmin) then |
---|
| 133 | enddo !DO i = ijb,ije |
---|
| 134 | enddo !do k=1,llm |
---|
[2281] | 135 | c$OMP END DO NOWAIT |
---|
[2270] | 136 | enddo ! do phase=1,nqo |
---|
| 137 | endif !if (use_iso(2)) then |
---|
| 138 | |
---|
[3331] | 139 | !write(*,*) 'check_isotopes 103' |
---|
| 140 | ! verifier que dexcess est raisonable |
---|
| 141 | if (use_iso(3).and.use_iso(2)) then |
---|
| 142 | ixt=indnum_fn_num(3) |
---|
| 143 | ixt2=indnum_fn_num(2) |
---|
| 144 | do phase=1,nqo |
---|
| 145 | iq=iqiso(ixt,phase) |
---|
| 146 | iq2=iqiso(ixt2,phase) |
---|
| 147 | c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) |
---|
| 148 | do k=1,llm |
---|
| 149 | DO i = ijb,ije |
---|
| 150 | if (q(i,k,iq).gt.qmin) then |
---|
| 151 | deltaD=(q(i,k,iq2)/q(i,k,phase)/tnat(2)-1)*1000 |
---|
| 152 | & -8.0*(q(i,k,iq)/q(i,k,phase)/tnat(3)-1)*1000 |
---|
| 153 | if ((deltaD.gt.dexcessmax).or.(deltaD.lt.dexcessmin)) then |
---|
| 154 | write(*,*) 'erreur detectee iso_verif_aberrant O18:' |
---|
| 155 | write(*,*) err_msg |
---|
| 156 | write(*,*) 'ixt,ixt2,phase=',ixt,ixt2,phase |
---|
| 157 | write(*,*) 'q,iq,i,k,=',q(i,k,phase),iq,i,k |
---|
| 158 | write(*,*) 'xt=',q(i,k,:) |
---|
| 159 | write(*,*) 'dexcess=',deltaD |
---|
| 160 | stop |
---|
| 161 | endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then |
---|
| 162 | endif !if (q(i,k,iq).gt.qmin) then |
---|
| 163 | enddo !DO i = ijb,ije |
---|
| 164 | enddo !do k=1,llm |
---|
| 165 | c$OMP END DO NOWAIT |
---|
| 166 | enddo ! do phase=1,nqo |
---|
| 167 | endif !if (use_iso(2)) then |
---|
[2270] | 168 | |
---|
[3331] | 169 | |
---|
| 170 | |
---|
[2281] | 171 | !write(*,*) 'check_isotopes 129' |
---|
[2270] | 172 | if (ok_isotrac) then |
---|
| 173 | |
---|
| 174 | if (use_iso(2).and.use_iso(1)) then |
---|
| 175 | do izone=1,ntraceurs_zone |
---|
| 176 | ixt=index_trac(izone,indnum_fn_num(2)) |
---|
| 177 | ieau=index_trac(izone,indnum_fn_num(1)) |
---|
| 178 | do phase=1,nqo |
---|
| 179 | iq=iqiso(ixt,phase) |
---|
| 180 | iqeau=iqiso(ieau,phase) |
---|
[2281] | 181 | c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) |
---|
[2270] | 182 | do k=1,llm |
---|
| 183 | DO i = ijb,ije |
---|
| 184 | if (q(i,k,iq).gt.qmin) then |
---|
| 185 | deltaD=(q(i,k,iq)/q(i,k,iqeau)/tnat(2)-1)*1000 |
---|
| 186 | if ((deltaD.gt.deltaDmax).or. |
---|
| 187 | & (deltaD.lt.deltaDmin)) then |
---|
| 188 | write(*,*) 'erreur dans iso_verif_aberrant trac:' |
---|
| 189 | write(*,*) err_msg |
---|
| 190 | write(*,*) 'izone,phase=',izone,phase |
---|
| 191 | write(*,*) 'ixt,ieau=',ixt,ieau |
---|
| 192 | write(*,*) 'q,iq,i,k,=',q(i,k,iq),iq,i,k |
---|
| 193 | write(*,*) 'deltaD=',deltaD |
---|
| 194 | stop |
---|
| 195 | endif !if ((deltaD.gt.deltaDmax).or. |
---|
| 196 | endif !if (q(i,k,iq).gt.qmin) then |
---|
| 197 | enddo !DO i = ijb,ije |
---|
| 198 | enddo ! do k=1,llm |
---|
[2281] | 199 | c$OMP END DO NOWAIT |
---|
[2270] | 200 | enddo ! do phase=1,nqo |
---|
| 201 | enddo !do izone=1,ntraceurs_zone |
---|
| 202 | endif !if (use_iso(2).and.use_iso(1)) then |
---|
| 203 | |
---|
| 204 | do iiso=1,niso |
---|
| 205 | do phase=1,nqo |
---|
| 206 | iq=iqiso(iiso,phase) |
---|
[2281] | 207 | c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) |
---|
[2270] | 208 | do k=1,llm |
---|
| 209 | DO i = ijb,ije |
---|
| 210 | xtractot=0.0 |
---|
| 211 | xiiso=q(i,k,iq) |
---|
| 212 | do izone=1,ntraceurs_zone |
---|
| 213 | iq=iqiso(index_trac(izone,iiso),phase) |
---|
| 214 | xtractot=xtractot+ q(i,k,iq) |
---|
| 215 | enddo !do izone=1,ntraceurs_zone |
---|
| 216 | if ((abs(xtractot-xiiso).gt.errmax).and. |
---|
| 217 | : (abs(xtractot-xiiso)/ |
---|
| 218 | : max(max(abs(xtractot),abs(xiiso)),1e-18) |
---|
| 219 | : .gt.errmaxrel)) then |
---|
| 220 | write(*,*) 'erreur detectee par iso_verif_traceurs:' |
---|
| 221 | write(*,*) err_msg |
---|
| 222 | write(*,*) 'iiso,phase=',iiso,phase |
---|
| 223 | write(*,*) 'i,k,=',i,k |
---|
| 224 | write(*,*) 'q(i,k,:)=',q(i,k,:) |
---|
| 225 | stop |
---|
| 226 | endif !if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and. |
---|
| 227 | |
---|
| 228 | ! bidouille pour éviter divergence: |
---|
| 229 | if (abs(xtractot).gt.ridicule) then |
---|
| 230 | do izone=1,ntraceurs_zone |
---|
| 231 | ixt=index_trac(izone,iiso) |
---|
| 232 | q(i,k,iq)=q(i,k,iq)/xtractot*xiiso |
---|
| 233 | enddo !do izone=1,ntraceurs_zone |
---|
| 234 | endif !if ((abs(xtractot).gt.ridicule) then |
---|
| 235 | enddo !DO i = ijb,ije |
---|
| 236 | enddo !do k=1,llm |
---|
[2281] | 237 | c$OMP END DO NOWAIT |
---|
[2270] | 238 | enddo !do phase=1,nqo |
---|
| 239 | enddo !do iiso=1,niso |
---|
| 240 | |
---|
| 241 | endif !if (ok_isotrac) then |
---|
| 242 | |
---|
| 243 | endif ! if (ok_isotopes) |
---|
[2281] | 244 | !write(*,*) 'check_isotopes 198' |
---|
[2270] | 245 | |
---|
| 246 | end |
---|
| 247 | |
---|
| 248 | |
---|