Ignore:
Timestamp:
Apr 3, 2023, 10:25:44 AM (15 months ago)
Author:
crisi
Message:

Bug corrections in LMDZiso, especially for water tagging

Location:
LMDZ6/trunk/libf/phylmdiso
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmdiso/add_phys_tend_mod.F90

    r4143 r4491  
    313313! Ajout des tendances sur le vent et l'eau liquide
    314314!======================================================================
     315
     316#ifdef ISO
     317#ifdef ISOVERIF
     318        if (iso_eau.gt.0) then
     319          call iso_verif_egalite_vect2D( &
     320                xt_seri,q_seri, &
     321                'add_phys_tend 321a: '//text,ntraciso,klon,klev)
     322          call iso_verif_egalite_vect2D( &
     323                xtl_seri,ql_seri, &
     324                'add_phys_tend 321b: '//text,ntraciso,klon,klev)
     325         endif !if (iso_eau.gt.0) then
     326#ifdef ISOTRAC
     327        call iso_verif_traceur_vect(xt_seri,klon,klev, &
     328     &           'add_phys_tend 328a: '//text)
     329        call iso_verif_traceur_vect(xtl_seri,klon,klev, &
     330     &           'add_phys_tend 328b: '//text)
     331#endif
     332#endif
     333#endif
    315334
    316335     u_seri(:,:)=u_seri(:,:)+zdu(:,:)
    … …  
    370389                'add_phys_tend 138: '//text,ntraciso,klon,klev)
    371390         endif !if (iso_eau.gt.0) then
     391#ifdef ISOTRAC
     392        call iso_verif_traceur_vect(xt_seri,klon,klev, &
     393     &           'add_phys_tend 374a: '//text)
     394        call iso_verif_traceur_vect(xtl_seri,klon,klev, &
     395     &           'add_phys_tend 374b: '//text)
     396#endif
    372397#endif
    373398#endif
    … …  
    495520                        'add_phys_tend 175'//text,ntraciso,klon,klev)
    496521      endif !if (iso_eau.gt.0) then
     522#ifdef ISOTRAC
     523        call iso_verif_traceur_vect(xt_seri,klon,klev, &
     524     &           'add_phys_tend 499'//text)
     525#endif
    497526#endif
    498527#endif

  • LMDZ6/trunk/libf/phylmdiso/cv30_routines.F90

    r4143 r4491  
    26472647     &          )
    26482648#ifdef ISO
    2649     use infotrac_phy, ONLY: ntraciso=>ntiso
     2649    use infotrac_phy, ONLY: ntraciso=>ntiso, niso
    26502650    use isotopes_mod, ONLY: essai_convergence, iso_eau,iso_HDO,ridicule
    2651     use isotopes_routines_mod, ONLY: appel_stewart_vectall
     2651    use isotopes_routines_mod, ONLY: appel_stewart_vectall,appel_stewart_debug
    26522652#ifdef ISOVERIF
    26532653    use isotopes_verif_mod, ONLY: errmax,errmaxrel, &
    … …  
    32083208       enddo
    32093209#endif
     3210
     3211        if (1.eq.0) then
    32103212        ! appel de appel_stewart_vectorise
    32113213        call appel_stewart_vectall(lwork,ncum, &
    … …  
    32183220     &          i,inb, & ! altitude: car cas particulier en INB
    32193221     &          na,nd,nloc,cvflag_grav,ginv,1e-16)
     3222
     3223        else !if (1.eq.0) then
     3224          ! truc simple sans fractionnement
     3225          ! juste pour debuggage
     3226          call appel_stewart_debug(lwork,nloc,inb,na,i, &
     3227                evap,water,rpprec,rr,wdtrain, &
     3228                xtevap,xtwater,xtp,xt,xtwdtrain)
     3229        endif ! if (1.eq.0) then
     3230
    32203231
    32213232#ifdef ISOVERIF

  • LMDZ6/trunk/libf/phylmdiso/cv3p_mixing.F90

    r4143 r4491  
    451451          call iso_verif_aberrant_choix( &
    452452                 xtelij(iso_HDO,il,i,j),elij(il,i,j), &
    453                  ridicule,deltalim,'cv3p_mixing 1993')
     453                 ridicule,deltalim_snow,'cv3p_mixing 1993')
    454454        endif !if (iso_hdo.gt.0) then
    455455#ifdef ISOTRAC 

  • LMDZ6/trunk/libf/phylmdiso/cv_driver.F90

    r4143 r4491  
    997997
    998998#ifdef ISO
    999        write(*,*) 'klev=',klev
     999!       write(*,*) 'klev=',klev
    10001000#ifdef ISOVERIF
    10011001       write(*,*) 'cv_driver 930: apres cv3_unsat'

  • LMDZ6/trunk/libf/phylmdiso/fisrtilp.F90

    r4380 r4491  
    509509        endif !if (iso_eau.gt.0) then
    510510        if (iso_HDO.gt.0) then 
    511          if (zrfl(i).gt.ridicule_rain) then
    512              call iso_verif_aberrant(zxtrfl(iso_HDO,i) &
    513      &           /zrfl(i),'il pleut 316')
    514          endif !if (zrfl(i).gt.ridicule_rain) then
     511             call iso_verif_aberrant_choix(zxtrfl(iso_HDO,i) &
     512     &           ,zrfl(i),ridicule_rain,deltalim_snow,'il pleut 316')
    515513        endif  !if (iso_HDO.gt.0) then 
    516514#ifdef ISOTRAC   
    … …  
    17321730        if (iso_HDO.gt.0) then
    17331731           do i=1,klon 
    1734               if (zcond(i).gt.ridicule) then
    1735                  call iso_verif_aberrant(zxtcond(iso_HDO,i) &
    1736      &                  /zcond(i), 'il pleut 637')               
    1737               endif  !if (zcond(i).gt.ridicule) then
    1738               IF ((t(i,1)) .LT. RTT) THEN
    17391732                  call iso_verif_aberrant_choix(zxtcond(iso_hdo,i), &
    17401733     &           zcond(i),ridicule_rain,deltalim_snow,'il pleut 1276')
    1741               endif
    17421734           enddo ! do i=1,klon 
    17431735        endif !if ((iso_eau.gt.0).and.(ixt.eq.iso_eau)) then
    … …  
    20312023           if (zoliq(i).gt.ridicule) then
    20322024            if (iso_HDO.gt.0) then 
    2033                 call iso_verif_aberrant(zxtoliq(iso_HDO,i)/zoliq(i), &
    2034      &           'il pleut 895a')
     2025                ! Camille 9 mars 2023: on est moins stricte avec le condensat
     2026                call iso_verif_aberrant_choix(zxtoliq(iso_HDO,i),zoliq(i), &
     2027     &           ridicule_rain,deltalim_snow, 'il pleut 895a')
    20352028              if (iso_O18.gt.0) then
    20362029                if (iso_verif_O18_aberrant_nostop(zxtoliq(iso_HDO,i)/zoliq(i), &
    20372030     &           zxtoliq(iso_O18,i)/zoliq(i),'il pleut 895b').eq.1) then
    20382031                   write(*,*) 'i,k,zoliq,zfice=',i,k,zoliq(i),zfice(i)
    2039                    stop
     2032                   !stop
    20402033                endif
    20412034              endif ! if (iso_HDO.gt.0) then
    … …  
    21712164          endif !if (iso_eau.gt.0) then
    21722165          if (iso_HDO.gt.0) then
    2173              if (zoliq(i).gt.ridicule) then
    2174                 call iso_verif_aberrant(zxtoliq(iso_HDO,i)/zoliq(i), &
    2175      &           'il pleut 963')
    2176              endif
     2166                call iso_verif_aberrant_choix(zxtoliq(iso_HDO,i),zoliq(i), &
     2167     &           ridicule_rain, deltalim_snow, 'il pleut 963')
    21772168          endif !if (iso_HDO.gt.0) then
    21782169        ! end cam verif
    … …  
    24822473            call iso_verif_aberrant(pxtrainfl(iso_HDO,i,k)/prfl(i,k), &
    24832474     &            'il pleut 1020')
    2484           endif !if (prfl(i,k).gt.ridicule_rain) then       
    2485           if (psfl(i,k).gt.ridicule_rain) then
    2486             call iso_verif_aberrant(pxtsnowfl(iso_HDO,i,k)/psfl(i,k), &
    2487      &           'il pleut 1024')
    2488           endif !if (psfl(i,k).gt.ridicule_rain) then
     2475          endif !if (prfl(i,k).gt.ridicule_rain) then           
     2476          call iso_verif_aberrant_choix(pxtsnowfl(iso_HDO,i,k),psfl(i,k), &
     2477     &           ridicule_rain, deltalim_snow, 'il pleut 1024')
    24892478          if (zq(i).gt.ridicule) then
    24902479              call iso_verif_aberrant_encadre(zxt(iso_HDO,i)/zq(i), &

  • LMDZ6/trunk/libf/phylmdiso/isotopes_mod.F90

    r4399 r4491  
    302302   IF(getKey('Rdefault',Rdefault,isoName)) CALL abort_physic(modname, 'can''t get Rdefault',1)
    303303   IF(getKey('Rmethox', Rmethox, isoName)) CALL abort_physic(modname, 'can''t get Rmethox', 1)
    304    IF(.NOT.Rdefault_smow) Rdefault(:) = 0.0
     304
     305   IF(.NOT.Rdefault_smow) then
     306        Rdefault(:) = 0.0
     307        if (iso_eau.gt.0) Rdefault(iso_eau) = 1.0 ! correction Camille 30 mars 2023
     308   ENDIF
     309   write(*,*) 'Rdefault=',Rdefault
    305310
    306311   !--- Sensitivity test: no kinetic effect in sfc evaporation

  • LMDZ6/trunk/libf/phylmdiso/isotopes_routines_mod.F90

    r4452 r4491  
    154154&                   'revap_ilp 131')
    155155       enddo
     156#ifdef ISOTRAC
     157     call iso_verif_traceur(zxtrfl_ancien(1,i), &
     158&           'iso_revap_fisrtilp 158: debut')
     159#endif
    156160#endif
    157161    endif !if (zrfln(i).gt.ridicule*1e-2) then 
    … …  
    193197        endif !if (zq(i).gt.ridicule) then
    194198     endif !if ((iso_HDO.gt.0.and.(iso_O18.gt.0) then
     199#ifdef ISOTRAC
     200     call iso_verif_traceur(zxtrfl_ancien(1,i), &
     201&           'iso_revap_fisrtilp 201: debut quand pas de precip')
     202#endif
    195203!     write(*,*) 'iso_routines tmp 184'
    196204#endif               
    … …  
    504512!             stop
    505513  endif
    506   enddo
     514  enddo !do iiso=1,niso
    507515enddo !do i=1,ncas_evap_liq
    508516#endif         
    … …  
    514522&       izone,zqevfl(1),Exi(1,1),fac_fluxtomixratio(1), &
    515523&       xtrevap_tag(1,1),1,hdiag(1))
     524        ! dans cette routine, zxtrfl reçoit zxtrfln_cas
    516525
    517526  enddo !do izone=1,ntraceurs_zone
     527
    518528#ifdef ISOVERIF
    519529do i=1,ncas_evap_liq
    … …  
    522532&              0.0,'revap_ilp 414')
    523533enddo
     534     call iso_verif_traceur(zxtrfl(1,cas_evap_liq(i)), &
     535&           'iso_revap_fisrtilp 470a: apres stewart_explicite_vectall')
    524536enddo !do i=1,ncas
    525537#endif           
    … …  
    536548     call iso_verif_traceur(zxt(1,cas_evap_liq(i)), &
    537549&           'iso_revap_fisrtilp 282')
     550     call iso_verif_traceur(zxtrfl(1,cas_evap_liq(i)), &
     551&           'iso_revap_fisrtilp 804a')
    538552     call iso_verif_traceur(zxtrfln(1,cas_evap_liq(i)), &
    539 &           'iso_revap_fisrtilp 804')
     553&           'iso_revap_fisrtilp 804b')
    540554     do ixt=1,ntraciso
    541555        call iso_verif_positif_choix(zxt(ixt,cas_evap_liq(i)), &
    … …  
    676690   if (iso_HDO.gt.0) then
    677691      do i=1,ncas_evap_glace
    678          if (zrfln_cas(i).gt.ridicule_rain) then
    679                 call iso_verif_aberrant( &
    680 &                  (zxtrfln_cas(iso_HDO,i) &
    681 &                 /zrfln_cas(i)), 'iso_revap_fisrtilp 4563')
    682          endif
     692                call iso_verif_aberrant_choix(zxtrfln_cas(iso_HDO,i), zrfln_cas(i), &
     693                  ridicule_rain,deltalim_snow, 'iso_revap_fisrtilp 4563')
    683694       enddo !do i=1,ncas_evap_glace
    684695   endif !if ((iso_HDO.gt.0).and.(ixt.eq.iso_HDO)) then
    … …  
    789800     call iso_verif_traceur(zxt(1,cas_evap_glace(i)), &
    790801&           'iso_revap_fisrtilp 1033')
     802     call iso_verif_traceur(zxtrfl(1,cas_evap_glace(i)), &
     803&           'iso_revap_fisrtilp 1035a')
    791804     call iso_verif_traceur(zxtrfln(1,cas_evap_glace(i)), &
    792 &           'iso_revap_fisrtilp 1035')
     805&           'iso_revap_fisrtilp 1035b')
    793806  enddo
    794807#endif         
    … …  
    830843    call iso_verif_traceur(zxt(1,i),'iso_revap_fisrtilp 532')
    831844    call iso_verif_traceur(zxtrfln(1,i), &
    832 &           'iso_revap_fisrtilp 533')
     845&           'iso_revap_fisrtilp 533a')
     846    call iso_verif_traceur(zxtrfl(1,i), &
     847&           'iso_revap_fisrtilp 533b')
    833848    do ixt=1,ntraciso
    834849        call iso_verif_positif_choix(zxt(ixt,i),0.0, &
    … …  
    836851    enddo
    837852enddo !do i=1,klon
    838         write(*,*) 'revap_ilp 814: sortie'
     853        !write(*,*) 'revap_ilp 814: sortie'
    839854#endif       
    840855#endif         
    … …  
    23842399    endif !if (iso_eau.gt.0) then
    23852400    if (abs(evap(i)).gt.ridicule_rain) then
    2386       if (iso_HDO.gt.0) then                 
    2387          if (iso_verif_aberrant_nostop(xtevap(iso_HDO,i)/evap(i), &
    2388 &           'iso_surf>iso_rosee_givre 3193').eq.1) then   
     2401      if (iso_HDO.gt.0) then                     
     2402         if (iso_verif_aberrant_choix_nostop(xtevap(iso_HDO,i),evap(i), &
     2403&           ridicule_rain,deltalim_snow,'iso_surf>iso_rosee_givre 3193').eq.1) then   
    23892404                write(*,*) 'zxtalphai(iso_HDO)=',zxtalphai(iso_HDO)
    23902405                write(*,*) 'deltaD1eff=',deltaD(xt1lay(iso_HDO,i)/q1lay(i))
    23912406                write(*,*) 'tsurf(i)=',tsurf(i)
    23922407                write(*,*) 'q1lay(i)=',q1lay(i)   
    2393                 stop 
     2408                !stop 
    23942409           endif !if (iso_verif_aberrant_nostop
    23952410      endif  !if (iso_HDO.gt.0) then 
    … …  
    25342549
    25352550! quelques verifs de bilan d'eau
    2536 #ifdef ISOVERIF      
     2551#ifdef ISOVERIF   
    25372552do il=1,ncas
    25382553  do ixt=1,niso
    … …  
    26472662&             'stewart_explicite_vectall 220')
    26482663    enddo
    2649 #endif
    2650 #ifdef ISOVERIF
    2651       if (iso_eau.gt.0) then
     2664    if (iso_eau.gt.0) then
    26522665        call iso_verif_egalite_choix( &
    26532666&              (Exi(iso_eau,il)*fac_ftmr(il)), &
    … …  
    26692682&                   *fac_ftmr(il))),'stewart_explicite 214')
    26702683       endif !if ((iso_HDO.gt.0).and.
    2671        if ((debug.eq.1).and.(il.eq.il_debug)) then
    2672            write(*,*) 'stewart_explicit 224: cas Pqisup<=0'
    2673            write(*,*) 'Eqi(il),deltaD=',Eqi(il), &
    2674 &                   deltaD((Exi(iso_HDO,il)/Eqi(il)))
    2675        endif
     2684
    26762685#endif             
    26772686else !if (Pqisup.eq.0) then
    … …  
    26792688h(il)=qeff(il)/qs(il)       
    26802689h(il)= MAX(MIN(h(il),1.0),0.0)
    2681 #ifdef ISOVERIF       
     2690#ifdef ISOVERIF
    26822691call iso_verif_positif(h(il)-thumxt1,'stewart_explicit 209')
    26832692#endif       
    … …  
    27392748      endif !if ((iso_HDO.gt.0).and.
    27402749    endif !if (iso_HDO.gt.0) then
    2741     if ((debug.eq.1).and.(il.eq.il_debug)) then
     2750    if ((debug.eq.1).and.(il.eq.il_debug).and.(Eqi(il).gt.0.)) then
    27422751           write(*,*) 'stewart_explicit 302: cas evap~0'
    27432752           write(*,*) 'deltaDv est inchangé:',deltaD( &
    … …  
    27792788f(il)=m(il)/m0(il)   
    27802789! verifs
    2781 #ifdef ISOVERIF        
     2790#ifdef ISOVERIF       
    27822791call iso_verif_positif((m(il)), &
    27832792&           'stewart_explicite 173')
    … …  
    28362845    ! rajout verif 4 sept 2009
    28372846    if (iso_HDO.gt.0) then
    2838      if (Pqisup(il).gt.ridicule) then
    2839         call iso_verif_aberrant((Rl0(iso_HDO,il)), &
    2840 &                   'stewart_explicite 368')
    2841      endif
    2842     endif
     2847        call iso_verif_aberrant_choix(Rl0(iso_HDO,il)*Pqisup(il),Pqisup(il), &
     2848&                   ridicule_rain,deltalim_snow,'stewart_explicite 368')
     2849    endif !if (iso_HDO.gt.0) then
    28432850  endif !(iso_eau.gt.0)
    28442851#endif
    … …  
    28682875&                   'stewart_explicite 271')
    28692876   enddo !do ixt=1,niso   
    2870 #endif
    2871 #ifdef ISOVERIF 
    28722877   if (iso_eau.gt.0) then
    28732878      call iso_verif_egalite_choix( &
    … …  
    29232928       endif !if (qp0(il)+Eqi(il)*fac_ftmr(il).gt.ridicule) then
    29242929     endif !if (iso_HDO.gt.0)
    2925      if ((debug.eq.1).and.(il.eq.il_debug)) then
     2930     if ((debug.eq.1).and.(il.eq.il_debug).and.(Eqi(il).gt.0.)) then
    29262931         write(*,*) 'stewart_explicit 442: tout se réévapore'
    29272932         write(*,*) 'Eqi(il),deltaD=',Eqi(il), &
    … …  
    29432948if ((h(il).gt.0.99).or. &
    29442949&           (h(il).gt.0.98).and.(f(il).lt.1e-3)) then
    2945 #ifdef ISOVERIF         
    2946 !            write(*,*) 'stewart_explicit 191: cas h=1: il=',il
    2947 #endif           
    29482950    do ixt=1,niso
    29492951    interm(ixt,il)=alphap(ixt,il)*(1.0-h(il)) &
    … …  
    29812983&           'stewart_explicite 261')
    29822984    enddo !do ixt=1,niso
    2983 #endif           
    2984 #ifdef ISOVERIF
    29852985    if (iso_eau.gt.0) then
    29862986        call iso_verif_egalite_choix( &
    … …  
    30073007     
    30083008     if (iso_HDO.gt.0) then
    3009        if (Pqiinf(il).gt.ridicule_rain) then
    3010            if (iso_verif_aberrant_nostop( &
    3011 &                  (Pxtiinf(iso_HDO,il)/Pqiinf(il)), &
    3012 &                   'stewart_explicie 248').eq.1) then
     3009       if (iso_verif_aberrant_choix_nostop(Pxtiinf(iso_HDO,il),Pqiinf(il), &
     3010&                ridicule_rain,deltalim_snow,'stewart_explicite 248').eq.1) then
    30133011             write(*,*) 'cas reeq totale, il=',il
    30143012             write(*,*) 'deltaDl0=',deltaD( &
    … …  
    30193017&                   (Rb(iso_hdo,il)))
    30203018             stop
    3021            endif
    3022        endif  !if ((iso_HDO.gt.0).and.(ixt.eq.iso_HDO).and.
     3019       endif !if (iso_verif_aberrant_choix_nostop
    30233020       if (qp0(il)+Eqi(il)*fac_ftmr(il).gt.ridicule) then
    30243021         call iso_verif_aberrant(( &
    … …  
    30273024       endif !if ((iso_HDO.gt.0).and.
    30283025     endif  !if (iso_HDO.gt.0) then
    3029     if ((debug.eq.1).and.(il.eq.il_debug)) then
     3026     if ((debug.eq.1).and.(il.eq.il_debug).and.(Eqi(il).gt.0.)) then
    30303027         write(*,*) 'stewart_explicit 526: cas h~1: rééq'
    30313028         write(*,*) 'Eqi(il),deltaD=',Eqi(il), &
    … …  
    30343031&                   (Rb0(iso_hdo,il))),deltaD( &
    30353032&                   (Rl0(iso_hdo,il)))
    3036     endif
     3033     endif !if ((debug.eq.1).and.(il.eq.il_debug)) then
    30373034#endif
    30383035    ! end verifs
    30393036
    30403037else if ((f(il).gt.0.998).and. &
    3041 &           (Eqi(il)*fac_ftmr(il).lt.1e-2*qp0(il))) then
     3038&           (Eqi(il)*fac_ftmr(il).lt.1e-2*qp0(il))) then ! if ((h(il).gt.0.99).or.
    30423039
    30433040!*** cas particulier pour éviter imprécisions numériques:
    … …  
    30833080&           'stewart_explicite 397')
    30843081   enddo !do ixt=1,niso
    3085 #endif
    3086 #ifdef ISOVERIF
    30873082   if (iso_eau.gt.0) then
    30883083        call iso_verif_egalite_choix( &
    … …  
    31673162           stop
    31683163         endif !if (iso_verif_aberrant_nostop((
    3169        endif !if ((iso_HDO.gt.0).and.
     3164       endif !if (qp0(il)+Eqi(il)*fac_ftmr(il).gt.ridicule) then
    31703165     endif !!if ((iso_HDO.gt.0)
    3171      if ((debug.eq.1).and.(il.eq.il_debug)) then
     3166     if ((debug.eq.1).and.(il.eq.il_debug).and.(Eqi(il).gt.0.)) then
    31723167         write(*,*) 'stewart_explicit 663: cas où réévap faible'
    31733168         write(*,*) 'ordre 1 pour la vapeur et le liquide'
    … …  
    32713266             endif !if (qp0(il)+Eqi(il)*fac_ftmr(il).gt.ridicule) then
    32723267          endif !if ((iso_HDO.gt.0)
    3273      if ((debug.eq.1).and.(il.eq.il_debug)) then
     3268     if ((debug.eq.1).and.(il.eq.il_debug).and.(Eqi(il).gt.0.)) then
    32743269         write(*,*) 'stewart_explicit 767: cas de réévap sèche'
    32753270         write(*,*) 'distill de Rayleigh'
    … …  
    33173312&                   'stewart_explicite 467')
    33183313    enddo !do ixt=1,niso   
    3319 #endif         
    3320 #ifdef ISOVERIF 
    33213314          if (iso_eau.gt.0) then
    33223315            call iso_verif_egalite_choix( &
    … …  
    33603353            endif !if (qp0(il)+Eqi(il)*fac_ftmr(il).gt.ridicule) then
    33613354          endif !if (iso_HDO.gt.0)
    3362       if ((debug.eq.1).and.(il.eq.il_debug)) then
     3355      if ((debug.eq.1).and.(il.eq.il_debug).and.(Eqi(il).gt.0.)) then
    33633356         write(*,*) 'stewart_explicit 831: flux de masse vap~0'
    33643357         write(*,*) 'Eqi(il),deltaD=',Eqi(il), &
    … …  
    34293422        enddo !do ixt=1,niso   
    34303423#endif
    3431 #ifdef ISOVERIF   
     3424#ifdef ISOVERIF
    34323425          if (iso_eau.gt.0) then
    34333426            call iso_verif_egalite_choix( &
    … …  
    35463539        enddo !do ixt=1,niso
    35473540#endif
    3548 #ifdef ISOVERIF      
     3541#ifdef ISOVERIF   
    35493542          if (iso_eau.gt.0) then
    35503543            call iso_verif_egalite_choix( &
    … …  
    37503743       ! pour meilleure convergence numérique:
    37513744       !xtnew=qp0+Eqi*fac_ftmr
    3752    endif   ! if (iso_eau.gt.0).and.(ixt.eq.iso_eau)   
    3753 
     3745   endif   ! if (iso_eau.gt.0).and.(ixt.eq.iso_eau) 
    37543746   if (iso_HDO.gt.0) then
    37553747     if (Pqiinf(il).gt.ridicule_rain) then
    … …  
    37783770&                   xtnew(iso_HDO,il)/(qp0(il)+Eqi(il) &
    37793771&                   *fac_ftmr(il))),'stewart_explicite 912')
    3780         endif !if (qp0(il)+Eqi(il)*fac_ftmr(il).gt.ridicule) then
    3781       endif !if (iso_HDO.gt.0)
    3782       if ((debug.eq.1).and.(il.eq.il_debug)) then
     3772       endif !if (qp0(il)+Eqi(il)*fac_ftmr(il).gt.ridicule) then
     3773       if ((debug.eq.1).and.(il.eq.il_debug).and.(Eqi(il).gt.0.)) then
    37833774         write(*,*) 'stewart_explicit 991: fcas général'
    37843775         write(*,*) 'mais avec formule simplifiée'
    3785          write(*,*) 'il,Eqi(il),deltaD=',il,Eqi(il), &
    3786 &                   deltaD((Exi(iso_HDO,il)/Eqi(il)))
    3787       endif
    3788 #endif     
     3776         write(*,*) 'il,Eqi(il)=',il,Eqi(il)
     3777         write(*,*) 'deltaD=',deltaD((Exi(iso_HDO,il)/Eqi(il)))
     3778       endif
     3779      endif !if (iso_HDO.gt.0) 
     3780#endif
    37893781  ! end verifs 
    37903782   
    … …  
    38383830! compression
    38393831if (ncas_Jsimple+ncas_rieman.gt.0) then
    3840 !#ifdef ISOVERIF          
     3832!#ifdef ISOVERIF         
    38413833!        write(*,*) 'stewart_explicite_vectall 873:compression_calculJ'
    38423834!#endif       
    … …  
    39273919    ! ******* traitement vectoriel du cas Rieman (=2533)
    39283920if (ncas_rieman.gt.0) then
    3929 !#ifdef ISOVERIF
    3930 !!          write(*,*) 'traitement vectoriel rieman: x',ncas_rieman
    3931 !          do icas_rieman=1+ncas_Jsimple,ncas_rieman+ncas_Jsimple
    3932 !!            write(*,*) 'ntot_cas(icas_rieman)=',ntot_cas(icas_rieman)
    3933 !          call iso_verif_positif(float(ntot_cas(icas_rieman))-1.0,
    3934 !     :           'stewart_expl 984: ntot faux')
    3935 !          enddo   !do icas_rieman=1,ncas_rieman       
    3936 !#endif           
     3921           
    39373922   icas_rieman=1+ncas_Jsimple
    3938 !           write(*,*) 'stewart_expl 988 tmp: icas_rieman=',icas_rieman
    3939 !           write(*,*) 'qp0_cas(1)=',qp0_cas(1)
    3940 !           write(*,*) 'A_cas(1)=',A_cas(1)
    3941 !           write(*,*) 'm0_cas(1)=',m0_cas(1)
    3942 !           write(*,*) 'm_cas(1)=',m_cas(1)
    3943 !           if (iso_eau.gt.0) then
    3944 !           write(*,*) 'beta_cas(iso_eau,1)=',beta_cas(iso_eau,1)
    3945 !           write(*,*) 'gama_cas(iso_eau,1)=',gama_cas(iso_eau,1)
    3946 !           endif
    3947 !           write(*,*) 'f_cas(1)=',f_cas(1)
    3948 !           write(*,*) 'g_cas(1)=',g_cas(1)
    3949 !#ifdef rieman           
    3950 !           call integrale_rieman_vectall
    3951 !     :       (ncas_rieman,m_cas(icas_rieman),
    3952 !     :       J(1,icas_rieman),e(1,icas_rieman),
    3953 !     :       qp0_cas(icas_rieman),A_cas(icas_rieman),
    3954 !     :       m0_cas(icas_rieman),beta_cas(1,icas_rieman),
    3955 !     :       gama_cas(1,icas_rieman),f_cas(icas_rieman),
    3956 !     :       g_cas(icas_rieman),ntot_cas(icas_rieman))
    3957 !#else
     3923
    39583924   call integrale_gauss_vectall &
    39593925&       (ncas_rieman,m_cas(icas_rieman), &
    … …  
    39643930!     :       g_cas(icas_rieman),ntot_cas(icas_rieman))
    39653931&       g_cas(icas_rieman))     
    3966 !#endif
    3967 !#ifdef ISOVERIF
    3968 !           do il=1+ncas_Jsimple,ncas_rieman+ncas_Jsimple
    3969 !             do ixt=1,niso
    3970 !               call integrale_Rieman_precision(m_cas(il),m0_cas(il),
    3971 !     :          Jtmp,etmp,ntot_cas(il)*1e2,
    3972 !     :           qp0_cas(il),A_cas(il),m0_cas(il),
    3973 !     :          beta_cas(ixt,il),gama_cas(ixt,il),f_cas(il),g_cas(il))
    3974 !               call iso_verif_egalite_choix((Jtmp),
    3975 !     :           (J(ixt,il)),
    3976 !     :           'stewart_exp 999: test rieman',errmax,errmaxrel)
    3977 !               if ((iso_eau.gt.0).and.(ixt.eq.iso_eau)) then
    3978 !                 write(*,*) 'stew exp tmp 1005: il,J(iso_eau,il),Jtmp=',
    3979 !     :           il,J(iso_eau,il),Jtmp
    3980 !               endif
    3981 !             enddo             
    3982 !           enddo
    3983 !#endif           
    39843932
    39853933  endif !if (ncas_rieman.gt.0) then
    … …  
    40083956&          -beta_cas(ixt,il)*gama_cas(ixt,il)*r_jqp0(ixt,il) &
    40093957&                   /g_cas(il)/g_cas(il))
    4010 !#ifdef ISOVERIF           
    4011 !            if ((iso_eau.gt.0).and.(iso_eau.eq.ixt)) then
    4012 !              write(*,*) 'stewart_explicite tmp 1071: il=',il
    4013 !              if (il.le.ncas_Jsimple) then
    4014 !                write(*,*) 'cas_Jsimple(il)=',cas_Jsimple(il)
    4015 !              else !if (il.le.ncas_Jsimple) then
    4016 !                write(*,*) 'cas_rieman(il)=',cas_rieman(il)
    4017 !              endif !if (il.le.ncas_Jsimple) then
    4018 !              write(*,*) 'f_cas(il),beta_cas(ixt,il),gama_cas(ixt,il)=',
    4019 !     :                   f_cas(il),beta_cas(ixt,il),gama_cas(ixt,il)
    4020 !              write(*,*) 'g_cas(il),r_jqp0(ixt,il),r_jl0(ixt,il)=',
    4021 !     ;          g_cas(il),r_jqp0(ixt,il),r_jl0(ixt,il)
    4022 !              write(*,*) 'A_cas(il)=',A_cas(il)
    4023 !              write(*,*) 'pond Rl0=',(f_cas(il)**beta_cas(ixt,il))
    4024 !     :          *(g_cas(il)**(-beta_cas(ixt,il)*gama_cas(ixt,il)))
    4025 !     :          +beta(ixt,il)*gama_cas(ixt,il)*r_jqp0(ixt,il)
    4026 !     :                   /f_cas(il)/g_cas(il)
    4027 !              write(*,*) 'pond Rb0=',gama_cas(ixt,il)*beta_cas(ixt,il)
    4028 !     :                   *r_jl0(ixt,il)/f_cas(il)/g_cas(il)
    4029 !              write(*,*) 'pondRl0=fac1*fac2+t3'
    4030 !              write(*,*) 'fac1=',
    4031 !     :           f_cas(il)**beta_cas(ixt,il)
    4032 !              write(*,*) 'fac2=',
    4033 !     :             g_cas(il)**(-beta_cas(ixt,il)*gama_cas(ixt,il))
    4034 !              write(*,*) 't3=',
    4035 !     :           beta_cas(ixt,il)*gama_cas(ixt,il)*r_jqp0(ixt,il)
    4036 !     :                   /f_cas(il)/g_cas(il)     
    4037 !            endif
    4038 !#endif
    40393958
    40403959    Pxtiinf_cas(ixt,il)=Pqiinf_cas(il)*Rl(ixt,il)
    … …  
    40974016    enddo
    40984017#endif
    4099 #ifdef ISOVERIF            
     4018#ifdef ISOVERIF         
    41004019    if (iso_eau.gt.0) then
    41014020      if (iso_verif_egalite_choix_nostop( &
    … …  
    44114330   if (iso_HDO.gt.0) then
    44124331    if (Pqiinf(il).gt.ridicule_rain) then
    4413       if (iso_verif_aberrant_nostop( &
    4414 &                  (Pxtiinf(iso_HDO,il)/Pqiinf(il)), &
    4415 &                   'stewart_explicie 871').eq.1)  then
     4332      if (iso_verif_aberrant_choix_nostop(Pxtiinf(iso_HDO,il),Pqiinf(il),ridicule_rain,deltalim_snow, &
     4333&                   'stewart_explicite 871').eq.1)  then
    44164334      write(*,*) 'deltaDl0=',deltaD( &
    44174335&           (Rl0(iso_HDO,il)))
    … …  
    44284346    endif !if (iso_HDO.gt.0)
    44294347
    4430     if ((debug.eq.1).and.(il.eq.il_debug)) then
     4348    if ((debug.eq.1).and.(il.eq.il_debug).and.(Eqi(il).gt.0.)) then
    44314349          write(*,*) 'stewart_explicit 1558: cas avec calcul J'
    44324350          write(*,*) 'Eqi(il),deltaD=',Eqi(il), &
    … …  
    44784396
    44794397#ifdef ISOVERIF
    4480 !          write(*,*) 'stewart_explicite vectall 1179: fin'
     4398          write(*,*) 'stewart_explicite vectall 1179: fin'
    44814399#endif         
    44824400
    … …  
    47024620&           'stewart_sublim_nofrac 39')             
    47034621endif !if ((iso_eau.gt.0).and.(ixt.eq.iso_eau)) then
    4704 if ((iso_HDO.gt.0).and. &
    4705 &           (Pqisup(il).gt.ridicule_rain)) then
    4706     call iso_verif_aberrant( &
    4707 &           (Pxtisup(iso_HDO,il)/Pqisup(il)), &
    4708 &           'stewart_sublim_nofrac 40')     
     4622if (iso_HDO.gt.0) then ! Camille 9 mars 2023: moins stricte pour condensat
     4623    call iso_verif_aberrant_choix(Pxtisup(iso_HDO,il),Pqisup(il), &
     4624&           ridicule_rain,deltalim_snow, 'stewart_sublim_nofrac 40')
    47094625endif !if ((iso_HDO.gt.0).and.(ixt.eq.iso_HDO).and.
    47104626enddo !do il=1,ncas
    … …  
    48344750        enddo ! do ixt=1,niso
    48354751        ! verif que deltaD(Pqiinf) raisonable
    4836         if ((iso_HDO.gt.0).and. &
    4837      &           (Pqiinf(il).gt.ridicule_rain)) then
    4838             call iso_verif_aberrant( &
    4839      &           (Pxtiinf(iso_HDO,il)/Pqiinf(il)), &
    4840      &           'stewart_sublim 175')
     4752        if (iso_HDO.gt.0) then
     4753            call iso_verif_aberrant_choix(Pxtiinf(iso_HDO,il),Pqiinf(il), &
     4754     &           ridicule_rain,deltalim_snow, 'stewart_sublim 175')
    48414755        endif !if ((iso_HDO.gt.0).and.(ixt.eq.iso_HDO).and.
    48424756        if (iso_eau.gt.0) then
    … …  
    49234837!         integer ntot_cas(ncas)
    49244838         integer il,ixt       
    4925          integer coeff_precision
    4926          parameter (coeff_precision=0.2)
    49274839
    49284840
    … …  
    1082710739     &          Exi_cas(1,1),Exi(1,1), &
    1082810740#endif       
    10829      &          xtp_avantevaptrac_cas,0,hdiag(1)) ! hdiag ne sera pas utilisé
     10741     &          xtp_avantevaptrac_cas,0,hdiag(1)) ! hdiag ne sera pas utilise
    1083010742
    1083110743       enddo ! do izone=1,ntraceurs_zone
    … …  
    1261512527                  if (iso_verif_aberrant_choix_nostop( &
    1261612528     &               zxtliq(iso_HDO,i),cond(i), &
    12617      &               ridicule,deltalim, &
     12529     &               ridicule,deltalim_snow, &
    1261812530     &               'condiso_liq_ice_vectall 32b').eq.1) then
    1261912531                    write(*,*) 'deltaDvap=',deltaD(xt(iso_hdo,i)/qt(i))
    … …  
    1263612548                        write(*,*) 'deltaO18(zxtliq/cond)=',deltaO(zxtliq(iso_O18,i)/cond(i))
    1263712549                        write(*,*) 'tcond(i)=',tcond(i)-t_coup,'°C'
    12638                         stop
     12550                        !stop ! Camille 9 mars 2023: trop strict
    1263912551                    endif  !if (iso_verif_O18_aberrant_nostop(     
    1264012552                    endif ! if (iso_O18.gt.0) then     
    … …  
    1271412626                        write(*,*) 'zxtalphai(iso_O18,i)=',zxtalphai(iso_O18,i)
    1271512627                        write(*,*) 'xt(1:niso,i)=',xt(1:niso,i)
    12716                         stop
     12628                        !stop ! Camille 9 mars 2023: trop strict
    1271712629                    endif  !if (iso_verif_O18_aberrant_nostop(     
    1271812630                    endif ! if (iso_O18.gt.0) then
    … …  
    1276512677          if (zfice(i).gt.0.9) then
    1276612678              if (iso_verif_aberrant_choix_nostop( &
    12767      &          zxtice(iso_HDO,i),cond(i),ridicule,deltalim_snow, &
     12679     &          zxtice(iso_HDO,i),cond(i),ridicule,deltalim_snow, &
     12680                ! Camille 9 mars 2023: pour le condensat, on laisse plus de
     12681                ! marge
    1276812682     &          'condiso_liq_ice_vect 412').eq.1) then
    1276912683                write(*,*) 'debug condiso_liq_ice_vect 449: i,zfice=', &
    … …  
    1347013384          endif !if (iso_eau.gt.0) then
    1347113385#ifdef ISOTRAC
    13472           call iso_verif_traceur(xtsnow_evap(1,i), &
     13386!          call iso_verif_traceur(xtsnow_evap(1,i), &
     13387!     &           'gestion neige 2146') ! attention car snow_evap parfois
     13388!     négatif -> il ne faut pas passer dans les verifs de positivité.
     13389          call iso_verif_traceur_justmass(xtsnow_evap(1,i), &
    1347313390     &           'gestion neige 2146')
    1347413391#endif         
    … …  
    1414114058        if (iso_HDO.gt.0) then
    1414214059            call iso_verif_aberrant_choix(-xtsol_evap(iso_HDO,i), &
    14143      &           sol_evap(i),ridicule_evap,deltalim, &
    14144      &           'calcul_iso_surf_sic 257')
     14060     &           sol_evap(i),ridicule_evap,deltalim_snow, &
     14061     &           'calcul_iso_surf_sic 257_sol_evap')
    1414514062        endif
    1414614063#endif           
    … …  
    1416214079              endif !if (iso_eau.gt.0) then
    1416314080              if (iso_HDO.gt.0) then
    14164                 if (evap(i).gt.ridicule_evap) then 
    14165                    call iso_verif_aberrant(xtevap(iso_HDO,i)/evap(i), &
    14166      &                  'calcul_iso_surf_sic 257')
    14167                endif !if (evap(i).gt.ridicule_evap) then 
     14081                   call iso_verif_aberrant_choix(xtevap(iso_HDO,i),evap(i), &
     14082     &                  ridicule_evap,deltalim_snow,'calcul_iso_surf_sic 257_evap')
    1416814083              endif !if (iso_eau.gt.0) then
    1416914084#ifdef ISOTRAC
    … …  
    1432314238        enddo !do i=1,knon
    1432414239        endif !if (iso_HDO.gt.0) then
    14325        
    14326         if (iso_eau.gt.0) then 
     14240           
    1432714241        do i=1,knon
     14242        if (iso_eau.gt.0) then
    1432814243          call iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i),  &
    14329      &         'calcul_iso_surf_lic_vectall 587',errmax,errmaxrel)
    14330         enddo
     14244     &         'calcul_iso_surf_lic_vectall 587a',errmax,errmaxrel)
    1433114245        endif
     14246        if (iso_HDO.gt.0) then
     14247          call iso_verif_aberrant_choix(xtsnow(iso_HDO,i), &
     14248     &               snow(i),ridicule_snow,deltalim_snow, &
     14249     &               'calcul_iso_surf_lic 587b')
     14250        endif
     14251        enddo !do i=1,knon
    1433214252#endif             
    1433314253       
    … …  
    1447614396     &                  'calcul_iso_surf_lic 363',errmax,errmaxrel) 
    1447714397              endif !if (iso_eau.gt.0) then
    14478               if (iso_HDO.gt.0) then 
    14479                 if (snow(i).gt.ridicule_snow) then 
    14480                  call iso_verif_aberrant(xtsnow(iso_HDO,i)/snow(i), &
    14481      &                  'calcul_iso_surf_lic 367')
    14482                 endif !if (snow(i).gt.ridicule_snow) then
     14398              if (iso_HDO.gt.0) then   
    1448314399                call iso_verif_aberrant_choix(xtsnow(iso_HDO,i), &
    14484      &               snow(i),ridicule,deltalim_snow, &
     14400     &               snow(i),ridicule_snow,deltalim_snow, &
    1448514401     &               'calcul_iso_surf_lic 797')
    14486                 if (abs(evap(i)).gt.ridicule_evap) then 
    14487                  call iso_verif_aberrant(xtevap(iso_HDO,i)/evap(i), &
    14488      &                  'calcul_iso_surf_lic 369')
    14489                 endif ! if (evap(i).gt.ridicule_evap) then
     14402                call iso_verif_aberrant_choix(xtevap(iso_HDO,i),evap(i), &
     14403     &                  ridicule_evap,deltalim_snow, 'calcul_iso_surf_lic 369')
    1449014404              endif !if (iso_eau.gt.0) then
    1449114405#ifdef ISOTRAC
    … …  
    1579515709              endif !if (iso_eau.gt.0) then
    1579615710              if (iso_HDO.gt.0) then 
    15797                  if (snow(i).gt.ridicule_snow) then
    15798                     call iso_verif_aberrant(xtsnow(iso_HDO,i)/snow(i), &
    15799      &                  'calcul_iso_surf_ter 749')
    15800                  endif !if (snow(i).gt.ridicule_snow) then
     15711                 call iso_verif_aberrant_choix(xtsnow(iso_HDO,i),snow(i), &
     15712     &                  ridicule_snow,deltalim_snow, 'calcul_iso_surf_ter 749')
    1580115713                 call iso_verif_aberrant_choix(xtsnow(iso_HDO,i), &
    1580215714     &               snow(i),ridicule,deltalim_snow, &
    … …  
    1580915721                     write(*,*) 'sol_evap,snow_evap=', &
    1581015722     &                   sol_evap(i),snow_evap(i)
    15811                      write(*,*) 'deltaDsol_evap=', &
     15723                     if (sol_evap(i).gt.ridicule_evap)write(*,*) 'deltaDsol_evap=', &
    1581215724     &                   deltaD(xtsol_evap(iso_hdo,i)/sol_evap(i))
    15813                      write(*,*) 'deltaDsnow_evap=', &
     15725                     if (snow_evap(i).gt.ridicule_evap)write(*,*) 'deltaDsnow_evap=', &
    1581415726     &                   deltaD(xtsnow_evap(iso_hdo,i)/snow_evap(i))
    1581515727                     write(*,*) 'deltaD1new=',deltaD( &
    … …  
    1609416006          endif !if (iso_eau.gt.0) then
    1609516007          if (iso_HDO.gt.0) then
    16096             if (snow(i,nsrf).gt.ridicule_snow) then
    16097             call iso_verif_aberrant(xtsnow(iso_hdo,i,nsrf)/snow(i,nsrf), &
    16098      &                'phyisoetat0 117')
    16099             endif
     16008            call iso_verif_aberrant_choix(xtsnow(iso_hdo,i,nsrf),snow(i,nsrf), &
     16009     &                ridicule_snow, deltalim_snow, 'phyisoetat0 117')
    1610016010          endif !if (iso_eau.gt.0) then
    1610116011        enddo !do nsrf=1,nbsrf
    … …  
    1618016090        real deltaD_snow_fall_O18,deltaD_rain_fall_O18
    1618116091        real alpha(niso),kcin(niso)
    16182         real q0,h0
    16183         parameter (q0=20e-3,h0=0.7)
    1618416092!      character*50 text
    1618516093       
    … …  
    1631416222        do k=1,klev
    1631516223          do ixt=1,niso
    16316            RMerlivat(ixt)=toce(ixt)/alpha(ixt) &
    16317      &                  *(1-kcin(ixt))/(1.0-kcin(ixt)*h0)
    16318            xt_ancien(ixt,i,k)=q_ancien(i,k)*RMerlivat(ixt) &
    16319      &                  *(min(q0,q_ancien(i,k))/q0)**(alpha(ixt)-1.0)
     16224           call iso_init_ideal(q_ancien(i,k),xt_ancien(ixt,i,k),ixt, &
     16225                alpha(ixt),kcin(ixt),toce(ixt))
     16226
    1632016227           if (q_ancien(i,k).gt.ridicule) then
    1632116228           xtl_ancien(ixt,i,k)=ql_ancien(i,k)*alpha(ixt) &
    … …  
    1633816245     &           'phyisoetat0 16067') 
    1633916246        enddo !do ixt=1,niso
    16340         if ((k.eq.1).and.(iso_HDO.gt.0).and.(iso_O18.gt.0) &
    16341                 .and.(abs(q_ancien(i,k)-q0).lt.1e-3)) then
    16342                 ! vérifier qu'on est proche de la fermeture de Merlivat
    16343             write(*,*) 'i,k=',i,k
    16344             write(*,*) 'q_ancien(i,k)=',q_ancien(i,k)
    16345             write(*,*) 'deltaD=',deltaD(xt_ancien(iso_HDO,i,k)/q_ancien(i,k))
    16346             write(*,*) 'deltaDM=',deltaD(RMerlivat(iso_HDO))
    16347             write(*,*) 'deltaO=',deltaO(xt_ancien(iso_O18,i,k)/q_ancien(i,k))
    16348             write(*,*) 'deltaOM=',deltaO(RMerlivat(iso_O18))
    16349             write(*,*) 'dexcess=',dexcess(xt_ancien(iso_HDO,i,k)/q_ancien(i,k), &
    16350                     xt_ancien(iso_O18,i,k)/q_ancien(i,k))
    16351             write(*,*) 'dexcessM=',dexcess(RMerlivat(iso_HDO),RMerlivat(iso_O18))
    16352             write(*,*) 'kcin=',kcin
    16353             write(*,*) 'toce=',toce
    16354             write(*,*) 'alpha=',alpha
    16355             call iso_verif_positif(20.0-abs(-80.0-deltaD(xt_ancien(iso_HDO,i,k)/q_ancien(i,k))), &
    16356                     'phyisoetat0 16398a')   
    16357             call iso_verif_positif(5.0-abs(10.0-dexcess(xt_ancien(iso_HDO,i,k)/q_ancien(i,k), &
    16358                     xt_ancien(iso_O18,i,k)/q_ancien(i,k))),'phyisoetat0 16398b')
    16359         endif
     16247
     16248        ! Camille 7 mars 2023: ajout d'un check
     16249        if ((i.eq.1).and.(k.eq.1).and.(iso_HDO.gt.0)) then
     16250        write(*,*) 'phyisoetat0 16362: q_ancien(1,1)=',q_ancien(1,1)
     16251        write(*,*) 'deltaD_ancien=',deltaD(xt_ancien(iso_HDO,i,k)/q_ancien(i,k))
     16252        write(*,*) 'xt_ancien(:,i,k)=',xt_ancien(:,i,k)
     16253        endif !if ((i.eq.1).and.(k.eq.1)) then
     16254
     16255        if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then
     16256            if (q_ancien(i,k).gt.ridicule) then 
     16257               if (iso_verif_o18_aberrant_nostop( &
     16258     &              xt_ancien(iso_HDO,i,k)/q_ancien(i,k), &
     16259     &              xt_ancien(iso_O18,i,k)/q_ancien(i,k), &
     16260     &              'phyisoetat0 16366 q_ancien').eq.1) then
     16261                  write(*,*) 'phyisoetat0 16367: i,k,q_ancien(i,k)=',i,k,q_ancien(i,k)
     16262                  write(*,*) 'xt_ancien(:,i,k)=',xt_ancien(:,i,k)
     16263                  stop
     16264              endif !  if (iso_verif_o18_aberrant_nostop
     16265            endif !if (q_seri(i,k).gt.errmax) then
     16266        endif !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then
    1636016267#endif
    1636116268
    … …  
    1654316450   CHARACTER(LEN=2) :: str2
    1654416451   CHARACTER(LEN=5) :: str5
    16545    CHARACTER(LEN=maxlen) :: outiso, oldIso, modname, nam(2)
     16452   CHARACTER(LEN=maxlen) :: outiso, oldIso, modname, nam(3), oldIso2
    1654616453   REAL :: xmin, xmax
    1654716454   LOGICAL :: found
    … …  
    1655916466      outiso = isoName(ixt)
    1656016467      oldIso = strTail(new2oldH2O(outiso), '_')            !--- Remove "H2O_" from "H2O_<iso>[_<tag>]"
     16468      i = INDEX(outiso, '_', .TRUE.)
     16469      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.
     16470!      write(*,*) 'tmp 16541:'
     16471!      write(*,*) 'outiso=',outiso
     16472!      write(*,*) 'oldIso=',oldIso
     16473!      write(*,*) 'oldIso2=',oldIso2
     16474
    1656116475      ! on lit seulement si ixt<=niso ou si on initialise les traceurs d'après fichier:
    1656216476#ifdef ISOTRAC
    1656316477      IF(ixt <= niso .OR. initialisation_isotrac == 0) THEN
    1656416478#endif
    16565       found = phyetat0iso_srf3(xtsnow,      "XTSNOW", "Surface snow", 0.)
    16566       if (.NOT.found) CALL abort_physic('isotopes_routines_mod', 'phyiso_etat0_fichier 16581: unfound isotopic variable',1)
    1656716479      found = phyetat0iso_srf3(fxtevap,     "XTEVAP", "evaporation",  0.)
     16480      if (.NOT.found) CALL abort_physic('isotopes_routines_mod', 'phyiso_etat0_fichier 16581a: unfound isotopic variable',1)
    1656816481      found = phyetat0iso_get2(xtrain_fall, "xtrain_f", "xrain fall", 0.)
    1656916482      found = phyetat0iso_get2(xtsnow_fall, "xtsnow_f", "xsnow fall", 0.)
    … …  
    1657116484      found = phyetat0iso_get3(xtl_ancien,  "XTLANCIEN", "QLANCIEN",  0.)
    1657216485      found = phyetat0iso_get3(xts_ancien,  "XTSANCIEN", "QSANCIEN",  0.)
    16573       found = phyetat0iso_get2(xtrun_off_lic_0, "XTRUNOFFLIC0", "RUNOFFLIC0", 0.)
    1657416486      found = phyetat0iso_get3(wake_deltaxt,  "WAKE_DELTAXT", "Delta hum. wake/env",  0.)
    1657516487#ifdef ISOVERIF
    … …  
    1658016492            DO nsrf = 1, nbsrf
    1658116493               CALL iso_verif_egalite(fxtevap(iso_eau,i,nsrf),fevap(i,nsrf),TRIM(modname)//' 231c')
    16582                CALL iso_verif_egalite( xtsnow(iso_eau,i,nsrf), snow(i,nsrf),TRIM(modname)//' 231d')
    1658316494            END DO
    1658416495         END DO
    … …  
    1659216503         END DO
    1659316504      END IF
    16594       IF(iso_eau > 0 .AND. ixt == iso_eau) THEN
    16595          DO i=1,klon
    16596             IF(iso_verif_egalite_nostop(run_off_lic_0(i),xtrun_off_lic_0(iso_eau,i),TRIM(modname)//' 326') == 1) THEN
    16597                WRITE(*,*) 'i=',i
    16598                STOP
    16599             END IF
    16600          END DO
    16601       END IF
    1660216505#endif
    1660316506      ! ces variables n'ont pas de traceurs:
    1660416507      IF(ixt <= niso) THEN
    1660516508         found = phyetat0iso_get2(xtsol, "XTSOL", "Surface humidity / bucket", 0.)
     16509         if (.NOT.found) CALL abort_physic('isotopes_routines_mod', 'phyiso_etat0_fichier 16581b: unfound isotopic variable',1)
    1660616510         found = phyetat0iso_get2(Rland_ice, "Rland_ice", "SR land ice", 0.)
    16607 #ifdef ISOVERIF
    16608 
     16511         found = phyetat0iso_srf3(xtsnow,      "XTSNOW", "Surface snow", 0.) ! CR avril 2023: deplacer ici
     16512         found = phyetat0iso_get2(xtrun_off_lic_0, "XTRUNOFFLIC0", "RUNOFFLIC0", 0.)
     16513#ifdef ISOVERIF
    1660916514         DO i=1,klon
    1661016515            IF(iso_verif_noNaN_nostop(xtsol(ixt,i),TRIM(modname)//' 95') == 1) THEN
    … …  
    1661216517               STOP
    1661316518            END IF
    16614          END DO
     16519            IF(ixt == iso_eau .AND. iso_eau > 0) THEN
     16520             DO nsrf = 1, nbsrf
     16521               CALL iso_verif_egalite(fxtevap(iso_eau,i,nsrf),fevap(i,nsrf),TRIM(modname)//' 231c')
     16522               CALL iso_verif_egalite( xtsnow(iso_eau,i,nsrf), snow(i,nsrf),TRIM(modname)//' 231d')
     16523             END DO
     16524             CALL iso_verif_egalite( xtrun_off_lic_0(iso_eau,i), run_off_lic_0(i),TRIM(modname)//' 231e')
     16525            ENDIF !IF(ixt == iso_eau .AND. iso_eau > 0) THEN
     16526         END DO !DO i=1,klon
    1661516527#endif
    1661616528      END IF
    … …  
    1668016592  nam(1) = TRIM(pref)//TRIM(outiso)
    1668116593  nam(2) = TRIM(pref)//TRIM(oldIso)
     16594  nam(3) = TRIM(pref)//TRIM(oldIso2)
    1668216595  lFound = phyetat0_get(iso_tmp, nam, descr, default)
    1668316596  field(ixt,:) = iso_tmp
    … …  
    1669216605  nam(1) = TRIM(pref)//TRIM(outiso)
    1669316606  nam(2) = TRIM(pref)//TRIM(oldIso)
     16607  nam(3) = TRIM(pref)//TRIM(oldIso2)
    1669416608  lFound = phyetat0_get(iso_tmp_lonlev, nam, descr, default)
    1669516609  field(ixt,:,:) = iso_tmp_lonlev(:,:)
    … …  
    1670316617  nam(1) = TRIM(pref)//TRIM(outiso)
    1670416618  nam(2) = TRIM(pref)//TRIM(oldIso)
     16619  nam(3) = TRIM(pref)//TRIM(oldIso2)
    1670516620  lFound = phyetat0_srf(iso_tmp_lonsrf, nam, descr, default)
    1670616621  field(ixt,:,:) = iso_tmp_lonsrf
    … …  
    1812918044
    1813018045      do iessai=1,nessai
    18131          day_nucl(iessai)   = 0.
    18132          month_nucl(iessai) = 0.
    18133          year_nucl(iessai)  = 0.
     18046         day_nucl(iessai)   = 0
     18047         month_nucl(iessai) = 0
     18048         year_nucl(iessai)  = 0
    1813418049         lat_nucl(iessai)   = 0.
    1813518050         lon_nucl(iessai)   = 0.
    … …  
    1881618731          do i=1,n
    1881718732            qttrac(i)=xt(ieau,i)
    18818 !            if (qt(i).gt.0.0) then ! modif C Risi juillet 2020
    18819             if ((qt(i).gt.0.0).and.(xt(ieau,i).gt.0.0)) then
    18820                zcondtrac(i)=(zcond(i)/qt(i))*xt(ieau,i)
     18733            if (qt(i).gt.0.0) then ! modif C Risi juillet 2020 ! remodif Camille 9 mars 2023
     18734!            if ((qt(i).gt.0.0).and.(xt(ieau,i).gt.0.0)) then
     18735               zcondtrac(i)=(zcond(i)/qt(i))*qttrac(i)
    1882118736            else !if (qt(i).eq.0) then
    1882218737#ifdef ISOVERIF             
    … …  
    1883618751            endif
    1883718752            if (iso_HDO.gt.0) then
    18838 !              if (qttrac(i).gt.ridicule_trac) then
    1883918753                call iso_verif_aberrant_choix(xttrac(iso_HDO,i), &
    1884018754     &           qttrac(i),ridicule_trac,deltalimtrac, &
    1884118755     &           'condisotrac 205')
    18842 !              endif
    1884318756            endif
    1884418757            call iso_verif_positif(qt(i)-cond(i), &
    … …  
    1885318766          call condiso_liq_ice_vectall(xttrac,qttrac,zcondtrac, &
    1885418767     &           tcond,zfice,zxticetrac,zxtliqtrac,n)
    18855 #ifdef ISOVERIF         
    18856           write(*,*) 'condisotrac 167: après condiso'
    18857 #endif           
     18768
    1885818769          do i=1,n
    1885918770          do iiso=1,niso         
    … …  
    1888118792#endif
    1888218793
     18794subroutine iso_init_ideal(q,xt,ixt,alpha,kcin,toce)
     18795
     18796        USE isotopes_mod, ONLY: iso_eau,iso_HDO,ridicule
     18797#ifdef ISOVERIF
     18798        USE isotopes_verif_mod
     18799#endif
     18800        implicit none
     18801
     18802        ! inputs
     18803        real q ! humidité spec
     18804        integer ixt ! indice isotopique
     18805        real alpha ! coef frac à l'eq
     18806        real kcin ! coef frac cinétique
     18807        real toce ! rapport iso ds ocean surface
     18808
     18809        ! outputs
     18810        real xt ! equivalent iso de l'humidité spec, même unité.
     18811
     18812        ! locals
     18813        real RMerlivat
     18814        real q0,h0 ! conditions initiales de la distill de Rayleigh
     18815        parameter (q0=20e-3,h0=0.7)
     18816
     18817        ! verifier que ixt est un isotope et pas un tagging
     18818        if (ixt.gt.niso) then
     18819          CALL abort_physic('isotopes_routines_mod', 'iso_init_ideal, ixt>niso', 1)
     18820        endif
     18821
     18822        ! R selon Merlivat:
     18823        RMerlivat=toce/alpha *(1.0-kcin)/(1.0-kcin*h0)
     18824
     18825        ! R d'après Rayleigh
     18826        xt=q*RMerlivat*(min(q0,q)/q0)**(alpha-1.0)
     18827
     18828#ifdef ISOVERIF
     18829        call iso_verif_noNaN(xt, 'isotopes_routines_mod 18930a: iso_init_ideal')
     18830        if ((iso_HDO.gt.0).and.(ixt.eq.iso_HDO)) then
     18831            if (q.gt.ridicule) then
     18832                call iso_verif_aberrant_encadre( xt/q, 'isotopes_routines_mod 18930b: iso_init_ideal')
     18833            endif
     18834        endif
     18835        if ((iso_eau.gt.0).and.(ixt.eq.iso_eau)) then
     18836             call iso_verif_egalite(xt,q, 'isotopes_routines_mod 18930c: iso_init_ideal')
     18837        endif
     18838#endif
     18839       
     18840
     18841end subroutine iso_init_ideal
     18842
     18843
     18844subroutine appel_stewart_debug(lwork,nloc,inb,na,i, &
     18845                evap,water,rpprec,rr,wdtrain, &
     18846                xtevap,xtwater,xtp,xt,xtwdtrain)
     18847USE isotopes_mod, ONLY: iso_eau, iso_HDO,thumxt1, &
     18848&       bidouille_anti_divergence,ridicule,Rdefault
     18849use infotrac_phy, ONLY: ntraciso=>ntiso, niso
     18850#ifdef ISOTRAC
     18851    use isotrac_mod, only: option_cond,izone_cond,index_iso,index_zone,izone_poubelle
     18852#endif
     18853#ifdef ISOVERIF
     18854USE isotopes_verif_mod
     18855#endif
     18856implicit none
     18857
     18858
     18859! inputs
     18860integer nloc,na,i ! dimension horiz effective
     18861logical lwork(nloc)
     18862real wdtrain(nloc),xtwdtrain(ntraciso,nloc)
     18863real xt(ntraciso,nloc,na)
     18864real evap(nloc,na),water(nloc,na),rpprec(nloc,na),rr(nloc,na)
     18865integer inb(nloc)
     18866
     18867! outputs
     18868real xtevap(ntraciso,nloc,na),xtwater(ntraciso,nloc,na),xtp(ntraciso,nloc,na)
     18869
     18870! locals
     18871integer il,ixt
     18872
     18873     do il=1,nloc
     18874       if (i.le.inb(il) .and. lwork(il)) then
     18875          if (wdtrain(il).gt.0.) then
     18876            do ixt=1,ntraciso
     18877             xtwater(ixt,il,i)= xtwdtrain(ixt,il)/wdtrain(il)*water(il,i)
     18878             xtevap(ixt,il,i)= xtwdtrain(ixt,il)/wdtrain(il)*evap(il,i)
     18879            enddo
     18880          else !if (wdtrain(il).gt.0.) then
     18881            do ixt=1,niso
     18882             xtwater(ixt,il,i)= Rdefault(ixt)*water(il,i)
     18883             xtevap(ixt,il,i)= Rdefault(ixt)*evap(il,i)
     18884            enddo
     18885#ifdef ISOTRAC
     18886            do ixt=1+niso,ntraciso
     18887             if (index_zone(ixt).eq.izone_poubelle) then
     18888               xtwater(ixt,il,i)= Rdefault(index_iso(ixt))*water(il,i)
     18889               xtevap(ixt,il,i)= Rdefault(index_iso(ixt))*evap(il,i)
     18890             else
     18891               xtwater(ixt,il,i)= 0.
     18892               xtevap(ixt,il,i)=0.
     18893             endif
     18894            enddo ! do ixt=1+niso,ntraciso
     18895#endif
     18896         endif !if (wdtrain(il).gt.0.) then
     18897         do ixt=1,ntraciso
     18898             xtp(ixt,il,i)= xt(ixt,il,i)/rr(il,i)*rpprec(il,i)
     18899         enddo !do ixt=1,ntraciso
     18900      endif
     18901    enddo ! do il=1,ncum
     18902end subroutine appel_stewart_debug
     18903
    1888318904END MODULE isotopes_routines_mod
    1888418905#endif

  • LMDZ6/trunk/libf/phylmdiso/isotopes_verif_mod.F90

    r4300 r4491  
    20282028
    20292029          if (iso_verif_egalite_choix_nostop(xtractot,x(iiso), &
    2030      &        err_msg//', verif trac egalite, iso '// &
     2030     &        err_msg//', verif trac egalite1, iso '// &
    20312031     &        TRIM(isoName(iiso)), &
    20322032     &        errmaxin,errmaxrelin).eq.1) then
    20332033            write(*,*) 'iso_verif_traceur 202: x=',x
    20342034!            write(*,*) 'xtractot=',xtractot
     2035            do izone=1,nzone 
     2036              ixt=itZonIso(izone,iiso)
     2037              write(*,*) 'izone,iiso,ixt=',izone,iiso,ixt
     2038            enddo
    20352039            iso_verif_tracm_choix_nostop=1
    20362040          endif
    … …  
    23982402        call iso_verif_egalite_std_vect( &
    23992403     &           xtractot,xiiso, &
    2400      &           err_msg//', verif trac egalite, iso ' &
     2404     &           err_msg//', verif trac egalite2, iso ' &
    24012405     &           //TRIM(isoName(iiso)), &
    24022406     &           n,m,errmax,errmaxrel)

  • LMDZ6/trunk/libf/phylmdiso/isotrac_routines_mod.F90

    r4325 r4491  
    508508!          write(*,*) 'compress 910: xtp_avantevap(iso_eau,cas(1))=',
    509509!     :           xtp_avantevap(iso_eau,cas(1))
     510!        write(*,*) 'compress_evap_liq_zone 510: ncas,ncum=',ncas,ncum
     511        ptrac(:)=0. ! CR 31 mars 2023: initialisation de ptrac
    510512
    511513          ieau=index_trac(izone,iso_eau)
    … …  
    519521            else
    520522#ifdef ISOVERIF
    521                 call iso_verif_egalite(( &
     523              call iso_verif_egalite(( &
    522524     &                   Eqi_prime(cas(il))),0.0, &
    523525     &                   'compress_stewart 979')
    … …  
    685687                call iso_verif_egalite(( &
    686688     &                   Eqi_prime(cas(il))),0.0, &
    687      &                   'compress_stewart 979')
     689     &                   'compress_stewart 979b')
    688690#endif               
    689691              Eqi_prime_cas(il)=0.0
    … …  
    870872            enddo ! do iiso=1,niso           
    871873          endif !if ((Eqi(il)*fac_ftmr(il).gt.evap_franche).and.
    872 #ifdef ISOVERIF
    873           if (il.eq.9) then           
    874             ixt=index_trac(2,2)
    875             write(*,*) 'il,ixt,xtrevap_tag(ixt,il)=', &
    876      &           il,xtrevap_tag(ixt,il)
    877           endif !if (il.eq.9) then
    878 #endif         
    879874         enddo !do il=1,ncas
    880875!         write(*,*) 'compress_stewart 1453 tmp: zxt=',
    … …  
    27912786       real xtractot
    27922787       
    2793         write(*,*) 'iso_verif_traceur tmp 822: xt(:,66,1)=',x(:,66,1)
    27942788        ! verif noNaN
    27952789        call iso_verif_traceur_noNaN_vect(x,n,m,err_msg)

  • LMDZ6/trunk/libf/phylmdiso/pbl_surface_mod.F90

    r4478 r4491  
    24072407            call iso_verif_noNaN(yxtevap(ixt,j), &
    24082408         &      'pbl_surface 1056a: apres surf_land')
     2409          enddo
     2410          do ixt=1,niso
    24092411            call iso_verif_noNaN(yxtsol(ixt,j), &
    24102412         &      'pbl_surface 1056b: apres surf_land')
    … …  
    24752477                   call iso_verif_noNaN(yxtevap(ixt,j), &
    24762478                        &      'pbl_surface 1095a: apres surf_landice')
     2479                enddo
     2480                do ixt=1,niso
    24772481                   call iso_verif_noNaN(yxtsol(ixt,j), &
    24782482                        &      'pbl_surface 1095b: apres surf_landice')
    … …  
    25702574            call iso_verif_noNaN(yxtevap(ixt,j), &
    25712575         &      'pbl_surface 1165a: apres surf_seaice')
     2576          enddo
     2577          do ixt=1,niso
    25722578            call iso_verif_noNaN(yxtsol(ixt,j), &
    25732579         &      'pbl_surface 1165b: apres surf_seaice')

  • LMDZ6/trunk/libf/phylmdiso/phyredem.F90

    r4389 r4491  
    581581      CALL put_field_srf1(pass, "XTEVAP"//TRIM(outiso), "Evaporation de surface",iso_tmp_lonsrf)
    582582
    583       iso_tmp_lonsrf(:,:)=xtsnow(ixt,:,:)
    584       CALL put_field_srf1(pass, "XTSNOW"//TRIM(outiso), "NEIGE",       iso_tmp_lonsrf)       
    585 
    586583      iso_tmp(:)=xtrain_fall(ixt,:)
    587584      CALL put_field(pass,    "xtrain_f"//TRIM(outiso), "precipitation liquide",iso_tmp)
    … …  
    602599      CALL put_field(pass,"WAKE_DELTAXT"//TRIM(outiso), "WAKE_DELTAQ", iso_tmp_lonlev)
    603600
     601      iso_tmp_lonlev(:,:)=wake_deltaxt(ixt,:,:)
     602      CALL put_field(pass,"WAKE_DELTAXT"//TRIM(outiso), "WAKE_DELTAXT",iso_tmp_lonlev)
     603
     604      ! variables seulement pour niso:
     605      if (ixt.le.niso) then
     606
     607      iso_tmp_lonsrf(:,:)=xtsnow(ixt,:,:)
     608      CALL put_field_srf1(pass, "XTSNOW"//TRIM(outiso), "NEIGE",       iso_tmp_lonsrf)
     609
     610      iso_tmp(:)=xtsol(ixt,:)
     611      CALL put_field(pass,      "XTSOL"//TRIM(outiso), "Eau dans le sol (mm)",iso_tmp)
     612
     613      iso_tmp(:)=Rland_ice(ixt,:)
     614      CALL put_field(pass,  "Rland_ice"//TRIM(outiso), "ratio land ice",      iso_tmp)
     615
    604616      iso_tmp(:)=xtrun_off_lic_0(ixt,:)
    605617      CALL put_field(pass,"XTRUNOFFLIC0"//TRIM(outiso), "Runofflic0",  iso_tmp)
    606618
    607       iso_tmp_lonlev(:,:)=wake_deltaxt(ixt,:,:)
    608       CALL put_field(pass,"WAKE_DELTAXT"//TRIM(outiso), "WAKE_DELTAXT",iso_tmp_lonlev)
    609 
    610       ! variables seulement pour niso:
    611       if (ixt.le.niso) then
    612 
    613       iso_tmp(:)=xtsol(ixt,:)
    614       CALL put_field(pass,      "XTSOL"//TRIM(outiso), "Eau dans le sol (mm)",iso_tmp)
    615 
    616       iso_tmp(:)=Rland_ice(ixt,:)
    617       CALL put_field(pass,  "Rland_ice"//TRIM(outiso), "ratio land ice",      iso_tmp)
    618 
    619619      endif ! if (ixt.le.niso) then
    620620

  • LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90

    r4479 r4491  
    157157        & iso_verif_aberrant_choix,iso_verif_positif, &
    158158        & iso_verif_positif_choix_vect,iso_verif_o18_aberrant_nostop, &
    159         & iso_verif_init,&
     159        & iso_verif_init,iso_verif_aberrant_enc_choix_nostop,&
    160160        & iso_verif_positif_strict_nostop,iso_verif_O18_aberrant_enc_vect2D
    161161#endif
    … …  
    24802480#ifdef ISO
    24812481#ifdef ISOVERIF
    2482 !    write(*,*) 'physiq 1847: qx(1,1,:)=',qx(1,1,:)
     2482    write(*,*) 'physiq 1847: qx(1,1,:)=',qx(1,1,:)
     2483    write(*,*) 'iqIsoPha(:,ivap)=',iqIsoPha(:,ivap)
    24832484    write(*,*) 'physiq 1846b: ok_isotopes,ntraciso,niso=',niso>0,ntraciso,niso
    24842485#endif
    … …  
    25962597     &              'physiq 2099 ql').eq.1) then
    25972598                  write(*,*) 'i,k,ql_seri(i,k)=',i,k,ql_seri(i,k)
    2598                   stop
     2599                  !stop
    25992600              endif !  if (iso_verif_o18_aberrant_nostop
    26002601            endif !if (q_seri(i,k).gt.errmax) then 
    … …  
    26052606     &              'physiq 2099 qs').eq.1) then
    26062607                  write(*,*) 'i,k,qs_seri(i,k)=',i,k,qs_seri(i,k)
    2607                   stop
     2608                  !stop
    26082609              endif !  if (iso_verif_o18_aberrant_nostop
    26092610            endif !if (q_seri(i,k).gt.errmax) then
    26102611          enddo !k=1,klev
    26112612         enddo  !i=1,klon
    2612         endif !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then
     2613        endif !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then         
     2614#ifdef ISOTRAC 
     2615         DO k = 1, klev
     2616           DO i = 1, klon   
     2617             call iso_verif_traceur(xt_seri(1,i,k),'physiq 2620a')
     2618             call iso_verif_traceur(xtl_seri(1,i,k),'physiq 2620b')
     2619             call iso_verif_traceur(xts_seri(1,i,k),'physiq 2620c')
     2620           enddo
     2621         enddo
     2622#endif
    26132623#endif
    26142624    !
    … …  
    28442854       ENDIF
    28452855    ENDIF
     2856
     2857#ifdef ISO
     2858#ifdef ISOVERIF
     2859#ifdef ISOTRAC
     2860         DO k = 1, klev
     2861           DO i = 1, klon         
     2862             call iso_verif_traceur(xt_seri(1,i,k), &
     2863     &           'physiq 2856: avant reevp')
     2864           enddo
     2865         enddo
     2866#endif
     2867#endif
     2868#endif
    28462869    !
    28472870    ! Re-evaporer l'eau liquide nuageuse
    … …  
    34583481                  call iso_verif_egalite(q_w(i,k),xt_w(iso_eau,i,k),'physiq 3338')
    34593482                endif
     3483                if (iso_HDO.gt.0) then
     3484                  if ((iso_verif_aberrant_enc_choix_nostop(xt_x(iso_hdo,i,k),q_x(i,k), &
     3485                        ridicule,deltalim,'physic 3462a xt_x').eq.1).or. &
     3486                        (iso_verif_aberrant_enc_choix_nostop(xt_w(iso_hdo,i,k),q_w(i,k), &
     3487                        ridicule,deltalim,'physic 3462b xt_w').eq.1)) then
     3488                     write(*,*) 'i,k=',i,k
     3489                     write(*,*) 'q_x(i,k),q_seri(i,k),wake_s(i),wake_deltaq(i,k)=', &
     3490                                q_x(i,k),q_seri(i,k),wake_s(i),wake_deltaq(i,k)
     3491                     write(*,*) 'xt_x(iso_hdo,i,k),xt_seri(iso_hdo,i,k),wake_s(i),wake_deltaxt(iso_hdo,i,k)=', &
     3492                                xt_x(iso_hdo,i,k),xt_seri(iso_hdo,i,k),wake_s(i),wake_deltaxt(iso_hdo,i,k)
     3493                     write(*,*) 'deltaD_seri,wake=',deltaD(xt_seri(iso_hdo,i,k)/q_seri(i,k)), &
     3494                                deltaD(wake_deltaxt(iso_hdo,i,k)/wake_deltaq(i,k))
     3495                     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))
     3496                     stop
     3497                  endif
     3498                endif
    34603499#endif
    34613500#endif
    … …  
    35903629              call iso_verif_aberrant_encadre( &
    35913630     &           xt_w(iso_hdo,i,k)/q_w(i,k), &
    3592      &          'physic 2657b')
     3631     &          'physic 2657c')
    35933632             endif !if (q_x(i,k).gt.ridicule) then
    35943633           endif !if (iso_HDO.gt.0) then
    … …  
    42554294             endif !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then   
    42564295             if ((iso_O17.gt.0).and.(iso_O18.gt.0)) then
    4257                if ((q_seri(i,k).gt.ridicule).and.(l.lt.nlevmaxO17)) then
     4296               if ((q_seri(i,k).gt.ridicule).and.(k.lt.nlevmaxO17)) then
    42584297                 call iso_verif_aberrant_o17(xt_seri(iso_o17,i,k) &
    42594298     &           /q_seri(i,k),xt_seri(iso_o18,i,k) &
    … …  
    46194658        do i=1,klon
    46204659           if ((iso_O17.gt.0).and.(iso_O18.gt.0)) then
    4621             if ((q_seri(i,k).gt.ridicule).and.(l.lt.nlevmaxO17)) then
     4660            if ((q_seri(i,k).gt.ridicule).and.(k.lt.nlevmaxO17)) then
    46224661             call iso_verif_aberrant_o17(xt_seri(iso_o17,i,k) &
    46234662     &           /q_seri(i,k),xt_seri(iso_o18,i,k) &
    … …  
    48414880            endif !if (q_seri(i,k).gt.errmax) then
    48424881            if (ql_seri(i,k).gt.ridicule) then
    4843                call iso_verif_aberrant( &
    4844      &           xtl_seri(iso_HDO,i,k)/ql_seri(i,k),'physiq 2871')
     4882               call iso_verif_aberrant_choix(xtl_seri(iso_HDO,i,k),ql_seri(i,k), &
     4883                        ridicule,deltalim_snow,'physiq 2871')
    48454884               if (iso_O18.gt.0) then 
    48464885                 if (iso_verif_o18_aberrant_nostop( &
    … …  
    48544893                        write(*,*) 'deltaO(d_ql_lsc(i,k))=',deltaO( &
    48554894     &                           d_xtl_lsc(iso_O18,i,k)/d_ql_lsc(i,k))
    4856                         stop
     4895                        !stop
    48574896                 endif
    48584897               endif ! if (iso_O18.gt.0) then 
    … …  
    48644903        do i=1,klon
    48654904          do k=1,nlev
    4866            if ((q_seri(i,k).gt.ridicule).and.(l.lt.nlevmaxO17)) then
     4905           if ((q_seri(i,k).gt.ridicule).and.(k.lt.nlevmaxO17)) then
    48674906            call iso_verif_aberrant_o17(xt_seri(iso_o17,i,k) &
    48684907     &           /q_seri(i,k),xt_seri(iso_o18,i,k) &
    … …  
    61496188     &              xt_seri(iso_HDO,i,k)/q_seri(i,k), &
    61506189     &              xt_seri(iso_O18,i,k)/q_seri(i,k), &
    6151      &              'physiq 5937, juste apres methox, qv').eq.1) then
     6190     &              'physiq 5937a, juste apres methox, qv').eq.1) then
    61526191                  write(*,*) 'physic 2444: i,k,q_seri(i,k)=',i,k,q_seri(i,k)
    61536192                  write(*,*) 'xt_seri(:,i,k)=',xt_seri(:,i,k)
    … …  
    61596198     &              xtl_seri(iso_HDO,i,k)/ql_seri(i,k), &
    61606199     &              xtl_seri(iso_O18,i,k)/ql_seri(i,k), &
    6161      &              'physiq 5937, juste apres methox, ql').eq.1) then
     6200     &              'physiq 5937b, juste apres methox, ql').eq.1) then
    61626201                  write(*,*) 'i,k,ql_seri(i,k)=',i,k,ql_seri(i,k)
    6163                   stop
     6202                  !stop
    61646203              endif !  if (iso_verif_o18_aberrant_nostop
    61656204            endif !if (q_seri(i,k).gt.errmax) then 
    … …  
    61686207     &              xts_seri(iso_HDO,i,k)/qs_seri(i,k), &
    61696208     &              xts_seri(iso_O18,i,k)/qs_seri(i,k), &
    6170      &              'physiq 5937, juste apres methox, qs').eq.1) then
     6209     &              'physiq 5937c, juste apres methox, qs').eq.1) then
    61716210                  write(*,*) 'i,k,qs_seri(i,k)=',i,k,qs_seri(i,k)
    6172                   stop
     6211                  !stop
    61736212              endif !  if (iso_verif_o18_aberrant_nostop
    61746213            endif !if (q_seri(i,k).gt.errmax) then

  • LMDZ6/trunk/libf/phylmdiso/reevap.F90

    r4143 r4491  
    8585      do ixt=1,ntiso
    8686        call iso_verif_noNaN(xt_seri(ixt,i,k), &
    87      &     'physiq 2417: apres evap tot')
     87     &     'reevap 2417: apres evap tot')
    8888      enddo
    8989      if (iso_eau.gt.0) then
    9090              call iso_verif_egalite_choix( &
    9191     &           xt_seri(iso_eau,i,k),q_seri(i,k), &
    92      &          'physiq 1891+, après reevap totale',errmax,errmaxrel)
     92     &          'reevap 1891+, après reevap totale',errmax,errmaxrel)
    9393              call iso_verif_egalite_choix( &
    9494     &           xtl_seri(iso_eau,i,k),ql_seri(i,k), &
    95      &          'physiq 2209+, après reevap totale',errmax,errmaxrel)
     95     &          'reevap 2209+, après reevap totale',errmax,errmaxrel)
    9696       endif !if (iso_eau.gt.0) then       
    9797      if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then 
    … …  
    100100     &              xt_seri(iso_HDO,i,k)/q_seri(i,k), &
    101101     &              xt_seri(iso_O18,i,k)/q_seri(i,k), &
    102      &              'physiq 2315: apres reevap totale').eq.1) then
     102     &              'reevap 2315: apres reevap totale').eq.1) then
    103103                  write(*,*) 'i,k,q_seri(i,k)=',i,k,q_seri(i,k)
    104104                  write(*,*) 'd_q_eva(i,k)=',d_q_eva(i,k)
    … …  
    111111#ifdef ISOTRAC     
    112112             call iso_verif_traceur(xt_seri(1,i,k), &
    113      &           'physiq 2165')
     113     &           'reevap 2165a')
    114114             call iso_verif_traceur_pbidouille(xt_seri(1,i,k), &
    115      &           'physiq 2165b')
     115     &           'reevap 2165b')
    116116#endif               
    117117         
    … …  
    146146      do ixt=1,ntiso
    147147      call iso_verif_noNaN(xt_seri(ixt,i,k), &
    148      &     'physiq 2417: apres evap tot')
     148     &     'reevap 2417: apres evap tot')
    149149      enddo
    150150      if (iso_eau.gt.0) then
    151151              call iso_verif_egalite_choix( &
    152152     &           xt_seri(iso_eau,i,k),q_seri(i,k), &
    153      &          'physiq 1891, après réévap totale',errmax,errmaxrel)
     153     &          'reevap 1891, après réévap totale',errmax,errmaxrel)
    154154              call iso_verif_egalite_choix( &
    155155     &           xtl_seri(iso_eau,i,k),ql_seri(i,k), &
    156      &          'physiq 2209, après réévap totale',errmax,errmaxrel)
     156     &          'reevap 2209, après réévap totale',errmax,errmaxrel)
    157157              call iso_verif_egalite_choix( &
    158158     &           xts_seri(iso_eau,i,k),qs_seri(i,k), &
    159      &          'physiq 2209b, après réévap totale',errmax,errmaxrel)
     159     &          'reevap 2209b, après réévap totale',errmax,errmaxrel)
    160160       endif !if (iso_eau.gt.0) then
    161161     
    … …  
    165165     &              xt_seri(iso_HDO,i,k)/q_seri(i,k), &
    166166     &              xt_seri(iso_O18,i,k)/q_seri(i,k), &
    167      &              'physiq 2408: apres reevap totale').eq.1) then
     167     &              'reevap 2408: apres reevap totale').eq.1) then
    168168                  write(*,*) 'i,k,q_seri(i,k)=',i,k,q_seri(i,k)                       
    169169                  stop
    … …  
    173173#ifdef ISOTRAC     
    174174             call iso_verif_traceur(xt_seri(1,i,k), &
    175      &           'physiq 2165')
     175     &           'reevap 2165c')
    176176             call iso_verif_traceur_pbidouille(xt_seri(1,i,k), &
    177      &           'physiq 2165b')
     177     &           'reevap 2165d')
    178178#endif                 
    179179#endif                 

  • LMDZ6/trunk/libf/phylmdiso/wake.F90

    r4374 r4491  
    3838  USE isotopes_verif_mod
    3939!, ONLY: errmax,errmaxrel
    40   USE isotopes_mod, ONLY: iso_eau,iso_hdo
     40  USE isotopes_mod, ONLY: iso_eau,iso_hdo,ridicule
    4141#endif
    4242#endif
    … …  
    714714          xtu(ixt,i,k)  =  xte(ixt,i,k) - deltaxtw(ixt,i,k)*sigmaw(i)
    715715        enddo !do ixt=1,ntraciso
    716 #endif
    717     END DO
    718   END DO
    719 
    720 
    721 #ifdef ISO
    722716#ifdef ISOVERIF
    723         ! TODO   
    724 #endif
    725 #endif
     717        if (iso_eau.gt.0) then
     718              call iso_verif_egalite(deltaqw(i,k),deltaxtw(iso_eau,i,k),'wake 723a')
     719              call iso_verif_egalite(qu(i,k),xtu(iso_eau,i,k),'wake 723b')
     720        endif
     721        if (iso_HDO.gt.0) then
     722               if (iso_verif_aberrant_enc_choix_nostop(xtu(iso_hdo,i,k),qu(i,k),ridicule,deltalim, &
     723                        'wake 723c xtu').eq.1) then
     724                     stop
     725               endif
     726         endif 
     727#endif
     728#endif
     729    END DO
     730  END DO
     731
     732
    726733
    727734  DO k = 1, klev - 1
    … …  
    18371844#ifdef ISO
    18381845#ifdef ISOVERIF
     1846      write(*,*) 'wake 1859'
    18391847      if (iso_eau.gt.0) then
    18401848        DO k= 1,klev
Note: See TracChangeset for help on using the changeset viewer.