Index: LMDZ6/branches/Ocean_skin/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/phys_output_write_mod.F90	(revision 3601)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/phys_output_write_mod.F90	(revision 3605)
@@ -17,5 +17,5 @@
   SUBROUTINE phys_output_write(itap, pdtphys, paprs, pphis, &
        pplay, lmax_th, aerosol_couple,         &
-       ok_ade, ok_aie, ivap, iliq, isol, new_aod, ok_sync, &
+       ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, new_aod, ok_sync, &
        ptconv, read_climoz, clevSTD, ptconvth, &
        d_u, d_t, qx, d_qx, zmasse, flag_aerosol, flag_aerosol_strat, ok_cdnc)
@@ -32,5 +32,6 @@
     USE phys_output_ctrlout_mod, ONLY: o_phis, o_aire, is_ter, is_lic, is_oce, &
          o_longitude, o_latitude, &
-         o_Ahyb, o_Bhyb,o_Ahyb_inter, o_Bhyb_inter, & 
+         o_Ahyb, o_Bhyb,o_Ahyb_bounds, o_Bhyb_bounds, & 
+         o_Ahyb_mid, o_Bhyb_mid,o_Ahyb_mid_bounds, o_Bhyb_mid_bounds, & 
          is_ave, is_sic, o_contfracATM, o_contfracOR, &
          o_aireTER, o_flat, o_slp, o_ptstar, o_pt0, o_tsol, &
@@ -72,5 +73,6 @@
          o_uwat, o_vwat, &
          o_ptop, o_fbase, o_plcl, o_plfc, &
-         o_wbeff, o_convoccur, o_cape_max, o_upwd, o_ep,o_epmax_diag, o_Ma, &
+         o_wbeff, o_convoccur, o_cape_max, o_upwd, o_ep,o_epmax_diag, &
+         o_Mipsh, o_Ma, &
          o_dnwd, o_dnwd0, o_ftime_deepcv, o_ftime_con, o_mc, &
          o_prw, o_prlw, o_prsw, o_s_pblh, o_s_pblt, o_s_lcl, &
@@ -87,5 +89,6 @@
          o_wake_s, o_wake_deltat, o_wake_deltaq, &
          o_wake_omg, o_dtwak, o_dqwak, o_dqwak2d, o_Vprecip, &
-         o_ftd, o_fqd, o_wdtrainA, o_wdtrainM, &
+         o_qtaa, o_Clwaa, &
+         o_ftd, o_fqd, o_wdtrainA, o_wdtrainS, o_wdtrainM, &
          o_n2, o_s2, o_proba_notrig, &
          o_random_notrig, o_ale_bl_stat, &
@@ -196,5 +199,7 @@
 ! Tropopause
          o_p_tropopause, o_z_tropopause, o_t_tropopause,  &
-         o_col_O3_strato, o_col_O3_tropo, &               ! Added ThL
+         o_col_O3_strato, o_col_O3_tropo,                 & 
+!--interactive CO2
+         o_flx_co2_ocean, o_flx_co2_land, o_flx_co2_ff, o_flx_co2_bb, &
          o_t_int, o_s_int, o_ds_ns, o_dt_ns, o_dter, o_dser, o_tkt, o_tks, &
          o_rf, o_taur
@@ -214,4 +219,7 @@
          o_vsed_aer, o_tau_strat_1020, o_ext_strat_1020, o_f_r_wet
 #endif
+    
+    USE phys_output_ctrlout_mod, ONLY: o_heat_volc, o_cool_volc !NL
+    USE phys_state_var_mod, ONLY: heat_volc, cool_volc !NL
 
     USE phys_state_var_mod, ONLY: pctsrf, rain_fall, snow_fall, &
@@ -227,5 +235,5 @@
          delta_tsurf, &
          wstar, cape, ema_pcb, ema_pct, &
-         ema_cbmf, Ma, fm_therm, ale_bl, alp_bl, ale, &
+         ema_cbmf, Mipsh, Ma, fm_therm, ale_bl, alp_bl, ale, &
          alp, cin, wake_pe, wake_dens, wake_s, wake_deltat, &
          wake_deltaq, ftd, fqd, ale_bl_trig, albsol1, &
@@ -244,4 +252,5 @@
 
     USE phys_local_var_mod, ONLY: zxfluxlat, slp, ptstar, pt0, zxtsol, zt2m, &
+         zt2m_cor,zq2m_cor,zu10m_cor,zv10m_cor, zrh2m_cor, zqsat2m_cor, &
          t2m_min_mon, t2m_max_mon, evap, &
          l_mixmin,l_mix, &
@@ -267,6 +276,6 @@
          kh         ,kh_x       ,kh_w       , &
          cv_gen, wake_h, &
-         wake_omg, d_t_wake, d_q_wake, Vprecip, &
-         wdtrainA, wdtrainM, n2, s2, proba_notrig, &
+         wake_omg, d_t_wake, d_q_wake, Vprecip, qtaa, Clw, &
+         wdtrainA, wdtrainS, wdtrainM, n2, s2, proba_notrig, &
          random_notrig, &
          alp_bl_det, alp_bl_fluct_m, alp_bl_conv, &
@@ -329,4 +338,6 @@
 #endif
 
+    USE carbon_cycle_mod, ONLY: fco2_ff, fco2_bb, fco2_land, fco2_ocean
+
     USE phys_output_var_mod, ONLY: vars_defined, snow_o, zfra_o, bils_diss, &
          bils_ec,bils_ech, bils_tke, bils_kinetic, bils_latent, bils_enthalp, &
@@ -368,5 +379,5 @@
     ! ug Pour les sorties XIOS
     USE xios
-    USE wxios, ONLY: wxios_closedef, missing_val
+    USE wxios, ONLY: wxios_closedef, missing_val, wxios_set_context
 #endif
     USE phys_cal_mod, ONLY : mth_len
@@ -390,5 +401,5 @@
     INTEGER, DIMENSION(klon) :: lmax_th
     LOGICAL :: aerosol_couple, ok_sync
-    LOGICAL :: ok_ade, ok_aie, new_aod
+    LOGICAL :: ok_ade, ok_aie, ok_volcan, new_aod
     LOGICAL, DIMENSION(klon, klev) :: ptconv, ptconvth
     REAL :: pdtphys
