Index: branches/UKMO/dev_r5107_xios_initialize_toyoce/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90
===================================================================
--- branches/UKMO/dev_r5107_xios_initialize_toyoce/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90	(revision 5277)
+++ branches/UKMO/dev_r5107_xios_initialize_toyoce/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90	(revision 5580)
@@ -24,4 +24,5 @@
    USE p4zsed          !  Sedimentation
    USE p4zint          !  time interpolation
+   USE p4zrem          !  remineralisation
    USE iom             !  I/O manager
    USE trd_oce         !  Ocean trends variables
@@ -36,7 +37,14 @@
    PUBLIC   p4z_sms        ! called in p4zsms.F90
 
-   REAL(wp) :: alkbudget, no3budget, silbudget, ferbudget
-   INTEGER ::  numco2, numnut  !: logical unit for co2 budget
-
+   REAL(wp) :: alkbudget, no3budget, silbudget, ferbudget, po4budget
+   REAL(wp) :: xfact1, xfact2
+   INTEGER ::  numco2, numnut, numnit  !: logical unit for co2 budget
+
+   !!* Array used to indicate negative tracer values
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xnegtr     !: ???
+
+
+   !! * Substitutions
+#  include "top_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
@@ -61,20 +69,17 @@
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index      
       !!
-      INTEGER ::   jnt, jn, jl
+      INTEGER ::   ji, jj, jk, jnt, jn, jl
+      REAL(wp) ::  ztra
+#if defined key_kriest
+      REAL(wp) ::  zcoef1, zcoef2
+#endif
       CHARACTER (len=25) :: charout
-      REAL(wp), POINTER, DIMENSION(:,:,:,:)  :: ztrdpis
       !!---------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('p4z_sms')
       !
-      IF( l_trdtrc )  THEN
-         CALL wrk_alloc( jpi, jpj, jpk, jp_pisces, ztrdpis ) 
-         DO jn = 1, jp_pisces
-            jl = jn + jp_pcs0 - 1
-            ztrdpis(:,:,:,jn) = trn(:,:,:,jl)
-         ENDDO
-      ENDIF
-      !
       IF( kt == nittrc000 ) THEN
+        !
+        ALLOCATE( xnegtr(jpi,jpj,jpk) )
         !
         CALL p4z_che                              ! initialize the chemical constants
@@ -88,4 +93,26 @@
       IF( ln_pisdmp .AND. MOD( kt - nn_dttrc, nn_pisdmp ) == 0 )   CALL p4z_dmp( kt )      ! Relaxation of some tracers
       !
+      !                                                                    !   set time step size (Euler/Leapfrog)
+      IF( ( neuler == 0 .AND. kt == nittrc000 ) .OR. ln_top_euler ) THEN   ;    rfact = rdttrc(1)     !  at nittrc000
+      ELSEIF( kt <= nittrc000 + nn_dttrc )                          THEN   ;    rfact = 2. * rdttrc(1)   ! at nittrc000 or nittrc000+nn_dttrc (Leapfrog)
+      ENDIF
+      !
+      IF( ( ln_top_euler .AND. kt == nittrc000 )  .OR. ( .NOT.ln_top_euler .AND. kt <= nittrc000 + nn_dttrc ) ) THEN
+         rfactr  = 1. / rfact
+         rfact2  = rfact / FLOAT( nrdttrc )
+         rfact2r = 1. / rfact2
+         xstep = rfact2 / rday         ! Time step duration for biology
+         IF(lwp) WRITE(numout,*) 
+         IF(lwp) WRITE(numout,*) '    Passive Tracer  time step    rfact  = ', rfact, ' rdt = ', rdttra(1)
+         IF(lwp) write(numout,*) '    PISCES  Biology time step    rfact2 = ', rfact2
+         IF(lwp) WRITE(numout,*)
+      ENDIF
+
+      IF( ( neuler == 0 .AND. kt == nittrc000 ) .OR. ln_top_euler ) THEN
+         DO jn = jp_pcs0, jp_pcs1              !   SMS on tracer without Asselin time-filter
+            trb(:,:,:,jn) = trn(:,:,:,jn)
+         END DO
+      ENDIF
+      !
       IF( ndayflxtr /= nday_year ) THEN      ! New days
          !
