Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14_biophys/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14_biophys/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 13649)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14_biophys/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 13677)
@@ -54,5 +54,5 @@
                                    zn_sed_c, zn_sed_ca, zn_sed_fe,      &
                                    zn_sed_n, zn_sed_si, zn_chl_srf,     &
-                                   scl_chl, chl_out
+                                   scl_chl,  chl_out,   xpar
       USE trc,               ONLY: med_diag, nittrc000, trn 
       USE trcnam_trp,        ONLY: ln_trcadv_cen2, ln_trcadv_tvd
@@ -749,4 +749,7 @@
              DEALLOCATE( remin3dn )
           ENDIF
+          IF( med_diag%MED_XPAR3%dgsave ) THEN
+             CALL iom_put( "MED_XPAR3" , xpar )
+          ENDIF  
 # if defined key_roam          
           IF( med_diag%PH3%dgsave ) THEN
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14_biophys/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14_biophys/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90	(revision 13649)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14_biophys/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90	(revision 13677)
@@ -317,4 +317,5 @@
    REAL(wp) ::   xlg        !: exposant for pigment absorption in green (NAMELIST)
    REAL(wp) ::   rpig       !: chla/chla+phea ratio                     (NAMELIST)
+   LOGICAL  ::   ln_rgb     !: use RGB light scheme rather than 2-band  (NAMELIST)
                                                         
    INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   neln    !: number of levels in the euphotic layer
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14_biophys/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14_biophys/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90	(revision 13649)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14_biophys/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90	(revision 13677)
@@ -94,5 +94,5 @@
       &    xthetarem,xo2min 
 #endif
-      NAMELIST/natopt/xkg0,xkr0,xkgp,xkrp,xlg,xlr,rpig
+      NAMELIST/natopt/xkg0,xkr0,xkgp,xkrp,xlg,xlr,rpig,ln_rgb
       INTEGER :: jl, jn
       INTEGER :: ios                 ! Local integer output status for namelist read
@@ -1097,4 +1097,5 @@
       xlr   = 0.
       rpig  = 0.
+      ln_rgb = .false.
 
       !READ(numnatm,natopt)
@@ -1120,4 +1121,5 @@
          WRITE(numout,*) ' red   chl exposant              xlr   = ',xlr
          WRITE(numout,*) ' chla/chla+phea ratio            rpig  = ',rpig
+         WRITE(numout,*) ' use RGB scheme                ln_rgb  = ',ln_rgb
          WRITE(numout,*) ' '
 
@@ -1533,4 +1535,9 @@
           med_diag%MED_XPAR%dgsave = .FALSE.
       ENDIF
+      IF  (iom_use("MED_XPAR3")) THEN 
+          med_diag%MED_XPAR3%dgsave = .TRUE.
+      ELSE 
+          med_diag%MED_XPAR3%dgsave = .FALSE.
+      ENDIF
       IF  (iom_use("INTFLX_N")) THEN 
           med_diag%INTFLX_N%dgsave = .TRUE.
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14_biophys/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcopt_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14_biophys/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcopt_medusa.F90	(revision 13649)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14_biophys/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcopt_medusa.F90	(revision 13677)
@@ -28,4 +28,6 @@
 
    PUBLIC   trc_opt_medusa   ! called in trcprg.F90
+
+   REAL(wp), DIMENSION(3,61) :: okrgb   !: tabulated attenuation coefficients for RGB absorption
 
    !!* Substitution
@@ -50,5 +52,5 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt   ! index of the time stepping
-      INTEGER  ::   ji, jj, jk
+      INTEGER  ::   ji, jj, jk, irgb
       REAL(wp) ::   zpig                                    ! total pigment
       REAL(wp) ::   zkr                                     ! total absorption coefficient in red
@@ -58,4 +60,6 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::   zpar0m          ! irradiance just below the surface
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zparr, zparg    ! red and green compound of par
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zekb, zekg, zekr
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze1, ze2, ze3
 
       CHARACTER (len=25) :: charout
@@ -113,30 +117,65 @@
       ! ---------------------
 
-      DO jk = 2, jpk                     ! determination of local par in w levels
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               totchl =trn(ji,jj,jk-1,jpchn)+trn(ji,jj,jk-1,jpchd)
-               zpig = MAX( TINY(0.), totchl/rpig) 
-               zkr  = xkr0 + xkrp * EXP( xlr * LOG( zpig ) )
-               zkg  = xkg0 + xkgp * EXP( xlg * LOG( zpig ) )
-               zparr(ji,jj,jk) = zparr(ji,jj,jk-1) * EXP( -zkr * fse3t(ji,jj,jk-1) )
-               zparg(ji,jj,jk) = zparg(ji,jj,jk-1) * EXP( -zkg * fse3t(ji,jj,jk-1) )
+      IF ( ln_rgb ) THEN
+         ! Mean PAR in T levels using RGB scheme
+         ! Should be same as in traqsr, but T rather than W levels
+         IF( kt == nittrc000 ) THEN
+            CALL trc_oce_rgb( okrgb )
+         ENDIF
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  totchl = trn(ji,jj,jk,jpchn) + trn(ji,jj,jk,jpchd)
+                  totchl = MIN( 10. , MAX( 0.05, totchl ) )
+                  irgb = NINT( 41 + 20.* LOG10( totchl ) + 1.e-15 )
+                  !                                                         
+                  zekb(ji,jj,jk) = okrgb(1,irgb) * fse3t(ji,jj,jk)
+                  zekg(ji,jj,jk) = okrgb(2,irgb) * fse3t(ji,jj,jk)
+                  zekr(ji,jj,jk) = okrgb(3,irgb) * fse3t(ji,jj,jk)
+               END DO
             END DO
