Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/CONFIG/makenemo
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/CONFIG/makenemo	(revision 10155)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/CONFIG/makenemo	(revision 10156)
@@ -360,4 +360,24 @@
     cd ${NEMO_TDIR}/${NEW_CONF} || cd -
 
+    #-
+    #- Get git version information if available
+    mkdir -p ${CONFIG_DIR}/${NEW_CONF}/BLD/inc
+    GITVERSIONHEADER=${CONFIG_DIR}/${NEW_CONF}/BLD/inc/gitversion.h90
+    echo "Extracting git commit info into" $GITVERSIONHEADER
+    NEMOCOMMIT=$( git describe --always --dirty )
+    if [ $? -ne 0 ]
+    then
+       echo "Unable to extract git commit info!"
+       NEMOCOMMIT="NOT AVAILABLE"
+    fi
+    echo "#define _NEMO_COMMIT_ID_ '${NEMOCOMMIT}'" > ${GITVERSIONHEADER}
+    NEMOBRANCH=$( git name-rev --name-only HEAD )
+    if [ $? -ne 0 ]
+    then
+       echo "Unable to extract git branch info!"
+       NEMOBRANCH="UNKNOWN"
+    fi
+    echo "#define _NEMO_BRANCH_ '${NEMOBRANCH}'" >> ${GITVERSIONHEADER}
+
 #if AGRIF we do a first preprocessing
     if [ ${#x_c} -eq 0 ]; then
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/inputs_fabm.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/inputs_fabm.F90	(revision 10156)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/inputs_fabm.F90	(revision 10156)
@@ -0,0 +1,254 @@
+MODULE inputs_fabm
+   !!======================================================================
+   !!                         ***  MODULE inputs_fabm  ***
+   !! TOP :   Input module of the FABM tracers
+   !!======================================================================
+
+#if defined key_fabm
+   !!----------------------------------------------------------------------
+   !!   'key_fabm'                                               FABM tracers
+   !!----------------------------------------------------------------------
+   !! initialize_inputs       : initialize input structures
+   !! update_inputs : update 2D input fields
+   !! trc_rnf_fabm : update river data
+   !!----------------------------------------------------------------------
+   USE par_trc
+   USE oce_trc
+   USE trc
+   USE iom
+   USE fldread
+   USE par_fabm
+   USE fabm
+
+   IMPLICIT NONE
+
+   PRIVATE
+
+   PUBLIC initialize_inputs
+   PUBLIC link_inputs
+   PUBLIC update_inputs
+   PUBLIC trc_rnf_fabm
+
+   #if defined key_trdtrc && defined key_iomput
+      REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:) :: tr_inp
+   #endif
+
+   TYPE, PUBLIC :: type_input_variable
+      TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf
+      INTEGER                              :: ntimes
+   END TYPE
+
+   TYPE, PUBLIC, EXTENDS(type_input_variable) :: type_input_data
+      TYPE(type_horizontal_variable_id)   :: horizontal_id
+      TYPE(type_input_data), POINTER   :: next => null()
+   END TYPE
+   TYPE (type_input_data), POINTER, PUBLIC :: first_input_data => NULL()
+
+   TYPE, PUBLIC, EXTENDS(type_input_variable):: type_river_data
+      INTEGER   :: jp_pos=0 ! position of linked state variable in trc fields
+      REAL(wp) :: rn_trrnfac=1._wp ! unit conversion factor
+      TYPE(type_river_data), POINTER   :: next => null()
+   END TYPE
+   TYPE (type_river_data), POINTER, PUBLIC :: first_river_data => NULL()
+
+   CONTAINS
+
+     SUBROUTINE initialize_inputs
+        TYPE(FLD_N)        :: sn, sn_empty
+        CHARACTER(LEN=256) :: name
+        REAL(wp) :: rfac
+        NAMELIST /variable/ name,sn
+        NAMELIST /riverdata/ name,sn,rfac
+        LOGICAL :: l_ext
+        INTEGER :: num, ierr, nmlunit
+        TYPE (type_input_data),POINTER :: input_data
+        TYPE (type_river_data),POINTER :: river_data
+        INTEGER :: jn
+        INTEGER , PARAMETER :: nbtimes = 366  !: maximum number of times record in a file
+        REAL(wp), DIMENSION(nbtimes) :: zsteps
+
+        ! Check if fabm_input.nml exists - if not, do nothing and return.
+        INQUIRE( FILE='fabm_input.nml', EXIST=l_ext )
+        IF (.NOT.l_ext) return
+
+        ! Open fabm_input.nml
+        CALL ctl_opn( nmlunit, 'fabm_input.nml', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, num, .FALSE. )
+
+        ! Read any number of "variable" namelists
+        DO
+           ! Initialize namelist variables
+           name = ''
+           sn = sn_empty
+
+           ! Read the namelist
+           READ(nmlunit,nml=variable,err=98,end=99)
+
+           ! Transfer namelist settings to new input_data object
+           ALLOCATE(input_data, STAT=ierr)
+           IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'inputs_fabm:initialize_inputs: unable to allocate input_data object for variable '//TRIM(name) )
+           input_data%horizontal_id = fabm_get_horizontal_variable_id(model,name)
+           IF (.NOT.fabm_is_variable_used(input_data%horizontal_id)) THEN
+              ! This variable was not found among FABM's horizontal variables (at least, those that are read by one or more FABM modules)
+              CALL ctl_stop('STOP', 'inputs_fabm:initialize_inputs: variable "'//TRIM(name)//'" was not found among horizontal FABM variables.')
+           END IF
+           ALLOCATE(input_data%sf(1), STAT=ierr)
+           IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'inputs_fabm:initialize_inputs: unable to allocate sf structure for variable '//TRIM(name) )
+           CALL fld_fill(input_data%sf, (/sn/), '', 'inputs_fabm:initialize_inputs', 'FABM variable '//TRIM(name), 'variable' )
+           ALLOCATE( input_data%sf(1)%fnow(jpi,jpj,1)   )
+           IF( sn%ln_tint ) ALLOCATE( input_data%sf(1)%fdta(jpi,jpj,1,2) )
+
+           ! Get number of record in file (if there is only one, we will read data
+           ! only at the very first time step)
+           CALL fld_clopn( input_data%sf(1) )
+           CALL iom_gettime( input_data%sf(1)%num, zsteps, kntime=input_data%ntimes)
+           CALL iom_close( input_data%sf(1)%num )
+
+           ! Prepend new input variable to list.
+           input_data%next => first_input_data
+           first_input_data => input_data
+        END DO
+
+  98    CALL ctl_stop('STOP', 'inputs_fabm:initialize_inputs: unable to read namelist "riverdata"')
+
+  99    REWIND(nmlunit)
+
+        ! Read any number of "riverdata" namelists
+        DO
+           ! Initialize namelist variables
+           name = ''
+           sn = sn_empty
+           rfac = 1._wp
+
+           ! Read the namelist
+           READ(nmlunit,nml=riverdata,err=198,end=199)
+
+           ! Transfer namelist settings to new river_data object
+           ALLOCATE(river_data, STAT=ierr)
+           IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'inputs_fabm:initialize_inputs: unable to allocate river_data object for variable '//TRIM(name) )
+           ! Check if river data name is in FABM states and
+           ! provide NEMO with position of the respective state variable
+           ! within tracer field
+           DO jn=1,jp_fabm
+             IF (TRIM(name) == TRIM(model%state_variables(jn)%name)) THEN
+               river_data%jp_pos = jp_fabm_m1+jn
+             END IF
+           END DO
+           IF (river_data%jp_pos == 0) THEN
+             ! This variable was not found among FABM's state variables
+             ! passed to NEMO!
+             CALL ctl_stop('STOP', 'inputs_fabm:initialize_inputs: variable "'//TRIM(name)//'" was not found among FABM state variables.')
+           END IF
+
+           ALLOCATE(river_data%sf(1), STAT=ierr)
+           IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'inputs_fabm:initialize_inputs: unable to allocate sf structure for variable '//TRIM(name) )
+           CALL fld_fill(river_data%sf, (/sn/), '', 'inputs_fabm:initialize_inputs', 'FABM variable '//TRIM(name), 'riverdata' )
+           ALLOCATE( river_data%sf(1)%fnow(jpi,jpj,1)   )
+           IF( sn%ln_tint ) ALLOCATE( river_data%sf(1)%fdta(jpi,jpj,1,2) )
+
+           ! Load unit conversion factor:
+           river_data%rn_trrnfac=rfac
+
+           ! Get number of record in file (if there is only one, we will read data
+           ! only at the very first time step)
+           CALL fld_clopn( river_data%sf(1) )
+           CALL iom_gettime( river_data%sf(1)%num, zsteps, kntime=river_data%ntimes)
+           CALL iom_close( river_data%sf(1)%num )
+
+           ! Prepend new input variable to list.
+           river_data%next => first_river_data
+           first_river_data => river_data
+        END DO
+
+  198   CALL ctl_stop('STOP', 'inputs_fabm:initialize_inputs: unable to read namelist "riverdata"')
+
+  199   RETURN
+
+     END SUBROUTINE initialize_inputs
+
+     SUBROUTINE link_inputs
+      TYPE (type_input_data),POINTER :: input_data
+
+      input_data => first_input_data
+      DO WHILE (ASSOCIATED(input_data))
+         ! Provide FABM with pointer to field that will receive prescribed data.
+         ! NB source=data_source_user guarantees that the prescribed data takes priority over any data FABM may already have for that variable.
+         CALL fabm_link_horizontal_data(model,input_data%horizontal_id,input_data%sf(1)%fnow(:,:,1),source=data_source_user)
+         input_data => input_data%next
+      END DO
+
+   END SUBROUTINE link_inputs
+
+   SUBROUTINE update_inputs( kt , l_write)
+      INTEGER, INTENT(IN) :: kt
+      LOGICAL, INTENT(IN), OPTIONAL :: l_write
+      TYPE (type_input_data),POINTER :: input_data
+      TYPE (type_river_data),POINTER :: river_data
+
+      input_data => first_input_data
+      DO WHILE (ASSOCIATED(input_data))
+         IF( kt == nit000 .OR. ( kt /= nit000 .AND. input_data%ntimes > 1 ) ) CALL fld_read( kt, 1, input_data%sf )
+#if defined key_trdtrc && defined key_iomput
+         IF ( .NOT.PRESENT(l_write).OR.l_write ) CALL iom_put( 'INP_'//TRIM(input_data%sf(1)%clvar), input_data%sf(1)%fnow(:,:,1)*tmask(:,:,1) )
+#endif
+         input_data => input_data%next
+      END DO
+
+      river_data => first_river_data
+      DO WHILE (ASSOCIATED(river_data))
+         IF( kt == nit000 .OR. ( kt /= nit000 .AND. river_data%ntimes > 1 ) ) CALL fld_read( kt, 1, river_data%sf )
+         river_data => river_data%next
+      END DO
+
+   END SUBROUTINE update_inputs
+
+   SUBROUTINE trc_rnf_fabm( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE trc_rnf_fabm  ***
+      !!
+      !! ** Purpose :   Add river loadings of biogeochemistry to states
+      !!
+      !! ** Action  :   tra (sms) updated with loadings at time-step kt
+      !!
+      !! This routines assumes river loadings to be given in
+      !! state variable units * m3 / sec
+      !!--------------------------------------------------------------------
+
+      INTEGER, INTENT(in) ::   kt          ! ocean time step
+      REAL(wp) :: zcoef
+      INTEGER :: ji,jj,jk
+      !
+      TYPE (type_river_data),POINTER :: river_data
+
+      river_data => first_river_data
+      DO WHILE (ASSOCIATED(river_data))
+#if defined key_trdtrc && defined key_iomput
+        tr_inp = 0.0_wp
+#endif
+        IF( kt == nit000 .OR. ( kt /= nit000 ) ) THEN
+            DO jj = 1, jpj
+              DO ji = 1, jpi
+                ! convert units and divide by surface area
+                ! loading / cell volume * vertical fraction of riverload
+                ! dtrc / dt (river) = riverload / e1e2t / e3t * e3t * h_rnf
+                !                    = riverload / e1e2t / h_rnf
+                zcoef = river_data%rn_trrnfac / e1e2t(ji,jj) / h_rnf(ji,jj)
+                DO jk = 1,nk_rnf(ji,jj)
+                  ! Add river loadings
+                  tra(ji,jj,jk,river_data%jp_pos) = tra(ji,jj,jk,river_data%jp_pos) + river_data%sf(1)%fnow(ji,jj,1)*zcoef
+#if defined key_trdtrc && defined key_iomput
+                  tr_inp(ji,jj,jk) = river_data%sf(1)%fnow(ji,jj,1)*zcoef
+#endif
+                END DO
+              END DO
+            END DO
+#if defined key_trdtrc && defined key_iomput
+            CALL iom_put( 'INP_'//TRIM(river_data%sf(1)%clvar), tr_inp(:,:,:) )
+#endif
+        END IF
+        river_data => river_data%next
+      END DO
+
+   END SUBROUTINE trc_rnf_fabm
+#endif
+   !!======================================================================
+END MODULE inputs_fabm
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/par_fabm.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/par_fabm.F90	(revision 10156)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/par_fabm.F90	(revision 10156)
@@ -0,0 +1,27 @@
+MODULE par_fabm
+
+   USE fabm
+
+   IMPLICIT NONE
+
+   TYPE (type_model) :: model !FABM model instance
+
+   INTEGER, PUBLIC :: jp_fabm0, jp_fabm1, jp_fabm, &
+                      jp_fabm_surface, jp_fabm_bottom, &
+                      jp_fabm_m1
+
+#if defined key_fabm
+   !!---------------------------------------------------------------------
+   !!   'key_fabm'                     FABM tracers
+   !!---------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER ::   lk_fabm     = .TRUE.   !: FABM flag 
+   LOGICAL, PUBLIC, ALLOCATABLE, DIMENSION(:) ::   lk_rad_fabm !: FABM negativity correction flag array 
+#else
+   !!---------------------------------------------------------------------
+   !!   Default                           No user defined tracers (FABM)
+   !!---------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER ::   lk_fabm     = .FALSE.  !: FABM flag 
+#endif
+
+   !!======================================================================
+END MODULE par_fabm
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/st2D_fabm.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/st2D_fabm.F90	(revision 10156)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/st2D_fabm.F90	(revision 10156)
@@ -0,0 +1,27 @@
+MODULE st2D_fabm
+   !!======================================================================
+   !!                      ***  MODULE  trc  ***
+   !! Passive tracers   :  module for tracers defined
+   !!======================================================================
+#if defined key_top
+   !!----------------------------------------------------------------------
+   !!   'key_top'                                                TOP models
+   !!----------------------------------------------------------------------
+   USE par_oce
+
+   IMPLICIT NONE
+   PUBLIC
+   !! passive tracers fields (before,now,after)
+   !! --------------------------------------------------
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)         ::  fabm_st2Dn           !: 2D state concentration for now time step
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)         ::  fabm_st2Da           !: 2D state concentration for next time step
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)         ::  fabm_st2Db           !: 2D state concentration for before time step
+#else
+   !!----------------------------------------------------------------------
+   !!  Empty module :                                     No passive tracer
+   !!----------------------------------------------------------------------
+#endif
+
+   !!======================================================================
+END MODULE st2d_fabm
+   
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/trcini_fabm.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/trcini_fabm.F90	(revision 10156)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/trcini_fabm.F90	(revision 10156)
@@ -0,0 +1,327 @@
+MODULE trcini_fabm
+   !!======================================================================
+   !!                         ***  MODULE trcini_fabm  ***
+   !! TOP :   initialisation of the FABM tracers
+   !!======================================================================
+   !! History :   2.0  !  2007-12  (C. Ethe, G. Madec) Original code
+   !!----------------------------------------------------------------------
+#if defined key_fabm
+   !!----------------------------------------------------------------------
+   !!   'key_fabm'                                               FABM tracers
+   !!----------------------------------------------------------------------
+   !! trc_ini_fabm   : FABM model initialisation
+   !!----------------------------------------------------------------------
+   USE par_trc         ! TOP parameters
+   USE oce_trc
+   USE trc
+   USE par_fabm
+   USE trcsms_fabm
+   USE fabm_config,ONLY: fabm_create_model_from_yaml_file
+   USE fabm,ONLY: type_external_variable, fabm_initialize_library
+   USE inputs_fabm,ONLY: initialize_inputs,link_inputs, &
+     type_input_variable,type_input_data,type_river_data, &
+     first_input_data,first_river_data
+#if defined key_git_version
+   USE fabm_version,ONLY: fabm_commit_id=>git_commit_id, &
+                          fabm_branch_name=>git_branch_name
+   USE fabm_types,ONLY: type_version,first_module_version
+#endif
+
+
+   IMPLICIT NONE
+   PRIVATE
+
+#if defined key_git_version
+#include "gitversion.h90"
+   CHARACTER(len=*),parameter :: git_commit_id = _NEMO_COMMIT_ID_
+   CHARACTER(len=*),parameter :: git_branch_name = _NEMO_BRANCH_
+#endif
+
+   PUBLIC   trc_ini_fabm   ! called by trcini.F90 module
+   PUBLIC   nemo_fabm_init
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE nemo_fabm_init()
+      INTEGER :: jn
+      INTEGER, PARAMETER :: xml_unit = 1979
+      TYPE (type_input_data),POINTER :: input_data
+      TYPE (type_river_data),POINTER :: river_data
+      CLASS (type_input_variable),POINTER :: input_pointer
+
+      ! Allow FABM to parse fabm.yaml. This ensures numbers of variables are known.
+      call fabm_create_model_from_yaml_file(model)
+
+      jp_fabm = size(model%state_variables)
+      jp_fabm_bottom = size(model%bottom_state_variables)
+      jp_fabm_surface = size(model%surface_state_variables)
+      jp_fabm0 = jptra + 1
+      jp_fabm1 = jptra + jp_fabm
+      jp_fabm_m1=jptra
+      jptra = jptra + jp_fabm
+      jpdia2d = jpdia2d + size(model%horizontal_diagnostic_variables)
+      jpdia3d = jpdia3d + size(model%diagnostic_variables)
+      jpdiabio = jpdiabio + jp_fabm
+
+      !Initialize input data structures.
+      call initialize_inputs
+
+      IF (lwp) THEN
+         ! write field_def_fabm.xml on lead process
+         OPEN(UNIT=xml_unit,FILE='field_def_fabm.xml',ACTION='WRITE',STATUS='REPLACE')
+
+         WRITE (xml_unit,1000) '<field_definition level="1" prec="4" operation="average" enabled=".TRUE." default_value="1.e20" >'
+
+         WRITE (xml_unit,1000) ' <field_group id="ptrc_T" grid_ref="grid_T_3D">'
+         DO jn=1,jp_fabm
+            CALL write_variable_xml(xml_unit,model%state_variables(jn))
+#if defined key_trdtrc
+            CALL write_trends_xml(xml_unit,model%state_variables(jn))
+#endif
+         END DO
+         WRITE (xml_unit,1000) ' </field_group>'
+
+         WRITE (xml_unit,1000) ' <field_group id="sf_T" grid_ref="grid_T_2D">'
+         DO jn=1,jp_fabm_surface
+            CALL write_variable_xml(xml_unit,model%surface_state_variables(jn))
+         END DO
+         DO jn=1,jp_fabm_bottom
+            CALL write_variable_xml(xml_unit,model%bottom_state_variables(jn))
+         END DO
+         WRITE (xml_unit,1000) ' </field_group>'
+
+         WRITE (xml_unit,1000) ' <field_group id="diad_T" grid_ref="grid_T_2D">'
+         DO jn=1,size(model%diagnostic_variables)
+            CALL write_variable_xml(xml_unit,model%diagnostic_variables(jn),3)
+         END DO
+         DO jn=1,size(model%horizontal_diagnostic_variables)
+            CALL write_variable_xml(xml_unit,model%horizontal_diagnostic_variables(jn))
+         END DO
+         WRITE (xml_unit,1000) ' </field_group>'
+
+         WRITE (xml_unit,1000) ' <field_group id="fabm_scalar" grid_ref="grid_0">'
+         DO jn=1,size(model%conserved_quantities)
+            CALL write_variable_xml(xml_unit,model%conserved_quantities(jn))
+         END DO
+         WRITE (xml_unit,1000) ' </field_group>'
+
+         WRITE (xml_unit,1000) ' <field_group id="fabm_input" grid_ref="grid_T_2D">'
+         input_data => first_input_data
+         DO WHILE (ASSOCIATED(input_data))
+           input_pointer => input_data
+           CALL write_input_xml(xml_unit,input_pointer)
+            input_data => input_data%next
+         END DO
+         river_data => first_river_data
+         DO WHILE (ASSOCIATED(river_data))
+           input_pointer => river_data
+           CALL write_input_xml(xml_unit,input_pointer,3)
+            river_data => river_data%next
+         END DO
+         WRITE (xml_unit,1000) ' </field_group>'
+
+         WRITE (xml_unit,1000) '</field_definition>'
+
+         CLOSE(xml_unit)
+      END IF
+      IF( lk_mpp )   CALL mppsync !Ensure field_def_fabm is ready.
+
+1000 FORMAT (A)
+
+   END SUBROUTINE nemo_fabm_init
+
+   SUBROUTINE write_variable_xml(xml_unit,variable,flag_grid_ref)
+      INTEGER,INTENT(IN) :: xml_unit
+      INTEGER,INTENT(IN),OPTIONAL :: flag_grid_ref
+      CLASS (type_external_variable),INTENT(IN) :: variable
+
+      CHARACTER(LEN=20) :: missing_value,string_dimensions
+      INTEGER :: number_dimensions
+
+      ! Check variable dimension for grid_ref specificaiton.
+      ! Default is to not specify the grid_ref in the field definition.
+      IF (present(flag_grid_ref)) THEN
+          number_dimensions=flag_grid_ref
+      ELSE
+          number_dimensions=-1 !default, don't specify grid_ref
+      ENDIF
+
+      WRITE (missing_value,'(E9.3)') variable%missing_value
+      WRITE (string_dimensions,'(I1)') number_dimensions
+      SELECT CASE (number_dimensions)
+      CASE (3)
+         WRITE (xml_unit,'(A)') '  <field id="'//TRIM(variable%name)//'" long_name="'//TRIM(variable%long_name)//'" unit="'//TRIM(variable%units)//'" default_value="'//TRIM(ADJUSTL(missing_value))//'" grid_ref="grid_T_3D" />'
+      CASE (2)
+         WRITE (xml_unit,'(A)') '  <field id="'//TRIM(variable%name)//'" long_name="'//TRIM(variable%long_name)//'" unit="'//TRIM(variable%units)//'" default_value="'//TRIM(ADJUSTL(missing_value))//'" grid_ref="grid_T_2D"/>'
+      CASE (0)
+         WRITE (xml_unit,'(A)') '  <field id="'//TRIM(variable%name)//'" long_name="'//TRIM(variable%long_name)//'" unit="'//TRIM(variable%units)//'" default_value="'//TRIM(ADJUSTL(missing_value))//'" grid_ref="1point"/>'
+      CASE (-1)
+         WRITE (xml_unit,'(A)') '  <field id="'//TRIM(variable%name)//'" long_name="'//TRIM(variable%long_name)//'" unit="'//TRIM(variable%units)//'" default_value="'//TRIM(ADJUSTL(missing_value))//'" />'
+      CASE default
+         IF(lwp) WRITE(numout,*) ' trc_ini_fabm: Failing to initialise output of variable '//TRIM(variable%name)//': Output of '//TRIM(ADJUSTL(string_dimensions))//'-dimensional variables not supported!!!'
+      END SELECT
+
+   END SUBROUTINE write_variable_xml
+
+   SUBROUTINE write_trends_xml(xml_unit,variable,flag_grid_ref)
+      INTEGER,INTENT(IN) :: xml_unit
+      INTEGER,INTENT(IN),OPTIONAL :: flag_grid_ref
+      CLASS (type_external_variable),INTENT(IN) :: variable
+
+      INTEGER :: number_dimensions,i
+      CHARACTER(LEN=20) :: missing_value,string_dimensions
+#if defined key_tracer_budget
+      CHARACTER(LEN=3),DIMENSION(10),PARAMETER :: trd_tags = (/ &
+        'LDF','BBL','FOR','ZDF','DMP','SMS','ATF', &
+        'RDB','RDN','VMV' /)
+      CHARACTER(LEN=3),DIMENSION(3),PARAMETER :: trd_e3t_tags = (/ &
+        'XAD','YAD','ZAD' /)
+#else
+      CHARACTER(LEN=3),DIMENSION(13),PARAMETER :: trd_tags = (/ &
+        'XAD','YAD','ZAD','LDF','BBL','FOR','ZDF','DMP','SMS','ATF', &
+        'RDB','RDN','VMV' /)
+#endif
+
+      ! Check variable dimension for grid_ref specificaiton.
+      ! Default is to not specify the grid_ref in the field definition.
+      IF (present(flag_grid_ref)) THEN
+          number_dimensions=flag_grid_ref
+      ELSE
+          number_dimensions=-1 !default, don't specify grid_ref
+      ENDIF
+
+      WRITE (missing_value,'(E9.3)') -2.E20
+      WRITE (string_dimensions,'(I1)') number_dimensions
+      SELECT CASE (number_dimensions)
+      CASE (3)
+        DO i=1,size(trd_tags)
+         WRITE (xml_unit,'(A)') '  <field id="'//TRIM(trd_tags(i))//'_'//TRIM(variable%name)//'" long_name="'//TRIM(variable%long_name)//' '//TRIM(trd_tags(i))//' trend" unit="'//TRIM(variable%units)//'/s" default_value="'//TRIM(ADJUSTL(missing_value))//'" grid_ref="grid_T_3D" />'
+        END DO
+#if defined key_tracer_budget
+        DO i=1,size(trd_e3t_tags)
+         WRITE (xml_unit,'(A)') '  <field id="'//TRIM(trd_e3t_tags(i))//'_'//TRIM(variable%name)//'_e3t" long_name="'//TRIM(variable%long_name)//' cell depth integrated '//TRIM(trd_e3t_tags(i))//' trend" unit="'//TRIM(variable%units)//'/s" default_value="'//TRIM(ADJUSTL(missing_value))//'" grid_ref="grid_T_3D" />'
+        END DO
+        WRITE (xml_unit,'(A)') '  <field id="'//TRIM(variable%name)//'_e3t" long_name="'//TRIM(variable%long_name)//' cell depth integrated" unit="'//TRIM(variable%units)//'/s" default_value="'//TRIM(ADJUSTL(missing_value))//'" grid_ref="grid_T_3D" />'
+#endif
+      CASE (-1)
+        DO i=1,size(trd_tags)
+         WRITE (xml_unit,'(A)') '  <field id="'//TRIM(trd_tags(i))//'_'//TRIM(variable%name)//'" long_name="'//TRIM(variable%long_name)//' '//TRIM(trd_tags(i))//' trend" unit="'//TRIM(variable%units)//'/s" default_value="'//TRIM(ADJUSTL(missing_value))//'" />'
+        END DO
+#if defined key_tracer_budget
+        DO i=1,size(trd_e3t_tags)
+         WRITE (xml_unit,'(A)') '  <field id="'//TRIM(trd_e3t_tags(i))//'_'//TRIM(variable%name)//'_e3t" long_name="'//TRIM(variable%long_name)//' cell depth integrated '//TRIM(trd_e3t_tags(i))//' trend" unit="'//TRIM(variable%units)//'/s" default_value="'//TRIM(ADJUSTL(missing_value))//'" />'
+        END DO
+        WRITE (xml_unit,'(A)') '  <field id="'//TRIM(variable%name)//'_e3t" long_name="'//TRIM(variable%long_name)//' cell depth integrated" unit="'//TRIM(variable%units)//'/s" default_value="'//TRIM(ADJUSTL(missing_value))//'" />'
+#endif
+      CASE default
+         IF(lwp) WRITE(numout,*) ' trc_ini_fabm: Failing to initialise trends of variable '//TRIM(variable%name)//': Output of '//TRIM(ADJUSTL(string_dimensions))//'-dimensional trends not supported!!!'
+      END SELECT
+
+   END SUBROUTINE write_trends_xml
+
+   SUBROUTINE write_input_xml(xml_unit,variable,flag_grid_ref)
+      INTEGER,INTENT(IN) :: xml_unit
+      INTEGER,INTENT(IN),OPTIONAL :: flag_grid_ref
+      CLASS(type_input_variable),POINTER,INTENT(IN) :: variable
+
+      INTEGER :: number_dimensions,i
+      CHARACTER(LEN=20) :: missing_value,string_dimensions
+
+      ! Check variable dimension for grid_ref specificaiton.
+      ! Default is to not specify the grid_ref in the field definition.
+      IF (present(flag_grid_ref)) THEN
+          number_dimensions=flag_grid_ref
+      ELSE
+          number_dimensions=-1 !default, don't specify grid_ref
+      ENDIF
+
+      WRITE (missing_value,'(E9.3)') -2.E20
+      WRITE (string_dimensions,'(I1)') number_dimensions
+      SELECT CASE (number_dimensions)
+      CASE (3)
+        WRITE (xml_unit,'(A)') '  <field id="'//'INP_'//TRIM(variable%sf(1)%clvar)//'" long_name="'//TRIM(variable%sf(1)%clvar)//' input" unit="" default_value="'//TRIM(ADJUSTL(missing_value))//'" grid_ref="grid_T_3D" />'
+      CASE (-1)
+        WRITE (xml_unit,'(A)') '  <field id="'//'INP_'//TRIM(variable%sf(1)%clvar)//'" long_name="'//TRIM(variable%sf(1)%clvar)//' input" unit="" default_value="'//TRIM(ADJUSTL(missing_value))//'" />'
+      CASE default
+         IF(lwp) WRITE(numout,*) ' trc_ini_fabm: Failing to initialise input diagnostic of variable '//TRIM(variable%sf(1)%clvar)//': Output of '//TRIM(ADJUSTL(string_dimensions))//'-dimensional diagnostic not supported!!!'
+      END SELECT
+
+   END SUBROUTINE write_input_xml
+
+   SUBROUTINE trc_ini_fabm
+      !!----------------------------------------------------------------------
+      !!                     ***  trc_ini_fabm  ***
+      !!
+      !! ** Purpose :   initialization for FABM model
+      !!
+      !! ** Method  : - Read the namcfc namelist and check the parameter values
+      !!----------------------------------------------------------------------
+#if defined key_git_version
+      TYPE (type_version),POINTER :: version
+#endif
+      INTEGER :: jn
+
+      !                       ! Allocate FABM arrays
+      IF( trc_sms_fabm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'trc_ini_fabm: unable to allocate FABM arrays' )
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' trc_ini_fabm: initialisation of FABM model'
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
+#if defined key_git_version
+      IF(lwp) WRITE(numout,*) ' NEMO version:   ',git_commit_id,' (',git_branch_name,' branch)'
+      IF(lwp) WRITE(numout,*) ' FABM version:   ',fabm_commit_id,' (',fabm_branch_name,' branch)'
+#endif
+
+      call fabm_initialize_library()
+#if defined key_git_version
+      version => first_module_version
+
+      do while (associated(version))
+         IF(lwp) WRITE(numout,*)  ' '//trim(version%module_name)//' version:   ',trim(version%version_string)
+         version => version%next
+      end do
+#endif
+
+      ! Log mapping of FABM states:
+      IF (lwp) THEN
+         IF (jp_fabm.gt.0) WRITE(numout,*) " FABM tracers:"
+         DO jn=1,jp_fabm
+            WRITE(numout,*) "   State",jn,":",trim(model%state_variables(jn)%name), &
+               " (",trim(model%state_variables(jn)%long_name), &
+               ") [",trim(model%state_variables(jn)%units),"]"
+         ENDDO
+         IF (jp_fabm_surface.gt.0) WRITE(numout,*) "FABM seasurface states:"
+         DO jn=1,jp_fabm_surface
+            WRITE(numout,*) "   State",jn,":",trim(model%surface_state_variables(jn)%name), &
+               " (",trim(model%surface_state_variables(jn)%long_name), &
+               ") [",trim(model%surface_state_variables(jn)%units),"]"
+         ENDDO
+         IF (jp_fabm_bottom.gt.0) WRITE(numout,*) "FABM seafloor states:"
+         DO jn=1,jp_fabm_bottom
+            WRITE(numout,*) "   State",jn,":",trim(model%bottom_state_variables(jn)%name), &
+               " (",trim(model%bottom_state_variables(jn)%long_name), &
+               ") [",trim(model%bottom_state_variables(jn)%units),"]"
+         ENDDO
+      ENDIF
+
+   END SUBROUTINE trc_ini_fabm
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Dummy module                                        No FABM model
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE nemo_fabm_init
+   END SUBROUTINE nemo_fabm_init
+
+   SUBROUTINE trc_ini_fabm            ! Empty routine
+   END SUBROUTINE trc_ini_fabm
+#endif
+
+   !!======================================================================
+END MODULE trcini_fabm
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/trcnam_fabm.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/trcnam_fabm.F90	(revision 10156)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/trcnam_fabm.F90	(revision 10156)
@@ -0,0 +1,63 @@
+MODULE trcnam_fabm
+   !!======================================================================
+   !!                      ***  MODULE trcnam_fabm  ***
+   !! TOP :   initialisation of some run parameters for FABM bio-model
+   !!======================================================================
+   !! History :   2.0  !  2007-12  (C. Ethe, G. Madec) Original code
+   !!----------------------------------------------------------------------
+#if defined key_fabm
+   !!----------------------------------------------------------------------
+   !!   'key_fabm'   :                                       FABM model
+   !!----------------------------------------------------------------------
+   !! trc_nam_fabm      : FABM initialisation
+   !!----------------------------------------------------------------------
+   USE oce_trc         ! Ocean variables
+   USE par_trc         ! TOP parameters
+   USE trc             ! TOP variables
+
+   USE par_fabm
+   USE trcsms_fabm
+
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   trc_nam_fabm   ! called by trcnam.F90 module
+   PUBLIC   trc_nam_fabm_override ! called by trcnam.F90 module
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
+   !! $Id$ 
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE trc_nam_fabm
+   END SUBROUTINE trc_nam_fabm
+
+   SUBROUTINE trc_nam_fabm_override
+      INTEGER :: jn
+
+      DO jn=1,jp_fabm
+         ctrcnm(jp_fabm_m1+jn) = model%state_variables(jn)%name
+         ctrcln(jp_fabm_m1+jn) = model%state_variables(jn)%long_name
+         ctrcun(jp_fabm_m1+jn) = model%state_variables(jn)%units
+         ln_trc_ini(jp_fabm_m1+jn) = .FALSE.
+      END DO
+   END SUBROUTINE trc_nam_fabm_override
+   
+#else
+   !!----------------------------------------------------------------------
+   !!  Dummy module :                                             No FABM
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_nam_fabm                      ! Empty routine
+   END  SUBROUTINE  trc_nam_fabm
+
+   SUBROUTINE trc_nam_fabm_override
+   END SUBROUTINE trc_nam_fabm_override
+#endif  
+
+   !!======================================================================
+END MODULE trcnam_fabm
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/trcrst_fabm.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/trcrst_fabm.F90	(revision 10156)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/trcrst_fabm.F90	(revision 10156)
@@ -0,0 +1,90 @@
+MODULE trcrst_fabm
+   !!======================================================================
+   !!                      ***  MODULE trcrst_fabm  ***
+   !! Read and write additional restart fields used by FABM
+   !!======================================================================
+   !! History :
+   !!----------------------------------------------------------------------
+#if defined key_fabm
+   !!----------------------------------------------------------------------
+   !!   'key_fabm'   :                                       FABM model
+   !!----------------------------------------------------------------------
+   !! trc_nam_fabm      : FABM initialisation
+   !!----------------------------------------------------------------------
+   USE oce_trc         ! Ocean variables
+   USE par_trc         ! TOP parameters
+   USE trc             ! TOP variables
+   USE iom
+
+   USE par_fabm
+   USE trcsms_fabm
+   USE st2D_fabm
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   trc_rst_read_fabm   ! called by trcrst.F90 module
+   PUBLIC   trc_rst_wri_fabm    ! called by trcrst.F90 module
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE trc_rst_read_fabm
+      INTEGER :: jn
+
+      DO jn=1,jp_fabm_surface
+         CALL iom_get( numrtr, jpdom_autoglo, 'fabm_st2Db'//TRIM(model%surface_state_variables(jn)%name), fabm_st2Db(:,:,jn) )
+      END DO
+      DO jn=1,jp_fabm_bottom
+         CALL iom_get( numrtr, jpdom_autoglo, 'fabm_st2Db'//TRIM(model%bottom_state_variables(jn)%name), fabm_st2Db(:,:,jp_fabm_surface+jn) )
+      END DO
+
+      DO jn=1,jp_fabm_surface
+         CALL iom_get( numrtr, jpdom_autoglo, 'fabm_st2Dn'//TRIM(model%surface_state_variables(jn)%name), fabm_st2Dn(:,:,jn) )
+      END DO
+      DO jn=1,jp_fabm_bottom
+         CALL iom_get( numrtr, jpdom_autoglo, 'fabm_st2Dn'//TRIM(model%bottom_state_variables(jn)%name), fabm_st2Dn(:,:,jp_fabm_surface+jn) )
+      END DO
+   END SUBROUTINE trc_rst_read_fabm
+
+   SUBROUTINE trc_rst_wri_fabm(kt)
+      INTEGER, INTENT( in ) ::   kt    ! ocean time-step index
+
+      INTEGER :: jn
+
+      DO jn=1,jp_fabm_surface
+         CALL iom_rstput( kt, nitrst, numrtw, 'fabm_st2Db'//TRIM(model%surface_state_variables(jn)%name), fabm_st2Db(:,:,jn) )
+      END DO
+      DO jn=1,jp_fabm_bottom
+         CALL iom_rstput( kt, nitrst, numrtw, 'fabm_st2Db'//TRIM(model%bottom_state_variables(jn)%name), fabm_st2Db(:,:,jp_fabm_surface+jn) )
+      END DO
+
+      DO jn=1,jp_fabm_surface
+         CALL iom_rstput( kt, nitrst, numrtw, 'fabm_st2Dn'//TRIM(model%surface_state_variables(jn)%name), fabm_st2Dn(:,:,jn) )
+      END DO
+      DO jn=1,jp_fabm_bottom
+         CALL iom_rstput( kt, nitrst, numrtw, 'fabm_st2Dn'//TRIM(model%bottom_state_variables(jn)%name), fabm_st2Dn(:,:,jp_fabm_surface+jn) )
+      END DO
+
+   END SUBROUTINE trc_rst_wri_fabm
+
+#else
+   !!----------------------------------------------------------------------
+   !!  Dummy module :                                             No FABM
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_rst_read_fabm
+   END  SUBROUTINE trc_rst_read_fabm
+
+   SUBROUTINE trc_rst_wri_fabm(kt)
+      INTEGER, INTENT( in ) ::   kt    ! ocean time-step index
+   END SUBROUTINE trc_rst_wri_fabm
+#endif
+
+   !!======================================================================
+END MODULE trcrst_fabm
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90	(revision 10156)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90	(revision 10156)
@@ -0,0 +1,498 @@
+MODULE trcsms_fabm
+   !!======================================================================
+   !!                         ***  MODULE trcsms_fabm  ***
+   !! TOP :   Main module of the FABM tracers
+   !!======================================================================
+   !! History :   1.0  !  2015-04  (PML) Original code
+   !!----------------------------------------------------------------------
+#if defined key_fabm
+   !!----------------------------------------------------------------------
+   !!   'key_fabm'                                               FABM tracers
+   !!----------------------------------------------------------------------
+   !! trc_sms_fabm       : FABM model main routine
+   !! trc_sms_fabm_alloc : allocate arrays specific to FABM sms
+   !!----------------------------------------------------------------------
+   USE par_trc         ! TOP parameters
+   USE oce_trc         ! Ocean variables
+   USE trc             ! TOP variables
+   USE trcbc
+   USE trd_oce
+   USE trdtrc
+#if defined key_trdtrc && defined key_iomput
+   USE trdtrc_oce
+#endif
+
+   USE oce, only: tsn  ! Needed?
+   USE sbc_oce, only: lk_oasis
+   USE dom_oce
+   USE zdf_oce
+   !USE iom
+   USE xios
+   USE cpl_oasis3
+   USE st2D_fabm
+   USE inputs_fabm
+   USE vertical_movement_fabm
+
+   !USE fldread         !  time interpolation
+
+   USE fabm
+
+   IMPLICIT NONE
+
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+
+   PRIVATE
+
+   PUBLIC   trc_sms_fabm       ! called by trcsms.F90 module
+   PUBLIC   trc_sms_fabm_alloc ! called by trcini_fabm.F90 module
+   PUBLIC   trc_sms_fabm_check_mass
+   PUBLIC   st2d_fabm_nxt ! 2D state intergration
+   PUBLIC   compute_fabm ! Compute FABM sources, sinks and diagnostics
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: flux    ! Cross-interface flux of pelagic variables (# m-2 s-1)
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   :: ext     ! Light extinction coefficient (m-1)
+
+   ! Work array for mass aggregation
+   REAL(wp), ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:)   :: current_total
+
+
+   ! Arrays for environmental variables
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:) :: prn,rho
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: taubot
+
+   ! repair counters
+   INTEGER :: repair_interior_count,repair_surface_count,repair_bottom_count
+
+   ! state check type
+   TYPE type_state
+      LOGICAL             :: valid
+      LOGICAL             :: repaired
+   END TYPE
+
+   REAL(wp), PUBLIC :: daynumber_in_year
+
+   TYPE (type_bulk_variable_id),SAVE :: swr_id
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE trc_sms_fabm( kt )
+      !!----------------------------------------------------------------------
+      !!                     ***  trc_sms_fabm  ***
+      !!
+      !! ** Purpose :   main routine of FABM model
+      !!
+      !! ** Method  : -
+      !!----------------------------------------------------------------------
+      !
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER :: jn
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrfabm
+
+!!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('trc_sms_fabm')
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,'(a,i0,a,i4.4,a,i2.2,a,i2.2,a,i5,a)') &
+          ' trc_sms_fabm:  FABM model, iteration ',kt,' ', &
+          nyear,'-',nmonth,'-',nday,' ',nsec_day," secs"
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
+
+      CALL update_inputs( kt )
+
+      CALL compute_fabm
+
+      CALL compute_vertical_movement( kt )
+
+      CALL st2d_fabm_nxt( kt )
+
+      IF( l_trdtrc )  CALL wrk_alloc( jpi, jpj, jpk, ztrfabm )
+
+      CALL trc_bc_read  ( kt )       ! tracers: surface and lateral Boundary Conditions
+      CALL trc_rnf_fabm ( kt ) ! River forcings
+
+      IF( l_trdtrc ) THEN      ! Save the trends in the mixed layer
+          DO jn = jp_fabm0, jp_fabm1
+            ztrfabm(:,:,:) = tra(:,:,:,jn)
+            CALL trd_trc( ztrfabm, jn, jptra_sms, kt )   ! save trends
+          END DO
+          CALL wrk_dealloc( jpi, jpj, jpk, ztrfabm )
+      END IF
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('trc_sms_fabm')
+
+   END SUBROUTINE trc_sms_fabm
+
+   SUBROUTINE compute_fabm()
+      INTEGER :: ji,jj,jk,jn
+      TYPE(type_state) :: valid_state
+      REAL(wp) :: zalfg,zztmpx,zztmpy
+
+      ! Validate current model state (setting argument to .TRUE. enables repair=clipping)
+      valid_state = check_state(.TRUE.)
+      IF (.NOT.valid_state%valid) THEN
+         WRITE(numout,*) "Invalid value in FABM encountered in area ",narea,"!!!"
+#if defined key_iomput
+         CALL xios_finalize                ! end mpp communications with xios
+         IF( lk_oasis ) CALL cpl_finalize    ! end coupling and mpp communications with OASIS
+#else
+         IF( lk_oasis ) THEN
+            CALL cpl_finalize              ! end coupling and mpp communications with OASIS
+         ELSE
+            IF( lk_mpp )   CALL mppstop    ! end mpp communications
+         ENDIF
+#endif
+      END IF
+      IF (valid_state%repaired) THEN
+         WRITE(numout,*) "Total interior repairs up to now on process",narea,":",repair_interior_count
+         WRITE(numout,*) "Total surface repairs up to now on process",narea,":",repair_surface_count
+         WRITE(numout,*) "Total bottom repairs up to now on process",narea,":",repair_bottom_count
+      ENDIF
+
+      ! Compute the now hydrostatic pressure
+      ! copied from istate.F90
+      ! ------------------------------------
+
+      zalfg = 0.5e-4 * grav ! FABM wants dbar, convert from Pa
+
+      rho = rau0 * ( 1. + rhd )
+
+      prn(:,:,1) = 10.1325 + zalfg * fse3t(:,:,1) * rho(:,:,1)
+
+      daynumber_in_year=(fjulday-fjulstartyear+1)*1._wp
+
+      DO jk = 2, jpk                                              ! Vertical integration from the surface
+         prn(:,:,jk) = prn(:,:,jk-1) + zalfg * ( &
+                     fse3t(:,:,jk-1) * rho(:,:,jk-1)  &
+                     + fse3t(:,:,jk) * rho(:,:,jk) )
+      END DO
+
+      ! Bottom stress
+      taubot(:,:) = 0._wp
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+               zztmpx = (  bfrua(ji  ,jj) * un(ji  ,jj,mbku(ji  ,jj))  &
+                      &  + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj))  )
+               zztmpy = (  bfrva(ji,  jj) * vn(ji,jj  ,mbkv(ji,jj  ))  &
+                      &  + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1))  )
+               taubot(ji,jj) = 0.5_wp * rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1)
+               !
+         ENDDO
+      ENDDO
+      ! Compute light extinction
+      DO jk=1,jpk
+          DO jj=1,jpj
+            call fabm_get_light_extinction(model,1,jpi,jj,jk,ext)
+         END DO
+      END DO
+
+      ! Compute light field (stored among FABM's internal diagnostics)
+      DO jj=1,jpj
+          DO ji=1,jpi
+            call fabm_get_light(model,1,jpk,ji,jj)
+         END DO
+      END DO
+
+      ! TODO: retrieve 3D shortwave and store in etot3
+
+      ! Zero rate array of interface-attached state variables
+      fabm_st2Da = 0._wp
+
+      ! Compute interfacial source terms and fluxes
+      DO jj=1,jpj
+         ! Process bottom (fabm_do_bottom increments rather than sets, so zero flux array first)
+         flux = 0._wp
+         CALL fabm_do_bottom(model,1,jpi,jj,flux,fabm_st2Da(:,jj,jp_fabm_surface+1:))
+         DO jn=1,jp_fabm
+             DO ji=1,jpi
+                 ! Divide bottom fluxes by height of bottom layer and add to source terms.
+                 ! TODO: is there perhaps an existing variable for fse3t(ji,jj,mbkt(ji,jj))??
+                 tra(ji,jj,mbkt(ji,jj),jp_fabm_m1+jn) = tra(ji,jj,mbkt(ji,jj),jp_fabm_m1+jn) + flux(ji,jn)/fse3t(ji,jj,mbkt(ji,jj))
+             END DO
+         END DO
+
+         ! Process surface (fabm_do_surface increments rather than sets, so zero flux array first)
+         flux = 0._wp
+         CALL fabm_do_surface(model,1,jpi,jj,flux,fabm_st2Da(:,jj,1:jp_fabm_surface))
+         DO jn=1,jp_fabm
+             ! Divide surface fluxes by height of surface layer and add to source terms.
+             tra(:,jj,1,jp_fabm_m1+jn) = tra(:,jj,1,jp_fabm_m1+jn) + flux(:,jn)/fse3t(:,jj,1)
+         END DO
+      END DO
+
+      ! Compute interior source terms (NB fabm_do increments rather than sets)
+      DO jk=1,jpk
+          DO jj=1,jpj
+              CALL fabm_do(model,1,jpi,jj,jk,tra(:,jj,jk,jp_fabm0:jp_fabm1))
+          END DO
+      END DO
+   END SUBROUTINE compute_fabm
+
+   FUNCTION check_state(repair) RESULT(exit_state)
+      LOGICAL, INTENT(IN) :: repair
+      TYPE(type_state) :: exit_state
+
+      INTEGER             :: jj,jk
+      LOGICAL             :: valid_int,valid_sf,valid_bt
+
+      exit_state%valid = .TRUE.
+      exit_state%repaired =.FALSE.
+      DO jk=1,jpk
+         DO jj=1,jpj
+            CALL fabm_check_state(model,1,jpi,jj,jk,repair,valid_int)
+            IF (repair.AND..NOT.valid_int) THEN
+               repair_interior_count = repair_interior_count + 1
+               exit_state%repaired = .TRUE.
+            END IF
+            IF (.NOT.(valid_int.OR.repair)) exit_state%valid = .FALSE.
+         END DO
+      END DO
+      DO jj=1,jpj
+         CALL fabm_check_surface_state(model,1,jpi,jj,repair,valid_sf)
+         IF (repair.AND..NOT.valid_sf) THEN
+            repair_surface_count = repair_surface_count + 1
+            exit_state%repaired = .TRUE.
+         END IF
+         IF (.NOT.(valid_sf.AND.valid_bt).AND..NOT.repair) exit_state%valid = .FALSE.
+         CALL fabm_check_bottom_state(model,1,jpi,jj,repair,valid_bt)
+         IF (repair.AND..NOT.valid_bt) THEN
+            repair_bottom_count = repair_bottom_count + 1
+            exit_state%repaired = .TRUE.
+         END IF
+         IF (.NOT.(valid_sf.AND.valid_bt).AND..NOT.repair) exit_state%valid = .FALSE.
+      END DO
+   END FUNCTION
+
+   SUBROUTINE trc_sms_fabm_check_mass()
+      REAL(wp) :: total(SIZE(model%conserved_quantities))
+      INTEGER :: jk,jj,jn
+
+      total = 0._wp
+
+      DO jk=1,jpk
+         DO jj=1,jpj
+            CALL fabm_get_conserved_quantities(model,1,jpi,jj,jk,current_total)
+            DO jn=1,SIZE(model%conserved_quantities)
+               total(jn) = total(jn) + SUM(cvol(:,jj,jk)*current_total(:,jn)*tmask_i(:,jj))
+            END DO
+         END DO
+      END DO
+
+      DO jj=1,jpj
+         CALL fabm_get_horizontal_conserved_quantities(model,1,jpi,jj,current_total)
+         DO jn=1,SIZE(model%conserved_quantities)
+            total(jn) = total(jn) + SUM(e1e2t(:,jj)*current_total(:,jn)*tmask_i(:,jj))
+         END DO
+      END DO
+
+      IF( lk_mpp ) CALL mpp_sum(total,SIZE(model%conserved_quantities))
+
+      DO jn=1,SIZE(model%conserved_quantities)
+         IF(lwp) WRITE(numout,*) 'FABM '//TRIM(model%conserved_quantities(jn)%name),total(jn),TRIM(model%conserved_quantities(jn)%units)//'*m3'
+      END DO
+
+   END SUBROUTINE trc_sms_fabm_check_mass
+
+   SUBROUTINE st2d_fabm_nxt( kt )
+      !!----------------------------------------------------------------------
+      !!                     ***  st2d_fabm_nxt  ***
+      !!
+      !! ** Purpose :   routine to integrate 2d states in time
+      !!
+      !! ** Method  :   based on integration of 3D passive tracer fields
+      !!                implemented in TOP_SRC/TRP/trcnxt.F90, plus
+      !!                tra_nxt_fix in OPA_SRC/TRA/tranxt.F90. Similar to
+      !!                time integration of sea surface height in
+      !!                OPA_SRC/DYN/sshwzv.F90.
+      !!----------------------------------------------------------------------
+      !
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      REAL(wp) :: z2dt
+      INTEGER :: jn
+
+!!----------------------------------------------------------------------
+      !
+      IF( neuler == 0 .AND. kt == nittrc000 ) THEN
+          z2dt = rdt                  ! set time step size (Euler)
+      ELSE
+          z2dt = 2._wp * rdt          ! set time step size (Leapfrog)
+      ENDIF
+
+      ! Forward Euler time step to compute "now"
+      DO jn=1,jp_fabm_surface+jp_fabm_bottom
+         fabm_st2Da(:,:,jn) = (fabm_st2db(:,:,jn) + z2dt * fabm_st2da(:,:,jn)) * tmask(:,:,1)
+      ENDDO
+
+      IF( neuler == 0 .AND. kt == nittrc000 ) THEN        ! Euler time-stepping at first time-step
+         !                                                ! (only swap)
+         fabm_st2Dn(:,:,:) = fabm_st2Da(:,:,:)
+         !
+      ELSE
+         ! Update now state + Asselin filter time stepping
+         fabm_st2Db(:,:,:) = (1._wp - 2._wp*atfp) * fabm_st2Dn(:,:,:) + &
+             atfp * ( fabm_st2Db(:,:,:) + fabm_st2Da(:,:,:) )
+         fabm_st2Dn(:,:,:) = fabm_st2Da(:,:,:)
+      ENDIF
+
+   END SUBROUTINE st2d_fabm_nxt
+
+   INTEGER FUNCTION trc_sms_fabm_alloc()
+      INTEGER :: jj,jk,jn
+      !!----------------------------------------------------------------------
+      !!              ***  ROUTINE trc_sms_fabm_alloc  ***
+      !!----------------------------------------------------------------------
+      !
+      ! ALLOCATE here the arrays specific to FABM
+      ALLOCATE( lk_rad_fabm(jp_fabm))
+      ALLOCATE( prn(jpi, jpj, jpk))
+      ALLOCATE( rho(jpi, jpj, jpk))
+      ALLOCATE( taubot(jpi, jpj))
+      ! ALLOCATE( tab(...) , STAT=trc_sms_fabm_alloc )
+
+      ! Allocate arrays to hold state for surface-attached and bottom-attached state variables
+      ALLOCATE(fabm_st2Dn(jpi, jpj, jp_fabm_surface+jp_fabm_bottom))
+      ALLOCATE(fabm_st2Da(jpi, jpj, jp_fabm_surface+jp_fabm_bottom))
+      ALLOCATE(fabm_st2Db(jpi, jpj, jp_fabm_surface+jp_fabm_bottom))
+
+      ! Work array to hold surface and bottom fluxes
+      ALLOCATE(flux(jpi,jp_fabm))
+
+      ! Work array to hold extinction coefficients
+      ALLOCATE(ext(jpi))
+      ext=0._wp
+
+      ! Allocate work arrays for vertical movement
+      ALLOCATE(w_ct(jpi,jpk,jp_fabm))
+      ALLOCATE(w_if(jpk,jp_fabm))
+      ALLOCATE(zwgt_if(jpk,jp_fabm))
+      ALLOCATE(flux_if(jpk,jp_fabm))
+      ALLOCATE(current_total(jpi,SIZE(model%conserved_quantities)))
+#if defined key_trdtrc && defined key_iomput
+      IF( lk_trdtrc ) ALLOCATE(tr_vmv(jpi,jpj,jpk,jp_fabm))
+      IF( lk_trdtrc ) ALLOCATE(tr_inp(jpi,jpj,jpk))
+#endif
+
+      trc_sms_fabm_alloc = 0      ! set to zero if no array to be allocated
+      !
+      IF( trc_sms_fabm_alloc /= 0 ) CALL ctl_warn('trc_sms_fabm_alloc : failed to allocate arrays')
+      !
+
+      ! Make FABM aware of diagnostics that are not needed [not included in output]
+      DO jn=1,size(model%diagnostic_variables)
+          !model%diagnostic_variables(jn)%save = iom_use(model%diagnostic_variables(jn)%name)
+      END DO
+      DO jn=1,size(model%horizontal_diagnostic_variables)
+          !model%horizontal_diagnostic_variables(jn)%save = iom_use(model%horizontal_diagnostic_variables(jn)%name)
+      END DO
+
+      ! Provide FABM with domain extents [after this, the save attribute of diagnostic variables can no longe change!]
+      call fabm_set_domain(model,jpi, jpj, jpk)
+
+      ! Provide FABM with the vertical indices of the surface and bottom, and the land-sea mask.
+      call model%set_bottom_index(mbkt)  ! NB mbkt extents should match dimension lengths provided to fabm_set_domain
+      call model%set_surface_index(1)
+      call fabm_set_mask(model,tmask,tmask(:,:,1)) ! NB tmask extents should match dimension lengths provided to fabm_set_domain
+
+      ! Send pointers to state data to FABM
+      do jn=1,jp_fabm
+         call fabm_link_bulk_state_data(model,jn,trn(:,:,:,jp_fabm_m1+jn))
+      end do
+      DO jn=1,jp_fabm_surface
+         CALL fabm_link_surface_state_data(model,jn,fabm_st2Dn(:,:,jn))
+      END DO
+      DO jn=1,jp_fabm_bottom
+         CALL fabm_link_bottom_state_data(model,jn,fabm_st2Dn(:,:,jp_fabm_surface+jn))
+      END DO
+
+      ! Send pointers to environmental data to FABM
+      call fabm_link_bulk_data(model,standard_variables%temperature,tsn(:,:,:,jp_tem))
+      call fabm_link_bulk_data(model,standard_variables%practical_salinity,tsn(:,:,:,jp_sal))
+      call fabm_link_bulk_data(model,standard_variables%density,rho(:,:,:))
+      call fabm_link_bulk_data(model,standard_variables%pressure,prn)
+      call fabm_link_horizontal_data(model,standard_variables%bottom_stress,taubot(:,:))
+      ! correct target for cell thickness depends on NEMO configuration:
+#ifdef key_vvl
+      call fabm_link_bulk_data(model,standard_variables%cell_thickness,e3t_n)
+#else
+      call fabm_link_bulk_data(model,standard_variables%cell_thickness,e3t_0)
+#endif
+      call fabm_link_horizontal_data(model,standard_variables%latitude,gphit)
+      call fabm_link_horizontal_data(model,standard_variables%longitude,glamt)
+      call fabm_link_scalar_data(model,standard_variables%number_of_days_since_start_of_the_year,daynumber_in_year)
+      call fabm_link_horizontal_data(model,standard_variables%wind_speed,wndm(:,:))
+      call fabm_link_horizontal_data(model,standard_variables%surface_downwelling_shortwave_flux,qsr(:,:))
+      call fabm_link_horizontal_data(model,standard_variables%bottom_depth_below_geoid,bathy(:,:))
+
+      swr_id = model%get_bulk_variable_id(standard_variables%downwelling_shortwave_flux)
+
+      ! Obtain user-specified input variables (read from NetCDF file)
+      call link_inputs
+      call update_inputs( nit000, .false. )
+
+      ! Check whether FABM has all required data
+      call fabm_check_ready(model)
+
+      ! Initialize state
+      DO jj=1,jpj
+         CALL fabm_initialize_surface_state(model,1,jpi,jj)
+         CALL fabm_initialize_bottom_state(model,1,jpi,jj)
+      END DO
+      DO jk=1,jpk
+         DO jj=1,jpj
+            CALL fabm_initialize_state(model,1,jpi,jj,jk)
+         END DO
+      END DO
+
+      ! Set mask for negativity corrections to the relevant states
+      lk_rad_fabm = .FALSE.
+      DO jn=1,jp_fabm
+        IF (model%state_variables(jn)%minimum.ge.0) THEN
+          lk_rad_fabm(jn)=.TRUE.
+          IF(lwp) WRITE(numout,*) 'FABM clipping for '//TRIM(model%state_variables(jn)%name)//' activated.'
+        END IF
+      END DO
+
+      ! Mask land points in states with zeros, not nice, but coherent
+      ! with NEMO "convention":
+      DO jn=jp_fabm0,jp_fabm1
+        WHERE (tmask==0._wp)
+          trn(:,:,:,jn)=0._wp
+        END WHERE
+      END DO
+      DO jn=1,jp_fabm_surface+jp_fabm_bottom
+        WHERE (tmask(:,:,1)==0._wp)
+          fabm_st2Dn(:,:,jn)=0._wp
+        END WHERE
+      END DO
+
+      ! Copy initial condition for interface-attached state variables to "previous" state field
+      ! NB NEMO does this itself for pelagic state variables (trb) in TOP_SRC/trcini.F90.
+      fabm_st2Db = fabm_st2Dn
+
+      ! Initialise repair counters
+      repair_interior_count = 0
+      repair_surface_count = 0
+      repair_bottom_count = 0
+
+   END FUNCTION trc_sms_fabm_alloc
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Dummy module                                        No FABM model
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_sms_fabm( kt )             ! Empty routine
+      INTEGER, INTENT( in ) ::   kt
+      WRITE(*,*) 'trc_sms_fabm: You should not have seen this print! error?', kt
+   END SUBROUTINE trc_sms_fabm
+#endif
+
+   !!======================================================================
+END MODULE trcsms_fabm
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/trcwri_fabm.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/trcwri_fabm.F90	(revision 10156)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/trcwri_fabm.F90	(revision 10156)
@@ -0,0 +1,134 @@
+MODULE trcwri_fabm
+   !!======================================================================
+   !!                       *** MODULE trcwri_fabm ***
+   !!    fabm :   Output of FABM tracers
+   !!======================================================================
+   !! History :   1.0  !  2009-05 (C. Ethe)  Original code
+   !!----------------------------------------------------------------------
+#if defined key_top && key_fabm && defined key_iomput
+   !!----------------------------------------------------------------------
+   !!   'key_fabm'                                           FABM model
+   !!----------------------------------------------------------------------
+   !! trc_wri_fabm   :  outputs of concentration fields
+   !!----------------------------------------------------------------------
+   USE trc         ! passive tracers common variables 
+   USE iom         ! I/O manager
+   USE trdtrc_oce
+   USE trcsms_fabm, only: trc_sms_fabm_check_mass
+   USE par_fabm
+   USE st2d_fabm
+   USE fabm, only: fabm_get_bulk_diagnostic_data, fabm_get_horizontal_diagnostic_data
+
+   IMPLICIT NONE
+   PRIVATE
+
+#if defined key_tracer_budget
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: tr_temp
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: fabm_st2d_temp
+#endif
+
+   INTERFACE trc_wri_fabm
+       MODULE PROCEDURE wri_fabm,wri_fabm_fl
+   END INTERFACE trc_wri_fabm
+
+
+   PUBLIC trc_wri_fabm 
+
+#  include "top_substitute.h90"
+CONTAINS
+
+   SUBROUTINE wri_fabm_fl(kt,fl)
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_wri_trc  ***
+      !!
+      !! ** Purpose :   output passive tracers fields 
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT( in )               :: fl
+      INTEGER, INTENT( in )               :: kt
+
+#if defined key_tracer_budget
+      INTEGER              :: jn
+      CHARACTER (len=20)   :: cltra
+      REAL(wp), DIMENSION(jpi,jpj,jpk)    :: trpool !temporary storage pool 3D
+      REAL(wp), DIMENSION(jpi,jpj)    :: st2dpool !temporary storage pool 2D
+      !!---------------------------------------------------------------------
+ 
+      ! write the tracer concentrations in the file
+      ! ---------------------------------------
+! depth integrated
+! for strict budgetting write this out at end of timestep as an average between 'now' and 'after' at kt
+      DO jn = 1, jp_fabm1
+        IF(ln_trdtrc (jn))THEN
+         trpool(:,:,:) = 0.5 * ( trn(:,:,:,jp_fabm0+jn-1)*fse3t_a(:,:,:) + &
+                             tr_temp(:,:,:,jn)*fse3t(:,:,:) )
+         cltra = TRIM( model%state_variables(jn)%name )//"_e3t"     ! depth integrated output
+         IF( kt == nittrc000 ) write(6,*)'output pool ',cltra
+         CALL iom_put( cltra, trpool)
+        ENDIF
+      END DO
+#else
+      CONTINUE
+#endif
+
+   END SUBROUTINE wri_fabm_fl
+
+   SUBROUTINE wri_fabm(kt)
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_wri_trc  ***
+      !!
+      !! ** Purpose :   output passive tracers fields 
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT( in )               :: kt
+      INTEGER              :: jn
+
+#if defined key_tracer_budget
+      IF( kt == nittrc000 ) THEN
+         ALLOCATE(tr_temp(jpi,jpj,jpk,jp_fabm),fabm_st2d_temp(jpi,jpj,jp_fabm_surface+jp_fabm_bottom))
+      ENDIF
+      tr_temp(:,:,:,:)=trn(:,:,:,jp_fabm0:jp_fabm1) ! slwa save for tracer budget (unfiltered trn)
+      fabm_st2d_temp(:,:,:)=fabm_st2dn(:,:,:)
+#endif
+      DO jn = 1, jp_fabm
+         CALL iom_put( model%state_variables(jn)%name, trn(:,:,:,jp_fabm0+jn-1) )
+      END DO
+      DO jn = 1, jp_fabm_surface
+         CALL iom_put( model%surface_state_variables(jn)%name, fabm_st2dn(:,:,jn) )
+      END DO
+      DO jn = 1, jp_fabm_bottom
+         CALL iom_put( model%bottom_state_variables(jn)%name, fabm_st2dn(:,:,jp_fabm_surface+jn) )
+      END DO
+
+      ! write 3D diagnostics in the file
+      ! ---------------------------------------
+      DO jn = 1, size(model%diagnostic_variables)
+         IF (model%diagnostic_variables(jn)%save) &
+             CALL iom_put( model%diagnostic_variables(jn)%name, fabm_get_bulk_diagnostic_data(model,jn))
+      END DO
+
+      ! write 2D diagnostics in the file
+      ! ---------------------------------------
+      DO jn = 1, size(model%horizontal_diagnostic_variables)
+         IF (model%horizontal_diagnostic_variables(jn)%save) &
+             CALL iom_put( model%horizontal_diagnostic_variables(jn)%name, fabm_get_horizontal_diagnostic_data(model,jn))
+      END DO
+      !
+      CALL trc_sms_fabm_check_mass
+
+   END SUBROUTINE wri_fabm
+
+#else
+   !!----------------------------------------------------------------------
+   !!  Dummy module :                                     No passive tracer
+   !!----------------------------------------------------------------------
+   PUBLIC trc_wri_fabm
+CONTAINS
+   SUBROUTINE trc_wri_fabm                     ! Empty routine  
+   END SUBROUTINE trc_wri_fabm
+#endif
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
+   !! $Id: trcwri_fabm.F90 3160 2011-11-20 14:27:18Z cetlod $ 
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!======================================================================
+END MODULE trcwri_fabm
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90	(revision 10156)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90	(revision 10156)
@@ -0,0 +1,223 @@
+MODULE vertical_movement_fabm
+   !!======================================================================
+   !!                         ***  MODULE vertical_movement_fabm  ***
+   !! TOP :   Module for the vertical movement of the FABM tracers
+   !!======================================================================
+
+#if defined key_fabm
+   !!----------------------------------------------------------------------
+   !!   'key_fabm'                                               FABM tracers
+   !!----------------------------------------------------------------------
+   !! compute_vertical_movement : compute vertical movement of FABM fields
+   !!----------------------------------------------------------------------
+   USE par_trc
+   USE oce_trc
+   USE trc
+   USE par_fabm
+   USE fabm
+   USE dom_oce
+#if defined key_trdtrc && defined key_iomput
+   USE iom
+   USE trdtrc_oce
+#endif
+
+   IMPLICIT NONE
+
+#  include "domzgr_substitute.h90"
+
+   PRIVATE
+
+   PUBLIC compute_vertical_movement
+
+   ! Work arrays for vertical advection (residual movement/sinking/floating)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:) :: w_ct
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:)   :: w_if
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:)   :: zwgt_if
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:)   :: flux_if
+#if defined key_trdtrc && defined key_iomput
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:,:) :: tr_vmv
+#endif
+
+   CONTAINS
+
+   SUBROUTINE compute_vertical_movement( kt )
+      !!----------------------------------------------------------------------
+      !!                     ***  compute_vertical_movement  ***
+      !!
+      !! ** Purpose :   compute vertical movement of FABM tracers
+      !!
+      !! ** Method  : Sets additional vertical velocity field and computes
+      !!              resulting advection using a conservative 3rd upwind
+      !!              scheme with QUICKEST TVD limiter, based on GOTM
+      !!              module adv_center.F90 (www.gotm.net). Currently assuming
+      !!              zero flux at sea surface and sea floor.
+      !!----------------------------------------------------------------------
+      !
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER :: ji,jj,jk,jn,k_floor,n_iter,n_count
+      INTEGER,PARAMETER :: n_itermax=100
+      REAL(wp) :: cmax_no,z2dt
+      REAL(wp),DIMENSION(jpk) :: tr_it,tr_u,tr_d,tr_c,tr_slope,c_no,flux_lim
+      REAL(wp),DIMENSION(jpk) :: phi_lim,x_fac
+#if defined key_trdtrc
+      CHARACTER (len=20) :: cltra
+#endif
+
+#if defined key_trdtrc && defined key_iomput
+      IF( lk_trdtrc ) tr_vmv = 0.0_wp
+#endif
+
+      IF( neuler == 0 .AND. kt == nittrc000 ) THEN
+          z2dt = rdt                  ! set time step size (Euler)
+      ELSE
+          z2dt = 2._wp * rdt          ! set time step size (Leapfrog)
+      ENDIF
+      ! Compute interior vertical velocities and include them in source array.
+      DO jj=1,jpj ! j-loop
+         ! Get vertical velocities at layer centres (entire 1:jpi,1:jpk slice).
+         DO jk=1,jpk
+            CALL fabm_get_vertical_movement(model,1,jpi,jj,jk,w_ct(:,jk,:))
+         END DO
+
+         DO ji=1,jpi ! i-loop
+            ! Only process this horizontal point (ji,jj) if number of layers exceeds 1
+            IF (mbkt(ji,jj)>1) THEN ! Level check
+               k_floor=mbkt(ji,jj)
+               ! Linearly interpolate to velocities at the interfaces between layers
+               ! Note:
+               !    - interface k sits between cell centre k and k-1,
+               !    - k [1,jpk] increases downwards
+               !    - upward velocity is positive, downward velocity is negative
+               zwgt_if(1,:)=0._wp ! surface
+               w_if(1,:)=0._wp ! surface
+               zwgt_if(2:k_floor,:)=spread(&
+                   fse3t(ji,jj,2:k_floor)/ (fse3t(ji,jj,1:k_floor-1)+fse3t(ji,jj,2:k_floor))&
+                   ,2,jp_fabm)
+               w_if(2:k_floor,:) = zwgt_if(2:k_floor,:)*w_ct(ji,1:k_floor-1,:)&
+                  +(1._wp-zwgt_if(1:k_floor-1,:))*w_ct(ji,2:k_floor,:)
+               zwgt_if(k_floor+1:,:)=0._wp ! sea floor and below
+               w_if(k_floor+1:,:)=0._wp ! sea floor and below
+
+               ! Advect:
+               DO jn=1,jp_fabm ! State loop
+                  ! get maximum Courant number:
+                  c_no(2:k_floor)=abs(w_if(2:k_floor,jn))*z2dt/ &
+                                ( 0.5_wp*(fse3t(ji,jj,2:k_floor) + &
+                                fse3t(ji,ji,1:k_floor-1)) )
+                  cmax_no=MAXVAL(c_no(2:k_floor))
+
+                  ! number of iterations:
+                  n_iter=min(n_itermax,int(cmax_no)+1)
+                  IF (ln_ctl.AND.(n_iter .gt. 1)) THEN
+                      WRITE(numout,*) 'vertical_movement_fabm():'
+                      WRITE(numout,*) '   Maximum Courant number is ',cmax_no,'.'
+                      WRITE(numout,*) '   ',n_iter,' iterations used for vertical advection.'
+                  ENDIF
+
+                  ! effective Courant number:
+                  c_no=c_no/n_iter
+
+                  tr_it(1:k_floor)=trb(ji,jj,1:k_floor,jp_fabm_m1+jn)
+                  DO n_count=1,n_iter ! Iterative loop
+                     !Compute slope ratio
+                     IF (k_floor.gt.2) THEN !More than 2 vertical wet points
+                        IF (k_floor.gt.3) THEN
+                          WHERE (w_if(3:k_floor-1,jn).ge.0._wp) !upward movement
+                           tr_u(3:k_floor-1)=tr_it(4:k_floor)
+                           tr_c(3:k_floor-1)=tr_it(3:k_floor-1)
+                           tr_d(3:k_floor-1)=tr_it(2:k_floor-2)
+                          ELSEWHERE !downward movement
+                           tr_u(3:k_floor-1)=tr_it(1:k_floor-3)
+                           tr_c(3:k_floor-1)=tr_it(2:k_floor-2)
+                           tr_d(3:k_floor-1)=tr_it(3:k_floor-1)
+                          ENDWHERE
+                        ENDIF
+                        IF (w_if(2,jn).ge.0._wp) THEN
+                           tr_u(2)=tr_it(3)
+                           tr_c(2)=tr_it(2)
+                           tr_d(2)=tr_it(1)
+                        ELSE
+                           tr_u(2)=tr_it(1)
+                           tr_c(2)=tr_it(1)
+                           tr_d(2)=tr_it(2)
+                        ENDIF
+                        IF (w_if(k_floor,jn).ge.0._wp) THEN
+                           tr_u(k_floor)=tr_it(k_floor)
+                           tr_c(k_floor)=tr_it(k_floor)
+                           tr_d(k_floor)=tr_it(k_floor-1)
+                        ELSE
+                           tr_u(k_floor)=tr_it(k_floor-2)
+                           tr_c(k_floor)=tr_it(k_floor-1)
+                           tr_d(k_floor)=tr_it(k_floor)
+                        ENDIF
+                     ELSE !only 2 vertical wet points, i.e. only 1 interface
+                        IF (w_if(k_floor,jn).ge.0._wp) THEN
+                           tr_u(2)=tr_it(2)
+                           tr_c(2)=tr_it(2)
+                           tr_d(2)=tr_it(1)
+                        ELSE
+                           tr_u(2)=tr_it(1)
+                           tr_c(2)=tr_it(1)
+                           tr_d(2)=tr_it(2)
+                        ENDIF
+                     ENDIF
+                     WHERE (abs(tr_d(2:k_floor)-tr_c(2:k_floor)).gt.1.e-10_wp)
+                        tr_slope(2:k_floor)= &
+                           (tr_c(2:k_floor)-tr_u(2:k_floor))/ &
+                           (tr_d(2:k_floor)-tr_c(2:k_floor))
+                     ELSEWHERE
+                        tr_slope(2:k_floor)=SIGN(1._wp,w_if(2:k_floor,jn))* &
+                              (tr_c(2:k_floor)-tr_u(2:k_floor))*1.e10_wp
+                     ENDWHERE
+
+                     !QUICKEST flux limiter:
+                     x_fac(2:k_floor)=(1._wp-2._wp*c_no(2:k_floor))/6._wp
+                     phi_lim(2:k_floor)=(0.5_wp+x_fac(2:k_floor)) + &
+                        (0.5_wp-x_Fac(2:k_floor))*tr_slope(2:k_floor)
+                     flux_lim(2:k_floor)=max( 0._wp, &
+                       min( phi_lim(2:k_floor),2._wp/(1._wp-c_no(2:k_floor)), &
+                         2._wp*tr_slope(2:k_floor)/(c_no(2:k_floor)+1.e-10_wp)) )
+
+                     ! Compute limited flux:
+                     flux_if(2:k_floor,jn) = w_if(2:k_floor,jn)* &
+                        ( tr_c(2:k_floor) + &
+                        0.5_wp*flux_lim(2:k_floor)*(1._wp-c_no(2:k_floor))* &
+                        (tr_d(2:k_floor)-tr_c(2:k_floor)) )
+
+                     ! Compute pseudo update for trend aggregation:
+                     tr_it(1:k_floor-1) = tr_it(1:k_floor-1) + &
+                        z2dt/float(n_iter)/fse3t(ji,jj,1:k_floor-1)* &
+                        flux_if(2:k_floor,jn)
+                     tr_it(2:k_floor) = tr_it(2:k_floor) - &
+                        z2dt/float(n_iter)/fse3t(ji,jj,2:k_floor)* &
+                        flux_if(2:k_floor,jn)
+
+                  ENDDO ! Iterative loop
+
+                  ! Estimate rate of change from pseudo state updates (source
+                  ! splitting):
+                  tra(ji,jj,1:k_floor,jp_fabm_m1+jn) = &
+                     tra(ji,jj,1:k_floor,jp_fabm_m1+jn) + &
+                     (tr_it(1:k_floor) - trb(ji,jj,1:k_floor,jp_fabm_m1+jn))/z2dt
+#if defined key_trdtrc && defined key_iomput
+                  IF( lk_trdtrc .AND. ln_trdtrc( jp_fabm_m1+jn ) ) THEN
+                    tr_vmv(ji,jj,1:k_floor,jn)=(tr_it(1:k_floor) - trb(ji,jj,1:k_floor,jn))/z2dt
+                  END IF
+#endif
+               ENDDO ! State loop
+            END IF ! Level check
+         END DO ! i-loop
+      END DO ! j-loop
+#if defined key_trdtrc && defined key_iomput
+      DO jn=1,jp_fabm ! State loop
+        IF( lk_trdtrc .AND. ln_trdtrc(jp_fabm_m1+jn) ) THEN
+          cltra = 'VMV_'//TRIM(ctrcnm(jp_fabm_m1+jn))
+          CALL iom_put( cltra,  tr_vmv(:,:,:,jn) )
+        END IF
+      ENDDO
+#endif
+
+   END SUBROUTINE compute_vertical_movement
+
+#endif
+END MODULE
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/MY_TRC/trcsms_my_trc.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/MY_TRC/trcsms_my_trc.F90	(revision 10155)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/MY_TRC/trcsms_my_trc.F90	(revision 10156)
@@ -18,4 +18,5 @@
    USE trd_oce
    USE trdtrc
+   USE trcbc, only : trc_bc_read
 
    IMPLICIT NONE
@@ -55,4 +56,6 @@
 
       IF( l_trdtrc )  CALL wrk_alloc( jpi, jpj, jpk, ztrmyt )
+
+      CALL trc_bc_read  ( kt )       ! tracers: surface and lateral Boundary Conditions
 
       IF( l_trdtrc ) THEN      ! Save the trends in the ixed layer
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/MY_TRC/trcwri_my_trc.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/MY_TRC/trcwri_my_trc.F90	(revision 10155)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/MY_TRC/trcwri_my_trc.F90	(revision 10156)
@@ -19,9 +19,17 @@
 
    PUBLIC trc_wri_my_trc 
+#if defined key_tracer_budget
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: trb_temp ! slwa
+#endif
+
 
 #  include "top_substitute.h90"
 CONTAINS
 
+#if defined key_tracer_budget
+   SUBROUTINE trc_wri_my_trc (kt, fl) ! slwa
+#else
    SUBROUTINE trc_wri_my_trc
+#endif
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE trc_wri_trc  ***
@@ -29,14 +37,52 @@
       !! ** Purpose :   output passive tracers fields 
       !!---------------------------------------------------------------------
+#if defined key_tracer_budget
+      INTEGER, INTENT( in ), OPTIONAL     :: fl 
+      INTEGER, INTENT( in )               :: kt
+      REAL(wp), DIMENSION(jpi,jpj,jpk)    :: trpool !tracer pool temporary output
+#else
+      INTEGER, INTENT( in )               :: kt
+#endif
       CHARACTER (len=20)   :: cltra
-      INTEGER              :: jn
+      INTEGER              :: jn,jk ! JC TODO jk defined here but may not be used
       !!---------------------------------------------------------------------
  
       ! write the tracer concentrations in the file
       ! ---------------------------------------
+
+
+#if defined key_tracer_budget
+      IF( PRESENT(fl)) THEN
+! depth integrated
+! for strict budgetting write this out at end of timestep as an average between 'now' and 'after' at kt
+         DO jn = jp_myt0, jp_myt1 
+          IF(ln_trdtrc (jn))THEN
+            trpool(:,:,:) = 0.5 * ( trn(:,:,:,jn) * fse3t_a(:,:,:) +  &
+                                        trb_temp(:,:,:,jn) * fse3t(:,:,:) )
+            cltra = TRIM( ctrcnm(jn) )//"e3t"     ! depth integrated output
+            IF( kt == nittrc000 ) write(6,*)'output pool ',cltra
+            DO jk = 1, jpk
+               trpool(:,:,jk) = trpool(:,:,jk)
+            END DO
+            CALL iom_put( cltra, trpool)
+
+          END IF
+         END DO
+
+      ELSE
+
+         IF( kt == nittrc000 ) THEN
+           ALLOCATE(trb_temp(jpi,jpj,jpk,jp_my_trc))  ! slwa
+         ENDIF
+         trb_temp(:,:,:,:)=trn(:,:,:,:) ! slwa save for tracer budget (unfiltered trn)
+
+
+      END IF
+#else
       DO jn = jp_myt0, jp_myt1
          cltra = TRIM( ctrcnm(jn) )                  ! short title for tracer
          CALL iom_put( cltra, trn(:,:,:,jn) )
       END DO
+#endif
       !
    END SUBROUTINE trc_wri_my_trc
@@ -48,5 +94,13 @@
    PUBLIC trc_wri_my_trc
 CONTAINS
-   SUBROUTINE trc_wri_my_trc                     ! Empty routine  
+#if defined key_tracer_budget
+   SUBROUTINE trc_wri_my_trc (kt, fl) ! slwa
+      INTEGER, INTENT( in ), OPTIONAL     :: fl 
+      INTEGER, INTENT( in )               :: kt
+#else
+   ! JC TODO Subroutine arguments (kt) inconsistent with earlier definition
+   SUBROUTINE trc_wri_my_trc (kt)
+      INTEGER, INTENT( in )               :: kt
+#endif
    END SUBROUTINE trc_wri_my_trc
 #endif
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90	(revision 10155)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90	(revision 10156)
@@ -56,5 +56,5 @@
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
       !!
-      INTEGER            :: jn
+      INTEGER            :: jn, jk
       CHARACTER (len=22) :: charout
       REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztrtrd
@@ -105,4 +105,9 @@
         DO jn = 1, jptra
            ztrtrd(:,:,:,jn) = tra(:,:,:,jn) - ztrtrd(:,:,:,jn)
+#if defined key_tracer_budget
+           DO jk = 1, jpkm1
+             ztrtrd(:,:,jk,jn) = ztrtrd(:,:,jk,jn) * e1t(:,:) * e2t(:,:) * fse3t(:,:,jk)  ! slwa
+           END DO
+#endif
            CALL trd_tra( kt, 'TRC', jn, jptra_ldf, ztrtrd(:,:,:,jn) )
         END DO
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/TRP/trcnxt.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/TRP/trcnxt.F90	(revision 10155)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/TRP/trcnxt.F90	(revision 10156)
@@ -33,4 +33,7 @@
    USE trdtra
    USE tranxt
+   USE trcbdy          ! BDY open boundaries
+   USE bdy_par, only: lk_bdy
+   USE iom
 # if defined key_agrif
    USE agrif_top_interp
@@ -93,4 +96,7 @@
       CHARACTER (len=22) :: charout
       REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  ztrdt 
+#if defined key_tracer_budget
+      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  ztrdt_m1 ! slwa
+#endif
       !!----------------------------------------------------------------------
       !
@@ -101,4 +107,19 @@
          WRITE(numout,*) 'trc_nxt : time stepping on passive tracers'
       ENDIF
+#if defined key_tracer_budget
+      IF( kt == nittrc000 .AND. l_trdtrc ) THEN
+         ALLOCATE( ztrdt_m1(jpi,jpj,jpk,jptra) )  ! slwa
+         IF( ln_rsttr .AND.    &                     ! Restart: read in restart  file
+            iom_varid( numrtr, 'atf_trend_'//TRIM(ctrcnm(1)), ldstop = .FALSE. ) > 0 ) THEN
+            IF(lwp) WRITE(numout,*) '          nittrc000-nn_dttrc ATF tracer trend read in the restart file'
+            DO jn = 1, jptra
+               CALL iom_get( numrtr, jpdom_autoglo, 'atf_trend_'//TRIM(ctrcnm(jn)), ztrdt_m1(:,:,:,jn) )   ! before tracer trend for atf
+            END DO
+         ELSE          
+           ztrdt_m1=0.0
+         ENDIF
+      ENDIF
+#endif
+
 
 #if defined key_agrif
@@ -111,7 +132,5 @@
 
 
-#if defined key_bdy
-!!      CALL bdy_trc( kt )               ! BDY open boundaries
-#endif
+      IF( lk_bdy )  CALL trc_bdy( kt )               ! BDY open boundaries
 
 
@@ -149,9 +168,37 @@
                zfact = 1.e0 / r2dt(jk)  
                ztrdt(:,:,jk,jn) = ( trb(:,:,jk,jn) - ztrdt(:,:,jk,jn) ) * zfact 
-               CALL trd_tra( kt, 'TRC', jn, jptra_atf, ztrdt )
+!slwa          CALL trd_tra( kt, 'TRC', jn, jptra_atf, ztrdt )
+#if defined key_tracer_budget
+               ztrdt(:,:,jk,jn) = ztrdt(:,:,jk,jn) * e1t(:,:) * e2t(:,:) * e3t_n(:,:,jk)  ! slwa vvl
+               !ztrdt(:,:,jk,jn) = ztrdt(:,:,jk,jn) * e1t(:,:) * e2t(:,:) * e3t_0(:,:,jk)  ! slwa CHANGE for vvl
+#endif
             END DO
+#if defined key_tracer_budget
+! slwa budget code
+              CALL trd_tra( kt, 'TRC', jn, jptra_atf, ztrdt_m1(:,:,:,jn) )
+#else
+              CALL trd_tra( kt, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
+#endif
          END DO
+#if defined key_tracer_budget
+        ztrdt_m1(:,:,:,:) = ztrdt(:,:,:,:)    ! need previous time step for budget slwa
+#endif
          CALL wrk_dealloc( jpi, jpj, jpk, jptra, ztrdt ) 
       END IF
+
+#if defined key_tracer_budget
+      !                                           Write in the tracer restart file
+      !                                          *******************************
+      IF( lrst_trc ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'trc : ATF trend at last time step for tracer budget written in tracer restart file ',   &
+            &                    'at it= ', kt,' date= ', ndastp
+         IF(lwp) WRITE(numout,*) '~~~~'
+         DO jn = 1, jptra
+            CALL iom_rstput( kt, nitrst, numrtw, 'atf_trend_'//TRIM(ctrcnm(jn)), ztrdt_m1(:,:,:,jn) )
+         END DO
+      ENDIF
+#endif
+
       !
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/TRP/trcrad.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/TRP/trcrad.F90	(revision 10155)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/TRP/trcrad.F90	(revision 10156)
@@ -18,4 +18,7 @@
    USE trdtra
    USE prtctl_trc          ! Print control for debbuging
+#if defined key_tracer_budget
+   USE iom
+#endif
 
    IMPLICIT NONE
@@ -51,4 +54,7 @@
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index      
       CHARACTER (len=22) :: charout
+      ! +++>>> FABM
+      INTEGER :: jn
+      ! FABM <<<+++
       !!----------------------------------------------------------------------
       !
@@ -65,5 +71,13 @@
       IF( lk_pisces  )   CALL trc_rad_sms( kt, trb, trn, jp_pcs0 , jp_pcs1, cpreserv='Y' )  ! PISCES model
       IF( lk_my_trc  )   CALL trc_rad_sms( kt, trb, trn, jp_myt0 , jp_myt1               )  ! MY_TRC model
-
+      ! +++>>> FABM
+      IF( lk_fabm  )   THEN
+        DO jn=1,jp_fabm ! state variable loop
+          IF (lk_rad_fabm(jn)) THEN
+           CALL trc_rad_sms( kt, trb, trn, jn+jp_fabm_m1 , jn+jp_fabm_m1 )
+          ENDIF
+        END DO
+      END IF
+      ! FABM <<<+++
       !
       IF(ln_ctl) THEN      ! print mean trends (used for debugging)
@@ -110,4 +124,7 @@
       REAL(wp) :: zcoef, ztrcorn, ztrmasn   !    "         "
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrtrdb, ztrtrdn   ! workspace arrays
+#if defined key_tracer_budget
+      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  ztrtrdb_m1 ! slwa 
+#endif
       REAL(wp) :: zs2rdt
       LOGICAL ::   lldebug = .FALSE.
@@ -116,4 +133,19 @@
  
       IF( l_trdtrc )  CALL wrk_alloc( jpi, jpj, jpk, ztrtrdb, ztrtrdn )
+#if defined key_tracer_budget
+      IF( kt == nittrc000 .AND. l_trdtrc) THEN
+         IF (.not. ALLOCATED(ztrtrdb_m1)) ALLOCATE( ztrtrdb_m1(jpi,jpj,jpk,jptra) )  ! slwa
+            DO jn = jp_sms0, jp_sms1
+               IF( ln_rsttr .AND.    &                     ! Restart: read in restart  file
+                  iom_varid( numrtr, 'rdb_trend_'//TRIM(ctrcnm(jn)), ldstop = .FALSE. ) > 0 ) THEN
+                  IF(lwp) WRITE(numout,*) '          nittrc000-nn_dttrc RDB tracer trend read for',TRIM(ctrcnm(jn))
+                  CALL iom_get( numrtr, jpdom_autoglo, 'rdb_trend_'//TRIM(ctrcnm(jn)), ztrtrdb_m1(:,:,:,jn) )   ! before tracer trend for rdb
+               ELSE
+                  IF(lwp) WRITE(numout,*) '          no nittrc000-nn_dttrc RDB tracer trend for',TRIM(ctrcnm(jn)),', setting to 0.'
+                  ztrtrdb_m1(:,:,:,jn)=0.0
+               ENDIF
+            END DO
+      ENDIF
+#endif
       
       IF( PRESENT( cpreserv )  ) THEN   !  total tracer concentration is preserved 
@@ -156,5 +188,15 @@
                ztrtrdb(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrdb(:,:,:) ) * zs2rdt
                ztrtrdn(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrdn(:,:,:) ) * zs2rdt 
+#if defined key_tracer_budget
+! slwa budget code
+               DO jk = 1, jpkm1
+                  ztrtrdb(:,:,jk) = ztrtrdb(:,:,jk) * e1t(:,:) * e2t(:,:) * fse3t(:,:,jk)
+                  ztrtrdn(:,:,jk) = ztrtrdn(:,:,jk) * e1t(:,:) * e2t(:,:) * fse3t(:,:,jk)
+               END DO
+               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb_m1(:,:,:,jn) )
+               ztrtrdb_m1(:,:,:,jn)=ztrtrdb(:,:,:)
+#else
                CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb )       ! Asselin-like trend handling
+#endif
                CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrdn )       ! standard     trend handling
               !
@@ -187,5 +229,15 @@
                ztrtrdb(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrdb(:,:,:) ) * zs2rdt
                ztrtrdn(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrdn(:,:,:) ) * zs2rdt 
+#if defined key_tracer_budget
+! slwa budget code
+               DO jk = 1, jpkm1
+                  ztrtrdb(:,:,jk) = ztrtrdb(:,:,jk) * e1t(:,:) * e2t(:,:) * fse3t(:,:,jk)
+                  ztrtrdn(:,:,jk) = ztrtrdn(:,:,jk) * e1t(:,:) * e2t(:,:) * fse3t(:,:,jk)
+               END DO
+               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb_m1(:,:,:,jn) )
+               ztrtrdb_m1(:,:,:,jn)=ztrtrdb(:,:,:)
+#else
                CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb )       ! Asselin-like trend handling
+#endif
                CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrdn )       ! standard     trend handling
               !
@@ -195,4 +247,18 @@
 
       ENDIF
+
+#if defined key_tracer_budget
+      !                                           Write in the tracer restart file
+      !                                          *******************************
+      IF( lrst_trc ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'trc : RDB trend at last time step for tracer budget written in tracer restart file ',   &
+            &                    'at it= ', kt,' date= ', ndastp
+         IF(lwp) WRITE(numout,*) '~~~~'
+         DO jn = jp_sms0, jp_sms1
+            CALL iom_rstput( kt, nitrst, numrtw, 'rdb_trend_'//TRIM(ctrcnm(jn)), ztrtrdb_m1(:,:,:,jn) )
+         END DO
+      ENDIF
+#endif
 
       IF( l_trdtrc )  CALL wrk_dealloc( jpi, jpj, jpk, ztrtrdb, ztrtrdn )
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90	(revision 10155)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90	(revision 10156)
@@ -113,4 +113,5 @@
            sbc_trc_b(:,:,:) = 0._wp
          ENDIF
+         sbc_trc(:,:,:) = 0._wp    !slwa initialise for vvl
       ELSE                                         ! Swap of forcing fields
          IF( ln_top_euler ) THEN
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 10155)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 10156)
@@ -27,4 +27,6 @@
    USE trcsbc          ! surface boundary condition          (trc_sbc routine)
    USE zpshde          ! partial step: hor. derivative       (zps_hde routine)
+   USE trcbdy          ! BDY open boundaries
+   USE bdy_par, only: lk_bdy
 
 #if defined key_agrif
@@ -68,4 +70,5 @@
          IF( ln_trcdmp )        CALL trc_dmp( kstp )            ! internal damping trends
          IF( ln_trcdmp_clo )    CALL trc_dmp_clo( kstp )        ! internal damping trends on closed seas only
+         IF( lk_bdy )           CALL trc_bdy_dmp( kstp )        ! BDY damping trends
                                 CALL trc_adv( kstp )            ! horizontal & vertical advection 
                                 CALL trc_ldf( kstp )            ! lateral mixing
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/TRP/trdtrc.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/TRP/trdtrc.F90	(revision 10155)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/TRP/trdtrc.F90	(revision 10156)
@@ -102,4 +102,13 @@
          END SELECT
                                           cltra = TRIM(cltra)//TRIM(ctrcnm(kjn))
+         ! +++>>>FABM
+#if defined key_tracer_budget
+! for outputting depth integrated
+         SELECT CASE( ktrd )
+         CASE( jptra_xad, jptra_yad, jptra_zad  ) 
+           cltra = TRIM(cltra)//"_e3t"
+         END SELECT
+#endif
+         ! FABM <<<+++
                                           CALL iom_put( cltra,  ptrtrd(:,:,:) )
          !
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/TRP/trdtrc_oce.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/TRP/trdtrc_oce.F90	(revision 10155)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/TRP/trdtrc_oce.F90	(revision 10156)
@@ -22,5 +22,10 @@
    CHARACTER(len=50) ::  cn_trdrst_trc_in     !: suffix of pass. tracer restart name (input)
    CHARACTER(len=50) ::  cn_trdrst_trc_out    !: suffix of pass. tracer restart name (output)
-   LOGICAL, DIMENSION(jptra) ::   ln_trdtrc   !: large trends diagnostic to write or not (namelist)
+   ! --->>> FABM
+   ! LOGICAL, DIMENSION(jptra) ::   ln_trdtrc   !: large trends diagnostic to write or not (namelist)
+   ! FABM <<<---
+   ! +++>>> FABM
+   LOGICAL, DIMENSION(jpmaxtrc) ::   ln_trdtrc   !: large trends diagnostic to write or not (namelist)
+   ! FABM <<<+++
 
 # if defined key_trdtrc && defined key_iomput
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/par_trc.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/par_trc.F90	(revision 10155)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/par_trc.F90	(revision 10156)
@@ -15,4 +15,7 @@
    USE par_cfc       ! CFC 11 and 12 tracers
    USE par_my_trc    ! user defined passive tracers
+   ! +++>>> FABM
+   USE par_fabm      ! FABM
+   ! FABM <<<+++
 
    IMPLICIT NONE
@@ -24,9 +27,21 @@
    ! Passive tracers : Total size
    ! ---------------               ! total number of passive tracers, of 2d and 3d output and trend arrays
-   INTEGER, PUBLIC,  PARAMETER ::   jptra    =  jp_pisces     + jp_cfc     + jp_c14b    + jp_my_trc
-   INTEGER, PUBLIC,  PARAMETER ::   jpdia2d  =  jp_pisces_2d  + jp_cfc_2d  + jp_c14b_2d + jp_my_trc_2d
-   INTEGER, PUBLIC,  PARAMETER ::   jpdia3d  =  jp_pisces_3d  + jp_cfc_3d  + jp_c14b_3d + jp_my_trc_3d
+   ! --->>> FABM
+   ! INTEGER, PUBLIC,  PARAMETER ::   jptra    =  jp_pisces     + jp_cfc     + jp_c14b    + jp_my_trc
+   ! INTEGER, PUBLIC,  PARAMETER ::   jpdia2d  =  jp_pisces_2d  + jp_cfc_2d  + jp_c14b_2d + jp_my_trc_2d
+   ! INTEGER, PUBLIC,  PARAMETER ::   jpdia3d  =  jp_pisces_3d  + jp_cfc_3d  + jp_c14b_3d + jp_my_trc_3d
+   ! FABM <<<---
+   ! +++>>> FABM
+   INTEGER, PUBLIC  ::   jptra    =  jp_pisces     + jp_cfc     + jp_c14b    + jp_my_trc
+   INTEGER, PUBLIC  ::   jpdia2d  =  jp_pisces_2d  + jp_cfc_2d  + jp_c14b_2d + jp_my_trc_2d
+   INTEGER, PUBLIC  ::   jpdia3d  =  jp_pisces_3d  + jp_cfc_3d  + jp_c14b_3d + jp_my_trc_3d
+   ! FABM <<<+++
    !                     ! total number of sms diagnostic arrays
-   INTEGER, PUBLIC,  PARAMETER ::   jpdiabio =  jp_pisces_trd + jp_cfc_trd + jp_c14b_trd + jp_my_trc_trd
+   ! --->>> FABM
+   ! INTEGER, PUBLIC,  PARAMETER ::   jpdiabio =  jp_pisces_trd + jp_cfc_trd + jp_c14b_trd + jp_my_trc_trd
+   ! FABM <<<---
+   ! +++>>> FABM
+   INTEGER, PUBLIC  ::   jpdiabio =  jp_pisces_trd + jp_cfc_trd + jp_c14b_trd + jp_my_trc_trd
+   ! FABM <<<+++
    
    !  1D configuration ("key_c1d")
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trc.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 10155)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 10156)
@@ -14,4 +14,7 @@
    USE par_oce
    USE par_trc
+#if defined key_bdy
+   USE bdy_oce, only: nb_bdy, OBC_DATA
+#endif
    
    IMPLICIT NONE
@@ -80,7 +83,14 @@
    END TYPE
 
-   REAL(wp), DIMENSION(jptra), PUBLIC         :: trc_ice_ratio, & ! ice-ocean tracer ratio
+   ! --->>> FABM 
+   !REAL(wp), DIMENSION(jptra), PUBLIC         :: trc_ice_ratio, & ! ice-ocean tracer ratio
+   !                                              trc_ice_prescr   ! prescribed ice trc cc
+   !CHARACTER(len=2), DIMENSION(jptra), PUBLIC :: cn_trc_o ! choice of ocean tracer cc
+   ! FABM <<<---
+   ! +++>>> FABM 
+   REAL(wp), DIMENSION(jpmaxtrc), PUBLIC         :: trc_ice_ratio, & ! ice-ocean tracer ratio
                                                  trc_ice_prescr   ! prescribed ice trc cc
-   CHARACTER(len=2), DIMENSION(jptra), PUBLIC :: cn_trc_o ! choice of ocean tracer cc
+   CHARACTER(len=2), DIMENSION(jpmaxtrc), PUBLIC :: cn_trc_o ! choice of ocean tracer cc
+   ! FABM <<<+++
 
    !! information for outputs
@@ -90,6 +100,22 @@
        CHARACTER(len = 80)  :: cllname  !: long name
        CHARACTER(len = 20)  :: clunit   !: unit
-       LOGICAL              :: llinit   !: read in a file or not
-       LOGICAL              :: llsave   !: save the tracer or not
+! --->>> FABM
+!       LOGICAL              :: llinit   !: read in a file or not
+!!#if defined  key_my_trc
+!       LOGICAL              :: llsbc   !: read in a file or not
+!       LOGICAL              :: llcbc   !: read in a file or not
+!       LOGICAL              :: llobc   !: read in a file or not
+!#endif
+!       LOGICAL              :: llsave   !: save the tracer or not
+! FABM <<<---
+! +++ FABM
+       LOGICAL              :: llinit=.FALSE.   !: read in a file or not
+#if defined  key_fabm
+       LOGICAL              :: llsbc=.FALSE.   !: read in a file or not
+       LOGICAL              :: llcbc=.FALSE.   !: read in a file or not
+       LOGICAL              :: llobc=.FALSE.   !: read in a file or not
+#endif
+       LOGICAL              :: llsave=.FALSE.   !: save the tracer or not
+! FABM <<<+++
    END TYPE PTRACER
    CHARACTER(len = 20), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)    ::  ctrcnm         !: tracer name 
@@ -191,4 +217,11 @@
 # endif
    !
+#if defined key_bdy
+   CHARACTER(len=20), PUBLIC, ALLOCATABLE,  SAVE,  DIMENSION(:)   ::  cn_trc_dflt          ! Default OBC condition for all tracers
+   CHARACTER(len=20), PUBLIC, ALLOCATABLE,  SAVE,  DIMENSION(:)   ::  cn_trc               ! Choice of boundary condition for tracers
+   INTEGER,           PUBLIC, ALLOCATABLE,  SAVE,  DIMENSION(:)   ::  nn_trcdmp_bdy        !: =T Tracer damping
+   ! External data structure of BDY for TOP. Available elements: cn_obc, ll_trc, trcnow, dmp
+   TYPE(OBC_DATA),    PUBLIC, ALLOCATABLE, DIMENSION(:,:), TARGET ::  trcdta_bdy           !: bdy external data (local process)
+#endif
 
    !!----------------------------------------------------------------------
@@ -213,4 +246,16 @@
          &      cvol(jpi,jpj,jpk)     , rdttrc(jpk)           , trai(jptra)           ,       &
          &      ctrcnm(jptra)         , ctrcln(jptra)         , ctrcun(jptra)         ,       & 
+! --->>> FABM
+!!#if defined key_my_trc
+! FABM <<<---
+! +++>>> FABM
+#if defined key_fabm
+! FABM <<<+++
+         &      ln_trc_sbc(jptra)     , ln_trc_cbc(jptra)     , ln_trc_obc(jptra)     ,       &
+#endif
+#if defined key_bdy
+         &      cn_trc_dflt(nb_bdy)   , cn_trc(nb_bdy)        , nn_trcdmp_bdy(nb_bdy) ,       &
+         &      trcdta_bdy(jptra,nb_bdy)                                              ,       &
+#endif
          &      ln_trc_ini(jptra)     , ln_trc_wri(jptra)     , qsr_mean(jpi,jpj)     ,  STAT = trc_alloc  )  
 
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcbc.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcbc.F90	(revision 10155)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcbc.F90	(revision 10156)
@@ -4,10 +4,12 @@
    !! TOP :  module for passive tracer boundary conditions
    !!=====================================================================
-   !!----------------------------------------------------------------------
-#if  defined key_top 
+   !! History :  3.5 !  2014-04  (M. Vichi, T. Lovato)  Original
+   !!            3.6 !  2015-03  (T . Lovato) Revision and BDY support
+   !!----------------------------------------------------------------------
+#if defined key_top
    !!----------------------------------------------------------------------
    !!   'key_top'                                                TOP model 
    !!----------------------------------------------------------------------
-   !!   trc_dta    : read and time interpolated passive tracer data
+   !!   trc_bc       : read and time interpolated tracer Boundary Conditions
    !!----------------------------------------------------------------------
    USE par_trc       !  passive tracers parameters
@@ -17,4 +19,7 @@
    USE lib_mpp       !  MPP library
    USE fldread       !  read input fields
+#if defined key_bdy
+   USE bdy_oce, only: nb_bdy , idx_bdy, ln_coords_file, rn_time_dmp, rn_time_dmp_out
+#endif
 
    IMPLICIT NONE
@@ -24,24 +29,22 @@
    PUBLIC   trc_bc_read    ! called in trcstp.F90 or within
 
-   INTEGER  , SAVE, PUBLIC                             :: nb_trcobc   ! number of tracers with open BC
-   INTEGER  , SAVE, PUBLIC                             :: nb_trcsbc   ! number of tracers with surface BC
-   INTEGER  , SAVE, PUBLIC                             :: nb_trccbc   ! number of tracers with coastal BC
+   INTEGER  , SAVE, PUBLIC                             :: nb_trcobc    ! number of tracers with open BC
+   INTEGER  , SAVE, PUBLIC                             :: nb_trcsbc    ! number of tracers with surface BC
+   INTEGER  , SAVE, PUBLIC                             :: nb_trccbc    ! number of tracers with coastal BC
    INTEGER  , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: n_trc_indobc ! index of tracer with OBC data
    INTEGER  , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: n_trc_indsbc ! index of tracer with SBC data
    INTEGER  , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: n_trc_indcbc ! index of tracer with CBC data
-   INTEGER  , SAVE, PUBLIC                             :: ntra_obc     ! MAX( 1, nb_trcxxx ) to avoid compilation error with bounds checking
-   INTEGER  , SAVE, PUBLIC                             :: ntra_sbc     ! MAX( 1, nb_trcxxx ) to avoid compilation error with bounds checking
-   INTEGER  , SAVE, PUBLIC                             :: ntra_cbc     ! MAX( 1, nb_trcxxx ) to avoid compilation error with bounds checking
-   REAL(wp) , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: rf_trofac   ! multiplicative factor for OBCtracer values
-   TYPE(FLD), SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: sf_trcobc   ! structure of data input OBC (file informations, fields read)
-   REAL(wp) , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: rf_trsfac   ! multiplicative factor for SBC tracer values
-   TYPE(FLD), SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: sf_trcsbc   ! structure of data input SBC (file informations, fields read)
-   REAL(wp) , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: rf_trcfac   ! multiplicative factor for CBC tracer values
-   TYPE(FLD), SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: sf_trccbc   ! structure of data input CBC (file informations, fields read)
+   REAL(wp) , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: rf_trsfac    ! multiplicative factor for SBC tracer values
+   TYPE(FLD), SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: sf_trcsbc    ! structure of data input SBC (file informations, fields read)
+   REAL(wp) , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: rf_trcfac    ! multiplicative factor for CBC tracer values
+   TYPE(FLD), SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: sf_trccbc    ! structure of data input CBC (file informations, fields read)
+   REAL(wp) , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:,:)  :: rf_trofac    ! multiplicative factor for OBCtracer values
+   TYPE(FLD), SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:,:), TARGET  :: sf_trcobc    ! structure of data input OBC (file informations, fields read)
+   TYPE(MAP_POINTER), ALLOCATABLE, DIMENSION(:,:) :: nbmap_ptr   ! array of pointers to nbmap
 
    !! * Substitutions
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 3.6 , NEMO Consortium (2015)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -60,12 +63,15 @@
       !
       INTEGER,INTENT(IN) :: ntrc                           ! number of tracers
