Index: /LMDZ6/trunk/DefLists/field_def_lmdz.xml
===================================================================
--- /LMDZ6/trunk/DefLists/field_def_lmdz.xml	(revision 5038)
+++ /LMDZ6/trunk/DefLists/field_def_lmdz.xml	(revision 5039)
@@ -636,4 +636,5 @@
 	<field id="tke_buoy" long_name="TKE Buoyancy term" unit="m2/s3" />
 	<field id="tke_shear" long_name="TKE Shear term" unit="m2/s3" />
+	<field id="tke_trans" long_name="TKE Transport term" unit="m2/s3" />
         <field id="tke_ter"    long_name="Max Turb. Kinetic Energy ter"    unit="m2/s2" />
         <field id="tke_lic"    long_name="Max Turb. Kinetic Energy lic"    unit="m2/s2" />
Index: /LMDZ6/trunk/libf/phylmd/lmdz_atke_exchange_coeff.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_atke_exchange_coeff.F90	(revision 5038)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_atke_exchange_coeff.F90	(revision 5039)
@@ -7,5 +7,5 @@
 subroutine atke_compute_km_kh(ngrid,nlay,dtime, &
                         wind_u,wind_v,temp,qvap,play,pinterf,cdrag_uv, &
-                        tke,eps,Km_out,Kh_out)
+                        tke,eps,tke_shear,tke_buoy,tke_trans,Km_out,Kh_out)
 
 !========================================================================
@@ -79,4 +79,7 @@
 
 REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT)    :: eps      ! output: TKE dissipation rate at interface between layers (m2/s3)
+REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT)    :: tke_shear! output: TKE shear production rate (m2/s3)
+REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT)    :: tke_buoy ! output: TKE buoyancy production rate (m2/s3)
+REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT)    :: tke_trans! output: TKE transport (diffusion) term (m2/s3)
 REAL, DIMENSION(ngrid,nlay), INTENT(OUT)      :: Km_out   ! output: Exchange coefficient for momentum at interface between layers (m2/s)
 REAL, DIMENSION(ngrid,nlay), INTENT(OUT)      :: Kh_out   ! output: Exchange coefficient for heat flux at interface between layers (m2/s)
@@ -261,4 +264,7 @@
                 shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))
                 eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay))
+                tke_shear(igrid,ilay)=l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay))*shear2(igrid,ilay) 
+                tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay))*shear2(igrid,ilay) &
+                                    *(Ri(igrid,ilay) / Prandtl(igrid,ilay))
             ENDDO
         ENDDO
@@ -278,5 +284,8 @@
             qq=max(0.,qq)
             tke(igrid,ilay)=0.5*(qq**2)
-            eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay)) 
+            eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay))
+            tke_shear(igrid,ilay)=l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay))*shear2(igrid,ilay)
+            tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay))*shear2(igrid,ilay) &
+                                *(Ri(igrid,ilay) / Prandtl(igrid,ilay)) 
             ENDDO
         ENDDO
@@ -293,5 +302,8 @@
             qq=(qq+l_exchange(igrid,ilay)*Sm(igrid,ilay)*dtime/sqrt(2.)      & 
                 *shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) &
-                /(1.+qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.))) 
+                /(1.+qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))
+            tke_shear(igrid,ilay)=l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay)
+            tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay) &
+                                *(Ri(igrid,ilay) / Prandtl(igrid,ilay)) 
             tke(igrid,ilay)=0.5*(qq**2) 
             eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay))
@@ -308,4 +320,7 @@
             eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay))
             qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
+            tke_shear(igrid,ilay)=l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay)
+            tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay) &
+                                *(Ri(igrid,ilay) / Prandtl(igrid,ilay))
             IF (Ri(igrid,ilay) .LT. 0.) THEN
                 netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))
@@ -327,4 +342,7 @@
             DO igrid=1,ngrid
             qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
