Index: branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90
===================================================================
--- branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90	(revision 6805)
+++ branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90	(revision 6808)
@@ -26,8 +26,7 @@
    USE trc_oce         ! share ocean/biogeo variables
    USE phycst          ! physical constants
+   USE ldftra          ! lateral diffusivity coefficients
    USE trabbl          ! active tracer: bottom boundary layer
    USE ldfslp          ! lateral diffusion: iso-neutral slopes
-   USE ldfeiv          ! eddy induced velocity coef. 
-   USE ldftra_oce      ! ocean tracer   lateral physics
    USE zdfmxl          ! vertical physics: mixed layer depth
    USE eosbn2          ! equation of state - Brunt Vaisala frequency
@@ -40,4 +39,5 @@
    USE fldread         ! read input fields 
    USE timing          ! Timing
+   USE wrk_nemo
 
    IMPLICIT NONE
@@ -50,8 +50,7 @@
    LOGICAL            ::   ln_dynwzv    !: vertical velocity read in a file (T) or computed from u/v (F)
    LOGICAL            ::   ln_dynbbl    !: bbl coef read in a file (T) or computed (F)
-   LOGICAL            ::   ln_degrad    !: degradation option enabled or not
    LOGICAL            ::   ln_dynrnf    !: read runoff data in file (T) or set to zero (F)
 
-   INTEGER  , PARAMETER ::   jpfld = 21     ! maximum number of fields to read
+   INTEGER  , PARAMETER ::   jpfld = 15     ! maximum number of fields to read
    INTEGER  , SAVE      ::   jf_tem         ! index of temperature
    INTEGER  , SAVE      ::   jf_sal         ! index of salinity
@@ -68,10 +67,4 @@
    INTEGER  , SAVE      ::   jf_ubl         ! index of u-bbl coef
    INTEGER  , SAVE      ::   jf_vbl         ! index of v-bbl coef
-   INTEGER  , SAVE      ::   jf_ahu         ! index of u-diffusivity coef
-   INTEGER  , SAVE      ::   jf_ahv         ! index of v-diffusivity coef 
-   INTEGER  , SAVE      ::   jf_ahw         ! index of w-diffusivity coef
-   INTEGER  , SAVE      ::   jf_eiu         ! index of u-eiv
-   INTEGER  , SAVE      ::   jf_eiv         ! index of v-eiv
-   INTEGER  , SAVE      ::   jf_eiw         ! index of w-eiv
    INTEGER  , SAVE      ::   jf_fmf         ! index of downward salt flux
 
@@ -92,5 +85,4 @@
 
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -112,11 +104,15 @@
       !!             - interpolates data if needed
       !!----------------------------------------------------------------------
-      !
-      USE oce, ONLY:  zts    => tsa 
+      USE oce, ONLY:  zts    => tsa
       USE oce, ONLY:  zuslp  => ua   , zvslp  => va
-      USE oce, ONLY:  zwslpi => rotb , zwslpj => rotn
-      USE oce, ONLY:  zu     => ub   , zv     => vb,  zw => hdivb
+      USE oce, ONLY:  zu     => ub   , zv     => vb,  zw => rke
       !
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !
+       REAL(wp), DIMENSION(jpi,jpj,jpk     )  :: zwslpi, zwslpj
+!      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts)  :: zts
+!      REAL(wp), DIMENSION(jpi,jpj,jpk     )  :: zuslp, zvslp, zwslpi, zwslpj
+!      REAL(wp), DIMENSION(jpi,jpj,jpk     )  :: zu, zv, zw
+      !
       !
       INTEGER  ::   ji, jj     ! dummy loop indices
@@ -138,5 +134,5 @@
          CALL fld_read( kt, 1, sf_dyn )      !==   read data at kt time step   ==!
          !
