Index: branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 10248)
+++ branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 10249)
@@ -1241,5 +1241,6 @@
       IF(lwm) WRITE( numond, nameos )
       !
-      rau0        = 1026._wp                 !: volumic mass of reference     [kg/m3]
+      rau0        = 1020._wp                 !: volumic mass of reference     [kg/m3]
+!     rau0        = 1026._wp                 !: volumic mass of reference     [kg/m3]
       rcp         = 3991.86795711963_wp      !: heat capacity     [J/K]
       !
Index: branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
===================================================================
--- branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 10248)
+++ branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 10249)
@@ -100,4 +100,8 @@
          IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
       ENDIF
+! slwa unless you use l_trdtra too, the above switches off trend calculations for l_trdtrc
+         l_trd = .FALSE.
+         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
+!slwa
       !
       IF( l_trd )  THEN
Index: branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/tradwl.F90
===================================================================
--- branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/tradwl.F90	(revision 10249)
+++ branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/tradwl.F90	(revision 10249)
@@ -0,0 +1,238 @@
+MODULE tradwl
+   !!======================================================================
+   !!                       ***  MODULE  tradwl  ***
+   !! Ocean physics: solar radiation penetration in the top ocean levels
+   !!======================================================================
+   !! History :  POLCOMS  !  1996-10  (J. Holt)  Original code
+   !!   NEMO     3.2  !  2010-03  (E. O'Dea)  Import to Nemo for use in Shelf Model
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   tra_dwl      : trend due to the solar radiation penetration
+   !!   tra_dwl_init : solar radiation penetration initialization
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and active tracers
+   USE dom_oce         ! ocean space and time domain
+   USE sbc_oce         ! surface boundary condition: ocean
+   USE trc_oce         ! share SMS/Ocean variables
+   USE trd_oce         ! trends: ocean variables
+   USE trdtra          ! ocean active tracers trends 
+   USE in_out_manager  ! I/O manager
+   USE phycst          ! physical constants
+   USE prtctl          ! Print control
+   USE iom             ! I/O manager
+   USE fldread         ! read input fields
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   tra_dwl        ! routine called by step.F90 (ln_tradwl=T)
+
+   !                                           !!* Namelist namtra_qsr: penetrative solar radiation
+   LOGICAL , PUBLIC ::   ln_tradwl  = .TRUE.    ! light absorption (dwl) flag
+   LOGICAL , PUBLIC ::   ln_vary_lambda  = .TRUE.    ! vary Lambda or not flag
+   
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
+   !! $Id$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE tra_dwl( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_qsr  ***
+      !!
+      !! ** Purpose :   Compute the temperature trend due to the solar radiation
+      !!      penetration and add it to the general temperature trend.
+      !!
+      !! ** Method  : The profile of the solar radiation within the ocean is defined
+      !!
+      !!          Jason Holt Oct 96
+      !!
+      !!          Calculates change in temperature due to penetrating
+      !!          radiation, with cooling at the surface layer
+      !!
+      !!          rad=rad0*exp(lambda*z)
+      !!
+      !!       Heat input into box is between z=K and z=K+1 is RAD(K)-RAD(K+1)
+      !!
+      !!
+      !! ** Action  : - update ta with the penetrative solar radiation trend
+      !!              - save the trend in ttrd ('key_trdtra')
+      !!
+      !!----------------------------------------------------------------------
+      USE oce, ONLY :   ztrdt => ua   ! use ua as 3D workspace   
+      USE oce, ONLY :   ztrds => va   ! use va as 3D workspace   
+      !!
+      INTEGER, INTENT(in) ::   kt     ! ocean time-step
+      !!
+      INTEGER  ::   ji, jj, jk           ! dummy loop indices
+      INTEGER  ::   irgb                 ! temporary integers
+      REAL(wp) ::   zchl, zcoef, zsi0r   ! temporary scalars
+      REAL(wp) ::   zc0, zc1, zc2, zc3   !    -         -
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zekb, zekg, zekr            ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze0, ze1 , ze2, ze3, zea    ! 3D workspace
+      !!----------------------------------------------------------------------
+      !! HERE GO VARIABLES USED IN POLCOMS CLEAN UP LATER
+
+      integer i,j,k
+!      real*8 dtmp(n-1)
+      real*8 dtmp(jpkm1)
+      real*8 z1,z2,Rad0,Rad1,Rad2,rD,SurfOut,cp
+      logical first
+      save first
+      data first/.true./
+      !!--------------------------End of POLCOMS variables Note instead of using saves
+      !!--------------------------Could shift this into initial code
+
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_dwl : penetration of the surface solar radiation'
+         IF(lwp) WRITE(numout,*) '~~~~~~~'
+         CALL tra_dwl_init
+         IF( .NOT.ln_tradwl )   RETURN
+      ENDIF
+
+      IF( l_trdtra ) THEN      ! Save ta and sa trends
+         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
+         ztrds(:,:,:) = 0.e0
+      ENDIF
+!--------------------------------------------------------------------
+!  Set transmissivity
+!--------------------------------------------------------------------
+!
+!  Normal value
+!
+!---------------------------------------------------------------------------
+!
+!  Convert Heat fluxes to units used in POL subroutine dwl
+!
+!---------------------------------------------------------------------------
+      cp=3986.0d0
+
+    DO jj = 2, jpj
+         DO ji = fs_2, fs_jpim1
+           qsr(ji,jj)  = qsr(ji,jj)  * (r1_rau0_rcp)
+         ENDDO       !ji
+    ENDDO            !jj
+!--------------------------------------------------------------------------------
+      if ( first ) then
+        rlambda2(:,:) = 0.0
+        first=.false.
+           if ( ln_vary_lambda ) then
+!        do j=1,jesub    ! Original Polcoms style Loop
+!          do i=1,iesub  ! Original Polcoms style Loop
+
+        do jj=2,jpjm1
+          do ji = fs_2, fs_jpim1   ! vector opt.
+
+!            if(ipexb(i,j).ne. 0) then  (Mask, use Tmask instead)
+
+              rlambda2(ji,jj)=-0.033*log(hbatt(ji,jj))+0.2583    ! JIAs formula
+              rlambda2(ji,jj)=max(0.05,rlambda2(ji,jj))     ! limit in deep water
+          enddo ! ji
+        enddo ! jj
+      rlambda = 0.0
+       else
+      rLambda=0.154
+       endif ! If vary lambda
+      endif ! If first
+
+!      do j=1,jesub      ! Original Polcoms Style Loop
+!        do i=1,iesub    ! Original Polcoms Style Loop
+              do jk=2,jpk
+        DO jj=2,jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+
+!--------------------------------------------------------------------
+! Calculate change in temperature
+!--------------------------------------------------------------------
+!
+!        rad0 = hfl_in(i,j)   ! change hfl_in to qsr I assume
+
+               rad0 = qsr(ji,jj)
+                     rD = rLambda2(ji,jj)  +rLambda      !  Transmissivity to be used here
+!       if rlambda 0 then rlambda2 not zer and vica versa 
+
+                z2=gdepw_0(ji,jj,jk-1)    ! grid box is from z=z1 to z=z2
+                z1=gdepw_0(ji,jj,jk)
+
+                Rad2=Rad0*(exp(-z2*rD)) ! radiation entering box
+                Rad1=Rad0*(exp(-z1*rD)) ! radiation leaving box
+
+
+                dtmp(jk)=1.0/(e3t_0(ji,jj,jk))*(Rad2-Rad1) !change in temperature
+                tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + dtmp(jk)
+         enddo  ! ji
+         enddo  ! jj
+             enddo !jk
+
+
+      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
+         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
+         !CEODCALL trd_mod( ztrdt, ztrds, jptra_trd_qsr, 'TRA', kt )
+      ENDIF
+      !                       ! print mean trends (used for debugging)
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
+      !
+   END SUBROUTINE tra_dwl
+
+
+   SUBROUTINE tra_dwl_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_dwl_init  ***
+      !!
+      !! ** Purpose :   Initialization for the penetrative solar radiation for Downwell routine
+      !!
+      !! ** Method  :   The profile of solar radiation within the ocean is set
+      !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
+      !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
+      !!      default values correspond to clear water (type I in Jerlov' 
+      !!      (1968) classification.
+      !!         called by tra_qsr at the first timestep (nit000)
+      !!
+      !! ** Action  : - initialize rn_si0, rn_si1 and rn_abs
+      !!
+      !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk            ! dummy loop indices
+      INTEGER  ::   ios                   ! Local integer output status for namelist read
+      INTEGER  ::   irgb, ierror          ! temporary integer
+      INTEGER  ::   ioptio, nqsr          ! temporary integer
+      REAL(wp) ::   zc0  , zc1            ! temporary scalars
+      REAL(wp) ::   zc2  , zc3  , zchl    !    -         -
+      REAL(wp) ::   zsi0r, zsi1r, zcoef   !    -         -
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zekb, zekg, zekr              ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze0 , ze1 , ze2 , ze3 , zea   ! 3D workspace
+      !!
+      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files
+      TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read
+      NAMELIST/namtra_dwl/  ln_tradwl, ln_vary_lambda
+      !!----------------------------------------------------------------------
+
+      REWIND( numnam_ref )            ! Read Namelist namtra_dwl in reference namelist :
+      READ  ( numnam_ref, namtra_dwl, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in reference namelist', lwp )
+      
+      REWIND( numnam_cfg )            ! Read Namelist namtra_dwl in configuration namelist :
+      READ  ( numnam_cfg, namtra_dwl, IOSTAT = ios, ERR = 902)
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist', lwp )
+      !
+      IF(lwp) THEN                ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'tra_dwl_init : '
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namtra_dwl : set the parameter of penetration'
+         WRITE(numout,*) '      Light penetration (T) or not (F)         ln_tradwl  = ', ln_tradwl
+         WRITE(numout,*) '      Vary Lambda  (T) or not (F))             ln_vary_lambda  = ', ln_vary_lambda
+      ENDIF
+
+   END SUBROUTINE tra_dwl_init
+
+   !!======================================================================
+END MODULE tradwl
Index: branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 10248)
+++ branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 10249)
@@ -58,4 +58,5 @@
 
    REAL(wp) ::   rbcp   ! Brown & Campana parameters for semi-implicit hpg
+   INTEGER  ::   warn_1, warn_2   ! indicators for warning statement
 
    !! * Substitutions
@@ -93,6 +94,6 @@
       INTEGER, INTENT(in) ::   kt    ! ocean time-step index
       !!
-      INTEGER  ::   jk, jn    ! dummy loop indices
-      REAL(wp) ::   zfact     ! local scalars
+      INTEGER  ::   jk, jn, ji, jj     ! dummy loop indices
+      REAL(wp) ::   zfact, zfreeze     ! local scalars
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
       !!----------------------------------------------------------------------
@@ -125,4 +126,33 @@
       ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dtra(:) = 2._wp* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
       ENDIF
+
+#if ( ! defined key_lim3 && ! defined key_lim2 && ! key_cice )
+      IF ( kt == nit000 ) warn_1=0
+      warn_2=0
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF ( tsa(ji,jj,jk,jp_tem) .lt. 0.0 ) THEN
+                  ! calculate freezing point
+                  zfreeze = ( -0.0575_wp + 1.710523E-3 * Sqrt(Abs(tsn(ji,jj,jk,jp_sal)))   & 
+                            - 2.154996E-4 * tsn(ji,jj,jk,jp_sal) ) * tsn(ji,jj,jk,jp_sal) - 7.53E-4 * ( 10.0_wp + fsdept(ji,jj,jk) )
+                  IF ( tsa(ji,jj,jk,jp_tem) .lt. zfreeze ) THEN
+                     tsa(ji,jj,jk,jp_tem)=zfreeze
+                     warn_2=1
+                  ENDIF
+               ENDIF
+            END DO
+         END DO
+      END DO
+      CALL mpp_max(warn_1)
+      CALL mpp_max(warn_2)
+      IF ( (warn_1 == 0) .and. (warn_2 /= 0) ) THEN
+         IF(lwp) THEN
+            CALL ctl_warn( ' Temperatures dropping below freezing point, ', &
+                      &    ' being forced to freezing point, no longer conservative' ) 
+         ENDIF
+         warn_1=1
+      ENDIF
+#endif
 
       ! trends computation initialisation
Index: branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 10248)
+++ branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 10249)
@@ -46,4 +46,5 @@
    LOGICAL , PUBLIC ::   ln_qsr_ice   !: light penetration for ice-model LIM3 (clem)
    INTEGER , PUBLIC ::   nn_chldta    !: use Chlorophyll data (=1) or not (=0)