+            tke_shear(igrid,ilay)=l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay)
+            tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay) &
+                                *(Ri(igrid,ilay) / Prandtl(igrid,ilay))
             qq=(l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) &
                 +qq*(1.+dtime*qq/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))) &
@@ -349,4 +367,6 @@
     tke(igrid,nlay+1)=0.
     eps(igrid,nlay+1)=0.
+    tke_shear(igrid,nlay+1)=0.
+    tke_buoy(igrid,nlay+1)=0.
 END DO
 
@@ -359,4 +379,6 @@
     tke(igrid,1)=ctkes*(ustar**2)
     eps(igrid,1)=0. ! arbitrary as TKE is not properly defined at the surface
+    tke_shear(igrid,1)=0.
+    tke_buoy(igrid,1)=0.
 END DO
 
@@ -364,6 +386,7 @@
 ! vertical diffusion of TKE 
 !==========================
+tke_trans(:,:)=0.
 IF (atke_ok_vdiff) THEN
-    CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)
+    CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke,tke_trans)
 ENDIF
 
@@ -387,5 +410,5 @@
 
 !===============================================================================================
-subroutine atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)
+subroutine atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke,tke_trans)
 
 ! routine that computes the vertical diffusion of TKE by the turbulence
@@ -408,5 +431,5 @@
 
 REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke    ! turbulent kinetic energy at interface between layers
-
+REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke_trans ! turbulent kinetic energy transport term (m2/s3)
 
 
@@ -480,4 +503,5 @@
 ! update TKE
 tke(:,:)=tke(:,:)+dtke(:,:)
+tke_trans(:,:)=dtke(:,:)/dtime
 
 
Index: /LMDZ6/trunk/libf/phylmd/lmdz_call_atke.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_call_atke.F90	(revision 5038)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_call_atke.F90	(revision 5039)
@@ -8,5 +8,5 @@
 contains
 
-subroutine call_atke(dtime,ngrid,nlay,cdrag_uv,cdrag_t,u_surf,v_surf,temp_surf, &
+subroutine call_atke(dtime,ngrid,nlay,nsrf,ni,cdrag_uv,cdrag_t,u_surf,v_surf,temp_surf, &
                         wind_u,wind_v,temp,qvap,play,pinterf, &
                         tke,eps,Km_out,Kh_out)
@@ -16,4 +16,5 @@
 
 USE lmdz_atke_turbulence_ini, ONLY : iflag_num_atke, rg, rd
+USE phys_local_var_mod, ONLY: tke_shear, tke_buoy, tke_trans
 
 implicit none
@@ -26,4 +27,6 @@
 INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
 INTEGER, INTENT(IN) :: nlay ! number of vertical index  
+INTEGER, INTENT(IN) :: nsrf ! surface tile index
+INTEGER, DIMENSION(ngrid), INTENT(IN) :: ni ! array of indices to move from knon to klon arrays
 
 
@@ -50,12 +53,13 @@
 
 
+REAL, DIMENSION(ngrid,nlay+1) :: tke_shear_term,tke_buoy_term,tke_trans_term
 REAL, DIMENSION(ngrid,nlay) :: wind_u_predict, wind_v_predict
 REAL, DIMENSION(ngrid) ::  wind1
-INTEGER i
+INTEGER i,j,k
 
 
 call atke_compute_km_kh(ngrid,nlay,dtime,&
                         wind_u,wind_v,temp,qvap,play,pinterf,cdrag_uv,&
-                        tke,eps,Km_out,Kh_out)
+                        tke,eps,tke_shear_term,tke_buoy_term,tke_trans_term,Km_out,Kh_out)
 
 
@@ -76,10 +80,19 @@
    call atke_compute_km_kh(ngrid,nlay,dtime,&
                         wind_u_predict,wind_v_predict,temp,qvap,play,pinterf,cdrag_uv, &
-                        tke,eps,Km_out,Kh_out)
+                        tke,eps,tke_shear_term,tke_buoy_term,tke_trans_term,Km_out,Kh_out)
 
 end if
 
 