-         IF( lk_ldfslp .AND. .NOT.lk_c1d .AND. sf_dyn(jf_tem)%ln_tint ) THEN    ! Computes slopes (here avt is used as workspace)                       
+         IF( l_ldfslp .AND. .NOT.lk_c1d .AND. sf_dyn(jf_tem)%ln_tint ) THEN    ! Computes slopes (here avt is used as workspace)                       
             zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:)   ! temperature
             zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:)   ! salinity 
@@ -162,5 +158,5 @@
       ENDIF
       ! 
-      IF( lk_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)                       
+      IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)                       
          iswap_tem = 0
          IF(  kt /= nit000 .AND. ( sf_dyn(jf_tem)%nrec_a(2) - nrecprev_tem ) /= 0 )  iswap_tem = 1
@@ -267,4 +263,6 @@
       rnf (:,:)        = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1)    ! river runoffs 
 
+      !                                               ! update eddy diffusivity coeff. and/or eiv coeff. at kt
+      IF( l_ldftra_time .OR. l_ldfeiv_time )   CALL ldf_tra( kt ) 
       !                                                      ! bbl diffusive coef
 #if defined key_trabbl && ! defined key_c1d
@@ -276,27 +274,4 @@
          CALL bbl( kt, nit000, 'TRC')
       END IF
-#endif
-#if ( ! defined key_degrad && defined key_traldf_c2d && defined key_traldf_eiv ) && ! defined key_c1d 
-      aeiw(:,:)        = sf_dyn(jf_eiw)%fnow(:,:,1) * tmask(:,:,1)    ! w-eiv
-      !                                                           ! Computes the horizontal values from the vertical value
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            aeiu(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji+1,jj  ) )  ! Average the diffusive coefficient at u- v- points
-            aeiv(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji  ,jj+1) )  ! at u- v- points
-         END DO
-      END DO
-      CALL lbc_lnk( aeiu, 'U', 1. )   ;   CALL lbc_lnk( aeiv, 'V', 1. )    ! lateral boundary condition
-#endif
-      
-#if defined key_degrad && ! defined key_c1d 
-      !                                          ! degrad option : diffusive and eiv coef are 3D
-      ahtu(:,:,:) = sf_dyn(jf_ahu)%fnow(:,:,:) * umask(:,:,:)
-      ahtv(:,:,:) = sf_dyn(jf_ahv)%fnow(:,:,:) * vmask(:,:,:)
-      ahtw(:,:,:) = sf_dyn(jf_ahw)%fnow(:,:,:) * tmask(:,:,:)
-#  if defined key_traldf_eiv 
-      aeiu(:,:,:) = sf_dyn(jf_eiu)%fnow(:,:,:) * umask(:,:,:)
-      aeiv(:,:,:) = sf_dyn(jf_eiv)%fnow(:,:,:) * vmask(:,:,:)
-      aeiw(:,:,:) = sf_dyn(jf_eiw)%fnow(:,:,:) * tmask(:,:,:)
-#  endif
 #endif
       !
@@ -339,12 +314,9 @@
       TYPE(FLD_N), DIMENSION(jpfld) ::  slf_d    ! array of namelist informations on the fields to read
       TYPE(FLD_N) :: sn_tem, sn_sal, sn_mld, sn_emp, sn_ice, sn_qsr, sn_wnd, sn_rnf  ! informations about the fields to be read
-      TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl          !   "                                 "
-      TYPE(FLD_N) :: sn_ahu, sn_ahv, sn_ahw, sn_eiu, sn_eiv, sn_eiw, sn_fmf  !   "                                 "
-      !!----------------------------------------------------------------------
-      !
-      NAMELIST/namdta_dyn/cn_dir, ln_dynwzv, ln_dynbbl, ln_degrad, ln_dynrnf,    &
+      TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl, sn_fmf          !   "                                 "
+      NAMELIST/namdta_dyn/cn_dir, ln_dynwzv, ln_dynbbl, ln_dynrnf,    &
          &                sn_tem, sn_sal, sn_mld, sn_emp, sn_ice, sn_qsr, sn_wnd, sn_rnf,  &
