Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 8442)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 8521)
@@ -6,4 +6,5 @@
    !! History :
    !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!   -   ! 2017-08 (A. Yool)            Amend bethic reservoir updating
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -62,159 +63,52 @@
       REAL(wp) :: fq0,fq1,fq2,fq3
 
+# if defined key_roam                     
+      !!----------------------------------------------------------------------
+      !! AXY (09/08/17): fix benthic submodel
       !!----------------------------------------------------------------------
       !! Process benthic in/out fluxes
       !! These can be handled outside of the 3D calculations since the
-      !! benthic pools (and fluxes) are 2D in nature; this code is
-      !! (shamelessly) borrowed from corresponding code in the LOBSTER
-      !! model
+      !! benthic pools (and fluxes) are 2D in nature; this code was
+      !! developed with help from George Nurser (NOC); it cannot be run
+      !! in a configuration with variable time-stepping with depth
       !!----------------------------------------------------------------------
       !!
-      !! IF(lwp) WRITE(numout,*) 'AXY: rdt = ', rdt
+      !! time-step calculation
       if (jorgben.eq.1) then
-         za_sed_n(:,:)  = zn_sed_n(:,:)  +                                 & 
-                          ( f_sbenin_n(:,:)  + f_fbenin_n(:,:)  -          &
-                            f_benout_n(:,:)  ) * (rdt / 86400.)
+         za_sed_n(:,:)  = zb_sed_n(:,:)  + ((2. * (rdt / 86400.)) * &
+                          ( f_sbenin_n(:,:)  + f_fbenin_n(:,:)  - f_benout_n(:,:)  ))
+         za_sed_fe(:,:) = zb_sed_fe(:,:) + ((2. * (rdt / 86400.)) * &
+                          ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) - f_benout_fe(:,:) ))
+         za_sed_c(:,:)  = zb_sed_c(:,:)  + ((2. * (rdt / 86400.)) * &
+                          ( f_sbenin_c(:,:)  + f_fbenin_c(:,:)  - f_benout_c(:,:)  ))
+      endif
+      if (jinorgben.eq.1) then
+         za_sed_si(:,:) = zb_sed_si(:,:) + ((2. * (rdt / 86400.)) * &
+                          ( f_fbenin_si(:,:) - f_benout_si(:,:) ))
+         za_sed_ca(:,:) = zb_sed_ca(:,:) + ((2. * (rdt / 86400.)) * &
+                          ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ))
+      endif
+      !!
+      !! time-level calculation
+      if (jorgben.eq.1) then
+         zb_sed_n(:,:)  = zn_sed_n(:,:)  + (atfp * &
+                          ( za_sed_n(:,:)  - (2. * zn_sed_n(:,:))  + zb_sed_n(:,:)  ))
          zn_sed_n(:,:)  = za_sed_n(:,:)
-         !!
-         za_sed_fe(:,:) = zn_sed_fe(:,:) +                                 &
-                          ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) -          &
-                            f_benout_fe(:,:) ) * (rdt / 86400.)
+         zb_sed_fe(:,:) = zn_sed_fe(:,:) + (atfp * &
+                          ( za_sed_fe(:,:) - (2. * zn_sed_fe(:,:)) + zb_sed_fe(:,:) ))
          zn_sed_fe(:,:) = za_sed_fe(:,:)
-         !!
-         za_sed_c(:,:)  = zn_sed_c(:,:)  +                                 &
-                          ( f_sbenin_c(:,:)  + f_fbenin_c(:,:)  -          &
-                            f_benout_c(:,:)  ) * (rdt / 86400.)
+         zb_sed_c(:,:)  = zn_sed_c(:,:)  + (atfp * &
+                          ( za_sed_c(:,:)  - (2. * zn_sed_c(:,:))  + zb_sed_c(:,:)  ))
          zn_sed_c(:,:)  = za_sed_c(:,:)
       endif
       if (jinorgben.eq.1) then