+! Diagnostics of tke loss/source terms
 
+ DO k=1,nlay+1
+    DO i=1,ngrid
+       j=ni(i)
+       tke_shear(j,k,nsrf)=tke_shear_term(i,k)
+       tke_buoy(j,k,nsrf)=tke_buoy_term(i,k)
+       tke_trans(j,k,nsrf)=tke_trans_term(i,k)
+    ENDDO
+ ENDDO
 
 
Index: /LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 5038)
+++ /LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 5039)
@@ -1996,5 +1996,5 @@
 
         IF (iflag_pbl>=50) THEN
-        CALL call_atke(dtime,knon,klev,ycdragm(1:knon), ycdragh(1:knon),yus0(1:knon),yvs0(1:knon),yts(1:knon), &
+        CALL call_atke(dtime,knon,klev,nsrf,ni,ycdragm(1:knon), ycdragh(1:knon),yus0(1:knon),yvs0(1:knon),yts(1:knon), &
                   yu(1:knon,:),yv(1:knon,:),yt(1:knon,:),yq(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:),       &
                   ytke(1:knon,:),yeps(1:knon,:), ycoefm(1:knon,:), ycoefh(1:knon,:))
@@ -2041,5 +2041,5 @@
         IF (iflag_pbl>=50) THEN
      
-        CALL call_atke(dtime,knon,klev,ycdragm_x(1:knon),ycdragh_x(1:knon),yus0(1:knon),yvs0(1:knon),yts_x(1:knon),    &
+        CALL call_atke(dtime,knon,klev,nsrf,ni,ycdragm_x(1:knon),ycdragh_x(1:knon),yus0(1:knon),yvs0(1:knon),yts_x(1:knon),    &
                        yu_x(1:knon,:),yv_x(1:knon,:),yt_x(1:knon,:),yq_x(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:),  &
                        ytke_x(1:knon,:),yeps_x(1:knon,:),ycoefm_x(1:knon,:), ycoefh_x(1:knon,:))
@@ -2081,5 +2081,5 @@
         IF (iflag_pbl>=50) THEN
         
-        CALL call_atke(dtime,knon,klev,ycdragm_w(1:knon),ycdragh_w(1:knon),yus0(1:knon),yvs0(1:knon),yts_w(1:knon), &
+        CALL call_atke(dtime,knon,klev,nsrf,ni,ycdragm_w(1:knon),ycdragh_w(1:knon),yus0(1:knon),yvs0(1:knon),yts_w(1:knon), &
                 yu_w(1:knon,:),yv_w(1:knon,:),yt_w(1:knon,:),yq_w(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:),      &
                 ytke_w(1:knon,:),yeps_w(1:knon,:),ycoefm_w(1:knon,:),ycoefh_w(1:knon,:))
Index: /LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 5038)
+++ /LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 5039)
@@ -26,4 +26,6 @@
       REAL, SAVE, ALLOCATABLE :: pbl_eps(:,:,:)
       !$OMP THREADPRIVATE(pbl_eps)
+      REAL, SAVE, ALLOCATABLE :: tke_shear(:,:,:), tke_buoy(:,:,:), tke_trans(:,:,:)
+      !$OMP THREADPRIVATE(tke_shear(:,:,:),tke_buoy(:,:,:),tke_trans(:,:,:))
       REAL, SAVE, ALLOCATABLE :: tr_seri(:,:,:)
       !$OMP THREADPRIVATE(tr_seri)
@@ -694,5 +696,7 @@
       ALLOCATE(l_mixmin(klon,klev+1,nbsrf),l_mix(klon,klev+1,nbsrf),wprime(klon,klev+1,nbsrf))
       ALLOCATE(pbl_eps(klon,klev+1,nbsrf+1))
+      ALLOCATE(tke_shear(klon,klev+1,nbsrf), tke_buoy(klon,klev+1,nbsrf), tke_trans(klon,klev+1,nbsrf))
       pbl_eps(:,:,:)=0.
+      tke_shear(:,:,:)=0.; tke_buoy(:,:,:)=0.; tke_trans(:,:,:)=0.
       l_mix(:,:,:)=0.;l_mixmin(:,:,:)=0.;wprime(:,:,:)=0. ! doit etre initialse car pas toujours remplis
       ALLOCATE(rhcl(klon,klev))
@@ -1055,4 +1059,5 @@
       DEALLOCATE(u_seri,v_seri)
       DEALLOCATE(l_mixmin,l_mix,wprime)
+      DEALLOCATE(tke_shear,tke_buoy,tke_trans)
       DEALLOCATE(pbl_eps)
       DEALLOCATE(rhcl)
Index: /LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 5038)
+++ /LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 5039)
@@ -1114,6 +1114,13 @@
   TYPE(ctrl_out), SAVE :: o_tke = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'tke ', 'TKE', 'm2/s2', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_tke_shear = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'tke_shear ', 'TKE shear term', 'm2/s3', (/ ('', i=1, 10) /))  
