Index: LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 3927)
+++ LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 3940)
@@ -1,4 +1,4 @@
 !
-! $Id: physiq_mod.F90 3666 2020-04-20 10:13:34Z lfalletti $
+! $Id: physiq_mod.F90 3908 2021-05-20 07:11:13Z idelkadi $
 !
 !#define IO_DEBUG
@@ -16,14 +16,108 @@
        d_u, d_v, d_t, d_qx, d_ps)
 
+! For clarity, the "USE" section is now arranged in alphabetical order,
+! with a separate section for CPP keys
+! PLEASE try to follow this rule 
+
+    USE ACAMA_GWD_rando_m, only: ACAMA_GWD_rando
+    USE aero_mod
+    USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, diag_phys_tend, prt_enerbil, &
+  &      fl_ebil, fl_cor_ebil
     USE assert_m, only: assert
+    USE change_srf_frac_mod
+    USE conf_phys_m, only: conf_phys
+    USE carbon_cycle_mod, ONLY : infocfields_init, RCO2_glo, carbon_cycle_rad
+    USE CFMIP_point_locations   ! IM stations CFMIP
+    USE cmp_seri_mod
+    USE dimphy
+    USE etat0_limit_unstruct_mod
+    USE FLOTT_GWD_rando_m, only: FLOTT_GWD_rando
+    USE fonte_neige_mod, ONLY  : fonte_neige_get_vars
+    USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg
     USE ioipsl, only: histbeg, histvert, histdef, histend, histsync, &
          histwrite, ju2ymds, ymds2ju, getin
-    USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg
+    USE ioipsl_getin_p_mod, ONLY : getin_p
+    USE indice_sol_mod
+    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, type_trac,ok_isotopes, &
+        nqtottr,itr_indice ! C Risi
+
+    USE iophy
+    USE limit_read_mod, ONLY : init_limit_read
+    USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo, grid1dTo2d_glo, grid_type, unstructured
+    USE mod_phys_lmdz_mpi_data, only: is_mpi_root
+    USE mod_phys_lmdz_para
+    USE netcdf95, only: nf95_close
+    USE netcdf, only: nf90_fill_real     ! IM for NMC files
+    USE open_climoz_m, only: open_climoz ! ozone climatology from a file
+    USE ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
+    USE pbl_surface_mod, ONLY : pbl_surface
+    USE phyaqua_mod, only: zenang_an
+    USE phystokenc_mod, ONLY: offline, phystokenc
     USE phys_cal_mod, only: year_len, mth_len, days_elapsed, jh_1jan, &
          year_cur, mth_cur,jD_cur, jH_cur, jD_ref, day_cur, hour
+!!  USE phys_local_var_mod, ONLY : a long list of variables
+!!              ==> see below, after "CPP Keys" section
+    USE phys_state_var_mod ! Variables sauvegardees de la physique
+    USE phys_output_mod
+    USE phys_output_ctrlout_mod
+    USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
+    USE readaerosol_mod, ONLY : init_aero_fromfile
+    USE readaerosolstrato_m, ONLY : init_readaerosolstrato
+    USE radlwsw_m, only: radlwsw
+    USE regr_horiz_time_climoz_m, ONLY: regr_horiz_time_climoz
+    USE regr_pr_time_av_m, only: regr_pr_time_av
+    USE surface_data,     ONLY : type_ocean, ok_veget, landice_opt
+    USE time_phylmdz_mod, only: annee_ref, current_time, day_ini, day_ref, &
+          day_step_phy, itau_phy, pdtphys, raz_date, start_time, update_time
+    USE tracinca_mod, ONLY: config_inca
+    USE tropopause_m,     ONLY: dyn_tropopause
+    USE vampir
+    USE VERTICAL_LAYERS_MOD, ONLY: aps,bps, ap, bp
     USE write_field_phy
