--- cs_convection_diffusion.c	2019-09-10 20:40:51.000000000 -0300
+++ cs_convection_diffusion_thermal_mod.c	2019-09-17 11:11:55.969320980 -0300
@@ -2644,6 +2644,7 @@
     /* Steady */
     if (idtvar < 0) {
 
+      bft_printf(_(" STEADY \n"));
       for (int g_id = 0; g_id < n_b_groups; g_id++) {
 #       pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
         for (int t_id = 0; t_id < n_b_threads; t_id++) {
@@ -2699,9 +2700,88 @@
         }
       }
 
-    /* Unsteady */
+      /* The scalar is internal_coupled and an implicit contribution
+       * is required */
+      if (icoupl > 0) {
+	bft_printf(_(" STEADY INTERNAL COUPLING \n"));
+	/* Prepare data for sending */
+	BFT_MALLOC(pvar_distant, n_distant, cs_real_t);
+	
+	for (cs_lnum_t ii = 0; ii < n_distant; ii++) {
+	  cs_lnum_t face_id = faces_distant[ii];
+	  cs_lnum_t jj = b_face_cells[face_id];
+	  cs_real_t pip;
+	  
+	  cs_real_t bldfrp = (cs_real_t) ircflp;
+	  /* Local limitation of the reconstruction */
+	  if (df_limiter != NULL && ircflp > 0)
+	    bldfrp = CS_MAX(df_limiter[jj], 0.);
+	  
+	  cs_b_cd_unsteady(bldfrp,
+			   diipb[face_id],
+			   grad[jj],
+			   _pvar[jj],
+			   &pip);
+	  pvar_distant[ii] = pip;
+	}
+	
+	/* Receive data */
+	BFT_MALLOC(pvar_local, n_local, cs_real_t);
+	cs_internal_coupling_exchange_var(cpl,
+					  1, /* Dimension */
+					  pvar_distant,
+					  pvar_local);
+	
+	/* Exchange diffusion limiter */
+	if (df_limiter != NULL) {
+	  BFT_MALLOC(df_limiter_local, n_local, cs_real_t);
+	  cs_internal_coupling_exchange_var(cpl,
+					    1, /* Dimension */
+					    df_limiter,
+					    df_limiter_local);
+	}
+	
+	/* Flux contribution */
+	assert(f != NULL);
+	cs_real_t *hintp = f->bc_coeffs->hint;
+	cs_real_t *hextp = f->bc_coeffs->hext;
+	for (cs_lnum_t ii = 0; ii < n_local; ii++) {
+	  cs_lnum_t face_id = faces_local[ii];
+	  cs_lnum_t jj = b_face_cells[face_id];
+	  cs_real_t pip, pjp;
+	  cs_real_t fluxi = 0.;
+	  
+	  cs_real_t bldfrp = (cs_real_t) ircflp;
+	  /* Local limitation of the reconstruction */
+	  if (df_limiter != NULL && ircflp > 0)
+	    bldfrp = CS_MAX(CS_MIN(df_limiter_local[ii], df_limiter[jj]), 0.);
+	  
+	  cs_b_cd_unsteady(bldfrp,
+			   diipb[face_id],
+			   grad[jj],
+			   _pvar[jj],
+			   &pip);
+	  
+	  pjp = pvar_local[ii];
+	  
+	  hint = hintp[face_id];
+	  hext = hextp[face_id];
+	  heq = hint * hext / (hint + hext);
+	  
+	  cs_b_diff_flux_coupling(idiffp,
+				  pip,
+				  pjp,
+				  heq,
+				  &fluxi);
+	  
+	  rhs[jj] -= thetap * fluxi;
+	  
+	  
+	}
+      }      
+      /* Unsteady */
     } else {
-
+      bft_printf(_(" UNSTEADY  \n"));
       for (int g_id = 0; g_id < n_b_groups; g_id++) {
 #       pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
         for (int t_id = 0; t_id < n_b_threads; t_id++) {
@@ -2757,6 +2837,7 @@
       /* The scalar is internal_coupled and an implicit contribution
        * is required */
       if (icoupl > 0) {
+	 bft_printf(_("UNSTEADY INTERNAL COUPLING \n"));
         /* Prepare data for sending */
         BFT_MALLOC(pvar_distant, n_distant, cs_real_t);
 
@@ -7758,10 +7839,86 @@
         }
       }
     }
-
+    /* The scalar is internal_coupled and an implicit contribution
+     * is required */
+    if (icoupl > 0) {
+      bft_printf(_(" STEADY INTERNAL COUPLING THERMAL \n"));
+      /* Prepare data for sending */
+      BFT_MALLOC(pvar_distant, n_distant, cs_real_t);
+      
+      for (cs_lnum_t ii = 0; ii < n_distant; ii++) {
+	cs_lnum_t face_id = faces_distant[ii];
+	cs_lnum_t jj = b_face_cells[face_id];
+	cs_real_t pip;
+	
+	cs_real_t bldfrp = (cs_real_t) ircflp;
+	/* Local limitation of the reconstruction */
+	if (df_limiter != NULL && ircflp > 0)
+	  bldfrp = CS_MAX(df_limiter[jj], 0.);
+	
+	cs_b_cd_unsteady(bldfrp,
+			 diipb[face_id],
+			 grad[jj],
+			 _pvar[jj],
+			 &pip);
+	pvar_distant[ii] = pip;
+      }
+      
+      /* Receive data */
+      BFT_MALLOC(pvar_local, n_local, cs_real_t);
+      cs_internal_coupling_exchange_var(cpl,
+					1, /* Dimension */
+					pvar_distant,
+					pvar_local);
+      
+      /* Exchange diffusion limiter */
+      if (df_limiter != NULL) {
+	BFT_MALLOC(df_limiter_local, n_local, cs_real_t);
+	cs_internal_coupling_exchange_var(cpl,
+					  1, /* Dimension */
+					  df_limiter,
+					  df_limiter_local);
+      }
+      
+      /* Flux contribution */
+      assert(f != NULL);
+      cs_real_t *hintp = f->bc_coeffs->hint;
+      cs_real_t *hextp = f->bc_coeffs->hext;
+      for (cs_lnum_t ii = 0; ii < n_local; ii++) {
+	cs_lnum_t face_id = faces_local[ii];
+	cs_lnum_t jj = b_face_cells[face_id];
+	cs_real_t pip, pjp;
+	cs_real_t fluxi = 0.;
+	
+	cs_real_t bldfrp = (cs_real_t) ircflp;
+	/* Local limitation of the reconstruction */
+	if (df_limiter != NULL && ircflp > 0)
+	  bldfrp = CS_MAX(CS_MIN(df_limiter_local[ii], df_limiter[jj]), 0.);
+	
+	cs_b_cd_unsteady(bldfrp,
+			 diipb[face_id],
+			 grad[jj],
+			 _pvar[jj],
+			 &pip);
+	
+	pjp = pvar_local[ii];
+	
+	hint = hintp[face_id];
+	hext = hextp[face_id];
+	heq = hint * hext / (hint + hext);
+	
+	cs_b_diff_flux_coupling(idiffp,
+				pip,
+				pjp,
+				heq,
+				&fluxi);
+	
+	rhs[jj] -= thetap * fluxi;
+      }
+    } 
+    
     /* Unsteady */
   } else {
-
     for (int g_id = 0; g_id < n_b_groups; g_id++) {
 #     pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
       for (int t_id = 0; t_id < n_b_threads; t_id++) {