@@ -417,5 +428,6 @@
     REAL, PARAMETER :: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
 !   REAL, PARAMETER :: missing_val=nf90_fill_real
-    REAL, DIMENSION(klev,2) :: Ahyb_bounds, Bhyb_bounds
+    REAL, DIMENSION(klev+1,2) :: Ahyb_bounds, Bhyb_bounds
+    REAL, DIMENSION(klev,2) :: Ahyb_mid_bounds, Bhyb_mid_bounds
     INTEGER :: ilev
 #ifndef CPP_XIOS
@@ -440,15 +452,34 @@
     CALL set_itau_iophy(itau_w)
 
-    IF (.NOT.vars_defined) THEN
-       iinitend = 2
-    ELSE
+ !   IF (.NOT.vars_defined) THEN
        iinitend = 1
-    ENDIF
-
-    DO ilev=1,klev
-      Ahyb_bounds(ilev,1) = ap(ilev)
-      Ahyb_bounds(ilev,2) = ap(ilev+1)
-      Bhyb_bounds(ilev,1) = bp(ilev)
-      Bhyb_bounds(ilev,2) = bp(ilev+1)
+ !   ELSE
+ !      iinitend = 1
+ !   ENDIF
+
+#ifdef CPP_XIOS
+    CALL wxios_set_context
+#endif
+
+    Ahyb_bounds(1,1) = 0.
+    Ahyb_bounds(1,2) = aps(1)
+    Bhyb_bounds(1,1) = 1.
+    Bhyb_bounds(1,2) = bps(1)    
+    DO ilev=2,klev
+      Ahyb_bounds(ilev,1) = aps(ilev-1)
+      Ahyb_bounds(ilev,2) = aps(ilev)
+      Bhyb_bounds(ilev,1) = bps(ilev-1)
+      Bhyb_bounds(ilev,2) = bps(ilev)
+    ENDDO
+     Ahyb_bounds(klev+1,1) = aps(klev)
+     Ahyb_bounds(klev+1,2) = 0.
+     Bhyb_bounds(klev+1,1) = bps(klev)
+     Bhyb_bounds(klev+1,2) = 0.
+
+    DO ilev=1, klev
+      Ahyb_mid_bounds(ilev,1) = ap(ilev)
+      Ahyb_mid_bounds(ilev,2) = ap(ilev+1)
+      Bhyb_mid_bounds(ilev,1) = bp(ilev)
+      Bhyb_mid_bounds(ilev,2) = bp(ilev+1)
     END DO
 
@@ -565,8 +596,12 @@
        CALL histwrite_phy("R_incl",R_incl)
        CALL histwrite_phy("solaire",solaire)
-       CALL histwrite_phy(o_Ahyb, aps)
-       CALL histwrite_phy(o_Bhyb, bps)
-       CALL histwrite_phy(o_Ahyb_inter, Ahyb_bounds)
-       CALL histwrite_phy(o_Bhyb_inter, Bhyb_bounds)
+       CALL histwrite_phy(o_Ahyb, ap)
+       CALL histwrite_phy(o_Bhyb, bp)
+       CALL histwrite_phy(o_Ahyb_bounds, Ahyb_bounds)
+       CALL histwrite_phy(o_Bhyb_bounds, Bhyb_bounds)
+       CALL histwrite_phy(o_Ahyb_mid, aps)
+       CALL histwrite_phy(o_Bhyb_mid, bps)
+       CALL histwrite_phy(o_Ahyb_mid_bounds, Ahyb_mid_bounds)
+       CALL histwrite_phy(o_Bhyb_mid_bounds, Bhyb_mid_bounds)
        CALL histwrite_phy(o_longitude, longitude_deg)
        CALL histwrite_phy(o_latitude, latitude_deg)
@@ -636,7 +671,7 @@
        CALL histwrite_phy(o_slp, slp)
        CALL histwrite_phy(o_tsol, zxtsol)
-       CALL histwrite_phy(o_t2m, zt2m)
-       CALL histwrite_phy(o_t2m_min, zt2m)
-       CALL histwrite_phy(o_t2m_max, zt2m)
+       CALL histwrite_phy(o_t2m, zt2m_cor)
+       CALL histwrite_phy(o_t2m_min, zt2m_cor)
+       CALL histwrite_phy(o_t2m_max, zt2m_cor)
        CALL histwrite_phy(o_t2m_max_mon, t2m_max_mon)
        CALL histwrite_phy(o_t2m_min_mon, t2m_min_mon)
@@ -644,5 +679,5 @@
        IF (vars_defined) THEN
           DO i=1, klon
-             zx_tmp_fi2d(i)=SQRT(zu10m(i)*zu10m(i)+zv10m(i)*zv10m(i))
+             zx_tmp_fi2d(i)=SQRT(zu10m_cor(i)*zu10m_cor(i)+zv10m_cor(i)*zv10m_cor(i))
           ENDDO
        ENDIF
@@ -651,5 +686,5 @@
        IF (vars_defined) THEN
           DO i=1, klon
-             zx_tmp_fi2d(i)=SQRT(zu10m(i)*zu10m(i)+zv10m(i)*zv10m(i))
+             zx_tmp_fi2d(i)=SQRT(zu10m_cor(i)*zu10m_cor(i)+zv10m_cor(i)*zv10m_cor(i))
           ENDDO
        ENDIF
@@ -664,8 +699,8 @@
        ENDIF
        CALL histwrite_phy(o_sicf, zx_tmp_fi2d)
-       CALL histwrite_phy(o_q2m, zq2m)
+       CALL histwrite_phy(o_q2m, zq2m_cor)
        CALL histwrite_phy(o_ustar, zustar)
-       CALL histwrite_phy(o_u10m, zu10m)
-       CALL histwrite_phy(o_v10m, zv10m)
+       CALL histwrite_phy(o_u10m, zu10m_cor)
+       CALL histwrite_phy(o_v10m, zv10m_cor)
 
        IF (vars_defined) THEN
@@ -716,6 +751,15 @@
        CALL histwrite_phy(o_fsnow, zfra_o)
        CALL histwrite_phy(o_evap, evap)