-      INTEGER            :: jl, jn                         ! dummy loop indices
+      INTEGER            :: jl, jn , ib, ibd, ii, ij, ik   ! dummy loop indices
       INTEGER            :: ierr0, ierr1, ierr2, ierr3     ! temporary integers
-      INTEGER            ::  ios                           ! Local integer output status for namelist read
+      INTEGER            :: ios                            ! Local integer output status for namelist read
+      INTEGER            :: nblen, igrd                    ! support arrays for BDY
       CHARACTER(len=100) :: clndta, clntrc
       !
-      CHARACTER(len=100) :: cn_dir
+      CHARACTER(len=100) :: cn_dir_sbc, cn_dir_cbc, cn_dir_obc
+
       TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: slf_i  ! local array of namelist informations on the fields to read
-      TYPE(FLD_N), DIMENSION(jpmaxtrc) :: sn_trcobc    ! open
+      TYPE(FLD_N), DIMENSION(jpmaxtrc,2) :: sn_trcobc  ! open
+      TYPE(FLD_N), DIMENSION(jpmaxtrc) :: sn_trcobc2   ! to read in multiple (2) open bdy
       TYPE(FLD_N), DIMENSION(jpmaxtrc) :: sn_trcsbc    ! surface
       TYPE(FLD_N), DIMENSION(jpmaxtrc) :: sn_trccbc    ! coastal
