Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils216/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbcssr.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils216/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbcssr.F90	(revision 10916)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils216/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbcssr.F90	(revision 10917)
@@ -5,4 +5,5 @@
    !!======================================================================
    !! History :  3.6  !  2017-06  (D. Ford)  Adapt from sbcssr.F90
+   !!                 !  2019-04  (D. Ford)  Add biomass balancing and error checking
    !!----------------------------------------------------------------------
 #if defined key_top
@@ -50,4 +51,5 @@
       !!                   surface only    (nn_chldmp = 1)
       !!                   mixed layer     (nn_chldmp = 2)
+      !!              - balance biomass variables
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in   ) ::   kt   ! ocean time step
@@ -55,10 +57,20 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
  
-      REAL(wp), DIMENSION(jpi,jpj) :: ztra, zchl
-      REAL(wp)                     :: zpft
+      REAL(wp), DIMENSION(jpi,jpj)           :: zdiff, zchl
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra) :: ztra_save, ztra_new
+      REAL(wp)                               :: zpft, zrat, zinc
+      REAL(wp)                               :: zmax_chl, zmax_c, zmax_n, zmax_p, zmax_s
+      INTEGER                                :: inum
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('trc_sbc_ssr')
       !
+      ! Maximum permitted (absolute) updates to phytoplankton variables
+      zmax_chl = 0.01
+      zmax_c   = 0.1
+      zmax_n   = 0.01
+      zmax_p   = 0.001
+      zmax_s   = 0.01
+      
       IF( kt == nittrc000 )  THEN
          !
@@ -68,4 +80,6 @@
             !
             IF (lwp) WRITE(numout,*) 'Damping chlorophyll on timestep ', kt
+            
+            ztra_save(:,:,:,:) = tra(:,:,:,:)
             !
             CALL fld_read( kt, 1, sf_chldmp )   ! Read Chl data and provides it at kt
@@ -76,5 +90,5 @@
                &        trb(:,:,1,jp_fabm_m1+jp_fabm_chl3) + &
                &        trb(:,:,1,jp_fabm_m1+jp_fabm_chl4)
-            ztra(:,:) = rn_chldmp * ( sf_chldmp(1)%fnow(:,:,1) - zchl(:,:) )
+            zdiff(:,:) = rn_chldmp * ( sf_chldmp(1)%fnow(:,:,1) - zchl(:,:) )
             !
             DO jj = 2, jpjm1
@@ -83,46 +97,237 @@
                      & ( sf_chldmp(1)%fnow(ji,jj,1) < 100.0 ) .AND. &
                      & ( zchl(ji,jj)                >   0.0 ) ) THEN
