source: LMDZ6/branches/DYNAMICO-conv/libf/dyn3d/check_isotopes.F @ 4976

Last change on this file since 4976 was 2281, checked in by crisi, 10 years ago

Camille Risi: corrections of bugs for the isotopic part

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