Index: LMDZ6/branches/DYNAMICO-conv/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/branches/DYNAMICO-conv/libf/phylmd/physiq_mod.F90	(revision 3355)
+++ LMDZ6/branches/DYNAMICO-conv/libf/phylmd/physiq_mod.F90	(revision 3356)
@@ -63,6 +63,4 @@
        !
        d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_t_diss, &
-       d_t_vdf_w,d_q_vdf_w, &
-       d_t_vdf_x,d_q_vdf_x, &
        d_ts, &
        !
@@ -127,4 +125,5 @@
        slab_wfbils, tpot, tpote,               &
        ue, uq, ve, vq, zxffonte,               &
+       uwat, vwat,                             &
        zxfqcalving, zxfluxlat,                 &
        zxrunofflic,                            &
@@ -135,6 +134,6 @@
        zxfluxlat_x, zxfluxlat_w, &
        !
-       dtvdf_x, dtvdf_w, &
-       dqvdf_x, dqvdf_w, &
+       d_t_vdf_x, d_t_vdf_w, &
+       d_q_vdf_x, d_q_vdf_w, &
        pbl_tke_input, &
        t_therm, q_therm, u_therm, v_therm, &
@@ -144,5 +143,5 @@
        !
        wake_k, &
-       ale_wake, alp_wake, & 
+       alp_wake, & 
        wake_h, wake_omg, &
                        ! tendencies of delta T and delta q:
@@ -153,20 +152,37 @@
        d_deltat_ajs_cv, d_deltaq_ajs_cv, & ! due to dry adjustment of (w) before convection
                        ! tendencies of wake fractional area and wake number per unit area:
-       d_s_wk,  d_dens_wk, &             ! due to wakes
-!!!       d_s_vdf, d_dens_vdf, &            ! due to vertical diffusion
-!!!       d_s_the, d_dens_the, &            ! due to thermals
+       d_s_wk,  d_dens_a_wk,  d_dens_wk, &  ! due to wakes
+!!!       d_s_vdf, d_dens_a_vdf, d_dens_vdf, & ! due to vertical diffusion
+!!!       d_s_the, d_dens_a_the, d_dens_the, & ! due to thermals
        !                                  
-       ptconv, &
+       ptconv, ratqsc, &
        wbeff, convoccur, zmax_th, &
        sens, flwp, fiwp,  &
-       ale_bl_stat,alp_bl_conv,alp_bl_det,  &
+       alp_bl_conv,alp_bl_det,  &
        alp_bl_fluct_m,alp_bl_fluct_tke,  &
        alp_bl_stat, n2, s2,  &
        proba_notrig, random_notrig,  &
-       !
-       dnwd, dnwd0,  &
-       upwd, omega,  &
+       cv_gen,  &
+       !
+       dnwd0,  &
+       omega,  &
        epmax_diag,  &
+
+       !    Deep convective variables used in phytrac
+       pmflxr, pmflxs,  &
+       wdtrainA, wdtrainM,  &
+       upwd, dnwd, &
        ep,  &
+       da, mp, &
+       phi, &
+       wght_cvfd, &
+       phi2, &
+       d1a, dam, &
+       ev, &
+       elij, &
+       clw, &
+       epmlmMm, eplaMm, &
+       sij, &
+
        cldemi,  &
        cldfra, cldtau, fiwc,  &
@@ -181,6 +197,5 @@
        fsolsw, wfbils, wfbilo,  &
        wfevap, wfrain, wfsnow,  &  
-       pmflxr, pmflxs, prfl,  &
-       psfl, fraca, Vprecip,  &
+       prfl, psfl, fraca, Vprecip,  &
        zw2,  &
        
@@ -192,5 +207,4 @@
        qwriteSTD, twriteSTD, rhwriteSTD, &    !pour calcul_STDlev.h
        
-       wdtrainA, wdtrainM,  &
        beta_prec,  &
        rneb,  &
@@ -333,8 +347,8 @@
     INTEGER igout
     !======================================================================
-    ! Clef controlant l'activation du cycle diurne:
+    ! Clef iflag_cycle_diurne controlant l'activation du cycle diurne:
     ! en attente du codage des cles par Fred
-    INTEGER iflag_cycle_diurne
-    PARAMETER (iflag_cycle_diurne=1)
+    ! iflag_cycle_diurne est initialise par conf_phys et se trouve 
+    ! dans clesphys.h (IM)
     !======================================================================
     ! Modele thermique du sol, a activer pour le cycle diurne:
@@ -446,18 +460,19 @@
     REAL dtadd(klon,klev)
 