-    USE dimphy
-USE infotrac_phy, ONLY: nqtot, nbtr, nqo, type_trac,ok_isotopes, &
-        nqtottr,itr_indice ! C Risi
+
+    !USE cmp_seri_mod
+!    USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, diag_phys_tend, prt_enerbil, &
+!  &      fl_ebil, fl_cor_ebil
+
+!!!!!!!!!!!!!!!!!! "USE" section for CPP keys !!!!!!!!!!!!!!!!!!!!!!!!
+!
+!
+#ifdef CPP_Dust
+    USE phytracr_spl_mod, ONLY: phytracr_spl, phytracr_spl_out_init
+    USE phys_output_write_spl_mod
+#else
+    USE phytrac_mod, ONLY : phytrac_init, phytrac
+    USE phys_output_write_mod
+#endif
+
+
+#ifdef REPROBUS
+    USE CHEM_REP, ONLY : Init_chem_rep_xjour, &
+         d_q_rep,d_ql_rep,d_qi_rep,ptrop,ttrop, &
+         ztrop, gravit,itroprep, Z1,Z2,fac,B
+#endif
+
+
+#ifdef CPP_RRTM
+    USE YOERAD, ONLY : NRADLP
+    USE YOESW, ONLY : RSUN
+#endif
+
+
+#ifdef CPP_StratAer
+    USE strataer_mod, ONLY: strataer_init
+#endif
+
+
+#ifdef CPP_XIOS
+    USE xios, ONLY: xios_update_calendar, xios_context_finalize, &
+            xios_get_field_attr, xios_field_is_active
+    USE wxios, ONLY: missing_val, missing_val_omp
+#endif
+#ifndef CPP_XIOS
+    USE paramLMDZ_phy_mod
+#endif
+
+
 #ifdef ISO
     USE infotrac_phy, ONLY:  &
@@ -64,22 +158,10 @@
 #endif
 #endif
-    USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo, grid1dTo2d_glo, grid_type, unstructured
-    USE mod_phys_lmdz_para
-    USE iophy
-    USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
-    USE phystokenc_mod, ONLY: offline, phystokenc
-    USE time_phylmdz_mod, only: raz_date, day_step_phy, update_time,current_time
-    USE vampir
-    USE pbl_surface_mod, ONLY : pbl_surface
-    USE change_srf_frac_mod
-    USE surface_data,     ONLY : type_ocean, ok_veget, ok_snow
-    USE tropopause_m,     ONLY: dyn_tropopause
-#ifdef CPP_Dust
-    USE phytracr_spl_mod, ONLY: phytracr_spl
-#endif
-#ifdef CPP_StratAer
-    USE strataer_mod, ONLY: strataer_init
-#endif
-    USE phys_local_var_mod, ONLY: phys_local_var_init, phys_local_var_end, &
+
+!
+!
+!!!!!!!!!!!!!!!!!!  END "USE" for CPP keys !!!!!!!!!!!!!!!!!!!!!!
+
+USE phys_local_var_mod, ONLY: phys_local_var_init, phys_local_var_end, &
        ! [Variables internes non sauvegardees de la physique]
        ! Variables locales pour effectuer les appels en serie
@@ -160,7 +242,5 @@
        cdragm, cdragh,                   &
        zustar, zu10m, zv10m, rh2m, qsat2m, &
-       zq2m, zt2m, weak_inversion, &
-       zq2m_cor,zt2m_cor,zu10m_cor,zv10m_cor, & ! pour corriger d'un bug
-       zrh2m_cor,zqsat2m_cor, &
+       zq2m, zt2m, zn2mout, weak_inversion, &
        zt2m_min_mon, zt2m_max_mon,   &         ! pour calcul_divers.h
        t2m_min_mon, t2m_max_mon,  &            ! pour calcul_divers.h
@@ -175,4 +255,5 @@
        zxrunofflic,                            &
        zxtsol, snow_lsc, zxfqfonte, zxqsurf,   &
+       delta_qsurf,                            &
        rain_lsc, rain_num,                     &
        !
@@ -236,5 +317,5 @@
        ref_liq, ref_ice, theta,  &
        ref_liq_pi, ref_ice_pi,  &
-       zphi, zx_rh,  &
+       zphi, zx_rh, zx_rhl, zx_rhi,  &
        pmfd, pmfu,  &
        !
@@ -284,69 +365,5 @@
 #endif
        !