+  TYPE(ctrl_out), SAVE :: o_tke_buoy = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'tke_buoy ', 'TKE buoyancy term', 'm2/s3', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_tke_trans = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'tke_trans ', 'TKE transport term', 'm2/s3', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_tke_dissip = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
-    'tke_dissip ', 'TKE DISSIPATION', 'm2/s3', (/ ('', i=1, 10) /))   
+    'tke_dissip ', 'TKE dissipation term', 'm2/s3', (/ ('', i=1, 10) /))
+
   TYPE(ctrl_out), SAVE :: o_tke_max = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'tke_max', 'TKE max', 'm2/s2',                                  &
Index: /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 5038)
+++ /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 5039)
@@ -146,5 +146,5 @@
          o_dqsphy, o_dqsphy2d, o_dqbsphy, o_dqbsphy2d, o_albe_srf, o_z0m_srf, o_z0h_srf, &
          o_ages_srf, o_snow_srf, o_alb1, o_alb2, o_tke, o_tke_dissip, &
-         o_tke_max, o_kz, o_kz_max, o_clwcon, &
+         o_tke_max, o_kz, o_kz_max, o_clwcon, o_tke_shear, o_tke_buoy, o_tke_trans,  &
          o_dtdyn, o_dqdyn, o_dqdyn2d, o_dqldyn, o_dqldyn2d, & 
          o_dqsdyn, o_dqsdyn2d, o_dqbsdyn, o_dqbsdyn2d, o_dudyn, o_dvdyn, &
@@ -314,5 +314,5 @@
          zn2mout, t2m_min_mon, t2m_max_mon, evap, &
          snowerosion, zxustartlic, zxrhoslic, zxqsaltlic, &
-         l_mixmin,l_mix, pbl_eps, &
+         l_mixmin,l_mix, pbl_eps, tke_shear, tke_buoy, tke_trans, &
          zu10m, zv10m, zq2m, zustar, zxqsurf, &
          rain_lsc, rain_num, snow_lsc, bils, sens, fder, &
@@ -1313,5 +1313,43 @@
           ENDIF
           