-       CALL histwrite_phy(o_tops, topsw*swradcorr)
-       CALL histwrite_phy(o_tops0, topsw0*swradcorr)
+
+       IF (vars_defined) THEN
+         zx_tmp_fi2d = topsw*swradcorr
+       ENDIF
+       CALL histwrite_phy(o_tops, zx_tmp_fi2d)
+
+       IF (vars_defined) THEN
+         zx_tmp_fi2d = topsw0*swradcorr
+       ENDIF
+       CALL histwrite_phy(o_tops0, zx_tmp_fi2d)
+
        CALL histwrite_phy(o_topl, toplw)
        CALL histwrite_phy(o_topl0, toplw0)
@@ -750,14 +794,41 @@
        ENDIF
        CALL histwrite_phy(o_nettop, zx_tmp_fi2d)
-       CALL histwrite_phy(o_SWup200, SWup200*swradcorr)
-       CALL histwrite_phy(o_SWup200clr, SWup200clr*swradcorr)
-       CALL histwrite_phy(o_SWdn200, SWdn200*swradcorr)
-       CALL histwrite_phy(o_SWdn200clr, SWdn200clr*swradcorr)
+       
+       IF (vars_defined) THEN
+          zx_tmp_fi2d = SWup200*swradcorr
+       ENDIF
+       CALL histwrite_phy(o_SWup200, zx_tmp_fi2d)
+       
+       IF (vars_defined) THEN
+          zx_tmp_fi2d = SWup200clr*swradcorr
+       ENDIF
+       CALL histwrite_phy(o_SWup200clr, zx_tmp_fi2d)
+       
+       IF (vars_defined) THEN
+          zx_tmp_fi2d = SWdn200*swradcorr
+       ENDIF
+       CALL histwrite_phy(o_SWdn200, zx_tmp_fi2d)
+       
+       
+       IF (vars_defined) THEN
+          zx_tmp_fi2d = SWdn200clr*swradcorr
+       ENDIF
+       CALL histwrite_phy(o_SWdn200clr, zx_tmp_fi2d)
+       
        CALL histwrite_phy(o_LWup200, LWup200)
        CALL histwrite_phy(o_LWup200clr, LWup200clr)
        CALL histwrite_phy(o_LWdn200, LWdn200)
        CALL histwrite_phy(o_LWdn200clr, LWdn200clr)
-       CALL histwrite_phy(o_sols, solsw*swradcorr)
-       CALL histwrite_phy(o_sols0, solsw0*swradcorr)
+       
+       IF (vars_defined) THEN
+          zx_tmp_fi2d = solsw*swradcorr
+       ENDIF
+       CALL histwrite_phy(o_sols, zx_tmp_fi2d)
+       
+       
+       IF (vars_defined) THEN
+          zx_tmp_fi2d = solsw0*swradcorr
+       ENDIF
+       CALL histwrite_phy(o_sols0, zx_tmp_fi2d)
        CALL histwrite_phy(o_soll, sollw)
        CALL histwrite_phy(o_soll0, sollw0)
@@ -909,5 +980,5 @@
           IF (iflag_pbl > 1) THEN
              CALL histwrite_phy(o_tke_srf(nsrf),  pbl_tke(:,1:klev,nsrf))
-             CALL histwrite_phy(o_l_mix(nsrf),  l_mix(:,1:klev,nsrf))
+             !CALL histwrite_phy(o_l_mix(nsrf),  l_mix(:,1:klev,nsrf))
              CALL histwrite_phy(o_l_mixmin(nsrf),  l_mixmin(:,1:klev,nsrf))
              CALL histwrite_phy(o_tke_max_srf(nsrf),  pbl_tke(:,1:klev,nsrf))
@@ -954,8 +1025,17 @@
        CALL histwrite_phy(o_cldt, cldt)
        CALL histwrite_phy(o_JrNt, JrNt)
-       CALL histwrite_phy(o_cldljn, cldl*JrNt)
-       CALL histwrite_phy(o_cldmjn, cldm*JrNt)
-       CALL histwrite_phy(o_cldhjn, cldh*JrNt)
-       CALL histwrite_phy(o_cldtjn, cldt*JrNt)
+       
+       IF (vars_defined)  zx_tmp_fi2d=cldl*JrNt      
+       CALL histwrite_phy(o_cldljn, zx_tmp_fi2d)
+       
+       IF (vars_defined)  zx_tmp_fi2d=cldm*JrNt      
+       CALL histwrite_phy(o_cldmjn, zx_tmp_fi2d)
+       
+       IF (vars_defined)  zx_tmp_fi2d=cldh*JrNt
+       CALL histwrite_phy(o_cldhjn, zx_tmp_fi2d)
+       
+       IF (vars_defined)  zx_tmp_fi2d=cldt*JrNt
+       CALL histwrite_phy(o_cldtjn, zx_tmp_fi2d)
+       
        CALL histwrite_phy(o_cldq, cldq)
        IF (vars_defined)       zx_tmp_fi2d(1:klon) = flwp(1:klon)
@@ -1108,4 +1188,5 @@
        ! Wakes
        IF (iflag_con.EQ.3) THEN
+          CALL histwrite_phy(o_Mipsh, Mipsh)
           IF (iflag_wake>=1) THEN
              CALL histwrite_phy(o_ale_wk, ale_wake)
@@ -1158,5 +1239,5 @@
              IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_q_wake(1:klon,1:klev)/pdtphys
              CALL histwrite_phy(o_dqwak, zx_tmp_fi3d)
-             CALL water_int(klon,klev,zx_tmp_fi3d,zmasse,zx_tmp_fi2d)
+             IF (vars_defined) CALL water_int(klon,klev,zx_tmp_fi3d,zmasse,zx_tmp_fi2d)
              CALL histwrite_phy(o_dqwak2d, zx_tmp_fi2d)
           ENDIF ! iflag_wake>=1
@@ -1168,5 +1249,8 @@
           ! etendue a iflag_con=3 (jyg)
           CALL histwrite_phy(o_Vprecip, Vprecip)
+          CALL histwrite_phy(o_qtaa, qtaa)
+          CALL histwrite_phy(o_clwaa, clw)
           CALL histwrite_phy(o_wdtrainA, wdtrainA)
+          CALL histwrite_phy(o_wdtrainS, wdtrainS)
           CALL histwrite_phy(o_wdtrainM, wdtrainM)
        ENDIF !(iflag_con.EQ.3.or.iflag_con.EQ.30)