+   INTEGER , PUBLIC ::   nn_kd490dta  !: use kd490dta data (=1) or not (=0)
    REAL(wp), PUBLIC ::   rn_abs       !: fraction absorbed in the very near surface (RGB & 2 bands)
    REAL(wp), PUBLIC ::   rn_si0       !: very near surface depth of extinction      (RGB & 2 bands)
@@ -54,4 +55,5 @@
    REAL(wp) ::   xsi1r                           !: inverse of rn_si1
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read)
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_kd490 ! structure of input kd490 (file informations, fields read)
    INTEGER, PUBLIC ::   nksr              ! levels below which the light cannot penetrate ( depth larger than 391 m)
    REAL(wp), DIMENSION(3,61) ::   rkrgb   !: tabulated attenuation coefficients for RGB absorption
@@ -306,4 +308,48 @@
             !
          ENDIF
+! slwa
+         IF( nn_kd490dta == 1 ) THEN                      !  use KD490 data read in   !
+            !                                             ! ------------------------- !
+               nksr = jpk - 1
+               !
+               CALL fld_read( kt, 1, sf_kd490 )     ! Read kd490 data and provide it at the current time step
+               !
+               zcoef  = ( 1. - rn_abs )
+               ze0(:,:,1) = rn_abs  * qsr(:,:)
+               ze1(:,:,1) = zcoef * qsr(:,:)
+               zea(:,:,1) =         qsr(:,:)
+               !
+               DO jk = 2, nksr+1
+!CDIR NOVERRCHK
+                  DO jj = 1, jpj
+!CDIR NOVERRCHK   
+                     DO ji = 1, jpi
+                        zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r     )
+                        zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * sf_kd490(1)%fnow(ji,jj,1) )
+                        ze0(ji,jj,jk) = zc0
+                        ze1(ji,jj,jk) = zc1
+                        zea(ji,jj,jk) = ( zc0 + zc1 ) * tmask(ji,jj,jk)
+                     END DO
+                  END DO
+               END DO
+               ! clem: store attenuation coefficient of the first ocean level
+               IF ( ln_qsr_ice ) THEN
+                  DO jj = 1, jpj
+                     DO ji = 1, jpi
+                        zzc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r     )
+                        zzc1 = zcoef  * EXP( - fse3t(ji,jj,1) * sf_kd490(1)%fnow(ji,jj,1) )
+                        fraqsr_1lev(ji,jj) = 1.0 - ( zzc0 + zzc1 ) * tmask(ji,jj,2) 
+                     END DO
+                  END DO
+               ENDIF
+               !
+               DO jk = 1, nksr                                        ! compute and add qsr trend to ta
+                  qsr_hc(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) )
+               END DO
+               zea(:,:,nksr+1:jpk) = 0.e0     ! 
+               CALL iom_put( 'qsr3d', zea )   ! Shortwave Radiation 3D distribution
+               !
+        ENDIF   ! use KD490 data
+!slwa
          !
          !                                        Add to the general trend
