Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90	(revision 10236)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90	(revision 10728)
@@ -33,4 +33,7 @@
    USE inputs_fabm
    USE vertical_movement_fabm
+   USE zdfmxl
+   USE asmbgc, ONLY: mld_choice_bgc
+   USE lbclnk
 
    !USE fldread         !  time interpolation
@@ -113,4 +116,6 @@
 
       CALL st2d_fabm_nxt( kt )
+      
+      CALL asmdiags_fabm( kt )
 
       IF( l_trdtrc )  CALL wrk_alloc( jpi, jpj, jpk, ztrfabm )
@@ -130,4 +135,88 @@
 
    END SUBROUTINE trc_sms_fabm
+   
+   SUBROUTINE asmdiags_fabm( kt )
+      INTEGER, INTENT(IN) :: kt
+      INTEGER :: ji,jj,jk,jkmax
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: pgrow_3d, ploss_3d, zmld
+	   
+      IF (kt == nittrc000) THEN
+         MLD_MAX(:,:) = 0.0
+      ENDIF
+      PGROW_AVG(:,:) = 0.0
+      PLOSS_AVG(:,:) = 0.0
+      PHYT_AVG(:,:)  = 0.0
+	     
+      pgrow_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_pgrow)
+      ploss_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_ploss)
+      
+      SELECT CASE( mld_choice_bgc )
+      CASE ( 1 )                   ! Turbocline/mixing depth [W points]
+         zmld(:,:) = hmld(:,:)
+      CASE ( 2 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
+         zmld(:,:) = hmlp(:,:)
+      CASE ( 3 )                   ! Kara MLD [Interpolated]
+#if defined key_karaml
+         IF ( ln_kara ) THEN
+            zmld(:,:) = hmld_kara(:,:)
+         ELSE
+            CALL ctl_stop( ' Kara mixed layer requested for BGC assimilation,', &
+               &           ' but ln_kara=.false.' )
+         ENDIF
+#else
+         CALL ctl_stop( ' Kara mixed layer requested for BGC assimilation,', &
+            &           ' but is not defined' )
+#endif
+      CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points]
+         zmld(:,:) = hmld_tref(:,:)
+      CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
+         zmld(:,:) = hmlpt(:,:)
+      END SELECT
+	
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+            !
+            jkmax = jpk-1
+            DO jk = jpk-1, 1, -1
+               IF ( ( zmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
+                  & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
+                  zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
+                  jkmax = jk
+               ENDIF
+            END DO
+            !
+            DO jk = 1, jkmax
+               PHYT_AVG(ji,jj) = PHYT_AVG(ji,jj) + &
+                  &              trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + &
+                  &              trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + &
+                  &              trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + &
+                  &              trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n)
+               IF ( pgrow_3d(ji,jj,jk) .GT. 0.0 ) THEN
+                  PGROW_AVG(ji,jj) = PGROW_AVG(ji,jj) + &
+                     &               pgrow_3d(ji,jj,jk)
+               ENDIF
+               IF ( ploss_3d(ji,jj,jk) .GT. 0.0 ) THEN
+                  PLOSS_AVG(ji,jj) = PLOSS_AVG(ji,jj) + &
+                     &               ploss_3d(ji,jj,jk)
+               ENDIF
+            END DO
+           
+            PHYT_AVG(ji,jj)  = PHYT_AVG(ji,jj)  / REAL(jkmax)
+            PGROW_AVG(ji,jj) = PGROW_AVG(ji,jj) / REAL(jkmax)
+            PLOSS_AVG(ji,jj) = PLOSS_AVG(ji,jj) / REAL(jkmax)
+	
+            IF ( zmld(ji,jj) .GT. MLD_MAX(ji,jj) ) THEN
+               MLD_MAX(ji,jj) = zmld(ji,jj)
+            ENDIF
+            !
+         END DO
+      END DO
+      
+      PHYT_AVG(:,:)  = PHYT_AVG(:,:)  * tmask(:,:,1)
+      PGROW_AVG(:,:) = PGROW_AVG(:,:) * tmask(:,:,1)
+      PLOSS_AVG(:,:) = PLOSS_AVG(:,:) * tmask(:,:,1)
+      MLD_MAX(:,:)   = MLD_MAX(:,:)   * tmask(:,:,1)
+	   
+   END SUBROUTINE asmdiags_fabm
 
    SUBROUTINE compute_fabm()