-    ! Variables pour le transport convectif
-    real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
-    real wght_cvfd(klon,klev)
 #ifndef CPP_XIOS
     REAL, SAVE :: missing_val=nf90_fill_real
 #endif
-    ! Variables pour le lessivage convectif
-    ! RomP >>> 
-    real phi2(klon,klev,klev)
-    real d1a(klon,klev),dam(klon,klev)
-    real ev(klon,klev)
-    real clw(klon,klev),elij(klon,klev,klev)
-    real epmlmMm(klon,klev,klev),eplaMm(klon,klev)
-    ! RomP <<<
+!!   Variables moved to phys_local_var_mod
+!!    ! Variables pour le transport convectif
+!!    real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
+!!    real wght_cvfd(klon,klev)
+!!    ! Variables pour le lessivage convectif
+!!    ! RomP >>> 
+!!    real phi2(klon,klev,klev)
+!!    real d1a(klon,klev),dam(klon,klev)
+!!    real ev(klon,klev)
+!!    real clw(klon,klev),elij(klon,klev,klev)
+!!    real epmlmMm(klon,klev,klev),eplaMm(klon,klev)
+!!    ! RomP <<<
     !IM definition dynamique o_trac dans phys_output_open
     !      type(ctrl_out) :: o_trac(nqtot)
@@ -534,4 +549,8 @@
     INTEGER k_upper_cv
     !------------------------------------------------------------------
+    ! Compteur de l'occurence de cvpas=1
+    INTEGER Ncvpaseq1
+    SAVE Ncvpaseq1
+    !$OMP THREADPRIVATE(Ncvpaseq1)
     !
     !==========================================================================
@@ -549,5 +568,6 @@
     ! variables supplementaires de concvl
     REAL Tconv(klon,klev)
-    REAL sij(klon,klev,klev)
+!!    variable moved to phys_local_var_mod
+!!    REAL sij(klon,klev,klev)
 !!    !
 !!    ! variables pour tester la conservation de l'energie dans concvl
@@ -624,4 +644,9 @@
     REAL, SAVE :: alp_offset
     !$OMP THREADPRIVATE(alp_offset)
+    REAL, SAVE :: dtcon_multistep_max=1.e6
+    !$OMP THREADPRIVATE(dtcon_multistep_max)
+    REAL, SAVE :: dqcon_multistep_max=1.e6
+    !$OMP THREADPRIVATE(dqcon_multistep_max)
+
   
     !
@@ -874,4 +899,14 @@
     INTEGER :: flag_inhib_tend = 0 !  0 is the default value
 !!    INTEGER :: flag_inhib_tend = 2
+    !
+    ! Logical switch to a bug : reseting to 0 convective variables at the 
+    ! begining of physiq.
+    LOGICAL, SAVE :: ok_bug_cv_trac = .TRUE.
+    !$OMP THREADPRIVATE(ok_bug_cv_trac)
+    !
+    ! Logical switch to a bug : changing wake_deltat when thermals are active
+    ! even when there are no wakes.
+    LOGICAL, SAVE :: ok_bug_split_th = .TRUE.
+    !$OMP THREADPRIVATE(ok_bug_split_th)
 
     !
@@ -886,5 +921,4 @@
     INTEGER kcbot(klon), kctop(klon), kdtop(klon)
     !
-    REAL ratqsc(klon,klev)
     real ratqsbas,ratqshaut,tau_ratqs
     save ratqsbas,ratqshaut,tau_ratqs
@@ -1111,4 +1145,8 @@
     REAL, dimension(klon, klev) :: cldfrarad   ! fraction nuageuse
 
+    !lwoff=y : offset LW CRE for radiation code and other schemes
+    REAL, SAVE :: betalwoff 
+    !OMP THREADPRIVATE(betalwoff)
+!
     INTEGER :: nbtr_tmp ! Number of tracer inside concvl
     REAL, dimension(klon,klev) :: sh_in ! Specific humidity entering in phytrac
@@ -1125,4 +1163,8 @@
     logical, parameter :: mass_fixer=.false.
     real qql1(klon),qql2(klon),corrqql
+
+    REAL pi
+
+    pi = 4. * ATAN(1.)
 
     ! Ehouarn: set value of jjmp1 since it is no longer a "fixed parameter"
@@ -1204,4 +1246,5 @@
        ENDIF
 
+       Ncvpaseq1 = 0
        dnwd0=0.0
        ftd=0.0
@@ -1211,5 +1254,4 @@
        pbase=0
        !IM 180608
-
 
        itau_con=0
@@ -1237,4 +1279,11 @@
     IF (debut) THEN
        CALL suphel ! initialiser constantes et parametres phys.