@@ -374,7 +420,8 @@
       CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files
       TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read
-      !!
-      NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice,  &
-         &                  nn_chldta, rn_abs, rn_si0, rn_si1
+      TYPE(FLD_N)        ::   sn_kd490 ! informations about the kd490 field to be read
+      !!
+      NAMELIST/namtra_qsr/  sn_chl, sn_kd490, cn_dir, ln_traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice,  &
+         &                  nn_chldta, rn_abs, rn_si0, rn_si1, nn_kd490dta
       !!----------------------------------------------------------------------
 
@@ -409,4 +456,5 @@
          WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0 = ', rn_si0
          WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1 = ', rn_si1
+         WRITE(numout,*) '      read in KD490 data                       nn_kd490dta  = ', nn_kd490dta
       ENDIF
 
@@ -422,4 +470,5 @@
          IF( ln_qsr_2bd  )   ioptio = ioptio + 1
          IF( ln_qsr_bio  )   ioptio = ioptio + 1
+         IF( nn_kd490dta == 1 )   ioptio = ioptio + 1
          !
          IF( ioptio /= 1 ) &
@@ -431,4 +480,5 @@
          IF( ln_qsr_2bd                      )   nqsr =  3
          IF( ln_qsr_bio                      )   nqsr =  4