@@ -1208,5 +1292,5 @@
           IF (slab_gm) THEN
              CALL histwrite_phy(o_slab_gm, dt_gm(:,1:nslay))
-          END IF
+          ENDIF
           IF (slab_hdiff) THEN
             IF (nslay.EQ.1) THEN
@@ -1257,5 +1341,5 @@
 !       CALL histwrite_phy(o_rh2m_max, zx_tmp_fi2d)
 
-       CALL histwrite_phy(o_qsat2m, qsat2m)
+       CALL histwrite_phy(o_qsat2m, zqsat2m_cor)
        CALL histwrite_phy(o_tpot, tpot)
        CALL histwrite_phy(o_tpote, tpote)
@@ -1322,13 +1406,18 @@
 
 ! ThL -- In the following, we assume read_climoz == 1
-       zx_tmp_fi2d = 0.0    ! Computation for strato, added ThL
-       DO k=1, klev
-          zx_tmp_fi2d(:) = zx_tmp_fi2d(:) + wo(:,k,1) * stratomask(:,k) * 1.e3
-       END DO
+       IF (vars_defined) THEN
+         zx_tmp_fi2d = 0.0    ! Computation for strato, added ThL
+         DO k=1, klev
+            zx_tmp_fi2d(:) = zx_tmp_fi2d(:) + wo(:,k,1) * stratomask(:,k) * 1.e3
+         END DO
+       ENDIF
        CALL histwrite_phy(o_col_O3_strato, zx_tmp_fi2d) ! Added ThL
-       zx_tmp_fi2d = 0.0    ! Computation for tropo, added ThL
-       DO k=1, klev
-          zx_tmp_fi2d(:) = zx_tmp_fi2d(:) + wo(:,k,1) * (1.0-stratomask(:,k)) * 1.e3
-       END DO
+
+       IF (vars_defined) THEN
+         zx_tmp_fi2d = 0.0    ! Computation for tropo, added ThL
+         DO k=1, klev
+            zx_tmp_fi2d(:) = zx_tmp_fi2d(:) + wo(:,k,1) * (1.0-stratomask(:,k)) * 1.e3
+         END DO
+       ENDIF
        CALL histwrite_phy(o_col_O3_tropo, zx_tmp_fi2d)   ! Added ThL
 ! end add ThL
@@ -1370,9 +1459,28 @@
        ENDIF
 #endif
+       !NL
+       IF (ok_volcan .AND. ok_ade) THEN
+          DO k=1, klev
+             IF (vars_defined) zx_tmp_fi3d(:,k)=heat_volc(:,k)*swradcorr(:)
+          ENDDO
+          CALL histwrite_phy(o_heat_volc, zx_tmp_fi3d)
+          DO k=1, klev
+             IF (vars_defined) zx_tmp_fi3d(:,k)=cool_volc(:,k)
+          ENDDO
+          CALL histwrite_phy(o_cool_volc, zx_tmp_fi3d)
+       ENDIF
        IF (ok_ade) THEN
-          CALL histwrite_phy(o_topswad, topswad_aero*swradcorr)
-          CALL histwrite_phy(o_topswad0, topswad0_aero*swradcorr)
-          CALL histwrite_phy(o_solswad, solswad_aero*swradcorr)
-          CALL histwrite_phy(o_solswad0, solswad0_aero*swradcorr)
+          IF (vars_defined) zx_tmp_fi2d(:)=topswad_aero*swradcorr
+          CALL histwrite_phy(o_topswad, zx_tmp_fi2d)
+          
+          IF (vars_defined) zx_tmp_fi2d(:)=topswad0_aero*swradcorr
+          CALL histwrite_phy(o_topswad0, zx_tmp_fi2d)
+                    
+          IF (vars_defined) zx_tmp_fi2d(:)=solswad_aero*swradcorr
+          CALL histwrite_phy(o_solswad, zx_tmp_fi2d)
+                    
+          IF (vars_defined) zx_tmp_fi2d(:)=solswad0_aero*swradcorr
+          CALL histwrite_phy(o_solswad0, zx_tmp_fi2d)
+          
           CALL histwrite_phy(o_toplwad, toplwad_aero)
           CALL histwrite_phy(o_toplwad0, toplwad0_aero)
@@ -1381,34 +1489,37 @@
           !====MS forcing diagnostics
           IF (new_aod) THEN
-             zx_tmp_fi2d(:)=topsw_aero(:,1)*swradcorr(:)
+          !ym warning : topsw_aero, solsw_aero, topsw0_aero, solsw0_aero are not defined by model
+          !ym => init to 0 in radlwsw_m.F90 ztopsw_aero, zsolsw_aero, ztopsw0_aero, zsolsw0_aero
+
+             IF (vars_defined) zx_tmp_fi2d(:)=topsw_aero(:,1)*swradcorr(:)
              CALL histwrite_phy(o_swtoaas_nat,zx_tmp_fi2d)
-             zx_tmp_fi2d(:)=solsw_aero(:,1)*swradcorr(:)
+             IF (vars_defined) zx_tmp_fi2d(:)=solsw_aero(:,1)*swradcorr(:)
              CALL histwrite_phy(o_swsrfas_nat,zx_tmp_fi2d)
-             zx_tmp_fi2d(:)=topsw0_aero(:,1)*swradcorr(:)
+             IF (vars_defined) zx_tmp_fi2d(:)=topsw0_aero(:,1)*swradcorr(:)
              CALL histwrite_phy(o_swtoacs_nat,zx_tmp_fi2d)
-             zx_tmp_fi2d(:)=solsw0_aero(:,1)*swradcorr(:)
+             IF (vars_defined) zx_tmp_fi2d(:)=solsw0_aero(:,1)*swradcorr(:)
              CALL histwrite_phy(o_swsrfcs_nat,zx_tmp_fi2d)
              !ant
-             zx_tmp_fi2d(:)=topsw_aero(:,2)*swradcorr(:)
+             IF (vars_defined) zx_tmp_fi2d(:)=topsw_aero(:,2)*swradcorr(:)
              CALL histwrite_phy(o_swtoaas_ant,zx_tmp_fi2d)
-             zx_tmp_fi2d(:)=solsw_aero(:,2)*swradcorr(:)
+             IF (vars_defined) zx_tmp_fi2d(:)=solsw_aero(:,2)*swradcorr(:)
              CALL histwrite_phy(o_swsrfas_ant,zx_tmp_fi2d)