-    USE phys_state_var_mod ! Variables sauvegardees de la physique
-#ifdef CPP_Dust
-    USE phys_output_write_spl_mod
-#else
-    USE phys_output_var_mod ! Variables pour les ecritures des sorties
-#endif
-
-    USE phys_output_write_mod
-    USE fonte_neige_mod, ONLY  : fonte_neige_get_vars
-    USE phys_output_mod
-    USE phys_output_ctrlout_mod
-    USE open_climoz_m, only: open_climoz ! ozone climatology from a file
-    USE regr_pr_time_av_m, only: regr_pr_time_av
-    USE netcdf95, only: nf95_close
-    !IM for NMC files
-    USE netcdf, only: nf90_fill_real
-    USE mod_phys_lmdz_mpi_data, only: is_mpi_root
-    USE aero_mod
-    USE ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
-    USE conf_phys_m, only: conf_phys
-    USE radlwsw_m, only: radlwsw
-    USE phyaqua_mod, only: zenang_an
-    USE time_phylmdz_mod, only: day_step_phy, annee_ref, day_ref, itau_phy, &
-         start_time, pdtphys, day_ini
-    USE tracinca_mod, ONLY: config_inca
-#ifdef CPP_XIOS
-    USE wxios, ONLY: missing_val, missing_val_omp
-    USE xios, ONLY: xios_get_field_attr, xios_field_is_active
-#endif
-#ifdef REPROBUS
-    USE CHEM_REP, ONLY : Init_chem_rep_xjour, &
-         d_q_rep,d_ql_rep,d_qi_rep,ptrop,ttrop, &
-         ztrop, gravit,itroprep, Z1,Z2,fac,B
-#endif
-    USE indice_sol_mod
-    USE phytrac_mod, ONLY : phytrac_init, phytrac
-    USE carbon_cycle_mod, ONLY : infocfields_init, RCO2_glo, carbon_cycle_rad
-
-#ifdef CPP_RRTM
-    USE YOERAD, ONLY : NRADLP
-    USE YOESW, ONLY : RSUN
-#endif
-    USE ioipsl_getin_p_mod, ONLY : getin_p
-
-#ifndef CPP_XIOS
-    USE paramLMDZ_phy_mod
-#endif
-
-    USE cmp_seri_mod
-    USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, diag_phys_tend, prt_enerbil, &
-  &      fl_ebil, fl_cor_ebil
-
-    !IM stations CFMIP
-    USE CFMIP_point_locations
-    USE FLOTT_GWD_rando_m, only: FLOTT_GWD_rando
-    USE ACAMA_GWD_rando_m, only: ACAMA_GWD_rando
-    USE VERTICAL_LAYERS_MOD, ONLY: aps,bps, ap, bp
-    USE etat0_limit_unstruct_mod
-#ifdef CPP_XIOS
-    USE xios, ONLY: xios_update_calendar, xios_context_finalize
-#endif
-    USE limit_read_mod, ONLY : init_limit_read
-    USE regr_horiz_time_climoz_m, ONLY: regr_horiz_time_climoz
-    USE readaerosol_mod, ONLY : init_aero_fromfile
-    USE readaerosolstrato_m, ONLY : init_readaerosolstrato
+
 
     IMPLICIT NONE
@@ -671,9 +688,4 @@
     !$OMP THREADPRIVATE(iflag_alp_wk_cond)
 
-    INTEGER,  SAVE               :: iflag_bug_t2m_ipslcm61=1 !
-    !$OMP THREADPRIVATE(iflag_bug_t2m_ipslcm61)
-    INTEGER,  SAVE               :: iflag_bug_t2m_stab_ipslcm61=-1 !
-    !$OMP THREADPRIVATE(iflag_bug_t2m_stab_ipslcm61)
-
     REAL t_w(klon,klev),q_w(klon,klev) ! temperature and moisture profiles in the wake region
     REAL t_x(klon,klev),q_x(klon,klev) ! temperature and moisture profiles in the off-wake region
@@ -1419,9 +1431,4 @@
        WRITE(lunout,*) 'debut physiq_mod tau_gl=',tau_gl
 
-       iflag_bug_t2m_ipslcm61 = 1
-       CALL getin_p('iflag_bug_t2m_ipslcm61', iflag_bug_t2m_ipslcm61)
-       iflag_bug_t2m_stab_ipslcm61 = -1
-       CALL getin_p('iflag_bug_t2m_stab_ipslcm61', iflag_bug_t2m_stab_ipslcm61)
-
        CALL getin_p('iflag_alp_wk_cond', iflag_alp_wk_cond) 
        CALL getin_p('random_notrig_max',random_notrig_max)
@@ -1451,4 +1458,10 @@
        iflag_phytrac = 1 ! by default we do want to call phytrac
        CALL getin_p('iflag_phytrac',iflag_phytrac)
+#ifdef CPP_Dust
+       IF (iflag_phytrac.EQ.0) THEN 
+         WRITE(lunout,*) 'In order to run with SPLA, iflag_phytrac will be forced to 1'
+         iflag_phytrac = 1
+       ENDIF
+#endif
        nvm_lmdz = 13
        CALL getin_p('NVM',nvm_lmdz)