+! tau_gl : constante de rappel de la temperature a la surface de la glace - en
+       tau_gl=5.
+       CALL getin_p('tau_gl', tau_gl)
+! tau_gl : constante de rappel de la temperature a la surface de la glace - en
+! secondes
+       tau_gl=86400.*tau_gl
+       print*,'debut physiq_mod tau_gl=',tau_gl
        CALL getin_p('iflag_alp_wk_cond', iflag_alp_wk_cond) 
        CALL getin_p('random_notrig_max',random_notrig_max)
@@ -1245,4 +1294,6 @@
                       ! 2 => convective adjustment and state variables are changed
        CALL getin_p('iflag_adjwk',iflag_adjwk)
+       CALL getin_p('dtcon_multistep_max',dtcon_multistep_max)
+       CALL getin_p('dqcon_multistep_max',dqcon_multistep_max)
        CALL getin_p('oliqmax',oliqmax)
        CALL getin_p('oicemax',oicemax)
@@ -1254,4 +1305,6 @@
                                   ! in rrtm/suphec.F90 (and rvtmp2 is set to 0).
        CALL getin_p('ok_bad_ecmwf_thermo',ok_bad_ecmwf_thermo)
+       CALL getin_p('ok_bug_cv_trac',ok_bug_cv_trac)
+       CALL getin_p('ok_bug_split_th',ok_bug_split_th)
        fl_ebil = 0 ! by default, conservation diagnostics are desactivated
        CALL getin_p('fl_ebil',fl_ebil)
@@ -1328,5 +1381,5 @@
        print*,'iflag_coupl,iflag_clos,iflag_wake', &
             iflag_coupl,iflag_clos,iflag_wake
-       print*,'iflag_CYCLE_DIURNE', iflag_cycle_diurne
+       print*,'iflag_cycle_diurne', iflag_cycle_diurne
        !
        IF (iflag_con.EQ.2.AND.iflag_cld_th.GT.-1) THEN
@@ -1379,5 +1432,6 @@
        print *,'physiq, nbapp_cv, nbapp_wk ',nbapp_cv,nbapp_wk
        IF (MOD(NINT(86400./phys_tstep),nbapp_cv).EQ.0) THEN
-          cvpas = NINT( 86400./phys_tstep)/nbapp_cv
+          cvpas_0 = NINT( 86400./phys_tstep)/nbapp_cv
+          cvpas = cvpas_0
        print *,'physiq, cvpas ',cvpas
        ELSE 
@@ -1510,11 +1564,15 @@
 !jyg<
        IF (klon_glo==1) THEN
-          pbl_tke(:,:,is_ave) = 0.
-          DO nsrf=1,nbsrf
-            DO k = 1,klev+1
-                 pbl_tke(:,k,is_ave) = pbl_tke(:,k,is_ave) &
-                     +pctsrf(:,nsrf)*pbl_tke(:,k,nsrf)
-            ENDDO
-          ENDDO
+          IF (iflag_pbl > 1) THEN          
+              pbl_tke(:,:,is_ave) = 0.
+              DO nsrf=1,nbsrf
+                DO k = 1,klev+1
+                     pbl_tke(:,k,is_ave) = pbl_tke(:,k,is_ave) &
+                         +pctsrf(:,nsrf)*pbl_tke(:,k,nsrf)
+                ENDDO
+              ENDDO
+          ELSE   ! (iflag_pbl > 1)
+              pbl_tke(:,:,:) = 0.
+          ENDIF  ! (iflag_pbl > 1)
         ELSE
           pbl_tke(:,:,is_ave) = 0. !ym missing init : maybe must be initialized in the same way that for klon_glo==1 ??
@@ -1686,27 +1744,12 @@
 
 #ifdef CPP_XIOS 