-             zx_tmp_fi2d(:)=topsw0_aero(:,2)*swradcorr(:)
+             IF (vars_defined) zx_tmp_fi2d(:)=topsw0_aero(:,2)*swradcorr(:)
              CALL histwrite_phy(o_swtoacs_ant,zx_tmp_fi2d)
-             zx_tmp_fi2d(:)=solsw0_aero(:,2)*swradcorr(:)
+             IF (vars_defined) zx_tmp_fi2d(:)=solsw0_aero(:,2)*swradcorr(:)
              CALL histwrite_phy(o_swsrfcs_ant,zx_tmp_fi2d)
              !cf
              IF (.not. aerosol_couple) THEN
-                zx_tmp_fi2d(:)=topswcf_aero(:,1)*swradcorr(:)
+                IF (vars_defined) zx_tmp_fi2d(:)=topswcf_aero(:,1)*swradcorr(:)
                 CALL histwrite_phy(o_swtoacf_nat,zx_tmp_fi2d)
-                zx_tmp_fi2d(:)=solswcf_aero(:,1)*swradcorr(:)
+                IF (vars_defined) zx_tmp_fi2d(:)=solswcf_aero(:,1)*swradcorr(:)
                 CALL histwrite_phy(o_swsrfcf_nat,zx_tmp_fi2d)
-                zx_tmp_fi2d(:)=topswcf_aero(:,2)*swradcorr(:)
+                IF (vars_defined) zx_tmp_fi2d(:)=topswcf_aero(:,2)*swradcorr(:)
                 CALL histwrite_phy(o_swtoacf_ant,zx_tmp_fi2d)
-                zx_tmp_fi2d(:)=solswcf_aero(:,2)*swradcorr(:)
+                IF (vars_defined) zx_tmp_fi2d(:)=solswcf_aero(:,2)*swradcorr(:)
                 CALL histwrite_phy(o_swsrfcf_ant,zx_tmp_fi2d)
-                zx_tmp_fi2d(:)=topswcf_aero(:,3)*swradcorr(:)
+                IF (vars_defined) zx_tmp_fi2d(:)=topswcf_aero(:,3)*swradcorr(:)
                 CALL histwrite_phy(o_swtoacf_zero,zx_tmp_fi2d)
-                zx_tmp_fi2d(:)=solswcf_aero(:,3)*swradcorr(:)
+                IF (vars_defined) zx_tmp_fi2d(:)=solswcf_aero(:,3)*swradcorr(:)
                 CALL histwrite_phy(o_swsrfcf_zero,zx_tmp_fi2d)
              ENDIF
@@ -1417,8 +1528,15 @@
        ENDIF
        IF (ok_aie) THEN
-          CALL histwrite_phy(o_topswai, topswai_aero*swradcorr)
-          CALL histwrite_phy(o_toplwai, toplwai_aero*swradcorr)
-          CALL histwrite_phy(o_solswai, solswai_aero*swradcorr)
-          CALL histwrite_phy(o_sollwai, sollwai_aero*swradcorr)
+          IF (vars_defined) zx_tmp_fi2d(:)= topswai_aero*swradcorr
+          CALL histwrite_phy(o_topswai, zx_tmp_fi2d)
+          
+          IF (vars_defined) zx_tmp_fi2d(:)=toplwai_aero*swradcorr
+          CALL histwrite_phy(o_toplwai, zx_tmp_fi2d)
+          
+          IF (vars_defined) zx_tmp_fi2d(:)=solswai_aero*swradcorr
+          CALL histwrite_phy(o_solswai, zx_tmp_fi2d)
+          
+          IF (vars_defined) zx_tmp_fi2d(:)=sollwai_aero*swradcorr
+          CALL histwrite_phy(o_sollwai, zx_tmp_fi2d)
        ENDIF
        IF (flag_aerosol.GT.0.AND.ok_cdnc) THEN
@@ -1435,5 +1553,5 @@
           CALL histwrite_phy(o_icc3dstra, icc3dstra)
           CALL histwrite_phy(o_cldicemxrat, zfice)
-          zx_tmp_fi3d(:,:)=1-zfice(:,:)
+          IF (vars_defined) zx_tmp_fi3d(:,:)=1-zfice(:,:)
           CALL histwrite_phy(o_cldwatmxrat, zx_tmp_fi3d)
           CALL histwrite_phy(o_reffclwtop, reffclwtop)
@@ -1450,5 +1568,8 @@
        CALL histwrite_phy(o_ovap, q_seri)
        CALL histwrite_phy(o_oliq, ql_seri)
-       CALL histwrite_phy(o_ocond, ql_seri+qs_seri)
+
+       IF (vars_defined) zx_tmp_fi3d = ql_seri+qs_seri
+       CALL histwrite_phy(o_ocond, zx_tmp_fi3d)
+
        CALL histwrite_phy(o_geop, zphi)
        CALL histwrite_phy(o_vitu, u_seri)
@@ -1457,5 +1578,7 @@
        CALL histwrite_phy(o_pres, pplay)
        CALL histwrite_phy(o_paprs, paprs(:,1:klev))
-       CALL histwrite_phy(o_zfull,zphi/RG)
+       
+       IF (vars_defined) zx_tmp_fi3d = zphi/RG
+       CALL histwrite_phy(o_zfull,zx_tmp_fi3d)
 
 #ifdef CPP_XIOS
@@ -1502,10 +1625,11 @@
        CALL histwrite_phy(o_rnebjn, zx_tmp_fi3d)
        CALL histwrite_phy(o_rhum, zx_rh)
-       CALL histwrite_phy(o_ozone, &
-            wo(:, :, 1) * dobson_u * 1e3 / zmasse / rmo3 * rmd)
+       
+       IF (vars_defined) zx_tmp_fi3d = wo(:, :, 1) * dobson_u * 1e3 / zmasse / rmo3 * rmd
+       CALL histwrite_phy(o_ozone, zx_tmp_fi3d)
 
        IF (read_climoz == 2) THEN
-          CALL histwrite_phy(o_ozone_light, &
-               wo(:, :, 2) * dobson_u * 1e3 / zmasse / rmo3 * rmd)
+         IF (vars_defined) zx_tmp_fi3d = wo(:, :, 2) * dobson_u * 1e3 / zmasse / rmo3 * rmd 
+         CALL histwrite_phy(o_ozone_light, zx_tmp_fi3d)
        ENDIF
 