@@ -74,8 +80,16 @@
       REAL(wp)   , DIMENSION(jpmaxtrc) :: rn_trcfac    ! multiplicative factor for tracer values
       !!
-      NAMELIST/namtrc_bc/ cn_dir, sn_trcobc, rn_trofac, sn_trcsbc, rn_trsfac, sn_trccbc, rn_trcfac 
+      NAMELIST/namtrc_bc/ cn_dir_sbc, cn_dir_cbc, cn_dir_obc, sn_trcobc2, rn_trofac, sn_trcsbc, rn_trsfac, sn_trccbc, rn_trcfac
+#if defined key_bdy
+      NAMELIST/namtrc_bdy/ cn_trc_dflt, cn_trc, nn_trcdmp_bdy
+#endif
       !!----------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('trc_bc_init')
       !
+      IF( lwp ) THEN
+         WRITE(numout,*) ' '
+         WRITE(numout,*) 'trc_bc_init : Tracers Boundary Conditions (BC)'
+         WRITE(numout,*) '~~~~~~~~~~~ '
+      ENDIF
       !  Initialisation and local array allocation
       ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0  
@@ -107,31 +121,5 @@
       n_trc_indcbc(:) = 0
       !
-      DO jn = 1, ntrc
-         IF( ln_trc_obc(jn) ) THEN
-             nb_trcobc       = nb_trcobc + 1 
-             n_trc_indobc(jn) = nb_trcobc 
-         ENDIF
-         IF( ln_trc_sbc(jn) ) THEN
-             nb_trcsbc       = nb_trcsbc + 1
-             n_trc_indsbc(jn) = nb_trcsbc
-         ENDIF
-         IF( ln_trc_cbc(jn) ) THEN
-             nb_trccbc       = nb_trccbc + 1
-             n_trc_indcbc(jn) = nb_trccbc
-         ENDIF
-      ENDDO
-      ntra_obc = MAX( 1, nb_trcobc )   ! To avoid compilation error with bounds checking
-      IF( lwp ) WRITE(numout,*) ' '
-      IF( lwp ) WRITE(numout,*) ' Number of passive tracers to be initialized with open boundary data :', nb_trcobc
-      IF( lwp ) WRITE(numout,*) ' '
-      ntra_sbc = MAX( 1, nb_trcsbc )   ! To avoid compilation error with bounds checking
-      IF( lwp ) WRITE(numout,*) ' '
-      IF( lwp ) WRITE(numout,*) ' Number of passive tracers to be initialized with surface boundary data :', nb_trcsbc
-      IF( lwp ) WRITE(numout,*) ' '
-      ntra_cbc = MAX( 1, nb_trccbc )   ! To avoid compilation error with bounds checking
-      IF( lwp ) WRITE(numout,*) ' '
-      IF( lwp ) WRITE(numout,*) ' Number of passive tracers to be initialized with coastal boundary data :', nb_trccbc
-      IF( lwp ) WRITE(numout,*) ' '
-
+      ! Read Boundary Conditions Namelists
       REWIND( numnat_ref )              ! Namelist namtrc_bc in reference namelist : Passive tracer data structure
       READ  ( numnat_ref, namtrc_bc, IOSTAT = ios, ERR = 901)
