| 1 | subroutine check_isotopes_seq(q,ip1jmp1,err_msg) |
|---|
| 2 | USE infotrac, ONLY: nqtot, nqo, niso, ntraciso, nzone, |
|---|
| 3 | & 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 (niso > 0) 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 (nzone > 0) then |
|---|
| 133 | |
|---|
| 134 | if (use_iso(2).and.use_iso(1)) then |
|---|
| 135 | do izone=1,nzone |
|---|
| 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,nzone |
|---|
| 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,nzone |
|---|
| 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,nzone |
|---|
| 188 | ixt=index_trac(izone,iiso) |
|---|
| 189 | q(i,k,iq)=q(i,k,iq)/xtractot*xiiso |
|---|
| 190 | enddo !do izone=1,nzone |
|---|
| 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 (nzone > 0) |
|---|
| 198 | |
|---|
| 199 | endif ! if (niso > 0) |
|---|
| 200 | !write(*,*) 'check_isotopes 198' |
|---|
| 201 | |
|---|
| 202 | end |
|---|
| 203 | |
|---|