@@ -1515,14 +1639,14 @@
 
        CALL histwrite_phy(o_dqphy,  d_qx(:,:,ivap))
-       CALL water_int(klon,klev,d_qx(:,:,ivap),zmasse,zx_tmp_fi2d)
+       IF (vars_defined) CALL water_int(klon,klev,d_qx(:,:,ivap),zmasse,zx_tmp_fi2d)
        CALL histwrite_phy(o_dqphy2d,  zx_tmp_fi2d)
 
        CALL histwrite_phy(o_dqlphy,  d_qx(:,:,iliq))
-       CALL water_int(klon,klev,d_qx(:,:,iliq),zmasse,zx_tmp_fi2d)
+       IF (vars_defined) CALL water_int(klon,klev,d_qx(:,:,iliq),zmasse,zx_tmp_fi2d)
        CALL histwrite_phy(o_dqlphy2d,  zx_tmp_fi2d)
 
        IF (nqo.EQ.3) THEN
        CALL histwrite_phy(o_dqsphy,  d_qx(:,:,isol))
-       CALL water_int(klon,klev,d_qx(:,:,isol),zmasse,zx_tmp_fi2d)
+       IF (vars_defined) CALL water_int(klon,klev,d_qx(:,:,isol),zmasse,zx_tmp_fi2d)
        CALL histwrite_phy(o_dqsphy2d,  zx_tmp_fi2d)
        ELSE 
@@ -1589,5 +1713,5 @@
        ENDIF
        CALL histwrite_phy(o_dtcon, zx_tmp_fi3d)
-       if(iflag_thermals.eq.0)then
+       IF (iflag_thermals.EQ.0) THEN
           IF (vars_defined) THEN
              zx_tmp_fi3d(1:klon,1:klev)=d_t_con(1:klon,1:klev)/pdtphys + &
@@ -1595,5 +1719,5 @@
           ENDIF
           CALL histwrite_phy(o_tntc, zx_tmp_fi3d)
-       else if(iflag_thermals.ge.1.and.iflag_wake.EQ.1)then
+       ELSE IF(iflag_thermals.GE.1.AND.iflag_wake.EQ.1) THEN
           IF (vars_defined) THEN
              zx_tmp_fi3d(1:klon,1:klev)=d_t_con(1:klon,1:klev)/pdtphys + &
@@ -1602,5 +1726,5 @@
           ENDIF
           CALL histwrite_phy(o_tntc, zx_tmp_fi3d)
-       endif
+       ENDIF
        IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_u_con(1:klon,1:klev)/pdtphys
        CALL histwrite_phy(o_ducon, zx_tmp_fi3d)
@@ -1609,5 +1733,5 @@
        IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_q_con(1:klon,1:klev)/pdtphys
        CALL histwrite_phy(o_dqcon, zx_tmp_fi3d)
-       CALL water_int(klon,klev,zx_tmp_fi3d,zmasse,zx_tmp_fi2d)
+       IF (vars_defined) CALL water_int(klon,klev,zx_tmp_fi3d,zmasse,zx_tmp_fi2d)
        CALL histwrite_phy(o_dqcon2d, zx_tmp_fi2d)
 
@@ -1631,5 +1755,5 @@
        IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_q_lsc(1:klon,1:klev)/pdtphys
        CALL histwrite_phy(o_dqlsc, zx_tmp_fi3d)
-       CALL water_int(klon,klev,zx_tmp_fi3d,zmasse,zx_tmp_fi2d)
+       IF (vars_defined) CALL water_int(klon,klev,zx_tmp_fi3d,zmasse,zx_tmp_fi2d)
        CALL histwrite_phy(o_dqlsc2d, zx_tmp_fi2d)
        IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=beta_prec(1:klon,1:klev)
@@ -1644,16 +1768,16 @@
           IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_q_lscth(1:klon,1:klev)/pdtphys
           CALL histwrite_phy(o_dqlscth, zx_tmp_fi3d)
-          CALL water_int(klon,klev,zx_tmp_fi3d,zmasse,zx_tmp_fi2d)
+          IF (vars_defined) CALL water_int(klon,klev,zx_tmp_fi3d,zmasse,zx_tmp_fi2d)
           CALL histwrite_phy(o_dqlscth2d, zx_tmp_fi2d)
           IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_q_lscst(1:klon,1:klev)/pdtphys
           CALL histwrite_phy(o_dqlscst, zx_tmp_fi3d)
-          CALL water_int(klon,klev,zx_tmp_fi3d,zmasse,zx_tmp_fi2d)
+          IF (vars_defined) CALL water_int(klon,klev,zx_tmp_fi3d,zmasse,zx_tmp_fi2d)
           CALL histwrite_phy(o_dqlscst2d, zx_tmp_fi2d)
           CALL histwrite_phy(o_plulth, plul_th)
           CALL histwrite_phy(o_plulst, plul_st)
           IF (vars_defined) THEN
-             do i=1,klon
+             DO i=1,klon
                 zx_tmp_fi2d(1:klon)=lmax_th(:)
-             enddo
+             ENDDO
           ENDIF
           CALL histwrite_phy(o_lmaxth, zx_tmp_fi2d)
@@ -1702,5 +1826,5 @@
        IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_q_vdf(1:klon,1:klev)/pdtphys
        CALL histwrite_phy(o_dqvdf, zx_tmp_fi3d)
-       CALL water_int(klon,klev,zx_tmp_fi3d,zmasse,zx_tmp_fi2d)
+       IF (vars_defined) CALL water_int(klon,klev,zx_tmp_fi3d,zmasse,zx_tmp_fi2d)
        CALL histwrite_phy(o_dqvdf2d, zx_tmp_fi2d)
        IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_eva(1:klon,1:klev)/pdtphys
@@ -1708,5 +1832,5 @@
        IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_q_eva(1:klon,1:klev)/pdtphys
        CALL histwrite_phy(o_dqeva, zx_tmp_fi3d)
