Index: /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90
===================================================================
--- /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90	(revision 11137)
+++ /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90	(revision 11137)
@@ -0,0 +1,502 @@
+MODULE dia25h 
+   !!======================================================================
+   !!                       ***  MODULE  diaharm  ***
+   !! Harmonic analysis of tidal constituents 
+   !!======================================================================
+   !! History :  3.6  !  2014  (E O'Dea)  Original code
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain
+   USE diainsitutem, ONLY: rinsitu_t, theta2t
+   USE in_out_manager  ! I/O units
+   USE iom             ! I/0 library
+   USE wrk_nemo        ! working arrays
+#if defined key_zdftke 
+   USE zdftke, ONLY: en
+#endif
+   USE zdf_oce, ONLY: avt, avm
+#if defined key_zdfgls
+   USE zdfgls, ONLY: mxln
+   USE zdf_oce, ONLY: en
+#endif
+   USE diatmb
+#if defined key_fabm
+   USE trc, ONLY: trn, visib
+   USE par_fabm
+   USE st2d_fabm, ONLY: fabm_st2dn
+   USE fabm, ONLY: fabm_get_bulk_diagnostic_data, &
+      &            fabm_get_horizontal_diagnostic_data
+#endif
+
+   IMPLICIT NONE
+   PRIVATE
+
+   LOGICAL , PUBLIC ::   ln_dia25h     !:  25h mean output
+   PUBLIC   dia_25h_init               ! routine called by nemogcm.F90
+   PUBLIC   dia_25h                    ! routine called by diawri.F90
+
+  !! * variables for calculating 25-hourly means
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   tn_25h  , sn_25h, rinsitu_t_25h  
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:)   ::   sshn_25h, insitu_bot_25h, temp_bot_25h 
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   un_25h  , vn_25h  , wn_25h
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   avt_25h , avm_25h
+#if defined key_zdfgls || key_zdftke
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   en_25h
+#endif
+#if defined key_zdfgls 
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   rmxln_25h
+#endif
+#if defined key_fabm
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:,:) :: fabm_25h
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:,:) :: fabm_3d_25h
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:)   :: fabm_surface_25h
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:)   :: fabm_bottom_25h
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:)   :: fabm_2d_25h
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:)   :: fabm_visib_25h
+#endif
+   INTEGER, SAVE :: cnt_25h     ! Counter for 25 hour means
+
+
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.6 , NEMO Consortium (2014)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dia_25h_init 
+      !!---------------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_25h_init  ***
+      !!     
+      !! ** Purpose: Initialization of 25h mean namelist 
+      !!        
+      !! ** Method : Read namelist
+      !!   History
+      !!   3.6  !  08-14  (E. O'Dea) Routine to initialize dia_25h
+      !!---------------------------------------------------------------------------
+      !!
+      INTEGER ::   ios                 ! Local integer output status for namelist read
+      INTEGER ::   ierror              ! Local integer for memory allocation
+      INTEGER ::   jn                  ! Loop counter
+      REAL(wp), DIMENSION(jpi,jpj,3)   ::   zwtmb                                 ! temporary workspace
+      !
+      NAMELIST/nam_dia25h/ ln_dia25h
+      !!----------------------------------------------------------------------
+      !
+      REWIND ( numnam_ref )              ! Read Namelist nam_dia25h in reference namelist : 25hour mean diagnostics
+      READ   ( numnam_ref, nam_dia25h, IOSTAT=ios, ERR= 901 )
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_dia25h in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist nam_dia25h in configuration namelist  25hour diagnostics
+      READ  ( numnam_cfg, nam_dia25h, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_dia25h in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, nam_dia25h )
+
+      IF(lwp) THEN                   ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'dia_25h_init : Output 25 hour mean diagnostics'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) 'Namelist nam_dia25h : set 25h outputs '
+         WRITE(numout,*) 'Switch for 25h diagnostics (T) or not (F)  ln_dia25h  = ', ln_dia25h
+      ENDIF
+      IF( .NOT. ln_dia25h )   RETURN
+      ! ------------------- !
+      ! 1 - Allocate memory !
+      ! ------------------- !
+      ALLOCATE( tn_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate tn_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( sn_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate sn_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( rinsitu_t_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate rinsitu_t_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( insitu_bot_25h(jpi,jpj), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate insitu_bot_25h' )   ;   RETURN
+      ENDIF      
+      ALLOCATE( temp_bot_25h(jpi,jpj), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate temp_bot_25h' )   ;   RETURN
+      ENDIF                           
+      ALLOCATE( un_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate un_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( vn_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate vn_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( wn_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate wn_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( avt_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate avt_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( avm_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate avm_25h' )   ;   RETURN
+      ENDIF
+# if defined key_zdfgls || defined key_zdftke
+      ALLOCATE( en_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate en_25h' )   ;   RETURN
+      ENDIF
+#endif
+# if defined key_zdfgls 
+      ALLOCATE( rmxln_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate rmxln_25h' )   ;   RETURN
+      ENDIF
+#endif
+      ALLOCATE( sshn_25h(jpi,jpj), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate sshn_25h' )   ;   RETURN
+      ENDIF
+#if defined key_fabm
+      ALLOCATE( fabm_25h(jpi,jpj,jpk,jp_fabm), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate fabm_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( fabm_3d_25h(jpi,jpj,jpk,jp_fabm_3d), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate fabm_3d_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( fabm_surface_25h(jpi,jpj,jp_fabm_surface), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate fabm_surface_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( fabm_bottom_25h(jpi,jpj,jp_fabm_bottom), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate fabm_bottom_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( fabm_2d_25h(jpi,jpj,jp_fabm_2d), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate fabm_2d_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( fabm_visib_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate fabm_visib_25h' )   ;   RETURN
+      ENDIF
+#endif 
+      ! ------------------------- !
+      ! 2 - Assign Initial Values !
+      ! ------------------------- !
+      cnt_25h = 1  ! sets the first value of sum at timestep 1 (note - should strictly be at timestep zero so before values used where possible) 
+      tn_25h(:,:,:) = tsb(:,:,:,jp_tem)
+      sn_25h(:,:,:) = tsb(:,:,:,jp_sal)
+      CALL theta2t
+      rinsitu_t_25h(:,:,:) = rinsitu_t(:,:,:)
+      CALL dia_calctmb( rinsitu_t(:,:,:),zwtmb )
+      insitu_bot_25h(:,:) = zwtmb(:,:,3)
+      CALL dia_calctmb( tn_25h(:,:,:),zwtmb )
+      temp_bot_25h(:,:) = zwtmb(:,:,3)
+      sshn_25h(:,:) = sshb(:,:)
+      un_25h(:,:,:) = ub(:,:,:)
+      vn_25h(:,:,:) = vb(:,:,:)
+      wn_25h(:,:,:) = wn(:,:,:)
+      avt_25h(:,:,:) = avt(:,:,:)
+      avm_25h(:,:,:) = avm(:,:,:)
+# if defined key_zdfgls || defined key_zdftke
+         en_25h(:,:,:) = en(:,:,:)
+#endif
+# if defined key_zdfgls
+         rmxln_25h(:,:,:) = mxln(:,:,:)
+#endif
+#if defined key_fabm
+      DO jn = 1, jp_fabm
+         fabm_25h(:,:,:,jn) = trn(:,:,:,jp_fabm_m1+jn)
+      END DO
+      DO jn = 1, jp_fabm_3d
+         fabm_3d_25h(:,:,:,jn) = fabm_get_bulk_diagnostic_data(model, jn)
+      END DO
+      DO jn = 1, jp_fabm_surface
+         fabm_surface_25h(:,:,jn) = fabm_st2dn(:,:,jn)
+      END DO
+      DO jn = 1, jp_fabm_bottom
+         fabm_bottom_25h(:,:,jn) = fabm_st2dn(:,:,jp_fabm_surface+jn)
+      END DO
+      DO jn = 1, jp_fabm_2d
+         fabm_2d_25h(:,:,jn) = fabm_get_horizontal_diagnostic_data(model,jn)
+      END DO
+      fabm_visib_25h(:,:,:) = visib(:,:,:)
+#endif
+#if defined key_lim3 || defined key_lim2
+         CALL ctl_stop('STOP', 'dia_25h not setup yet to do tidemean ice')
+#endif 
+
+      ! -------------------------- !
+      ! 3 - Return to dia_wri      !
+      ! -------------------------- !
+
+
+   END SUBROUTINE dia_25h_init
+
+
+   SUBROUTINE dia_25h( kt )  
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE dia_25h  ***
+      !!         
+      !!
+      !!--------------------------------------------------------------------
+      !!                   
+      !! ** Purpose :   Write diagnostics with M2/S2 tide removed
+      !!
+      !! ** Method  :   
+      !!      25hr mean outputs for shelf seas
+      !!
+      !! History :
+      !!   ?.0  !  07-04  (A. Hines) New routine, developed from dia_wri_foam
+      !!   3.4  !  02-13  (J. Siddorn) Routine taken from old dia_wri_foam
+      !!   3.6  !  08-14  (E. O'Dea) adapted for VN3.6
+      !!----------------------------------------------------------------------
+      !! * Modules used
+
+      IMPLICIT NONE
+
+      !! * Arguments
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+
+
+      !! * Local declarations
+      INTEGER ::   ji, jj, jk, jn
+
+      LOGICAL ::   ll_print = .FALSE.    ! =T print and flush numout
+      REAL(wp)                         ::   zsto, zout, zmax, zjulian, zdt, zmdi  ! temporary reals
+      INTEGER                          ::   i_steps                               ! no of timesteps per hour
+      REAL(wp), DIMENSION(jpi,jpj    ) ::   zw2d, un_dm, vn_dm                    ! temporary workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d                                  ! temporary workspace
+      REAL(wp), DIMENSION(jpi,jpj,3)   ::   zwtmb                                 ! temporary workspace
+      INTEGER                          ::   iyear0, nimonth0,iday0                ! start year,imonth,day
+
+      !!----------------------------------------------------------------------
+
+      ! 0. Initialisation
+      ! -----------------
+      ! Define frequency of summing to create 25 h mean
+      zdt = rdt
+      IF( nacc == 1 ) zdt = rdtmin
+
+      IF( MOD( 3600,INT(zdt) ) == 0 ) THEN
+         i_steps = 3600/INT(zdt)
+      ELSE
+         CALL ctl_stop('STOP', 'dia_wri_tide: timestep must give MOD(3600,rdt) = 0 otherwise no hourly values are possible')
+      ENDIF
+
+#if defined key_lim3 || defined key_lim2
+      CALL ctl_stop('STOP', 'dia_wri_tide not setup yet to do tidemean ice')
+#endif
+
+      ! local variable for debugging
+      ll_print = ll_print .AND. lwp
+
+      ! Sum of 25 hourly instantaneous values to give a 25h mean from 24hours
+      ! every day
+      IF( MOD( kt, i_steps ) == 0  .and. kt .ne. nn_it000 ) THEN
+
+         IF (lwp) THEN
+              WRITE(numout,*) 'dia_wri_tide : Summing instantaneous hourly diagnostics at timestep ',kt
+              WRITE(numout,*) '~~~~~~~~~~~~ '
+         ENDIF
+
+         tn_25h(:,:,:)        = tn_25h(:,:,:) + tsn(:,:,:,jp_tem)
+         sn_25h(:,:,:)        = sn_25h(:,:,:) + tsn(:,:,:,jp_sal)
+         CALL theta2t
+         rinsitu_t_25h(:,:,:)  = rinsitu_t_25h(:,:,:) + rinsitu_t(:,:,:)
+         CALL dia_calctmb( rinsitu_t(:,:,:),zwtmb )
+         insitu_bot_25h(:,:)  = insitu_bot_25h(:,:) + zwtmb(:,:,3)
+         zw3d(:,:,:)          = tsn(:,:,:,jp_tem)
+         CALL dia_calctmb( zw3d,zwtmb )
+         temp_bot_25h(:,:)    = temp_bot_25h(:,:) + zwtmb(:,:,3)
+         sshn_25h(:,:)        = sshn_25h(:,:) + sshn (:,:)
+         un_25h(:,:,:)        = un_25h(:,:,:) + un(:,:,:)
+         vn_25h(:,:,:)        = vn_25h(:,:,:) + vn(:,:,:)
+         wn_25h(:,:,:)        = wn_25h(:,:,:) + wn(:,:,:)
+         avt_25h(:,:,:)       = avt_25h(:,:,:) + avt(:,:,:)
+         avm_25h(:,:,:)       = avm_25h(:,:,:) + avm(:,:,:)
+# if defined key_zdfgls || defined key_zdftke
+         en_25h(:,:,:)        = en_25h(:,:,:) + en(:,:,:)
+#endif
+# if defined key_zdfgls
+         rmxln_25h(:,:,:)      = rmxln_25h(:,:,:) + mxln(:,:,:)
+#endif
+#if defined key_fabm
+      DO jn = 1, jp_fabm
+         fabm_25h(:,:,:,jn) = fabm_25h(:,:,:,jn) + trn(:,:,:,jp_fabm_m1+jn)
+      END DO
+      DO jn = 1, jp_fabm_3d
+         fabm_3d_25h(:,:,:,jn) = fabm_3d_25h(:,:,:,jn) + fabm_get_bulk_diagnostic_data(model, jn)
+      END DO
+      DO jn = 1, jp_fabm_surface
+         fabm_surface_25h(:,:,jn) = fabm_surface_25h(:,:,jn) + fabm_st2dn(:,:,jn)
+      END DO
+      DO jn = 1, jp_fabm_bottom
+         fabm_bottom_25h(:,:,jn) = fabm_bottom_25h(:,:,jn) + fabm_st2dn(:,:,jp_fabm_surface+jn)
+      END DO
+      DO jn = 1, jp_fabm_2d
+         fabm_2d_25h(:,:,jn) = fabm_2d_25h(:,:,jn) + fabm_get_horizontal_diagnostic_data(model,jn)
+      END DO
+      fabm_visib_25h(:,:,:) = fabm_visib_25h(:,:,:) + visib(:,:,:)
+#endif
+         cnt_25h = cnt_25h + 1
+
+         IF (lwp) THEN
+            WRITE(numout,*) 'dia_tide : Summed the following number of hourly values so far',cnt_25h
+         ENDIF
+
+      ENDIF ! MOD( kt, i_steps ) == 0
+
+         ! Write data for 25 hour mean output streams
+      IF( cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps*24) == 0 .AND. kt .NE. nn_it000 ) THEN
+
+            IF(lwp) THEN
+               WRITE(numout,*) 'dia_wri_tide : Writing 25 hour mean tide diagnostics at timestep', kt
+               WRITE(numout,*) '~~~~~~~~~~~~ '
+            ENDIF
+
+            tn_25h(:,:,:)        = tn_25h(:,:,:) / 25.0_wp
+            sn_25h(:,:,:)        = sn_25h(:,:,:) / 25.0_wp
+            rinsitu_t_25h(:,:,:)  = rinsitu_t_25h(:,:,:) / 25.0_wp
+            insitu_bot_25h(:,:)  = insitu_bot_25h(:,:) / 25.0_wp 
+            temp_bot_25h(:,:)    = temp_bot_25h(:,:) /25.0_wp
+            sshn_25h(:,:)        = sshn_25h(:,:) / 25.0_wp
+            un_25h(:,:,:)        = un_25h(:,:,:) / 25.0_wp
+            vn_25h(:,:,:)        = vn_25h(:,:,:) / 25.0_wp
+            wn_25h(:,:,:)        = wn_25h(:,:,:) / 25.0_wp
+            avt_25h(:,:,:)       = avt_25h(:,:,:) / 25.0_wp
+            avm_25h(:,:,:)       = avm_25h(:,:,:) / 25.0_wp
+# if defined key_zdfgls || defined key_zdftke
+            en_25h(:,:,:)        = en_25h(:,:,:) / 25.0_wp
+#endif
+# if defined key_zdfgls
+            rmxln_25h(:,:,:)       = rmxln_25h(:,:,:) / 25.0_wp
+#endif
+#if defined key_fabm
+            fabm_25h(:,:,:,:)       = fabm_25h(:,:,:,:)       / 25.0_wp
+            fabm_3d_25h(:,:,:,:)    = fabm_3d_25h(:,:,:,:)    / 25.0_wp
+            fabm_surface_25h(:,:,:) = fabm_surface_25h(:,:,:) / 25.0_wp
+            fabm_bottom_25h(:,:,:)  = fabm_bottom_25h(:,:,:)  / 25.0_wp
+            fabm_2d_25h(:,:,:)      = fabm_2d_25h(:,:,:)      / 25.0_wp
+            fabm_visib_25h(:,:,:)   = fabm_visib_25h(:,:,:)   / 25.0_wp
+#endif
+
+            IF (lwp)  WRITE(numout,*) 'dia_wri_tide : Mean calculated by dividing 25 hour sums and writing output'
+            zmdi=1.e+20 !missing data indicator for masking
+            ! write tracers (instantaneous)
+            zw3d(:,:,:) = tn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            CALL iom_put("temper25h", zw3d)   ! potential temperature
+            CALL theta2t                                                                    ! calculate insitu temp
+            zw3d(:,:,:) = rinsitu_t_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            CALL iom_put("tempis25h", zw3d)   ! in-situ temperature
+            zw2d(:,:) = insitu_bot_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
+            CALL iom_put("tempisbot25h", zw2d) ! bottom in-situ temperature
+            zw2d(:,:) = temp_bot_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
+            CALL iom_put("temperbot25h",zw2d) ! bottom potential temperature
+            zw3d(:,:,:) = sn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            CALL iom_put( "salin25h", zw3d  )   ! salinity
+            zw2d(:,:) = sshn_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
+            CALL iom_put( "ssh25h", zw2d )   ! sea surface 
+
+#if defined key_fabm
+            ! Write ERSEM variables
+            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  )
+            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  )
+            END DO
+            DO jn = 1, jp_fabm_surface
+               zw2d(:,:) = fabm_surface_25h(:,:,jn)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
+               CALL iom_put( TRIM(model%surface_state_variables(jn)%name)//"25h", zw2d  )
+            END DO
+            DO jn = 1, jp_fabm_bottom
+               zw2d(:,:) = fabm_bottom_25h(:,:,jn)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
+               CALL iom_put( TRIM(model%bottom_state_variables(jn)%name)//"25h", zw2d  )
+            END DO
+            DO jn = 1, jp_fabm_2d
+               zw2d(:,:) = fabm_2d_25h(:,:,jn)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
+               CALL iom_put( TRIM(model%horizontal_diagnostic_variables(jn)%name)//"25h", zw2d  )
+            END DO
+            zw3d(:,:,:) = fabm_visib_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            CALL iom_put( "visib25h", zw3d  )
+#endif
+
+            ! Write velocities (instantaneous)
+            zw3d(:,:,:) = un_25h(:,:,:)*umask(:,:,:) + zmdi*(1.0-umask(:,:,:))
+            CALL iom_put("vozocrtx25h", zw3d)    ! i-current
+            zw3d(:,:,:) = vn_25h(:,:,:)*vmask(:,:,:) + zmdi*(1.0-vmask(:,:,:))
+            CALL iom_put("vomecrty25h", zw3d  )   ! j-current
+
+            zw3d(:,:,:) = wn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            CALL iom_put("vomecrtz25h", zw3d )   ! k-current
+            zw3d(:,:,:) = avt_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            CALL iom_put("avt25h", zw3d )   ! diffusivity
+            zw3d(:,:,:) = avm_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            CALL iom_put("avm25h", zw3d)   ! viscosity
+#if defined key_zdftke || defined key_zdfgls 
+!           zw3d(:,:,:) = en_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+!           CALL iom_put("tke25h", zw3d)   ! tke
+#endif
+#if defined key_zdfgls 
+!           zw3d(:,:,:) = rmxln_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+!           CALL iom_put( "mxln25h",zw3d)
+#endif
+
+            ! After the write reset the values to cnt=1 and sum values equal current value 
+            tn_25h(:,:,:) = tsn(:,:,:,jp_tem)
+            sn_25h(:,:,:) = tsn(:,:,:,jp_sal)
+            CALL theta2t
+            rinsitu_t_25h(:,:,:) = rinsitu_t(:,:,:)
+            CALL dia_calctmb( rinsitu_t(:,:,:),zwtmb )
+            insitu_bot_25h(:,:) = zwtmb(:,:,3)
+            CALL dia_calctmb( tn_25h(:,:,:),zwtmb)
+            temp_bot_25h(:,:) = zwtmb(:,:,3)
+            sshn_25h(:,:) = sshn (:,:)
+            un_25h(:,:,:) = un(:,:,:)
+            vn_25h(:,:,:) = vn(:,:,:)
+            wn_25h(:,:,:) = wn(:,:,:)
+            avt_25h(:,:,:) = avt(:,:,:)
+            avm_25h(:,:,:) = avm(:,:,:)
+# if defined key_zdfgls || defined key_zdftke
+            en_25h(:,:,:) = en(:,:,:)
+#endif
+# if defined key_zdfgls
+            rmxln_25h(:,:,:) = mxln(:,:,:)
+#endif
+#if defined key_fabm
+      DO jn = 1, jp_fabm
+         fabm_25h(:,:,:,jn) = trn(:,:,:,jp_fabm_m1+jn)
+      END DO
+      DO jn = 1, jp_fabm_3d
+         fabm_3d_25h(:,:,:,jn) = fabm_get_bulk_diagnostic_data(model, jn)
+      END DO
+      DO jn = 1, jp_fabm_surface
+         fabm_surface_25h(:,:,jn) = fabm_st2dn(:,:,jn)
+      END DO
+      DO jn = 1, jp_fabm_bottom
+         fabm_bottom_25h(:,:,jn) = fabm_st2dn(:,:,jp_fabm_surface+jn)
+      END DO
+      DO jn = 1, jp_fabm_2d
+         fabm_2d_25h(:,:,jn) = fabm_get_horizontal_diagnostic_data(model,jn)
+      END DO
+      fabm_visib_25h(:,:,:) = visib(:,:,:)
+#endif
+            cnt_25h = 1
+            IF (lwp)  WRITE(numout,*) 'dia_wri_tide : After 25hr mean write, reset sum to current value and cnt_25h to one for overlapping average',cnt_25h
+
+      ENDIF !  cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps * 24) == 0 .AND. kt .NE. nn_it000
+
+
+   END SUBROUTINE dia_25h 
+
+   !!======================================================================
+END MODULE dia25h
Index: /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaharmana.F90
===================================================================
--- /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaharmana.F90	(revision 11137)
+++ /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaharmana.F90	(revision 11137)
@@ -0,0 +1,667 @@
+MODULE harmonic_analysis
+   !!======================================================================
+   !!                       ***  MODULE  example  ***
+   !! Ocean physics:  On line harmonic analyser
+   !!                 
+   !!=====================================================================
+
+#if defined key_harm_ana
+   !!----------------------------------------------------------------------
+   !!   'key_harm_ana'  :                Calculate harmonic analysis
+   !!----------------------------------------------------------------------
+   !!   harm_ana        :
+   !!   harm_ana_init   :
+   !!----------------------------------------------------------------------
+
+   USE oce             ! ocean dynamics and tracers 
+   USE dom_oce         ! ocean space and time domain
+   USE iom
+   USE in_out_manager  ! I/O units
+   USE phycst          ! physical constants
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE bdy_par         ! Unstructured boundary parameters
+   USE bdy_oce         ! ocean open boundary conditions
+   USE bdytides        ! tidal bdy forcing
+   USE daymod          ! calendar
+   USE tideini
+   USE restart
+   USE ioipsl, ONLY : ju2ymds    ! for calendar
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! *  Routine accessibility
+   PUBLIC harm_ana    ! routine called in step.F90 module
+
+   !! * Module variables
+   INTEGER, PARAMETER ::  nharm_max  =   jptides_max ! max number of harmonics to be analysed 
+   INTEGER, PARAMETER ::  nhm_max    =   2*jptides_max+1 
+   INTEGER, PARAMETER ::  nvab       = 2 ! number of 3D variables
+   INTEGER            ::  nharm
+   INTEGER            ::  nhm 
+   INTEGER ::                 & !!! ** toto namelist (namtoto) **
+      nflag  =  1                ! default value of nflag 
+   REAL(wp), DIMENSION(nharm_max) ::                & 
+      om_tide                     ! tidal frequencies ( rads/sec)
+   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:)   ::                & 
+      bzz,c,x    ! work arrays
+   REAL(wp) :: cca,ssa,zm,bt
+   REAL(wp) :: ccau,ccav,ssau,ssav
+   REAL(wp), PUBLIC, ALLOCATABLE,DIMENSION(:) :: anau, anav, anaf   ! nodel/phase corrections used by diaharmana
+!
+   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:) ::   &
+      bssh,bubar,bvbar        ! work array for ssh anaylsis
+   REAL(wp), ALLOCATABLE,SAVE, DIMENSION(:,:,:) ::   &
+      cosampz, cosampu,cosampv,      &
+      sinampz, sinampu,sinampv
+   REAL(WP), ALLOCATABLE,SAVE,DIMENSION(:,:) :: cc,a
+!
+   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:) ::   &
+      gout,hout        ! arrays for output
+   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:) ::   &
+      huout,hvout,guout,gvout        ! arrays for output
+   REAL(wp), PUBLIC ::   fjulday_startharm       !: Julian Day since start of harmonic analysis
+
+   !! * Substitutions
+
+   !!----------------------------------------------------------------------
+   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! or LIM 2.0 , UCL-LOCEAN-IPSL (2005)
+   !! or  TOP 1.0 , LOCEAN-IPSL (2005)
+   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/module_example,v 1.3 2005/03/27 18:34:47 opalod Exp $ 
+   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE harm_ana( kt )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE harm_ana  ***
+      !!
+      !! ** Purpose :   Harmonic analyser
+      !!
+      !! ** Method  :   
+      !!
+      !! ** Action  : - first action (share memory array/varible modified
+      !!                in this routine
+      !!              - second action .....
+      !!              - .....
+      !!
+      !! References :
+      !!   Give references if exist otherwise suppress these lines
+      !!
+      !! History :
+      !!   9.0  !  03-08  (Autor Names)  Original code
+      !!        !  02-08  (Author names)  brief description of modifications
+      !!----------------------------------------------------------------------
+      !! * Modules used
+      
+      !! * arguments
+      INTEGER, INTENT( in  ) ::   &  
+         kt                          ! describe it!!!
+
+      !! * local declarations
+      INTEGER ::   ji, jj, jk,  &        ! dummy loop arguments
+                   ih,i1,i2,iv,jgrid
+
+
+      !!--------------------------------------------------------------------
+
+
+
+      IF( kt == nit000  )   CALL harm_ana_init    ! Initialization (first time-step only)
+
+! this bit done every time step
+        !nharm=ncpt_anal
+        nharm=nb_harmo
+        nhm=2*nb_harmo+1
+       
+        c(1)=1.0
+
+        do ih=1,nb_harmo
+             c(2*ih)=    cos((fjulday-fjulday_startharm)*86400._wp*om_tide(ih)  )
+             c(2*ih+1)=  sin((fjulday-fjulday_startharm)*86400._wp*om_tide(ih)  )
+        enddo
+
+
+         do ji=1,jpi
+          do jj=1,jpj
+            do ih=1,nhm
+            bssh (ih,ji,jj)=bssh (ih,ji,jj)+c(ih)*sshn(ji,jj)
+            bubar(ih,ji,jj)=bubar(ih,ji,jj)+c(ih)*un_b(ji,jj)
+            bvbar(ih,ji,jj)=bvbar(ih,ji,jj)+c(ih)*vn_b(ji,jj)
+            enddo
+          enddo
+        enddo
+!  
+!        IF(lwp) WRITE(666,'(4E15.5)') adatrj*86400._wp, sshn(10,10),bssh(2,10,10),bssh(3,10,10)
+!
+         do i1=1,nhm
+           do i2=1,nhm
+             cc(i1,i2)=cc(i1,i2)+c(i1)*c(i2)
+            enddo
+         enddo
+!
+! At End of run
+
+      IF (kt ==  nitend  ) then
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'harm_ana : harmonic analysis of tides at end of run'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
+      CALL harm_rst_write(kt)     ! Dump out data for a restarted run
+
+
+      if(ln_harm_ana_compute) then
+               WRITE(numout,*) "Computing harmonics at last step"
+!
+        cosampu=0.0_wp
+        sinampu=0.0_wp
+        cosampv=0.0_wp
+        sinampv=0.0_wp
+        cosampz=0.0_wp
+        sinampz=0.0_wp
+!
+     do jgrid=1,3 ! elevation, Ubar,Vbar
+        do ji=1,jpi
+            do jj=1,jpj
+               bt=1.0
+               do ih=1,nhm
+                   if(jgrid .eq. 1) then
+                     bzz(ih)=bssh(ih,ji,jj)
+                   endif
+                   if(jgrid .eq. 2) then
+                     bzz(ih)=bubar(ih,ji,jj)
+                   endif
+                   if(jgrid .eq. 3) then
+                     bzz(ih)=bvbar(ih,ji,jj)
+                   endif
+                   bt=bt*bzz(ih)
+               enddo
+
+               do i1=1,nhm
+                   do i2=1,nhm
+                       a(i1,i2)=cc(i1,i2)
+                   enddo
+               enddo
+!	now do gaussian elimination of the system
+!	a * x = b
+! 	the matrix x is (a0,a1,b1,a2,b2 ...)
+!	the matrix a and rhs b solved here for x
+              x=0.0d0
+              if(bt.ne.0.) then
+                  if(lwp .and. ji == 10 .and. ji == 10) then
+              endif
+              call gelim(a,bzz,x,nhm)
+
+              do ih=1,nb_harmo
+                if(jgrid .eq. 1) then
+                   cosampz(ih,ji,jj)=x(ih*2)
+                   sinampz(ih,ji,jj)=x(ih*2+1)
+                endif
+                if(jgrid .eq. 2) then
+                   cosampu(ih,ji,jj)=x(ih*2)
+                   sinampu(ih,ji,jj)=x(ih*2+1)
+                endif
+                if(jgrid .eq. 3) then
+                   cosampv(ih,ji,jj)=x(ih*2)
+                   sinampv(ih,ji,jj)=x(ih*2+1)
+                endif
+              enddo
+                if(jgrid .eq. 1) then
+                  cosampz(0,ji,jj)=x(1)
+                  sinampz(0,ji,jj)=0.0_wp
+                endif
+                if(jgrid .eq. 2) then
+                  cosampu(0,ji,jj)=x(1)
+                  sinampu(0,ji,jj)=0.0_wp
+                endif
+                if(jgrid .eq. 3) then
+                  cosampv(0,ji,jj)=x(1)
+                  sinampv(0,ji,jj)=0.0_wp
+                endif
+              endif
+           enddo
+        enddo
+      enddo ! jgrid
+!
+!
+        CALL harm_ana_out     ! output analysis (last time step)
+     ELSE ! ln_harmana_compute 
+         WRITE(numout,*) "Skipping Computing harmonics at last step"
+     ENDIF
+       
+    ENDIF
+   END SUBROUTINE harm_ana
+
+   SUBROUTINE harm_ana_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE harm_ana_init  ***
+      !!                   
+      !! ** Purpose :   initialization of ....
+      !!
+      !! ** Method  :   blah blah blah ...
+      !!
+      !! ** input   :   Namlist namexa
+      !!
+      !! ** Action  :   ...  
+      !!
+      !! history :
+      !!   9.0  !  03-08  (Autor Names)  Original code
+      !!----------------------------------------------------------------------
+      !! * local declarations
+      INTEGER ::   ji, jj, jk, jit,ih   ! dummy loop indices
+
+      !!----------------------------------------------------------------------
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'harm_init : initialization of harmonic analysis of tides'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
+     ! DO ALLOCATIONS
+
+      ALLOCATE( bubar(nb_harmo*2+1,jpi,jpj) )
+      ALLOCATE( bvbar(nb_harmo*2+1,jpi,jpj) )
+      ALLOCATE( bssh (nb_harmo*2+1,jpi,jpj) )
+
+      ALLOCATE( cosampz(0:nb_harmo*2+1,jpi,jpj))
+      ALLOCATE( cosampu(0:nb_harmo*2+1,jpi,jpj))
+      ALLOCATE( cosampv(0:nb_harmo*2+1,jpi,jpj))
+
+      ALLOCATE( sinampz(0:nb_harmo*2+1,jpi,jpj))
+      ALLOCATE( sinampu(0:nb_harmo*2+1,jpi,jpj))
+      ALLOCATE( sinampv(0:nb_harmo*2+1,jpi,jpj))
+
+      ALLOCATE( cc(nb_harmo*2+1,nb_harmo*2+1) )
+
+      ALLOCATE( a(nb_harmo*2+1,nb_harmo*2+1) )
+
+      ALLOCATE( anau(nb_harmo) )
+      ALLOCATE( anav(nb_harmo) )
+      ALLOCATE( anaf(nb_harmo) )
+
+      ALLOCATE( bzz(nb_harmo*2+1) )
+      ALLOCATE( x(nb_harmo*2+1) )
+      ALLOCATE( c(nb_harmo*2+1) )
+
+      ALLOCATE( gout(jpi,jpj))
+      ALLOCATE( hout(jpi,jpj))
+
+      ALLOCATE( guout(jpi,jpj))
+      ALLOCATE( huout(jpi,jpj))
+
+      ALLOCATE( gvout(jpi,jpj))
+      ALLOCATE( hvout(jpi,jpj))
+
+        ! END ALLOCATE 
+
+                 do ih=1,nb_harmo
+                  om_tide(ih)= omega_tide(ih) 
+                 enddo
+       if(ln_harmana_read) then
+               WRITE(numout,*) "Reading previous harmonic data from previous run"
+               ! Need to read in  bssh bz, cc anau anav and anaf 
+               call harm_rst_read  ! This reads in from the previous day
+                                   ! Currrently the data in in assci format
+       else
+               WRITE(numout,*) "Starting harmonic analysis from Fresh "
+         bssh(:,:,:)  = 0.0_wp
+         bubar(:,:,:) = 0.0_wp
+         bvbar(:,:,:) = 0.0_wp
+
+         cc=0.0_wp
+
+         anau(:) = 0.0_wp
+         anav(:) = 0.0_wp
+         anaf(:) = 0.0_wp
+         bzz(:) = 0.0_wp
+         x(:) = 0.0_wp
+         c(:) = 0.0_wp
+
+         DO ih = 1, nb_harmo
+            anau(ih) = utide(ih)
+            anav(ih) = v0tide(ih)
+            anaf(ih) = ftide(ih)
+
+         END DO
+         fjulday_startharm=fjulday !Set this at very start and store
+      endif
+
+   END SUBROUTINE harm_ana_init
+!
+   SUBROUTINE gelim (a,b,x,n)
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE harm_ana  ***
+      !!
+      !! ** Purpose :   Guassian elimination
+      !!
+      !! ** Method  :   
+      !!
+      !! ** Action  : - first action (share memory array/varible modified
+      !!                in this routine
+      !!              - second action .....
+      !!              - .....
+      !!
+      !! References :
+      !!   Give references if exist otherwise suppress these lines
+      !!
+      !! History :
+        implicit none
+!
+        integer  :: n
+   REAL(WP) :: b(nb_harmo*2+1),a(nb_harmo*2+1,nb_harmo*2+1)
+   REAL(WP) :: x(nb_harmo*2+1)
+        INTEGER  :: row,col,prow,pivrow,rrow,ntemp
+   REAL(WP) ::  atemp
+   REAL(WP) ::  pivot
+   REAL(WP) ::  m
+   REAL(WP) ::  tempsol(nb_harmo*2+1)
+
+
+        do row=1,n-1
+           pivrow=row
+           pivot=a(row,n-row+1)
+           do prow=row+1,n
+              if (abs(a(prow,n-row+1)).gt.abs(pivot)  ) then
+                 pivot=a(prow,n-row+1)
+                 pivrow=prow
+              endif
+           enddo
+!	swap row and prow
+           if ( pivrow .ne. row ) then
+              atemp=b(pivrow)
+              b(pivrow)=b(row)
+              b(row)=atemp
+              do col=1,n
+                 atemp=a(pivrow,col)
+                 a(pivrow,col)=a(row,col)
+                 a(row,col)=atemp
+              enddo
+           endif
+
+           do rrow=row+1,n
+              if (a(row,row).ne.0) then
+   
+                 m=-a(rrow,n-row+1)/a(row,n-row+1)
+                 do col=1,n
+                    a(rrow,col)=m*a(row,col)+a(rrow,col)
+                 enddo
+                 b(rrow)=m*b(row)+b(rrow)
+              endif
+           enddo
+        enddo
+!	back substitution now
+
+        x(1)=b(n)/a(n,1)
+        do row=n-1,1,-1
+           x(n-row+1)=b(row)
+           do col=1,(n-row)
+              x(n-row+1)=(x(n-row+1)-a(row,col)*x(col)) 
+           enddo
+
+           x(n-row+1)=(x(n-row+1)/a(row,(n-row)+1))
+        enddo
+
+        return
+        END SUBROUTINE gelim
+
+   SUBROUTINE harm_ana_out
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE harm_ana_init  ***
+      !!                   
+      !! ** Purpose :   initialization of ....
+      !!
+      !! ** Method  :   blah blah blah ...
+      !!
+      !! ** input   :   Namlist namexa
+      !!
+      !! ** Action  :   ...  
+      !!
+      !! history :
+      !!   9.0  !  03-08  (Autor Names)  Original code
+      !!----------------------------------------------------------------------
+        USE dianam          ! build name of file (routine)
+ 
+      !! * local declarations
+      INTEGER ::   ji, jj, jk, jgrid,jit, ih, jjk   ! dummy loop indices
+      INTEGER :: nh_T
+      INTEGER :: nid_harm
+      CHARACTER (len=40) ::           &
+         clhstnamt, clop1, clop2          ! temporary names 
+      CHARACTER (len=40) ::           &
+         clhstnamu, clhstnamv   ! temporary names 
+      REAL(wp) ::   &
+         zsto1, zsto2, zout, zmax,            &  ! temporary scalars
+         zjulian, zdt, zmdi  
+         
+        do jgrid=1,3
+          do ih=1,nb_harmo
+            hout = 0.0
+            gout = 0.0
+            do jj=1,nlcj
+                do ji=1,nlci
+                    if(jgrid .eq. 1) then ! SSH 
+                       cca=cosampz(ih,ji,jj)
+                       ssa=sinampz(ih,ji,jj)
+                    endif
+                    if(jgrid .eq. 2) then ! UBAR
+                       cca=cosampu(ih,ji,jj)
+                       ssa=sinampu(ih,ji,jj)
+                    endif
+                    if(jgrid .eq. 3) then ! VBAR
+                       cca=cosampv(ih,ji,jj)
+                       ssa=sinampv(ih,ji,jj)
+                    endif
+
+                    hout(ji,jj)=sqrt(cca**2+ssa**2)
+        
+                    if (cca.ne.0.0) then
+                         gout(ji,jj)=(180.0/rpi)*atan(ssa/cca)
+                    else
+                        if (ssa.ne.0.0) then
+                            if (ssa.gt. 0.0) then
+                                 gout(ji,jj)=90.0
+                            else
+                                 gout(ji,jj)=270.0
+                            endif
+                        else
+                            gout(ji,jj)=0.0
+                        endif
+                    endif
+
+                    if (gout(ji,jj).lt.0) then
+                        gout(ji,jj)=gout(ji,jj)+180.0
+                    endif
+                    if (ssa.lt.0) then
+                        gout(ji,jj)=gout(ji,jj)+180.0
+                    endif
+                    
+                    if (hout(ji,jj).ne.0) then
+                        hout(ji,jj)=hout(ji,jj)/anaf(ih)
+                    endif
+                    if (gout(ji,jj).ne.0) then
+                                                !Convert Rad to degree and take
+                                                !modulus
+                         gout(ji,jj)=gout(ji,jj)+MOD( (anau(ih)+anav(ih))/rad , 360.0)
+                        if (gout(ji,jj).gt.360.0) then
+                            gout(ji,jj)=gout(ji,jj)-360.0
+                        else if (gout(ji,jj).lt.0) then
+                            gout(ji,jj)=gout(ji,jj)+360.0
+                        endif
+                    endif
+                enddo
+            enddo
+!NETCDF OUTPUT
+            if(jgrid==1) then
+             WRITE(numout,*) TRIM(Wave(ntide(ih))%cname_tide)//'amp'
+             WRITE(numout,*) TRIM(Wave(ntide(ih))%cname_tide)//'phase'
+             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'amp', hout(:,:) )
+             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'phase', gout(:,:) )
+            endif
+            if(jgrid==2) then
+             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'amp_u', hout(:,:) )
+             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'phase_u', gout(:,:) )
+            endif
+            if(jgrid==3) then
+             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'amp_v', hout(:,:) )
+             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'phase_v', gout(:,:) )
+            endif
+          enddo
+        enddo
+ 
+!
+   END SUBROUTINE harm_ana_out
+
+   SUBROUTINE harm_rst_write(kt)
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE harm_ana_init  ***
+      !!                   
+      !! ** Purpose :  To write out cummulated Tidal Harmomnic data to file for
+      !!               restarting
+      !!
+      !! ** Method  :   restart files will be dated by default
+      !!
+      !! ** input   :   
+      !!
+      !! ** Action  :   ...  
+      !!
+      !! history :
+      !!   0.0  !  01-16  (Enda O'Dea)  Original code
+      !! ASSUMES  dated file for rose  , can change later to be more generic
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt     ! ocean time-step
+      INTEGER             ::   iyear, imonth, iday,ih
+      REAL (wp)           ::   zsec
+      !!
+      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step define as a character
+      CHARACTER(LEN=50)   ::   clname   ! ocean output restart file name
+      CHARACTER(LEN=150)  ::   clpath   ! full path to ocean output restart file
+      CHARACTER(LEN=250)  ::   clfinal   ! full name
+
+      CALL ju2ymds( fjulday , iyear, imonth, iday, zsec)
+
+      WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday
+      ! create the file
+      clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_harm_ana"
+      clpath = TRIM(cn_ocerst_outdir)
+      IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
+      WRITE(numout,*) 'Open tidal harmonics restart file for writing: ',TRIM(clpath)//clname
+
+
+
+
+!!      clhstnam=TRIM(cexper)//'.restart_harm_ana'
+      ! Open a file to write the data into
+         write (clfinal,'(a,''_'',i4.4)') trim(clpath)//trim(clname),nproc
+         open(66,file=TRIM(clfinal))
+         write(66,'(1e20.9)') cc
+         !These Three which contain the most data can be moved to a regular
+         !restart file
+
+            CALL iom_rstput( kt, nitrst, numrow, 'Mean_bssh'     , bssh(1,:,:)       )
+            CALL iom_rstput( kt, nitrst, numrow, 'Mean_bubar'    , bubar(1,:,:)       )
+            CALL iom_rstput( kt, nitrst, numrow, 'Mean_bvbar'    , bvbar(1,:,:)       )
+         do ih=1,nb_harmo
+            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bssh_cos'    , bssh (ih*2  , : , : )     )
+            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bssh_sin'    , bssh (ih*2+1, : , : )     )
+
+            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bubar_cos'   , bubar(ih*2  , : , : )    )
+            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bubar_sin'   , bubar(ih*2+1, : , : )    )
+
+            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bvbar_cos'   , bvbar(ih*2  , : , : )    )
+            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bvbar_sin'   , bvbar(ih*2+1, : , : )    )
+         enddo
+
+         write(66,'(1e20.9)') anau
+         write(66,'(1e20.9)') anav
+         write(66,'(1e20.9)') anaf
+         write(66,'(1e25.20)') fjulday_startharm
+         close(66)
+ 
+   END SUBROUTINE harm_rst_write
+
+   SUBROUTINE harm_rst_read
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE harm_ana_init  ***
+      !!                   
+      !! ** Purpose :  To read in  cummulated Tidal Harmomnic data to file for
+      !!               restarting
+      !!
+      !! ** Method  :   
+      !!
+      !! ** input   :   
+      !!
+      !! ** Action  :   ...  
+      !!
+      !! history :
+      !!   0.0  !  01-16  (Enda O'Dea)  Original code
+      !! ASSUMES  dated file for rose  , can change later to be more generic
+      !!----------------------------------------------------------------------
+      CHARACTER(LEN=50)   ::   clname   ! ocean output restart file name
+      CHARACTER(LEN=150)  ::   clpath   ! full path to ocean output restart file
+      CHARACTER(LEN=250)  ::   clfinal   ! full name
+      INTEGER             ::   ih
+
+      ! create the file
+      clname = "restart_harm_ana"
+      clpath = TRIM(cn_ocerst_indir)
+      IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
+      WRITE(numout,*) 'Open tidal harmonics restart file for reading: ',TRIM(clpath)//clname
+
+
+
+
+!!      clhstnam=TRIM(cexper)//'.restart_harm_ana'
+      ! Open a file to read the data ifrom
+         write (clfinal,'(a,''_'',i4.4)') trim(clpath)//trim(clname),nproc
+         open(66,file=TRIM(clfinal),status='old')
+         read(66,'(1e20.9)') cc
+!2D regular fields can be read from normal restart this saves space and handy to
+!view in netcdf format also.
+         CALL iom_get( numror,jpdom_autoglo, 'Mean_bssh'     , bssh(1,:,:)       )
+         CALL iom_get( numror,jpdom_autoglo, 'Mean_bubar'    , bubar(1,:,:)       )
+         CALL iom_get( numror,jpdom_autoglo, 'Mean_bvbar'    , bvbar(1,:,:)       )
+          
+         do ih=1,nb_harmo
+
+            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bssh_cos'    , bssh (ih*2  , : , : )     )
+            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bssh_sin'    , bssh (ih*2+1, : , : )     )
+
+            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bubar_cos'   , bubar(ih*2  , : , : )    )
+            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bubar_sin'   , bubar(ih*2+1, : , : )    )
+
+            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bvbar_cos'   , bvbar(ih*2  , : , : )    )
+            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bvbar_sin'   , bvbar(ih*2+1, : , : )    )
+         enddo
+         read(66,'(1e20.9)') anau
+         read(66,'(1e20.9)') anav
+         read(66,'(1e20.9)') anaf
+         read(66,'(1e25.20)') fjulday_startharm
+      WRITE(numout,*) 'Checking anaf is correct'
+      WRITE(numout,*) anaf
+      WRITE(numout,*) '---end of anaf checking----'
+
+         close(66)
+
+   END SUBROUTINE harm_rst_read
+
+   !!======================================================================
+#else
+!!---------------------------------------------------------------------------------
+!!   Dummy module                                   NO harmonic Analysis
+!!---------------------------------------------------------------------------------
+        CONTAINS
+           SUBROUTINE harm_rst_write(kt)     ! Dummy routine
+           END SUBROUTINE harm_rst_write
+           SUBROUTINE harm_rst_read    ! Dummy routine
+           END SUBROUTINE harm_rst_read
+           SUBROUTINE harm_ana_out      ! Dummy routine
+           END SUBROUTINE harm_ana_out
+           SUBROUTINE harm_ana_init
+           END SUBROUTINE harm_ana_init
+           SUBROUTINE harm_ana( kt )
+           END SUBROUTINE harm_ana
+           SUBROUTINE gelim (a,b,x,n)
+           END SUBROUTINE gelim
+           
+#endif
+
+END MODULE harmonic_analysis
Index: /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diainsitutem.F90
===================================================================
--- /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diainsitutem.F90	(revision 11137)
+++ /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diainsitutem.F90	(revision 11137)
@@ -0,0 +1,124 @@
+MODULE diainsitutem
+   !!======================================================================
+   !!                       ***  MODULE  diaharm  ***
+   !! Harmonic analysis of tidal constituents
+   !!======================================================================
+   !! Initial version: D. Acreman June 2006 Original Code
+   !! History :  3.6  !  2014  (E O'Dea) for 3.6 25hour mean output
+   !!----------------------------------------------------------------------
+
+  USE dom_oce  ! ocean space and time domain
+  USE oce,            ONLY: tsn
+  USE par_oce,        ONLY: jpi, jpj, jpk, jpkm1
+  USE lbclnk,         ONLY: lbc_lnk
+  USE lib_mpp         ! MPP library
+
+  IMPLICIT NONE
+  PRIVATE
+
+  PUBLIC theta2t
+
+  !! * Accessibility
+  PUBLIC insitu_tem_alloc          ! routines called by step.F90
+
+  REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: rinsitu_t
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+
+CONTAINS
+
+   REAL FUNCTION insitu_tem_alloc()
+      !!----------------------------------------------------------------------
+      INTEGER, DIMENSION(2) :: ierr
+      !!----------------------------------------------------------------------
+      !
+      ierr = 0
+      !
+      ALLOCATE( rinsitu_t(jpi,jpj,jpk), STAT=ierr(1) )
+         !
+      insitu_tem_alloc = MAXVAL(ierr)
+      IF( lk_mpp )   CALL mpp_sum( insitu_tem_alloc )
+      !
+  END FUNCTION insitu_tem_alloc
+
+!-----------------------------------------------------------------------------------
+!
+! Calculate the insitu temperature by integrating the adiabatic lapse rate from zero 
+! to pressure at depth of tracer level. Based on UM subroutine POTTEM and function ATG
+!
+! Initial version: D. Acreman June 2006
+!
+!-----------------------------------------------------------------------------------
+
+  SUBROUTINE theta2t()
+
+    INTEGER, PARAMETER :: inum_steps=10                ! number of steps in integration
+    INTEGER            :: jstep                        ! iteration counter
+    INTEGER            :: ji, jj, jk                  ! loop indices
+    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zP            ! pressure (decibars)
+    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zT            ! temperature at pressure p
+    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zTB           ! temperature at p-dp
+    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zTA           ! temperature at p+dp
+    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zDP           ! pressure step
+    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zSS           ! salinity(PSU) - 35.0
+    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zLAPSE        ! adiabatic lapse rate
+
+!CDIR IEXPAND (ATG)
+
+! Set the pressure interval for the integration. The integration is carried out from 
+! zero (pressure at the surface) to the pressure at the depth of the tracer point. The
+! pressure in decibars is represented by the depth in metres. Pressures are "Oceanographic"
+! pressures equal to absolute pressure minus one atmosphere
+     zDP(:,:,:) = 0.0
+     DO jk = 1, jpkm1
+        DO jj = 1, jpj
+           DO ji = 1, jpi
+              ! These loops expanded for case where fsdept may be 1D
+              zDP(ji,jj,jk) = fsdept(ji,jj,jk) / real(inum_steps)
+           END DO
+        END DO
+     END DO
+
+! Salinity at each point
+     zSS(:,:,:)      = tsn(:,:,:,jp_sal) - 35.0 
+
+! Set initial values of temperature and pressure. zT is the temperature at pressure zP, 
+! zTB is the temperature at pressure zP-zdP and zTA is the temperature at pressure zP+zdP
+     zT(:,:,:)  = tsn(:,:,:,jp_tem)
+     zP(:,:,:)  = 0.0         ! Pressure at surface
+     CALL ATG(zP, zT, zSS, zLAPSE)
+     zTB(:,:,:) = zT(:,:,:)  - zLAPSE(:,:,:) * zDP(:,:,:)
+
+     interation: DO jstep=1, inum_steps
+        ! Calculate lapse rate (dT/dP) and hence TA
+        CALL ATG(zP, zT, zSS, zLAPSE)
+        zTA(:,:,:) = zTB(:,:,:) + 2.0 * zLAPSE(:,:,:) * zDP(:,:,:)
+        ! Have calculated TB, T and TA for this pressure, now advance solution by dP
+        zP(:,:,:)  = zP(:,:,:)  + zDP
+        zTB(:,:,:) = zT(:,:,:)
+        zT(:,:,:)  = zTA(:,:,:)
+
+     END DO interation
+
+     rinsitu_t(:,:,:) = zT(:,:,:) * tmask(:,:,:)
+     CALL lbc_lnk( rinsitu_t,  'T', 1.0)
+
+   END SUBROUTINE theta2t
+
+   SUBROUTINE ATG(P,T,DS,zLAPSE)
+
+     REAL, INTENT(IN   ) :: P(jpi,jpj,jpk)        ! PRESSURE (DECIBARS)
+     REAL, INTENT(IN   ) :: T(jpi,jpj,jpk)        ! TEMPERATURE (DEG C)
+     REAL, INTENT(IN   ) :: DS(jpi,jpj,jpk)       ! SALINITY (PSU) -35.0
+     REAL, INTENT(  OUT) :: zLAPSE(jpi,jpj,jpk)   ! LAPSE RATE
+
+       zLAPSE = ((( -2.1687E-16*T+1.8676E-14)*T-4.6206E-13)*P &
+           + ((2.7759E-12*T-1.1351E-10)*DS+((-5.4481E-14*T    &
+           + 8.733E-12)*T-6.7795E-10)*T+1.8741E-8))*P         &
+           + (-4.2393E-8*T+1.8932E-6)*DS                      &
+           + ((6.6228E-10*T-6.836E-8)*T+8.5258E-6)*T+3.5803E-5 
+
+   END SUBROUTINE ATG
+
+END MODULE diainsitutem
Index: /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaopfoam.F90
===================================================================
--- /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaopfoam.F90	(revision 11137)
+++ /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaopfoam.F90	(revision 11137)
@@ -0,0 +1,182 @@
+MODULE diaopfoam
+   !!======================================================================
+   !!                       ***  MODULE  diaopfoam  ***
+   !! Output stream for operational use
+   !!======================================================================
+   !! History :  3.6  !  2016  (P Sykes)  Original code
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain
+   USE diainsitutem, ONLY: rinsitu_t, theta2t
+   USE in_out_manager  ! I/O units
+   USE iom             ! I/0 library
+   USE wrk_nemo        ! working arrays
+   USE diatmb
+   USE diurnal_bulk
+   USE cool_skin
+
+
+   IMPLICIT NONE
+   PRIVATE
+
+   LOGICAL , PUBLIC ::   ln_diaopfoam     !: Diaopfoam output 
+   PUBLIC   dia_diaopfoam_init            ! routine called by nemogcm.F90
+   PUBLIC   dia_diaopfoam                 ! routine called by diawri.F90
+   PUBLIC   calc_max_cur                  ! routine called by diaopfoam.F90
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.6 , NEMO Consortium (2014)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dia_diaopfoam_init
+      !!---------------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_wri_diaop_init  ***
+      !!     
+      !! ** Purpose: Initialization of diaopfoam namelist 
+      !!        
+      !! ** Method : Read namelist
+      !!   History
+      !!   3.4  !  03-14  (P. Sykes) Routine to initialize dia_wri_diaop
+      !!---------------------------------------------------------------------------
+      !!
+      INTEGER            ::   ios      ! Local integer output status for namelist read
+      INTEGER            ::   ierror   ! local integer
+      !!
+      NAMELIST/nam_diadiaop/ ln_diaopfoam
+      !!
+      !!----------------------------------------------------------------------
+      !
+      ln_diaopfoam = .false.         ! default value for diaopfoam stream
+      REWIND ( numnam_ref )              ! Read Namelist nam_diadiaop in reference namelist : 3D hourly diagnostics
+      READ   ( numnam_ref, nam_diadiaop, IOSTAT=ios, ERR= 901 )
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diadiaop in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist nam_diadiaop in configuration namelist  3D hourly diagnostics
+      READ  ( numnam_cfg, nam_diadiaop, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diadiaop in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, nam_diadiaop )
+      !
+      IF(lwp) THEN                   ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'dia_diaopfoam_init : Output Diaopfoam Diagnostics'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist nam_diadiaop : set diaopfoam outputs '
+         WRITE(numout,*) '      Switch for diaopfoam diagnostics (T) or not (F)  ln_diaopfoam  = ', ln_diaopfoam
+      ENDIF
+    END SUBROUTINE dia_diaopfoam_init
+
+    SUBROUTINE dia_diaopfoam
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE dia_diaopfoam  ***
+      !! ** Purpose :   Write 3D hourly diagnostics for operational use
+      !!
+      !!
+      !! History :
+      !!   3.6  !  11-16  (P. Sykes) 
+      !!         
+      !!--------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace
+      REAL(wp) :: zmdi
+      REAL(wp), POINTER, DIMENSION(:,:)   :: zwu
+      REAL(wp), POINTER, DIMENSION(:,:)   :: zwv
+      REAL(wp), POINTER, DIMENSION(:,:)   :: zwz
+
+      CALL wrk_alloc( jpi , jpj      , zwu )
+      CALL wrk_alloc( jpi , jpj      , zwv )
+      CALL wrk_alloc( jpi , jpj      , zwz )
+
+      zmdi=1.e+20                               !  missing data indicator for masking
+      ! Diaopfoam stream if needed
+      IF (ln_diaopfoam) THEN
+         CALL theta2t
+         CALL iom_put( "insitut_op" , rinsitu_t(:,:,:)                      )    ! insitu temperature
+         CALL iom_put( "toce_op"   , tsn(:,:,:,jp_tem)                     )    ! temperature
+         CALL iom_put( "soce_op"   , tsn(:,:,:,jp_sal)                     )    ! salinity
+         IF (ln_diurnal) THEN
+            CALL iom_put( "sst_wl_op"   , x_dsst                    )    ! warm layer 
+            CALL iom_put( "sst_cs_op"   , x_csdsst                  )    ! cool skin 
+         ENDIF
+         zw2d(:,:)=sshn(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
+         CALL iom_put( "ssh_op"  , zw2d(:,:)                          ) ! sea surface height
+         CALL iom_put( "uoce_op"   , un                                    )    ! i-current      
+         CALL iom_put( "voce_op"   , vn                                    )    ! j-current
+         !CALL iom_put( "woce_op"   , wn                                    )    ! k-current
+!#if defined key_spm
+!         cltra = TRIM(ctrc3d(5))//"_op"
+!         zw3d(:,:,:) = trc3d(:,:,:,5)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) ! Visibility
+!         CALL iom_put( cltra, zw3d  )
+!#endif
+         CALL calc_max_cur(zwu,zwv,zwz,zmdi)
+         CALL iom_put( "maxu" , zwu                                     ) ! max u current
+         CALL iom_put( "maxv" , zwv                                     ) ! max v current
+         CALL iom_put( "maxz" , zwz                                     ) ! max current depth
+      ENDIF
+   END SUBROUTINE dia_diaopfoam
+
+   SUBROUTINE calc_max_cur(zmax_u, zmax_v, zmax_z, inmdi)
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE calc_max_cur  ***
+      !!
+      !! ** Purpose :   To locate within the water column the magnitude and
+      !!                vertical location of the strongest horizontal
+      !!                current. The vertical component is ignored since it
+      !!                is an order of magnitude smaller than the horizontal
+      !!                flow, in general. 
+      !!
+      !! ** Method :    A. Map U,V to T-grid.
+      !!                B. Calculate the magnitude of the current for every
+      !!                   grid cell.
+      !!                C. Locate the vertical index using FORTRAN's builtin
+      !!                   MAXLOC function.
+      !!                D. Copy the U,V,Z components of the relevant
+      !!                   indices.
+      !!
+      !! ** Returns :   Value of u, v component and depth of maximum
+      !!                horizontal current on T-grid.
+      !!
+      !!---------------------------------------------------------------------
+      IMPLICIT NONE
+      REAL(wp), DIMENSION(jpi, jpj), INTENT(  OUT) :: zmax_u, zmax_v, zmax_z
+      REAL(wp), DIMENSION(jpi, jpj, jpk)           :: zmax_u_t, zmax_v_t
+      REAL(wp), DIMENSION(jpi, jpj, jpk)           :: zcmag
+      REAL(wp), DIMENSION(jpi, jpj)                :: zmaxk
+      REAL(wp),                      INTENT(IN   ) :: inmdi
+      INTEGER :: ji, jj, jk
+
+      ! Initialise output arrays
+      zmax_u(:, :) = inmdi
+      zmax_v(:, :) = inmdi
+      zmax_z(:, :) = inmdi
+      ! Map to T-grid
+      zmax_u_t(:, :, :) = 0._wp
+      zmax_v_t(:, :, :) = 0._wp
+      DO jk = 1,jpk
+         DO jj = 2,jpj
+            DO ji = 2,jpi
+               zmax_u_t(ji, jj, jk) = 0.5 * (un(ji, jj, jk) + un(ji-1, jj, jk))
+               zmax_v_t(ji, jj, jk) = 0.5 * (vn(ji, jj, jk) + vn(ji, jj-1, jk))
+            END DO
+         END DO
+      END DO
+      ! Calculate absolute velocity
+      zcmag = sqrt((zmax_u_t)**2 + (zmax_v_t)**2)
+      ! Find max. current
+      zmaxk = maxloc(zcmag, dim=3)
+      ! Output values
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+             zmax_u(ji, jj) = zmax_u_t(ji, jj, INT(zmaxk(ji, jj)))
+             zmax_v(ji, jj) = zmax_v_t(ji, jj, INT(zmaxk(ji, jj)))
+             zmax_z(ji, jj) = fsdept(ji, jj, INT(zmaxk(ji, jj)))
+         END DO
+      END DO      
+   END SUBROUTINE calc_max_cur
+
+END MODULE diaopfoam
Index: /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90
===================================================================
--- /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90	(revision 11137)
+++ /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90	(revision 11137)
@@ -0,0 +1,195 @@
+MODULE diatmb 
+   !!======================================================================
+   !!                       ***  MODULE  diaharm  ***
+   !! Harmonic analysis of tidal constituents 
+   !!======================================================================
+   !! History :  3.6  !  2014  (E O'Dea)  Original code
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain
+   USE in_out_manager  ! I/O units
+   USE iom             ! I/0 library
+   USE wrk_nemo        ! working arrays
+#if defined key_fabm
+   USE trc, ONLY: trn, visib
+   USE par_fabm
+   USE fabm, ONLY: fabm_get_bulk_diagnostic_data
+#endif
+
+
+   IMPLICIT NONE
+   PRIVATE
+
+   LOGICAL , PUBLIC ::   ln_diatmb     !: Top Middle and Bottom output 
+   PUBLIC   dia_tmb_init            ! routine called by nemogcm.F90
+   PUBLIC   dia_tmb                 ! routine called by diawri.F90
+   PUBLIC   dia_calctmb             ! routine called by dia25h.F90
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.6 , NEMO Consortium (2014)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dia_tmb_init 
+      !!---------------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_tmb_init  ***
+      !!     
+      !! ** Purpose: Initialization of tmb namelist 
+      !!        
+      !! ** Method : Read namelist
+      !!   History
+      !!   3.6  !  08-14  (E. O'Dea) Routine to initialize dia_tmb
+      !!---------------------------------------------------------------------------
+      !!
+      INTEGER ::   ios                 ! Local integer output status for namelist read
+      !
+      NAMELIST/nam_diatmb/ ln_diatmb
+      !!----------------------------------------------------------------------
+      !
+      REWIND ( numnam_ref )              ! Read Namelist nam_diatmb in reference namelist : TMB diagnostics
+      READ   ( numnam_ref, nam_diatmb, IOSTAT=ios, ERR= 901 )
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diatmb in reference namelist', lwp )
+ 
+      REWIND( numnam_cfg )              ! Namelist nam_diatmb in configuration namelist  TMB diagnostics
+      READ  ( numnam_cfg, nam_diatmb, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diatmb in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, nam_diatmb )
+
+      IF(lwp) THEN                   ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'dia_tmb_init : Output Top, Middle, Bottom Diagnostics'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) 'Namelist nam_diatmb : set tmb outputs '
+         WRITE(numout,*) 'Switch for TMB diagnostics (T) or not (F)  ln_diatmb  = ', ln_diatmb
+      ENDIF
+
+   END SUBROUTINE dia_tmb_init
+
+   SUBROUTINE dia_calctmb( pinfield,pouttmb )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_tmb  ***
+      !!                   
+      !! ** Purpose :    Find the Top, Mid and Bottom fields of water Column
+      !!
+      !! ** Method  :   
+      !!      use mbathy to find surface, mid and bottom of model levels
+      !!
+      !! History :
+      !!   3.6  !  08-14  (E. O'Dea) Routine based on dia_wri_foam
+      !!----------------------------------------------------------------------
+      !! * Modules used
+
+      ! Routine to map 3d field to top, middle, bottom
+      IMPLICIT NONE
+
+
+      ! Routine arguments
+      REAL(wp), DIMENSION(jpi, jpj, jpk), INTENT(IN   ) :: pinfield    ! Input 3d field and mask
+      REAL(wp), DIMENSION(jpi, jpj, 3  ), INTENT(  OUT) :: pouttmb     ! Output top, middle, bottom
+
+
+
+      ! Local variables
+      INTEGER :: ji,jj,jk  ! Dummy loop indices
+
+      ! Local Real
+      REAL(wp)                         ::   zmdi  !  set masked values
+
+      zmdi=1.e+20 !missing data indicator for masking
+
+      ! Calculate top
+      pouttmb(:,:,1) = pinfield(:,:,1)*tmask(:,:,1)  + zmdi*(1.0-tmask(:,:,1))
+
+      ! Calculate middle
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            jk              = max(1,mbathy(ji,jj)/2)
+            pouttmb(ji,jj,2) = pinfield(ji,jj,jk)*tmask(ji,jj,jk)  + zmdi*(1.0-tmask(ji,jj,jk))
+         END DO
+      END DO
+
+      ! Calculate bottom
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            jk              = max(1,mbathy(ji,jj) )
+            pouttmb(ji,jj,3) = pinfield(ji,jj,jk)*tmask(ji,jj,jk)  + zmdi*(1.0-tmask(ji,jj,jk))
+         END DO
+      END DO
+
+   END SUBROUTINE dia_calctmb
+
+
+
+   SUBROUTINE dia_tmb
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE dia_tmb  ***
+      !! ** Purpose :   Write diagnostics for Top, Mid and Bottom of water Column
+      !!
+      !! ** Method  :   
+      !!      use mbathy to find surface, mid and bottom of model levels
+      !!      calls calctmb to retrieve TMB values before sending to iom_put
+      !!
+      !! History :
+      !!   3.6  !  08-14  (E. O'Dea) 
+      !!         
+      !!--------------------------------------------------------------------
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwtmb    ! temporary workspace 
+      REAL(wp)                         ::   zmdi      ! set masked values
+      INTEGER                          ::   jn        ! loop counter
+
+      zmdi=1.e+20 !missing data indicator for maskin
+
+      IF (ln_diatmb) THEN
+         CALL wrk_alloc( jpi , jpj, 3 , zwtmb )
+         CALL dia_calctmb(  tsn(:,:,:,jp_tem),zwtmb )
+         !ssh already output but here we output it masked
+         CALL iom_put( "sshnmasked" , sshn(:,:)*tmask(:,:,1) + zmdi*(1.0 - tmask(:,:,1)) )   ! tmb Temperature
+         CALL iom_put( "top_temp" , zwtmb(:,:,1) )    ! tmb Temperature
+         CALL iom_put( "mid_temp" , zwtmb(:,:,2) )    ! tmb Temperature
+         CALL iom_put( "bot_temp" , zwtmb(:,:,3) )    ! tmb Temperature
+!         CALL iom_put( "sotrefml" , hmld_tref(:,:) )    ! "T criterion Mixed Layer Depth
+
+         CALL dia_calctmb(  tsn(:,:,:,jp_sal),zwtmb )
+         CALL iom_put( "top_sal" , zwtmb(:,:,1) )    ! tmb Salinity 
+         CALL iom_put( "mid_sal" , zwtmb(:,:,2) )    ! tmb Salinity
+         CALL iom_put( "bot_sal" , zwtmb(:,:,3) )    ! tmb Salinity
+
+         CALL dia_calctmb(  un(:,:,:),zwtmb )
+         CALL iom_put( "top_u" , zwtmb(:,:,1) )    ! tmb  U Velocity
+         CALL iom_put( "mid_u" , zwtmb(:,:,2) )    ! tmb  U Velocity
+         CALL iom_put( "bot_u" , zwtmb(:,:,3) )    ! tmb  U Velocity
+!Called in  dynspg_ts.F90        CALL iom_put( "baro_u" , un_b )    ! Barotropic  U Velocity
+
+         CALL dia_calctmb(  vn(:,:,:),zwtmb )
+         CALL iom_put( "top_v" , zwtmb(:,:,1) )    ! tmb  V Velocity
+         CALL iom_put( "mid_v" , zwtmb(:,:,2) )    ! tmb  V Velocity
+         CALL iom_put( "bot_v" , zwtmb(:,:,3) )    ! tmb  V Velocity
+!Called in  dynspg_ts.F90       CALL iom_put( "baro_v" , vn_b )    ! Barotropic  V Velocity
+
+#if defined key_fabm
+         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) )
+         END DO
+         DO jn = 1, jp_fabm_3d
+            CALL dia_calctmb( fabm_get_bulk_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) )
+         END DO
+         CALL dia_calctmb( visib(:,:,:), zwtmb )
+         CALL iom_put( "top_visib" , zwtmb(:,:,1) )
+         CALL iom_put( "mid_visib" , zwtmb(:,:,2) )
+         CALL iom_put( "bot_visib" , zwtmb(:,:,3) )
+#endif
+      ELSE
+         CALL ctl_warn('dia_tmb: tmb diagnostic is set to false you should not have seen this')
+      ENDIF
+
+   END SUBROUTINE dia_tmb
+   !!======================================================================
+END MODULE diatmb
Index: /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIU/cool_skin.F90
===================================================================
--- /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIU/cool_skin.F90	(revision 11137)
+++ /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIU/cool_skin.F90	(revision 11137)
@@ -0,0 +1,149 @@
+MODULE cool_skin
+   !!======================================================================
+   !!                    ***  MODULE  cool_skin  ***
+   !!     Cool skin thickness and delta T correction using Artele et al. (2002)
+   !!     [see also Tu and Tsuang (2005)]
+   !!
+   !!=====================================================================
+   !! History :        !  2012-01  (P. Sykes)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   diurnal_sst_coolskin_step  : timestep the cool skin corrections
+   !!----------------------------------------------------------------------
+   USE par_kind
+   USE phycst
+   USE dom_oce
+   USE in_out_manager
+   USE sbc_oce
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   
+   IMPLICIT NONE
+   
+   !Namelist parameters
+
+   !Parameters
+   REAL(wp), PRIVATE, PARAMETER :: pp_k = 0.596_wp          !Thermal conductivity of seawater
+   REAL(wp), PRIVATE, PARAMETER :: pp_v = 1.05e-6_wp        !Kinematic viscosity of seawater
+   REAL(wp), PRIVATE, PARAMETER :: pp_C = 86400             !seconds [see Tu and Tsuang (2005)]
+   REAL(wp), PRIVATE, PARAMETER :: pp_cw = 3993._wp         !specific heat capacity of seawater
+   REAL(wp), PRIVATE, PARAMETER :: pp_h = 10._wp            !reference depth [using 10m from Artale et al. (2002)]
+   REAL(wp), PRIVATE, PARAMETER :: pp_rhoa = 1.20421_wp     !density of air (at 20C)
+   REAL(wp), PRIVATE, PARAMETER :: pp_cda = 1.45e-3_wp      !assumed air-sea drag coefficient for calculating wind speed
+   
+   !Key variables
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_csdsst     !Cool skin delta SST
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_csthick  !Cool skin thickness
+
+   PRIVATE
+   PUBLIC diurnal_sst_coolskin_step, diurnal_sst_coolskin_init
+
+      !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   
+   CONTAINS 
+   
+   SUBROUTINE diurnal_sst_coolskin_init
+      !!----------------------------------------------------------------------
+      !! *** ROUTINE diurnal_sst_coolskin_init ***
+      !!
+      !! ** Purpose :   initalise the coolskin model
+      !!
+      !! ** Method : 
+      !!
+      !! ** Reference :
+      !! 
+      !!----------------------------------------------------------------------
+      
+      IMPLICIT NONE
+      
+      ALLOCATE( x_csdsst(jpi,jpj), x_csthick(jpi,jpj) )
+      x_csdsst = 0.
+      x_csthick = 0.
+      
+   END SUBROUTINE diurnal_sst_coolskin_init
+ 
+   SUBROUTINE diurnal_sst_coolskin_step(psqflux, pstauflux, psrho, rdt)
+      !!----------------------------------------------------------------------
+      !! *** ROUTINE diurnal_sst_takaya_step ***
+      !!
+      !! ** Purpose :   Timestep the Artale cool skin model
+      !!
+      !! ** Method : 
+      !!
+      !! ** Reference : 
+      !!----------------------------------------------------------------------
+     
+      IMPLICIT NONE
+     
+      !Dummy variables
+      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: psqflux     !Heat (non-solar)(Watts)
+      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: pstauflux   !Wind stress (kg/ m s^2)
+      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: psrho       !Water density (kg/m^3)
+      REAL(wp), INTENT(IN) :: rdt                             !Timestep
+     
+      !Local variables
+      REAL(wp), DIMENSION(jpi,jpj) :: z_fv                    !Friction velocity     
+      REAL(wp), DIMENSION(jpi,jpj) :: z_gamma                 !Dimensionless function of windspeed
+      REAL(wp), DIMENSION(jpi,jpj) :: z_lamda                 !Sauders (dimensionless) proportionality constant
+      REAL(wp), DIMENSION(jpi,jpj) :: z_wspd                  !Windspeed (m/s)
+      REAL(wp) :: z_ztx                                       !Temporary u wind stress
+      REAL(wp) :: z_zty                                       !Temporary v wind stress
+      REAL(wp) :: z_zmod                                      !Temporary total wind stress
+     
+      INTEGER :: ji,jj
+      
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            
+            !Calcualte wind speed from wind stress and friction velocity
+            IF( tmask(ji,jj,1) == 1. .AND. pstauflux(ji,jj) /= 0 .AND. psrho(ji,jj) /=0 ) THEN
+               z_fv(ji,jj) = SQRT( pstauflux(ji,jj) / psrho(ji,jj) )
+               z_wspd(ji,jj) = SQRT( pstauflux(ji,jj) / ( pp_cda * pp_rhoa ) )
+            ELSE
+               z_fv(ji,jj) = 0.
+               z_wspd(ji,jj) = 0.     
+            ENDIF
+
+ 
+            !Calculate gamma function which is dependent upon wind speed
+            IF( tmask(ji,jj,1) == 1. ) THEN
+               IF( ( z_wspd(ji,jj) <= 7.5 ) ) z_gamma(ji,jj) = ( 0.2 * z_wspd(ji,jj) ) + 0.5
+               IF( ( z_wspd(ji,jj) > 7.5 ) .AND. ( z_wspd(ji,jj) < 10. ) ) z_gamma(ji,jj) = ( 1.6 * z_wspd(ji,jj) ) - 10.
+               IF( ( z_wspd(ji,jj) >= 10. ) ) z_gamma(ji,jj) = 6.
+            ENDIF
+
+
+            !Calulate lamda function
+            IF( tmask(ji,jj,1) == 1. .AND. z_fv(ji,jj) /= 0 ) THEN
+               z_lamda(ji,jj) = ( z_fv(ji,jj) * pp_k * pp_C ) / ( z_gamma(ji,jj) * psrho(ji,jj) * pp_cw * pp_h * pp_v )
+            ELSE
+               z_lamda(ji,jj) = 0.
+            ENDIF
+
+
+
+            !Calculate the cool skin thickness - only when heat flux is out of the ocean
+            IF( tmask(ji,jj,1) == 1. .AND. z_fv(ji,jj) /= 0 .AND. psqflux(ji,jj) < 0 ) THEN
+                x_csthick(ji,jj) = ( z_lamda(ji,jj) * pp_v ) / z_fv(ji,jj)
+            ELSE
+                x_csthick(ji,jj) = 0.
+            ENDIF
+
+
+
+            !Calculate the cool skin correction - only when the heat flux is out of the ocean
+            IF( tmask(ji,jj,1) == 1. .AND. x_csthick(ji,jj) /= 0. .AND. psqflux(ji,jj) < 0. ) THEN
+               x_csdsst(ji,jj) = ( psqflux(ji,jj) * x_csthick(ji,jj) ) / pp_k
+             ELSE
+               x_csdsst(ji,jj) = 0.
+            ENDIF
+
+         ENDDO
+      ENDDO
+
+   END SUBROUTINE diurnal_sst_coolskin_step
+
+
+END MODULE cool_skin
Index: /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIU/diurnal_bulk.F90
===================================================================
--- /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIU/diurnal_bulk.F90	(revision 11137)
+++ /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIU/diurnal_bulk.F90	(revision 11137)
@@ -0,0 +1,312 @@
+MODULE diurnal_bulk
+   !!======================================================================
+   !!                    ***  MODULE  diurnal_bulk  ***
+   !!     Takaya model of diurnal warming (Takaya, 2010)
+   !!=====================================================================
+   !! History :        !  11-10  (J. While)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   diurnal_sst_bulk_init  : initalise diurnal model
+   !!   diurnal_sst_bulk_step  : timestep the diurnal model
+   !!----------------------------------------------------------------------
+   USE par_kind
+   USE phycst
+   USE dom_oce
+   USE lib_mpp
+   USE solfrac_mod
+   USE in_out_manager
+   
+   IMPLICIT NONE
+   
+   !Namelist parameters
+   LOGICAL, PUBLIC :: ln_diurnal = .FALSE.
+   LOGICAL, PUBLIC :: ln_diurnal_only = .FALSE.
+
+   !Parameters
+   REAL(wp), PRIVATE, PARAMETER :: pp_alpha = 2.0e-4_wp
+   REAL(wp), PRIVATE, PARAMETER :: pp_veltol = 0._wp
+   REAL(wp), PRIVATE, PARAMETER :: pp_min_fvel = 1.e-10_wp 
+   
+   !Key variables
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_dsst     !Delta SST
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_solfrac  !Fraction of 
+                                                      !absorbed 
+                                                      !radiation
+
+   PRIVATE
+   PUBLIC diurnal_sst_bulk_init, diurnal_sst_takaya_step
+   
+   CONTAINS 
+   
+   SUBROUTINE diurnal_sst_bulk_init
+      !!----------------------------------------------------------------------
+      !! *** ROUTINE diurnal_sst_init ***
+      !!
+      !! ** Purpose : Initalise the Takaya diurnal model
+      
+      !!----------------------------------------------------------------------
+      
+      IMPLICIT NONE
+
+      INTEGER ::   ios                 ! Local integer output status for namelist read
+      
+      NAMELIST /namdiu/ ln_diurnal, ln_diurnal_only
+      
+      !Read the namelist
+      REWIND ( numnam_ref )              ! Read Namelist namdiu in reference namelist
+      READ   ( numnam_ref, namdiu, IOSTAT=ios, ERR= 901 )
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdiu in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namdiu in configuration namelist
+      READ  ( numnam_cfg, namdiu, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdiu in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namdiu )
+     
+      IF(lwp) THEN                   ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'namdiu_init : Diurnal Model settings are:'
+         WRITE(numout,*)    'ln_diurnal = ', ln_diurnal
+         WRITE(numout,*)    'ln_diurnal_only = ', ln_diurnal_only
+      ENDIF
+
+      IF ( ln_diurnal_only .AND. ( .NOT. ln_diurnal ) ) THEN
+         CALL ctl_stop( "ln_diurnal_only set, but ln_diurnal = FALSE !" )
+      ENDIF
+      
+      IF ( ln_diurnal ) THEN      
+         
+         !Allocate arrays
+         ALLOCATE( x_dsst(jpi,jpj), x_solfrac(jpi,jpj) )
+         
+         !Initalise the solar fraction
+         x_solfrac=0._wp
+
+         IF ( ln_diurnal_only ) THEN
+            CALL ctl_warn( "ln_diurnal_only set; only the diurnal component of SST will be calculated" )
+         ENDIF
+      ENDIF
+      
+   END SUBROUTINE diurnal_sst_bulk_init
+   
+   SUBROUTINE diurnal_sst_takaya_step(psolflux, pqflux, ptauflux, prho, p_rdt,&
+            &                  pla, pthick, pcoolthick, pmu, ld_calcfrac, &
+            &                  p_fvel_bkginc, p_hflux_bkginc)
+      !!----------------------------------------------------------------------
+      !! *** ROUTINE diurnal_sst_takaya_step ***
+      !!
+      !! ** Purpose :   Timestep the Takaya diurnal model
+      !!
+      !! ** Method :    1) Calculate the Obukhov length
+      !!                2) Calculate the Simularity function
+      !!                2) Calculate the increment to dsst
+      !!                3) Apply the increment
+      !! ** Reference : Refinements to a prognostic scheme of skin sea surface
+      !!                temperature, Takaya et al, JGR, 2010 
+      !!----------------------------------------------------------------------
+      
+      IMPLICIT NONE
+      
+      !Dummy variables
+      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: psolflux  !solar flux (Watts)
+      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: pqflux    !heat (non-solar) 
+                                                            !flux (Watts)
+      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: ptauflux  !wind stress 
+                                                            !(kg/ m s^2)
+      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: prho      !water density 
+                                                            !(kg/m^3)
+      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: pLa    
+                                                     !Langmuir number
+      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: pthick 
+                                                     !warm layer thickness (m)
+      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: pcoolthick 
+                                                     !cool skin thickness (m)
+      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: pmu    
+                                                     !mu parameter
+      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: p_hflux_bkginc  
+                                                     !optional increment to the
+                                                     !heat flux
+      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: p_fvel_bkginc  
+                                                     !optional increment to the
+                                                     !friction velocity
+      REAL(wp), INTENT(IN) :: p_rdt              !timestep
+      
+      LOGICAL, OPTIONAL, &
+      &     INTENT(IN)  :: ld_calcfrac !Set TRUE to recalculate the 
+                                               !solar fraction
+      
+      !Local variables 
+      REAL(wp), DIMENSION(jpi,jpj) :: z_fvel              ! friction velocity     
+      REAL(wp), DIMENSION(jpi,jpj) :: zthick, zcoolthick, zmu, zla
+      REAL(wp), DIMENSION(jpi,jpj) :: z_abflux            ! absorbed flux           
+      REAL(wp), DIMENSION(jpi,jpj) :: z_fla               !Langmuir fuction value 
+      
+      LOGICAL  :: ll_calcfrac
+      
+      INTEGER :: ji,jj
+
+      !Set optional arguments to their defualts
+      IF ( .NOT. PRESENT(pthick) ) THEN 
+         zthick(:,:) = 3._wp
+      ELSE
+         zthick(:,:) = pthick(:,:)
+      ENDIF
+      IF ( .NOT. PRESENT(pcoolthick) ) THEN 
+         zcoolthick(:,:) = 0._wp
+      ELSE
+         zcoolthick(:,:) = pcoolthick(:,:)
+      ENDIF
+      IF ( .NOT. PRESENT(pmu) ) THEN
+         zmu(:,:) = 0.3_wp
+      ELSE
+         zmu(:,:) = pmu(:,:)
+      ENDIF
+      IF ( .NOT. PRESENT(pla) ) THEN
+         zla(:,:) = 0.3_wp
+      ELSE
+         zla(:,:) = pla(:,:)
+      ENDIF
+      IF ( .NOT. PRESENT(ld_calcfrac) ) THEN
+         ll_calcfrac = .FALSE.
+      ELSE
+         ll_calcfrac = ld_calcfrac
+      ENDIF     
+      
+      !If not done already, calculate the solar fraction
+      IF (ll_calcfrac ) THEN
+         DO jj = 1,jpj
+            DO ji = 1, jpi
+               IF(  ( x_solfrac(ji,jj) == 0. ) .AND. ( tmask(ji,jj,1) == 1. ) ) &
+               &   x_solfrac(ji,jj) = solfrac( zcoolthick(ji,jj),zthick(ji,jj) )
+            END DO
+         END DO   
+      ENDIF   
+
+      !convert solar flux and heat flux to absorbed flux   
+      WHERE ( tmask(:,:,1) == 1.) 
+         z_abflux(:,:) = (  x_solfrac(:,:) *  psolflux (:,:)) + pqflux(:,:)     
+      ELSEWHERE
+         z_abflux(:,:) = 0._wp
+      ENDWHERE
+      IF( PRESENT(p_hflux_bkginc) ) z_abflux(:,:) = z_abflux(:,:) + p_hflux_bkginc !Optional increment
+      WHERE ( ABS( z_abflux(:,:) ) < rsmall )
+         z_abflux(:,:) = rsmall
+      ENDWHERE 
+      
+      !Calculate the friction velocity
+      WHERE ( (ptauflux /= 0) .AND. ( tmask(:,:,1) == 1.) )   
+         z_fvel(:,:) = SQRT( ptauflux(:,:) / prho(:,:) )
+      ELSEWHERE
+         z_fvel(:,:) = 0.
+      ENDWHERE
+      IF( PRESENT(p_hflux_bkginc) ) z_fvel(:,:) = z_fvel(:,:) + p_fvel_bkginc !Optional increment
+      
+       
+       
+      !Calculate the Langmuir function value
+      WHERE ( tmask(:,:,1) == 1.)
+         z_fla(:,:) = MAX( 1._wp, zla(:,:)**( -2._wp / 3._wp ) )  
+      ENDWHERE     
+      
+      !!Increment the temperature using the implicit solution
+      x_dsst(:,:) = t_imp( x_dsst(:,:), p_rdt, z_abflux(:,:), z_fvel(:,:), &
+                                 z_fla(:,:), zmu(:,:), zthick(:,:), prho(:,:))
+      
+                                
+      
+   END SUBROUTINE diurnal_sst_takaya_step
+
+   
+   FUNCTION t_imp(p_dsst, p_rdt, p_abflux, p_fvel, &
+                          p_fla, pmu, pthick, prho )
+                          
+      IMPLICIT NONE
+      
+      !Function definition
+      REAL(wp), DIMENSION(jpi,jpj) :: t_imp
+      !Dummy variables
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_dsst    !Delta SST
+      REAL(wp), INTENT(IN) ::                 p_rdt           !Timestep
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_abflux  !Heat forcing
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fvel    !Friction velocity
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fla     !Langmuir number
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pmu       !Structure parameter
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pthick    !Layer thickness
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: prho      !Water density
+   
+      !Local variables
+      REAL(wp) :: z_olength          ! Obukhov length
+      REAL(wp) :: z_sigma, z_sigma2
+      REAL(wp) :: z_term1, z_term2     
+      REAL(wp) :: z_stabfunc         ! stability function value
+      REAL(wp) :: z_fvel      
+      
+      CHARACTER(200) :: warn_string
+       
+      INTEGER :: ji,jj
+                     
+      DO jj = 1, jpj
+         DO ji = 1, jpi   
+            
+            !Only calculate outside tmask
+            IF ( tmask(ji,jj,1) /= 1. ) THEN
+               t_imp(ji,jj) = 0._wp
+               CYCLE   
+            END IF
+            
+            IF (p_fvel(ji,jj) < pp_min_fvel) THEN
+               z_fvel = pp_min_fvel
+               WRITE(warn_string,*) "diurnal_sst_takaya step: "&
+               &//"friction velocity < minimum\n" &
+               &//"Setting friction velocity =",pp_min_fvel
+               CALL ctl_warn(warn_string)
+               
+            ELSE
+               z_fvel = p_fvel(ji,jj)
+            ENDIF
+                 
+            !Calculate the Obukhov length
+            IF ( (z_fvel < pp_veltol ) .AND. &
+            &    (p_dsst(ji,jj) > 0.) ) THEN
+               z_olength =  z_fvel  / &
+               &     SQRT( p_dsst(ji,jj) * vkarmn * grav * &
+               &             pp_alpha / ( 5._wp * pthick(ji,jj) ) )
+            ELSE
+               z_olength = &
+               &   ( prho(ji,jj) * rcp * z_fvel**3._wp ) / &
+               &   ( vkarmn * grav * pp_alpha *&
+               &     p_abflux(ji,jj) )
+            ENDIF
+             
+            !Calculate the stability function          
+            z_sigma = pthick(ji,jj) / z_olength
+            z_sigma2 = z_sigma * z_sigma
+  
+            IF ( z_sigma >= 0 ) THEN
+               z_stabfunc = 1._wp + &
+               &  ( ( 5._wp * z_sigma + 4._wp * z_sigma2 ) / &
+               &  ( 1._wp + 3._wp * z_sigma + 0.25_wp * &
+               &    z_sigma2 ) )
+            ELSE
+               z_stabfunc = 1._wp / &
+               &                   SQRT( 1._wp - 16._wp * z_sigma )
+            ENDIF
+
+            !Calculate the T increment
+            z_term1 = ( p_abflux(ji,jj) * ( pmu(ji,jj) + 1._wp)  / &
+            & ( pmu(ji,jj) * pthick(ji,jj) * prho(ji,jj) * rcp ) )
+            
+             
+            z_term2 = -( ( pmu(ji,jj) + 1._wp) * &
+            &                       ( vkarmn * z_fvel * p_fla(ji,jj) ) / &
+            &      ( pthick(ji,jj) * z_stabfunc ) )     
+          
+            t_imp(ji,jj) = ( p_dsst(ji,jj) + p_rdt * z_term1 ) / &
+                           ( 1._wp - p_rdt * z_term2 )
+
+         END DO
+      END DO
+      
+      END FUNCTION t_imp
+
+END MODULE diurnal_bulk
Index: /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIU/solfrac_mod.F90
===================================================================
--- /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIU/solfrac_mod.F90	(revision 11137)
+++ /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIU/solfrac_mod.F90	(revision 11137)
@@ -0,0 +1,56 @@
+MODULE solfrac_mod
+   !!======================================================================
+   !!                    ***  MODULE  solfrac  ***
+   !!     POSH representation of solar absorbtion (Gntermann, 2009)
+   !!=====================================================================
+   !! History :        !  11-10  (J. While)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   solfrac  : function to calculate the solar fraction
+   !!----------------------------------------------------------------------
+   
+   USE par_kind
+   IMPLICIT NONE
+     
+   !Parameters 
+   REAL(wp), PRIVATE, PARAMETER, DIMENSION(9) :: &
+   &                                     pp_wgt = (/0.2370, 0.36,  0.1790, &
+   &                                                0.087,  0.08,  0.025,  &
+   &                                                0.025,  0.007, 0.0004/)
+   REAL(wp), PRIVATE, PARAMETER, DIMENSION(9) :: &
+   &                                    pp_len = (/34.84,   2.266,   0.0315,  &
+   &                                               0.0055,  8.32e-4, 1.26e-4, &
+   &                                                3.13e-4, 7.82e-4, 1.44e-5/)
+   
+   PUBLIC solfrac
+                                                   
+CONTAINS
+
+   REAL(dp) FUNCTION solfrac(ptop,pbottom)
+       !!----------------------------------------------------------------------
+      !! *** ROUTINE solfrac ***
+      !!
+      !! ** Purpose :   Calculate the solar fraction absorbed between two 
+      !!                layers
+      !!
+      !! ** Reference : POSH a model of diurnal warming, Gentemann et al, 
+      !!                 JGR, 2009 
+      !!----------------------------------------------------------------------
+      
+      !Dummy variabes
+      REAL(wp), INTENT(IN) :: ptop, pbottom  !Top and bottom of layer
+      
+      !local variables
+      INTEGER :: jt
+      
+      !Calculate the solar fraction absorbed between the two layers
+      solfrac = 0._wp
+      DO jt = 1, 9
+           solfrac = solfrac + pp_wgt(jt) * ( exp ( -ptop / pp_len(jt) ) &
+            &                                 - exp ( -pbottom / pp_len(jt) ) )
+      END DO
+      
+   END FUNCTION
+   
+END MODULE solfrac_mod
Index: /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/inputs_fabm.F90
===================================================================
--- /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/inputs_fabm.F90	(revision 11137)
+++ /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/inputs_fabm.F90	(revision 11137)
@@ -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/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/par_fabm.F90
===================================================================
--- /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/par_fabm.F90	(revision 11137)
+++ /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/par_fabm.F90	(revision 11137)
@@ -0,0 +1,46 @@
+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, jp_fabm_2d, jp_fabm_3d
+
+   ! Variables needed for OBS/ASM
+   INTEGER, PUBLIC :: jp_fabm_chl1, jp_fabm_chl2, &
+                      jp_fabm_chl3, jp_fabm_chl4, &
+                      jp_fabm_p1c,  jp_fabm_p1n,  &
+                      jp_fabm_p1p,  jp_fabm_p1s,  &
+                      jp_fabm_p2c,  jp_fabm_p2n,  &
+                      jp_fabm_p2p,  jp_fabm_p3c,  &
+                      jp_fabm_p3n,  jp_fabm_p3p,  &
+                      jp_fabm_p4c,  jp_fabm_p4n,  &
+                      jp_fabm_p4p,  jp_fabm_z4c,  &
+                      jp_fabm_z5c,  jp_fabm_z5n,  &
+                      jp_fabm_z5p,  jp_fabm_z6c,  &
+                      jp_fabm_z6n,  jp_fabm_z6p,  &
+                      jp_fabm_n1p,  jp_fabm_n3n,  &
+                      jp_fabm_n4n,  jp_fabm_n5s,  &
+                      jp_fabm_o2o,  jp_fabm_o3c,  &
+                      jp_fabm_o3a,  jp_fabm_o3ph, &
+                      jp_fabm_o3pc, jp_fabm_xeps
+
+#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/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/st2D_fabm.F90
===================================================================
--- /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/st2D_fabm.F90	(revision 11137)
+++ /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/st2D_fabm.F90	(revision 11137)
@@ -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/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/trcini_fabm.F90
===================================================================
--- /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/trcini_fabm.F90	(revision 11137)
+++ /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/trcini_fabm.F90	(revision 11137)
@@ -0,0 +1,509 @@
+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, &
+                  fabm_get_bulk_diagnostic_data
+   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
+      jp_fabm_2d = size(model%horizontal_diagnostic_variables)
+      jp_fabm_3d = size(model%diagnostic_variables)
+      jpdia2d = jpdia2d + jp_fabm_2d
+      jpdia3d = jpdia3d + jp_fabm_3d
+      jpdiabio = jpdiabio + jp_fabm
+
+      !Initialize input data structures.
+      call initialize_inputs
+
+      ! Get indexes for select state variables
+      jp_fabm_chl1 = fabm_state_index( 'P1_Chl' )
+      jp_fabm_chl2 = fabm_state_index( 'P2_Chl' )
+      jp_fabm_chl3 = fabm_state_index( 'P3_Chl' )
+      jp_fabm_chl4 = fabm_state_index( 'P4_Chl' )
+      jp_fabm_p1c  = fabm_state_index( 'P1_c' )
+      jp_fabm_p1n  = fabm_state_index( 'P1_n' )
+      jp_fabm_p1p  = fabm_state_index( 'P1_p' )
+      jp_fabm_p1s  = fabm_state_index( 'P1_s' )
+      jp_fabm_p2c  = fabm_state_index( 'P2_c' )
+      jp_fabm_p2n  = fabm_state_index( 'P2_n' )
+      jp_fabm_p2p  = fabm_state_index( 'P2_p' )
+      jp_fabm_p3c  = fabm_state_index( 'P3_c' )
+      jp_fabm_p3n  = fabm_state_index( 'P3_n' )
+      jp_fabm_p3p  = fabm_state_index( 'P3_p' )
+      jp_fabm_p4c  = fabm_state_index( 'P4_c' )
+      jp_fabm_p4n  = fabm_state_index( 'P4_n' )
+      jp_fabm_p4p  = fabm_state_index( 'P4_p' )
+      jp_fabm_z4c  = fabm_state_index( 'Z4_c' )
+      jp_fabm_z5c  = fabm_state_index( 'Z5_c' )
+      jp_fabm_z5n  = fabm_state_index( 'Z5_n' )
+      jp_fabm_z5p  = fabm_state_index( 'Z5_p' )
+      jp_fabm_z6c  = fabm_state_index( 'Z6_c' )
+      jp_fabm_z6n  = fabm_state_index( 'Z6_n' )
+      jp_fabm_z6p  = fabm_state_index( 'Z6_p' )
+      jp_fabm_n1p  = fabm_state_index( 'N1_p' )
+      jp_fabm_n3n  = fabm_state_index( 'N3_n' )
+      jp_fabm_n4n  = fabm_state_index( 'N4_n' )
+      jp_fabm_n5s  = fabm_state_index( 'N5_s' )
+      jp_fabm_o2o  = fabm_state_index( 'O2_o' )
+      jp_fabm_o3c  = fabm_state_index( 'O3_c' )
+      jp_fabm_o3a  = fabm_state_index( 'O3_bioalk' )
+
+      ! Get indexes for select diagnostic variables
+      jp_fabm_o3ph = fabm_diag_index( 'O3_pH' )
+      jp_fabm_o3pc = fabm_diag_index( 'O3_pCO2' )
+      jp_fabm_xeps = fabm_diag_index( 'light_xEPS')
+
+      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
+            CALL write_25hourm_xml(xml_unit,model%state_variables(jn))
+            CALL write_tmb_xml(xml_unit,model%state_variables(jn))
+         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))
+            CALL write_25hourm_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))
+            CALL write_25hourm_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)
+            CALL write_25hourm_xml(xml_unit,model%diagnostic_variables(jn),3)
+            CALL write_tmb_xml(xml_unit,model%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
+         WRITE (xml_unit,'(A)') '  <field id="visib" long_name="visibility" unit="m" default_value="1.e+20" grid_ref="grid_T_3D" />'
+         WRITE (xml_unit,'(A)') '  <field id="visib25h" long_name="visibility 25-hour mean" unit="m" default_value="1.e+20" grid_ref="grid_T_3D" />'
+         WRITE (xml_unit,'(A)') '  <field id="top_visib" long_name="Top-level visibility" unit="m" default_value="1.e+20" grid_ref="grid_T_2D" />'
+         WRITE (xml_unit,'(A)') '  <field id="mid_visib" long_name="Middle-level visibility" unit="m" default_value="1.e+20" grid_ref="grid_T_2D" />'
+         WRITE (xml_unit,'(A)') '  <field id="bot_visib" long_name="Bottom-level visibility" unit="m" default_value="1.e+20" grid_ref="grid_T_2D" />'
+         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_25hourm_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)') 1.e+20
+      WRITE (string_dimensions,'(I1)') number_dimensions
+      SELECT CASE (number_dimensions)
+      CASE (3)
+         WRITE (xml_unit,'(A)') '  <field id="'//TRIM(variable%name)//'25h'//'" long_name="'//TRIM(variable%long_name)//' 25-hour mean'//'" 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)//'25h'//'" long_name="'//TRIM(variable%long_name)//' 25-hour mean'//'" 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)//'25h'//'" long_name="'//TRIM(variable%long_name)//' 25-hour mean'//'" unit="'//TRIM(variable%units)//'" default_value="'//TRIM(ADJUSTL(missing_value))//'" grid_ref="1point"/>'
+      CASE (-1)
+         WRITE (xml_unit,'(A)') '  <field id="'//TRIM(variable%name)//'25h'//'" long_name="'//TRIM(variable%long_name)//' 25-hour mean'//'" 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)//'25h'//': Output of '//TRIM(ADJUSTL(string_dimensions))//'-dimensional variables not supported!!!'
+      END SELECT
+
+   END SUBROUTINE write_25hourm_xml
+
+   SUBROUTINE write_tmb_xml(xml_unit,variable)
+      INTEGER,INTENT(IN) :: xml_unit
+      CLASS (type_external_variable),INTENT(IN) :: variable
+
+      CHARACTER(LEN=20) :: missing_value
+
+      WRITE (missing_value,'(E9.3)') 1.e+20
+      WRITE (xml_unit,'(A)') '  <field id="top_'//TRIM(variable%name)//'" long_name="Top-level '//TRIM(variable%long_name)//'" unit="'//TRIM(variable%units)//'" default_value="'//TRIM(ADJUSTL(missing_value))//'" grid_ref="grid_T_2D"/>'
+      WRITE (xml_unit,'(A)') '  <field id="mid_'//TRIM(variable%name)//'" long_name="Middle-level '//TRIM(variable%long_name)//'" unit="'//TRIM(variable%units)//'" default_value="'//TRIM(ADJUSTL(missing_value))//'" grid_ref="grid_T_2D"/>'
+      WRITE (xml_unit,'(A)') '  <field id="bot_'//TRIM(variable%name)//'" long_name="Bottom-level '//TRIM(variable%long_name)//'" unit="'//TRIM(variable%units)//'" default_value="'//TRIM(ADJUSTL(missing_value))//'" grid_ref="grid_T_2D"/>'
+
+   END SUBROUTINE write_tmb_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
+
+      ! Initialise visibility
+      visib(:,:,:) = 1.7 / fabm_get_bulk_diagnostic_data(model, jp_fabm_xeps)
+
+      ! 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
+
+   INTEGER FUNCTION fabm_state_index( state_name )
+      !!----------------------------------------------------------------------
+      !!                     ***  fabm_state_index  ***  
+      !!
+      !! ** Purpose :   return index of a given FABM state variable
+      !!
+      !! ** Method  : - loop through state variables until found
+      !!----------------------------------------------------------------------
+      
+      IMPLICIT NONE
+      
+      CHARACTER(LEN=256), INTENT(IN) :: state_name
+      
+      INTEGER                        :: jn
+
+      !!----------------------------------------------------------------------
+      
+      fabm_state_index = -1
+      DO jn=1,jp_fabm
+         IF (TRIM(model%state_variables(jn)%name) == TRIM(state_name)) THEN
+            fabm_state_index = jn
+            EXIT
+         ENDIF
+      END DO
+      IF (fabm_state_index == -1) THEN
+         CALL ctl_stop( 'Could not find '//TRIM(state_name)//' state variable' )
+      ELSE
+         IF (lwp) WRITE(numout,*) 'Index for '//TRIM(state_name)//' is: ', fabm_state_index
+      ENDIF
+   
+   END FUNCTION fabm_state_index
+
+   INTEGER FUNCTION fabm_diag_index( diag_name )
+      !!----------------------------------------------------------------------
+      !!                     ***  fabm_state_index  ***  
+      !!
+      !! ** Purpose :   return index of a given FABM diagnostic variable
+      !!
+      !! ** Method  : - loop through diagnostic variables until found
+      !!----------------------------------------------------------------------
+      
+      IMPLICIT NONE
+      
+      CHARACTER(LEN=256), INTENT(IN) :: diag_name
+      
+      INTEGER                        :: jn
+
+      !!----------------------------------------------------------------------
+      
+      fabm_diag_index = -1
+      DO jn = 1, SIZE(model%diagnostic_variables)
+         IF (TRIM(model%diagnostic_variables(jn)%name) == TRIM(diag_name)) THEN
+            fabm_diag_index = jn
+            EXIT
+         ENDIF
+      END DO
+      IF (fabm_diag_index == -1) THEN
+         CALL ctl_stop( 'Could not find '//TRIM(diag_name)//' diagnostic' )
+      ELSE
+         IF (lwp) WRITE(numout,*) 'Index for '//TRIM(diag_name)//' is: ', fabm_diag_index
+      ENDIF
+   
+   END FUNCTION fabm_diag_index
+
+#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/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/trcnam_fabm.F90
===================================================================
--- /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/trcnam_fabm.F90	(revision 11137)
+++ /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/trcnam_fabm.F90	(revision 11137)
@@ -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/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/trcrst_fabm.F90
===================================================================
--- /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/trcrst_fabm.F90	(revision 11137)
+++ /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/trcrst_fabm.F90	(revision 11137)
@@ -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/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90
===================================================================
--- /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90	(revision 11137)
+++ /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90	(revision 11137)
@@ -0,0 +1,501 @@
+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
+   USE par_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
+      
+      visib(:,:,:) = 1.7 / fabm_get_bulk_diagnostic_data(model, jp_fabm_xeps)
+
+      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/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/trcwri_fabm.F90
===================================================================
--- /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/trcwri_fabm.F90	(revision 11137)
+++ /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/trcwri_fabm.F90	(revision 11137)
@@ -0,0 +1,138 @@
+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
+      
+      REAL(wp), DIMENSION(jpi,jpj,jpk)    :: zw3d
+
+#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
+      zw3d(:,:,:) = visib(:,:,:)*tmask(:,:,:) + 1.e+20*(1.0-tmask(:,:,:))
+      CALL iom_put( 'visib', zw3d )
+
+      ! 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$
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!======================================================================
+END MODULE trcwri_fabm
Index: /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90
===================================================================
--- /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90	(revision 11137)
+++ /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90	(revision 11137)
@@ -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/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/trcbdy.F90
===================================================================
--- /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/trcbdy.F90	(revision 11137)
+++ /branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/trcbdy.F90	(revision 11137)
@@ -0,0 +1,305 @@
+MODULE trcbdy
+   !!======================================================================
+   !!                       ***  MODULE  bdytrc  ***
+   !! Ocean tracers:   Apply boundary conditions for tracers in TOP component
+   !!======================================================================
+   !! History :  1.0  !  2005-01  (J. Chanut, A. Sellar)  Original code
+   !!            3.0  !  2008-04  (NEMO team)  add in the reference version
+   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge
+   !!            3.5  !  2012     (S. Mocavero, I. Epicoco) Optimization of BDY communications
+   !!            3.6  !  2015     (T. Lovato) Adapt BDY for tracers in TOP component
+   !!----------------------------------------------------------------------
+#if defined key_bdy && key_top
+   !!----------------------------------------------------------------------
+   !!   'key_bdy'                     Unstructured Open Boundary Conditions
+   !!----------------------------------------------------------------------
+   !!   trc_bdy            : Apply open boundary conditions to T and S
+   !!   trc_bdy_frs        : Apply Flow Relaxation Scheme
+   !!----------------------------------------------------------------------
+   USE timing                       ! Timing
+   USE oce_trc                      ! ocean dynamics and tracers variables
+   USE par_trc
+   USE trc                          ! ocean space and time domain variables 
+   USE bdylib                       ! for orlanski library routines
+   USE lbclnk                       ! ocean lateral boundary conditions (or mpp link)
+   USE in_out_manager               ! I/O manager
+   USE bdy_oce, only: idx_bdy, OBC_INDEX, BDYTMASK, lk_bdy       ! ocean open boundary conditions
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC trc_bdy      ! routine called in trcnxt.F90 
+   PUBLIC trc_bdy_dmp  ! routine called in trcstp.F90 
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.6 , NEMO Consortium (2015)
+   !! $Id$
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE trc_bdy( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  SUBROUTINE trc_bdy  ***
+      !!
+      !! ** Purpose : - Apply open boundary conditions for tracers in TOP component
+      !!                and scale the tracer data
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) :: kt     ! Main time step counter
+      !!
+      INTEGER               :: ib_bdy, jn ! Loop indeces
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('trc_bdy')
+      !
+      DO jn = 1, jptra
+         DO ib_bdy=1, nb_bdy
+
+            SELECT CASE( trcdta_bdy(jn,ib_bdy)%cn_obc )
+            CASE('none')
+               CYCLE
+            CASE('frs')
+               CALL bdy_trc_frs( jn, idx_bdy(ib_bdy), trcdta_bdy(jn,ib_bdy), kt )
+            CASE('specified')
+               CALL bdy_trc_spe( jn, idx_bdy(ib_bdy), trcdta_bdy(jn,ib_bdy), kt )
+            CASE('neumann')
+               CALL bdy_trc_nmn( jn, idx_bdy(ib_bdy), trcdta_bdy(jn,ib_bdy), kt )
+            CASE('orlanski')
+               CALL bdy_trc_orlanski( jn, idx_bdy(ib_bdy), trcdta_bdy(jn,ib_bdy), ll_npo=.false. )
+            CASE('orlanski_npo')
+               CALL bdy_trc_orlanski( jn, idx_bdy(ib_bdy), trcdta_bdy(jn,ib_bdy), ll_npo=.true. )
+            CASE DEFAULT
+               CALL ctl_stop( 'bdy_tra : unrecognised option for open boundaries for T and S' )
+            END SELECT
+
+            ! Boundary points should be updated
+            CALL lbc_bdy_lnk( tra(:,:,:,jn), 'T', 1., ib_bdy )
+
+         ENDDO
+      ENDDO
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('trc_bdy')
+
+   END SUBROUTINE trc_bdy
+
+   SUBROUTINE bdy_trc_frs( jn, idx, dta, kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_trc_frs  ***
+      !!                    
+      !! ** Purpose : Apply the Flow Relaxation Scheme for tracers at open boundaries.
+      !! 
+      !! Reference : Engedahl H., 1995, Tellus, 365-382.
+      !!----------------------------------------------------------------------
+      INTEGER,         INTENT(in) ::   kt
+      INTEGER,         INTENT(in) ::   jn   ! Tracer index
+      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
+      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
+      !! 
+      REAL(wp) ::   zwgt           ! boundary weight
+      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
+      INTEGER  ::   ii, ij         ! 2D addresses
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_trc_frs')
+      ! 
+      igrd = 1                       ! Everything is at T-points here
+      DO ib = 1, idx%nblen(igrd)
+         DO ik = 1, jpkm1
+            ii = idx%nbi(ib,igrd)
+            ij = idx%nbj(ib,igrd)
+            zwgt = idx%nbw(ib,igrd)
+            tra(ii,ij,ik,jn) = ( tra(ii,ij,ik,jn) + zwgt * ( ( dta%trc(ib,ik) * dta%rn_fac)  & 
+                        &  - tra(ii,ij,ik,jn) ) ) * tmask(ii,ij,ik)
+         END DO
+      END DO 
+      !
+      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_trc_frs')
+      !
+   END SUBROUTINE bdy_trc_frs
+  
+   SUBROUTINE bdy_trc_spe( jn, idx, dta, kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_trc_frs  ***
+      !!                    
+      !! ** Purpose : Apply a specified value for tracers at open boundaries.
+      !! 
+      !!----------------------------------------------------------------------
+      INTEGER,         INTENT(in) ::   kt
+      INTEGER,         INTENT(in) ::   jn   ! Tracer index
+      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
+      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
+      !! 
+      REAL(wp) ::   zwgt           ! boundary weight
+      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
+      INTEGER  ::   ii, ij         ! 2D addresses
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_trc_spe')
+      !
+      igrd = 1                       ! Everything is at T-points here
+      DO ib = 1, idx%nblenrim(igrd)
+         ii = idx%nbi(ib,igrd)
+         ij = idx%nbj(ib,igrd)
+         DO ik = 1, jpkm1
+            tra(ii,ij,ik,jn) = dta%trc(ib,ik) * dta%rn_fac * tmask(ii,ij,ik)
+         END DO
+      END DO
+      !
+      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_trc_spe')
+      !
+   END SUBROUTINE bdy_trc_spe
+
+   SUBROUTINE bdy_trc_nmn( jn, idx, dta, kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_trc_nmn  ***
+      !!                    
+      !! ** Purpose : Duplicate the value for tracers at open boundaries.
+      !! 
+      !!----------------------------------------------------------------------
+      INTEGER,         INTENT(in) ::   kt
+      INTEGER,         INTENT(in) ::   jn   ! Tracer index
+      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
+      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
+      !! 
+      REAL(wp) ::   zcoef, zcoef1, zcoef2           ! boundary weight
+      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
+      INTEGER  ::   ii, ij, ip, jp   ! 2D addresses
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_trc_nmn')
+      !
+      igrd = 1                       ! Everything is at T-points here
+      DO ib = 1, idx%nblenrim(igrd)
+         ii = idx%nbi(ib,igrd)
+         ij = idx%nbj(ib,igrd)
+         DO ik = 1, jpkm1
+            ! search the sense of the gradient
+            zcoef1 = bdytmask(ii-1,ij  )*tmask(ii-1,ij,ik) +  bdytmask(ii+1,ij  )*tmask(ii+1,ij,ik)
+            zcoef2 = bdytmask(ii  ,ij-1)*tmask(ii,ij-1,ik) +  bdytmask(ii  ,ij+1)*tmask(ii,ij+1,ik)
+            IF ( nint(zcoef1+zcoef2) == 0) THEN
+               ! corner
+               zcoef = tmask(ii-1,ij,ik) + tmask(ii+1,ij,ik) +  tmask(ii,ij-1,ik) +  tmask(ii,ij+1,ik)
+               IF (zcoef > .5_wp) THEN ! Only set not isolated points.
+                 tra(ii,ij,ik,jn) = tra(ii-1,ij  ,ik,jn) * tmask(ii-1,ij  ,ik) + &
+                   &              tra(ii+1,ij  ,ik,jn) * tmask(ii+1,ij  ,ik) + &
+                   &              tra(ii  ,ij-1,ik,jn) * tmask(ii  ,ij-1,ik) + &
+                   &              tra(ii  ,ij+1,ik,jn) * tmask(ii  ,ij+1,ik)
+                 tra(ii,ij,ik,jn) = ( tra(ii,ij,ik,jn) / zcoef ) * tmask(ii,ij,ik)
+               ENDIF
+            ELSEIF ( nint(zcoef1+zcoef2) == 2) THEN
+               ! oblique corner
+               zcoef = tmask(ii-1,ij,ik)*bdytmask(ii-1,ij  ) + tmask(ii+1,ij,ik)*bdytmask(ii+1,ij  ) + &
+                  &  tmask(ii,ij-1,ik)*bdytmask(ii,ij -1 ) +  tmask(ii,ij+1,ik)*bdytmask(ii,ij+1  )
+               tra(ii,ij,ik,jn) = tra(ii-1,ij  ,ik,jn) * tmask(ii-1,ij  ,ik)*bdytmask(ii-1,ij  ) + &
+                  &              tra(ii+1,ij  ,ik,jn) * tmask(ii+1,ij  ,ik)*bdytmask(ii+1,ij  )  + &
+                  &              tra(ii  ,ij-1,ik,jn) * tmask(ii  ,ij-1,ik)*bdytmask(ii,ij -1 ) + &
+                  &              tra(ii  ,ij+1,ik,jn) * tmask(ii  ,ij+1,ik)*bdytmask(ii,ij+1  )
+ 
+               tra(ii,ij,ik,jn) = ( tra(ii,ij,ik,jn) / MAX(1._wp, zcoef) ) * tmask(ii,ij,ik)
+            ELSE
+               ip = nint(bdytmask(ii+1,ij  )*tmask(ii+1,ij,ik) - bdytmask(ii-1,ij  )*tmask(ii-1,ij,ik))
+               jp = nint(bdytmask(ii  ,ij+1)*tmask(ii,ij+1,ik) - bdytmask(ii  ,ij-1)*tmask(ii,ij-1,ik))
+               tra(ii,ij,ik,jn) = tra(ii+ip,ij+jp,ik,jn) * tmask(ii+ip,ij+jp,ik) * tmask(ii,ij,ik)
+            ENDIF
+         END DO
+      END DO
+      !
+      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_trc_nmn')
+      !
+   END SUBROUTINE bdy_trc_nmn
+ 
+
+   SUBROUTINE bdy_trc_orlanski( jn, idx, dta, ll_npo )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_trc_orlanski  ***
+      !!             
+      !!              - Apply Orlanski radiation to tracers of TOP component. 
+      !!              - Wrapper routine for bdy_orlanski_3d
+      !! 
+      !!
+      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)    
+      !!----------------------------------------------------------------------
+      INTEGER,                      INTENT(in) ::   jn      ! Tracer index
+      TYPE(OBC_INDEX),              INTENT(in) ::   idx     ! OBC indices
+      TYPE(OBC_DATA),               INTENT(in) ::   dta     ! OBC external data
+      LOGICAL,                      INTENT(in) ::   ll_npo  ! switch for NPO version
+
+      INTEGER  ::   igrd                                    ! grid index
+      !!----------------------------------------------------------------------
+
+      IF( nn_timing == 1 ) CALL timing_start('bdy_trc_orlanski')
+      !
+      igrd = 1      ! Orlanski bc on tracers; 
+      !            
+      CALL bdy_orlanski_3d( idx, igrd, trb(:,:,:,jn), tra(:,:,:,jn), (dta%trc * dta%rn_fac), ll_npo )
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_trc_orlanski')
+      !
+
+   END SUBROUTINE bdy_trc_orlanski
+
+   SUBROUTINE trc_bdy_dmp( kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE trc_bdy_dmp  ***
+      !!                    
+      !! ** Purpose : Apply damping for tracers at open boundaries.
+      !!             It currently applies the damping to all tracers!!!
+      !! 
+      !!----------------------------------------------------------------------
+      INTEGER,         INTENT(in) ::   kt
+      !! 
+      INTEGER  ::   jn             ! Tracer index
+      REAL(wp) ::   zwgt           ! boundary weight
+      REAL(wp) ::   zta, zsa, ztime
+      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
+      INTEGER  ::   ii, ij         ! 2D addresses
+      INTEGER  ::   ib_bdy         ! Loop index
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('trc_bdy_dmp')
+      !
+      DO jn = 1, jptra
+         DO ib_bdy=1, nb_bdy
+            IF ( trcdta_bdy(jn, ib_bdy)%dmp ) THEN
+               igrd = 1                       ! Everything is at T-points here
+               DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
+                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                  zwgt = idx_bdy(ib_bdy)%nbd(ib,igrd)
+                  DO ik = 1, jpkm1
+                     zta = zwgt * ( trcdta_bdy(jn, ib_bdy)%trc(ib,ik) - trb(ii,ij,ik,jn) ) * tmask(ii,ij,ik)
+                     tra(ii,ij,ik,jn) = tra(ii,ij,ik,jn) + zta
+                  END DO
+               END DO
+            ENDIF
+         ENDDO
+      ENDDO
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('trc_bdy_dmp')
+      !
+   END SUBROUTINE trc_bdy_dmp
+ 
+#else
+   !!----------------------------------------------------------------------
+   !!   Dummy module                   NO Unstruct Open Boundary Conditions
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_bdy(kt)      ! Empty routine
+      WRITE(*,*) 'trc_bdy: You should not have seen this print! error?', kt
+   END SUBROUTINE trc_bdy
+
+   SUBROUTINE trc_bdy_dmp(kt)      ! Empty routine
+      WRITE(*,*) 'trc_bdy_dmp: You should not have seen this print! error?', kt
+   END SUBROUTINE trc_bdy_dmp
+
+#endif
+
+   !!======================================================================
+END MODULE trcbdy