-                     zpft = ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) / zchl(ji,jj) ) * ztra(ji,jj)
-                     tra(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) = tra(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + zpft
-                     IF( nn_chldmp == 2 ) THEN
-                        DO jk = 2, jpkm1
-                           IF( fsdept(ji,jj,jk) < hmlp (ji,jj) ) THEN
-                              tra(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) = tra(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) + zpft
-                           ENDIF
-                        END DO
+                     ! P1
+                     IF ( ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) > 0.0 ) .AND. &
+                        & ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_p1c)  > 0.0 ) .AND. &
+                        & ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_p1n)  > 0.0 ) .AND. &
+                        & ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_p1p)  > 0.0 ) .AND. &
+                        & ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_p1s)  > 0.0 ) ) THEN
+                        ! Chlorophyll
+                        zpft = ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) / zchl(ji,jj) ) * zdiff(ji,jj)
+                        zpft = MAX( -1.0 * zmax_chl, MIN( zpft, zmax_chl ) )
+                        tra(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) = tra(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + zpft
+                        ! Carbon
+                        zrat = trb(ji,jj,1,jp_fabm_m1+jp_fabm_p1c) / trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl1)
+                        zinc = MAX( -1.0 * zmax_c, MIN( (zpft * zrat), zmax_c ) )
+                        tra(ji,jj,1,jp_fabm_m1+jp_fabm_p1c) = tra(ji,jj,1,jp_fabm_m1+jp_fabm_p1c) + zinc
+                        ! Nitrogen
+                        zrat = trb(ji,jj,1,jp_fabm_m1+jp_fabm_p1n) / trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl1)
+                        zinc = MAX( -1.0 * zmax_n, MIN( (zpft * zrat), zmax_n ) )
+                        tra(ji,jj,1,jp_fabm_m1+jp_fabm_p1n) = tra(ji,jj,1,jp_fabm_m1+jp_fabm_p1n) + zinc
+                        ! Phosphorus
+                        zrat = trb(ji,jj,1,jp_fabm_m1+jp_fabm_p1p) / trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl1)
+                        zinc = MAX( -1.0 * zmax_p, MIN( (zpft * zrat), zmax_p ) )
+                        tra(ji,jj,1,jp_fabm_m1+jp_fabm_p1p) = tra(ji,jj,1,jp_fabm_m1+jp_fabm_p1p) + zinc
+                        ! Silicon
+                        zrat = trb(ji,jj,1,jp_fabm_m1+jp_fabm_p1s) / trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl1)
+                        zinc = MAX( -1.0 * zmax_s, MIN( (zpft * zrat), zmax_s ) )
+                        tra(ji,jj,1,jp_fabm_m1+jp_fabm_p1s) = tra(ji,jj,1,jp_fabm_m1+jp_fabm_p1s) + zinc
+                        !
+                        IF( nn_chldmp == 2 ) THEN
+                           ! Mixed layer
+                           DO jk = 2, jpkm1
+                              IF( fsdept(ji,jj,jk) < hmlp (ji,jj) ) THEN
+                                 IF ( ( trb(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) > 0.0 ) .AND. &
+                                    & ( trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p1c)  > 0.0 ) .AND. &
+                                    & ( trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n)  > 0.0 ) .AND. &
+                                    & ( trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p1p)  > 0.0 ) .AND. &
+                                    & ( trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p1s)  > 0.0 ) ) THEN
+                                    ! Chlorophyll
+                                    tra(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) = tra(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) + zpft
+                                    ! Carbon
+                                    zrat = trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p1c) / trb(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1)
+                                    zinc = MAX( -1.0 * zmax_c, MIN( (zpft * zrat), zmax_c ) )
+                                    tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p1c) = tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p1c) + zinc
+                                    ! Nitrogen
+                                    zrat = trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) / trb(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1)
+                                    zinc = MAX( -1.0 * zmax_n, MIN( (zpft * zrat), zmax_n ) )
+                                    tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) = tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + zinc
+                                    ! Phosphorus
+                                    zrat = trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p1p) / trb(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1)
+                                    zinc = MAX( -1.0 * zmax_p, MIN( (zpft * zrat), zmax_p ) )
+                                    tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p1p) = tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p1p) + zinc
+                                    ! Silicon
+                                    zrat = trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p1s) / trb(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1)
+                                    zinc = MAX( -1.0 * zmax_s, MIN( (zpft * zrat), zmax_s ) )
+                                    tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p1s) = tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p1s) + zinc
+                                 ENDIF
+                              ENDIF
+                           END DO
+                        ENDIF
                      ENDIF