-       CALL water_int(klon,klev,zx_tmp_fi3d,zmasse,zx_tmp_fi2d)
+       IF (vars_defined) CALL water_int(klon,klev,zx_tmp_fi3d,zmasse,zx_tmp_fi2d)
        CALL histwrite_phy(o_dqeva2d, zx_tmp_fi2d)
        CALL histwrite_phy(o_ratqs, ratqs)
@@ -1747,5 +1871,5 @@
           ENDIF
           CALL histwrite_phy(o_dqthe, zx_tmp_fi3d)
-          CALL water_int(klon,klev,zx_tmp_fi3d,zmasse,zx_tmp_fi2d)
+          IF (vars_defined) CALL water_int(klon,klev,zx_tmp_fi3d,zmasse,zx_tmp_fi2d)
           CALL histwrite_phy(o_dqthe2d, zx_tmp_fi2d)
        ENDIF !iflag_thermals
@@ -1754,5 +1878,5 @@
        IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_q_ajsb(1:klon,1:klev)/pdtphys
        CALL histwrite_phy(o_dqajs, zx_tmp_fi3d)
-       CALL water_int(klon,klev,zx_tmp_fi3d,zmasse,zx_tmp_fi2d)
+       IF (vars_defined) CALL water_int(klon,klev,zx_tmp_fi3d,zmasse,zx_tmp_fi2d)
        CALL histwrite_phy(o_dqajs2d, zx_tmp_fi2d)
        IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_swr(1:klon,1:klev)/pdtphys
@@ -1790,6 +1914,10 @@
 
        IF (ok_hines) THEN
-          CALL histwrite_phy(o_du_gwd_hines, du_gwd_hines/pdtphys)
-          CALL histwrite_phy(o_dv_gwd_hines, dv_gwd_hines/pdtphys)
+          IF (vars_defined) zx_tmp_fi3d=du_gwd_hines/pdtphys
+          CALL histwrite_phy(o_du_gwd_hines, zx_tmp_fi3d)
+
+          IF (vars_defined) zx_tmp_fi3d= dv_gwd_hines/pdtphys          
+          CALL histwrite_phy(o_dv_gwd_hines, zx_tmp_fi3d)
+          
           IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_hin(1:klon,1:klev)/pdtphys
           CALL histwrite_phy(o_dthin, zx_tmp_fi3d)
@@ -1799,6 +1927,10 @@
 
        IF (.not. ok_hines .and. ok_gwd_rando) THEN
-          CALL histwrite_phy(o_du_gwd_front, du_gwd_front / pdtphys)
-          CALL histwrite_phy(o_dv_gwd_front, dv_gwd_front / pdtphys)
+          IF (vars_defined)  zx_tmp_fi3d=du_gwd_front / pdtphys
+          CALL histwrite_phy(o_du_gwd_front, zx_tmp_fi3d)
+          
+          IF (vars_defined)  zx_tmp_fi3d=dv_gwd_front / pdtphys
+          CALL histwrite_phy(o_dv_gwd_front, zx_tmp_fi3d)
+          
           CALL histwrite_phy(o_ustr_gwd_front, zustr_gwd_front)
           CALL histwrite_phy(o_vstr_gwd_front, zvstr_gwd_front)
@@ -1806,6 +1938,9 @@
 
        IF (ok_gwd_rando) THEN
-          CALL histwrite_phy(o_du_gwd_rando, du_gwd_rando / pdtphys)
-          CALL histwrite_phy(o_dv_gwd_rando, dv_gwd_rando / pdtphys)
+          IF (vars_defined)  zx_tmp_fi3d=du_gwd_rando / pdtphys
+          CALL histwrite_phy(o_du_gwd_rando, zx_tmp_fi3d)
+          
+          IF (vars_defined)  zx_tmp_fi3d=dv_gwd_rando / pdtphys
+          CALL histwrite_phy(o_dv_gwd_rando, zx_tmp_fi3d)
           CALL histwrite_phy(o_ustr_gwd_rando, zustr_gwd_rando)
           CALL histwrite_phy(o_vstr_gwd_rando, zvstr_gwd_rando)
@@ -1815,30 +1950,53 @@
 
        IF (ok_qch4) THEN
-          CALL histwrite_phy(o_dqch4, d_q_ch4 / pdtphys)
-       ENDIF
-
-       DO k=1, klevp1
-         zx_tmp_fi3d1(:,k)=swup(:,k)*swradcorr(:)
-       ENDDO
+          IF (vars_defined) zx_tmp_fi3d=d_q_ch4 / pdtphys
+          CALL histwrite_phy(o_dqch4, zx_tmp_fi3d)
+       ENDIF
+       
+       IF (vars_defined) THEN
+         DO k=1, klevp1
+           zx_tmp_fi3d1(:,k)=swup(:,k)*swradcorr(:)
+         ENDDO
+       ENDIF
+       
        CALL histwrite_phy(o_rsu, zx_tmp_fi3d1)
-       DO k=1, klevp1
-         zx_tmp_fi3d1(:,k)=swdn(:,k)*swradcorr(:)
-       ENDDO
+
+
+       IF (vars_defined) THEN
+         DO k=1, klevp1
+           zx_tmp_fi3d1(:,k)=swdn(:,k)*swradcorr(:)
+         ENDDO
+       ENDIF
+       
        CALL histwrite_phy(o_rsd, zx_tmp_fi3d1)
-       DO k=1, klevp1
-         zx_tmp_fi3d1(:,k)=swup0(:,k)*swradcorr(:)
-       ENDDO
+
+       IF (vars_defined) THEN
+         DO k=1, klevp1
+           zx_tmp_fi3d1(:,k)=swup0(:,k)*swradcorr(:)
+         ENDDO
+       ENDIF
+       
        CALL histwrite_phy(o_rsucs, zx_tmp_fi3d1)
-       DO k=1, klevp1
-         zx_tmp_fi3d1(:,k)=swupc0(:,k)*swradcorr(:)
-       ENDDO
+
+       IF (vars_defined) THEN
+         DO k=1, klevp1
+           zx_tmp_fi3d1(:,k)=swupc0(:,k)*swradcorr(:)
+         ENDDO
+       ENDIF
        CALL histwrite_phy(o_rsucsaf, zx_tmp_fi3d1)