@@ -1508,5 +1521,5 @@
        tau_overturning_th(:)=0.
 
-       IF (type_trac == 'inca') THEN
+       IF (type_trac == 'inca' .OR. type_trac == 'inco') THEN
           ! jg : initialisation jusqu'au ces variables sont dans restart
           ccm(:,:,:) = 0.
@@ -1763,9 +1776,9 @@
 #ifdef CPP_COSP
       IF (ok_cosp) THEN
-           DO k = 1, klev
-             DO i = 1, klon
-               phicosp(i,k) = pphi(i,k) + pphis(i)
-             ENDDO
-           ENDDO
+!           DO k = 1, klev
+!             DO i = 1, klon
+!               phicosp(i,k) = pphi(i,k) + pphis(i)
+!             ENDDO
+!           ENDDO
         CALL phys_cosp(itap,phys_tstep,freq_cosp, &
                ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
@@ -1785,9 +1798,9 @@
 #ifdef CPP_COSP2
         IF (ok_cosp) THEN
-           DO k = 1, klev
-             DO i = 1, klon
-               phicosp(i,k) = pphi(i,k) + pphis(i)
-             ENDDO
-           ENDDO
+!           DO k = 1, klev
+!             DO i = 1, klon
+!               phicosp(i,k) = pphi(i,k) + pphis(i)
+!             ENDDO
+!           ENDDO
           CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
                ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
@@ -1835,9 +1848,21 @@
 
        CALL iniradia(klon,klev,paprs(1,1:klev+1))
-       ! Initialisation des champs dans phytrac qui sont utilisés par phys_output_write
+
+       ! Initialisation des champs dans phytrac* qui sont utilisés par phys_output_write*
+#ifdef CPP_Dust
+       ! Quand on utilise SPLA, on force iflag_phytrac=1
+       CALL phytracr_spl_out_init()
+       CALL phys_output_write_spl(itap, pdtphys, paprs, pphis,                  &
+                                pplay, lmax_th, aerosol_couple,                 &
+                                ok_ade, ok_aie, ivap, ok_sync,                  &
+                                ptconv, read_climoz, clevSTD,                   &
+                                ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse,      &
+                                flag_aerosol, flag_aerosol_strat, ok_cdnc)
+#else
+       ! phys_output_write écrit des variables traceurs seulement si iflag_phytrac == 1 
+       ! donc seulement dans ce cas on doit appeler phytrac_init()
        IF (iflag_phytrac == 1 ) THEN
           CALL phytrac_init()
-        ENDIF
-
+       ENDIF
        CALL phys_output_write(itap, pdtphys, paprs, pphis,                    &
                               pplay, lmax_th, aerosol_couple,                 &
@@ -1846,4 +1871,6 @@
                               ptconvth, d_u, d_t, qx, d_qx, zmasse,           &
                               flag_aerosol, flag_aerosol_strat, ok_cdnc)
+#endif
+
 
 #ifdef CPP_XIOS
@@ -2070,5 +2097,5 @@
        !c         ENDDO
        !
-       IF (type_trac == 'inca') THEN
+       IF (type_trac == 'inca' .OR. type_trac == 'inco') THEN			! ModThL
 #ifdef INCA
           CALL VTe(VTphysiq)
@@ -2087,5 +2114,5 @@
                klon, &
                nqtot, &
-               nqo, &
+               nqo+nqCO2, &
                pdtphys, &
                annee_ref, &
@@ -2840,5 +2867,6 @@
     !   s_therm,   s_trmb1,   s_trmb2, s_trmb3,
     !   zu10m,     zv10m,   fder,
-    !   zxqsurf,   rh2m,      zxfluxu, zxfluxv,
+    !   zxqsurf,   delta_qsurf,
+    !   rh2m,      zxfluxu, zxfluxv,
     !   frugs,     agesno,    fsollw,  fsolsw,
     !   d_ts,      fevap,     fluxlat, t2m,
@@ -2890,6 +2918,6 @@
             debut,     lafin, &
             longitude_deg, latitude_deg, rugoro,  zrmu0,      &
-            zsig,      sollwdown, pphi,    cldt,      &
-            rain_fall, snow_fall, solsw,   sollw,     &
+            sollwdown,    cldt,      &
+            rain_fall, snow_fall, solsw,   solswfdiff, sollw,     &
             gustiness,                                &
             t_seri,    q_seri,    u_seri,  v_seri,    &
@@ -2901,4 +2929,5 @@
                                 !albedo SB <<<
             cdragh,    cdragm,  u1,    v1,            &
+            beta_aridity, &
                                 !albedo SB >>>
                                 ! albsol1,   albsol2,   sens,    evap,      &
@@ -2906,5 +2935,5 @@
                                 !albedo SB <<<
             albsol3_lic,runoff,   snowhgt,   qsnow, to_ice, sissnow, &
-            zxtsol,    zxfluxlat, zt2m,    qsat2m,  &
+            zxtsol,    zxfluxlat, zt2m,    qsat2m,  zn2mout, &
             d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf, d_t_diss, &
                                 !nrlmd<
@@ -2927,5 +2956,5 @@
             s_therm,   s_trmb1,   s_trmb2, s_trmb3, &
             zustar, zu10m,     zv10m,   fder, &
-            zxqsurf,   rh2m,      zxfluxu, zxfluxv, &
+            zxqsurf, delta_qsurf,   rh2m,      zxfluxu, zxfluxv, &
             z0m, z0h,     agesno,    fsollw,  fsolsw, &
             d_ts,      fevap,     fluxlat, t2m, &
@@ -3041,21 +3070,4 @@
 #endif
        ENDIF
-
-!add limitation for t,q at and wind at 10m
-        if ( iflag_bug_t2m_ipslcm61 == 0 ) THEN
-          CALL borne_var_surf( klon,klev,nbsrf,                 &
-            iflag_bug_t2m_stab_ipslcm61,                        &
-            t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),    &
-            ftsol,zxqsurf,pctsrf,paprs,                         &
-            t2m, q2m, u10m, v10m,                               &
-            zt2m_cor, zq2m_cor, zu10m_cor, zv10m_cor,           &
-            zrh2m_cor, zqsat2m_cor)
-        ELSE
-          zt2m_cor(:)=zt2m(:)
-          zq2m_cor(:)=zq2m(:)
-          zu10m_cor(:)=zu10m(:)
-          zv10m_cor(:)=zv10m(:)
-          zqsat2m_cor=999.999
-        ENDIF
 
        !---------------------------------------------------------------------