-                     !
-                     zpft = ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) / zchl(ji,jj) ) * ztra(ji,jj)
-                     tra(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) = tra(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + zpft
-                     IF( nn_chldmp == 2 ) THEN
-                        DO jk = 2, jpkm1
-                           IF( fsdept(ji,jj,jk) < hmlp (ji,jj) ) THEN
-                              tra(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) = tra(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) + zpft
-                           ENDIF
-                        END DO
+                     ! P2
+                     IF ( ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) > 0.0 ) .AND. &
+                        & ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_p2c)  > 0.0 ) .AND. &
+                        & ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_p2n)  > 0.0 ) .AND. &
+                        & ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_p2p)  > 0.0 ) ) THEN
+                        ! Chlorophyll
+                        zpft = ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) / zchl(ji,jj) ) * zdiff(ji,jj)
+                        zpft = MAX( -1.0 * zmax_chl, MIN( zpft, zmax_chl ) )
+                        tra(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) = tra(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + zpft
+                        ! Carbon
+                        zrat = trb(ji,jj,1,jp_fabm_m1+jp_fabm_p2c) / trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl2)
+                        zinc = MAX( -1.0 * zmax_c, MIN( (zpft * zrat), zmax_c ) )
+                        tra(ji,jj,1,jp_fabm_m1+jp_fabm_p2c) = tra(ji,jj,1,jp_fabm_m1+jp_fabm_p2c) + zinc
+                        ! Nitrogen
+                        zrat = trb(ji,jj,1,jp_fabm_m1+jp_fabm_p2n) / trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl2)
+                        zinc = MAX( -1.0 * zmax_n, MIN( (zpft * zrat), zmax_n ) )
+                        tra(ji,jj,1,jp_fabm_m1+jp_fabm_p2n) = tra(ji,jj,1,jp_fabm_m1+jp_fabm_p2n) + zinc
+                        ! Phosphorus
+                        zrat = trb(ji,jj,1,jp_fabm_m1+jp_fabm_p2p) / trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl2)
+                        zinc = MAX( -1.0 * zmax_p, MIN( (zpft * zrat), zmax_p ) )
+                        tra(ji,jj,1,jp_fabm_m1+jp_fabm_p2p) = tra(ji,jj,1,jp_fabm_m1+jp_fabm_p2p) + zinc
+                        !
+                        IF( nn_chldmp == 2 ) THEN
+                           ! Mixed layer
+                           DO jk = 2, jpkm1
+                              IF( fsdept(ji,jj,jk) < hmlp (ji,jj) ) THEN
+                                 IF ( ( trb(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) > 0.0 ) .AND. &
+                                    & ( trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p2c)  > 0.0 ) .AND. &
+                                    & ( trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n)  > 0.0 ) .AND. &
+                                    & ( trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p2p)  > 0.0 ) ) THEN
+                                    ! Chlorophyll
+                                    tra(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) = tra(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) + zpft
+                                    ! Carbon
+                                    zrat = trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p2c) / trb(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2)
+                                    zinc = MAX( -1.0 * zmax_c, MIN( (zpft * zrat), zmax_c ) )
+                                    tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p2c) = tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p2c) + zinc
+                                    ! Nitrogen
+                                    zrat = trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) / trb(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2)
+                                    zinc = MAX( -1.0 * zmax_n, MIN( (zpft * zrat), zmax_n ) )
+                                    tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) = tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + zinc
+                                    ! Phosphorus
+                                    zrat = trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p2p) / trb(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2)
+                                    zinc = MAX( -1.0 * zmax_p, MIN( (zpft * zrat), zmax_p ) )
+                                    tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p2p) = tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p2p) + zinc
+                                 ENDIF
+                              ENDIF
+                           END DO
+                        ENDIF
                      ENDIF
-                     !
-                     zpft = ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) / zchl(ji,jj) ) * ztra(ji,jj)
-                     tra(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) = tra(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + zpft
-                     IF( nn_chldmp == 2 ) THEN
-                        DO jk = 2, jpkm1
-                           IF( fsdept(ji,jj,jk) < hmlp (ji,jj) ) THEN
-                              tra(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) = tra(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) + zpft
-                           ENDIF
-                        END DO
+                     ! P3
+                     IF ( ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) > 0.0 ) .AND. &
+                        & ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_p3c)  > 0.0 ) .AND. &
+                        & ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_p3n)  > 0.0 ) .AND. &
+                        & ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_p3p)  > 0.0 ) ) THEN
+                        ! Chlorophyll
+                        zpft = ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) / zchl(ji,jj) ) * zdiff(ji,jj)
+                        zpft = MAX( -1.0 * zmax_chl, MIN( zpft, zmax_chl ) )
+                        tra(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) = tra(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + zpft
+                        ! Carbon
+                        zrat = trb(ji,jj,1,jp_fabm_m1+jp_fabm_p3c) / trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl3)
+                        zinc = MAX( -1.0 * zmax_c, MIN( (zpft * zrat), zmax_c ) )
+                        tra(ji,jj,1,jp_fabm_m1+jp_fabm_p3c) = tra(ji,jj,1,jp_fabm_m1+jp_fabm_p3c) + zinc
+                        ! Nitrogen
+                        zrat = trb(ji,jj,1,jp_fabm_m1+jp_fabm_p3n) / trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl3)
+                        zinc = MAX( -1.0 * zmax_n, MIN( (zpft * zrat), zmax_n ) )
+                        tra(ji,jj,1,jp_fabm_m1+jp_fabm_p3n) = tra(ji,jj,1,jp_fabm_m1+jp_fabm_p3n) + zinc
+                        ! Phosphorus
+                        zrat = trb(ji,jj,1,jp_fabm_m1+jp_fabm_p3p) / trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl3)
+                        zinc = MAX( -1.0 * zmax_p, MIN( (zpft * zrat), zmax_p ) )
+                        tra(ji,jj,1,jp_fabm_m1+jp_fabm_p3p) = tra(ji,jj,1,jp_fabm_m1+jp_fabm_p3p) + zinc
+                        !
+                        IF( nn_chldmp == 2 ) THEN
+                           ! Mixed layer
+                           DO jk = 2, jpkm1
+                              IF( fsdept(ji,jj,jk) < hmlp (ji,jj) ) THEN
+                                 IF ( ( trb(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) > 0.0 ) .AND. &
+                                    & ( trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p3c)  > 0.0 ) .AND. &
+                                    & ( trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n)  > 0.0 ) .AND. &
+                                    & ( trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p3p)  > 0.0 ) ) THEN
+                                    ! Chlorophyll
+                                    tra(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) = tra(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) + zpft
+                                    ! Carbon
+                                    zrat = trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p3c) / trb(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3)
+                                    zinc = MAX( -1.0 * zmax_c, MIN( (zpft * zrat), zmax_c ) )
+                                    tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p3c) = tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p3c) + zinc
+                                    ! Nitrogen
+                                    zrat = trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) / trb(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3)
+                                    zinc = MAX( -1.0 * zmax_n, MIN( (zpft * zrat), zmax_n ) )
+                                    tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) = tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + zinc
+                                    ! Phosphorus
+                                    zrat = trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p3p) / trb(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3)
+                                    zinc = MAX( -1.0 * zmax_p, MIN( (zpft * zrat), zmax_p ) )
+                                    tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p3p) = tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p3p) + zinc
+                                 ENDIF
+                              ENDIF
+                           END DO
+                        ENDIF
                      ENDIF