-       DO k=1, klevp1
-         zx_tmp_fi3d1(:,k)=swdn0(:,k)*swradcorr(:)
-       ENDDO
+
+       IF (vars_defined) THEN
+         DO k=1, klevp1
+           zx_tmp_fi3d1(:,k)=swdn0(:,k)*swradcorr(:)
+         ENDDO
+       ENDIF
        CALL histwrite_phy(o_rsdcs, zx_tmp_fi3d1)
-       DO k=1, klevp1
-         zx_tmp_fi3d1(:,k)=swdnc0(:,k)*swradcorr(:)
-       ENDDO
+
+
+       IF (vars_defined) THEN
+         DO k=1, klevp1
+           zx_tmp_fi3d1(:,k)=swdnc0(:,k)*swradcorr(:)
+         ENDDO
+       ENDIF
        CALL histwrite_phy(o_rsdcsaf, zx_tmp_fi3d1)
 
@@ -1894,5 +2052,6 @@
        ELSE IF (iflag_con == 2) THEN
           CALL histwrite_phy(o_mcd,  pmfd)
-          CALL histwrite_phy(o_dmc,  pmfu + pmfd)
+          IF (vars_defined) zx_tmp_fi3d = pmfu + pmfd
+          CALL histwrite_phy(o_dmc,  zx_tmp_fi3d)
        ENDIF
        CALL histwrite_phy(o_ref_liq, ref_liq)
@@ -1908,19 +2067,27 @@
           IF (vars_defined) zx_tmp_fi2d(:) = lwup0p(:,klevp1)
           CALL histwrite_phy(o_rlutcs4co2, zx_tmp_fi2d)
-          DO k=1, klevp1
-            zx_tmp_fi3d1(:,k)=swupp(:,k)*swradcorr(:)
-          ENDDO
+          IF (vars_defined) THEN
+            DO k=1, klevp1
+              zx_tmp_fi3d1(:,k)=swupp(:,k)*swradcorr(:)
+            ENDDO
+          ENDIF
           CALL histwrite_phy(o_rsu4co2, zx_tmp_fi3d1)
-          DO k=1, klevp1
-            zx_tmp_fi3d1(:,k)=swup0p(:,k)*swradcorr(:)
-          ENDDO
+          IF (vars_defined) THEN
+            DO k=1, klevp1
+              zx_tmp_fi3d1(:,k)=swup0p(:,k)*swradcorr(:)
+            ENDDO
+          ENDIF
           CALL histwrite_phy(o_rsucs4co2, zx_tmp_fi3d1)
-          DO k=1, klevp1
-            zx_tmp_fi3d1(:,k)=swdnp(:,k)*swradcorr(:)
-          ENDDO
+          IF (vars_defined) THEN
+            DO k=1, klevp1
+              zx_tmp_fi3d1(:,k)=swdnp(:,k)*swradcorr(:)
+            ENDDO
+          ENDIF
           CALL histwrite_phy(o_rsd4co2, zx_tmp_fi3d1)
-          DO k=1, klevp1
-            zx_tmp_fi3d1(:,k)=swdn0p(:,k)*swradcorr(:)
-          ENDDO
+          IF (vars_defined) THEN
+            DO k=1, klevp1
+              zx_tmp_fi3d1(:,k)=swdn0p(:,k)*swradcorr(:)
+            ENDDO
+          ENDIF
           CALL histwrite_phy(o_rsdcs4co2, zx_tmp_fi3d1)
           CALL histwrite_phy(o_rlu4co2, lwupp)
@@ -2071,7 +2238,6 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        IF (iflag_phytrac == 1 ) then
-       IF (nqtot.GE.nqo+1) THEN
-          DO iq=nqo+1, nqtot
-            IF (type_trac == 'lmdz' .OR. type_trac == 'repr' .OR. type_trac == 'coag') THEN
+         IF (type_trac == 'lmdz' .OR. type_trac == 'repr' .OR. type_trac == 'coag') THEN
+           DO iq=nqo+1, nqtot
              !--3D fields
              CALL histwrite_phy(o_trac(iq-nqo), tr_seri(:,:,iq-nqo))
@@ -2089,5 +2255,5 @@
              CALL histwrite_phy(o_dtr_sat(iq-nqo),d_tr_sat(:,:,iq-nqo))
              CALL histwrite_phy(o_dtr_uscav(iq-nqo),d_tr_uscav(:,:,iq-nqo))
-             !--2D fields
+            !--2D fields
              CALL histwrite_phy(o_dtr_dry(iq-nqo), flux_tr_dry(:,iq-nqo))
              zx_tmp_fi2d=0.
@@ -2100,9 +2266,32 @@
              CALL histwrite_phy(o_trac_cum(iq-nqo), zx_tmp_fi2d)
 #endif
-            ENDIF
-          ENDDO
-       ENDIF
-
-       IF (type_trac == 'repr') THEN
+           ENDDO !--iq
+         ENDIF   !--type_trac
+!
+         IF (type_trac == 'co2i') THEN
+           DO iq=nqo+1, nqtot
+             !--3D fields
+             CALL histwrite_phy(o_trac(iq-nqo), tr_seri(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_vdf(iq-nqo),d_tr_cl(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_the(iq-nqo),d_tr_th(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_con(iq-nqo),d_tr_cv(:,:,iq-nqo))
+             !--2D fields
+             !--CO2 burden
+             zx_tmp_fi2d=0.
+             IF (vars_defined) THEN
+                DO k=1,klev
+                   zx_tmp_fi2d(:)=zx_tmp_fi2d(:)+zmasse(:,k)*tr_seri(:,k,iq-nqo)
+                ENDDO
+             ENDIF
+             CALL histwrite_phy(o_trac_cum(iq-nqo), zx_tmp_fi2d)
+           ENDDO !--iq
+           !--CO2 net fluxes
+           CALL histwrite_phy(o_flx_co2_land,  fco2_land)
+           CALL histwrite_phy(o_flx_co2_ocean, fco2_ocean)
+           CALL histwrite_phy(o_flx_co2_ff,    fco2_ff)
+           CALL histwrite_phy(o_flx_co2_bb,    fco2_bb)
+         ENDIF !--type_trac co2i
+
+         IF (type_trac == 'repr') THEN
 #ifdef REPROBUS
            DO iq=1,nbnas
@@ -2110,5 +2299,5 @@
            ENDDO
 #endif
-       ENDIF
+         ENDIF
 
        ENDIF   !(iflag_phytrac==1)