@@ -3977,5 +3989,5 @@
              ENDDO
           ENDIF
-        
+         
 #ifdef ISOVERIF
         write(*,*) 'physiq 3977: verif des inputs de calwake'
@@ -4312,5 +4324,4 @@
           !
           ENDIF  ! (mod(iflag_pbl_split/10,10) .GE. 1)
-
           !
           CALL add_phys_tend(d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs,  &
@@ -4456,9 +4467,8 @@
          iflag_ratqs,iflag_con,iflag_cld_th,pdtphys,  &
          ratqsbas,ratqshaut,ratqsp0, ratqsdp, &
-         tau_ratqs,fact_cldcon,   &
+         tau_ratqs,fact_cldcon,wake_s, wake_deltaq,   &
          ptconv,ptconvth,clwcon0th, rnebcon0th,     &
          paprs,pplay,q_seri,zqsat,fm_therm, &
-         ratqs,ratqsc)
-
+         ratqs,ratqsc,ratqs_inter)
 
     !
@@ -5013,4 +5023,8 @@
           ENDIF
           zx_rh(i,k) = q_seri(i,k)/zx_qs
+            IF (iflag_ice_thermo .GT. 0) THEN
+          zx_rhl(i,k) = q_seri(i,k)/(qsatl(zx_t)/pplay(i,k))
+          zx_rhi(i,k) = q_seri(i,k)/(qsats(zx_t)/pplay(i,k))
+            ENDIF
           zqsat(i,k)=zx_qs
        ENDDO
@@ -5039,5 +5053,5 @@
     ENDDO
 
-    IF (type_trac == 'inca') THEN
+    IF (type_trac == 'inca' .OR. type_trac == 'inco') THEN	! ModThL
 #ifdef INCA
        CALL VTe(VTphysiq)
@@ -5082,5 +5096,5 @@
             nbp_lon, &
             nbp_lat-1, &
-            tr_seri, &
+            tr_seri(:,:,1+nqCO2:nbtr), &
             ftsol, &
             paprs, &
@@ -5093,6 +5107,6 @@
        CALL VTe(VTinca)
        CALL VTb(VTphysiq)
-#endif 
-    ENDIF !type_trac = inca
+#endif
+    ENDIF !type_trac = inca or inco
     IF (type_trac == 'repr') THEN
 #ifdef REPROBUS
@@ -5265,5 +5279,6 @@
 
        IF (ok_newmicro) then