-         za_sed_si(:,:) = zn_sed_si(:,:) +                                 & 
-                          ( f_fbenin_si(:,:) - f_benout_si(:,:) ) *        &
-                          (rdt / 86400.)
+         zb_sed_si(:,:) = zn_sed_si(:,:) + (atfp * &
+                          ( za_sed_si(:,:) - (2. * zn_sed_si(:,:)) + zb_sed_si(:,:) ))
          zn_sed_si(:,:) = za_sed_si(:,:)
-         !!
-         za_sed_ca(:,:) = zn_sed_ca(:,:) +                                 &
-                          ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ) *        &
-                          (rdt / 86400.)
+         zb_sed_ca(:,:) = zn_sed_ca(:,:) + (atfp * &
+                          ( za_sed_ca(:,:) - (2. * zn_sed_ca(:,:)) + zb_sed_ca(:,:) ))
          zn_sed_ca(:,:) = za_sed_ca(:,:)
       endif
-      !!
-      if (ibenthic.eq.2) then
-         !! The code below (in this if ... then ... endif loop) is
-         !! effectively commented out because it does not work as 
-         !! anticipated; it can be deleted at a later date
-         if (jorgben.eq.1) then
-            za_sed_n(:,:)  = ( f_sbenin_n(:,:)  + f_fbenin_n(:,:)  -       &
-                               f_benout_n(:,:)  ) * rdt
-            za_sed_fe(:,:) = ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) -       &
-                               f_benout_fe(:,:) ) * rdt
-            za_sed_c(:,:)  = ( f_sbenin_c(:,:)  + f_fbenin_c(:,:)  -       &
-                               f_benout_c(:,:)  ) * rdt
-         endif
-         if (jinorgben.eq.1) then
-            za_sed_si(:,:) = ( f_fbenin_si(:,:) - f_benout_si(:,:) ) * rdt
-            za_sed_ca(:,:) = ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ) * rdt
-         endif
-         !!
-         !! Leap-frog scheme - only in explicit case, otherwise the 
-         !! time stepping is already being done in trczdf
-         !! IF( l_trczdf_exp .AND. (ln_trcadv_cen2 .OR. ln_trcadv_tvd) ) THEN
-         !!    zfact = 2. * rdttra(jk) * FLOAT( ndttrc )
-         !!    IF( neuler == 0 .AND. kt == nittrc000 )  zfact = rdttra(jk) * 
-         !!                                             FLOAT(ndttrc)
-         !!    if (jorgben.eq.1) then
-         !!       za_sed_n(:,:)  = zb_sed_n(:,:)  + ( zfact * za_sed_n(:,:)  )
-         !!      za_sed_fe(:,:) = zb_sed_fe(:,:) + ( zfact * za_sed_fe(:,:) )
-         !!       za_sed_c(:,:)  = zb_sed_c(:,:)  + ( zfact * za_sed_c(:,:)  )
-         !!    endif
-         !!    if (jinorgben.eq.1) then
-         !!       za_sed_si(:,:) = zb_sed_si(:,:) + ( zfact * za_sed_si(:,:) )
-         !!       za_sed_ca(:,:) = zb_sed_ca(:,:) + ( zfact * za_sed_ca(:,:) )
-         !!    endif
-         !! ENDIF
-         !! 
-         !! Time filter and swap of arrays
-         IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd  ) THEN ! centred or tvd scheme
-            IF( neuler == 0 .AND. kt == nittrc000 ) THEN
-               if (jorgben.eq.1) then
-                  zb_sed_n(:,:)  = zn_sed_n(:,:)
-                  zn_sed_n(:,:)  = za_sed_n(:,:)
-                  za_sed_n(:,:)  = 0.0
-                  !!
-                  zb_sed_fe(:,:) = zn_sed_fe(:,:)
-                  zn_sed_fe(:,:) = za_sed_fe(:,:)
-                  za_sed_fe(:,:) = 0.0
-                  !!
-                  zb_sed_c(:,:)  = zn_sed_c(:,:)
-                  zn_sed_c(:,:)  = za_sed_c(:,:)
-                  za_sed_c(:,:)  = 0.0
-               endif
-               if (jinorgben.eq.1) then
-                  zb_sed_si(:,:) = zn_sed_si(:,:)
-                  zn_sed_si(:,:) = za_sed_si(:,:)
-                  za_sed_si(:,:) = 0.0
-                  !!
-                  zb_sed_ca(:,:) = zn_sed_ca(:,:)
-                  zn_sed_ca(:,:) = za_sed_ca(:,:)
-                  za_sed_ca(:,:) = 0.0
-               endif
-            ELSE
-               if (jorgben.eq.1) then
-                  zb_sed_n(:,:)  = (atfp  *                                 &
-                                    ( zb_sed_n(:,:)  + za_sed_n(:,:)  )) +  &
-                                      (atfp1 * zn_sed_n(:,:) )
-                  zn_sed_n(:,:)  = za_sed_n(:,:)
-                  za_sed_n(:,:)  = 0.0
-                  !!
-                  zb_sed_fe(:,:) = (atfp  *                                 &
-                                    ( zb_sed_fe(:,:) + za_sed_fe(:,:) )) +  &
-                                      (atfp1 * zn_sed_fe(:,:))
-                  zn_sed_fe(:,:) = za_sed_fe(:,:)
-                  za_sed_fe(:,:) = 0.0
-                  !!
-                  zb_sed_c(:,:)  = (atfp  *                                 &
-                                    ( zb_sed_c(:,:)  + za_sed_c(:,:)  )) +  &
-                                      (atfp1 * zn_sed_c(:,:) )
-                  zn_sed_c(:,:)  = za_sed_c(:,:)
-                  za_sed_c(:,:)  = 0.0
-               endif
-               if (jinorgben.eq.1) then
-                  zb_sed_si(:,:) = (atfp  *                                 &
-                                    ( zb_sed_si(:,:) + za_sed_si(:,:) )) +  &
-                                      (atfp1 * zn_sed_si(:,:))
-                  zn_sed_si(:,:) = za_sed_si(:,:)
-                  za_sed_si(:,:) = 0.0
-                  !!
-                  zb_sed_ca(:,:) = (atfp  *                                 &
-                                    ( zb_sed_ca(:,:) + za_sed_ca(:,:) )) +  &
-                                      (atfp1 * zn_sed_ca(:,:))
-                  zn_sed_ca(:,:) = za_sed_ca(:,:)
-                  za_sed_ca(:,:) = 0.0
-               endif
-            ENDIF
-         ELSE                   !  case of smolar scheme or muscl
-            if (jorgben.eq.1) then
-               zb_sed_n(:,:)  = za_sed_n(:,:)
-               zn_sed_n(:,:)  = za_sed_n(:,:)
-               za_sed_n(:,:)  = 0.0
-               !!
-               zb_sed_fe(:,:) = za_sed_fe(:,:)
-               zn_sed_fe(:,:) = za_sed_fe(:,:)
-               za_sed_fe(:,:) = 0.0
-               !!
-               zb_sed_c(:,:)  = za_sed_c(:,:)
-               zn_sed_c(:,:)  = za_sed_c(:,:)
-               za_sed_c(:,:)  = 0.0
-            endif
-            if (jinorgben.eq.1) then
-               zb_sed_si(:,:) = za_sed_si(:,:)
-               zn_sed_si(:,:) = za_sed_si(:,:)
-               za_sed_si(:,:) = 0.0
-               !!
-               zb_sed_ca(:,:) = za_sed_ca(:,:)
-               zn_sed_ca(:,:) = za_sed_ca(:,:)
-               za_sed_ca(:,:) = 0.0
-            endif
-         ENDIF
-      endif
+# endif      
 
 #  if defined key_debug_medusa
