Changeset 4491

Timestamp:

Apr 3, 2023, 10:25:44 AM (20 months ago)

Author:

crisi

Message:

Bug corrections in LMDZiso, especially for water tagging

Location:

LMDZ6/trunk/libf/phylmdiso

Files:

• add_phys_tend_mod.F90 (modified) (3 diffs)
• cv30_routines.F90 (modified) (3 diffs)
• cv3p_mixing.F90 (modified) (1 diff)
• cv_driver.F90 (modified) (1 diff)
• fisrtilp.F90 (modified) (5 diffs)
• isotopes_mod.F90 (modified) (1 diff)
• isotopes_routines_mod.F90 (modified) (76 diffs)
• isotopes_verif_mod.F90 (modified) (2 diffs)
• isotrac_routines_mod.F90 (modified) (5 diffs)
• pbl_surface_mod.F90 (modified) (3 diffs)
• phyredem.F90 (modified) (2 diffs)
• physiq_mod.F90 (modified) (15 diffs)
• reevap.F90 (modified) (6 diffs)
• wake.F90 (modified) (3 diffs)

LMDZ6/trunk/libf/phylmdiso/add_phys_tend_mod.F90

@@ -313 +313 @@
 ! Ajout des tendances sur le vent et l'eau liquide
 !======================================================================