@@ -139,58 +127,175 @@
 
       REWIND( numnat_cfg )              ! Namelist namtrc_bc in configuration namelist : Passive tracer data structure
-      READ  ( numnat_cfg, namtrc_bc, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_bc in configuration namelist', lwp )
-      IF(lwm) WRITE ( numont, namtrc_bc )
-
-      ! print some information for each 
+#if defined key_bdy
+      DO ib = 1, nb_bdy
+#endif
+        READ  ( numnat_cfg, namtrc_bc, IOSTAT = ios, ERR = 902 )
+902     IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_bc in configuration namelist', lwp )
+        IF(lwm) WRITE ( numont, namtrc_bc )
+#if defined key_bdy
+        sn_trcobc(:,ib)=sn_trcobc2(:)
+      ENDDO
+#endif
+
+#if defined key_bdy
+      REWIND( numnat_ref )              ! Namelist namtrc_bc in reference namelist : Passive tracer data structure
+      READ  ( numnat_ref, namtrc_bdy, IOSTAT = ios, ERR = 903)
+903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_bdy in reference namelist', lwp )
+
+      REWIND( numnat_cfg )              ! Namelist namtrc_bc in configuration namelist : Passive tracer data structure
+      READ  ( numnat_cfg, namtrc_bdy, IOSTAT = ios, ERR = 904 )
+904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_bdy in configuration namelist', lwp )
+      IF(lwm) WRITE ( numont, namtrc_bdy )
+      ! setup up preliminary informations for BDY structure
+      DO jn = 1, ntrc
+         DO ib = 1, nb_bdy
+            ! Set type of obc in BDY data structure (around here we may plug user override of obc type from nml)
+            IF ( ln_trc_obc(jn) ) THEN
+               trcdta_bdy(jn,ib)%cn_obc = TRIM( cn_trc(ib) )
+            ELSE
+               trcdta_bdy(jn,ib)%cn_obc = TRIM( cn_trc_dflt(ib) )
+            ENDIF
+            ! set damping use in BDY data structure
+            trcdta_bdy(jn,ib)%dmp = .false.
+            IF(nn_trcdmp_bdy(ib) .EQ. 1 .AND. ln_trc_obc(jn) ) trcdta_bdy(jn,ib)%dmp = .true.
+            IF(nn_trcdmp_bdy(ib) .EQ. 2 ) trcdta_bdy(jn,ib)%dmp = .true.
+            IF(trcdta_bdy(jn,ib)%cn_obc == 'frs' .AND. nn_trcdmp_bdy(ib) .NE. 0 )  &
+                & CALL ctl_stop( 'Use FRS OR relaxation' )
+            IF (nn_trcdmp_bdy(ib) .LT. 0 .OR. nn_trcdmp_bdy(ib) .GT. 2)            &
+                & CALL ctl_stop( 'Not a valid option for nn_trcdmp_bdy. Allowed: 0,1,2.' )
+         ENDDO
+      ENDDO
+
+#else
+      ! Force all tracers OBC to false if bdy not used
+      ln_trc_obc = .false.
+#endif
+      ! compose BC data indexes
+      DO jn = 1, ntrc
+         IF( ln_trc_obc(jn) ) THEN
+             nb_trcobc       = nb_trcobc + 1  ; n_trc_indobc(jn) = nb_trcobc
+         ENDIF
+         IF( ln_trc_sbc(jn) ) THEN
+             nb_trcsbc       = nb_trcsbc + 1  ; n_trc_indsbc(jn) = nb_trcsbc
+         ENDIF
+         IF( ln_trc_cbc(jn) ) THEN
+             nb_trccbc       = nb_trccbc + 1  ; n_trc_indcbc(jn) = nb_trccbc
+         ENDIF
+      ENDDO
+      ntra_obc = MAX( 1, nb_trcobc )   ! To avoid compilation error with bounds checking
+      IF( lwp ) WRITE(numout,*) ' '
+      IF( lwp ) WRITE(numout,*) ' Number of passive tracers to be initialized with open boundary data :', nb_trcobc
+      IF( lwp ) WRITE(numout,*) ' '
+      ntra_sbc = MAX( 1, nb_trcsbc )   ! To avoid compilation error with bounds checking
+      IF( lwp ) WRITE(numout,*) ' '
+      IF( lwp ) WRITE(numout,*) ' Number of passive tracers to be initialized with surface boundary data :', nb_trcsbc
+      IF( lwp ) WRITE(numout,*) ' '
+      ntra_cbc = MAX( 1, nb_trccbc )   ! To avoid compilation error with bounds checking
+      IF( lwp ) WRITE(numout,*) ' '
+      IF( lwp ) WRITE(numout,*) ' Number of passive tracers to be initialized with coastal boundary data :', nb_trccbc
+      IF( lwp ) WRITE(numout,*) ' '
+
+      ! Print summmary of Boundary Conditions
       IF( lwp ) THEN
-         DO jn = 1, ntrc
-            IF( ln_trc_obc(jn) )  THEN    
-               clndta = TRIM( sn_trcobc(jn)%clvar ) 
-               IF(lwp) WRITE(numout,*) 'Preparing to read OBC data file for passive tracer number :', jn, ' name : ', clndta, & 
-               &               ' multiplicative factor : ', rn_trofac(jn)
-            ENDIF
-            IF( ln_trc_sbc(jn) )  THEN    
-               clndta = TRIM( sn_trcsbc(jn)%clvar ) 
-               IF(lwp) WRITE(numout,*) 'Preparing to read SBC data file for passive tracer number :', jn, ' name : ', clndta, & 
-               &               ' multiplicative factor : ', rn_trsfac(jn)
-            ENDIF
-            IF( ln_trc_cbc(jn) )  THEN    
-               clndta = TRIM( sn_trccbc(jn)%clvar ) 
-               IF(lwp) WRITE(numout,*) 'Preparing to read CBC data file for passive tracer number :', jn, ' name : ', clndta, & 
-               &               ' multiplicative factor : ', rn_trcfac(jn)
-            ENDIF
-         END DO
-      ENDIF
-      !
-      ! The following code is written this way to reduce memory usage and repeated for each boundary data
-      ! MAV: note that this is just a placeholder and the dimensions must be changed according to 
-      !      what will be done with BDY. A new structure will probably need to be included
-      !
+         WRITE(numout,*) ' '
+         WRITE(numout,'(a,i3)') '   Total tracers to be initialized with SURFACE BCs data:', nb_trcsbc
+         IF ( nb_trcsbc > 0 ) THEN
+            WRITE(numout,*) '   #trc        NAME        Boundary     Mult.Fact. '
+            DO jn = 1, ntrc
+               IF ( ln_trc_sbc(jn) ) WRITE(numout,9001) jn, TRIM( sn_trcsbc(jn)%clvar ), 'SBC', rn_trsfac(jn)
+            ENDDO
+         ENDIF
+         WRITE(numout,'(2a)') '   SURFACE BC data repository : ', TRIM(cn_dir_sbc)
+
+         WRITE(numout,*) ' '
+         WRITE(numout,'(a,i3)') '   Total tracers to be initialized with COASTAL BCs data:', nb_trccbc
+         IF ( nb_trccbc > 0 ) THEN
+            WRITE(numout,*) '   #trc        NAME        Boundary     Mult.Fact. '
+            DO jn = 1, ntrc
+               IF ( ln_trc_cbc(jn) ) WRITE(numout, 9001) jn, TRIM( sn_trccbc(jn)%clvar ), 'CBC', rn_trcfac(jn)
+            ENDDO
+         ENDIF
+         WRITE(numout,'(2a)') '   COASTAL BC data repository : ', TRIM(cn_dir_cbc)
+
+         WRITE(numout,*) ' '
+         WRITE(numout,'(a,i3)') '   Total tracers to be initialized with OPEN BCs data:', nb_trcobc
+#if defined key_bdy
+         IF ( nb_trcobc > 0 ) THEN
+            WRITE(numout,*) '   #trc        NAME        Boundary     Mult.Fact.   OBC Settings'
+            DO jn = 1, ntrc
+               DO ib = 1, nb_bdy
+                 IF ( ln_trc_obc(jn) )  WRITE(numout, 9001) jn, TRIM( sn_trcobc(jn,ib)%clvar ), 'OBC', rn_trofac(jn), (trcdta_bdy(jn,ib)%cn_obc)
+               ENDDO
+               !IF ( ln_trc_obc(jn) )  WRITE(numout, 9001) jn, TRIM( sn_trcobc(jn,ib)%clvar ), 'OBC', rn_trofac(jn), (trcdta_bdy(jn,ib)%cn_obc,ib=1,nb_bdy)
+               IF ( .NOT. ln_trc_obc(jn) )  WRITE(numout, 9002) jn, 'Set data to IC and use default condition', (trcdta_bdy(jn,ib)%cn_obc,ib=1,nb_bdy)
+            ENDDO
+            WRITE(numout,*) ' '
+            DO ib = 1, nb_bdy
+                IF (nn_trcdmp_bdy(ib) .EQ. 0) WRITE(numout,9003) '   Boundary ',ib,' -> NO damping of tracers'
+                IF (nn_trcdmp_bdy(ib) .EQ. 1) WRITE(numout,9003) '   Boundary ',ib,' -> damping ONLY for tracers with external data provided'
+                IF (nn_trcdmp_bdy(ib) .EQ. 2) WRITE(numout,9003) '   Boundary ',ib,' -> damping of ALL tracers'
+                IF (nn_trcdmp_bdy(ib) .GT. 0) THEN
+                   WRITE(numout,9003) '     USE damping parameters from nambdy for boundary ', ib,' : '
+                   WRITE(numout,'(a,f10.2,a)') '     - Inflow damping time scale  : ',rn_time_dmp(ib),' days'
+                   WRITE(numout,'(a,f10.2,a)') '     - Outflow damping time scale : ',rn_time_dmp_out(ib),' days'
+                ENDIF
+            ENDDO
+         ENDIF
+#endif
+         WRITE(numout,'(2a)') '   OPEN BC data repository : ', TRIM(cn_dir_obc)
+      ENDIF
+9001  FORMAT(2x,i5, 3x, a15, 3x, a5, 6x, e11.3, 4x, 10a13)
+9002  FORMAT(2x,i5, 3x, a41, 3x, 10a13)
+9003  FORMAT(a, i5, a)
+
+      !
+#if defined key_bdy
       ! OPEN Lateral boundary conditions
-      IF( nb_trcobc > 0 ) THEN       !  allocate only if the number of tracer to initialise is greater than zero
-         ALLOCATE( sf_trcobc(nb_trcobc), rf_trofac(nb_trcobc), STAT=ierr1 )
+      IF( nb_trcobc > 0 ) THEN 
+         ALLOCATE ( sf_trcobc(nb_trcobc,nb_bdy), rf_trofac(nb_trcobc,nb_bdy), nbmap_ptr(nb_trcobc,nb_bdy), STAT=ierr1 )
          IF( ierr1 > 0 ) THEN
-            CALL ctl_stop( 'trc_bc_init: unable to allocate  sf_trcobc structure' )   ;   RETURN
-         ENDIF
-         !
-         DO jn = 1, ntrc
-            IF( ln_trc_obc(jn) ) THEN      ! update passive tracers arrays with input data read from file
-               jl = n_trc_indobc(jn)
-               slf_i(jl)    = sn_trcobc(jn)
-               rf_trofac(jl) = rn_trofac(jn)
-                                            ALLOCATE( sf_trcobc(jl)%fnow(jpi,jpj,jpk)   , STAT=ierr2 )
-               IF( sn_trcobc(jn)%ln_tint )  ALLOCATE( sf_trcobc(jl)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 )
-               IF( ierr2 + ierr3 > 0 ) THEN
-                 CALL ctl_stop( 'trc_bc_init : unable to allocate passive tracer OBC data arrays' )   ;   RETURN
+            CALL ctl_stop( 'trc_bc_init: unable to allocate sf_trcobc structure' )   ;   RETURN
+         ENDIF
+
+         igrd = 1                       ! Everything is at T-points here
+
+         DO ib = 1, nb_bdy
+            DO jn = 1, ntrc
+
+               nblen = idx_bdy(ib)%nblen(igrd)
+
+               IF ( ln_trc_obc(jn) ) THEN
+               ! Initialise from external data
+                  jl = n_trc_indobc(jn)
+                  slf_i(jl)    = sn_trcobc(jn,ib)
+                  rf_trofac(jl,ib) = rn_trofac(jn)
+                                               ALLOCATE( sf_trcobc(jl,ib)%fnow(nblen,1,jpk)   , STAT=ierr2 )
+                  IF( sn_trcobc(jn,ib)%ln_tint )  ALLOCATE( sf_trcobc(jl,ib)%fdta(nblen,1,jpk,2) , STAT=ierr3 )
+                  IF( ierr2 + ierr3 > 0 ) THEN
+                    CALL ctl_stop( 'trc_bc_init : unable to allocate passive tracer OBC data arrays' )   ;   RETURN
+                  ENDIF
+                  trcdta_bdy(jn,ib)%trc => sf_trcobc(jl,ib)%fnow(:,1,:)
+                  trcdta_bdy(jn,ib)%rn_fac = rf_trofac(jl,ib)
+                  ! create OBC mapping array
+                  nbmap_ptr(jl,ib)%ptr => idx_bdy(ib)%nbmap(:,igrd)
+                  nbmap_ptr(jl,ib)%ll_unstruc = ln_coords_file(igrd)
+               ELSE
+               ! Initialise obc arrays from initial conditions
+                  ALLOCATE ( trcdta_bdy(jn,ib)%trc(nblen,jpk) )
+                  DO ibd = 1, nblen
+                     DO ik = 1, jpkm1
+                        ii = idx_bdy(ib)%nbi(ibd,igrd)
+                        ij = idx_bdy(ib)%nbj(ibd,igrd)
+                        trcdta_bdy(jn,ib)%trc(ibd,ik) = trn(ii,ij,ik,jn) * tmask(ii,ij,ik)
+                     END DO
+                  END DO
+                  trcdta_bdy(jn,ib)%rn_fac = 1._wp
                ENDIF
-            ENDIF
-            !   
+            ENDDO
+            CALL fld_fill( sf_trcobc(:,ib), slf_i, cn_dir_obc, 'trc_bc_init', 'Passive tracer OBC data', 'namtrc_bc' )
          ENDDO
-         !                         ! fill sf_trcdta with slf_i and control print
-         CALL fld_fill( sf_trcobc, slf_i, cn_dir, 'trc_bc_init', 'Passive tracer OBC data', 'namtrc_bc' )
-         !
-      ENDIF
-      !
+
+      ENDIF
+#endif
       ! SURFACE Boundary conditions
       IF( nb_trcsbc > 0 ) THEN       !  allocate only if the number of tracer to initialise is greater than zero
@@ -214,5 +319,5 @@
          ENDDO
          !                         ! fill sf_trcsbc with slf_i and control print
-         CALL fld_fill( sf_trcsbc, slf_i, cn_dir, 'trc_bc_init', 'Passive tracer SBC data', 'namtrc_bc' )
+         CALL fld_fill( sf_trcsbc, slf_i, cn_dir_sbc, 'trc_bc_init', 'Passive tracer SBC data', 'namtrc_bc' )
          !
       ENDIF
@@ -239,5 +344,5 @@
          ENDDO
          !                         ! fill sf_trccbc with slf_i and control print
-         CALL fld_fill( sf_trccbc, slf_i, cn_dir, 'trc_bc_init', 'Passive tracer CBC data', 'namtrc_bc' )
+         CALL fld_fill( sf_trccbc, slf_i, cn_dir_cbc, 'trc_bc_init', 'Passive tracer CBC data', 'namtrc_bc' )
          !
       ENDIF
@@ -249,5 +354,5 @@
 
 
-   SUBROUTINE trc_bc_read(kt)
+   SUBROUTINE trc_bc_read(kt, jit)
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE trc_bc_init  ***
@@ -264,26 +369,64 @@
       !! * Arguments
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
-
+      INTEGER, INTENT( in ), OPTIONAL ::   jit   ! subcycle time-step index (for timesplitting option)
+      INTEGER :: ib
       !!---------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('trc_bc_read')
 
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'trc_bc_read : Surface boundary conditions for passive tracers.'
-         IF(lwp) WRITE(numout,*) '~~~~~~~ '
-      ENDIF
-
-      ! OPEN boundary conditions: DOES NOT WORK. Waiting for stable BDY
-      IF( nb_trcobc > 0 ) THEN
-        if (lwp) write(numout,'(a,i5,a,i5)') '   reading OBC data for ', nb_trcobc ,' variables at step ', kt
-        CALL fld_read(kt,1,sf_trcobc)
-        ! vertical interpolation on s-grid and partial step to be added
-      ENDIF
-
-      ! SURFACE boundary conditions       
-      IF( nb_trcsbc > 0 ) THEN
-        if (lwp) write(numout,'(a,i5,a,i5)') '   reading SBC data for ', nb_trcsbc ,' variables at step ', kt
-        CALL fld_read(kt,1,sf_trcsbc)
+      IF( kt == nit000 .AND. lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'trc_bc_read : Surface boundary conditions for passive tracers.'
+         WRITE(numout,*) '~~~~~~~~~~~ '
+      ENDIF
+
+      IF ( PRESENT(jit) ) THEN 
+
+#ifdef key_bdy
+         ! OPEN boundary conditions (use time_offset=+1 as they are applied at the end of the step)
+         IF( nb_trcobc > 0 ) THEN
+           if (lwp) write(numout,'(a,i5,a,i10)') '   reading OBC data for ', nb_trcobc ,' variable(s) at step ', kt
+           DO ib = 1,nb_bdy
+             CALL fld_read(kt=kt, kn_fsbc=1, sd=sf_trcobc(:,ib), map=nbmap_ptr(:,ib), kit=jit, kt_offset=+1)
+           ENDDO
+         ENDIF
+#endif
+
+         ! SURFACE boundary conditions
+         IF( nb_trcsbc > 0 ) THEN
+           if (lwp) write(numout,'(a,i5,a,i10)') '   reading SBC data for ', nb_trcsbc ,' variable(s) at step ', kt
+           CALL fld_read(kt=kt, kn_fsbc=1, sd=sf_trcsbc, kit=jit)
+         ENDIF
+
+         ! COASTAL boundary conditions
+         IF( nb_trccbc > 0 ) THEN
+           if (lwp) write(numout,'(a,i5,a,i10)') '   reading CBC data for ', nb_trccbc ,' variable(s) at step ', kt
+           CALL fld_read(kt=kt, kn_fsbc=1, sd=sf_trccbc, kit=jit)
+         ENDIF
+
+      ELSE
+
+#ifdef key_bdy
+         ! OPEN boundary conditions (use time_offset=+1 as they are applied at the end of the step)
+         IF( nb_trcobc > 0 ) THEN
+           if (lwp) write(numout,'(a,i5,a,i10)') '   reading OBC data for ', nb_trcobc ,' variable(s) at step ', kt
+           DO ib = 1,nb_bdy
+             CALL fld_read(kt=kt, kn_fsbc=1, sd=sf_trcobc(:,ib), map=nbmap_ptr(:,ib), kt_offset=+1)
+           ENDDO
+         ENDIF
+#endif
+
+         ! SURFACE boundary conditions
+         IF( nb_trcsbc > 0 ) THEN
+           if (lwp) write(numout,'(a,i5,a,i10)') '   reading SBC data for ', nb_trcsbc ,' variable(s) at step ', kt
+           CALL fld_read(kt=kt, kn_fsbc=1, sd=sf_trcsbc)
+         ENDIF
+
+         ! COASTAL boundary conditions
+         IF( nb_trccbc > 0 ) THEN
+           if (lwp) write(numout,'(a,i5,a,i10)') '   reading CBC data for ', nb_trccbc ,' variable(s) at step ', kt
+           CALL fld_read(kt=kt, kn_fsbc=1, sd=sf_trccbc)
+         ENDIF
+
       ENDIF
 
@@ -303,4 +446,10 @@
    !!----------------------------------------------------------------------
 CONTAINS
+
+   SUBROUTINE trc_bc_init( ntrc )        ! Empty routine
+      INTEGER,INTENT(IN) :: ntrc                           ! number of tracers
+      WRITE(*,*) 'trc_bc_init: You should not have seen this print! error?', kt
+   END SUBROUTINE trc_bc_init
+
    SUBROUTINE trc_bc_read( kt )        ! Empty routine
       WRITE(*,*) 'trc_bc_read: You should not have seen this print! error?', kt
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcdta.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcdta.F90	(revision 10155)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcdta.F90	(revision 10156)
@@ -9,4 +9,5 @@
    !!            3.4   !  2010-11  (C. Ethe, G. Madec)  use of fldread + dynamical allocation 
    !!            3.5   !  2013-08  (M. Vichi)  generalization for other BGC models
+   !!            3.6   !  2015-03  (T. Lovato) revision of code log info
    !!----------------------------------------------------------------------
 #if  defined key_top 
@@ -72,4 +73,10 @@
       IF( nn_timing == 1 )  CALL timing_start('trc_dta_init')
       !
+      IF( lwp ) THEN
+         WRITE(numout,*) ' '
+         WRITE(numout,*) '  trc_dta_init : Tracers Initial Conditions (IC)'
+         WRITE(numout,*) '  ~~~~~~~~~~~ '
+      ENDIF
+      !
       !  Initialisation
       ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0  
@@ -77,5 +84,5 @@
       ALLOCATE( n_trc_index(ntrc), slf_i(ntrc), STAT=ierr0 )
       IF( ierr0 > 0 ) THEN
-         CALL ctl_stop( 'trc_nam: unable to allocate n_trc_index' )   ;   RETURN
+         CALL ctl_stop( 'trc_dta_init: unable to allocate n_trc_index' )   ;   RETURN
       ENDIF
       nb_trcdta      = 0
@@ -97,9 +104,9 @@
       REWIND( numnat_ref )              ! Namelist namtrc_dta in reference namelist : Passive tracer input data
       READ  ( numnat_ref, namtrc_dta, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_dta in reference namelist', lwp )
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_dta_init in reference namelist', lwp )
 
       REWIND( numnat_cfg )              ! Namelist namtrc_dta in configuration namelist : Passive tracer input data
       READ  ( numnat_cfg, namtrc_dta, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_dta in configuration namelist', lwp )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_dta_init in configuration namelist', lwp )
       IF(lwm) WRITE ( numont, namtrc_dta )
 
@@ -109,12 +116,14 @@
                clndta = TRIM( sn_trcdta(jn)%clvar ) 
                clntrc = TRIM( ctrcnm   (jn)       ) 
+               if (jn > jptra) clntrc='Dummy' ! By pass weird formats in ocean.output if ntrc > jptra
                zfact  = rn_trfac(jn)
                IF( clndta /=  clntrc ) THEN 
-                  CALL ctl_warn( 'trc_dta_init: passive tracer data initialisation :  ',   &
-                  &              'the variable name in the data file : '//clndta//   & 
-                  &              '  must be the same than the name of the passive tracer : '//clntrc//' ')
+                  CALL ctl_warn( 'trc_dta_init: passive tracer data initialisation    ',   &
+                  &              'Input name of data file : '//TRIM(clndta)//   &
+                  &              ' differs from that of tracer : '//TRIM(clntrc)//' ')
                ENDIF
-               WRITE(numout,*) ' read an initial file for passive tracer number :', jn, ' name : ', clndta, & 
-               &               ' multiplicative factor : ', zfact
+               WRITE(numout,*) ' '
+               WRITE(numout,'(a, i3,3a,e11.3)') ' Read IC file for tracer number :', &
+               &            jn, ', name : ', TRIM(clndta), ', Multiplicative Scaling factor : ', zfact
             ENDIF
          END DO
@@ -124,5 +133,5 @@
          ALLOCATE( sf_trcdta(nb_trcdta), rf_trfac(nb_trcdta), STAT=ierr1 )
          IF( ierr1 > 0 ) THEN
-            CALL ctl_stop( 'trc_dta_ini: unable to allocate  sf_trcdta structure' )   ;   RETURN
+            CALL ctl_stop( 'trc_dta_init: unable to allocate  sf_trcdta structure' )   ;   RETURN
          ENDIF
          !
@@ -135,5 +144,5 @@
                IF( sn_trcdta(jn)%ln_tint )  ALLOCATE( sf_trcdta(jl)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 )
                IF( ierr2 + ierr3 > 0 ) THEN
-                 CALL ctl_stop( 'trc_dta : unable to allocate passive tracer data arrays' )   ;   RETURN
+                 CALL ctl_stop( 'trc_dta_init : unable to allocate passive tracer data arrays' )   ;   RETURN
                ENDIF
             ENDIF
@@ -141,5 +150,5 @@
          ENDDO
          !                         ! fill sf_trcdta with slf_i and control print
-         CALL fld_fill( sf_trcdta, slf_i, cn_dir, 'trc_dta', 'Passive tracer data', 'namtrc' )
+         CALL fld_fill( sf_trcdta, slf_i, cn_dir, 'trc_dta_init', 'Passive tracer data', 'namtrc' )
          !
       ENDIF
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcini.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 10155)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 10156)
@@ -24,4 +24,8 @@
    USE trcini_c14b     ! C14 bomb initialisation
    USE trcini_my_trc   ! MY_TRC   initialisation
+   ! +++>>> FABM
+   USE trcsms_fabm     ! FABM initialisation
+   USE trcini_fabm     ! FABM initialisation
+   ! FABM <<<FABM
    USE trcdta          ! initialisation from files
    USE daymod          ! calendar manager
@@ -32,4 +36,5 @@
    USE sbc_oce
    USE trcice          ! tracers in sea ice
+   USE trcbc,   only : trc_bc_init ! generalized Boundary Conditions
  
    IMPLICIT NONE
@@ -69,4 +74,8 @@
       IF(lwp) WRITE(numout,*) 'trc_init : initial set up of the passive tracers'
       IF(lwp) WRITE(numout,*) '~~~~~~~'
+      ! +++>>> FABM
+      ! Allow FABM to update numbers of biogeochemical tracers, diagnostics (jptra etc.)
+      IF( lk_fabm ) CALL nemo_fabm_init
+      ! FABM <<<+++
 
       CALL top_alloc()              ! allocate TOP arrays
@@ -101,4 +110,7 @@
       IF( lk_c14b    )       CALL trc_ini_c14b         ! C14 bomb  tracer
       IF( lk_my_trc  )       CALL trc_ini_my_trc       ! MY_TRC  tracers
+      ! +++>>> FABM
+      IF( lk_fabm    )       CALL trc_ini_fabm         ! FABM    tracers
+      ! FABM <<<+++
 
       CALL trc_ice_ini                                 ! Tracers in sea ice
@@ -110,4 +122,5 @@
       ENDIF
 
+      ! Initialisation of tracers Initial Conditions
       IF( ln_trcdta )      CALL trc_dta_init(jptra)
 
@@ -144,4 +157,17 @@
         ! 
       ENDIF
+      ! --->>> FABM
+! Initialisation of tracers Boundary Conditions  - here so that you can use initial condition as boundary
+      !IF( lk_my_trc )     CALL trc_bc_init(jptra)
+      ! FABM <<<---
+      ! FABM +++>>>
+! Initialisation of FABM diagnostics and tracer boundary conditions (so that you can use initial condition as boundary)
+      IF( lk_fabm )     THEN
+          wndm=0._wp !uninitiased field at this point
+          qsr=0._wp !uninitiased field at this point
+          CALL compute_fabm ! only needed to set-up diagnostics
+          CALL trc_bc_init(jptra)
+      ENDIF
+      ! FABM <<<+++
  
       tra(:,:,:,:) = 0._wp
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcnam.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcnam.F90	(revision 10155)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcnam.F90	(revision 10156)
@@ -25,4 +25,7 @@
    USE trcnam_c14b       ! C14 SMS namelist
    USE trcnam_my_trc     ! MY_TRC SMS namelist
+   ! +++>>> FABM
+   USE trcnam_fabm       ! FABM SMS namelist
+   ! FABM <<<+++
    USE trd_oce       
    USE trdtrc_oce
@@ -34,4 +37,10 @@
    PUBLIC trc_nam_run  ! called in trcini
    PUBLIC trc_nam      ! called in trcini
+   PUBLIC trc_nam_dia
+#if defined key_trdmxl_trc  || defined key_trdtrc
+   NAMELIST/namtrc_trd/ nn_trd_trc, nn_ctls_trc, rn_ucf_trc, &
+      &                ln_trdmxl_trc_restart, ln_trdmxl_trc_instant, &
+      &                cn_trdrst_trc_in, cn_trdrst_trc_out, ln_trdtrc
+#endif
 
    !! * Substitutions
@@ -57,4 +66,8 @@
       !!---------------------------------------------------------------------
       INTEGER  ::   jn                  ! dummy loop indice
+#if defined key_trdmxl_trc  || defined key_trdtrc
+      INTEGER :: ios
+#endif
+
       !                                        !   Parameters of the run 
       IF( .NOT. lk_offline ) CALL trc_nam_run
@@ -168,4 +181,10 @@
       ELSE                    ;   IF(lwp) WRITE(numout,*) '          MY_TRC not used'
       ENDIF
+
+      ! +++>>> FABM
+      IF( lk_fabm    ) THEN   ;   CALL trc_nam_fabm        ! FABM tracers
+      ELSE                    ;   IF(lwp) WRITE(numout,*) '          FABM not used'
+      ENDIF
+      ! FABM <<<+++
       !
    END SUBROUTINE trc_nam
@@ -187,5 +206,5 @@
 
 
-      IF(lwp) WRITE(numout,*) 'trc_nam : read the passive tracer namelists'
+      IF(lwp) WRITE(numout,*) 'trc_nam_run : read the passive tracer namelists'
       IF(lwp) WRITE(numout,*) '~~~~~~~'
 
@@ -234,5 +253,10 @@
 
       ! --- Namelist declarations --- !
-      TYPE(TRC_I_NML), DIMENSION(jptra) :: sn_tri_tracer
+      ! --->>> FABM
+      !TYPE(TRC_I_NML), DIMENSION(jptra) :: sn_tri_tracer
+      ! FABM <<<--- 
+      ! +++>>> FABM
+      TYPE(TRC_I_NML), DIMENSION(jpmaxtrc) :: sn_tri_tracer
+      ! FABM <<<+++ 
       NAMELIST/namtrc_ice/ nn_ice_tr, sn_tri_tracer
 
@@ -278,7 +302,12 @@
       !!
       !!---------------------------------------------------------------------
-      TYPE(PTRACER), DIMENSION(jptra) :: sn_tracer  ! type of tracer for saving if not key_iomput
-      !!
-      NAMELIST/namtrc/ sn_tracer, ln_trcdta,ln_trcdmp, ln_trcdmp_clo
+      ! --->>> FABM
+      !TYPE(PTRACER), DIMENSION(jptra) :: sn_tracer  ! type of tracer for saving if not key_iomput
+      ! FABM <<<---
+      ! +++>>> FABM
+      TYPE(PTRACER), DIMENSION(jpmaxtrc) :: sn_tracer  ! type of tracer for saving if not key_iomput
+      ! FABM <<<+++
+      !!
+      NAMELIST/namtrc/ sn_tracer, ln_trcdta, ln_trcdmp, ln_trcdmp_clo
   
       INTEGER  ::   ios                 ! Local integer output status for namelist read
@@ -286,7 +315,15 @@
       !!---------------------------------------------------------------------
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'trc_nam : read the passive tracer namelists'
+      IF(lwp) WRITE(numout,*) 'trc_nam_trc : read the passive tracer namelists'
       IF(lwp) WRITE(numout,*) '~~~~~~~'
 
+      ! Initialise logical flags to .FALSE.:
+      sn_tracer(:)%llinit = .FALSE.
+      sn_tracer(:)%llsave = .FALSE.
+#ifdef key_fabm
+      sn_tracer(:)%llsbc = .FALSE.
+      sn_tracer(:)%llcbc = .FALSE.
+      sn_tracer(:)%llcbc = .FALSE.
+#endif
 
       REWIND( numnat_ref )              ! Namelist namtrc in reference namelist : Passive tracer variables
@@ -304,7 +341,20 @@
          ctrcun    (jn) = TRIM( sn_tracer(jn)%clunit  )
          ln_trc_ini(jn) =       sn_tracer(jn)%llinit
+! --->>> FABM
+!!#if defined key_my_trc
+! FABM <<<---
+! +++>>> FABM
+#if defined key_fabm
+! FABM <<<+++
+         ln_trc_sbc(jn) =       sn_tracer(jn)%llsbc
+         ln_trc_cbc(jn) =       sn_tracer(jn)%llcbc
+         ln_trc_obc(jn) =       sn_tracer(jn)%llobc
+#endif
          ln_trc_wri(jn) =       sn_tracer(jn)%llsave
       END DO
-      
+     
+      ! +++>>> FABM
+      if (lk_fabm) CALL trc_nam_fabm_override
+      ! FABM <<<+++
     END SUBROUTINE trc_nam_trc
 
@@ -322,7 +372,7 @@
       INTEGER ::  ierr
 #if defined key_trdmxl_trc  || defined key_trdtrc
-      NAMELIST/namtrc_trd/ nn_trd_trc, nn_ctls_trc, rn_ucf_trc, &
-         &                ln_trdmxl_trc_restart, ln_trdmxl_trc_instant, &
-         &                cn_trdrst_trc_in, cn_trdrst_trc_out, ln_trdtrc
+!     NAMELIST/namtrc_trd/ nn_trd_trc, nn_ctls_trc, rn_ucf_trc, &
+!        &                ln_trdmxl_trc_restart, ln_trdmxl_trc_instant, &
+!        &                cn_trdrst_trc_in, cn_trdrst_trc_out, ln_trdtrc
 #endif
       NAMELIST/namtrc_dia/ ln_diatrc, ln_diabio, nn_writedia, nn_writebio
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcrst.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 10155)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 10156)
@@ -28,4 +28,7 @@
    USE iom
    USE daymod
+   ! +++>>> FABM
+   USE trcrst_fabm
+   ! FABM <<<+++
    IMPLICIT NONE
    PRIVATE
@@ -117,4 +120,8 @@
          CALL iom_get( numrtr, jpdom_autoglo, 'TRB'//ctrcnm(jn), trb(:,:,:,jn) )
       END DO
+      ! +++>>> FABM
+
+      IF (lk_fabm) CALL trc_rst_read_fabm
+      ! FABM <<<+++
       !
    END SUBROUTINE trc_rst_read
@@ -142,4 +149,7 @@
          CALL iom_rstput( kt, nitrst, numrtw, 'TRB'//ctrcnm(jn), trb(:,:,:,jn) )
       END DO
+      ! +++>>> FABM
+      IF (lk_fabm) CALL trc_rst_wri_fabm(kt)
+      ! FABM <<<+++
       !
       IF( kt == nitrst ) THEN
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcsms.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcsms.F90	(revision 10155)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcsms.F90	(revision 10156)
@@ -19,4 +19,7 @@
    USE trcsms_c14b        ! C14b tracer 
    USE trcsms_my_trc      ! MY_TRC  tracers
+   ! +++>>>> FABM
+   USE trcsms_fabm        ! FABM tracers
+   ! FABM <<<+++
    USE prtctl_trc         ! Print control for debbuging
 
@@ -52,4 +55,7 @@
       IF( lk_c14b    )   CALL trc_sms_c14b   ( kt )    ! surface fluxes of C14
       IF( lk_my_trc  )   CALL trc_sms_my_trc ( kt )    ! MY_TRC  tracers
+      ! +++>>> FABM
+      IF( lk_fabm    )   CALL trc_sms_fabm ( kt )      ! FABM tracers
+      ! FABM <<<+++
 
       IF(ln_ctl) THEN      ! print mean trends (used for debugging)
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcstp.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcstp.F90	(revision 10155)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcstp.F90	(revision 10156)
@@ -100,4 +100,8 @@
          IF( lrst_trc )            CALL trc_rst_wri  ( kt )       ! write tracer restart file
          IF( lk_trdmxl_trc  )      CALL trd_mxl_trc  ( kt )       ! trends: Mixed-layer
+#if defined key_tracer_budget
+!slwa tracer budget
+         IF( lk_iomput ) CALL trc_wri (kt, 2) 
+#endif
          !
          IF( nn_dttrc /= 1   )     CALL trc_sub_reset( kt )       ! resetting physical variables when sub-stepping
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcwri.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcwri.F90	(revision 10155)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_FABM/NEMOGCM/NEMO/TOP_SRC/trcwri.F90	(revision 10156)
@@ -21,4 +21,7 @@
    USE trcwri_c14b
    USE trcwri_my_trc
+   ! +++>>> FABM
+   USE trcwri_fabm
+   ! FABM <<<+++
 
    IMPLICIT NONE
@@ -32,5 +35,9 @@
 CONTAINS
 
+#if defined key_tracer_budget
+   SUBROUTINE trc_wri( kt , fl)  !slwa
+#else
    SUBROUTINE trc_wri( kt )
+#endif
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE trc_wri  ***
@@ -39,4 +46,9 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT( in )     :: kt
+      ! +++>>>FABM
+#if defined key_tracer_budget
+      INTEGER, INTENT( in ), OPTIONAL     :: fl  ! slwa
+#endif
+      ! FABM <<<+++
       !
       INTEGER                   :: jn
@@ -59,6 +71,20 @@
       IF( lk_cfc     )   CALL trc_wri_cfc        ! surface fluxes of CFC
       IF( lk_c14b    )   CALL trc_wri_c14b       ! surface fluxes of C14
-      IF( lk_my_trc  )   CALL trc_wri_my_trc     ! MY_TRC  tracers
+      ! +++>>>FABM
+#if defined key_tracer_budget
+      IF( PRESENT(fl) ) THEN
+         IF( lk_fabm    )   CALL trc_wri_fabm (kt, fl) ! MY_TRC  tracers for budget
+         IF( lk_my_trc ) CALL trc_wri_my_trc (kt, fl)    ! MY_TRC  tracers for budget
+      ELSE
+         IF( lk_fabm    )   CALL trc_wri_fabm (kt) ! FABM  tracers for budget
+         IF( lk_my_trc  )   CALL trc_wri_my_trc (kt) ! MY_TRC  tracers
+      ENDIF
+#else
+      IF( lk_fabm  )   CALL trc_wri_fabm (kt)     ! FABM  tracers
+      IF( lk_my_trc  )   CALL trc_wri_my_trc(kt)     ! MY_TRC  tracers
+#endif
+      ! FABM <<<+++
       !
+
       IF( nn_timing == 1 )  CALL timing_stop('trc_wri')
       !