-       !--setting up swaero_diag to TRUE in XIOS case 
-       IF (xios_field_is_active("topswad").OR.xios_field_is_active("topswad0").OR. & 
-           xios_field_is_active("solswad").OR.xios_field_is_active("solswad0").OR. & 
-           xios_field_is_active("topswai").OR.xios_field_is_active("solswai").OR.  & 
-             (iflag_rrtm==1.AND.(xios_field_is_active("toplwad").OR.xios_field_is_active("toplwad0").OR. & 
-                                 xios_field_is_active("sollwad").OR.xios_field_is_active("sollwad0"))))  & 
-           !!!--for now these fields are not in the XML files so they are omitted 
-           !!!  xios_field_is_active("toplwai").OR.xios_field_is_active("sollwai") !))) & 
-           swaero_diag=.TRUE. 
-
-       !--setting up dryaod_diag to TRUE in XIOS case 
-       DO naero = 1, naero_tot-1
-         IF (xios_field_is_active("dryod550_"//name_aero_tau(naero))) dryaod_diag=.TRUE. 
-       ENDDO
-       !
-       !--setting up ok_4xCO2atm to TRUE in XIOS case 
-       IF (xios_field_is_active("rsut4co2").OR.xios_field_is_active("rlut4co2").OR. & 
-           xios_field_is_active("rsutcs4co2").OR.xios_field_is_active("rlutcs4co2").OR. &
-           xios_field_is_active("rsu4co2").OR.xios_field_is_active("rsucs4co2").OR. &
-           xios_field_is_active("rsd4co2").OR.xios_field_is_active("rsdcs4co2").OR. &
-           xios_field_is_active("rlu4co2").OR.xios_field_is_active("rlucs4co2").OR. &
-           xios_field_is_active("rld4co2").OR.xios_field_is_active("rldcs4co2")) &
-           ok_4xCO2atm=.TRUE. 
+! Need to put this initialisation after phyetat0 as in the coupled model the XIOS context is only
+! initialised at that moment
+       ! Get "missing_val" value from XML files (from temperature variable)
+       !$OMP MASTER
+       CALL xios_get_field_attr("temp",default_value=missing_val_omp)
+       !$OMP END MASTER
+       !$OMP BARRIER
+       missing_val=missing_val_omp
 #endif 
 
@@ -1824,4 +1867,17 @@
        WRITE(*,*)'lat2_beta=',lat2_beta
        WRITE(*,*)'mskocean_beta=',mskocean_beta
+
+      !lwoff=y : offset LW CRE for radiation code and other schemes
+      !lwoff=y : betalwoff=1.
+      betalwoff=0.
+      IF (ok_lwoff) THEN
+         betalwoff=1.
+      ENDIF
+      WRITE(*,*)'ok_lwoff=',ok_lwoff
+      ! 
+      !lwoff=y to begin only sollw and sollwdown are set up to CS values
+      sollw = sollw + betalwoff * (sollw0 - sollw)
+      sollwdown(:)= sollwdown(:) + betalwoff *(-1.*ZFLDN0(:,1) - &
+                    sollwdown(:))
     ENDIF
     !
@@ -1890,18 +1946,23 @@
        ENDDO
     ENDDO
-    da(:,:)=0.
-    mp(:,:)=0.
-    phi(:,:,:)=0.
-    ! RomP >>>
-    phi2(:,:,:)=0.
     beta_prec_fisrt(:,:)=0.
     beta_prec(:,:)=0.
-    epmlmMm(:,:,:)=0.
-    eplaMm(:,:)=0.
-    d1a(:,:)=0.
-    dam(:,:)=0.
-    pmflxr=0.
-    pmflxs=0.
-    ! RomP <<<
+    !
+    !   Output variables from the convective scheme should not be set to 0 
+    !   since convection is not always called at every time step.
+    IF (ok_bug_cv_trac) THEN
+      da(:,:)=0.
+      mp(:,:)=0.
+      phi(:,:,:)=0.
+      ! RomP >>>
+      phi2(:,:,:)=0.
+      epmlmMm(:,:,:)=0.
+      eplaMm(:,:)=0.
+      d1a(:,:)=0.
+      dam(:,:)=0.
+      pmflxr(:,:)=0.
+      pmflxs(:,:)=0.
+      ! RomP <<<
+    ENDIF
 
     !
@@ -2080,10 +2141,11 @@
           IF(adjust_tropopause) THEN
              CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),   &
-                      ro3i,  press_edg_climoz,  paprs,   wo, time_climoz,    &
-                      longitude_deg,   latitude_deg,    press_cen_climoz,    &
+                      ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1),    &
+                      time_climoz ,  longitude_deg,   latitude_deg,          &
                       dyn_tropopause(t_seri, ztsol, paprs, pplay, rot))
           ELSE
-             CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),   &
-                      ro3i,  press_edg_climoz,  paprs,  wo,  time_climoz)
+             CALL regr_pr_time_av(ncid_climoz,  vars_climoz(1:read_climoz),  &
+                      ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1),    &
+                      time_climoz )
           END IF
           ! Convert from mole fraction of ozone to column density of ozone in a
@@ -2197,4 +2259,5 @@
        END SELECT
     ENDIF
+    sza_o = ACOS (rmu0) *180./pi
 
     IF (mydebug) THEN
@@ -2234,5 +2297,7 @@
 !!jyg       IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,2) .eq. 1) THEN
        IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,10) .ge. 1) THEN