+
+#ifdef ISO
+#ifdef ISOVERIF
+        if (iso_eau.gt.0) then
+          call iso_verif_egalite_vect2D( &
+                xt_seri,q_seri, &
+                'add_phys_tend 321a: '//text,ntraciso,klon,klev)
+          call iso_verif_egalite_vect2D( &
+                xtl_seri,ql_seri, &
+                'add_phys_tend 321b: '//text,ntraciso,klon,klev)
+         endif !if (iso_eau.gt.0) then
+#ifdef ISOTRAC
+        call iso_verif_traceur_vect(xt_seri,klon,klev, &
+     &           'add_phys_tend 328a: '//text)
+        call iso_verif_traceur_vect(xtl_seri,klon,klev, &
+     &           'add_phys_tend 328b: '//text)
+#endif
+#endif
+#endif
 
      u_seri(:,:)=u_seri(:,:)+zdu(:,:)
@@ -370 +389 @@
                 'add_phys_tend 138: '//text,ntraciso,klon,klev)
          endif !if (iso_eau.gt.0) then
+#ifdef ISOTRAC
+        call iso_verif_traceur_vect(xt_seri,klon,klev, &
+     &           'add_phys_tend 374a: '//text)
+        call iso_verif_traceur_vect(xtl_seri,klon,klev, &
+     &           'add_phys_tend 374b: '//text)
+#endif
 #endif
 #endif
@@ -495 +520 @@
                         'add_phys_tend 175'//text,ntraciso,klon,klev)
       endif !if (iso_eau.gt.0) then
+#ifdef ISOTRAC
+        call iso_verif_traceur_vect(xt_seri,klon,klev, &
+     &           'add_phys_tend 499'//text)
+#endif
 #endif
 #endif

LMDZ6/trunk/libf/phylmdiso/cv30_routines.F90

@@ -2647 +2647 @@
      &          )
 #ifdef ISO
-    use infotrac_phy, ONLY: ntraciso=>ntiso
+    use infotrac_phy, ONLY: ntraciso=>ntiso, niso
     use isotopes_mod, ONLY: essai_convergence, iso_eau,iso_HDO,ridicule
-    use isotopes_routines_mod, ONLY: appel_stewart_vectall
+    use isotopes_routines_mod, ONLY: appel_stewart_vectall,appel_stewart_debug
 #ifdef ISOVERIF
     use isotopes_verif_mod, ONLY: errmax,errmaxrel, &
@@ -3208 +3208 @@
        enddo
 #endif
+
+        if (1.eq.0) then
         ! appel de appel_stewart_vectorise
         call appel_stewart_vectall(lwork,ncum, &
@@ -3218 +3220 @@
      &          i,inb, & ! altitude: car cas particulier en INB 
      &          na,nd,nloc,cvflag_grav,ginv,1e-16)
+
+        else !if (1.eq.0) then
+          ! truc simple sans fractionnement
+          ! juste pour debuggage
+          call appel_stewart_debug(lwork,nloc,inb,na,i, &
+                evap,water,rpprec,rr,wdtrain, &
+                xtevap,xtwater,xtp,xt,xtwdtrain)
+        endif ! if (1.eq.0) then
+
 
 #ifdef ISOVERIF

LMDZ6/trunk/libf/phylmdiso/cv3p_mixing.F90

@@ -451 +451 @@
           call iso_verif_aberrant_choix( &
                  xtelij(iso_HDO,il,i,j),elij(il,i,j), &
-                 ridicule,deltalim,'cv3p_mixing 1993')
+                 ridicule,deltalim_snow,'cv3p_mixing 1993')
         endif !if (iso_hdo.gt.0) then
 #ifdef ISOTRAC  

LMDZ6/trunk/libf/phylmdiso/cv_driver.F90

@@ -997 +997 @@
 
 #ifdef ISO
-       write(*,*) 'klev=',klev
+!       write(*,*) 'klev=',klev
 #ifdef ISOVERIF
        write(*,*) 'cv_driver 930: apres cv3_unsat'

LMDZ6/trunk/libf/phylmdiso/fisrtilp.F90

@@ -509 +509 @@
         endif !if (iso_eau.gt.0) then
         if (iso_HDO.gt.0) then  
-         if (zrfl(i).gt.ridicule_rain) then
-             call iso_verif_aberrant(zxtrfl(iso_HDO,i) &
-     &           /zrfl(i),'il pleut 316')
-         endif !if (zrfl(i).gt.ridicule_rain) then
+             call iso_verif_aberrant_choix(zxtrfl(iso_HDO,i) &
+     &           ,zrfl(i),ridicule_rain,deltalim_snow,'il pleut 316')
         endif  !if (iso_HDO.gt.0) then  
 #ifdef ISOTRAC   
@@ -1732 +1730 @@
         if (iso_HDO.gt.0) then
            do i=1,klon  
-              if (zcond(i).gt.ridicule) then
-                 call iso_verif_aberrant(zxtcond(iso_HDO,i) &
-     &                  /zcond(i), 'il pleut 637')                
-              endif  !if (zcond(i).gt.ridicule) then
-              IF ((t(i,1)) .LT. RTT) THEN
                   call iso_verif_aberrant_choix(zxtcond(iso_hdo,i), &
      &           zcond(i),ridicule_rain,deltalim_snow,'il pleut 1276')
-              endif
            enddo ! do i=1,klon  
         endif !if ((iso_eau.gt.0).and.(ixt.eq.iso_eau)) then
@@ -2031 +2023 @@
            if (zoliq(i).gt.ridicule) then
             if (iso_HDO.gt.0) then  
-                call iso_verif_aberrant(zxtoliq(iso_HDO,i)/zoliq(i), &
-     &           'il pleut 895a')
+                ! Camille 9 mars 2023: on est moins stricte avec le condensat
+                call iso_verif_aberrant_choix(zxtoliq(iso_HDO,i),zoliq(i), &
+     &           ridicule_rain,deltalim_snow, 'il pleut 895a')
               if (iso_O18.gt.0) then
                 if (iso_verif_O18_aberrant_nostop(zxtoliq(iso_HDO,i)/zoliq(i), &
      &           zxtoliq(iso_O18,i)/zoliq(i),'il pleut 895b').eq.1) then
                    write(*,*) 'i,k,zoliq,zfice=',i,k,zoliq(i),zfice(i)
-                   stop
+                   !stop
                 endif
               endif ! if (iso_HDO.gt.0) then
@@ -2171 +2164 @@
           endif !if (iso_eau.gt.0) then
           if (iso_HDO.gt.0) then
-             if (zoliq(i).gt.ridicule) then
-                call iso_verif_aberrant(zxtoliq(iso_HDO,i)/zoliq(i), &
-     &           'il pleut 963')
-             endif
+                call iso_verif_aberrant_choix(zxtoliq(iso_HDO,i),zoliq(i), &
+     &           ridicule_rain, deltalim_snow, 'il pleut 963')
           endif !if (iso_HDO.gt.0) then
         ! end cam verif
@@ -2482 +2473 @@
             call iso_verif_aberrant(pxtrainfl(iso_HDO,i,k)/prfl(i,k), &
      &            'il pleut 1020')
-          endif !if (prfl(i,k).gt.ridicule_rain) then        
-          if (psfl(i,k).gt.ridicule_rain) then
-            call iso_verif_aberrant(pxtsnowfl(iso_HDO,i,k)/psfl(i,k), &
-     &           'il pleut 1024')
-          endif !if (psfl(i,k).gt.ridicule_rain) then
+          endif !if (prfl(i,k).gt.ridicule_rain) then           
+          call iso_verif_aberrant_choix(pxtsnowfl(iso_HDO,i,k),psfl(i,k), &
+     &           ridicule_rain, deltalim_snow, 'il pleut 1024')
           if (zq(i).gt.ridicule) then
               call iso_verif_aberrant_encadre(zxt(iso_HDO,i)/zq(i), &

LMDZ6/trunk/libf/phylmdiso/isotopes_mod.F90

@@ -302 +302 @@
    IF(getKey('Rdefault',Rdefault,isoName)) CALL abort_physic(modname, 'can''t get Rdefault',1)
    IF(getKey('Rmethox', Rmethox, isoName)) CALL abort_physic(modname, 'can''t get Rmethox', 1)
-   IF(.NOT.Rdefault_smow) Rdefault(:) = 0.0
+
+   IF(.NOT.Rdefault_smow) then
+        Rdefault(:) = 0.0
+        if (iso_eau.gt.0) Rdefault(iso_eau) = 1.0 ! correction Camille 30 mars 2023
+   ENDIF
+   write(*,*) 'Rdefault=',Rdefault
 
    !--- Sensitivity test: no kinetic effect in sfc evaporation

LMDZ6/trunk/libf/phylmdiso/isotopes_routines_mod.F90

@@ -154 +154 @@
 &                   'revap_ilp 131')
        enddo
+#ifdef ISOTRAC
+     call iso_verif_traceur(zxtrfl_ancien(1,i), &
+&           'iso_revap_fisrtilp 158: debut')
+#endif
 #endif
     endif !if (zrfln(i).gt.ridicule*1e-2) then  
@@ -193 +197 @@
         endif !if (zq(i).gt.ridicule) then
      endif !if ((iso_HDO.gt.0.and.(iso_O18.gt.0) then
+#ifdef ISOTRAC
+     call iso_verif_traceur(zxtrfl_ancien(1,i), &
+&           'iso_revap_fisrtilp 201: debut quand pas de precip')
+#endif
 !     write(*,*) 'iso_routines tmp 184'
 #endif               
@@ -504 +512 @@
 !             stop
   endif
-  enddo
+  enddo !do iiso=1,niso
 enddo !do i=1,ncas_evap_liq
 #endif          
@@ -514 +522 @@
 &       izone,zqevfl(1),Exi(1,1),fac_fluxtomixratio(1), &
 &       xtrevap_tag(1,1),1,hdiag(1))
+        ! dans cette routine, zxtrfl reçoit zxtrfln_cas
 
   enddo !do izone=1,ntraceurs_zone
+
 #ifdef ISOVERIF
 do i=1,ncas_evap_liq
@@ -522 +532 @@
 &              0.0,'revap_ilp 414')
 enddo
+     call iso_verif_traceur(zxtrfl(1,cas_evap_liq(i)), &
+&           'iso_revap_fisrtilp 470a: apres stewart_explicite_vectall')
 enddo !do i=1,ncas
 #endif            
@@ -536 +548 @@
      call iso_verif_traceur(zxt(1,cas_evap_liq(i)), &
 &           'iso_revap_fisrtilp 282')
+     call iso_verif_traceur(zxtrfl(1,cas_evap_liq(i)), &
+&           'iso_revap_fisrtilp 804a')
      call iso_verif_traceur(zxtrfln(1,cas_evap_liq(i)), &
-&           'iso_revap_fisrtilp 804')
+&           'iso_revap_fisrtilp 804b')
      do ixt=1,ntraciso
         call iso_verif_positif_choix(zxt(ixt,cas_evap_liq(i)), &
@@ -676 +690 @@
    if (iso_HDO.gt.0) then
       do i=1,ncas_evap_glace
-         if (zrfln_cas(i).gt.ridicule_rain) then
-                call iso_verif_aberrant( &
-&                  (zxtrfln_cas(iso_HDO,i) &
-&                 /zrfln_cas(i)), 'iso_revap_fisrtilp 4563')
-         endif
+                call iso_verif_aberrant_choix(zxtrfln_cas(iso_HDO,i), zrfln_cas(i), &
+                  ridicule_rain,deltalim_snow, 'iso_revap_fisrtilp 4563')
        enddo !do i=1,ncas_evap_glace
    endif !if ((iso_HDO.gt.0).and.(ixt.eq.iso_HDO)) then
@@ -789 +800 @@
      call iso_verif_traceur(zxt(1,cas_evap_glace(i)), &
 &           'iso_revap_fisrtilp 1033')
+     call iso_verif_traceur(zxtrfl(1,cas_evap_glace(i)), &
+&           'iso_revap_fisrtilp 1035a')
      call iso_verif_traceur(zxtrfln(1,cas_evap_glace(i)), &
-&           'iso_revap_fisrtilp 1035')
+&           'iso_revap_fisrtilp 1035b')
   enddo
 #endif          
@@ -830 +843 @@
     call iso_verif_traceur(zxt(1,i),'iso_revap_fisrtilp 532')
     call iso_verif_traceur(zxtrfln(1,i), &
-&           'iso_revap_fisrtilp 533') 
+&           'iso_revap_fisrtilp 533a') 
+    call iso_verif_traceur(zxtrfl(1,i), &
+&           'iso_revap_fisrtilp 533b')
     do ixt=1,ntraciso
         call iso_verif_positif_choix(zxt(ixt,i),0.0, &
@@ -836 +851 @@
     enddo
 enddo !do i=1,klon
-        write(*,*) 'revap_ilp 814: sortie'
+        !write(*,*) 'revap_ilp 814: sortie'
 #endif       
 #endif          
@@ -2384 +2399 @@
     endif !if (iso_eau.gt.0) then
     if (abs(evap(i)).gt.ridicule_rain) then
-      if (iso_HDO.gt.0) then                  
-         if (iso_verif_aberrant_nostop(xtevap(iso_HDO,i)/evap(i), &
-&           'iso_surf>iso_rosee_givre 3193').eq.1) then    
+      if (iso_HDO.gt.0) then                      
+         if (iso_verif_aberrant_choix_nostop(xtevap(iso_HDO,i),evap(i), &
+&           ridicule_rain,deltalim_snow,'iso_surf>iso_rosee_givre 3193').eq.1) then    
                 write(*,*) 'zxtalphai(iso_HDO)=',zxtalphai(iso_HDO)
                 write(*,*) 'deltaD1eff=',deltaD(xt1lay(iso_HDO,i)/q1lay(i))
                 write(*,*) 'tsurf(i)=',tsurf(i)
                 write(*,*) 'q1lay(i)=',q1lay(i)   
-                stop  
+                !stop  
            endif !if (iso_verif_aberrant_nostop
       endif  !if (iso_HDO.gt.0) then  
@@ -2534 +2549 @@
 
 ! quelques verifs de bilan d'eau
-#ifdef ISOVERIF      
+#ifdef ISOVERIF   
 do il=1,ncas
   do ixt=1,niso
@@ -2647 +2662 @@
 &             'stewart_explicite_vectall 220')
     enddo
-#endif
-#ifdef ISOVERIF
-      if (iso_eau.gt.0) then
+    if (iso_eau.gt.0) then
         call iso_verif_egalite_choix( &
 &              (Exi(iso_eau,il)*fac_ftmr(il)), &
@@ -2669 +2682 @@
 &                   *fac_ftmr(il))),'stewart_explicite 214')
        endif !if ((iso_HDO.gt.0).and.
-       if ((debug.eq.1).and.(il.eq.il_debug)) then
-           write(*,*) 'stewart_explicit 224: cas Pqisup<=0'
-           write(*,*) 'Eqi(il),deltaD=',Eqi(il), &
-&                   deltaD((Exi(iso_HDO,il)/Eqi(il)))
-       endif
+
 #endif              
 else !if (Pqisup.eq.0) then
@@ -2679 +2688 @@
 h(il)=qeff(il)/qs(il)       
 h(il)= MAX(MIN(h(il),1.0),0.0)
-#ifdef ISOVERIF        
+#ifdef ISOVERIF
 call iso_verif_positif(h(il)-thumxt1,'stewart_explicit 209') 
 #endif        
@@ -2739 +2748 @@
       endif !if ((iso_HDO.gt.0).and.
     endif !if (iso_HDO.gt.0) then
-    if ((debug.eq.1).and.(il.eq.il_debug)) then
+    if ((debug.eq.1).and.(il.eq.il_debug).and.(Eqi(il).gt.0.)) then
            write(*,*) 'stewart_explicit 302: cas evap~0'
            write(*,*) 'deltaDv est inchangé:',deltaD( &
@@ -2779 +2788 @@
 f(il)=m(il)/m0(il)    
 ! verifs
-#ifdef ISOVERIF        
+#ifdef ISOVERIF       
 call iso_verif_positif((m(il)), &
 &           'stewart_explicite 173')
@@ -2836 +2845 @@
     ! rajout verif 4 sept 2009
     if (iso_HDO.gt.0) then
-     if (Pqisup(il).gt.ridicule) then
-        call iso_verif_aberrant((Rl0(iso_HDO,il)), &
-&                   'stewart_explicite 368')
-     endif
-    endif
+        call iso_verif_aberrant_choix(Rl0(iso_HDO,il)*Pqisup(il),Pqisup(il), &
+&                   ridicule_rain,deltalim_snow,'stewart_explicite 368')
+    endif !if (iso_HDO.gt.0) then
   endif !(iso_eau.gt.0)
 #endif
@@ -2868 +2875 @@
 &                   'stewart_explicite 271')
    enddo !do ixt=1,niso   
-#endif
-#ifdef ISOVERIF  
    if (iso_eau.gt.0) then
       call iso_verif_egalite_choix( &
@@ -2923 +2928 @@
        endif !if (qp0(il)+Eqi(il)*fac_ftmr(il).gt.ridicule) then
      endif !if (iso_HDO.gt.0)
-     if ((debug.eq.1).and.(il.eq.il_debug)) then
+     if ((debug.eq.1).and.(il.eq.il_debug).and.(Eqi(il).gt.0.)) then
          write(*,*) 'stewart_explicit 442: tout se réévapore'
          write(*,*) 'Eqi(il),deltaD=',Eqi(il), &
@@ -2943 +2948 @@
 if ((h(il).gt.0.99).or. &
 &           (h(il).gt.0.98).and.(f(il).lt.1e-3)) then
-#ifdef ISOVERIF          
-!            write(*,*) 'stewart_explicit 191: cas h=1: il=',il
-#endif            
     do ixt=1,niso
     interm(ixt,il)=alphap(ixt,il)*(1.0-h(il)) &
@@ -2981 +2983 @@
 &           'stewart_explicite 261')
     enddo !do ixt=1,niso
-#endif           
-#ifdef ISOVERIF
     if (iso_eau.gt.0) then
         call iso_verif_egalite_choix( &
@@ -3007 +3007 @@
       
      if (iso_HDO.gt.0) then 
-       if (Pqiinf(il).gt.ridicule_rain) then
-           if (iso_verif_aberrant_nostop( &
-&                  (Pxtiinf(iso_HDO,il)/Pqiinf(il)), &
-&                   'stewart_explicie 248').eq.1) then
+       if (iso_verif_aberrant_choix_nostop(Pxtiinf(iso_HDO,il),Pqiinf(il), &
+&                ridicule_rain,deltalim_snow,'stewart_explicite 248').eq.1) then
              write(*,*) 'cas reeq totale, il=',il
              write(*,*) 'deltaDl0=',deltaD( &
@@ -3019 +3017 @@
 &                   (Rb(iso_hdo,il)))
              stop
-           endif 
-       endif  !if ((iso_HDO.gt.0).and.(ixt.eq.iso_HDO).and.
+       endif !if (iso_verif_aberrant_choix_nostop
        if (qp0(il)+Eqi(il)*fac_ftmr(il).gt.ridicule) then
          call iso_verif_aberrant(( &
@@ -3027 +3024 @@
        endif !if ((iso_HDO.gt.0).and.
      endif  !if (iso_HDO.gt.0) then
-    if ((debug.eq.1).and.(il.eq.il_debug)) then
+     if ((debug.eq.1).and.(il.eq.il_debug).and.(Eqi(il).gt.0.)) then
          write(*,*) 'stewart_explicit 526: cas h~1: rééq'
          write(*,*) 'Eqi(il),deltaD=',Eqi(il), &
@@ -3034 +3031 @@
 &                   (Rb0(iso_hdo,il))),deltaD( &
 &                   (Rl0(iso_hdo,il)))
-    endif 
+     endif !if ((debug.eq.1).and.(il.eq.il_debug)) then
 #endif
     ! end verifs
 
 else if ((f(il).gt.0.998).and. &
-&           (Eqi(il)*fac_ftmr(il).lt.1e-2*qp0(il))) then
+&           (Eqi(il)*fac_ftmr(il).lt.1e-2*qp0(il))) then ! if ((h(il).gt.0.99).or. 
 
 !*** cas particulier pour éviter imprécisions numériques:
@@ -3083 +3080 @@
 &           'stewart_explicite 397')
    enddo !do ixt=1,niso 
-#endif
-#ifdef ISOVERIF
    if (iso_eau.gt.0) then
         call iso_verif_egalite_choix( &
@@ -3167 +3162 @@
            stop
          endif !if (iso_verif_aberrant_nostop((
-       endif !if ((iso_HDO.gt.0).and.
+       endif !if (qp0(il)+Eqi(il)*fac_ftmr(il).gt.ridicule) then
      endif !!if ((iso_HDO.gt.0)
-     if ((debug.eq.1).and.(il.eq.il_debug)) then
+     if ((debug.eq.1).and.(il.eq.il_debug).and.(Eqi(il).gt.0.)) then
          write(*,*) 'stewart_explicit 663: cas où réévap faible'
          write(*,*) 'ordre 1 pour la vapeur et le liquide'
@@ -3271 +3266 @@
              endif !if (qp0(il)+Eqi(il)*fac_ftmr(il).gt.ridicule) then
           endif !if ((iso_HDO.gt.0)
-     if ((debug.eq.1).and.(il.eq.il_debug)) then
+     if ((debug.eq.1).and.(il.eq.il_debug).and.(Eqi(il).gt.0.)) then
          write(*,*) 'stewart_explicit 767: cas de réévap sèche'
          write(*,*) 'distill de Rayleigh'
@@ -3317 +3312 @@
 &                   'stewart_explicite 467')
     enddo !do ixt=1,niso    
-#endif          
-#ifdef ISOVERIF  
           if (iso_eau.gt.0) then
             call iso_verif_egalite_choix( &
@@ -3360 +3353 @@
             endif !if (qp0(il)+Eqi(il)*fac_ftmr(il).gt.ridicule) then
           endif !if (iso_HDO.gt.0)
-      if ((debug.eq.1).and.(il.eq.il_debug)) then
+      if ((debug.eq.1).and.(il.eq.il_debug).and.(Eqi(il).gt.0.)) then
          write(*,*) 'stewart_explicit 831: flux de masse vap~0'
          write(*,*) 'Eqi(il),deltaD=',Eqi(il), &
@@ -3429 +3422 @@
         enddo !do ixt=1,niso    
 #endif
-#ifdef ISOVERIF   
+#ifdef ISOVERIF
           if (iso_eau.gt.0) then
             call iso_verif_egalite_choix( &
@@ -3546 +3539 @@
         enddo !do ixt=1,niso 
 #endif
-#ifdef ISOVERIF      
+#ifdef ISOVERIF   
           if (iso_eau.gt.0) then
             call iso_verif_egalite_choix( &
@@ -3750 +3743 @@
        ! pour meilleure convergence numérique:
        !xtnew=qp0+Eqi*fac_ftmr
-   endif   ! if (iso_eau.gt.0).and.(ixt.eq.iso_eau)    
-
+   endif   ! if (iso_eau.gt.0).and.(ixt.eq.iso_eau)  
    if (iso_HDO.gt.0) then
      if (Pqiinf(il).gt.ridicule_rain) then
@@ -3778 +3770 @@
 &                   xtnew(iso_HDO,il)/(qp0(il)+Eqi(il) &
 &                   *fac_ftmr(il))),'stewart_explicite 912')
-        endif !if (qp0(il)+Eqi(il)*fac_ftmr(il).gt.ridicule) then
-      endif !if (iso_HDO.gt.0)
-      if ((debug.eq.1).and.(il.eq.il_debug)) then
+       endif !if (qp0(il)+Eqi(il)*fac_ftmr(il).gt.ridicule) then
+       if ((debug.eq.1).and.(il.eq.il_debug).and.(Eqi(il).gt.0.)) then
          write(*,*) 'stewart_explicit 991: fcas général'
          write(*,*) 'mais avec formule simplifiée'
-         write(*,*) 'il,Eqi(il),deltaD=',il,Eqi(il), &
-&                   deltaD((Exi(iso_HDO,il)/Eqi(il)))
-      endif
-#endif     
+         write(*,*) 'il,Eqi(il)=',il,Eqi(il)
+         write(*,*) 'deltaD=',deltaD((Exi(iso_HDO,il)/Eqi(il)))
+       endif
+      endif !if (iso_HDO.gt.0)  
+#endif
   ! end verifs  
     
@@ -3838 +3830 @@
 ! compression
 if (ncas_Jsimple+ncas_rieman.gt.0) then
-!#ifdef ISOVERIF          
+!#ifdef ISOVERIF         
 !        write(*,*) 'stewart_explicite_vectall 873:compression_calculJ' 
 !#endif        
@@ -3927 +3919 @@
     ! ******* traitement vectoriel du cas Rieman (=2533)
 if (ncas_rieman.gt.0) then
-!#ifdef ISOVERIF
-!!          write(*,*) 'traitement vectoriel rieman: x',ncas_rieman
-!          do icas_rieman=1+ncas_Jsimple,ncas_rieman+ncas_Jsimple
-!!            write(*,*) 'ntot_cas(icas_rieman)=',ntot_cas(icas_rieman)
-!          call iso_verif_positif(float(ntot_cas(icas_rieman))-1.0,
-!     :           'stewart_expl 984: ntot faux') 
-!          enddo   !do icas_rieman=1,ncas_rieman       
-!#endif            
+            
    icas_rieman=1+ncas_Jsimple
-!           write(*,*) 'stewart_expl 988 tmp: icas_rieman=',icas_rieman
-!           write(*,*) 'qp0_cas(1)=',qp0_cas(1)
-!           write(*,*) 'A_cas(1)=',A_cas(1)
-!           write(*,*) 'm0_cas(1)=',m0_cas(1)
-!           write(*,*) 'm_cas(1)=',m_cas(1)
-!           if (iso_eau.gt.0) then
-!           write(*,*) 'beta_cas(iso_eau,1)=',beta_cas(iso_eau,1)
-!           write(*,*) 'gama_cas(iso_eau,1)=',gama_cas(iso_eau,1)
-!           endif
-!           write(*,*) 'f_cas(1)=',f_cas(1)
-!           write(*,*) 'g_cas(1)=',g_cas(1)
-!#ifdef rieman           
-!           call integrale_rieman_vectall
-!     :       (ncas_rieman,m_cas(icas_rieman),
-!     :       J(1,icas_rieman),e(1,icas_rieman),
-!     :       qp0_cas(icas_rieman),A_cas(icas_rieman),
-!     :       m0_cas(icas_rieman),beta_cas(1,icas_rieman),
-!     :       gama_cas(1,icas_rieman),f_cas(icas_rieman),
-!     :       g_cas(icas_rieman),ntot_cas(icas_rieman))
-!#else
+
    call integrale_gauss_vectall &
 &       (ncas_rieman,m_cas(icas_rieman), &
@@ -3964 +3930 @@
 !     :       g_cas(icas_rieman),ntot_cas(icas_rieman))
 &       g_cas(icas_rieman))      
-!#endif
-!#ifdef ISOVERIF
-!           do il=1+ncas_Jsimple,ncas_rieman+ncas_Jsimple
-!             do ixt=1,niso
-!               call integrale_Rieman_precision(m_cas(il),m0_cas(il),
-!     :          Jtmp,etmp,ntot_cas(il)*1e2,
-!     :           qp0_cas(il),A_cas(il),m0_cas(il),
-!     :          beta_cas(ixt,il),gama_cas(ixt,il),f_cas(il),g_cas(il))
-!               call iso_verif_egalite_choix((Jtmp),
-!     :           (J(ixt,il)),
-!     :           'stewart_exp 999: test rieman',errmax,errmaxrel)
-!               if ((iso_eau.gt.0).and.(ixt.eq.iso_eau)) then
-!                 write(*,*) 'stew exp tmp 1005: il,J(iso_eau,il),Jtmp=',
-!     :           il,J(iso_eau,il),Jtmp
-!               endif
-!             enddo             
-!           enddo
-!#endif           
 
   endif !if (ncas_rieman.gt.0) then
@@ -4008 +3956 @@
 &          -beta_cas(ixt,il)*gama_cas(ixt,il)*r_jqp0(ixt,il) &
 &                   /g_cas(il)/g_cas(il))
-!#ifdef ISOVERIF            
-!            if ((iso_eau.gt.0).and.(iso_eau.eq.ixt)) then
-!              write(*,*) 'stewart_explicite tmp 1071: il=',il
-!              if (il.le.ncas_Jsimple) then
-!                write(*,*) 'cas_Jsimple(il)=',cas_Jsimple(il)
-!              else !if (il.le.ncas_Jsimple) then
-!                write(*,*) 'cas_rieman(il)=',cas_rieman(il)
-!              endif !if (il.le.ncas_Jsimple) then
-!              write(*,*) 'f_cas(il),beta_cas(ixt,il),gama_cas(ixt,il)=',
-!     :                   f_cas(il),beta_cas(ixt,il),gama_cas(ixt,il)
-!              write(*,*) 'g_cas(il),r_jqp0(ixt,il),r_jl0(ixt,il)=',
-!     ;          g_cas(il),r_jqp0(ixt,il),r_jl0(ixt,il)
-!              write(*,*) 'A_cas(il)=',A_cas(il)
-!              write(*,*) 'pond Rl0=',(f_cas(il)**beta_cas(ixt,il))
-!     :          *(g_cas(il)**(-beta_cas(ixt,il)*gama_cas(ixt,il)))
-!     :          +beta(ixt,il)*gama_cas(ixt,il)*r_jqp0(ixt,il)
-!     :                   /f_cas(il)/g_cas(il)
-!              write(*,*) 'pond Rb0=',gama_cas(ixt,il)*beta_cas(ixt,il)
-!     :                   *r_jl0(ixt,il)/f_cas(il)/g_cas(il)
-!              write(*,*) 'pondRl0=fac1*fac2+t3'
-!              write(*,*) 'fac1=',
-!     :           f_cas(il)**beta_cas(ixt,il)
-!              write(*,*) 'fac2=',
-!     :             g_cas(il)**(-beta_cas(ixt,il)*gama_cas(ixt,il)) 
-!              write(*,*) 't3=',
-!     :           beta_cas(ixt,il)*gama_cas(ixt,il)*r_jqp0(ixt,il)
-!     :                   /f_cas(il)/g_cas(il)     
-!            endif
-!#endif
 
     Pxtiinf_cas(ixt,il)=Pqiinf_cas(il)*Rl(ixt,il) 
@@ -4097 +4016 @@
     enddo
 #endif
-#ifdef ISOVERIF            
+#ifdef ISOVERIF         
     if (iso_eau.gt.0) then
       if (iso_verif_egalite_choix_nostop( &
@@ -4411 +4330 @@
    if (iso_HDO.gt.0) then
     if (Pqiinf(il).gt.ridicule_rain) then
-      if (iso_verif_aberrant_nostop( &
-&                  (Pxtiinf(iso_HDO,il)/Pqiinf(il)), &
-&                   'stewart_explicie 871').eq.1)  then
+      if (iso_verif_aberrant_choix_nostop(Pxtiinf(iso_HDO,il),Pqiinf(il),ridicule_rain,deltalim_snow, &
+&                   'stewart_explicite 871').eq.1)  then
       write(*,*) 'deltaDl0=',deltaD( &
 &           (Rl0(iso_HDO,il)))
@@ -4428 +4346 @@
     endif !if (iso_HDO.gt.0)
 
-    if ((debug.eq.1).and.(il.eq.il_debug)) then
+    if ((debug.eq.1).and.(il.eq.il_debug).and.(Eqi(il).gt.0.)) then
           write(*,*) 'stewart_explicit 1558: cas avec calcul J'
           write(*,*) 'Eqi(il),deltaD=',Eqi(il), &
@@ -4478 +4396 @@
 
 #ifdef ISOVERIF
-!          write(*,*) 'stewart_explicite vectall 1179: fin'
+          write(*,*) 'stewart_explicite vectall 1179: fin'
 #endif         
 
@@ -4702 +4620 @@
 &           'stewart_sublim_nofrac 39')              
 endif !if ((iso_eau.gt.0).and.(ixt.eq.iso_eau)) then
-if ((iso_HDO.gt.0).and. &
-&           (Pqisup(il).gt.ridicule_rain)) then
-    call iso_verif_aberrant( &
-&           (Pxtisup(iso_HDO,il)/Pqisup(il)), &
-&           'stewart_sublim_nofrac 40')     
+if (iso_HDO.gt.0) then ! Camille 9 mars 2023: moins stricte pour condensat
+    call iso_verif_aberrant_choix(Pxtisup(iso_HDO,il),Pqisup(il), &
+&           ridicule_rain,deltalim_snow, 'stewart_sublim_nofrac 40') 
 endif !if ((iso_HDO.gt.0).and.(ixt.eq.iso_HDO).and.
 enddo !do il=1,ncas
@@ -4834 +4750 @@
         enddo ! do ixt=1,niso
         ! verif que deltaD(Pqiinf) raisonable
-        if ((iso_HDO.gt.0).and. &
-     &           (Pqiinf(il).gt.ridicule_rain)) then
-            call iso_verif_aberrant( &
-     &           (Pxtiinf(iso_HDO,il)/Pqiinf(il)), &
-     &           'stewart_sublim 175')
+        if (iso_HDO.gt.0) then
+            call iso_verif_aberrant_choix(Pxtiinf(iso_HDO,il),Pqiinf(il), &
+     &           ridicule_rain,deltalim_snow, 'stewart_sublim 175')
         endif !if ((iso_HDO.gt.0).and.(ixt.eq.iso_HDO).and.
         if (iso_eau.gt.0) then
@@ -4923 +4837 @@
 !         integer ntot_cas(ncas)
          integer il,ixt        
-         integer coeff_precision
-         parameter (coeff_precision=0.2)
 
 
@@ -10827 +10739 @@
      &          Exi_cas(1,1),Exi(1,1), &
 #endif       
-     &          xtp_avantevaptrac_cas,0,hdiag(1)) ! hdiag ne sera pas utilisé
+     &          xtp_avantevaptrac_cas,0,hdiag(1)) ! hdiag ne sera pas utilise
 
        enddo ! do izone=1,ntraceurs_zone
@@ -12615 +12527 @@
                   if (iso_verif_aberrant_choix_nostop( &
      &               zxtliq(iso_HDO,i),cond(i), &
-     &               ridicule,deltalim, &
+     &               ridicule,deltalim_snow, &
      &               'condiso_liq_ice_vectall 32b').eq.1) then
                     write(*,*) 'deltaDvap=',deltaD(xt(iso_hdo,i)/qt(i))
@@ -12636 +12548 @@
                         write(*,*) 'deltaO18(zxtliq/cond)=',deltaO(zxtliq(iso_O18,i)/cond(i))
                         write(*,*) 'tcond(i)=',tcond(i)-t_coup,'°C'
-                        stop
+                        !stop ! Camille 9 mars 2023: trop strict
                     endif  !if (iso_verif_O18_aberrant_nostop(     
                     endif ! if (iso_O18.gt.0) then      
@@ -12714 +12626 @@
                         write(*,*) 'zxtalphai(iso_O18,i)=',zxtalphai(iso_O18,i)
                         write(*,*) 'xt(1:niso,i)=',xt(1:niso,i)
-                        stop
+                        !stop ! Camille 9 mars 2023: trop strict
                     endif  !if (iso_verif_O18_aberrant_nostop(     
                     endif ! if (iso_O18.gt.0) then
@@ -12765 +12677 @@
           if (zfice(i).gt.0.9) then
               if (iso_verif_aberrant_choix_nostop( &
-     &          zxtice(iso_HDO,i),cond(i),ridicule,deltalim_snow, &
+     &          zxtice(iso_HDO,i),cond(i),ridicule,deltalim_snow, & 
+                ! Camille 9 mars 2023: pour le condensat, on laisse plus de
+                ! marge
      &          'condiso_liq_ice_vect 412').eq.1) then
                 write(*,*) 'debug condiso_liq_ice_vect 449: i,zfice=', &
@@ -13470 +13384 @@
           endif !if (iso_eau.gt.0) then
 #ifdef ISOTRAC
-          call iso_verif_traceur(xtsnow_evap(1,i), &
+!          call iso_verif_traceur(xtsnow_evap(1,i), &
+!     &           'gestion neige 2146') ! attention car snow_evap parfois
+!     négatif -> il ne faut pas passer dans les verifs de positivité.
+          call iso_verif_traceur_justmass(xtsnow_evap(1,i), &
      &           'gestion neige 2146')
 #endif          
@@ -14141 +14058 @@
         if (iso_HDO.gt.0) then
             call iso_verif_aberrant_choix(-xtsol_evap(iso_HDO,i), &
-     &           sol_evap(i),ridicule_evap,deltalim, &
-     &           'calcul_iso_surf_sic 257')
+     &           sol_evap(i),ridicule_evap,deltalim_snow, &
+     &           'calcul_iso_surf_sic 257_sol_evap')
         endif
 #endif            
@@ -14162 +14079 @@
               endif !if (iso_eau.gt.0) then
               if (iso_HDO.gt.0) then
-                if (evap(i).gt.ridicule_evap) then  
-                   call iso_verif_aberrant(xtevap(iso_HDO,i)/evap(i), &
-     &                  'calcul_iso_surf_sic 257')
-               endif !if (evap(i).gt.ridicule_evap) then  
+                   call iso_verif_aberrant_choix(xtevap(iso_HDO,i),evap(i), &
+     &                  ridicule_evap,deltalim_snow,'calcul_iso_surf_sic 257_evap')
               endif !if (iso_eau.gt.0) then
 #ifdef ISOTRAC
@@ -14323 +14238 @@
         enddo !do i=1,knon
         endif !if (iso_HDO.gt.0) then
-        
-        if (iso_eau.gt.0) then  
+            
         do i=1,knon
+        if (iso_eau.gt.0) then
           call iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i),  &
-     &         'calcul_iso_surf_lic_vectall 587',errmax,errmaxrel)
-        enddo
+     &         'calcul_iso_surf_lic_vectall 587a',errmax,errmaxrel)
         endif
+        if (iso_HDO.gt.0) then
+          call iso_verif_aberrant_choix(xtsnow(iso_HDO,i), &
+     &               snow(i),ridicule_snow,deltalim_snow, &
+     &               'calcul_iso_surf_lic 587b')
+        endif
+        enddo !do i=1,knon
 #endif             
         
@@ -14476 +14396 @@
      &                  'calcul_iso_surf_lic 363',errmax,errmaxrel)  
               endif !if (iso_eau.gt.0) then
-              if (iso_HDO.gt.0) then  
-                if (snow(i).gt.ridicule_snow) then  
-                 call iso_verif_aberrant(xtsnow(iso_HDO,i)/snow(i), &
-     &                  'calcul_iso_surf_lic 367')
-                endif !if (snow(i).gt.ridicule_snow) then 
+              if (iso_HDO.gt.0) then   
                 call iso_verif_aberrant_choix(xtsnow(iso_HDO,i), &
-     &               snow(i),ridicule,deltalim_snow, &
+     &               snow(i),ridicule_snow,deltalim_snow, &
      &               'calcul_iso_surf_lic 797') 
-                if (abs(evap(i)).gt.ridicule_evap) then  
-                 call iso_verif_aberrant(xtevap(iso_HDO,i)/evap(i), &
-     &                  'calcul_iso_surf_lic 369')
-                endif ! if (evap(i).gt.ridicule_evap) then 
+                call iso_verif_aberrant_choix(xtevap(iso_HDO,i),evap(i), &
+     &                  ridicule_evap,deltalim_snow, 'calcul_iso_surf_lic 369')
               endif !if (iso_eau.gt.0) then
 #ifdef ISOTRAC
@@ -15795 +15709 @@
               endif !if (iso_eau.gt.0) then
               if (iso_HDO.gt.0) then  
-                 if (snow(i).gt.ridicule_snow) then 
-                    call iso_verif_aberrant(xtsnow(iso_HDO,i)/snow(i), &
-     &                  'calcul_iso_surf_ter 749')
-                 endif !if (snow(i).gt.ridicule_snow) then
+                 call iso_verif_aberrant_choix(xtsnow(iso_HDO,i),snow(i), &
+     &                  ridicule_snow,deltalim_snow, 'calcul_iso_surf_ter 749')
                  call iso_verif_aberrant_choix(xtsnow(iso_HDO,i), &
      &               snow(i),ridicule,deltalim_snow, &
@@ -15809 +15721 @@
                      write(*,*) 'sol_evap,snow_evap=', &
      &                   sol_evap(i),snow_evap(i)
-                     write(*,*) 'deltaDsol_evap=', &
+                     if (sol_evap(i).gt.ridicule_evap)write(*,*) 'deltaDsol_evap=', &
      &                   deltaD(xtsol_evap(iso_hdo,i)/sol_evap(i))
-                     write(*,*) 'deltaDsnow_evap=', &
+                     if (snow_evap(i).gt.ridicule_evap)write(*,*) 'deltaDsnow_evap=', &
      &                   deltaD(xtsnow_evap(iso_hdo,i)/snow_evap(i))
                      write(*,*) 'deltaD1new=',deltaD( &
@@ -16094 +16006 @@
           endif !if (iso_eau.gt.0) then
           if (iso_HDO.gt.0) then
-            if (snow(i,nsrf).gt.ridicule_snow) then
-            call iso_verif_aberrant(xtsnow(iso_hdo,i,nsrf)/snow(i,nsrf), &
-     &                'phyisoetat0 117')
-            endif
+            call iso_verif_aberrant_choix(xtsnow(iso_hdo,i,nsrf),snow(i,nsrf), &
+     &                ridicule_snow, deltalim_snow, 'phyisoetat0 117')
           endif !if (iso_eau.gt.0) then
         enddo !do nsrf=1,nbsrf
@@ -16180 +16090 @@
         real deltaD_snow_fall_O18,deltaD_rain_fall_O18
         real alpha(niso),kcin(niso)
-        real q0,h0
-        parameter (q0=20e-3,h0=0.7)
 !      character*50 text
         
@@ -16314 +16222 @@
         do k=1,klev
           do ixt=1,niso
-           RMerlivat(ixt)=toce(ixt)/alpha(ixt) &
-     &                  *(1-kcin(ixt))/(1.0-kcin(ixt)*h0)
-           xt_ancien(ixt,i,k)=q_ancien(i,k)*RMerlivat(ixt) &
-     &                  *(min(q0,q_ancien(i,k))/q0)**(alpha(ixt)-1.0)
+           call iso_init_ideal(q_ancien(i,k),xt_ancien(ixt,i,k),ixt, &
+                alpha(ixt),kcin(ixt),toce(ixt))
+
            if (q_ancien(i,k).gt.ridicule) then
            xtl_ancien(ixt,i,k)=ql_ancien(i,k)*alpha(ixt) &
@@ -16338 +16245 @@
      &           'phyisoetat0 16067')  
         enddo !do ixt=1,niso
-        if ((k.eq.1).and.(iso_HDO.gt.0).and.(iso_O18.gt.0) &
-                .and.(abs(q_ancien(i,k)-q0).lt.1e-3)) then
-                ! vérifier qu'on est proche de la fermeture de Merlivat
-            write(*,*) 'i,k=',i,k
-            write(*,*) 'q_ancien(i,k)=',q_ancien(i,k)
-            write(*,*) 'deltaD=',deltaD(xt_ancien(iso_HDO,i,k)/q_ancien(i,k))
-            write(*,*) 'deltaDM=',deltaD(RMerlivat(iso_HDO))
-            write(*,*) 'deltaO=',deltaO(xt_ancien(iso_O18,i,k)/q_ancien(i,k))
-            write(*,*) 'deltaOM=',deltaO(RMerlivat(iso_O18))
-            write(*,*) 'dexcess=',dexcess(xt_ancien(iso_HDO,i,k)/q_ancien(i,k), &
-                    xt_ancien(iso_O18,i,k)/q_ancien(i,k))
-            write(*,*) 'dexcessM=',dexcess(RMerlivat(iso_HDO),RMerlivat(iso_O18))
-            write(*,*) 'kcin=',kcin
-            write(*,*) 'toce=',toce
-            write(*,*) 'alpha=',alpha
-            call iso_verif_positif(20.0-abs(-80.0-deltaD(xt_ancien(iso_HDO,i,k)/q_ancien(i,k))), &
-                    'phyisoetat0 16398a')   
-            call iso_verif_positif(5.0-abs(10.0-dexcess(xt_ancien(iso_HDO,i,k)/q_ancien(i,k), &
-                    xt_ancien(iso_O18,i,k)/q_ancien(i,k))),'phyisoetat0 16398b') 
-        endif 
+
+        ! Camille 7 mars 2023: ajout d'un check
+        if ((i.eq.1).and.(k.eq.1).and.(iso_HDO.gt.0)) then
+        write(*,*) 'phyisoetat0 16362: q_ancien(1,1)=',q_ancien(1,1) 
+        write(*,*) 'deltaD_ancien=',deltaD(xt_ancien(iso_HDO,i,k)/q_ancien(i,k))
+        write(*,*) 'xt_ancien(:,i,k)=',xt_ancien(:,i,k)
+        endif !if ((i.eq.1).and.(k.eq.1)) then
+
+        if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then 
+            if (q_ancien(i,k).gt.ridicule) then  
+               if (iso_verif_o18_aberrant_nostop( &
+     &              xt_ancien(iso_HDO,i,k)/q_ancien(i,k), &
+     &              xt_ancien(iso_O18,i,k)/q_ancien(i,k), &
+     &              'phyisoetat0 16366 q_ancien').eq.1) then
+                  write(*,*) 'phyisoetat0 16367: i,k,q_ancien(i,k)=',i,k,q_ancien(i,k) 
+                  write(*,*) 'xt_ancien(:,i,k)=',xt_ancien(:,i,k)
+                  stop
+              endif !  if (iso_verif_o18_aberrant_nostop 
+            endif !if (q_seri(i,k).gt.errmax) then  
+        endif !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then 
 #endif 
 
@@ -16543 +16450 @@
    CHARACTER(LEN=2) :: str2
    CHARACTER(LEN=5) :: str5
-   CHARACTER(LEN=maxlen) :: outiso, oldIso, modname, nam(2)
+   CHARACTER(LEN=maxlen) :: outiso, oldIso, modname, nam(3), oldIso2
    REAL :: xmin, xmax
    LOGICAL :: found
@@ -16559 +16466 @@
       outiso = isoName(ixt)
       oldIso = strTail(new2oldH2O(outiso), '_')            !--- Remove "H2O_" from "H2O_<iso>[_<tag>]"
+      i = INDEX(outiso, '_', .TRUE.)
+      oldIso2 = outiso(1:i-1)//outiso(i+1:LEN_TRIM(outiso)) ! CR 2023: on ajoute cette possibilité aussi, elle correspond au cas le plus récent.
+!      write(*,*) 'tmp 16541:'
+!      write(*,*) 'outiso=',outiso
+!      write(*,*) 'oldIso=',oldIso
+!      write(*,*) 'oldIso2=',oldIso2
+
       ! on lit seulement si ixt<=niso ou si on initialise les traceurs d'après fichier:
 #ifdef ISOTRAC
       IF(ixt <= niso .OR. initialisation_isotrac == 0) THEN
 #endif
-      found = phyetat0iso_srf3(xtsnow,      "XTSNOW", "Surface snow", 0.)
-      if (.NOT.found) CALL abort_physic('isotopes_routines_mod', 'phyiso_etat0_fichier 16581: unfound isotopic variable',1)
       found = phyetat0iso_srf3(fxtevap,     "XTEVAP", "evaporation",  0.)
+      if (.NOT.found) CALL abort_physic('isotopes_routines_mod', 'phyiso_etat0_fichier 16581a: unfound isotopic variable',1)
       found = phyetat0iso_get2(xtrain_fall, "xtrain_f", "xrain fall", 0.)
       found = phyetat0iso_get2(xtsnow_fall, "xtsnow_f", "xsnow fall", 0.)
@@ -16571 +16484 @@
       found = phyetat0iso_get3(xtl_ancien,  "XTLANCIEN", "QLANCIEN",  0.)
       found = phyetat0iso_get3(xts_ancien,  "XTSANCIEN", "QSANCIEN",  0.)
-      found = phyetat0iso_get2(xtrun_off_lic_0, "XTRUNOFFLIC0", "RUNOFFLIC0", 0.)
       found = phyetat0iso_get3(wake_deltaxt,  "WAKE_DELTAXT", "Delta hum. wake/env",  0.)
 #ifdef ISOVERIF
@@ -16580 +16492 @@
             DO nsrf = 1, nbsrf
                CALL iso_verif_egalite(fxtevap(iso_eau,i,nsrf),fevap(i,nsrf),TRIM(modname)//' 231c')
-               CALL iso_verif_egalite( xtsnow(iso_eau,i,nsrf), snow(i,nsrf),TRIM(modname)//' 231d')
             END DO
          END DO
@@ -16592 +16503 @@
          END DO
       END IF
-      IF(iso_eau > 0 .AND. ixt == iso_eau) THEN
-         DO i=1,klon
-            IF(iso_verif_egalite_nostop(run_off_lic_0(i),xtrun_off_lic_0(iso_eau,i),TRIM(modname)//' 326') == 1) THEN
-               WRITE(*,*) 'i=',i
-               STOP
-            END IF
-         END DO
-      END IF
 #endif
       ! ces variables n'ont pas de traceurs:
       IF(ixt <= niso) THEN
          found = phyetat0iso_get2(xtsol, "XTSOL", "Surface humidity / bucket", 0.)
+         if (.NOT.found) CALL abort_physic('isotopes_routines_mod', 'phyiso_etat0_fichier 16581b: unfound isotopic variable',1)
          found = phyetat0iso_get2(Rland_ice, "Rland_ice", "SR land ice", 0.)
-#ifdef ISOVERIF
-
+         found = phyetat0iso_srf3(xtsnow,      "XTSNOW", "Surface snow", 0.) ! CR avril 2023: deplacer ici 
+         found = phyetat0iso_get2(xtrun_off_lic_0, "XTRUNOFFLIC0", "RUNOFFLIC0", 0.)
+#ifdef ISOVERIF
          DO i=1,klon
             IF(iso_verif_noNaN_nostop(xtsol(ixt,i),TRIM(modname)//' 95') == 1) THEN
@@ -16612 +16517 @@
                STOP
             END IF
-         END DO
+            IF(ixt == iso_eau .AND. iso_eau > 0) THEN
+             DO nsrf = 1, nbsrf
+               CALL iso_verif_egalite(fxtevap(iso_eau,i,nsrf),fevap(i,nsrf),TRIM(modname)//' 231c')
+               CALL iso_verif_egalite( xtsnow(iso_eau,i,nsrf), snow(i,nsrf),TRIM(modname)//' 231d')
+             END DO
+             CALL iso_verif_egalite( xtrun_off_lic_0(iso_eau,i), run_off_lic_0(i),TRIM(modname)//' 231e')
+            ENDIF !IF(ixt == iso_eau .AND. iso_eau > 0) THEN
+         END DO !DO i=1,klon
 #endif
       END IF
@@ -16680 +16592 @@
   nam(1) = TRIM(pref)//TRIM(outiso)
   nam(2) = TRIM(pref)//TRIM(oldIso)
+  nam(3) = TRIM(pref)//TRIM(oldIso2)
   lFound = phyetat0_get(iso_tmp, nam, descr, default)
   field(ixt,:) = iso_tmp
@@ -16692 +16605 @@
   nam(1) = TRIM(pref)//TRIM(outiso)
   nam(2) = TRIM(pref)//TRIM(oldIso)
+  nam(3) = TRIM(pref)//TRIM(oldIso2)
   lFound = phyetat0_get(iso_tmp_lonlev, nam, descr, default)
   field(ixt,:,:) = iso_tmp_lonlev(:,:)
@@ -16703 +16617 @@
   nam(1) = TRIM(pref)//TRIM(outiso)
   nam(2) = TRIM(pref)//TRIM(oldIso)
+  nam(3) = TRIM(pref)//TRIM(oldIso2)
   lFound = phyetat0_srf(iso_tmp_lonsrf, nam, descr, default)
   field(ixt,:,:) = iso_tmp_lonsrf
@@ -18129 +18044 @@
 
       do iessai=1,nessai
-         day_nucl(iessai)   = 0.
-         month_nucl(iessai) = 0.
-         year_nucl(iessai)  = 0. 
+         day_nucl(iessai)   = 0
+         month_nucl(iessai) = 0
+         year_nucl(iessai)  = 0
          lat_nucl(iessai)   = 0.
          lon_nucl(iessai)   = 0.
@@ -18816 +18731 @@
           do i=1,n
             qttrac(i)=xt(ieau,i)
-!            if (qt(i).gt.0.0) then ! modif C Risi juillet 2020
-            if ((qt(i).gt.0.0).and.(xt(ieau,i).gt.0.0)) then
-               zcondtrac(i)=(zcond(i)/qt(i))*xt(ieau,i)
+            if (qt(i).gt.0.0) then ! modif C Risi juillet 2020 ! remodif Camille 9 mars 2023
+!            if ((qt(i).gt.0.0).and.(xt(ieau,i).gt.0.0)) then
+               zcondtrac(i)=(zcond(i)/qt(i))*qttrac(i)
             else !if (qt(i).eq.0) then
 #ifdef ISOVERIF              
@@ -18836 +18751 @@
             endif
             if (iso_HDO.gt.0) then
-!              if (qttrac(i).gt.ridicule_trac) then
                 call iso_verif_aberrant_choix(xttrac(iso_HDO,i), &
      &           qttrac(i),ridicule_trac,deltalimtrac, &
      &           'condisotrac 205')
-!              endif
             endif
             call iso_verif_positif(qt(i)-cond(i), &
@@ -18853 +18766 @@
           call condiso_liq_ice_vectall(xttrac,qttrac,zcondtrac, &
      &           tcond,zfice,zxticetrac,zxtliqtrac,n)
-#ifdef ISOVERIF          
-          write(*,*) 'condisotrac 167: après condiso'
-#endif           
+
           do i=1,n
           do iiso=1,niso          
@@ -18881 +18792 @@
 #endif
 
+subroutine iso_init_ideal(q,xt,ixt,alpha,kcin,toce)
+
+        USE isotopes_mod, ONLY: iso_eau,iso_HDO,ridicule
+#ifdef ISOVERIF
+        USE isotopes_verif_mod
+#endif
+        implicit none
+
+        ! inputs
+        real q ! humidité spec
+        integer ixt ! indice isotopique
+        real alpha ! coef frac à l'eq
+        real kcin ! coef frac cinétique
+        real toce ! rapport iso ds ocean surface
+
+        ! outputs
+        real xt ! equivalent iso de l'humidité spec, même unité.
+
+        ! locals
+        real RMerlivat
+        real q0,h0 ! conditions initiales de la distill de Rayleigh
+        parameter (q0=20e-3,h0=0.7)
+
+        ! verifier que ixt est un isotope et pas un tagging
+        if (ixt.gt.niso) then
+          CALL abort_physic('isotopes_routines_mod', 'iso_init_ideal, ixt>niso', 1)
+        endif
+
+        ! R selon Merlivat:
+        RMerlivat=toce/alpha *(1.0-kcin)/(1.0-kcin*h0)
+
+        ! R d'après Rayleigh
+        xt=q*RMerlivat*(min(q0,q)/q0)**(alpha-1.0)
+
+#ifdef ISOVERIF
+        call iso_verif_noNaN(xt, 'isotopes_routines_mod 18930a: iso_init_ideal')
+        if ((iso_HDO.gt.0).and.(ixt.eq.iso_HDO)) then
+            if (q.gt.ridicule) then
+                call iso_verif_aberrant_encadre( xt/q, 'isotopes_routines_mod 18930b: iso_init_ideal')
+            endif
+        endif
+        if ((iso_eau.gt.0).and.(ixt.eq.iso_eau)) then
+             call iso_verif_egalite(xt,q, 'isotopes_routines_mod 18930c: iso_init_ideal')
+        endif
+#endif
+        
+
+end subroutine iso_init_ideal
+
+
+subroutine appel_stewart_debug(lwork,nloc,inb,na,i, &
+                evap,water,rpprec,rr,wdtrain, &
+                xtevap,xtwater,xtp,xt,xtwdtrain)
+USE isotopes_mod, ONLY: iso_eau, iso_HDO,thumxt1, &
+&       bidouille_anti_divergence,ridicule,Rdefault
+use infotrac_phy, ONLY: ntraciso=>ntiso, niso
+#ifdef ISOTRAC
+    use isotrac_mod, only: option_cond,izone_cond,index_iso,index_zone,izone_poubelle
+#endif
+#ifdef ISOVERIF
+USE isotopes_verif_mod
+#endif
+implicit none
+
+
+! inputs
+integer nloc,na,i ! dimension horiz effective
+logical lwork(nloc)
+real wdtrain(nloc),xtwdtrain(ntraciso,nloc)
+real xt(ntraciso,nloc,na)
+real evap(nloc,na),water(nloc,na),rpprec(nloc,na),rr(nloc,na)
+integer inb(nloc)
+
+! outputs
+real xtevap(ntraciso,nloc,na),xtwater(ntraciso,nloc,na),xtp(ntraciso,nloc,na)
+
+! locals
+integer il,ixt
+
+     do il=1,nloc
+       if (i.le.inb(il) .and. lwork(il)) then
+          if (wdtrain(il).gt.0.) then
+            do ixt=1,ntraciso
+             xtwater(ixt,il,i)= xtwdtrain(ixt,il)/wdtrain(il)*water(il,i)
+             xtevap(ixt,il,i)= xtwdtrain(ixt,il)/wdtrain(il)*evap(il,i)
+            enddo
+          else !if (wdtrain(il).gt.0.) then
+            do ixt=1,niso
+             xtwater(ixt,il,i)= Rdefault(ixt)*water(il,i)
+             xtevap(ixt,il,i)= Rdefault(ixt)*evap(il,i)
+            enddo
+#ifdef ISOTRAC
+            do ixt=1+niso,ntraciso
+             if (index_zone(ixt).eq.izone_poubelle) then
+               xtwater(ixt,il,i)= Rdefault(index_iso(ixt))*water(il,i)
+               xtevap(ixt,il,i)= Rdefault(index_iso(ixt))*evap(il,i)
+             else
+               xtwater(ixt,il,i)= 0.
+               xtevap(ixt,il,i)=0.
+             endif
+            enddo ! do ixt=1+niso,ntraciso
+#endif 
+         endif !if (wdtrain(il).gt.0.) then
+         do ixt=1,ntraciso
+             xtp(ixt,il,i)= xt(ixt,il,i)/rr(il,i)*rpprec(il,i)
+         enddo !do ixt=1,ntraciso
+      endif
+    enddo ! do il=1,ncum
+end subroutine appel_stewart_debug
+
 END MODULE isotopes_routines_mod
 #endif

LMDZ6/trunk/libf/phylmdiso/isotopes_verif_mod.F90

@@ -2028 +2028 @@
 
           if (iso_verif_egalite_choix_nostop(xtractot,x(iiso), &
-     &        err_msg//', verif trac egalite, iso '// &
+     &        err_msg//', verif trac egalite1, iso '// &
      &        TRIM(isoName(iiso)), &
      &        errmaxin,errmaxrelin).eq.1) then
             write(*,*) 'iso_verif_traceur 202: x=',x
 !            write(*,*) 'xtractot=',xtractot
+            do izone=1,nzone  
+              ixt=itZonIso(izone,iiso)
+              write(*,*) 'izone,iiso,ixt=',izone,iiso,ixt
+            enddo
             iso_verif_tracm_choix_nostop=1
           endif
@@ -2398 +2402 @@
         call iso_verif_egalite_std_vect( &
      &           xtractot,xiiso, &
-     &           err_msg//', verif trac egalite, iso ' &
+     &           err_msg//', verif trac egalite2, iso ' &
      &           //TRIM(isoName(iiso)), &
      &           n,m,errmax,errmaxrel)

LMDZ6/trunk/libf/phylmdiso/isotrac_routines_mod.F90

@@ -508 +508 @@
 !          write(*,*) 'compress 910: xtp_avantevap(iso_eau,cas(1))=',
 !     :           xtp_avantevap(iso_eau,cas(1))
+!        write(*,*) 'compress_evap_liq_zone 510: ncas,ncum=',ncas,ncum
+        ptrac(:)=0. ! CR 31 mars 2023: initialisation de ptrac
 
           ieau=index_trac(izone,iso_eau)
@@ -519 +521 @@
             else
 #ifdef ISOVERIF 
-                call iso_verif_egalite(( &
+              call iso_verif_egalite(( &
      &                   Eqi_prime(cas(il))),0.0, &
      &                   'compress_stewart 979')
@@ -685 +687 @@
                 call iso_verif_egalite(( &
      &                   Eqi_prime(cas(il))),0.0, &
-     &                   'compress_stewart 979')
+     &                   'compress_stewart 979b')
 #endif                
               Eqi_prime_cas(il)=0.0 
@@ -870 +872 @@
             enddo ! do iiso=1,niso           
           endif !if ((Eqi(il)*fac_ftmr(il).gt.evap_franche).and.
-#ifdef ISOVERIF
-          if (il.eq.9) then            
-            ixt=index_trac(2,2)
-            write(*,*) 'il,ixt,xtrevap_tag(ixt,il)=', &
-     &           il,xtrevap_tag(ixt,il)
-          endif !if (il.eq.9) then
-#endif          
          enddo !do il=1,ncas
 !         write(*,*) 'compress_stewart 1453 tmp: zxt=',
@@ -2791 +2786 @@
        real xtractot
         
-        write(*,*) 'iso_verif_traceur tmp 822: xt(:,66,1)=',x(:,66,1)
         ! verif noNaN
         call iso_verif_traceur_noNaN_vect(x,n,m,err_msg)

LMDZ6/trunk/libf/phylmdiso/pbl_surface_mod.F90

@@ -2407 +2407 @@
             call iso_verif_noNaN(yxtevap(ixt,j), &
          &      'pbl_surface 1056a: apres surf_land')
+          enddo
+          do ixt=1,niso
             call iso_verif_noNaN(yxtsol(ixt,j), &
          &      'pbl_surface 1056b: apres surf_land')
@@ -2475 +2477 @@
                    call iso_verif_noNaN(yxtevap(ixt,j), &
                         &      'pbl_surface 1095a: apres surf_landice')
+                enddo
+                do ixt=1,niso
                    call iso_verif_noNaN(yxtsol(ixt,j), &
                         &      'pbl_surface 1095b: apres surf_landice')
@@ -2570 +2574 @@
             call iso_verif_noNaN(yxtevap(ixt,j), &
          &      'pbl_surface 1165a: apres surf_seaice')
+          enddo
+          do ixt=1,niso
             call iso_verif_noNaN(yxtsol(ixt,j), &
          &      'pbl_surface 1165b: apres surf_seaice')

LMDZ6/trunk/libf/phylmdiso/phyredem.F90

@@ -581 +581 @@
       CALL put_field_srf1(pass, "XTEVAP"//TRIM(outiso), "Evaporation de surface",iso_tmp_lonsrf)
 
-      iso_tmp_lonsrf(:,:)=xtsnow(ixt,:,:)
-      CALL put_field_srf1(pass, "XTSNOW"//TRIM(outiso), "NEIGE",       iso_tmp_lonsrf)       
-
       iso_tmp(:)=xtrain_fall(ixt,:)
       CALL put_field(pass,    "xtrain_f"//TRIM(outiso), "precipitation liquide",iso_tmp)
@@ -602 +599 @@
       CALL put_field(pass,"WAKE_DELTAXT"//TRIM(outiso), "WAKE_DELTAQ", iso_tmp_lonlev)
 
+      iso_tmp_lonlev(:,:)=wake_deltaxt(ixt,:,:)
+      CALL put_field(pass,"WAKE_DELTAXT"//TRIM(outiso), "WAKE_DELTAXT",iso_tmp_lonlev)
+
+      ! variables seulement pour niso:
+      if (ixt.le.niso) then
+
+      iso_tmp_lonsrf(:,:)=xtsnow(ixt,:,:)
+      CALL put_field_srf1(pass, "XTSNOW"//TRIM(outiso), "NEIGE",       iso_tmp_lonsrf)
+
+      iso_tmp(:)=xtsol(ixt,:)
+      CALL put_field(pass,      "XTSOL"//TRIM(outiso), "Eau dans le sol (mm)",iso_tmp)
+
+      iso_tmp(:)=Rland_ice(ixt,:)
+      CALL put_field(pass,  "Rland_ice"//TRIM(outiso), "ratio land ice",      iso_tmp)
+
       iso_tmp(:)=xtrun_off_lic_0(ixt,:)
       CALL put_field(pass,"XTRUNOFFLIC0"//TRIM(outiso), "Runofflic0",  iso_tmp)
 
-      iso_tmp_lonlev(:,:)=wake_deltaxt(ixt,:,:)
-      CALL put_field(pass,"WAKE_DELTAXT"//TRIM(outiso), "WAKE_DELTAXT",iso_tmp_lonlev)
-
-      ! variables seulement pour niso:
-      if (ixt.le.niso) then
-
-      iso_tmp(:)=xtsol(ixt,:)
-      CALL put_field(pass,      "XTSOL"//TRIM(outiso), "Eau dans le sol (mm)",iso_tmp)
-
-      iso_tmp(:)=Rland_ice(ixt,:)
-      CALL put_field(pass,  "Rland_ice"//TRIM(outiso), "ratio land ice",      iso_tmp)
-
       endif ! if (ixt.le.niso) then
 

LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90

@@ -157 +157 @@
         & iso_verif_aberrant_choix,iso_verif_positif, &
         & iso_verif_positif_choix_vect,iso_verif_o18_aberrant_nostop, &
-        & iso_verif_init,&
+        & iso_verif_init,iso_verif_aberrant_enc_choix_nostop,&
         & iso_verif_positif_strict_nostop,iso_verif_O18_aberrant_enc_vect2D
 #endif
@@ -2480 +2480 @@
 #ifdef ISO
 #ifdef ISOVERIF
-!    write(*,*) 'physiq 1847: qx(1,1,:)=',qx(1,1,:)
+    write(*,*) 'physiq 1847: qx(1,1,:)=',qx(1,1,:)
+    write(*,*) 'iqIsoPha(:,ivap)=',iqIsoPha(:,ivap)
     write(*,*) 'physiq 1846b: ok_isotopes,ntraciso,niso=',niso>0,ntraciso,niso
 #endif
@@ -2596 +2597 @@
      &              'physiq 2099 ql').eq.1) then
                   write(*,*) 'i,k,ql_seri(i,k)=',i,k,ql_seri(i,k) 
-                  stop
+                  !stop
               endif !  if (iso_verif_o18_aberrant_nostop 
             endif !if (q_seri(i,k).gt.errmax) then  
@@ -2605 +2606 @@
      &              'physiq 2099 qs').eq.1) then
                   write(*,*) 'i,k,qs_seri(i,k)=',i,k,qs_seri(i,k) 
-                  stop
+                  !stop
               endif !  if (iso_verif_o18_aberrant_nostop 
             endif !if (q_seri(i,k).gt.errmax) then
           enddo !k=1,klev
          enddo  !i=1,klon
-        endif !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then
+        endif !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then         
+#ifdef ISOTRAC  
+         DO k = 1, klev
+           DO i = 1, klon    
+             call iso_verif_traceur(xt_seri(1,i,k),'physiq 2620a')
+             call iso_verif_traceur(xtl_seri(1,i,k),'physiq 2620b')
+             call iso_verif_traceur(xts_seri(1,i,k),'physiq 2620c') 
+           enddo
+         enddo
+#endif
 #endif
     !
@@ -2844 +2854 @@
        ENDIF
     ENDIF
+
+#ifdef ISO
+#ifdef ISOVERIF 
+#ifdef ISOTRAC 
+         DO k = 1, klev
+           DO i = 1, klon         
+             call iso_verif_traceur(xt_seri(1,i,k), &
+     &           'physiq 2856: avant reevp')
+           enddo
+         enddo
+#endif
+#endif
+#endif
     !
     ! Re-evaporer l'eau liquide nuageuse
@@ -3458 +3481 @@
                   call iso_verif_egalite(q_w(i,k),xt_w(iso_eau,i,k),'physiq 3338')
                 endif
+                if (iso_HDO.gt.0) then
+                  if ((iso_verif_aberrant_enc_choix_nostop(xt_x(iso_hdo,i,k),q_x(i,k), &
+                        ridicule,deltalim,'physic 3462a xt_x').eq.1).or. &
+                        (iso_verif_aberrant_enc_choix_nostop(xt_w(iso_hdo,i,k),q_w(i,k), &
+                        ridicule,deltalim,'physic 3462b xt_w').eq.1)) then
+                     write(*,*) 'i,k=',i,k
+                     write(*,*) 'q_x(i,k),q_seri(i,k),wake_s(i),wake_deltaq(i,k)=', &
+                                q_x(i,k),q_seri(i,k),wake_s(i),wake_deltaq(i,k)
+                     write(*,*) 'xt_x(iso_hdo,i,k),xt_seri(iso_hdo,i,k),wake_s(i),wake_deltaxt(iso_hdo,i,k)=', &
+                                xt_x(iso_hdo,i,k),xt_seri(iso_hdo,i,k),wake_s(i),wake_deltaxt(iso_hdo,i,k)
+                     write(*,*) 'deltaD_seri,wake=',deltaD(xt_seri(iso_hdo,i,k)/q_seri(i,k)), &
+                                deltaD(wake_deltaxt(iso_hdo,i,k)/wake_deltaq(i,k))
+                     write(*,*) 'deltaD_x,deltaD_w=',deltaD(xt_x(iso_hdo,i,k)/q_x(i,k)),deltaD(xt_w(iso_hdo,i,k)/q_w(i,k))
+                     stop
+                  endif
+                endif
 #endif
 #endif
@@ -3590 +3629 @@
               call iso_verif_aberrant_encadre( &
      &           xt_w(iso_hdo,i,k)/q_w(i,k), &
-     &          'physic 2657b')
+     &          'physic 2657c')
              endif !if (q_x(i,k).gt.ridicule) then
            endif !if (iso_HDO.gt.0) then
@@ -4255 +4294 @@
              endif !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then    
              if ((iso_O17.gt.0).and.(iso_O18.gt.0)) then
-               if ((q_seri(i,k).gt.ridicule).and.(l.lt.nlevmaxO17)) then
+               if ((q_seri(i,k).gt.ridicule).and.(k.lt.nlevmaxO17)) then
                  call iso_verif_aberrant_o17(xt_seri(iso_o17,i,k) &
      &           /q_seri(i,k),xt_seri(iso_o18,i,k) &
@@ -4619 +4658 @@
         do i=1,klon
            if ((iso_O17.gt.0).and.(iso_O18.gt.0)) then
-            if ((q_seri(i,k).gt.ridicule).and.(l.lt.nlevmaxO17)) then
+            if ((q_seri(i,k).gt.ridicule).and.(k.lt.nlevmaxO17)) then
              call iso_verif_aberrant_o17(xt_seri(iso_o17,i,k) &
      &           /q_seri(i,k),xt_seri(iso_o18,i,k) &
@@ -4841 +4880 @@
             endif !if (q_seri(i,k).gt.errmax) then
             if (ql_seri(i,k).gt.ridicule) then
-               call iso_verif_aberrant( &
-     &           xtl_seri(iso_HDO,i,k)/ql_seri(i,k),'physiq 2871')
+               call iso_verif_aberrant_choix(xtl_seri(iso_HDO,i,k),ql_seri(i,k), &
+                        ridicule,deltalim_snow,'physiq 2871')
                if (iso_O18.gt.0) then  
                  if (iso_verif_o18_aberrant_nostop( &
@@ -4854 +4893 @@
                         write(*,*) 'deltaO(d_ql_lsc(i,k))=',deltaO( &
      &                           d_xtl_lsc(iso_O18,i,k)/d_ql_lsc(i,k))
-                        stop
+                        !stop
                  endif
                endif ! if (iso_O18.gt.0) then  
@@ -4864 +4903 @@
         do i=1,klon 
           do k=1,nlev
-           if ((q_seri(i,k).gt.ridicule).and.(l.lt.nlevmaxO17)) then
+           if ((q_seri(i,k).gt.ridicule).and.(k.lt.nlevmaxO17)) then
             call iso_verif_aberrant_o17(xt_seri(iso_o17,i,k) &
      &           /q_seri(i,k),xt_seri(iso_o18,i,k) &
@@ -6149 +6188 @@
      &              xt_seri(iso_HDO,i,k)/q_seri(i,k), &
      &              xt_seri(iso_O18,i,k)/q_seri(i,k), &
-     &              'physiq 5937, juste apres methox, qv').eq.1) then
+     &              'physiq 5937a, juste apres methox, qv').eq.1) then
                   write(*,*) 'physic 2444: i,k,q_seri(i,k)=',i,k,q_seri(i,k) 
                   write(*,*) 'xt_seri(:,i,k)=',xt_seri(:,i,k)
@@ -6159 +6198 @@
      &              xtl_seri(iso_HDO,i,k)/ql_seri(i,k), &
      &              xtl_seri(iso_O18,i,k)/ql_seri(i,k), &
-     &              'physiq 5937, juste apres methox, ql').eq.1) then
+     &              'physiq 5937b, juste apres methox, ql').eq.1) then
                   write(*,*) 'i,k,ql_seri(i,k)=',i,k,ql_seri(i,k) 
-                  stop
+                  !stop
               endif !  if (iso_verif_o18_aberrant_nostop 
             endif !if (q_seri(i,k).gt.errmax) then  
@@ -6168 +6207 @@
      &              xts_seri(iso_HDO,i,k)/qs_seri(i,k), &
      &              xts_seri(iso_O18,i,k)/qs_seri(i,k), &
-     &              'physiq 5937, juste apres methox, qs').eq.1) then
+     &              'physiq 5937c, juste apres methox, qs').eq.1) then
                   write(*,*) 'i,k,qs_seri(i,k)=',i,k,qs_seri(i,k) 
-                  stop
+                  !stop
               endif !  if (iso_verif_o18_aberrant_nostop 
             endif !if (q_seri(i,k).gt.errmax) then

LMDZ6/trunk/libf/phylmdiso/reevap.F90

@@ -85 +85 @@
       do ixt=1,ntiso
         call iso_verif_noNaN(xt_seri(ixt,i,k), &
-     &     'physiq 2417: apres evap tot')
+     &     'reevap 2417: apres evap tot')
       enddo
       if (iso_eau.gt.0) then
               call iso_verif_egalite_choix( &
      &           xt_seri(iso_eau,i,k),q_seri(i,k), &
-     &          'physiq 1891+, après reevap totale',errmax,errmaxrel)
+     &          'reevap 1891+, après reevap totale',errmax,errmaxrel)
               call iso_verif_egalite_choix( &
      &           xtl_seri(iso_eau,i,k),ql_seri(i,k), &
-     &          'physiq 2209+, après reevap totale',errmax,errmaxrel)
+     &          'reevap 2209+, après reevap totale',errmax,errmaxrel)
        endif !if (iso_eau.gt.0) then       
       if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then  
@@ -100 +100 @@
      &              xt_seri(iso_HDO,i,k)/q_seri(i,k), &
      &              xt_seri(iso_O18,i,k)/q_seri(i,k), &
-     &              'physiq 2315: apres reevap totale').eq.1) then
+     &              'reevap 2315: apres reevap totale').eq.1) then
                   write(*,*) 'i,k,q_seri(i,k)=',i,k,q_seri(i,k)
                   write(*,*) 'd_q_eva(i,k)=',d_q_eva(i,k)
@@ -111 +111 @@
 #ifdef ISOTRAC      
              call iso_verif_traceur(xt_seri(1,i,k), &
-     &           'physiq 2165')
+     &           'reevap 2165a')
              call iso_verif_traceur_pbidouille(xt_seri(1,i,k), &
-     &           'physiq 2165b')
+     &           'reevap 2165b')
 #endif               
          
@@ -146 +146 @@
       do ixt=1,ntiso
       call iso_verif_noNaN(xt_seri(ixt,i,k), &
-     &     'physiq 2417: apres evap tot')
+     &     'reevap 2417: apres evap tot')
       enddo
       if (iso_eau.gt.0) then
               call iso_verif_egalite_choix( &
      &           xt_seri(iso_eau,i,k),q_seri(i,k), &
-     &          'physiq 1891, après réévap totale',errmax,errmaxrel)
+     &          'reevap 1891, après réévap totale',errmax,errmaxrel)
               call iso_verif_egalite_choix( &
      &           xtl_seri(iso_eau,i,k),ql_seri(i,k), &
-     &          'physiq 2209, après réévap totale',errmax,errmaxrel)
+     &          'reevap 2209, après réévap totale',errmax,errmaxrel)
               call iso_verif_egalite_choix( &
      &           xts_seri(iso_eau,i,k),qs_seri(i,k), &
-     &          'physiq 2209b, après réévap totale',errmax,errmaxrel)
+     &          'reevap 2209b, après réévap totale',errmax,errmaxrel)
        endif !if (iso_eau.gt.0) then
      
@@ -165 +165 @@
      &              xt_seri(iso_HDO,i,k)/q_seri(i,k), &
      &              xt_seri(iso_O18,i,k)/q_seri(i,k), &
-     &              'physiq 2408: apres reevap totale').eq.1) then
+     &              'reevap 2408: apres reevap totale').eq.1) then
                   write(*,*) 'i,k,q_seri(i,k)=',i,k,q_seri(i,k)                       
                   stop
@@ -173 +173 @@
 #ifdef ISOTRAC      
              call iso_verif_traceur(xt_seri(1,i,k), &
-     &           'physiq 2165')
+     &           'reevap 2165c')
              call iso_verif_traceur_pbidouille(xt_seri(1,i,k), &
-     &           'physiq 2165b')
+     &           'reevap 2165d')
 #endif                  
 #endif                  

LMDZ6/trunk/libf/phylmdiso/wake.F90

@@ -38 +38 @@
   USE isotopes_verif_mod
 !, ONLY: errmax,errmaxrel
-  USE isotopes_mod, ONLY: iso_eau,iso_hdo
+  USE isotopes_mod, ONLY: iso_eau,iso_hdo,ridicule
 #endif
 #endif
@@ -714 +714 @@
           xtu(ixt,i,k)  =  xte(ixt,i,k) - deltaxtw(ixt,i,k)*sigmaw(i)
         enddo !do ixt=1,ntraciso
-#endif
-    END DO
-  END DO
-
-
-#ifdef ISO
 #ifdef ISOVERIF
-        ! TODO   
-#endif
-#endif
+        if (iso_eau.gt.0) then
+              call iso_verif_egalite(deltaqw(i,k),deltaxtw(iso_eau,i,k),'wake 723a')
+              call iso_verif_egalite(qu(i,k),xtu(iso_eau,i,k),'wake 723b')
+        endif
+        if (iso_HDO.gt.0) then
+               if (iso_verif_aberrant_enc_choix_nostop(xtu(iso_hdo,i,k),qu(i,k),ridicule,deltalim, &
+                        'wake 723c xtu').eq.1) then
+                     stop
+               endif
+         endif  
+#endif
+#endif
+    END DO
+  END DO
+
+
 
   DO k = 1, klev - 1
@@ -1837 +1844 @@
 #ifdef ISO
 #ifdef ISOVERIF
+      write(*,*) 'wake 1859'
       if (iso_eau.gt.0) then
         DO k= 1,klev