@@ -253,7 +147,8 @@
                   fq1 = f_sbenin_n(ji,jj) + f_fbenin_n(ji,jj)
                   fq2 = fq0 + fq1
-                  IF (lwp) write (numout,'(a,2i3,a,3f15.10)')               &
-                     'AXY N   cons: (i,j)=',ji,jj,', (flx,ben,err)=',      &
-                     fq0,fq1,fq2
+                  fq3 = f_benout_n(ji,jj)
+                  if (lwp) write (numout,'a,2i3,a,4f15,5)')                   &
+                     'AXY N   cons: (i,j)=',ji,jj,', (flx,ben,err,out)=',      &
+                     fq0,fq1,fq2,fq3
                ENDIF
             ENDDO
@@ -266,7 +161,8 @@
                   fq1 = f_fbenin_si(ji,jj)
                   fq2 = fq0 + fq1
-                  IF (lwp) write (numout,'(a,2i3,a,3f15.10)')               &
-                     'AXY Si  cons: (i,j)=',ji,jj,', (flx,ben,err)=',     &
-                     fq0,fq1,fq2
+                  fq3 = f_benout_si(ji,jj)
+                  if (lwp) write (numout,'a,2i3,a,4f15,5)')                   &
+                     'AXY Si  cons: (i,j)=',ji,jj,', (flx,ben,err,out)=',     &
+                     fq0,fq1,fq2,fq3
                ENDIF
             ENDDO