@@ -105,28 +132,57 @@
       DO jnt = 1, nrdttrc          ! Potential time splitting if requested
          !
-         CALL p4z_bio (kt, jnt)    ! Biology
-         CALL p4z_sed (kt, jnt)    ! Sedimentation
-         !
+         CALL p4z_bio( kt, jnt )   ! Biology
+         CALL p4z_sed( kt, jnt )   ! Sedimentation
+         CALL p4z_lys( kt, jnt )   ! Compute CaCO3 saturation
+         CALL p4z_flx( kt, jnt )   ! Compute surface fluxes
+         !
+         xnegtr(:,:,:) = 1.e0
          DO jn = jp_pcs0, jp_pcs1
-            trb(:,:,:,jn) = trn(:,:,:,jn)
-         ENDDO
-         !
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     IF( ( trb(ji,jj,jk,jn) + tra(ji,jj,jk,jn) ) < 0.e0 ) THEN
+                        ztra             = ABS( trb(ji,jj,jk,jn) ) / ( ABS( tra(ji,jj,jk,jn) ) + rtrn )
+                        xnegtr(ji,jj,jk) = MIN( xnegtr(ji,jj,jk),  ztra )
+                     ENDIF
+                 END DO
+               END DO
+            END DO
+         END DO
+         !                                ! where at least 1 tracer concentration becomes negative
+         !                                ! 
+         DO jn = jp_pcs0, jp_pcs1
+           trb(:,:,:,jn) = trb(:,:,:,jn) + xnegtr(:,:,:) * tra(:,:,:,jn)
+         END DO
+        !
+         DO jn = jp_pcs0, jp_pcs1
+            tra(:,:,:,jn) = 0._wp
+         END DO
+         !
+         IF( ln_top_euler ) THEN
+            DO jn = jp_pcs0, jp_pcs1
+               trn(:,:,:,jn) = trb(:,:,:,jn)
+            END DO
+         ENDIF
       END DO
 
-      IF( l_trdtrc )  THEN
-         DO jn = 1, jp_pisces
-            jl = jn + jp_pcs0 - 1
-            ztrdpis(:,:,:,jn) = ( ztrdpis(:,:,:,jn) - trn(:,:,:,jl) ) * rfact2r
-         ENDDO
-      ENDIF
-      CALL p4z_lys( kt )             ! Compute CaCO3 saturation
-      CALL p4z_flx( kt )             ! Compute surface fluxes
-
-      DO jn = jp_pcs0, jp_pcs1
-        CALL lbc_lnk( trn(:,:,:,jn), 'T', 1. )
-        CALL lbc_lnk( trb(:,:,:,jn), 'T', 1. )
-        CALL lbc_lnk( tra(:,:,:,jn), 'T', 1. )
+#if defined key_kriest
+      ! 
+      zcoef1 = 1.e0 / xkr_massp 
+      zcoef2 = 1.e0 / xkr_massp / 1.1
+      DO jk = 1,jpkm1
+         trb(:,:,jk,jpnum) = MAX(  trb(:,:,jk,jpnum), trb(:,:,jk,jppoc) * zcoef1 / xnumm(jk)  )
+         trb(:,:,jk,jpnum) = MIN(  trb(:,:,jk,jpnum), trb(:,:,jk,jppoc) * zcoef2              )
       END DO
       !
+#endif
+      !
+      !
+      IF( l_trdtrc ) THEN
+         DO jn = jp_pcs0, jp_pcs1
+           CALL trd_trc( tra(:,:,:,jn), jn, jptra_sms, kt )   ! save trends
+         END DO
+      END IF
+      !
       IF( lk_sed ) THEN 
          !
@@ -134,5 +190,5 @@
          !
          DO jn = jp_pcs0, jp_pcs1
-           CALL lbc_lnk( trn(:,:,:,jn), 'T', 1. )
+           CALL lbc_lnk( trb(:,:,:,jn), 'T', 1. )
          END DO
          !
@@ -141,13 +197,5 @@
       IF( lrst_trc )  CALL p4z_rst( kt, 'WRITE' )  !* Write PISCES informations in restart file 
       !