-         &                sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl,          &
-         &                sn_ahu, sn_ahv, sn_ahw, sn_eiu, sn_eiv, sn_eiw, sn_fmf
+         &                sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl, sn_fmf  
+      !!----------------------------------------------------------------------
       !
       REWIND( numnam_ref )              ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data
@@ -365,13 +337,8 @@
          WRITE(numout,*) '      vertical velocity read from file (T) or computed (F) ln_dynwzv  = ', ln_dynwzv
          WRITE(numout,*) '      bbl coef read from file (T) or computed (F)          ln_dynbbl  = ', ln_dynbbl
-         WRITE(numout,*) '      degradation option enabled (T) or not (F)            ln_degrad  = ', ln_degrad
          WRITE(numout,*) '      river runoff option enabled (T) or not (F)           ln_dynrnf  = ', ln_dynrnf
          WRITE(numout,*)
       ENDIF
       ! 
-      IF( ln_degrad .AND. .NOT.lk_degrad ) THEN
-         CALL ctl_warn( 'dta_dyn_init: degradation option requires key_degrad activated ; force ln_degrad to false' )
-         ln_degrad = .FALSE.
-      ENDIF
       IF( ln_dynbbl .AND. ( .NOT.lk_trabbl .OR. lk_c1d ) ) THEN
          CALL ctl_warn( 'dta_dyn_init: bbl option requires key_trabbl activated ; force ln_dynbbl to false' )
@@ -391,50 +358,31 @@
                 jf_rnf = jfld + 1  ;  jfld  = jf_rnf
          slf_d(jf_rnf) = sn_rnf
+         ! Activate runoff key of sbc_oce
+         ln_rnf = .true.
+         WRITE(numout,*) 'dta_dyn : Activate the runoff data structure from ocean core ( force ln_rnf = .true.) '
+         WRITE(numout,*)
       ELSE
          rnf (:,:) = 0._wp
       ENDIF
 
-      !
-      IF( .NOT.ln_degrad ) THEN     ! no degrad option
-         IF( lk_traldf_eiv .AND. ln_dynbbl ) THEN        ! eiv & bbl
-                 jf_ubl  = jfld + 1 ;        jf_vbl  = jfld + 2 ;        jf_eiw  = jfld + 3   ;   jfld = jf_eiw
-           slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl  ;   slf_d(jf_eiw) = sn_eiw
-         ENDIF
-         IF( .NOT.lk_traldf_eiv .AND. ln_dynbbl ) THEN   ! no eiv & bbl
+      IF( ln_dynbbl ) THEN         ! eiv & bbl
                  jf_ubl  = jfld + 1 ;        jf_vbl  = jfld + 2 ;  jfld = jf_vbl
            slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl
-         ENDIF
-         IF( lk_traldf_eiv .AND. .NOT.ln_dynbbl ) THEN   ! eiv & no bbl
-           jf_eiw = jfld + 1 ; jfld = jf_eiw ; slf_d(jf_eiw) = sn_eiw
-         ENDIF
-      ELSE
-              jf_ahu  = jfld + 1 ;        jf_ahv  = jfld + 2 ;        jf_ahw  = jfld + 3  ;  jfld = jf_ahw
-        slf_d(jf_ahu) = sn_ahu  ;   slf_d(jf_ahv) = sn_ahv  ;   slf_d(jf_ahw) = sn_ahw
-        IF( lk_traldf_eiv .AND. ln_dynbbl ) THEN         ! eiv & bbl
-                 jf_ubl  = jfld + 1 ;        jf_vbl  = jfld + 2 ;
-           slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl
-                 jf_eiu  = jfld + 3 ;        jf_eiv  = jfld + 4 ;    jf_eiw  = jfld + 5   ;  jfld = jf_eiw 
-           slf_d(jf_eiu) = sn_eiu  ;   slf_d(jf_eiv) = sn_eiv  ;    slf_d(jf_eiw) = sn_eiw
-        ENDIF
-        IF( .NOT.lk_traldf_eiv .AND. ln_dynbbl ) THEN    ! no eiv & bbl
-                 jf_ubl  = jfld + 1 ;        jf_vbl  = jfld + 2 ;  jfld = jf_vbl
-           slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl
-        ENDIF
-        IF( lk_traldf_eiv .AND. .NOT.ln_dynbbl ) THEN    ! eiv & no bbl
-                 jf_eiu  = jfld + 1 ;         jf_eiv  = jfld + 2 ;    jf_eiw  = jfld + 3   ; jfld = jf_eiw 
-           slf_d(jf_eiu) = sn_eiu  ;   slf_d(jf_eiv) = sn_eiv  ;   slf_d(jf_eiw) = sn_eiw
-        ENDIF
-      ENDIF
-  
+      ENDIF
+
+
       ALLOCATE( sf_dyn(jfld), STAT=ierr )         ! set sf structure
       IF( ierr > 0 ) THEN
          CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' )   ;   RETURN
       ENDIF
