Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90	(revision 13175)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90	(revision 13241)
@@ -1082,5 +1082,5 @@
 #elif defined key_fabm
       totalk_tavg(:,:,:)  = totalk_tavg(:,:,:) + &
-         &                  fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta) / pnumtimes_tavg
+         &                  model%get_interior_diagnostic_data(jp_fabm_o3ta) / pnumtimes_tavg
       totalk_tavg(:,:,:)  = totalk_tavg(:,:,:) * tmask(:,:,:)
 #endif
@@ -1127,5 +1127,5 @@
          cchl_p_tavg(:,:,:)     = cchl_p(:,:,:)
 #elif defined key_fabm
-         totalk_tavg(:,:,:)  = fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta)
+         totalk_tavg(:,:,:)  = model%get_interior_diagnostic_data(jp_fabm_o3ta)
          totalk_tavg(:,:,:)  = totalk_tavg(:,:,:) * tmask(:,:,:)
 #endif
Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90	(revision 13175)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90	(revision 13241)
@@ -25,6 +25,4 @@
    USE par_fabm
    USE st2d_fabm, ONLY: fabm_st2dn
-   USE fabm, ONLY: fabm_get_interior_diagnostic_data, &
-      &            fabm_get_horizontal_diagnostic_data
 #endif
 
@@ -211,5 +209,5 @@
       END DO
       DO jn = 1, jp_fabm_3d
-         fabm_3d_25h(:,:,:,jn) = fabm_get_interior_diagnostic_data(model, jn)
+         fabm_3d_25h(:,:,:,jn) = model%get_interior_diagnostic_data(jn)
       END DO
       DO jn = 1, jp_fabm_surface
@@ -220,5 +218,5 @@
       END DO
       DO jn = 1, jp_fabm_2d
-         fabm_2d_25h(:,:,jn) = fabm_get_horizontal_diagnostic_data(model,jn)
+         fabm_2d_25h(:,:,jn) = model%get_horizontal_diagnostic_data(jn)
       END DO
 #endif
@@ -327,5 +325,5 @@
       END DO
       DO jn = 1, jp_fabm_3d
-         fabm_3d_25h(:,:,:,jn) = fabm_3d_25h(:,:,:,jn) + fabm_get_interior_diagnostic_data(model, jn)
+         fabm_3d_25h(:,:,:,jn) = fabm_3d_25h(:,:,:,jn) + model%get_interior_diagnostic_data(jn)
       END DO
       DO jn = 1, jp_fabm_surface
@@ -336,5 +334,5 @@
       END DO
       DO jn = 1, jp_fabm_2d
-         fabm_2d_25h(:,:,jn) = fabm_2d_25h(:,:,jn) + fabm_get_horizontal_diagnostic_data(model,jn)
+         fabm_2d_25h(:,:,jn) = fabm_2d_25h(:,:,jn) + model%get_horizontal_diagnostic_data(jn)
       END DO
 #endif
@@ -401,9 +399,9 @@
             DO jn = 1, jp_fabm
                zw3d(:,:,:) = fabm_25h(:,:,:,jn)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