-        END DO
-      END DO
-
-      DO jk = 1, jpkm1                   ! mean par in t levels
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               totchl =trn(ji,jj,jk  ,jpchn)+trn(ji,jj,jk  ,jpchd)
-               zpig = MAX( TINY(0.), totchl/rpig) 
-               zkr  = xkr0 + xkrp * EXP( xlr * LOG( zpig ) )
-               zkg  = xkg0 + xkgp * EXP( xlg * LOG( zpig ) )
-               zparr(ji,jj,jk)    = zparr(ji,jj,jk) / zkr / fse3t(ji,jj,jk) * ( 1 - EXP( -zkr*fse3t(ji,jj,jk) ) )
-               zparg(ji,jj,jk)    = zparg(ji,jj,jk) / zkg / fse3t(ji,jj,jk) * ( 1 - EXP( -zkg*fse3t(ji,jj,jk) ) )
-               xpar (ji,jj,jk) = MAX( zparr(ji,jj,jk) + zparg(ji,jj,jk), 1.e-15 )
+         END DO
+         ze1(:,:,1) = zpar0m(:,:) * EXP( -0.5 * zekb(:,:,1) ) / 3.0
+         ze2(:,:,1) = zpar0m(:,:) * EXP( -0.5 * zekg(:,:,1) ) / 3.0
+         ze3(:,:,1) = zpar0m(:,:) * EXP( -0.5 * zekr(:,:,1) ) / 3.0
+         !
+         DO jk = 2, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ze1(ji,jj,jk) = ze1(ji,jj,jk-1) * EXP( -0.5 * ( zekb(ji,jj,jk-1) + zekb(ji,jj,jk) ) )
+                  ze2(ji,jj,jk) = ze2(ji,jj,jk-1) * EXP( -0.5 * ( zekg(ji,jj,jk-1) + zekg(ji,jj,jk) ) )
+                  ze3(ji,jj,jk) = ze3(ji,jj,jk-1) * EXP( -0.5 * ( zekr(ji,jj,jk-1) + zekr(ji,jj,jk) ) )
+                  xpar(ji,jj,jk) = MAX( ze1(ji,jj,jk) + ze2(ji,jj,jk) + ze3(ji,jj,jk), 1.e-15 )
+               END DO
             END DO
          END DO
-      END DO
+      ELSE
+         DO jk = 2, jpk                     ! determination of local par in w levels
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  totchl =trn(ji,jj,jk-1,jpchn)+trn(ji,jj,jk-1,jpchd)
+                  zpig = MAX( TINY(0.), totchl/rpig) 
+                  zkr  = xkr0 + xkrp * EXP( xlr * LOG( zpig ) )
+                  zkg  = xkg0 + xkgp * EXP( xlg * LOG( zpig ) )
+                  zparr(ji,jj,jk) = zparr(ji,jj,jk-1) * EXP( -zkr * fse3t(ji,jj,jk-1) )
+                  zparg(ji,jj,jk) = zparg(ji,jj,jk-1) * EXP( -zkg * fse3t(ji,jj,jk-1) )
+               END DO
+           END DO
+         END DO
+
+         DO jk = 1, jpkm1                   ! mean par in t levels
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  totchl =trn(ji,jj,jk  ,jpchn)+trn(ji,jj,jk  ,jpchd)
+                  zpig = MAX( TINY(0.), totchl/rpig) 
+                  zkr  = xkr0 + xkrp * EXP( xlr * LOG( zpig ) )
+                  zkg  = xkg0 + xkgp * EXP( xlg * LOG( zpig ) )
+                  zparr(ji,jj,jk)    = zparr(ji,jj,jk) / zkr / fse3t(ji,jj,jk) * ( 1 - EXP( -zkr*fse3t(ji,jj,jk) ) )
+                  zparg(ji,jj,jk)    = zparg(ji,jj,jk) / zkg / fse3t(ji,jj,jk) * ( 1 - EXP( -zkg*fse3t(ji,jj,jk) ) )
+                  xpar (ji,jj,jk) = MAX( zparr(ji,jj,jk) + zparg(ji,jj,jk), 1.e-15 )
+               END DO
+            END DO
+         END DO
+      ENDIF
 
       ! 3. Determination of euphotic layer depth
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14_biophys/NEMOGCM/NEMO/TOP_SRC/trc.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14_biophys/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 13649)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14_biophys/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 13677)
@@ -154,5 +154,7 @@
                   MIGRAZP3, MIGRAZD3, MEGRAZP3, MEGRAZD3, MEGRAZZ3,                                  &
                   O2SAT3, PBSI3, PCAL3, REMOC3,                                                      &
-                  PNLIMJ3, PNLIMN3, PNLIMFE3, PDLIMJ3, PDLIMN3, PDLIMFE3, PDLIMSI3       
+                  PNLIMJ3, PNLIMN3, PNLIMFE3, PDLIMJ3, PDLIMN3, PDLIMFE3, PDLIMSI3,                  &
+! DAF (22/10/20): some more diagnostics
+                  MED_XPAR3
                   !!
                   !! list of all MEDUSA diagnostics that could be called by iom_use