+         IF( nn_kd490dta == 1                )   nqsr =  5
          !
          IF(lwp) THEN                   ! Print the choice
@@ -438,4 +488,5 @@
             IF( nqsr ==  3 )   WRITE(numout,*) '         2 bands light penetration'
             IF( nqsr ==  4 )   WRITE(numout,*) '         bio-model light penetration'
+            IF( nqsr ==  5 )   WRITE(numout,*) '         KD490 light penetration'
          ENDIF
          !
@@ -447,6 +498,18 @@
          xsi0r = 1.e0 / rn_si0
          xsi1r = 1.e0 / rn_si1
+         IF( nn_kd490dta == 1 ) THEN           !* KD490 data : set sf_kd490 structure
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '        KD490 read in a file'
+            ALLOCATE( sf_kd490(1), STAT=ierror )
+            IF( ierror > 0 ) THEN
+               CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_kd490 structure' )   ;   RETURN
+            ENDIF
+            ALLOCATE( sf_kd490(1)%fnow(jpi,jpj,1)   )
+            IF( sn_kd490%ln_tint )ALLOCATE( sf_kd490(1)%fdta(jpi,jpj,1,2) )
+            !                                        ! fill sf_kd490 with sn_kd490 and control print
+            CALL fld_fill( sf_kd490, (/ sn_kd490 /), cn_dir, 'tra_qsr_init',   &
+               &                                         'Solar penetration function of read KD490', 'namtra_qsr' )
          !                                ! ---------------------------------- !