+      !                                         ! fill sf with slf_i and control print
+      CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' )
       ! Open file for each variable to get his number of dimension
       DO ifpr = 1, jfld
-         CALL iom_open( TRIM( cn_dir )//TRIM( slf_d(ifpr)%clname ), inum )
-         idv   = iom_varid( inum , slf_d(ifpr)%clvar )  ! id of the variable sdjf%clvar
-         idimv = iom_file ( inum )%ndims(idv)             ! number of dimension for variable sdjf%clvar
-         IF( inum /= 0 )   CALL iom_close( inum )       ! close file if already open
+         CALL fld_clopn( sf_dyn(ifpr), nyear, nmonth, nday )
+         idv   = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar )        ! id of the variable sdjf%clvar
+         idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv)                 ! number of dimension for variable sdjf%clvar
+         IF( sf_dyn(ifpr)%num /= 0 )   CALL iom_close( sf_dyn(ifpr)%num ) ! close file if already open
+         ierr1=0
          IF( idimv == 3 ) THEN    ! 2D variable
                                       ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,1)    , STAT=ierr0 )
@@ -448,8 +396,6 @@
          ENDIF
       END DO
-      !                                         ! fill sf with slf_i and control print
-      CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' )
-      !
-      IF( lk_ldfslp .AND. .NOT.lk_c1d ) THEN                  ! slopes 
+      !
+      IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN                  ! slopes 
          IF( sf_dyn(jf_tem)%ln_tint ) THEN      ! time interpolation
             ALLOCATE( uslpdta (jpi,jpj,jpk,2), vslpdta (jpi,jpj,jpk,2),    &
@@ -478,4 +424,5 @@
    END SUBROUTINE dta_dyn_init
 
+
    SUBROUTINE dta_dyn_wzv( pu, pv, pw )
       !!----------------------------------------------------------------------
@@ -506,19 +453,20 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zu  = pu(ji  ,jj  ,jk) * umask(ji  ,jj  ,jk) * e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk)
-               zu1 = pu(ji-1,jj  ,jk) * umask(ji-1,jj  ,jk) * e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk)
-               zv  = pv(ji  ,jj  ,jk) * vmask(ji  ,jj  ,jk) * e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk)
-               zv1 = pv(ji  ,jj-1,jk) * vmask(ji  ,jj-1,jk) * e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk)
-               zet = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
-               zhdiv(ji,jj,jk) = ( zu - zu1 + zv - zv1 ) * zet 
+               zu  = pu(ji  ,jj  ,jk) * umask(ji  ,jj  ,jk) * e2u(ji  ,jj  ) * e3u_n(ji  ,jj  ,jk)
+               zu1 = pu(ji-1,jj  ,jk) * umask(ji-1,jj  ,jk) * e2u(ji-1,jj  ) * e3u_n(ji-1,jj  ,jk)
+               zv  = pv(ji  ,jj  ,jk) * vmask(ji  ,jj  ,jk) * e1v(ji  ,jj  ) * e3v_n(ji  ,jj  ,jk)
+               zv1 = pv(ji  ,jj-1,jk) * vmask(ji  ,jj-1,jk) * e1v(ji  ,jj-1) * e3v_n(ji  ,jj-1,jk)
+               zhdiv(ji,jj,jk) = ( zu - zu1 + zv - zv1 ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
             END DO
          END DO
       END DO
+      !                              !  update the horizontal divergence with the runoff inflow
+      IF( ln_dynrnf )   zhdiv(:,:,1) = zhdiv(:,:,1) - rnf(:,:) * r1_rau0 / e3t_n(:,:,1)
+      !
       CALL lbc_lnk( zhdiv, 'T', 1. )      ! Lateral boundary conditions on zhdiv
-      !
       ! computation of vertical velocity from the bottom
       pw(:,:,jpk) = 0._wp
       DO jk = jpkm1, 1, -1
-         pw(:,:,jk) = pw(:,:,jk+1) - fse3t(:,:,jk) * zhdiv(:,:,jk)
+         pw(:,:,jk) = pw(:,:,jk+1) - e3t_n(:,:,jk) * zhdiv(:,:,jk)
       END DO
       !
@@ -539,8 +487,8 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpj   ! meridional diapycnal slopes
       !!---------------------------------------------------------------------
-#if defined key_ldfslp && ! defined key_c1d
-      CALL eos    ( pts, rhd, rhop, gdept_0(:,:,:) )
-      CALL eos_rab( pts, rab_n )       ! now local thermal/haline expension ratio at T-points
-      CALL bn2    ( pts, rab_n, rn2  ) ! now    Brunt-Vaisala
+      IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)                       
+         CALL eos    ( pts, rhd, rhop, gdept_0(:,:,:) )
+         CALL eos_rab( pts, rab_n )       ! now local thermal/haline expension ratio at T-points
+         CALL bn2    ( pts, rab_n, rn2  ) ! now    Brunt-Vaisala
 
       ! Partial steps: before Horizontal DErivative