-          print *,'debut du splitting de la PBL'
+          print *,'debut du splitting de la PBL, wake_s = ', wake_s(:)
+          print *,'debut du splitting de la PBL, wake_deltat = ', wake_deltat(:,1)
+          print *,'debut du splitting de la PBL, wake_deltaq = ', wake_deltaq(:,1)
        ENDIF
        ! !!
@@ -2246,7 +2311,7 @@
           gustiness(1:klon)=0
        ELSE IF (iflag_gusts==1) THEN
-          do i = 1, klon
-             gustiness(i)=f_gust_bl*ale_bl(i)+f_gust_wk*ale_wake(i)
-          enddo
+          gustiness(1:klon)=f_gust_bl*ale_bl(1:klon)+f_gust_wk*ale_wake(1:klon)
+       ELSE IF (iflag_gusts==2) THEN
+          gustiness(1:klon)=f_gust_bl*ale_bl_stat(1:klon)+f_gust_wk*ale_wake(1:klon)
           ! ELSE IF (iflag_gusts==2) THEN
           !    do i = 1, klon
@@ -2323,5 +2388,5 @@
           d_deltaq_vdf(:,:) = d_q_vdf_w(:,:)-d_q_vdf_x(:,:)
           CALL add_wake_tend &
-             (d_deltat_vdf, d_deltaq_vdf, dsig0, ddens0, wkoccur1, 'vdf', abortphy) 
+             (d_deltat_vdf, d_deltaq_vdf, dsig0, ddens0, ddens0, wkoccur1, 'vdf', abortphy) 
        ELSE
           d_deltat_vdf(:,:) = 0.
@@ -2465,6 +2530,33 @@
     ! Appel de la convection tous les "cvpas"
     !
-    IF (MOD(itapcv,cvpas).EQ.0) THEN
-
+!!jyg    IF (MOD(itapcv,cvpas).EQ.0) THEN
+!!    print *,' physiq : itapcv, cvpas, itap-1, cvpas_0 ', &
+!!                       itapcv, cvpas, itap-1, cvpas_0
+    IF (MOD(itapcv,cvpas).EQ.0 .OR. MOD(itap-1,cvpas_0).EQ.0) THEN
+
+    !
+    ! Mettre a zero des variables de sortie (pour securite)
+    !
+    pmflxr(:,:) = 0.
+    pmflxs(:,:) = 0.
+    wdtrainA(:,:) = 0.
+    wdtrainM(:,:) = 0.
+    upwd(:,:) = 0.
+    dnwd(:,:) = 0.
+    ep(:,:) = 0.
+    da(:,:)=0.
+    mp(:,:)=0.
+    wght_cvfd(:,:)=0.
+    phi(:,:,:)=0.
+    phi2(:,:,:)=0.
+    epmlmMm(:,:,:)=0.
+    eplaMm(:,:)=0.
+    d1a(:,:)=0.
+    dam(:,:)=0.
+    elij(:,:,:)=0.
+    ev(:,:)=0.
+    clw(:,:)=0.
+    sij(:,:,:)=0.
+    !
     IF (iflag_con.EQ.1) THEN
        abort_message ='reactiver le call conlmd dans physiq.F'
@@ -2539,5 +2631,5 @@
              IF (iflag_adjwk == 2) THEN
                CALL add_wake_tend &
-                 (d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, ddens0, wkoccur1, 'ajs_cv', abortphy) 
+                 (d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, ddens0, ddens0, wkoccur1, 'ajs_cv', abortphy) 
              ENDIF  ! (iflag_adjwk == 2)
           ENDIF  ! (iflag_adjwk >= 1)
@@ -2568,8 +2660,13 @@
           ! Calculate the upmost level of deep convection loops: k_upper_cv
           !  (near 22 km)
-          izero = klon/2+1/klon
           k_upper_cv = klev
+          !izero = klon/2+1/klon
+          !DO k = klev,1,-1
+          !   IF (pphi(izero,k) > 22.e4) k_upper_cv = k
+          !ENDDO
+          ! FH : nouveau calcul base sur un profil global sans quoi
+          ! le modele etait sensible au decoupage de domaines
           DO k = klev,1,-1
-             IF (pphi(izero,k) > 22.e4) k_upper_cv = k
+             IF (-7*log(presnivs(k)/presnivs(1)) > 25.) k_upper_cv = k
           ENDDO
           IF (prt_level .ge. 5) THEN
@@ -2615,7 +2712,38 @@
           clwcon0=qcondc
           pmfu(:,:)=upwd(:,:)+dnwd(:,:)