-      IF( l_trdtrc ) THEN
-         DO jn = 1, jp_pisces
-            jl = jn + jp_pcs0 - 1
-             ztrdpis(:,:,:,jn) = ztrdpis(:,:,:,jn) + tra(:,:,:,jl)
-             CALL trd_trc( ztrdpis(:,:,:,jn), jn, jptra_sms, kt )   ! save trends
-          END DO
-          CALL wrk_dealloc( jpi, jpj, jpk, jp_pisces, ztrdpis ) 
-      END IF
-      !
+
       IF( lk_iomput .OR. ln_check_mass )  CALL p4z_chk_mass( kt ) ! Mass conservation checking
 
@@ -280,7 +328,7 @@
                ztmas   = tmask(ji,jj,jk)
                ztmas1  = 1. - tmask(ji,jj,jk)
-               zcaralk = trn(ji,jj,jk,jptal) - borat(ji,jj,jk) / (  1. + 1.E-8 / ( rtrn + akb3(ji,jj,jk) )  )
-               zco3    = ( zcaralk - trn(ji,jj,jk,jpdic) ) * ztmas + 0.5e-3 * ztmas1
-               zbicarb = ( 2. * trn(ji,jj,jk,jpdic) - zcaralk )
+               zcaralk = trb(ji,jj,jk,jptal) - borat(ji,jj,jk) / (  1. + 1.E-8 / ( rtrn + akb3(ji,jj,jk) )  )
+               zco3    = ( zcaralk - trb(ji,jj,jk,jpdic) ) * ztmas + 0.5e-3 * ztmas1
+               zbicarb = ( 2. * trb(ji,jj,jk,jpdic) - zcaralk )
                hi(ji,jj,jk) = ( ak23(ji,jj,jk) * zbicarb / zco3 ) * ztmas + 1.e-9 * ztmas1
             END DO
@@ -361,5 +409,6 @@
       REAL(wp) ::  silmean = 91.51     ! mean value of silicate
       !
-      REAL(wp) :: zarea, zalksum, zpo4sum, zno3sum, zsilsum
+      REAL(wp) :: zarea, zalksumn, zpo4sumn, zno3sumn, zsilsumn
+      REAL(wp) :: zalksumb, zpo4sumb, zno3sumb, zsilsumb
       !!---------------------------------------------------------------------
 
@@ -374,23 +423,44 @@
          zarea          = 1._wp / glob_sum( cvol(:,:,:) ) * 1e6              
 