@@ -549,21 +497,20 @@
          &                                        rhd, gru , grv    )  ! of t, s, rd at the last ocean level
       IF( ln_zps .AND.        ln_isfcav)                            &
-         &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
-         &                                        rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
-         &                                 gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the first ocean level
-
-      rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl
-      CALL zdf_mxl( kt )            ! mixed layer depth
-      CALL ldf_slp( kt, rhd, rn2 )  ! slopes
-      puslp (:,:,:) = uslp (:,:,:) 
-      pvslp (:,:,:) = vslp (:,:,:) 
-      pwslpi(:,:,:) = wslpi(:,:,:) 
-      pwslpj(:,:,:) = wslpj(:,:,:) 
-#else
-      puslp (:,:,:) = 0.            ! to avoid warning when compiling
-      pvslp (:,:,:) = 0.
-      pwslpi(:,:,:) = 0.
-      pwslpj(:,:,:) = 0.
-#endif
+         &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
+         &                                        rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level
+
+         rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl
+         CALL zdf_mxl( kt )            ! mixed layer depth
+         CALL ldf_slp( kt, rhd, rn2 )  ! slopes
+         puslp (:,:,:) = uslp (:,:,:) 
+         pvslp (:,:,:) = vslp (:,:,:) 
+         pwslpi(:,:,:) = wslpi(:,:,:) 
+         pwslpj(:,:,:) = wslpj(:,:,:) 
+     ELSE
+         puslp (:,:,:) = 0.            ! to avoid warning when compiling
+         pvslp (:,:,:) = 0.
+         pwslpi(:,:,:) = 0.
+         pwslpj(:,:,:) = 0.
+     ENDIF
       !
    END SUBROUTINE dta_dyn_slp