-
+          !
+          !jyg<
+          ! If convective tendencies are too large, then call convection 
+          !  every time step
+          cvpas = cvpas_0
+          DO k=1,k_upper_cv
+             DO i=1,klon
+               IF (d_t_con(i,k) > 6.721 .AND. d_t_con(i,k) < 6.722 .AND.&
+                   d_q_con(i,k) > -.0002171 .AND. d_q_con(i,k) < -.0002170) THEN
+                     dtcon_multistep_max = 3.
+                     dqcon_multistep_max = 0.02
+               ENDIF
+             ENDDO
+          ENDDO
+!
+          DO k=1,k_upper_cv
+             DO i=1,klon
+!!               IF (abs(d_t_con(i,k)) > 0.24 .OR. &
+!!                   abs(d_q_con(i,k)) > 2.e-2) THEN
+               IF (abs(d_t_con(i,k)) > dtcon_multistep_max .OR. &
+                   abs(d_q_con(i,k)) > dqcon_multistep_max) THEN
+                 cvpas = 1
+!!                 print *,'physiq1, i,k,d_t_con(i,k),d_q_con(i,k) ', &
+!!                                   i,k,d_t_con(i,k),d_q_con(i,k)
+               ENDIF
+             ENDDO
+          ENDDO
+!!!   Ligne a ne surtout pas remettre sans avoir murement reflechi (jyg)
+!!!          call bcast(cvpas)
+!!!   ------------------------------------------------------------
+          !>jyg
+          !
           DO i = 1, klon
-             IF (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
+             IF (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+cvpas
           ENDDO
           !
@@ -2725,7 +2853,15 @@
     proba_notrig(:) = 1.
     itapcv = 0
-    ENDIF !  (MOD(itapcv,cvpas).EQ.0)
+    ENDIF !  (MOD(itapcv,cvpas).EQ.0 .OR. MOD(itapcv,cvpas_0).EQ.0)
 !
     itapcv = itapcv+1
+    !
+    ! Compter les steps ou cvpas=1
+    IF (cvpas == 1) THEN
+      Ncvpaseq1 = Ncvpaseq1+1
+    ENDIF
+    IF (mod(itap,1000) == 0) THEN
+      print *,' physiq, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
+    ENDIF
 
 !!!jyg  Appel diagnostique a add_phys_tend pour tester la conservation de 
@@ -2849,7 +2985,7 @@
                t_seri, q_seri, omega,  &
                dt_dwn, dq_dwn, M_dwn, M_up,  &
-               dt_a, dq_a,  &
-               sigd,  &
-               wake_deltat, wake_deltaq, wake_s, wake_dens,  &
+               dt_a, dq_a, cv_gen,  &
+               sigd, cin,  &
+               wake_deltat, wake_deltaq, wake_s, awake_dens, wake_dens,  &
                wake_dth, wake_h,  &
 !!               wake_pe, wake_fip, wake_gfl,  &
@@ -2861,5 +2997,5 @@
                wake_omg, wake_dp_deltomg,  &
                wake_spread, wake_Cstar, d_deltat_wk_gw,  &
-               d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_wk)
+               d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_a_wk, d_dens_wk)
           !
           !jyg    Reinitialize itapwk when wakes have been called
@@ -2880,8 +3016,14 @@
 
          CALL add_wake_tend &
-            (d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_wk, wake_k, &
+            (d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_a_wk, d_dens_wk, wake_k, &
              'wake', abortphy) 
           call prt_enerbil('wake',itap)
        ENDIF   ! (iflag_wake_tend .GT. 0.)
+       !
+       IF (prt_level .GE. 10) THEN
+         print *,' physiq, after calwake, wake_s: ',wake_s(:)
+         print *,' physiq, after calwake, wake_deltat: ',wake_deltat(:,1)
+         print *,' physiq, after calwake, wake_deltaq: ',wake_deltaq(:,1)
+       ENDIF
 
        IF (iflag_alp_wk_cond .GT. 0.) THEN
@@ -3011,6 +3153,11 @@
              ENDDO
           !
-             CALL add_wake_tend &
-                 (d_deltat_the, d_deltaq_the, dsig0, ddens0, wkoccur1, 'the', abortphy) 
+             IF (ok_bug_split_th) THEN
+               CALL add_wake_tend &
+                   (d_deltat_the, d_deltaq_the, dsig0, ddens0, ddens0, wkoccur1, 'the', abortphy) 
+             ELSE
+               CALL add_wake_tend &
+                   (d_deltat_the, d_deltaq_the, dsig0, ddens0, ddens0, wake_k, 'the', abortphy) 
+             ENDIF
              call prt_enerbil('the',itap)
           !
@@ -3026,5 +3173,5 @@
                           ale_bl_trig, ale_bl_stat, ale_bl,  &
                           alp_bl, alp_bl_stat, &
-                          proba_notrig, random_notrig)
+                          proba_notrig, random_notrig, cv_gen)
           !>jyg
 