-          IF (iflag_rrtm.NE.0) THEN 
+! AI          IF (iflag_rrtm.NE.0) THEN 
+          IF (iflag_rrtm.EQ.1) THEN
 #ifdef CPP_RRTM
              IF (ok_cdnc.AND.NRADLP.NE.3) THEN
@@ -5452,5 +5467,5 @@
                heat,heat0,cool,cool0,albpla, &
                heat_volc,cool_volc, &
-               topsw,toplw,solsw,sollw, &
+               topsw,toplw,solsw,solswfdiff,sollw, &
                sollwdown, &
                topsw0,toplw0,solsw0,sollw0, &
@@ -5539,5 +5554,5 @@
                      heatp,heat0p,coolp,cool0p,albplap, &
                      heat_volc,cool_volc, &
-                     topswp,toplwp,solswp,sollwp, &
+                     topswp,toplwp,solswp,solswfdiffp,sollwp, &
                      sollwdownp, &
                      topsw0p,toplw0p,solsw0p,sollw0p, &
@@ -5976,5 +5991,8 @@
 
     CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pctsrf,pbl_tke)
-
+   !
+   ! Prevent pbl_tke_w from becoming negative
+    wake_delta_pbl_tke(:,:,:) = max(wake_delta_pbl_tke(:,:,:), -pbl_tke(:,:,:))
+   !
 
        ENDIF
@@ -6090,8 +6108,14 @@
 #ifdef CPP_COSPV2
        IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
+!        IF (MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
 
           IF (prt_level .GE.10) THEN
              print*,'freq_cosp',freq_cosp
           ENDIF
+           DO k = 1, klev
+             DO i = 1, klon
+               phicosp(i,k) = pphi(i,k) + pphis(i)
+             ENDDO
+           ENDDO
           mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
                  print*,'Dans physiq.F avant appel '
@@ -6194,11 +6218,14 @@
     ELSE
        sh_in(:,:) = qx(:,:,ivap)
-       ch_in(:,:) = qx(:,:,iliq)
+       IF (nqo .EQ. 3) THEN 
+          ch_in(:,:) = qx(:,:,iliq) + qx(:,:,isol)
+       ELSE 
+          ch_in(:,:) = qx(:,:,iliq)
+       ENDIF
     ENDIF
 
-    IF (iflag_phytrac == 1 ) THEN
-
 #ifdef CPP_Dust
-      CALL       phytracr_spl ( debut,lafin , jD_cur,jH_cur,iflag_con,       &  ! I
+    !  Avec SPLA, iflag_phytrac est forcé =1 
+    CALL       phytracr_spl ( debut,lafin , jD_cur,jH_cur,iflag_con,       &  ! I
                       pdtphys,ftsol,                                   &  ! I
                       t,q_seri,paprs,pplay,RHcl,                  &  ! I
@@ -6216,6 +6243,6 @@
 
 #else
-
-    CALL phytrac ( &
+    IF (iflag_phytrac == 1 ) THEN
+      CALL phytrac ( &
          itap,     days_elapsed+1,    jH_cur,   debut, &
          lafin,    phys_tstep,     u, v,     t, &
@@ -6254,7 +6281,8 @@
 
 #endif
-
-#endif
     ENDIF    ! (iflag_phytrac=1)
+
+#endif
+    !ENDIF    ! (iflag_phytrac=1)
 
     IF (offline) THEN
@@ -6345,5 +6373,5 @@
 #endif
     !
-    IF (type_trac == 'inca') THEN
+    IF (type_trac == 'inca' .OR. type_trac == 'inco') THEN
 #ifdef INCA
        CALL VTe(VTphysiq)
@@ -6354,5 +6382,5 @@
             pplay, &
             t_seri, &
-            tr_seri, &
+            tr_seri(:,:,1+nqCO2:nbtr), &
             nbtr, &
             paprs, &
@@ -6620,5 +6648,5 @@
   CALL phys_output_write_spl(itap, pdtphys, paprs, pphis,  &
        pplay, lmax_th, aerosol_couple,                 &
-       ok_ade, ok_aie, ivap, ok_sync,         &
+       ok_ade, ok_aie, ivap, ok_sync,                  &
        ptconv, read_climoz, clevSTD,                   &
        ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse,      &
@@ -6639,5 +6667,5 @@
 #endif
 
-! On remet des variables a .false. apres un premier appel
+! Pour XIOS : On remet des variables a .false. apres un premier appel
     IF (debut) THEN
 #ifdef CPP_XIOS