-          CALL histwrite_phy(o_tke_dissip, zx_tmp_fi3d)    
+          CALL histwrite_phy(o_tke_dissip, zx_tmp_fi3d)   
+
+          zx_tmp_fi3d=0.
+          IF (vars_defined) THEN
+             DO nsrf=1,nbsrf
+                DO k=1,klev
+                   zx_tmp_fi3d(:,k)=zx_tmp_fi3d(:,k) &
+                        +pctsrf(:,nsrf)*tke_shear(:,k,nsrf)
+                ENDDO
+             ENDDO
+          ENDIF
+
+          CALL histwrite_phy(o_tke_shear, zx_tmp_fi3d)
+
+          zx_tmp_fi3d=0.
+          IF (vars_defined) THEN
+             DO nsrf=1,nbsrf
+                DO k=1,klev
+                   zx_tmp_fi3d(:,k)=zx_tmp_fi3d(:,k) &
+                        +pctsrf(:,nsrf)*tke_buoy(:,k,nsrf)
+                ENDDO
+             ENDDO
+          ENDIF
+
+          CALL histwrite_phy(o_tke_buoy, zx_tmp_fi3d)
+
+
+          zx_tmp_fi3d=0.
+          IF (vars_defined) THEN
+             DO nsrf=1,nbsrf
+                DO k=1,klev
+                   zx_tmp_fi3d(:,k)=zx_tmp_fi3d(:,k) &
+                        +pctsrf(:,nsrf)*tke_trans(:,k,nsrf)
+                ENDDO
+             ENDDO
+          ENDIF
+
+          CALL histwrite_phy(o_tke_trans, zx_tmp_fi3d)
+
        ENDIF
 
Index: /LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90	(revision 5038)
+++ /LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90	(revision 5039)
@@ -28,4 +28,6 @@
       REAL, SAVE, ALLOCATABLE :: pbl_eps(:,:,:)
       !$OMP THREADPRIVATE(pbl_eps)
+      REAL, SAVE, ALLOCATABLE :: tke_shear(:,:,:), tke_buoy(:,:,:), tke_trans(:,:,:)
+      !$OMP THREADPRIVATE(tke_shear(:,:,:),tke_buoy(:,:,:),tke_trans(:,:,:))
       REAL, SAVE, ALLOCATABLE :: tr_seri(:,:,:)
       !$OMP THREADPRIVATE(tr_seri)
@@ -815,5 +817,7 @@
       ALLOCATE(l_mixmin(klon,klev+1,nbsrf),l_mix(klon,klev+1,nbsrf),wprime(klon,klev+1,nbsrf))
       ALLOCATE(pbl_eps(klon,klev+1,nbsrf+1))
+      ALLOCATE(tke_shear(klon,klev+1,nbsrf), tke_buoy(klon,klev+1,nbsrf), tke_trans(klon,klev+1,nbsrf))
       pbl_eps(:,:,:)=0.
+      tke_shear(:,:,:)=0.; tke_buoy(:,:,:)=0.; tke_trans(:,:,:)=0.
       l_mix(:,:,:)=0.;l_mixmin(:,:,:)=0.;wprime(:,:,:)=0. ! doit etre initialse car pas toujours remplis
       ALLOCATE(rhcl(klon,klev))
@@ -1253,4 +1257,5 @@
       DEALLOCATE(u_seri,v_seri)
       DEALLOCATE(l_mixmin,l_mix,wprime)
+      DEALLOCATE(tke_shear,tke_buoy,tke_trans)
       DEALLOCATE(pbl_eps)
       DEALLOCATE(rhcl)
Index: /LMDZ6/trunk/libf/phylmdiso/phys_output_ctrlout_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/phys_output_ctrlout_mod.F90	(revision 5038)
+++ /LMDZ6/trunk/libf/phylmdiso/phys_output_ctrlout_mod.F90	(revision 5039)
@@ -1114,6 +1114,13 @@
   TYPE(ctrl_out), SAVE :: o_tke = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'tke ', 'TKE', 'm2/s2', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_tke_shear = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'tke_shear ', 'TKE shear term', 'm2/s3', (/ ('', i=1, 10) /))  
+  TYPE(ctrl_out), SAVE :: o_tke_buoy = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'tke_buoy ', 'TKE buoyancy term', 'm2/s3', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_tke_trans = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'tke_trans ', 'TKE transport term', 'm2/s3', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_tke_dissip = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
-    'tke_dissip ', 'TKE DISSIPATION', 'm2/s3', (/ ('', i=1, 10) /))   
+    'tke_dissip ', 'TKE dissipation term', 'm2/s3', (/ ('', i=1, 10) /))
+
   TYPE(ctrl_out), SAVE :: o_tke_max = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'tke_max', 'TKE max', 'm2/s2',                                  &
