source: trunk/LMDZ.COMMON/libf/dyn3d/check_isotopes.F

Last change on this file was 1508, checked in by emillour, 9 years ago

Common dynamics:
Updates in the dynamics (seq and ) to keep up with updates
in LMDZ5 (up to LMDZ5 trunk, rev 2325):
IMPORTANT: Modifications for isotopes are only done in dyn3d, not in dyn3dpar

as in LMDZ5 these modifications were done in dyn3dmem.
Related LMDZ5 revisions are r2270 and r2281

  • in dynlonlat_phylonlat:
  • add module "grid_atob_m.F90" (a regridding utility so far only used by phylmd/ce0l.F90, used to be dyn3d_common/grid_atob.F)
  • in misc:
  • follow up updates on wxios.F (add missing_val module variable)
  • in dyn3d_common:
  • pression.F => pression.F90
  • misc_mod.F90: moved from misc to dyn3d_common
  • added new iso_verif_dyn.F
  • covcont.F => covcont.F90
  • infotrac.F90 : add handling of isotopes (reading of corresponding traceur.def for planets not implemented)
  • dynetat0.F => dynetat0.F90 with some code factorization
  • dynredem.F => dynredem.F90 with some code factorization
  • added dynredem_mod.F90: routines used by dynredem
  • iniacademic.F90 : added isotopes-related initialization for Earth case
  • in dyn3d:
  • added check_isotopes.F
  • modified (isotopes) advtrac.F90, caladvtrac.F
  • guide_mod.F90: ported updates
  • leapfrog.F : (isotopes) updates (NB: call integrd with nqtot tracers)
  • qminimium.F : adaptations for isotopes (copied over, except that #include comvert.h is not needed).
  • vlsplt.F: adaptations for isotopes (copied over, except than #include logic.h, comvert.h not needed, and replace "include comconst.h" with use comconst_mod, ONLY: pi)
  • vlspltqs.F : same as vlsplt.F, but also keeping added modification for CP(T)
  • in dyn3dpar:
  • leapfrog_p.F: remove unecessary #ifdef CPP_EARTH cpp flag. and call integrd_p with nqtot tracers (only important for Earth)
  • dynredem_p.F => dynredem_p.F90 and some code factorization
  • and no isotopes-relates changes in dyn3dpar (since these changes have been made in LMDZ5 dyn3dmem).

EM
File size: 7.7 KB
Line 
1        subroutine check_isotopes_seq(q,ip1jmp1,err_msg)
2        USE infotrac, ONLY: nqtot,niso,nqo,ntraceurs_zone,ntraciso,
3     &                      ok_isotopes,ok_isotrac,iqiso,use_iso,
4     &                      indnum_fn_num,tnat,index_trac
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=200.0,deltaDmin=-999.9)
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                  stop
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                  stop
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                  stop
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                  stop
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                  stop
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                  stop
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
Note: See TracBrowser for help on using the repository browser.