-                     !
-                     zpft = ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) / zchl(ji,jj) ) * ztra(ji,jj)
-                     tra(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) = tra(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) + zpft
-                     IF( nn_chldmp == 2 ) THEN
-                        DO jk = 2, jpkm1
-                           IF( fsdept(ji,jj,jk) < hmlp (ji,jj) ) THEN
-                              tra(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) = tra(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) + zpft
-                           ENDIF
-                        END DO
+                     ! P4
+                     IF ( ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) > 0.0 ) .AND. &
+                        & ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_p4c)  > 0.0 ) .AND. &
+                        & ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_p4n)  > 0.0 ) .AND. &
+                        & ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_p4p)  > 0.0 ) ) THEN
+                        ! Chlorophyll
+                        zpft = ( trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) / zchl(ji,jj) ) * zdiff(ji,jj)
+                        zpft = MAX( -1.0 * zmax_chl, MIN( zpft, zmax_chl ) )
+                        tra(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) = tra(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) + zpft
+                        ! Carbon
+                        zrat = trb(ji,jj,1,jp_fabm_m1+jp_fabm_p4c) / trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
+                        zinc = MAX( -1.0 * zmax_c, MIN( (zpft * zrat), zmax_c ) )
+                        tra(ji,jj,1,jp_fabm_m1+jp_fabm_p4c) = tra(ji,jj,1,jp_fabm_m1+jp_fabm_p4c) + zinc
+                        ! Nitrogen
+                        zrat = trb(ji,jj,1,jp_fabm_m1+jp_fabm_p4n) / trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
+                        zinc = MAX( -1.0 * zmax_n, MIN( (zpft * zrat), zmax_n ) )
+                        tra(ji,jj,1,jp_fabm_m1+jp_fabm_p4n) = tra(ji,jj,1,jp_fabm_m1+jp_fabm_p4n) + zinc
+                        ! Phosphorus
+                        zrat = trb(ji,jj,1,jp_fabm_m1+jp_fabm_p4p) / trb(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
+                        zinc = MAX( -1.0 * zmax_p, MIN( (zpft * zrat), zmax_p ) )
+                        tra(ji,jj,1,jp_fabm_m1+jp_fabm_p4p) = tra(ji,jj,1,jp_fabm_m1+jp_fabm_p4p) + zinc
+                        !
+                        IF( nn_chldmp == 2 ) THEN
+                           ! Mixed layer
+                           DO jk = 2, jpkm1
+                              IF( fsdept(ji,jj,jk) < hmlp (ji,jj) ) THEN
+                                 IF ( ( trb(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) > 0.0 ) .AND. &
+                                    & ( trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p4c)  > 0.0 ) .AND. &
+                                    & ( trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n)  > 0.0 ) .AND. &
+                                    & ( trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p4p)  > 0.0 ) ) THEN
+                                    ! Chlorophyll
+                                    tra(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) = tra(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) + zpft
+                                    ! Carbon
+                                    zrat = trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p4c) / trb(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4)
+                                    zinc = MAX( -1.0 * zmax_c, MIN( (zpft * zrat), zmax_c ) )
+                                    tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p4c) = tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p4c) + zinc
+                                    ! Nitrogen
+                                    zrat = trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) / trb(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4)
+                                    zinc = MAX( -1.0 * zmax_n, MIN( (zpft * zrat), zmax_n ) )
+                                    tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) = tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) + zinc
+                                    ! Phosphorus
+                                    zrat = trb(ji,jj,jk,jp_fabm_m1+jp_fabm_p4p) / trb(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4)
+                                    zinc = MAX( -1.0 * zmax_p, MIN( (zpft * zrat), zmax_p ) )
+                                    tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p4p) = tra(ji,jj,jk,jp_fabm_m1+jp_fabm_p4p) + zinc
+                                 ENDIF
+                              ENDIF
+                           END DO
+                        ENDIF
                      ENDIF
                   ENDIF
                END DO
             END DO