@@ -3488,4 +3635,10 @@
                         tausum_aero, drytausum_aero, tau3d_aero)
 #endif
+
+                   IF (flag_aerosol .EQ. 7) THEN
+                      CALL MACv2SP(pphis,pplay,paprs,longitude_deg,latitude_deg,  &
+                                   tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm)
+                   ENDIF
+
                    !
                 ELSE IF (NSW.EQ.2) THEN 
@@ -3542,5 +3695,5 @@
        !
        !--WMO criterion to determine tropopause
-       CALL stratosphere_mask(missing_val, t_seri, pplay, latitude_deg)
+       CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
        !
        !--STRAT AEROSOL
@@ -3582,5 +3735,5 @@
 #ifdef CPP_StratAer
        !--compute stratospheric mask
-       CALL stratosphere_mask(missing_val, t_seri, pplay, latitude_deg)
+       CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
        !--interactive strat aerosols
        CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut)
@@ -3612,5 +3765,5 @@
 #endif
           ENDIF
-          CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, &
+          CALL newmicro (flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, &
                paprs, pplay, t_seri, cldliq, cldfra, &
                cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
@@ -3642,4 +3795,7 @@
           ! global
           !
+!IM 251017 begin
+!               print*,'physiq betaCRF global zdtime=',zdtime
+!IM 251017 end
           DO k=1, klev
              DO i=1, klon
@@ -3756,4 +3912,5 @@
              print *,' ->radlwsw, number 1 '
           ENDIF
+
           !
           CALL radlwsw &
@@ -3772,4 +3929,5 @@
                tau_aero_lw_rrtm, & 
                cldtaupirad,new_aod, &
+!              zqsat, flwcrad, fiwcrad, &
                zqsat, flwc, fiwc, &
                ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
@@ -3778,6 +3936,6 @@
                sollwdown, &
                topsw0,toplw0,solsw0,sollw0, &
-               lwdn0, lwdn, lwup0, lwup,  &
-               swdn0, swdn, swup0, swup, &
+               lwdnc0, lwdn0, lwdn, lwupc0, lwup0, lwup,  &
+               swdnc0, swdn0, swdn, swupc0, swup0, swup, &
                topswad_aero, solswad_aero, &
                topswai_aero, solswai_aero, &
@@ -3794,4 +3952,14 @@
                ZSWFT0_i, ZFSDN0, ZFSUP0)
 
+          !lwoff=y, betalwoff=1. : offset LW CRE for radiation code and other
+          !schemes
+          toplw = toplw + betalwoff * (toplw0 - toplw)
+          sollw = sollw + betalwoff * (sollw0 - sollw)
+          lwdn = lwdn + betalwoff * (lwdn0 - lwdn)
+          lwup = lwup + betalwoff * (lwup0 - lwup)
+          sollwdown(:)= sollwdown(:) + betalwoff *(-1.*ZFLDN0(:,1) - &
+                        sollwdown(:))
+          cool = cool + betalwoff * (cool0 - cool)
+ 
 #ifndef CPP_XIOS
           !--OB 30/05/2016 modified 21/10/2016
@@ -3800,4 +3968,5 @@
           !--this is necessary to get the right swaero at first step
           !--but only in the case of no XIOS as XIOS is covered elsewhere 
+          IF (debut) swaerofree_diag = .FALSE.
           IF (debut) swaero_diag = .FALSE.
           IF (debut) dryaod_diag = .FALSE.
@@ -3805,5 +3974,5 @@
           !--as for swaero_diag, see above
           IF (debut) ok_4xCO2atm = .FALSE.
-#endif
+
           !
           !IM 2eme calcul radiatif pour le cas perturbe ou au moins un
@@ -3816,4 +3985,5 @@
               RCFC11_per.NE.RCFC11_act.OR. &
               RCFC12_per.NE.RCFC12_act) ok_4xCO2atm =.TRUE. 
+#endif
    ! 
           IF (ok_4xCO2atm) THEN
@@ -3844,4 +4014,5 @@
                      tau_aero_lw_rrtm, &
                      cldtaupi,new_aod, &
+!                    zqsat, flwcrad, fiwcrad, &
                      zqsat, flwc, fiwc, &
                      ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
@@ -3850,6 +4021,6 @@
                      sollwdownp, &
                      topsw0p,toplw0p,solsw0p,sollw0p, &
