Index: branches/DEV_R1821_Rivers/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- branches/DEV_R1821_Rivers/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 1937)
+++ branches/DEV_R1821_Rivers/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 1938)
@@ -21,4 +21,6 @@
    USE in_out_manager  ! I/O manager
    USE prtctl          ! Print control
+   USE sbcrnf          ! River runoff  
+   USE sbcmod          ! ln_rnf  
 
    IMPLICIT NONE
@@ -103,6 +105,10 @@
       INTEGER, INTENT(in) ::   kt     ! ocean time-step index
       !!
-      INTEGER  ::   ji, jj                   ! dummy loop indices
-      REAL(wp) ::   zta, zsa, zsrau, zse3t   ! temporary scalars
+      INTEGER  ::   ji, jj, jk           ! dummy loop indices  
+      REAL(wp) ::   zta, zsa             ! temporary scalars, adjustment to temperature and salinity  
+      REAL(wp) ::   azta, azsa           ! temporary scalars, calculations of automatic change to temp & sal due to vvl (done elsewhere)  
+      REAL(wp) ::   zsrau, zse3t, zdep   ! temporary scalars, 1/density, 1/height of box, 1/height of effected water column  
+      REAL(wp) ::   dheat, dsalt          ! total change of temperature and salinity  
+      REAL(wp) ::   tot_sal1, tot_tmp1  
       !!----------------------------------------------------------------------
 
@@ -125,5 +131,5 @@
       IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration
 
-      ! Concentration dillution effect on (t,s)
+      ! Concentration dillution effect on (t,s) due to evapouration, precipitation and qns, but not river runoff  
       DO jj = 2, jpj
          DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -132,10 +138,10 @@
 #endif
             IF( lk_vvl) THEN
-               zta = ro0cpr * qns(ji,jj) * zse3t &                   ! temperature : heat flux
-                &    - emp(ji,jj) * zsrau * tn(ji,jj,1)  * zse3t     ! & cooling/heating effet of EMP flux
+               zta =  ro0cpr * qns(ji,jj) * zse3t &                  ! temperature : heat flux 
+                &    - emp(ji,jj) * zsrau * tn(ji,jj,1) * zse3t      ! & cooling/heating effet of EMP flux 
                zsa = 0.e0                                            ! No salinity concent./dilut. effect
             ELSE
-               zta = ro0cpr * qns(ji,jj) * zse3t     ! temperature : heat flux
-               zsa = emps(ji,jj) * zsrau * sn(ji,jj,1)   * zse3t     ! salinity :  concent./dilut. effect
+               zta =  ro0cpr * qns(ji,jj) * zse3t                    ! temperature : heat flux 
+               zsa =  emps(ji,jj) * zsrau * sn(ji,jj,1) * zse3t      ! salinity :  concent./dilut. effect 
             ENDIF
             ta(ji,jj,1) = ta(ji,jj,1) + zta                          ! add the trend to the general tracer trend
@@ -143,4 +149,73 @@
          END DO
       END DO
+
+      IF ( ln_rnf ) THEN  
+      ! Concentration / dilution effect on (t,s) due to river runoff  
+        DO jj=1,jpj  
+           DO ji=1,jpi  
+              rnf_dep(ji,jj)=0  
+              DO jk=1,rnf_mod_dep(ji,jj)                          ! recalculates rnf_dep to be the depth  
+                rnf_dep(ji,jj)=rnf_dep(ji,jj)+fse3t(ji,jj,jk)    ! in metres to the bottom of the relevant grid box  
+              ENDDO  
+              zdep = 1. / rnf_dep(ji,jj)  
+              zse3t= 1. / fse3t(ji,jj,1)  
+              IF ( rnf_tmp(ji,jj) == -999 )   rnf_tmp(ji,jj)=tn(ji,jj,1)        ! if not specified set runoff temp to be sst  
+  
+              IF ( rnf(ji,jj) .gt. 0.0 ) THEN  
+  
+                IF( lk_vvl) THEN  
+                  !!!indirect flux, concentration or dilution effect  
+                  !!!force a dilution effect in all levels;  
+                  dheat=0.0  
+                  dsalt=0.0  
+                  DO jk=1, rnf_mod_dep(ji,jj)  
+                    zta = -tn(ji,jj,jk) * rnf(ji,jj) * zsrau * zdep  
+                    zsa = -sn(ji,jj,jk) * rnf(ji,jj) * zsrau * zdep  
+                    ta(ji,jj,jk)=ta(ji,jj,jk)+zta  
+                    sa(ji,jj,jk)=sa(ji,jj,jk)+zsa  
+                    dheat=dheat+zta*fse3t(ji,jj,jk)  
+                    dsalt=dsalt+zsa*fse3t(ji,jj,jk)  
+                  ENDDO  
+                  !!!negate this total change in heat and salt content from top level  
+                  zta=-dheat*zse3t  
+                  zsa=-dsalt*zse3t  
+                  ta(ji,jj,1)=ta(ji,jj,1)+zta  
+                  sa(ji,jj,1)=sa(ji,jj,1)+zsa  
+    
+                  !!!direct flux  
+                  zta = rnf_tmp(ji,jj) * rnf(ji,jj) * zsrau * zdep  
+                  zsa = rnf_sal(ji,jj) * rnf(ji,jj) * zsrau * zdep  
+    
+                  DO jk=1, rnf_mod_dep(ji,jj)  
+                    ta(ji,jj,jk) = ta(ji,jj,jk) + zta  
+                    sa(ji,jj,jk) = sa(ji,jj,jk) + zsa  
+                  ENDDO  
+   
+                ELSE  
+                  DO jk=1, rnf_mod_dep(ji,jj)  
+                    zta = ( rnf_tmp(ji,jj)-tn(ji,jj,jk) ) * rnf(ji,jj) * zsrau * zdep  
+                    zsa = ( rnf_sal(ji,jj)-sn(ji,jj,jk) ) * rnf(ji,jj) * zsrau * zdep  
+                    ta(ji,jj,jk) = ta(ji,jj,jk) + zta  
+                    sa(ji,jj,jk) = sa(ji,jj,jk) + zsa  
+                  ENDDO  
+                ENDIF  
+  
+              ELSEIF (rnf(ji,jj) .lt. 0.) THEN   !! for use in baltic when flow is out of domain, want no change in temp and sal  
+  
+                IF( lk_vvl) THEN  
+                  !calculate automatic adjustment to sal and temp due to dilution/concentraion effect   
+                  azta = -tn(ji,jj,1) * rnf(ji,jj) * zsrau * zse3t  
+                  azsa = -sn(ji,jj,1) * rnf(ji,jj) * zsrau * zse3t  
+                  !!!negate this change in sal and temp   
+                  ta(ji,jj,1)=ta(ji,jj,1)-azta  
+                  sa(ji,jj,1)=sa(ji,jj,1)-azsa  
+                ENDIF  
+  
+              ENDIF  
+  
+           ENDDO  
+        ENDDO  
+       
+      ENDIF  
 
       IF( l_trdtra ) THEN      ! save the sbc trends for diagnostic