@@ -278,22 +174,24 @@
                   fq0 = fflx_c(ji,jj)
                   fq1 = f_sbenin_c(ji,jj) + f_fbenin_c(ji,jj) + f_fbenin_ca(ji,jj)
-                  fq2 = f_co2flux(ji,jj) * e3t_n(ji,jj,1)
+                  fq2 = f_co2flux(ji,jj) * fse3t(ji,jj,1)
                   fq3 = fq0 + fq1
-                  IF (lwp) write (numout,'(a,2i3,a,4f15.10)')               &
-                    'AXY C   cons: (i,j)=',ji,jj,', (flx,ben,asf,err)=',  &
-                    fq0,fq1,fq2,fq3
-               ENDIF
-            ENDDO
-         ENDDO   
-         !! alkalinity
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               if (tmask(ji,jj,1) == 1) then
-                  fq0 = fflx_a(ji,jj)
-                  fq1 = 2.0 * f_fbenin_ca(ji,jj)
-                  fq2 = fq0 + fq1
-                  IF (lwp) write (numout,'(a,2i3,a,3f15.10)')               &
-                     'AXY alk cons: (i,j)=',ji,jj,', (flx,ben,err)=',     &
-                     fq0,fq1,fq2
+                  fq4 = f_benout_c(ji,jj) + f_benout_ca(ji,jj)
+                  if (lwp) write (numout,'a,2i3,a,5f15,5)')                   &
+                     'AXY C   cons: (i,j)=',ji,jj,', (flx,ben,asf,err,out)=', &
+                     fq0,fq1,fq2,fq3,fq4
+                ENDIF
+             ENDDO
+          ENDDO   
+          !! alkalinity
+          DO jj = 2,jpjm1
+             DO ji = 2,jpim1
+                if (tmask(ji,jj,1) == 1) then
+                   fq0 = fflx_a(ji,jj)
+                   fq1 = 2.0 * f_fbenin_ca(ji,jj)
+                   fq2 = fq0 + fq1
+                   fq3 = 2.0 * f_benout_ca(ji,jj)
+                   if (lwp) write (numout,'a,2i3,a,4f15,5)')                   &
+                      'AXY alk cons: (i,j)=',ji,jj,', (flx,ben,err,out)=',     &
+                      fq0,fq1,fq2,fq3
                ENDIF
             ENDDO
@@ -360,4 +258,7 @@
             CALL iom_put( "OCAL_LVL"  , fccd )
          ENDIF
+         IF ( med_diag%CHL_MLD%dgsave ) THEN
+            CALL iom_put( "CHL_MLD"  , fchl_ml )
+         ENDIF
          IF ( med_diag%PN_JLIM%dgsave ) THEN
             CALL iom_put( "PN_JLIM"  , fjln2d )