-         zalksum = glob_sum( trn(:,:,:,jptal) * cvol(:,:,:)  ) * zarea
-         zpo4sum = glob_sum( trn(:,:,:,jppo4) * cvol(:,:,:)  ) * zarea * po4r
-         zno3sum = glob_sum( trn(:,:,:,jpno3) * cvol(:,:,:)  ) * zarea * rno3
-         zsilsum = glob_sum( trn(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea
+         zalksumn = glob_sum( trn(:,:,:,jptal) * cvol(:,:,:)  ) * zarea
+         zpo4sumn = glob_sum( trn(:,:,:,jppo4) * cvol(:,:,:)  ) * zarea * po4r
+         zno3sumn = glob_sum( trn(:,:,:,jpno3) * cvol(:,:,:)  ) * zarea * rno3
+         zsilsumn = glob_sum( trn(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea
  
-         IF(lwp) WRITE(numout,*) '       TALK mean : ', zalksum
-         trn(:,:,:,jptal) = trn(:,:,:,jptal) * alkmean / zalksum
-
-         IF(lwp) WRITE(numout,*) '       PO4  mean : ', zpo4sum
-         trn(:,:,:,jppo4) = trn(:,:,:,jppo4) * po4mean / zpo4sum
-
-         IF(lwp) WRITE(numout,*) '       NO3  mean : ', zno3sum
-         trn(:,:,:,jpno3) = trn(:,:,:,jpno3) * no3mean / zno3sum
-
-         IF(lwp) WRITE(numout,*) '       SiO3 mean : ', zsilsum
-         trn(:,:,:,jpsil) = MIN( 400.e-6,trn(:,:,:,jpsil) * silmean / zsilsum )
-         !
-      ENDIF
-
+         IF(lwp) WRITE(numout,*) '       TALKN mean : ', zalksumn
+         trn(:,:,:,jptal) = trn(:,:,:,jptal) * alkmean / zalksumn
+
+         IF(lwp) WRITE(numout,*) '       PO4N  mean : ', zpo4sumn
+         trn(:,:,:,jppo4) = trn(:,:,:,jppo4) * po4mean / zpo4sumn
+
+         IF(lwp) WRITE(numout,*) '       NO3N  mean : ', zno3sumn
+         trn(:,:,:,jpno3) = trn(:,:,:,jpno3) * no3mean / zno3sumn
+
+         IF(lwp) WRITE(numout,*) '       SiO3N mean : ', zsilsumn
+         trn(:,:,:,jpsil) = MIN( 400.e-6,trn(:,:,:,jpsil) * silmean / zsilsumn )
+         !
+         !
+         IF( .NOT. ln_top_euler ) THEN
+            zalksumb = glob_sum( trb(:,:,:,jptal) * cvol(:,:,:)  ) * zarea
+            zpo4sumb = glob_sum( trb(:,:,:,jppo4) * cvol(:,:,:)  ) * zarea * po4r
+            zno3sumb = glob_sum( trb(:,:,:,jpno3) * cvol(:,:,:)  ) * zarea * rno3
+            zsilsumb = glob_sum( trb(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea
+ 
+            IF(lwp) WRITE(numout,*) ' '
+            IF(lwp) WRITE(numout,*) '       TALKB mean : ', zalksumb
+            trb(:,:,:,jptal) = trb(:,:,:,jptal) * alkmean / zalksumb
+
+            IF(lwp) WRITE(numout,*) '       PO4B  mean : ', zpo4sumb
+            trb(:,:,:,jppo4) = trb(:,:,:,jppo4) * po4mean / zpo4sumb
+
+            IF(lwp) WRITE(numout,*) '       NO3B  mean : ', zno3sumb
+            trb(:,:,:,jpno3) = trb(:,:,:,jpno3) * no3mean / zno3sumb
+
+            IF(lwp) WRITE(numout,*) '       SiO3B mean : ', zsilsumb
+            trb(:,:,:,jpsil) = MIN( 400.e-6,trb(:,:,:,jpsil) * silmean / zsilsumb )
+        ENDIF
+        !
+      ENDIF
+        !
    END SUBROUTINE p4z_dmp
 
@@ -406,5 +476,11 @@
       INTEGER , INTENT( in ) ::   kt      ! ocean time-step index      
       REAL(wp)               ::  zfact       
-      !!
+      REAL(wp) ::  zrdenittot, zsdenittot, znitrpottot
+      CHARACTER(LEN=100)   ::   cltxt
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zvol
+      INTEGER :: jk
+      !!----------------------------------------------------------------------
+
+      !
       !!---------------------------------------------------------------------
 
@@ -413,7 +489,14 @@
             CALL ctl_opn( numco2, 'carbon.budget'  , 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
             CALL ctl_opn( numnut, 'nutrient.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
+            CALL ctl_opn( numnit, 'nitrogen.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
+            xfact1 = rfact2r * 12. / 1.e15 * ryyss    ! conversion molC/kt --> PgC/yr
+            xfact2 = 1.e+3 * rno3 * 14. / 1.e12 * ryyss   ! conversion molC/l/s ----> TgN/m3/yr
+            cltxt='time-step   Alkalinity        Nitrate        Phosphorus         Silicate           Iron'
+            IF( lwp ) WRITE(numnut,*)  TRIM(cltxt)
+            IF( lwp ) WRITE(numnut,*) 
          ENDIF
       ENDIF
 
+      !
       IF( iom_use( "pno3tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
          !   Compute the budget of NO3, ALK, Si, Fer
@@ -431,5 +514,18 @@
       ENDIF
       !
-      IF( iom_use( "psiltot" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
+      IF( iom_use( "ppo4tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
+         po4budget = glob_sum( (   trn(:,:,:,jppo4)                     &
+            &                    + trn(:,:,:,jpphy) + trn(:,:,:,jpdia)  &
+            &                    + trn(:,:,:,jpzoo) + trn(:,:,:,jpmes)  &
+            &                    + trn(:,:,:,jppoc)                     &
+#if ! defined key_kriest
+            &                    + trn(:,:,:,jpgoc)                     &
+#endif
+            &                    + trn(:,:,:,jpdoc)                     ) * cvol(:,:,:)  )
+         po4budget = po4budget / areatot
+         CALL iom_put( "ppo4tot", po4budget )
+      ENDIF
+      !
+      IF( iom_use( "psiltot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
          silbudget = glob_sum( (   trn(:,:,:,jpsil) + trn(:,:,:,jpgsi)  &
             &                    + trn(:,:,:,jpdsi)                     ) * cvol(:,:,:)  )
@@ -439,5 +535,5 @@
       ENDIF
       !
-      IF( iom_use( "palktot" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
+      IF( iom_use( "palktot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
          alkbudget = glob_sum( (   trn(:,:,:,jpno3) * rno3              &
             &                    + trn(:,:,:,jptal)                     &
@@ -448,5 +544,5 @@
       ENDIF
       !
-      IF( iom_use( "pfertot" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
+      IF( iom_use( "pfertot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
          ferbudget = glob_sum( (   trn(:,:,:,jpfer) + trn(:,:,:,jpnfe)  &
             &                    + trn(:,:,:,jpdfe)                     &
@@ -462,19 +558,45 @@
       ENDIF
       !
+
+      ! Global budget of N SMS : denitrification in the water column and in the sediment
+      !                          nitrogen fixation by the diazotrophs
+      ! --------------------------------------------------------------------------------
+      IF( iom_use( "tnfix" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
+         znitrpottot  = glob_sum ( nitrpot(:,:,:) * nitrfix * cvol(:,:,:) )
+         CALL iom_put( "tnfix"  , znitrpottot * 1.e+3 * rno3 )  ! Global  nitrogen fixation molC/l  to molN/m3 
+      ENDIF
+      !
+      IF( iom_use( "tdenit" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
+         zrdenittot   = glob_sum ( denitr(:,:,:) * rdenit * xnegtr(:,:,:) * cvol(:,:,:) )
+         CALL iom_put( "tdenit"  , zrdenittot * 1.e+3 * rno3 )  ! Total denitrification molC/l to molN/m3 
+      ENDIF
+      !
+      IF( iom_use( "Sdenit" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
+         zsdenittot   = glob_sum ( sdenit(:,:) * e1e2t(:,:) )
+         CALL iom_put( "Sdenit", sdenit(:,:) * zfact * tmask(:,:,1) )  ! Nitrate reduction in the sediments
+      ENDIF
+
       IF( ln_check_mass .AND. kt == nitend ) THEN   ! Compute the budget of NO3, ALK, Si, Fer
-         zfact = rfact2r * 12. / 1.e15 * ryyss    ! conversion molC/kt --> PgC/year
          t_atm_co2_flx  = t_atm_co2_flx / glob_sum( e1e2t(:,:) )
-         t_oce_co2_flx  = t_oce_co2_flx         * zfact * (-1 )
-         tpp            = tpp           * 1000. * zfact
-         t_oce_co2_exp  = t_oce_co2_exp * 1000. * zfact
+         t_oce_co2_flx  = t_oce_co2_flx         * xfact1 * (-1 )
+         tpp            = tpp           * 1000. * xfact1
+         t_oce_co2_exp  = t_oce_co2_exp * 1000. * xfact1
          IF( lwp ) WRITE(numco2,9000) ndastp, t_atm_co2_flx, t_oce_co2_flx, tpp, t_oce_co2_exp
-         IF( lwp ) WRITE(numnut,9500) ndastp, alkbudget        * 1.e+06, &
+         IF( lwp ) WRITE(numnut,9100) ndastp, alkbudget        * 1.e+06, &
              &                                no3budget * rno3 * 1.e+06, &
+             &                                po4budget * po4r * 1.e+06, &
              &                                silbudget        * 1.e+06, &
              &                                ferbudget        * 1.e+09
+         !
+         IF( lwp ) WRITE(numnit,9200) ndastp, znitrpottot * xfact2  , &
+         &                             zrdenittot  * xfact2  , &
+         &                             zsdenittot  * xfact2
+
       ENDIF
       !
  9000  FORMAT(i8,f10.5,e18.10,f10.5,f10.5)
- 9500  FORMAT(i8,4e18.10)
+ 9100  FORMAT(i8,5e18.10)
+ 9200  FORMAT(i8,3f10.5)
+
        !
    END SUBROUTINE p4z_chk_mass