-                     lwdn0p, lwdnp, lwup0p, lwupp,  &
-                     swdn0p, swdnp, swup0p, swupp, &
+                     lwdnc0p, lwdn0p, lwdnp, lwupc0p, lwup0p, lwupp,  &
+                     swdnc0p, swdn0p, swdnp, swupc0p, swup0p, swupp, &
                      topswad_aerop, solswad_aerop, &
                      topswai_aerop, solswai_aerop, &
@@ -3865,6 +4036,5 @@
                      ZLWFT0_i, ZFLDN0, ZFLUP0, &
                      ZSWFT0_i, ZFSDN0, ZFSUP0)
-          endif
-          !
+          endif !ok_4xCO2atm
        ENDIF ! aerosol_couple
        itaprad = 0
@@ -4257,5 +4427,5 @@
 
 
-    CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pbl_tke)
+    CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pctsrf,pbl_tke)
 
 
@@ -4412,6 +4582,6 @@
     !
     CALL transp (paprs,zxtsol, &
-         t_seri, q_seri, u_seri, v_seri, zphi, &
-         ve, vq, ue, uq)
+         t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zphi, &
+         ve, vq, ue, uq, vwat, uwat)
     !
     !IM global posePB BEG
@@ -4709,4 +4879,53 @@
 #endif
 
+! On remet des variables a .false. apres un premier appel
+    if (debut) then
+#ifdef CPP_XIOS
+      swaero_diag=.FALSE.
+      swaerofree_diag=.FALSE.
+      dryaod_diag=.FALSE.
+      ok_4xCO2atm= .FALSE.
+!      write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
+
+      IF (is_master) then
+        !--setting up swaero_diag to TRUE in XIOS case 
+        IF (xios_field_is_active("topswad").OR.xios_field_is_active("topswad0").OR. & 
+           xios_field_is_active("solswad").OR.xios_field_is_active("solswad0").OR. & 
+           xios_field_is_active("topswai").OR.xios_field_is_active("solswai").OR.  & 
+             (iflag_rrtm==1.AND.(xios_field_is_active("toplwad").OR.xios_field_is_active("toplwad0").OR. & 
+                                 xios_field_is_active("sollwad").OR.xios_field_is_active("sollwad0"))))  & 
+           !!!--for now these fields are not in the XML files so they are omitted 
+           !!!  xios_field_is_active("toplwai").OR.xios_field_is_active("sollwai") !))) & 
+           swaero_diag=.TRUE. 
+
+        !--setting up swaerofree_diag to TRUE in XIOS case 
+        IF (xios_field_is_active("SWdnSFCcleanclr").OR.xios_field_is_active("SWupSFCcleanclr").OR. &
+           xios_field_is_active("SWupTOAcleanclr").OR.xios_field_is_active("rsucsaf").OR.   &
+           xios_field_is_active("rsdcsaf") .OR. xios_field_is_active("LWdnSFCcleanclr").OR. &
+           xios_field_is_active("LWupTOAcleanclr")) &
+           swaerofree_diag=.TRUE. 
+
+        !--setting up dryaod_diag to TRUE in XIOS case 
+        DO naero = 1, naero_tot-1
+         IF (xios_field_is_active("dryod550_"//name_aero_tau(naero))) dryaod_diag=.TRUE. 
+        ENDDO
+        !
+        !--setting up ok_4xCO2atm to TRUE in XIOS case 
+        IF (xios_field_is_active("rsut4co2").OR.xios_field_is_active("rlut4co2").OR. & 
+           xios_field_is_active("rsutcs4co2").OR.xios_field_is_active("rlutcs4co2").OR. &
+           xios_field_is_active("rsu4co2").OR.xios_field_is_active("rsucs4co2").OR. &
+           xios_field_is_active("rsd4co2").OR.xios_field_is_active("rsdcs4co2").OR. &
+           xios_field_is_active("rlu4co2").OR.xios_field_is_active("rlucs4co2").OR. &
+           xios_field_is_active("rld4co2").OR.xios_field_is_active("rldcs4co2")) &
+           ok_4xCO2atm=.TRUE. 
+      endif
+      !$OMP BARRIER
+      call bcast(swaero_diag)
+      call bcast(swaerofree_diag)
+      call bcast(dryaod_diag)
+      call bcast(ok_4xCO2atm)
+!      write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
+#endif
+    endif
 
     !====================================================================
@@ -4745,4 +4964,5 @@
        IF (is_omp_master) CALL xios_context_finalize
 #endif
+       print *,' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
     ENDIF
 