-               CALL iom_put( TRIM(model%state_variables(jn)%name)//"25h", zw3d  )
+               CALL iom_put( TRIM(model%interior_state_variables(jn)%name)//"25h", zw3d  )
             END DO
             DO jn = 1, jp_fabm_3d
                zw3d(:,:,:) = fabm_3d_25h(:,:,:,jn)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
-               CALL iom_put( TRIM(model%diagnostic_variables(jn)%name)//"25h", zw3d  )
+               CALL iom_put( TRIM(model%interior_diagnostic_variables(jn)%name)//"25h", zw3d  )
             END DO
             DO jn = 1, jp_fabm_surface
@@ -468,5 +466,5 @@
       END DO
       DO jn = 1, jp_fabm_3d
-         fabm_3d_25h(:,:,:,jn) = fabm_get_interior_diagnostic_data(model, jn)
+         fabm_3d_25h(:,:,:,jn) = model%get_interior_diagnostic_data(jn)
       END DO
       DO jn = 1, jp_fabm_surface
@@ -477,5 +475,5 @@
       END DO
       DO jn = 1, jp_fabm_2d
-         fabm_2d_25h(:,:,jn) = fabm_get_horizontal_diagnostic_data(model,jn)
+         fabm_2d_25h(:,:,jn) = model%get_horizontal_diagnostic_data(jn)
       END DO
 #endif
Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90	(revision 13175)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90	(revision 13241)
@@ -14,5 +14,4 @@
    USE trc, ONLY: trn
    USE par_fabm
-   USE fabm, ONLY: fabm_get_interior_diagnostic_data
 #endif
 
@@ -172,13 +171,13 @@
          DO jn = 1, jp_fabm
             CALL dia_calctmb( trn(:,:,:,jp_fabm_m1+jn), zwtmb )
-            CALL iom_put( "top_"//TRIM(model%state_variables(jn)%name) , zwtmb(:,:,1) )
-            CALL iom_put( "mid_"//TRIM(model%state_variables(jn)%name) , zwtmb(:,:,2) )
-            CALL iom_put( "bot_"//TRIM(model%state_variables(jn)%name) , zwtmb(:,:,3) )
+            CALL iom_put( "top_"//TRIM(model%interior_state_variables(jn)%name) , zwtmb(:,:,1) )
+            CALL iom_put( "mid_"//TRIM(model%interior_state_variables(jn)%name) , zwtmb(:,:,2) )
+            CALL iom_put( "bot_"//TRIM(model%interior_state_variables(jn)%name) , zwtmb(:,:,3) )
          END DO
          DO jn = 1, jp_fabm_3d
-            CALL dia_calctmb( fabm_get_interior_diagnostic_data(model, jn), zwtmb )
-            CALL iom_put( "top_"//TRIM(model%diagnostic_variables(jn)%name) , zwtmb(:,:,1) )
-            CALL iom_put( "mid_"//TRIM(model%diagnostic_variables(jn)%name) , zwtmb(:,:,2) )
-            CALL iom_put( "bot_"//TRIM(model%diagnostic_variables(jn)%name) , zwtmb(:,:,3) )
+            CALL dia_calctmb( model%get_interior_diagnostic_data(jn), zwtmb )
+            CALL iom_put( "top_"//TRIM(model%interior_diagnostic_variables(jn)%name) , zwtmb(:,:,1) )
+            CALL iom_put( "mid_"//TRIM(model%interior_diagnostic_variables(jn)%name) , zwtmb(:,:,2) )
+            CALL iom_put( "bot_"//TRIM(model%interior_diagnostic_variables(jn)%name) , zwtmb(:,:,3) )
          END DO
 #endif
Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 13175)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 13241)
@@ -45,4 +45,5 @@
    LOGICAL , PUBLIC ::   ln_qsr_bio   !: bio-model      light absorption flag
    LOGICAL , PUBLIC ::   ln_qsr_ice   !: light penetration for ice-model LIM3 (clem)
+   LOGICAL , PUBLIC ::   ln_qsr_spec  !: spectral model heating from ERSEM
    INTEGER , PUBLIC ::   nn_chldta    !: use Chlorophyll data (=1) or not (=0)
    INTEGER , PUBLIC ::   nn_kd490dta  !: use kd490dta data (=1) or not (=0)
@@ -150,5 +151,34 @@
       
       !                                           ! ============================================== !
-      IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN      !  bio-model fluxes  : all vertical coordinates  !
+      IF( ln_qsr_spec ) THEN                      !  ERSEM spectral heating                        !
+         !                                        ! ============================================== !
+         DO jk = 1, jpkm1
+           qsr_hc(:,:,jk) = r1_rau0_rcp * ( etot3(:,:,jk) )
+         END DO
+         !                                        Add to the general trend
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z1_e3t = zfact / fse3t(ji,jj,jk)
+                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) * z1_e3t
+               END DO
+            END DO
+         END DO
+         CALL iom_put( 'qsr3d', etot3 )   ! Shortwave Radiation 3D distribution
+         IF ( ln_qsr_ice ) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF ( qsr(ji,jj) /= 0._wp ) THEN
+                     fraqsr_1lev(ji,jj) = ( qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) ) )
+                  ELSE
+                     fraqsr_1lev(ji,jj) = 1.
+                  ENDIF
+               END DO
+            END DO
+         ENDIF
+         ! 
+
+         !                                        ! ============================================== !
+      ELSEIF( lk_qsr_bio .AND. ln_qsr_bio ) THEN      !  bio-model fluxes  : all vertical coordinates  !
          !                                        ! ============================================== !
          DO jk = 1, jpkm1
@@ -423,5 +453,5 @@
       !!
       NAMELIST/namtra_qsr/  sn_chl, sn_kd490, cn_dir, ln_traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice,  &
-         &                  nn_chldta, rn_abs, rn_si0, rn_si1, nn_kd490dta
+         &                  ln_qsr_spec, nn_chldta, rn_abs, rn_si0, rn_si1, nn_kd490dta
       !!----------------------------------------------------------------------
 
@@ -451,4 +481,5 @@
          WRITE(numout,*) '      2 band               light penetration   ln_qsr_2bd = ', ln_qsr_2bd
          WRITE(numout,*) '      bio-model            light penetration   ln_qsr_bio = ', ln_qsr_bio
+         WRITE(numout,*) '      ERSEM spectral heating model             ln_qsr_spec= ', ln_qsr_spec
          WRITE(numout,*) '      light penetration for ice-model LIM3     ln_qsr_ice = ', ln_qsr_ice
          WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)    nn_chldta  = ', nn_chldta
@@ -470,4 +501,5 @@
          IF( ln_qsr_2bd  )   ioptio = ioptio + 1
          IF( ln_qsr_bio  )   ioptio = ioptio + 1
+         IF( ln_qsr_spec )   ioptio = ioptio + 1
          IF( nn_kd490dta == 1 )   ioptio = ioptio + 1
          !
@@ -481,4 +513,5 @@
          IF( ln_qsr_bio                      )   nqsr =  4
          IF( nn_kd490dta == 1                )   nqsr =  5
+         IF( ln_qsr_spec                     )   nqsr =  6
          !
          IF(lwp) THEN                   ! Print the choice
@@ -489,5 +522,13 @@
             IF( nqsr ==  4 )   WRITE(numout,*) '         bio-model light penetration'
             IF( nqsr ==  5 )   WRITE(numout,*) '         KD490 light penetration'
-         ENDIF
+            IF( nqsr ==  6 )   WRITE(numout,*) '         ERSEM spectral light penetration'
+         ENDIF
+#if ! defined key_fabm
+         !
+         IF( nqsr ==  6 ) THEN
+            CALL ctl_stop( 'ln_qsr_spec=.true. so trying to use ERSEM spectral light penetration', &
+               &           'but not running with ERSEM' )
+         ENDIF
+#endif
          !
       ENDIF
Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/inputs_fabm.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/inputs_fabm.F90	(revision 13175)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/inputs_fabm.F90	(revision 13241)
@@ -19,7 +19,9 @@
    USE fldread
    USE par_fabm
-   USE fabm
+   USE fabm, only: type_fabm_horizontal_variable_id
 
    IMPLICIT NONE
+
+#  include "vectopt_loop_substitute.h90"
 
    PRIVATE
@@ -40,6 +42,6 @@
 
    TYPE, PUBLIC, EXTENDS(type_input_variable) :: type_input_data
-      TYPE(type_horizontal_variable_id)   :: horizontal_id
-      TYPE(type_input_data), POINTER   :: next => null()
+      TYPE(type_fabm_horizontal_variable_id) :: horizontal_id
+      TYPE(type_input_data), POINTER         :: next => null()
    END TYPE
    TYPE (type_input_data), POINTER, PUBLIC :: first_input_data => NULL()
@@ -87,6 +89,6 @@
            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
+           input_data%horizontal_id = model%get_horizontal_variable_id(name)
+           IF (.NOT.model%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.')
@@ -130,5 +132,5 @@
            ! within tracer field
            DO jn=1,jp_fabm
-             IF (TRIM(name) == TRIM(model%state_variables(jn)%name)) THEN
+             IF (TRIM(name) == TRIM(model%interior_state_variables(jn)%name)) THEN
                river_data%jp_pos = jp_fabm_m1+jn
              END IF
@@ -173,5 +175,5 @@
          ! 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)
+         CALL model%link_horizontal_data(input_data%horizontal_id,input_data%sf(1)%fnow(:,:,1),source=data_source_user)
          input_data => input_data%next
       END DO
@@ -226,6 +228,6 @@
 #endif
         IF( kt == nit000 .OR. ( kt /= nit000 ) ) THEN
-            DO jj = 1, jpj
-              DO ji = 1, jpi
+            DO jj = 2, jpjm1
+              DO ji = fs_2, fs_jpim1
                 ! convert units and divide by surface area
                 ! loading / cell volume * vertical fraction of riverload
Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/par_fabm.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/par_fabm.F90	(revision 13175)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/par_fabm.F90	(revision 13241)
@@ -1,9 +1,13 @@
 MODULE par_fabm
 
+#if defined key_fabm
+#  include "fabm_version.h"
+#  if _FABM_API_VERSION_ < 1
+#    error You need FABM 1.0 or later
+#  endif
    USE fabm
+#endif
 
    IMPLICIT NONE
-
-   TYPE (type_model) :: model !FABM model instance
 
    INTEGER, PUBLIC :: jp_fabm0, jp_fabm1, jp_fabm, &
@@ -37,10 +41,13 @@
                       jp_fabm_pgrow, jp_fabm_ploss
 
+   LOGICAL, PUBLIC, ALLOCATABLE, DIMENSION(:) ::   lk_rad_fabm !: FABM negativity correction flag array
+
 #if defined key_fabm
+   CLASS (type_fabm_model), POINTER :: model !FABM model instance
+
    !!---------------------------------------------------------------------
    !!   '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
    !!---------------------------------------------------------------------
Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcini_fabm.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcini_fabm.F90	(revision 13175)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcini_fabm.F90	(revision 13241)
@@ -4,5 +4,6 @@
    !! TOP :   initialisation of the FABM tracers
    !!======================================================================
-   !! History :   2.0  !  2007-12  (C. Ethe, G. Madec) Original code
+   !! History :   1.0  !  2015-04  (PML) Original code
+   !! History :   1.1  !  2020-06  (PML) Update to FABM 1.0, improved performance
    !!----------------------------------------------------------------------
 #if defined key_fabm
@@ -17,7 +18,7 @@
    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, &
+   USE fabm, only: fabm_create_model, type_fabm_variable
+   USE fabm_driver
+   USE inputs_fabm,ONLY: initialize_inputs, link_inputs, &
      type_input_variable,type_input_data,type_river_data, &
      first_input_data,first_river_data
@@ -33,5 +34,5 @@
 
 #if defined key_git_version
-!#include "gitversion.h90"
+#  include "gitversion.h90"
    CHARACTER(len=*),parameter :: git_commit_id = _NEMO_COMMIT_ID_
    CHARACTER(len=*),parameter :: git_branch_name = _NEMO_BRANCH_
@@ -39,5 +40,11 @@
 
    PUBLIC   trc_ini_fabm   ! called by trcini.F90 module
-   PUBLIC   nemo_fabm_init
+   PUBLIC   nemo_fabm_configure
+
+   TYPE,extends(type_base_driver) :: type_nemo_fabm_driver
+   contains
+      procedure :: fatal_error => nemo_fabm_driver_fatal_error
+      procedure :: log_message => nemo_fabm_driver_log_message
+   end type
 
    !!----------------------------------------------------------------------
@@ -48,5 +55,5 @@
 CONTAINS
 
-   SUBROUTINE nemo_fabm_init()
+   SUBROUTINE nemo_fabm_configure()
       INTEGER :: jn
       INTEGER, PARAMETER :: xml_unit = 1979
@@ -55,8 +62,10 @@
       CLASS (type_input_variable),POINTER :: input_pointer
 
+      ALLOCATE(type_nemo_fabm_driver::driver)
+
       ! 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)
+      model => fabm_create_model()
+
+      jp_fabm = size(model%interior_state_variables)
       jp_fabm_bottom = size(model%bottom_state_variables)
       jp_fabm_surface = size(model%surface_state_variables)
@@ -66,10 +75,12 @@
       jptra = jptra + jp_fabm
       jp_fabm_2d = size(model%horizontal_diagnostic_variables)
-      jp_fabm_3d = size(model%diagnostic_variables)
+      jp_fabm_3d = size(model%interior_diagnostic_variables)
       jpdia2d = jpdia2d + jp_fabm_2d
       jpdia3d = jpdia3d + jp_fabm_3d
       jpdiabio = jpdiabio + jp_fabm
 
-      !Initialize input data structures.
+      ! Read inputs (river and additional 2D forcing) from fabm_input.nml
+      ! This must be done before writing field_def_fabm.xml, as that file
+      ! also describes the additional input variables.
       call initialize_inputs
 
@@ -128,5 +139,5 @@
       IF (lwp) THEN
          ! write field_def_fabm.xml on lead process
-         OPEN(UNIT=xml_unit,FILE='field_def_fabm.xml',ACTION='WRITE',STATUS='REPLACE')
+         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" >'
@@ -134,10 +145,10 @@
          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))
+            CALL write_variable_xml(xml_unit,model%interior_state_variables(jn))
 #if defined key_trdtrc
-            CALL write_trends_xml(xml_unit,model%state_variables(jn))
-#endif
-            CALL write_25hourm_xml(xml_unit,model%state_variables(jn))
-            CALL write_tmb_xml(xml_unit,model%state_variables(jn))
+            CALL write_trends_xml(xml_unit,model%interior_state_variables(jn))
+#endif
+            CALL write_25hourm_xml(xml_unit,model%interior_state_variables(jn))
+            CALL write_tmb_xml(xml_unit,model%interior_state_variables(jn))
          END DO
          WRITE (xml_unit,1000) ' </field_group>'
@@ -155,12 +166,22 @@
 
          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)
-            CALL write_25hourm_xml(xml_unit,model%diagnostic_variables(jn),3)
-            CALL write_tmb_xml(xml_unit,model%diagnostic_variables(jn))
+         DO jn=1,size(model%interior_diagnostic_variables)
+            CALL write_variable_xml(xml_unit,model%interior_diagnostic_variables(jn),3)
+            CALL write_25hourm_xml(xml_unit,model%interior_diagnostic_variables(jn),3)
+            CALL write_tmb_xml(xml_unit,model%interior_diagnostic_variables(jn))
          END DO
          DO jn=1,size(model%horizontal_diagnostic_variables)
             CALL write_variable_xml(xml_unit,model%horizontal_diagnostic_variables(jn))
             CALL write_25hourm_xml(xml_unit,model%horizontal_diagnostic_variables(jn))
+         END DO
+         DO jn=1,size(model%interior_state_variables)
+            WRITE (xml_unit,'(A)') '  <field id="'//TRIM(model%interior_state_variables(jn)%name)// &
+               &                   '_VINT" long_name="depth-integrated '//TRIM(model%interior_state_variables(jn)%long_name)// &
+               &                   '" unit="'//TRIM(model%interior_state_variables(jn)%units)//'*m" default_value="0.0" />'
+         END DO
+         DO jn=1,size(model%interior_diagnostic_variables)
+            WRITE (xml_unit,'(A)') '  <field id="'//TRIM(model%interior_diagnostic_variables(jn)%name)// &
+               &                   '_VINT" long_name="depth-integrated '//TRIM(model%interior_diagnostic_variables(jn)%long_name)// &
+               &                   '" unit="'//TRIM(model%interior_diagnostic_variables(jn)%units)//'*m" default_value="0.0" />'
          END DO
          WRITE (xml_unit,1000) ' </field_group>'
@@ -195,10 +216,10 @@
 1000 FORMAT (A)
 
-   END SUBROUTINE nemo_fabm_init
+   END SUBROUTINE nemo_fabm_configure
 
    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
+      CLASS (type_fabm_variable),INTENT(IN) :: variable
 
       CHARACTER(LEN=20) :: missing_value,string_dimensions
@@ -233,5 +254,5 @@
       INTEGER,INTENT(IN) :: xml_unit
       INTEGER,INTENT(IN),OPTIONAL :: flag_grid_ref
-      CLASS (type_external_variable),INTENT(IN) :: variable
+      CLASS (type_fabm_variable),INTENT(IN) :: variable
 
       CHARACTER(LEN=20) :: missing_value,string_dimensions
@@ -265,5 +286,5 @@
    SUBROUTINE write_tmb_xml(xml_unit,variable)
       INTEGER,INTENT(IN) :: xml_unit
-      CLASS (type_external_variable),INTENT(IN) :: variable
+      CLASS (type_fabm_variable),INTENT(IN) :: variable
 
       CHARACTER(LEN=20) :: missing_value
@@ -279,5 +300,5 @@
       INTEGER,INTENT(IN) :: xml_unit
       INTEGER,INTENT(IN),OPTIONAL :: flag_grid_ref
-      CLASS (type_external_variable),INTENT(IN) :: variable
+      CLASS (type_fabm_variable),INTENT(IN) :: variable
 
       INTEGER :: number_dimensions,i
@@ -383,13 +404,10 @@
       !! ** Purpose :   initialization for FABM model
       !!
-      !! ** Method  : - Read the namcfc namelist and check the parameter values
+      !! ** Method  : - Allocate FABM arrays, configure domain, send data
       !!----------------------------------------------------------------------
 #if defined key_git_version
-      TYPE (type_version),POINTER :: 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,*)
@@ -399,10 +417,5 @@
       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)
@@ -411,25 +424,28 @@
 #endif
 
+      ! Allocate FABM arrays
+      IF(trc_sms_fabm_alloc() /= 0) CALL ctl_stop( 'STOP', 'trc_ini_fabm: unable to allocate FABM arrays' )
+
       ! Log mapping of FABM states:
       IF (lwp) THEN
-         IF (jp_fabm.gt.0) WRITE(numout,*) " FABM tracers:"
+         IF (jp_fabm > 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:"
+            WRITE(numout,*) "   State",jn,":",trim(model%interior_state_variables(jn)%name), &
+               " (",trim(model%interior_state_variables(jn)%long_name), &
+               ") [",trim(model%interior_state_variables(jn)%units),"]"
+         END DO
+         IF (jp_fabm_surface > 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:"
+         END DO
+         IF (jp_fabm_bottom > 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 DO
+      END IF
       
       ! Initialise variables required for Hemmings et al. (2008)
@@ -442,4 +458,19 @@
    END SUBROUTINE trc_ini_fabm
 
+   SUBROUTINE nemo_fabm_driver_fatal_error(self, location, message)
+      CLASS (type_nemo_fabm_driver), INTENT(INOUT) :: self
+      CHARACTER(len=*),              INTENT(IN)    :: location, message
+
+      CALL ctl_stop('STOP', TRIM(location)//': '//TRIM(message))
+      STOP
+   END SUBROUTINE
+
+   SUBROUTINE nemo_fabm_driver_log_message(self, message)
+      CLASS (type_nemo_fabm_driver), INTENT(INOUT) :: self
+      CHARACTER(len=*),              INTENT(IN)    :: message
+
+      IF(lwp) WRITE (numout,*) TRIM(message)
+   END SUBROUTINE
+
    INTEGER FUNCTION fabm_state_index( state_name )
       !!----------------------------------------------------------------------
@@ -461,5 +492,5 @@
       fabm_state_index = -1
       DO jn=1,jp_fabm
-         IF (TRIM(model%state_variables(jn)%name) == TRIM(state_name)) THEN
+         IF (TRIM(model%interior_state_variables(jn)%name) == TRIM(state_name)) THEN
             fabm_state_index = jn
             EXIT
@@ -492,6 +523,6 @@
       
       fabm_diag_index = -1
-      DO jn = 1, SIZE(model%diagnostic_variables)
-         IF (TRIM(model%diagnostic_variables(jn)%name) == TRIM(diag_name)) THEN
+      DO jn = 1, SIZE(model%interior_diagnostic_variables)
+         IF (TRIM(model%interior_diagnostic_variables(jn)%name) == TRIM(diag_name)) THEN
             fabm_diag_index = jn
             EXIT
Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcnam_fabm.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcnam_fabm.F90	(revision 13175)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcnam_fabm.F90	(revision 13241)
@@ -4,5 +4,6 @@
    !! TOP :   initialisation of some run parameters for FABM bio-model
    !!======================================================================
-   !! History :   2.0  !  2007-12  (C. Ethe, G. Madec) Original code
+   !! History :   1.0  !  2015-04  (PML) Original code
+   !! History :   1.1  !  2020-06  (PML) Update to FABM 1.0, improved performance
    !!----------------------------------------------------------------------
 #if defined key_fabm
@@ -18,5 +19,4 @@
    USE par_fabm
    USE trcsms_fabm
-
 
    IMPLICIT NONE
@@ -35,17 +35,43 @@
 
    SUBROUTINE trc_nam_fabm
+      LOGICAL :: l_ext
+      INTEGER :: nmlunit, ios
+      NAMELIST/namfabm/ nn_adv
+
+      ! Read NEMO-FABM coupler settings from namfabm
+      nn_adv = 3
+      INQUIRE( FILE='namelist_fabm_ref', EXIST=l_ext )
+      IF (l_ext) then
+         CALL ctl_opn( nmlunit, 'namelist_fabm_ref', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE.)
+         READ(nmlunit, nml=namfabm, iostat=ios)
+         IF( ios /= 0 ) CALL ctl_nam ( ios , 'namfabm in namelist_fabm_ref', .TRUE. )
+      END IF
+      INQUIRE( FILE='namelist_fabm_cfg', EXIST=l_ext )
+      IF (l_ext) then
+         CALL ctl_opn( nmlunit, 'namelist_fabm_cfg', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE.)
+         READ(nmlunit, nml=namfabm, iostat=ios)
+         IF( ios /= 0 ) CALL ctl_nam ( ios , 'namfabm in namelist_fabm_cfg', .TRUE. )
+      END IF
+      IF (nn_adv /= 1 .and. nn_adv /= 3) CALL ctl_stop( 'STOP', 'trc_ini_fabm: nn_adv must be 1 or 3.' )
    END SUBROUTINE trc_nam_fabm
 
-   SUBROUTINE trc_nam_fabm_override
+   SUBROUTINE trc_nam_fabm_override(sn_tracer)
+      TYPE(PTRACER), DIMENSION(jpmaxtrc), INTENT(INOUT) :: sn_tracer
+
       INTEGER :: jn
+      CHARACTER(LEN=3) :: index
 
       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.
+         IF (sn_tracer(jn)%clsname /= 'NONAME' .AND. sn_tracer(jn)%clsname /= model%interior_state_variables(jn)%name) THEN
+            WRITE (index,'(i0)') jn
+            CALL ctl_stop('Tracer name mismatch in namtrc: '//TRIM(sn_tracer(jn)%clsname)//' found at sn_tracer('//TRIM(index)//') where '//TRIM(model%interior_state_variables(jn)%name)//' was expected.')
+         END IF
+         sn_tracer(jn)%clsname = TRIM(model%interior_state_variables(jn)%name)
+         sn_tracer(jn)%cllname = TRIM(model%interior_state_variables(jn)%long_name)
+         sn_tracer(jn)%clunit = TRIM(model%interior_state_variables(jn)%units)
+         sn_tracer(jn)%llinit = .FALSE.
       END DO
    END SUBROUTINE trc_nam_fabm_override
-   
+
 #else
    !!----------------------------------------------------------------------
Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcrst_fabm.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcrst_fabm.F90	(revision 13175)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcrst_fabm.F90	(revision 13241)
@@ -4,5 +4,6 @@
    !! Read and write additional restart fields used by FABM
    !!======================================================================
-   !! History :
+   !! History :   1.0  !  2015-04  (PML) Original code
+   !! History :   1.1  !  2020-06  (PML) Update to FABM 1.0, improved performance
    !!----------------------------------------------------------------------
 #if defined key_fabm
Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90	(revision 13175)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90	(revision 13241)
@@ -5,4 +5,5 @@
    !!======================================================================
    !! History :   1.0  !  2015-04  (PML) Original code
+   !! History :   1.1  !  2020-06  (PML) Update to FABM 1.0, improved performance
    !!----------------------------------------------------------------------
 #if defined key_fabm
@@ -19,4 +20,5 @@
    USE trd_oce
    USE trdtrc
+   USE traqsr,only: ln_qsr_spec
 #if defined key_trdtrc && defined key_iomput
    USE trdtrc_oce
@@ -24,8 +26,8 @@
 
    USE oce, only: tsn  ! Needed?
-   USE sbc_oce, only: lk_oasis
+   USE sbc_oce, only: lk_oasis,fr_i
    USE dom_oce
    USE zdf_oce
-   !USE iom
+   USE iom
    USE xios
    USE cpl_oasis3
@@ -36,8 +38,7 @@
    USE asmbgc, ONLY: mld_choice_bgc
    USE lbclnk
+   USE fabm_types, ONLY: standard_variables
 
    !USE fldread         !  time interpolation
-
-   USE fabm
 
    IMPLICIT NONE
@@ -48,34 +49,29 @@
    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
+   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 ! called by trcwri_fabm.F90
+
+   ! Work arrays
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: flux          ! Cross-interface flux of pelagic variables (# m-2 s-1)
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: current_total ! Totals of conserved quantities
+
+   ! Arrays for environmental variables that are computed by the coupler
+   REAL(wp), PUBLIC, TARGET, ALLOCATABLE, DIMENSION(:,:,:) :: prn,rho
+   REAL(wp), PUBLIC, TARGET, ALLOCATABLE, DIMENSION(:,:) :: taubot
+   REAL(wp), PUBLIC, TARGET :: daynumber_in_year
+   REAL(wp), PUBLIC, POINTER :: swrad(:,:,:)  ! pointer to ERSEM spectral heating term
+
+   ! State repair counters
+   INTEGER, SAVE :: repair_interior_count = 0
+   INTEGER, SAVE :: repair_surface_count  = 0
+   INTEGER, SAVE :: repair_bottom_count   = 0
+
+   ! Coupler parameters
+   INTEGER, PUBLIC :: nn_adv  ! Vertical advection scheme for sinking/floating/movement
+                              ! (1: 1st order upwind, 3: 3rd order TVD)
+
+   ! Flag indicating whether model%start has been called (will be done on-demand)
+   LOGICAL, SAVE :: started = .false.
 
    !!----------------------------------------------------------------------
@@ -96,6 +92,7 @@
       !
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      INTEGER :: jn
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrfabm
+      INTEGER :: jn, jk
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrfabm, pdat
+      REAL(wp), DIMENSION(jpi,jpj)    :: vint
 
 !!----------------------------------------------------------------------
@@ -109,9 +106,11 @@
       IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
 
+      IF (.NOT. started) CALL nemo_fabm_start
+
       CALL update_inputs( kt )
 
-      CALL compute_fabm
-
-      CALL compute_vertical_movement( kt )
+      CALL compute_fabm( kt )
+
+      CALL compute_vertical_movement( kt, nn_adv )
 
       CALL st2d_fabm_nxt( kt )
@@ -123,4 +122,28 @@
       CALL trc_bc_read  ( kt )       ! tracers: surface and lateral Boundary Conditions
       CALL trc_rnf_fabm ( kt ) ! River forcings
+
+      ! Send 3D diagnostics to output (these apply to time "n")
+      DO jn = 1, size(model%interior_diagnostic_variables)
+         IF (model%interior_diagnostic_variables(jn)%save) THEN
+            ! Save 3D field
+            pdat => model%get_interior_diagnostic_data(jn)
+            CALL iom_put(model%interior_diagnostic_variables(jn)%name, pdat)
+
+            ! Save depth integral if selected for output in XIOS
+            IF (iom_use(TRIM(model%interior_diagnostic_variables(jn)%name)//'_VINT')) THEN
+               vint = 0._wp
+               DO jk = 1, jpkm1
+                  vint = vint + pdat(:,:,jk) * fse3t(:,:,jk) * tmask(:,:,jk)
+               END DO
+               CALL iom_put(TRIM(model%interior_diagnostic_variables(jn)%name)//'_VINT', vint)
+            END IF
+         END IF
+      END DO
+
+      ! Send 2D diagnostics to output (these apply to time "n")
+      DO jn = 1, size(model%horizontal_diagnostic_variables)
+         IF (model%horizontal_diagnostic_variables(jn)%save) &
+             CALL iom_put( model%horizontal_diagnostic_variables(jn)%name, model%get_horizontal_diagnostic_data(jn))
+      END DO
 
       IF( l_trdtrc ) THEN      ! Save the trends in the mixed layer
@@ -139,5 +162,6 @@
       INTEGER, INTENT(IN) :: kt
       INTEGER :: ji,jj,jk,jkmax
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: pgrow_3d, ploss_3d, zmld
+      REAL(wp), DIMENSION(jpi,jpj) :: zmld
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: pgrow_3d, ploss_3d
 	   
       IF (kt == nittrc000) THEN
@@ -148,6 +172,6 @@
       PHYT_AVG(:,:)  = 0.0
 	     
-      pgrow_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_pgrow)
-      ploss_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_ploss)
+      pgrow_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_pgrow)
+      ploss_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_ploss)
       
       SELECT CASE( mld_choice_bgc )
@@ -220,12 +244,14 @@
    END SUBROUTINE asmdiags_fabm
 
-   SUBROUTINE compute_fabm()
+   SUBROUTINE compute_fabm( kt )
+      INTEGER, INTENT(in) :: kt   ! ocean time-step index
+
       INTEGER :: ji,jj,jk,jn
-      TYPE(type_state) :: valid_state
+      LOGICAL :: valid, repaired
       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
+      CALL check_state(.TRUE., valid, repaired)
+      IF (.NOT. valid) THEN
          WRITE(numout,*) "Invalid value in FABM encountered in area ",narea,"!!!"
 #if defined key_iomput
@@ -240,5 +266,5 @@
 #endif
       END IF
-      IF (valid_state%repaired) THEN
+      IF (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
@@ -246,49 +272,44 @@
       ENDIF
 
+      daynumber_in_year = fjulday - fjulstartyear + 1
+
       ! 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
+      IF (ALLOCATED(rho)) rho = rau0 * ( 1._wp + rhd )
+
+      IF (ALLOCATED(prn)) THEN
+         zalfg = 0.5e-4_wp * grav ! FABM wants dbar, convert from Pa (and multiply with 0.5 to average 2 cell thicknesses below)
+         prn(:,:,1) = 10.1325_wp + zalfg * fse3t(:,:,1) * rho(:,:,1)
+         DO jk = 2, jpkm1                                              ! Vertical integration from the surface
+            prn(:,:,jk) = prn(:,:,jk-1) + zalfg * ( &
+                        fse3t(:,:,jk-1) * rho(:,:,jk-1)  &
+                        + fse3t(:,:,jk) * rho(:,:,jk) )
+         END DO
+      END IF
+
+      ! Compute the bottom stress
+      ! copied from diawri.F90
+      ! ------------------------------------
+
+      IF (ALLOCATED(taubot)) THEN
+         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)
+                  !
+            END DO
+         END DO
+      END IF
+
+      CALL model%prepare_inputs(real(kt, wp),nyear,nmonth,nday,REAL(nsec_day,wp))
+
+      ! Retrieve 3D shortwave and store in etot3
+      IF (ln_qsr_spec) etot3(:,:,:) = swrad(:,:,:)
 
       ! Zero rate array of interface-attached state variables
@@ -296,87 +317,96 @@
 
       ! Compute interfacial source terms and fluxes
-      DO jj=1,jpj
-         ! Process bottom (fabm_do_bottom increments rather than sets, so zero flux array first)
+      DO jj=2,jpjm1
+         ! Process bottom (get_bottom_sources 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:))
+         CALL model%get_bottom_sources(fs_2,fs_jpim1,jj,flux,fabm_st2Da(fs_2:fs_jpim1,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)
+            ! Divide bottom fluxes by height of bottom layer and add to source terms.
+            DO ji=fs_2,fs_jpim1
+               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 (get_surface_sources 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))
+         CALL model%get_surface_sources(fs_2,fs_jpim1,jj,flux,fabm_st2Da(fs_2:fs_jpim1,jj,1:jp_fabm_surface))
+         ! Divide surface fluxes by height of surface layer and add to source terms.
          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
+            DO ji=fs_2,fs_jpim1
+               tra(ji,jj,1,jp_fabm_m1+jn) = tra(ji,jj,1,jp_fabm_m1+jn) + flux(ji,jn)/fse3t(ji,jj,1)
+            END DO
+         END DO
+      END DO
+
+      ! Compute interior source terms (NB get_interior_sources increments rather than sets)
+      DO jk=1,jpkm1
+         DO jj=2,jpjm1
+            CALL model%get_interior_sources(fs_2,fs_jpim1,jj,jk,tra(fs_2:fs_jpim1,jj,jk,jp_fabm0:jp_fabm1))
+         END DO
+      END DO
+
+      CALL model%finalize_outputs()
    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
+   SUBROUTINE check_state(repair, valid, repaired)
+      LOGICAL, INTENT(IN)  :: repair
+      LOGICAL, INTENT(OUT) :: valid, repaired
+
+      INTEGER :: jj, jk
+      LOGICAL :: valid_int, valid_sf, valid_bt
+
+      valid = .TRUE.     ! Whether the model state is valid after this subroutine returns
+      repaired = .FALSE. ! Whether the model state has been repaired by this subroutine
+      DO jk=1,jpkm1
+         DO jj=2,jpjm1
+            CALL model%check_interior_state(fs_2, fs_jpim1, jj, jk, repair, valid_int)
+            IF (repair .AND. .NOT. valid_int) THEN
                repair_interior_count = repair_interior_count + 1
-               exit_state%repaired = .TRUE.
+               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
+            IF (.NOT. (valid_int .OR. repair)) valid = .FALSE.
+         END DO
+      END DO
+      DO jj=2,jpjm1
+         CALL model%check_surface_state(fs_2, fs_jpim1, jj, repair, valid_sf)
+         IF (repair .AND. .NOT. valid_sf) THEN
             repair_surface_count = repair_surface_count + 1
-            exit_state%repaired = .TRUE.
+            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
+         IF (.NOT. (valid_sf .AND. valid_bt) .AND. .NOT. repair) valid = .FALSE.
+         CALL model%check_bottom_state(fs_2, fs_jpim1, jj, repair, valid_bt)
+         IF (repair .AND. .NOT. valid_bt) THEN
             repair_bottom_count = repair_bottom_count + 1
-            exit_state%repaired = .TRUE.
+            repaired = .TRUE.
          END IF
-         IF (.NOT.(valid_sf.AND.valid_bt).AND..NOT.repair) exit_state%valid = .FALSE.
-      END DO
-   END FUNCTION
+         IF (.NOT. (valid_sf .AND. valid_bt) .AND. .NOT. repair) valid = .FALSE.
+      END DO
+   END SUBROUTINE
 
    SUBROUTINE trc_sms_fabm_check_mass()
       REAL(wp) :: total(SIZE(model%conserved_quantities))
-      INTEGER :: jk,jj,jn
+      INTEGER :: ji,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)
+      IF (.NOT. started) CALL nemo_fabm_start
+
+      DO jk=1,jpkm1
+         DO jj=2,jpjm1
+            CALL model%get_interior_conserved_quantities(fs_2,fs_jpim1,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 ji=fs_2,fs_jpim1
+                  total(jn) = total(jn) + cvol(ji,jj,jk) * current_total(ji,jn) * tmask_i(ji,jj)
+               END DO
+            END DO
+         END DO
+      END DO
+
+      DO jj=2,jpjm1
+         CALL model%get_horizontal_conserved_quantities(fs_2,fs_jpim1,jj,current_total)
          DO jn=1,SIZE(model%conserved_quantities)
-            total(jn) = total(jn) + SUM(e1e2t(:,jj)*current_total(:,jn)*tmask_i(:,jj))
+            DO ji=fs_2,fs_jpim1
+               total(jn) = total(jn) + e1e2t(ji,jj) * current_total(ji,jn) * tmask_i(ji,jj)
+            END DO
          END DO
       END DO
@@ -434,5 +464,5 @@
 
    INTEGER FUNCTION trc_sms_fabm_alloc()
-      INTEGER :: jj,jk,jn
+      INTEGER :: jn
       !!----------------------------------------------------------------------
       !!              ***  ROUTINE trc_sms_fabm_alloc  ***
@@ -441,7 +471,7 @@
       ! 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))
+      IF (model%variable_needs_values(fabm_standard_variables%pressure)) ALLOCATE(prn(jpi, jpj, jpk))
+      IF (ALLOCATED(prn) .or. model%variable_needs_values(fabm_standard_variables%density)) ALLOCATE(rho(jpi, jpj, jpk))
+      IF (model%variable_needs_values(fabm_standard_variables%bottom_stress)) ALLOCATE(taubot(jpi, jpj))
       ! ALLOCATE( tab(...) , STAT=trc_sms_fabm_alloc )
 
@@ -452,16 +482,9 @@
 
       ! 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(flux(fs_2:fs_jpim1,jp_fabm))
 
       ! 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)))
+      ALLOCATE(w_ct(fs_2:fs_jpim1,1:jpkm1,jp_fabm))
+      ALLOCATE(current_total(fs_2:fs_jpim1,SIZE(model%conserved_quantities)))
 #if defined key_trdtrc && defined key_iomput
       IF( lk_trdtrc ) ALLOCATE(tr_vmv(jpi,jpj,jpk,jp_fabm))
@@ -474,91 +497,58 @@
       !
 
-      ! 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
+      ! Provide FABM with domain extents
+      CALL model%set_domain(jpi, jpj, jpk, rdt)
+      CALL model%set_domain_start(fs_2, 2, 1)
+      CALL model%set_domain_stop(fs_jpim1, jpjm1, jpkm1)
+
+      ! Provide FABM with the vertical indices of the bottom, and the land-sea mask.
+      CALL model%set_bottom_index(mbkt)  ! NB mbkt extents should match dimension lengths provided to set_domain
+      CALL model%set_mask(tmask,tmask(:,:,1)) ! NB tmask extents should match dimension lengths provided to set_domain
+
+      ! Initialize state and send pointers to state data to FABM
+      ! We mask land points in states with zeros, as per with NEMO "convention"
+      ! NB we cannot call model%initialize_*_state at this point, because model%start has not been called yet.
+      DO jn=1,jp_fabm
+         trn(:,:,:,jp_fabm_m1+jn) = model%interior_state_variables(jn)%initial_value * tmask
+         CALL model%link_interior_state_data(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))
+         fabm_st2Dn(:,:,jn) = model%surface_state_variables(jn)%initial_value * tmask(:,:,1)
+         CALL model%link_surface_state_data(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))
+         fabm_st2Dn(:,:,jp_fabm_surface+jn) = model%bottom_state_variables(jn)%initial_value * tmask(:,:,1)
+         CALL model%link_bottom_state_data(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)
+      CALL model%link_interior_data(fabm_standard_variables%depth, fsdept(:,:,:))
+      CALL model%link_interior_data(fabm_standard_variables%temperature, tsn(:,:,:,jp_tem))
+      CALL model%link_interior_data(fabm_standard_variables%practical_salinity, tsn(:,:,:,jp_sal))
+      IF (ALLOCATED(rho)) CALL model%link_interior_data(fabm_standard_variables%density, rho(:,:,:))
+      IF (ALLOCATED(prn)) CALL model%link_interior_data(fabm_standard_variables%pressure, prn)
+      IF (ALLOCATED(taubot)) CALL model%link_horizontal_data(fabm_standard_variables%bottom_stress, taubot(:,:))
+      CALL model%link_interior_data(fabm_standard_variables%cell_thickness, fse3t(:,:,:))
+      CALL model%link_horizontal_data(fabm_standard_variables%latitude, gphit)
+      CALL model%link_horizontal_data(fabm_standard_variables%longitude, glamt)
+      CALL model%link_scalar(fabm_standard_variables%number_of_days_since_start_of_the_year, daynumber_in_year)
+      CALL model%link_horizontal_data(fabm_standard_variables%wind_speed, wndm(:,:))
+      CALL model%link_horizontal_data(fabm_standard_variables%surface_downwelling_shortwave_flux, qsr(:,:))
+      CALL model%link_horizontal_data(fabm_standard_variables%bottom_depth_below_geoid, bathy(:,:))
+      CALL model%link_horizontal_data(fabm_standard_variables%ice_area_fraction, fr_i(:,:))
 
       ! 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
+      CALL link_inputs
+      CALL update_inputs(nit000, .FALSE.)
 
       ! Set mask for negativity corrections to the relevant states
-      lk_rad_fabm = .FALSE.
+      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.'
+        IF (model%interior_state_variables(jn)%minimum >= 0._wp) THEN
+          lk_rad_fabm(jn) = .TRUE.
+          IF(lwp) WRITE(numout,*) 'FABM clipping for '//TRIM(model%interior_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
@@ -566,10 +556,29 @@
       fabm_st2Db = fabm_st2Dn
 
-      ! Initialise repair counters
-      repair_interior_count = 0
-      repair_surface_count = 0
-      repair_bottom_count = 0
+      ! Pointer to spectral heating term
+      swrad => model%get_data(model%get_interior_variable_id(standard_variables%net_rate_of_absorption_of_shortwave_energy_in_layer))
+
 
    END FUNCTION trc_sms_fabm_alloc
+
+   SUBROUTINE nemo_fabm_start()
+      INTEGER :: jn
+
+      ! Make FABM aware of diagnostics that are not needed [not included in output]
+      ! This works only after iom has completely initialised, because it depends on iom_use
+      DO jn=1,size(model%interior_diagnostic_variables)
+         model%interior_diagnostic_variables(jn)%save = iom_use(model%interior_diagnostic_variables(jn)%name) &
+            .or. iom_use(TRIM(model%interior_diagnostic_variables(jn)%name)//'_VINT')
+      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
+
+      ! Check whether FABM has all required data
+      ! [after this, the save attribute of diagnostic variables can no longer change!]
+      CALL model%start()
+
+      started = .TRUE.
+   END SUBROUTINE
 
 #else
Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcwri_fabm.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcwri_fabm.F90	(revision 13175)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcwri_fabm.F90	(revision 13241)
@@ -4,5 +4,6 @@
    !!    fabm :   Output of FABM tracers
    !!======================================================================
-   !! History :   1.0  !  2009-05 (C. Ethe)  Original code
+   !! History :   1.0  !  2015-04  (PML) Original code
+   !! History :   1.1  !  2020-06  (PML) Update to FABM 1.0, improved performance
    !!----------------------------------------------------------------------
 #if defined key_top && key_fabm && defined key_iomput
@@ -18,5 +19,4 @@
    USE par_fabm
    USE st2d_fabm
-   USE fabm, only: fabm_get_bulk_diagnostic_data, fabm_get_horizontal_diagnostic_data
 
    IMPLICIT NONE
@@ -31,5 +31,4 @@
        MODULE PROCEDURE wri_fabm,wri_fabm_fl
    END INTERFACE trc_wri_fabm
-
 
    PUBLIC trc_wri_fabm 
@@ -58,9 +57,9 @@
 ! 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(:,:,:) + &
+      DO jn = 1, jp_fabm
+        IF(ln_trdtrc (jp_fabm_m1+jn))THEN
+         trpool(:,:,:) = 0.5 * ( trn(:,:,:,jp_fabm_m1+jn)*fse3t_a(:,:,:) + &
                              tr_temp(:,:,:,jn)*fse3t(:,:,:) )
-         cltra = TRIM( model%state_variables(jn)%name )//"_e3t"     ! depth integrated output
+         cltra = TRIM( model%interior_state_variables(jn)%name )//"_e3t"     ! depth integrated output
          IF( kt == nittrc000 ) write(6,*)'output pool ',cltra
          CALL iom_put( cltra, trpool)
@@ -80,5 +79,6 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT( in )               :: kt
-      INTEGER              :: jn
+      INTEGER              :: jn, jk
+      REAL(wp), DIMENSION(jpi,jpj)    :: vint
 
 #if defined key_tracer_budget
@@ -90,5 +90,15 @@
 #endif
       DO jn = 1, jp_fabm
-         CALL iom_put( model%state_variables(jn)%name, trn(:,:,:,jp_fabm0+jn-1) )
+         ! Save 3D field
+         CALL iom_put(model%interior_state_variables(jn)%name, trn(:,:,:,jp_fabm_m1+jn))
+
+         ! Save depth integral if selected for output in XIOS
+         IF (iom_use(TRIM(model%interior_state_variables(jn)%name)//'_VINT')) THEN
+            vint = 0._wp
+            DO jk = 1, jpkm1
+               vint = vint + trn(:,:,jk,jp_fabm_m1+jn) * fse3t(:,:,jk) * tmask(:,:,jk)
+            END DO
+            CALL iom_put(TRIM(model%interior_state_variables(jn)%name)//'_VINT', vint)
+         END IF
       END DO
       DO jn = 1, jp_fabm_surface
@@ -99,18 +109,4 @@
       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
 
@@ -121,8 +117,21 @@
    !!  Dummy module :                                     No passive tracer
    !!----------------------------------------------------------------------
+   INTERFACE trc_wri_fabm
+       MODULE PROCEDURE wri_fabm,wri_fabm_fl
+   END INTERFACE trc_wri_fabm
+
    PUBLIC trc_wri_fabm
-CONTAINS
-   SUBROUTINE trc_wri_fabm                     ! Empty routine  
-   END SUBROUTINE trc_wri_fabm
+
+   CONTAINS
+
+   SUBROUTINE wri_fabm_fl(kt,fl)
+      INTEGER, INTENT( in )               :: fl
+      INTEGER, INTENT( in )               :: kt
+   END SUBROUTINE wri_fabm_fl
+
+   SUBROUTINE wri_fabm(kt)                 ! Empty routine  
+      INTEGER, INTENT( in )               :: kt
+   END SUBROUTINE wri_fabm
+
 #endif
 
Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90	(revision 13175)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90	(revision 13241)
@@ -15,5 +15,4 @@
    USE trc
    USE par_fabm
-   USE fabm
    USE dom_oce
 #if defined key_trdtrc && defined key_iomput
@@ -25,4 +24,5 @@
 
 #  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
 
    PRIVATE
@@ -32,7 +32,4 @@
    ! 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
@@ -41,23 +38,20 @@
    CONTAINS
 
-   SUBROUTINE compute_vertical_movement( kt )
+   SUBROUTINE compute_vertical_movement( kt, method )
       !!----------------------------------------------------------------------
       !!                     ***  compute_vertical_movement  ***
       !!
-      !! ** Purpose :   compute vertical movement of FABM tracers
+      !! ** Purpose : compute vertical movement of FABM tracers through the water
+      !!              (sinking/floating/active movement)
       !!
-      !! ** 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.
+      !! ** Method  : Retrieves additional vertical velocity field and applies
+      !!              advection scheme.
       !!----------------------------------------------------------------------
       !
-      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
+      INTEGER, INTENT(in) ::   kt     ! ocean time-step index
+      INTEGER, INTENT(in) ::   method ! advection method (1: 1st order upstream, 3: 3rd order TVD with QUICKEST limiter)
+
+      INTEGER :: ji,jj,jk,jn,k_floor
+      REAL(wp) :: zwgt_if(1:jpkm1-1), dc(1:jpkm1), w_if(1:jpkm1-1), z2dt, h(1:jpkm1)
 #if defined key_trdtrc
       CHARACTER (len=20) :: cltra
@@ -69,144 +63,53 @@
 
       IF( neuler == 0 .AND. kt == nittrc000 ) THEN
-          z2dt = rdt                  ! set time step size (Euler)
+         z2dt = rdt                  ! set time step size (Euler)
       ELSE
-          z2dt = 2._wp * rdt          ! set time step size (Leapfrog)
+         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,:))
+      DO jj=2,jpjm1 ! j-loop
+         ! Get vertical velocities at layer centres (entire i-k slice).
+         DO jk=1,jpkm1
+            CALL model%get_vertical_movement(fs_2,fs_jpim1,jj,jk,w_ct(:,jk,:))
          END DO
-
-         DO ji=1,jpi ! i-loop
+         DO ji=fs_2,fs_jpim1 ! 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)
+            k_floor = mbkt(ji,jj)
+            IF (k_floor > 1) THEN
                ! 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
+               !    - interface k sits between cell centre k and k+1 (k=0 for surface)
+               !    - k [1,jpkm1] 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
+               h(1:k_floor) = fse3t(ji,jj,1:k_floor)
+               zwgt_if(1:k_floor-1) = h(2:k_floor) / (h(1:k_floor-1) + h(2:k_floor))
 
                ! 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,jj,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
+                  IF (ALL(w_ct(ji,1:k_floor,jn) == 0._wp)) CYCLE
+
+                  ! Compute velocities at interfaces
+                  w_if(1:k_floor-1) = zwgt_if(1:k_floor-1) * w_ct(ji,1:k_floor-1,jn) + (1._wp - zwgt_if(1:k_floor-1)) * w_ct(ji,2:k_floor,jn)
+
+                  ! Compute change (per volume) due to vertical movement per layer
+                  IF (method == 1) THEN
+                     CALL advect_1(k_floor, trn(ji,jj,1:k_floor,jp_fabm_m1+jn), w_if(1:k_floor-1), h(1:k_floor), z2dt, dc(1:k_floor))
+                  ELSE
+                     CALL advect_3(k_floor, trb(ji,jj,1:k_floor,jp_fabm_m1+jn), w_if(1:k_floor-1), h(1:k_floor), z2dt, dc(1:k_floor))
                   END IF
-#endif
-               ENDDO ! State loop
+
+                  ! Incorporate change due to vertical movement in sources-sinks
+                  tra(ji,jj,1:k_floor,jp_fabm_m1+jn) = tra(ji,jj,1:k_floor,jp_fabm_m1+jn) + dc(1:k_floor)
+
+#if defined key_trdtrc && defined key_iomput
+                  ! Store change due to vertical movement as diagnostic
+                  IF( lk_trdtrc .AND. ln_trdtrc( jp_fabm_m1+jn)) tr_vmv(ji,jj,1:k_floor,jn) = dc(1:k_floor)
+#endif
+               END DO ! 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
@@ -220,4 +123,158 @@
    END SUBROUTINE compute_vertical_movement
 
+   SUBROUTINE advect_1(nk, c, w, h, dt, trend)
+      INTEGER,  INTENT(IN)  :: nk
+      REAL(wp), INTENT(IN)  :: c(1:nk)
+      REAL(wp), INTENT(IN)  :: w(1:nk-1)
+      REAL(wp), INTENT(IN)  :: h(1:nk)
+      REAL(wp), INTENT(IN)  :: dt
+      REAL(wp), INTENT(OUT) :: trend(1:nk)
+
+      REAL(wp) :: flux(0:nk)
+      INTEGER  :: jk
+      ! Compute fluxes (per surface area) over at interfaces (remember: positive for upwards)
+      flux(0) = 0._wp
+      DO jk=1,nk-1 ! k-loop
+         IF (w(jk) > 0) THEN
+            ! Upward movement (source layer is jk+1)
+            flux(jk) = min(w(jk), h(jk+1)/dt) * c(jk+1)
+         ELSE
+            ! Downward movement (source layer is jk)
+            flux(jk) = max(w(jk), -h(jk)/dt) * c(jk)
+         END IF
+      END DO
+      flux(nk) = 0._wp
+      trend = (flux(1:nk) - flux(0:nk-1)) / h
+   END SUBROUTINE
+
+   SUBROUTINE advect_3(nk, c_old, w, h, dt, trend)
+      INTEGER,  INTENT(IN)  :: nk
+      REAL(wp), INTENT(IN)  :: c_old(1:nk)
+      REAL(wp), INTENT(IN)  :: w(1:nk-1)
+      REAL(wp), INTENT(IN)  :: h(1:nk)
+      REAL(wp), INTENT(IN)  :: dt
+      REAL(wp), INTENT(OUT) :: trend(1:nk)
+
+      INTEGER, PARAMETER :: n_itermax=100
+      REAL(wp) :: cmax_no
+      REAL(wp) :: cfl(1:nk-1)
+      INTEGER  :: n_iter, n_count, jk
+      REAL(wp) :: c(1:nk)
+      REAL(wp) :: tr_u(1:nk-1)
+      REAL(wp) :: tr_c(1:nk-1)
+      REAL(wp) :: tr_d(1:nk-1)
+      REAL(wp) :: delta_tr_u(1:nk-1)
+      REAL(wp) :: delta_tr(1:nk-1)
+      REAL(wp) :: ratio(1:nk-1)
+      REAL(wp) :: x_fac(1:nk-1)
+      REAL(wp) :: phi_lim(1:nk-1)
+      REAL(wp) :: limiter(1:nk-1)
+      REAL(wp) :: flux_if(1:nk-1)
+
+      c(:) = c_old(:)
+
+      ! get maximum Courant number:
+      cfl = ABS(w) * dt / (0.5_wp * (h(2:nk) + h(1:nk-1)))
+      cmax_no = MAXVAL(cfl)
+
+      ! number of iterations:
+      n_iter = MIN(n_itermax, INT(cmax_no) + 1)
+      IF (ln_ctl.AND.(n_iter .gt. 1)) THEN
+         WRITE(numout,*) 'compute_vertical_movement::advect_3():'
+         WRITE(numout,*) '   Maximum Courant number is ',cmax_no,'.'
+         WRITE(numout,*) '   ',n_iter,' iterations used for vertical advection.'
+      ENDIF
+
+      ! effective Courant number:
+      cfl = cfl/n_iter
+
+      DO n_count=1,n_iter ! Iterative loop
+         ! Determine tracer concentration at 1.5 upstream (tr_u), 0.5 upstream (tr_c), 0.5 downstream (tr_d) from interface
+         IF (nk.gt.2) THEN
+            ! More than 2 vertical wet points
+            IF (nk.gt.3) THEN
+               WHERE (w(2:nk-2).ge.0._wp)
+                  !upward movement
+                  tr_u(2:nk-2)=c(4:nk)
+                  tr_c(2:nk-2)=c(3:nk-1)
+                  tr_d(2:nk-2)=c(2:nk-2)
+               ELSEWHERE
+                  ! downward movement
+                  tr_u(2:nk-2)=c(1:nk-3)
+                  tr_c(2:nk-2)=c(2:nk-2)
+                  tr_d(2:nk-2)=c(3:nk-1)
+               ENDWHERE
+            ENDIF
+
+            ! Interface between surface layer and the next
+            IF (w(1).ge.0._wp) THEN
+               ! upward movement
+               tr_u(1)=c(3)
+               tr_c(1)=c(2)
+               tr_d(1)=c(1)
+            ELSE
+               ! downward movement
+               tr_u(1)=c(1)
+               tr_c(1)=c(1)
+               tr_d(1)=c(2)
+            ENDIF
+
+            ! Interface between bottom layer and the previous
+            IF (w(nk-1).ge.0._wp) THEN
+               ! upward movement
+               tr_u(nk-1)=c(nk)
+               tr_c(nk-1)=c(nk)
+               tr_d(nk-1)=c(nk-1)
+            ELSE
+               ! downward movement
+               tr_u(nk-1)=c(nk-2)
+               tr_c(nk-1)=c(nk-1)
+               tr_d(nk-1)=c(nk)
+            ENDIF
+         ELSE
+            ! only 2 vertical wet points, i.e. only 1 interface
+            IF (w(1).ge.0._wp) THEN
+               ! upward movement
+               tr_u(1)=c(2)
+               tr_c(1)=c(2)
+               tr_d(1)=c(1)
+            ELSE
+               ! downward movement
+               tr_u(1)=c(1)
+               tr_c(1)=c(1)
+               tr_d(1)=c(2)
+            ENDIF
+         ENDIF
+
+         delta_tr_u = tr_c - tr_u
+         delta_tr = tr_d - tr_c
+         WHERE (delta_tr * delta_tr_u > 0._wp)
+            ! Monotonic function over tr_u, tr_c, r_d
+
+            ! Compute slope ratio
+            ratio = delta_tr_u / delta_tr
+
+            ! QUICKEST flux limiter
+            x_fac = (1._wp - 2._wp * cfl) / 6._wp
+            phi_lim = (0.5_wp + x_fac) + (0.5_wp - x_fac) * ratio
+            limiter = MIN(phi_lim, 2._wp / (1._wp - cfl), 2._wp * ratio / (cfl + 1.e-10_wp))
+
+            ! Compute limited flux
+            flux_if = w * (tr_c + 0.5_wp * limiter * (1._wp - cfl) * delta_tr)
+         ELSEWHERE
+            ! Non-monotonic, use 1st order upstream
+            flux_if = w * tr_c
+         ENDWHERE
+
+         ! Compute pseudo update for trend aggregation:
+         c(1:nk-1) = c(1:nk-1) + dt / real(n_iter, wp) / h(1:nk-1) * flux_if
+         c(2:nk)   = c(2:nk)   - dt / real(n_iter, wp) / h(2:nk)   * flux_if
+
+      ENDDO ! Iterative loop
+
+      ! Estimate rate of change from pseudo state updates (source splitting):
+      trend = (c - c_old) / dt
+   END SUBROUTINE
+
 #endif
 END MODULE
Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/trcini.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 13175)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 13241)
@@ -76,5 +76,5 @@
       ! +++>>> FABM
       ! Allow FABM to update numbers of biogeochemical tracers, diagnostics (jptra etc.)
