source: LMDZ6/branches/LMDZ-INCA-Dyn/libf/dyn3dmem/check_isotopes_loc.F @ 4484

Last change on this file since 4484 was 2281, checked in by crisi, 9 years ago

Camille Risi: corrections of bugs for the isotopic part
File size: 7.9 KB
Line 
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
14        integer ixt,phase,k,i,iq,iiso,izone,ieau,iqeau
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)
26        parameter (deltaDmax=200.0,deltaDmin=-999.9)
27        parameter (ridicule=1e-12)
28        real deltaD
29
30        if (ok_isotopes) then
31
32        !write(*,*) 'check_isotopes 31: err_msg=',err_msg
33        ! verifier que rien n'est NaN
34        do ixt=1,ntraciso
35          do phase=1,nqo
36            iq=iqiso(ixt,phase)
37c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
38            do k=1,llm
39              DO i = ijb,ije
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                  stop
48                endif  !if ((x(ixt,i,j).gt.-borne).and.
49              enddo !DO i = ijb,ije
50            enddo !do k=1,llm
51c$OMP END DO NOWAIT
52          enddo !do phase=1,nqo
53        enddo !do ixt=1,ntraciso
54
55        !write(*,*) 'check_isotopes 52'
56        ! verifier que l'eau normale est OK
57        if (use_iso(1)) then
58          ixt=indnum_fn_num(1)
59          do phase=1,nqo
60            iq=iqiso(ixt,phase)
61c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
62            do k=1,llm
63            DO i = ijb,ije 
64              if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and.
65     :          (abs((q(i,k,phase)-q(i,k,iq))/
66     :           max(max(abs(q(i,k,phase)),abs(q(i,k,iq))),1e-18))
67     :           .gt.errmaxrel)) then
68                  write(*,*) 'erreur detectee par iso_verif_egalite:'
69                  write(*,*) err_msg
70                  write(*,*) 'ixt,phase,ijb=',ixt,phase,ijb
71                  write(*,*) 'q,iq,i,k=',q(i,k,iq),iq,i,k
72                  write(*,*) 'q(i,k,phase)=',q(i,k,phase)
73                  stop
74              endif !if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and.
75              ! bidouille pour éviter divergence:
76              q(i,k,iq)= q(i,k,phase)
77            enddo ! DO i = ijb,ije
78            enddo !do k=1,llm
79c$OMP END DO NOWAIT
80          enddo ! do phase=1,nqo
81        endif !if (use_iso(1)) then
82       
83        !write(*,*) 'check_isotopes 78'
84        ! verifier que HDO est raisonable
85        if (use_iso(2)) then
86          ixt=indnum_fn_num(2)
87          do phase=1,nqo
88            iq=iqiso(ixt,phase)
89c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
90            do k=1,llm
91            DO i = ijb,ije
92            if (q(i,k,iq).gt.qmin) then
93             deltaD=(q(i,k,iq)/q(i,k,phase)/tnat(2)-1)*1000
94             if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
95                  write(*,*) 'erreur detectee par iso_verif_aberrant:'
96                  write(*,*) err_msg
97                  write(*,*) 'ixt,phase=',ixt,phase
98                  write(*,*) 'q,iq,i,k,=',q(i,k,iq),iq,i,k
99                  write(*,*) 'q=',q(i,k,:)
100                  write(*,*) 'deltaD=',deltaD
101                  stop
102             endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
103            endif !if (q(i,k,iq).gt.qmin) then
104            enddo !DO i = ijb,ije
105            enddo !do k=1,llm
106c$OMP END DO NOWAIT
107          enddo ! do phase=1,nqo
108        endif !if (use_iso(2)) then
109
110        !write(*,*) 'check_isotopes 103'
111        ! verifier que O18 est raisonable
112        if (use_iso(3)) then
113          ixt=indnum_fn_num(3)
114          do phase=1,nqo
115            iq=iqiso(ixt,phase)
116c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
117            do k=1,llm
118            DO i = ijb,ije
119            if (q(i,k,iq).gt.qmin) then
120             deltaD=(q(i,k,iq)/q(i,k,phase)/tnat(3)-1)*1000
121             if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
122                  write(*,*) 'erreur detectee iso_verif_aberrant O18:'
123                  write(*,*) err_msg
124                  write(*,*) 'ixt,phase=',ixt,phase
125                  write(*,*) 'q,iq,i,k,=',q(i,k,phase),iq,i,k
126                  write(*,*) 'xt=',q(i,k,:)
127                  write(*,*) 'deltaO18=',deltaD
128                  stop
129             endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
130            endif !if (q(i,k,iq).gt.qmin) then
131            enddo !DO i = ijb,ije
132            enddo !do k=1,llm
133c$OMP END DO NOWAIT
134          enddo ! do phase=1,nqo
135        endif !if (use_iso(2)) then
136
137
138        !write(*,*) 'check_isotopes 129'
139        if (ok_isotrac) then
140
141          if (use_iso(2).and.use_iso(1)) then
142            do izone=1,ntraceurs_zone
143             ixt=index_trac(izone,indnum_fn_num(2))
144             ieau=index_trac(izone,indnum_fn_num(1))
145             do phase=1,nqo
146               iq=iqiso(ixt,phase)
147               iqeau=iqiso(ieau,phase)
148c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
149               do k=1,llm
150                DO i = ijb,ije
151                if (q(i,k,iq).gt.qmin) then
152                 deltaD=(q(i,k,iq)/q(i,k,iqeau)/tnat(2)-1)*1000
153                 if ((deltaD.gt.deltaDmax).or.
154     &                   (deltaD.lt.deltaDmin)) then
155                  write(*,*) 'erreur dans iso_verif_aberrant trac:'
156                  write(*,*) err_msg
157                  write(*,*) 'izone,phase=',izone,phase
158                  write(*,*) 'ixt,ieau=',ixt,ieau
159                  write(*,*) 'q,iq,i,k,=',q(i,k,iq),iq,i,k
160                  write(*,*) 'deltaD=',deltaD
161                  stop
162                 endif !if ((deltaD.gt.deltaDmax).or.
163                endif !if (q(i,k,iq).gt.qmin) then
164                enddo !DO i = ijb,ije
165                enddo  ! do k=1,llm
166c$OMP END DO NOWAIT
167              enddo ! do phase=1,nqo   
168            enddo !do izone=1,ntraceurs_zone
169          endif !if (use_iso(2).and.use_iso(1)) then
170
171          do iiso=1,niso
172           do phase=1,nqo
173              iq=iqiso(iiso,phase)
174c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
175              do k=1,llm
176                DO i = ijb,ije
177                   xtractot=0.0
178                   xiiso=q(i,k,iq)
179                   do izone=1,ntraceurs_zone
180                      iq=iqiso(index_trac(izone,iiso),phase)
181                      xtractot=xtractot+ q(i,k,iq)
182                   enddo !do izone=1,ntraceurs_zone
183                   if ((abs(xtractot-xiiso).gt.errmax).and.
184     :                  (abs(xtractot-xiiso)/
185     :                  max(max(abs(xtractot),abs(xiiso)),1e-18)
186     :                  .gt.errmaxrel)) then
187                  write(*,*) 'erreur detectee par iso_verif_traceurs:'
188                  write(*,*) err_msg
189                  write(*,*) 'iiso,phase=',iiso,phase
190                  write(*,*) 'i,k,=',i,k
191                  write(*,*) 'q(i,k,:)=',q(i,k,:)
192                  stop
193                 endif !if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and.
194                 
195                 ! bidouille pour éviter divergence:
196                 if (abs(xtractot).gt.ridicule) then
197                   do izone=1,ntraceurs_zone
198                     ixt=index_trac(izone,iiso)
199                     q(i,k,iq)=q(i,k,iq)/xtractot*xiiso
200                   enddo !do izone=1,ntraceurs_zone               
201                  endif !if ((abs(xtractot).gt.ridicule) then
202                enddo !DO i = ijb,ije
203              enddo !do k=1,llm
204c$OMP END DO NOWAIT
205           enddo !do phase=1,nqo
206          enddo !do iiso=1,niso
207
208        endif !if (ok_isotrac) then
209
210        endif ! if (ok_isotopes)
211        !write(*,*) 'check_isotopes 198'
212       
213        end
214
215
Note: See TracBrowser for help on using the repository browser.