-         IF( ln_qsr_rgb ) THEN            !  Red-Green-Blue light penetration  !
+         ELSEIF( ln_qsr_rgb ) THEN            !  Red-Green-Blue light penetration  !
             !                             ! ---------------------------------- !
             !
Index: branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 10248)
+++ branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 10249)
@@ -25,4 +25,5 @@
    USE trd_oce         ! trends: ocean variables
    USE trdtra          ! trends manager: tracers 
+   USE tradwl          ! solar radiation penetration (downwell method)
    !
    USE in_out_manager  ! I/O manager
@@ -33,4 +34,7 @@
    USE timing          ! Timing
    USE eosbn2
+#if defined key_asminc   
+   USE asminc          ! Assimilation increment
+#endif
 
    IMPLICIT NONE
@@ -138,5 +142,5 @@
 
 !!gm      IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration
-      IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration
+      IF( .NOT.ln_traqsr .and. .NOT.ln_tradwl ) THEN     ! no solar radiation penetration
          qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns
          qsr(:,:) = 0.e0                     ! qsr set to zero
@@ -278,4 +282,23 @@
          END DO  
       ENDIF
+
+#if defined key_asminc
+! WARNING: THIS MAY WELL NOT BE REQUIRED - WE DON'T WANT TO CHANGE T&S BUT THIS MAY COMPENSATE ANOTHER TERM...
+! Rate of change in e3t for each level is ssh_iau*e3t_0/ht_0
+! Contribution to tsa should be rate of change in level / per m of ocean? (hence the division by fse3t_n)
+      IF( ln_sshinc ) THEN         ! input of heat and salt due to assimilation
+         DO jj = 2, jpj 
+            DO ji = fs_2, fs_jpim1
+               zdep = ssh_iau(ji,jj) / ( ht_0(ji,jj) + 1.0 - ssmask(ji, jj) )
+               DO jk = 1, jpkm1
+                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
+                                        &            + tsn(ji,jj,jk,jp_tem) * zdep * ( e3t_0(ji,jj,jk) / fse3t_n(ji,jj,jk) )
+                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   &
+                                        &            + tsn(ji,jj,jk,jp_sal) * zdep * ( e3t_0(ji,jj,jk) / fse3t_n(ji,jj,jk) )
+               END DO
+            END DO  
+         END DO  
+      ENDIF
+#endif
  
       IF( l_trdtra )   THEN                      ! send trends for further diagnostics