+            
+            ! Save updates
+            ztra_new(:,:,:,:) = tra(:,:,:,:) - ztra_save(:,:,:,:)
+            CALL iom_open( 'chldmp.balincs', inum, ldwrt = .TRUE., kiolib = jprstlib)
+            zdiff(:,:) = zdiff(:,:) * tmask(:,:,1)
+            CALL iom_rstput( kt, kt, inum, 'ChlTot', zdiff(:,:) )
+            CALL iom_rstput( kt, kt, inum, 'P1_Chl', ztra_new(:,:,:,jp_fabm_chl1) )
+            CALL iom_rstput( kt, kt, inum, 'P1_c',   ztra_new(:,:,:,jp_fabm_p1c)  )
+            CALL iom_rstput( kt, kt, inum, 'P1_n',   ztra_new(:,:,:,jp_fabm_p1n)  )
+            CALL iom_rstput( kt, kt, inum, 'P1_p',   ztra_new(:,:,:,jp_fabm_p1p)  )
+            CALL iom_rstput( kt, kt, inum, 'P1_s',   ztra_new(:,:,:,jp_fabm_p1s)  )
+            CALL iom_rstput( kt, kt, inum, 'P2_Chl', ztra_new(:,:,:,jp_fabm_chl2) )
+            CALL iom_rstput( kt, kt, inum, 'P2_c',   ztra_new(:,:,:,jp_fabm_p2c)  )
+            CALL iom_rstput( kt, kt, inum, 'P2_n',   ztra_new(:,:,:,jp_fabm_p2n)  )
+            CALL iom_rstput( kt, kt, inum, 'P2_p',   ztra_new(:,:,:,jp_fabm_p2p)  )
+            CALL iom_rstput( kt, kt, inum, 'P3_Chl', ztra_new(:,:,:,jp_fabm_chl3) )
+            CALL iom_rstput( kt, kt, inum, 'P3_c',   ztra_new(:,:,:,jp_fabm_p3c)  )
+            CALL iom_rstput( kt, kt, inum, 'P3_n',   ztra_new(:,:,:,jp_fabm_p3n)  )
+            CALL iom_rstput( kt, kt, inum, 'P3_p',   ztra_new(:,:,:,jp_fabm_p3p)  )
+            CALL iom_rstput( kt, kt, inum, 'P4_Chl', ztra_new(:,:,:,jp_fabm_chl4) )
+            CALL iom_rstput( kt, kt, inum, 'P4_c',   ztra_new(:,:,:,jp_fabm_p4c)  )
+            CALL iom_rstput( kt, kt, inum, 'P4_n',   ztra_new(:,:,:,jp_fabm_p4n)  )
+            CALL iom_rstput( kt, kt, inum, 'P4_p',   ztra_new(:,:,:,jp_fabm_p4p)  )
+            CALL iom_close( inum )
 #else
             CALL ctl_stop( 'STOP', 'trc_sbc_ssr: only works with FABM-ERSEM' )