-      IF( lk_fabm ) CALL nemo_fabm_init
+      IF( lk_fabm ) CALL nemo_fabm_configure
       ! FABM <<<+++
 
@@ -123,5 +123,4 @@
       ! Initialisation of tracers Initial Conditions
       IF( ln_trcdta )      CALL trc_dta_init(jptra)
-
 
       IF( ln_rsttr ) THEN
@@ -162,12 +161,7 @@
       ! 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
+      IF( lk_fabm )     CALL trc_bc_init(jptra) 
       ! FABM <<<+++
- 
+
       tra(:,:,:,:) = 0._wp
       IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav )   &              ! Partial steps: before horizontal gradient of passive
Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/trcnam.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/trcnam.F90	(revision 13175)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/trcnam.F90	(revision 13241)
@@ -325,4 +325,5 @@
       sn_tracer(:)%llcbc = .FALSE.
       sn_tracer(:)%llcbc = .FALSE.
+      sn_tracer(:)%clsname = 'NONAME'
 #endif
 
@@ -335,4 +336,8 @@
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc in configuration namelist', lwp )
       IF(lwm) WRITE ( numont, namtrc )
+
+      ! +++>>> FABM
+      if (lk_fabm) CALL trc_nam_fabm_override(sn_tracer)
+      ! FABM <<<+++
 
       DO jn = 1, jptra
@@ -354,7 +359,4 @@
       END DO
      
-      ! +++>>> FABM
-      if (lk_fabm) CALL trc_nam_fabm_override
-      ! FABM <<<+++
     END SUBROUTINE trc_nam_trc
 
@@ -380,8 +382,4 @@
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       !!---------------------------------------------------------------------
-
-      IF(lwp) WRITE(numout,*) 
-      IF(lwp) WRITE(numout,*) 'trc_nam_dia : read the passive tracer diagnostics options'
-      IF(lwp) WRITE(numout,*) '~~~~~~~'
 
       IF(lwp) WRITE(numout,*)
