Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90	(revision 11393)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90	(revision 11394)
@@ -30,4 +30,5 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE stopack
 
    IMPLICIT NONE
@@ -38,4 +39,11 @@
 
    INTEGER ::   nldf = -2   ! type of lateral diffusion used defined from ln_dynldf_... namlist logicals)
+
+#if defined key_dynldf_c3d
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahm10,ahm20,ahm30,ahm40
+#endif
+#if defined key_dynldf_c2d
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) :: ahm10,ahm20,ahm30,ahm40
+#endif
 
    !! * Substitutions
@@ -67,4 +75,48 @@
          ztrdv(:,:,:) = va(:,:,:) 
       ENDIF
+
+#if defined key_dynldf_c3d
+         IF( kt .eq. nit000 .and. (nn_spp_ahm1+nn_spp_ahm2) .gt. 0 ) THEN
+             ALLOCATE ( ahm10(jpi,jpj,jpk), ahm20(jpi,jpj,jpk) )
+             ALLOCATE ( ahm30(jpi,jpj,jpk), ahm40(jpi,jpj,jpk) )
+             ahm10 = ahm1
+             ahm20 = ahm2
+             ahm30 = ahm3
+             ahm40 = ahm4
+         ENDIF
+#endif
+#if defined key_dynldf_c2d
+         IF( kt .eq. nit000 .and. (nn_spp_ahm1+nn_spp_ahm2) .gt. 0 ) THEN
+             ALLOCATE ( ahm10(jpi,jpj), ahm20(jpi,jpj) )
+             ALLOCATE ( ahm30(jpi,jpj), ahm40(jpi,jpj) )
+             ahm10 = ahm1
+             ahm20 = ahm2
+             ahm30 = ahm3
+             ahm40 = ahm4
+         ENDIF
+#endif
+
+#if defined key_traldf_c3d || defined key_traldf_c2d
+         IF( nn_spp_ahm1 .GT. 0) THEN
+           IF( ln_dynldf_lap ) THEN
+             ahm1 = ahm10
+             CALL spp_ahm(kt,ahm1,nn_spp_ahm1,rn_ahm1_sd,jk_spp_ahm1)
+           ENDIF
+           IF( ln_dynldf_bilap ) THEN
+             ahm3 = ahm30
+             CALL spp_ahm(kt,ahm3,nn_spp_ahm1,rn_ahm1_sd,jk_spp_ahm3)
+           ENDIF
+         ENDIF
+         IF( nn_spp_ahm2 .GT. 0) THEN
+           IF( ln_dynldf_lap ) THEN
+             ahm2 = ahm20
+             CALL spp_ahm(kt,ahm2,nn_spp_ahm2,rn_ahm2_sd,jk_spp_ahm2)
+           ENDIF
+           IF( ln_dynldf_bilap ) THEN
+             ahm4 = ahm40
+             CALL spp_ahm(kt,ahm4,nn_spp_ahm2,rn_ahm2_sd,jk_spp_ahm4)
+           ENDIF
+         ENDIF
+#endif
 
       SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 11393)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 11394)
@@ -22,4 +22,5 @@
    USE wrk_nemo       ! work arrays
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE stopack
 
    IMPLICIT NONE
@@ -154,4 +155,7 @@
                END DO
             END DO
+            
+            IF ( nn_spp_icealb > 0 ) CALL spp_gen( 1, pa_ice_cs(:,:,jl), nn_spp_icealb, rn_icealb_sd, jk_spp_alb, jl )
+                        
          END DO
 
@@ -209,4 +213,7 @@
               END DO
             END DO
+            
+            IF ( nn_spp_icealb > 0 ) CALL spp_gen( 1, pa_ice_os(:,:,jl), nn_spp_icealb, rn_icealb_sd, jk_spp_alb, jl )
+            
          END DO
          ! Effect of the clouds (2d order polynomial)
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 11393)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 11394)
@@ -51,4 +51,5 @@
    USE par_ice_2
 #endif
+   USE stopack
 
    IMPLICIT NONE
@@ -89,4 +90,5 @@
    REAL(wp) ::   rn_efac     ! multiplication factor for evaporation (clem)
    REAL(wp) ::   rn_vfac     ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem)
+   REAL(wp), ALLOCATABLE, SAVE ::   rn_vfac0(:,:) ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem)
    REAL(wp) ::   rn_zqt      ! z(q,t) : height of humidity and temperature measurements
    REAL(wp) ::   rn_zu       ! z(u)   : height of wind measurements
@@ -198,4 +200,12 @@
          sfx(:,:) = 0._wp                          ! salt flux; zero unless ice is present (computed in limsbc(_2).F90)
          !
+         ALLOCATE ( rn_vfac0(jpi,jpj) )
+         rn_vfac0(:,:) = rn_vfac
+         !
+      ENDIF
+
+      IF( nn_spp_relw > 0 ) THEN
+         rn_vfac0 = rn_vfac
+         CALL spp_gen(kt, rn_vfac0, nn_spp_relw, rn_relw_sd, jk_spp_relw )
       ENDIF
 
@@ -291,6 +301,6 @@
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vect. opt.
-            zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  )
-            zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  )
+            zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  )
+            zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  )
          END DO
       END DO
@@ -467,8 +477,10 @@
             DO ji = 2, jpim1   ! B grid : NO vector opt
                ! ... scalar wind at I-point (fld being at T-point)
-               zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)   &
-                  &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - rn_vfac * u_ice(ji,jj)
-               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   &
-                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj)
+               zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)    &
+                  &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) &
+                  &      - rn_vfac0(ji,jj) * u_ice(ji,jj)
+               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)    &
+                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) &
+                  &      - rn_vfac0(ji,jj) * v_ice(ji,jj)
                zwnorm_f = zcoef_wnorm * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f )
                ! ... ice stress at I-point
@@ -476,8 +488,10 @@
                vtau_ice(ji,jj) = zwnorm_f * zwndj_f
                ! ... scalar wind at T-point (fld being at T-point)
-               zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  u_ice(ji,jj+1) + u_ice(ji+1,jj+1)   &
-                  &                                                    + u_ice(ji,jj  ) + u_ice(ji+1,jj  )  )
-               zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  v_ice(ji,jj+1) + v_ice(ji+1,jj+1)   &
-                  &                                                    + v_ice(ji,jj  ) + v_ice(ji+1,jj  )  )
+               zwndi_t = sf(jp_wndi)%fnow(ji,jj,1)                                         &
+                  &      - rn_vfac0(ji,jj) * 0.25 * (  u_ice(ji,jj+1) + u_ice(ji+1,jj+1)   &
+                  &                                  + u_ice(ji,jj  ) + u_ice(ji+1,jj  )  )
+               zwndj_t = sf(jp_wndj)%fnow(ji,jj,1)                                         &
+                  &      - rn_vfac0(ji,jj) * 0.25 * (  v_ice(ji,jj+1) + v_ice(ji+1,jj+1)   &
+                  &                                  + v_ice(ji,jj  ) + v_ice(ji+1,jj  )  )
                wndm_ice(ji,jj)  = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
             END DO
@@ -490,6 +504,6 @@
          DO jj = 2, jpj
             DO ji = fs_2, jpi   ! vect. opt.
-               zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  )
-               zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  )
+               zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  )
+               zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  )
                wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
             END DO
@@ -497,8 +511,10 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vect. opt.
-               utau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          &
-                  &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) )
-               vtau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          &
-                  &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) )
+               utau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )             &
+                  &              * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) &
+                  &                  - rn_vfac0(ji,jj) * u_ice(ji,jj) )
+               vtau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )             &
+                  &              * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) &
+                  &                  - rn_vfac0(ji,jj) * v_ice(ji,jj) )
             END DO
          END DO
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 11393)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 11394)
@@ -53,4 +53,5 @@
    REAL(wp)          , PUBLIC ::   rn_avt_rnf      !: runoffs, value of the additional vertical mixing coef. [m2/s]
    REAL(wp)          , PUBLIC ::   rn_rfact        !: multiplicative factor for runoff
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rn_avt_rnf0
 
    LOGICAL           , PUBLIC ::   l_rnfcpl = .false.       ! runoffs recieved from oasis
@@ -447,4 +448,7 @@
          !                                      !    - mixed upstream-centered (ln_traadv_cen2=T)
          !
+         ALLOCATE ( rn_avt_rnf0(jpi,jpj) )
+         rn_avt_rnf0(:,:) = rn_avt_rnf
+         !
          IF ( ln_rnf_depth )   CALL ctl_warn( 'sbc_rnf_init: increased mixing turned on but effects may already',   &
             &                                              'be spread through depth by ln_rnf_depth'               )
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 11393)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 11394)
@@ -25,4 +25,6 @@
    USE timing         ! Timing
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE stopack
+   USE wrk_nemo       ! Memory Allocation
 
    IMPLICIT NONE
@@ -75,6 +77,6 @@
       REAL(wp) ::   zerp     ! local scalar for evaporation damping
       REAL(wp) ::   zqrp     ! local scalar for heat flux damping
-      REAL(wp) ::   zsrp     ! local scalar for unit conversion of rn_deds factor
       REAL(wp) ::   zerp_bnd ! local scalar for unit conversion of rn_epr_max factor
+      REAL(wp), POINTER, DIMENSION(:,:) :: rn_dqdt_s, zsrp
       INTEGER  ::   ierror   ! return error code
       !!
@@ -95,7 +97,12 @@
             !
             IF( nn_sstr == 1 ) THEN                                   !* Temperature restoring term
+
+               CALL wrk_alloc( jpi, jpj, rn_dqdt_s)
+               rn_dqdt_s=rn_dqdt
+
+               IF( nn_spp_dqdt > 0 ) CALL spp_gen(kt, rn_dqdt_s,nn_spp_dqdt,rn_dqdt_sd,jk_spp_dqdt )
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     zqrp = rn_dqdt * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) )
+                     zqrp = rn_dqdt_s(ji,jj) * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) )
                      qns(ji,jj) = qns(ji,jj) + zqrp
                      qrp(ji,jj) = zqrp
@@ -103,12 +110,15 @@
                END DO
                CALL iom_put( "qrp", qrp )                             ! heat flux damping
+               CALL wrk_dealloc( jpi, jpj, rn_dqdt_s )
             ENDIF
             !
             IF( nn_sssr == 1 ) THEN                                   !* Salinity damping term (salt flux only (sfx))
-               zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s]
+               CALL wrk_alloc( jpi, jpj, zsrp)
+               zsrp = rn_deds
+               IF( nn_spp_dedt > 0 ) CALL spp_gen(kt, zsrp, nn_spp_dedt, rn_dedt_sd, jk_spp_deds )
 !CDIR COLLAPSE
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths
+                     zerp = (zsrp(ji,jj)/rday) * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths
                         &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) ) 
                      sfx(ji,jj) = sfx(ji,jj) + zerp                 ! salt flux
@@ -117,12 +127,15 @@
                END DO
                CALL iom_put( "erp", erp )                             ! freshwater flux damping
+               CALL wrk_dealloc( jpi,jpj, zsrp )
                !
             ELSEIF( nn_sssr == 2 ) THEN                               !* Salinity damping term (volume flux (emp) and associated heat flux (qns)
-               zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s]
+               CALL wrk_alloc( jpi, jpj, zsrp)
+               zsrp = rn_deds
+               IF( nn_spp_dedt > 0 ) CALL spp_gen(kt, zsrp, nn_spp_dedt, rn_dedt_sd, jk_spp_deds )
                zerp_bnd = rn_sssr_bnd / rday                          !       -              -    
 !CDIR COLLAPSE
                DO jj = 1, jpj
                   DO ji = 1, jpi                            
-                     zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths
+                     zerp = (zsrp(ji,jj)/rday) * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths
                         &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) )   &
                         &        / MAX(  sss_m(ji,jj), 1.e-20   )
@@ -134,4 +147,5 @@
                END DO
                CALL iom_put( "erp", erp )                             ! freshwater flux damping
+               CALL wrk_dealloc( jpi,jpj,zsrp )
             ENDIF
             !
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/STO/stopack.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/STO/stopack.F90	(revision 11394)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/STO/stopack.F90	(revision 11394)
@@ -0,0 +1,2812 @@
+!#define NEMO_V34
+
+! Portability to V34 is achieved through (aggressive) precompiler directives:
+! Uncomment the first line for NEMO v3.4
+
+#if defined NEMO_V34
+
+! Change name of trend indexes
+#define  jptra_xad jptra_trd_xad
+#define  jptra_yad jptra_trd_yad
+#define  jptra_zad jptra_trd_zad
+#define  jptra_ldf jptra_trd_ldf
+#define  jptra_zdf jptra_trd_zdf
+#define  jptra_bbc jptra_trd_bbc
+#define  jptra_bbl jptra_trd_bbl
+#define  jptra_npc jptra_trd_npc
+#define  jptra_dmp jptra_trd_dmp
+#define  jptra_qsr jptra_trd_qsr
+#define  jptra_nsr jptra_trd_nsr
+#define  jptra_atf jptra_trd_atf
+
+#define  jptra_sad  -999
+#define  jptra_zdfp -999
+#define  jptra_evd  -999
+
+#define  jpdyn_hpg jpdyn_trd_hpg
+#define  jpdyn_keg jpdyn_trd_keg
+#define  jpdyn_rvo jpdyn_trd_rvo
+#define  jpdyn_pvo jpdyn_trd_pvo
+#define  jpdyn_ldf jpdyn_trd_ldf
+#define  jpdyn_had jpdyn_trd_had
+#define  jpdyn_zad jpdyn_trd_zad
+#define  jpdyn_zdf jpdyn_trd_zdf
+#define  jpdyn_spg jpdyn_trd_spg
+#define  jpdyn_dat jpdyn_trd_dat
+#define  jpdyn_swf jpdyn_trd_swf
+#define  jpdyn_bfr jpdyn_trd_bfr
+
+#define jpdyn_spgflt  -999
+#define jpdyn_spgexp  -999
+#define jpdyn_atf     -999
+#define jpdyn_tau     -999
+#define jpdyn_bfri    -999
+
+! Change name of namelist units
+#define numnam_ref numnam
+#define numnam_cfg numnam
+#define lwm        lwp   
+#define numond     numout
+
+#define wmask      tmask  
+
+#endif
+
+MODULE stopack
+   !!======================================================================
+   !!                       ***  MODULE  stopack ***
+   !!        Calculate and Apply sotchastic physics perturbations
+   !!======================================================================
+   !! History :  1.0  !  2018-02  (A. Storto), Original SPPT code
+   !!                                          for NEMO 3.6
+   !!            2.0  !  2019-05  (A. Storto), upgrades and updates:
+   !!                                          (SPP, SKEB and sea-ice)
+   !!----------------------------------------------------------------------
+   !!   
+   !!   stopack       : Generate stochastic physics perturbations
+   !!   
+   !!                   Method
+   !!                   ======
+   !!                   The module allows users to activate:
+   !!   		- SPPT (Stochastically perturbed parameterization
+   !!			  tendencies )scheme for user-selected trends for
+   !!  			  tracers, momentum and sea-ice 
+   !!   		- SPP (Schastically perturbed parameters) scheme
+   !!			  for some (namelist) parameters
+   !!   		- SEB (Stochastic Energy backscatter)
+   !!			  backscatter energy dissipated numerically or
+   !! 			  through deep convection.
+   !!   
+   !!   
+   !!                   Acknowledgements: C3S funded ERGO project
+   !!   
+   !!----------------------------------------------------------------------
+   USE par_kind
+   USE timing          ! Timing
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain variables 
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE in_out_manager  ! I/O manager
+#ifdef NEMO_V34
+   USE trdmod_oce
+   USE restart
+#else
+   USE trd_oce
+#endif
+   USE iom             ! I/O manager library
+   USE lib_mpp         ! MPP library
+#if defined key_lim3
+   USE ice             ! LIM-3 variables
+#endif
+#if defined key_lim2
+   USE ice_2           ! LIM-2 variables
+#endif
+   USE wrk_nemo
+   USE diaptr
+   USE zdf_oce         
+   USE phycst
+
+   IMPLICIT NONE
+   PRIVATE
+
+   ! Accessible routines
+   PUBLIC tra_sppt_collect
+   PUBLIC dyn_sppt_collect
+   PUBLIC tra_sppt_apply 
+   PUBLIC dyn_sppt_apply 
+   PUBLIC stopack_rst
+   PUBLIC stopack_init
+   PUBLIC stopack_pert
+   PUBLIC sppt_ice
+   PUBLIC spp_aht
+   PUBLIC spp_ahm
+   PUBLIC spp_gen
+   PUBLIC skeb_comp
+   PUBLIC skeb_apply
+
+   ! Main logical switch
+   LOGICAL, PUBLIC :: ln_stopack
+
+   ! Internal switches for SPPT
+   LOGICAL, PUBLIC, SAVE :: ln_sppt_tra = .FALSE.
+   LOGICAL, PUBLIC, SAVE :: ln_sppt_dyn = .FALSE.
+   LOGICAL, PUBLIC, SAVE :: ln_sppt_ice = .FALSE.
+
+   ! SPPT Options
+   INTEGER :: sppt_filter_pass, sppt_step, nn_vwei, &
+   & nn_sppt_step_bound, nn_rndm_freq, nn_deftau
+   REAL(wp), SAVE :: rn_sppt_tau, rn_sppt_bound, rn_distcoast, rn_sppt_stdev
+   REAL(wp), SAVE :: rn_skeb_tau, rn_skeb_stdev
+   REAL(wp), SAVE :: rn_spp_tau, rn_spp_stdev
+   INTEGER :: skeb_filter_pass, spp_filter_pass
+ 
+   ! SPPT Logical switches for individual tendencies
+   LOGICAL :: ln_sppt_taumap, ln_stopack_restart, ln_distcoast, &
+   ln_sppt_traxad, ln_sppt_trayad, ln_sppt_trazad, ln_sppt_trasad, ln_sppt_traldf, &
+   ln_sppt_trazdf, ln_sppt_trazdfp,ln_sppt_traevd, ln_sppt_trabbc, ln_sppt_trabbl, &
+   ln_sppt_tranpc, ln_sppt_tradmp, ln_sppt_traqsr, ln_sppt_transr, ln_sppt_traatf    
+   LOGICAL :: &
+   ln_sppt_dynhpg, ln_sppt_dynspg, ln_sppt_dynkeg, ln_sppt_dynrvo, ln_sppt_dynpvo, ln_sppt_dynzad,&
+   ln_sppt_dynldf, ln_sppt_dynzdf, ln_sppt_dynbfr, ln_sppt_dynatf, ln_sppt_dyntau, ln_sppt_glocon
+   LOGICAL, PUBLIC :: &
+   ln_sppt_icehdf, ln_sppt_icelat, ln_sppt_icezdf
+
+   ! Arrays
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   gauss_n_2d
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   gauss_n_2d_p
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   gauss_n_2d_k
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   coeff_bfr
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gauss_n, zdc
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gauss_nu, gauss_nv
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   spptt, sppts, spptu, spptv
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   gauss_b, sppt_tau, sppt_a, sppt_b
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   skeb_a, skeb_b
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   spp_a, spp_b
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   gauss_bp, gauss_bk
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   g2d_save
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:    ) ::   gauss_w
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zicewrk
+   REAL(wp),              SAVE                   ::   flt_fac,flt_fac_k,flt_fac_p
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   spp_tau
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   skeb_tau
+
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_alb    = 1
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_rhg    = 2
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_relw   = 3
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_dqdt   = 4
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_deds   = 5
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_arnf   = 6
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_geot   = 7
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_qsi0   = 8
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_bfr    = 9
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_aevd   = 10 
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_avt    = 11
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_avm    = 12
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_tkelc  = 13
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_tkedf  = 14
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_tkeds  = 15
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_tkebb  = 16
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_tkefr  = 17
+
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_ahtu   = 18
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_ahtv   = 19
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_ahtw   = 20
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_ahtt   = 21
+
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_ahubbl = 22
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_ahvbbl = 23
+
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_ahm1   = 24
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_ahm2   = 25
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_ahm3   = 26
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_ahm4   = 27
+
+   INTEGER, PARAMETER, PUBLIC :: jk_skeb_dnum  = 31
+   INTEGER, PARAMETER, PUBLIC :: jk_skeb_dcon  = 32
+   INTEGER, PARAMETER, PUBLIC :: jk_skeb_tot   = 33
+
+   INTEGER, PARAMETER, PUBLIC :: jk_sppt_tem   = 41
+   INTEGER, PARAMETER, PUBLIC :: jk_sppt_sal   = 42
+   INTEGER, PARAMETER, PUBLIC :: jk_sppt_uvl   = 43
+   INTEGER, PARAMETER, PUBLIC :: jk_sppt_vvl   = 44
+
+   INTEGER, PARAMETER, PUBLIC :: jk_spp        = 27
+   INTEGER, PARAMETER, PUBLIC :: jk_stpk_tot   = 44
+   REAL(wp), SAVE :: rn_mmax( jk_stpk_tot ) = 0._wp
+
+   INTEGER, SAVE :: numdiag      = 601
+   INTEGER, SAVE :: numrep        = 602
+   INTEGER, SAVE :: lkt
+   
+   ! Randome generator seed
+   INTEGER, SAVE   :: nn_stopack_seed(4)
+   LOGICAL,   SAVE, PUBLIC :: ln_stopack_repr = .TRUE.
+
+   ! SPP switches
+   INTEGER, PUBLIC, SAVE :: nn_spp_bfr,nn_spp_dqdt,nn_spp_dedt,nn_spp_avt,nn_spp_avm,nn_spp_qsi0,nn_spp_relw, nn_spp_arnf,nn_spp_geot,nn_spp_aevd,nn_spp_ahubbl,nn_spp_ahvbbl
+   INTEGER, PUBLIC, SAVE :: nn_spp_tkelc,nn_spp_tkedf,nn_spp_tkeds,nn_spp_tkebb,nn_spp_tkefr
+   INTEGER, PUBLIC, SAVE :: nn_spp_ahtu, nn_spp_ahtv, nn_spp_ahtw, nn_spp_ahtt
+   INTEGER, PUBLIC, SAVE :: nn_spp_ahm1, nn_spp_ahm2
+   ! SPP Sea-ice switches
+   INTEGER, PUBLIC, SAVE :: nn_spp_icealb,nn_spp_icestr
+
+   ! SPP parameters
+   REAL(wp), PUBLIC, SAVE :: rn_bfr_sd, rn_dqdt_sd, rn_dedt_sd, rn_avt_sd, rn_avm_sd, rn_qsi0_sd, rn_relw_sd, rn_arnf_sd, rn_geot_sd, rn_aevd_sd, rn_ahubbl_sd, rn_ahvbbl_sd
+   REAL(wp), PUBLIC, SAVE :: rn_tkelc_sd,rn_tkedf_sd,rn_tkeds_sd,rn_tkebb_sd,rn_tkefr_sd
+   REAL(wp), PUBLIC, SAVE :: rn_ahtu_sd, rn_ahtv_sd, rn_ahtw_sd, rn_ahtt_sd
+   REAL(wp), PUBLIC, SAVE :: rn_ahm1_sd, rn_ahm2_sd
+   REAL(wp), PUBLIC, SAVE :: rn_icestr_sd, rn_icealb_sd
+
+   ! Internal switches for SPP
+   LOGICAL,   SAVE :: ln_spp_own_gauss
+   INTEGER         :: nn_spp
+
+   LOGICAL,  SAVE :: ln_stopack_diags
+
+   LOGICAL,  SAVE :: ln_skeb_own_gauss
+   LOGICAL,  SAVE, PUBLIC :: ln_skeb
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) :: dpsiv, dpsiu
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dnum , dcon
+   INTEGER,  SAVE :: jpri,jprj
+   INTEGER,  SAVE :: nn_dcom_freq = 1
+   REAL(wp), SAVE :: rn_skeb, rn_kh, rn_kc
+   LOGICAL,  SAVE, PUBLIC :: ln_dpsiv = .FALSE.
+   LOGICAL,  SAVE, PUBLIC :: ln_dpsiu = .FALSE.
+   REAL(wp), SAVE :: rn_beta_num, rn_beta_con
+   INTEGER,  SAVE :: nn_dconv = 1
+
+   ! Default/initial seeds
+   INTEGER(KIND=i8) :: x=1234567890987654321_8
+   INTEGER(KIND=i8) :: y=362436362436362436_8
+   INTEGER(KIND=i8) :: z=1066149217761810_8
+   INTEGER(KIND=i8) :: w=123456123456123456_8
+   ! Parameters to generate real random variates
+   REAL(KIND=wp), PARAMETER :: huge64=9223372036854775808.0  ! +1
+   REAL(KIND=wp), PARAMETER :: zero=0.0, half=0.5, one=1.0, two=2.0
+   ! Variables to store 2 Gaussian random numbers with current index (ig)
+   INTEGER(KIND=i8), SAVE :: ig=1
+   REAL(KIND=wp), SAVE :: gran1, gran2
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: bdytra.F90 4292 2013-11-20 16:28:04Z cetlod $ 
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+#if defined NEMO_V34
+   SUBROUTINE tra_sppt_collect( ptrdx, ptrdy, ktrd )
+#else
+   SUBROUTINE tra_sppt_collect( ptrdx, ptrdy, ktrd, kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_sppt_collect ***
+      !!
+      !! ** Purpose :   Collect tracer tendencies (additive)
+      !!                This function is called by the tendency diagnostics 
+      !!                module
+      !!----------------------------------------------------------------------
+      INTEGER                   , INTENT(in   ) ::   kt      ! time step
+#endif
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptrdx   ! Temperature 
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptrdy   ! Salinity
+      INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
+
+      IF( ktrd .eq. jptra_xad .AND. ln_sppt_traxad ) THEN
+           spptt = spptt + ptrdx ; sppts = sppts + ptrdy
+      ENDIF
+      IF( ktrd .eq. jptra_yad .AND. ln_sppt_trayad ) THEN
+           spptt = spptt + ptrdx ; sppts = sppts + ptrdy
+      ENDIF
+      IF( ktrd .eq. jptra_zad .AND. ln_sppt_trazad ) THEN
+           spptt = spptt + ptrdx ; sppts = sppts + ptrdy
+      ENDIF
+      IF( ktrd .eq. jptra_sad .AND. ln_sppt_trasad ) THEN
+           spptt = spptt + ptrdx ; sppts = sppts + ptrdy
+      ENDIF
+      IF( ktrd .eq. jptra_ldf .AND. ln_sppt_traldf ) THEN
+           spptt = spptt + ptrdx ; sppts = sppts + ptrdy
+      ENDIF
+      IF( ktrd .eq. jptra_zdf .AND. ln_sppt_trazdf ) THEN
+           spptt = spptt + ptrdx ; sppts = sppts + ptrdy
+      ENDIF
+      IF( ktrd .eq. jptra_zdfp.AND. ln_sppt_trazdfp) THEN
+           spptt = spptt + ptrdx ; sppts = sppts + ptrdy
+      ENDIF
+      IF( ktrd .eq. jptra_evd .AND. ln_sppt_traevd ) THEN
+           spptt = spptt + ptrdx ; sppts = sppts + ptrdy
+      ENDIF
+      IF( ktrd .eq. jptra_bbc .AND. ln_sppt_trabbc ) THEN
+           spptt = spptt + ptrdx ; sppts = sppts + ptrdy
+      ENDIF
+      IF( ktrd .eq. jptra_bbl .AND. ln_sppt_trabbl ) THEN
+           spptt = spptt + ptrdx ; sppts = sppts + ptrdy
+      ENDIF
+      IF( ktrd .eq. jptra_npc .AND. ln_sppt_tranpc ) THEN
+           spptt = spptt + ptrdx ; sppts = sppts + ptrdy
+      ENDIF
+      IF( ktrd .eq. jptra_dmp .AND. ln_sppt_tradmp ) THEN
+           spptt = spptt + ptrdx ; sppts = sppts + ptrdy
+      ENDIF
+      IF( ktrd .eq. jptra_qsr .AND. ln_sppt_traqsr ) THEN
+           spptt = spptt + ptrdx ; sppts = sppts + ptrdy
+      ENDIF
+      IF( ktrd .eq. jptra_nsr .AND. ln_sppt_transr ) THEN
+           spptt = spptt + ptrdx ; sppts = sppts + ptrdy
+      ENDIF
+      IF( ktrd .eq. jptra_atf .AND. ln_sppt_traatf ) THEN
+           spptt = spptt + ptrdx ; sppts = sppts + ptrdy
+      ENDIF
+
+   END SUBROUTINE tra_sppt_collect
+
+#if defined NEMO_V34
+   SUBROUTINE dyn_sppt_collect( ptrdx, ptrdy, ktrd )
+#else
+   SUBROUTINE dyn_sppt_collect( ptrdx, ptrdy, ktrd, kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dyn_sppt_collect ***
+      !!
+      !! ** Purpose :   Collect momentum tendencies (additive)
+      !!                This function is called by the tendency diagnostics 
+      !!                module
+      !!----------------------------------------------------------------------
+      INTEGER                   , INTENT(in   ) ::   kt      ! time step
+#endif
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptrdx   ! Temperature 
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptrdy   ! Salinity
+      INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
+
+      IF( ktrd .eq. jpdyn_hpg .AND. ln_sppt_dynhpg ) THEN
+           spptu = spptu + ptrdx ; spptv = spptv + ptrdy
+      ENDIF
+      IF( ktrd .eq. jpdyn_spg .AND. ln_sppt_dynspg ) THEN
+           spptu = spptu + ptrdx ; spptv = spptv + ptrdy
+      ENDIF
+      IF( ktrd .eq. jpdyn_spgflt .AND. ln_sppt_dynspg ) THEN
+           spptu = spptu + ptrdx ; spptv = spptv + ptrdy
+      ENDIF
+      IF( ktrd .eq. jpdyn_spgexp .AND. ln_sppt_dynspg ) THEN
+           spptu = spptu + ptrdx ; spptv = spptv + ptrdy
+      ENDIF
+      IF( ktrd .eq. jpdyn_keg .AND. ln_sppt_dynkeg ) THEN
+           spptu = spptu + ptrdx ; spptv = spptv + ptrdy
+      ENDIF
+      IF( ktrd .eq. jpdyn_rvo .AND. ln_sppt_dynrvo ) THEN
+           spptu = spptu + ptrdx ; spptv = spptv + ptrdy
+      ENDIF
+      IF( ktrd .eq. jpdyn_pvo .AND. ln_sppt_dynpvo ) THEN
+           spptu = spptu + ptrdx ; spptv = spptv + ptrdy
+      ENDIF
+      IF( ktrd .eq. jpdyn_zad .AND. ln_sppt_dynzad ) THEN
+           spptu = spptu + ptrdx ; spptv = spptv + ptrdy
+      ENDIF
+      IF( ktrd .eq. jpdyn_ldf .AND. ln_sppt_dynldf ) THEN
+           spptu = spptu + ptrdx ; spptv = spptv + ptrdy
+      ENDIF
+      IF( ktrd .eq. jpdyn_zdf .AND. ln_sppt_dynzdf ) THEN
+           spptu = spptu + ptrdx ; spptv = spptv + ptrdy
+      ENDIF
+      IF( ktrd .eq. jpdyn_bfr .AND. ln_sppt_dynbfr ) THEN
+           spptu = spptu + ptrdx ; spptv = spptv + ptrdy
+      ENDIF
+      IF( ktrd .eq. jpdyn_atf .AND. ln_sppt_dynatf ) THEN
+           spptu = spptu + ptrdx ; spptv = spptv + ptrdy
+      ENDIF
+      IF( ktrd .eq. jpdyn_tau .AND. ln_sppt_dyntau ) THEN
+           spptu = spptu + ptrdx ; spptv = spptv + ptrdy
+      ENDIF
+      IF( ktrd .eq. jpdyn_bfri.AND. ln_sppt_dynbfr ) THEN
+           spptu = spptu + ptrdx ; spptv = spptv + ptrdy
+      ENDIF
+
+   END SUBROUTINE dyn_sppt_collect
+
+   SUBROUTINE stopack_pert( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE stopack_pert ***
+      !!
+      !! ** Purpose :   Update perturbation every nn_rndm_freq timestep
+      !!                Calculate for SPPT, eventually for SPP and SKEB
+      !!                If they are activeated and with different error scales
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) :: kt
+      INTEGER :: ji,jj
+
+      IF( ln_stopack_diags ) lkt = kt
+      IF( MOD( kt - 1, nn_rndm_freq ) == 0 .OR. kt == nit000 ) THEN
+
+       IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) ' Calculating perturbation at timestep ', kt
+         WRITE(numout,*)
+       ENDIF
+
+       ! Generate red noise
+       CALL gaussian_ar1_field ( gauss_n , sppt_filter_pass, gauss_w , gauss_b,&
+        sppt_a, sppt_b, gauss_n_2d,flt_fac,1)
+
+       ! Generate red noise for SPP, SKEB if required
+       IF( ln_spp_own_gauss ) CALL gaussian_ar1_field ( gauss_n , spp_filter_pass, gauss_w , gauss_bp,&
+        spp_a, spp_b, gauss_n_2d_p,flt_fac_p,2)
+       IF( ln_skeb_own_gauss ) CALL gaussian_ar1_field ( gauss_n ,skeb_filter_pass, gauss_w , gauss_bk,&
+        skeb_a, skeb_b, gauss_n_2d_k,flt_fac_k,3)
+
+      ! Staggering on U-,V- grids
+      ! (later should account also for possibly different bottom levels)
+       DO jj=1,jpj-1
+         DO ji=1,jpi-1
+          gauss_nu(ji,jj,:) = 0.5_wp * (gauss_n(ji,jj,:)+gauss_n(ji+1,jj,:))
+          gauss_nv(ji,jj,:) = 0.5_wp * (gauss_n(ji,jj,:)+gauss_n(ji,jj+1,:))
+         ENDDO
+       ENDDO
+       CALL lbc_lnk( gauss_nu, 'U', -1._wp)
+       CALL lbc_lnk( gauss_nv, 'V', -1._wp)
+
+      ENDIF
+
+#ifdef key_iomput
+      ! Output the perturbation field
+      CALL iom_put( "sppt_ar1" , gauss_n )
+      IF(ln_spp_own_gauss)  CALL iom_put( "spp_ar1"  , gauss_n_2d_p)
+      IF(ln_skeb_own_gauss) CALL iom_put( "skeb_ar1" , gauss_n_2d_k)
+#endif
+
+   END SUBROUTINE stopack_pert
+
+#if defined key_lim2
+   SUBROUTINE sppt_ice( kstep, kopt, z1,z2,z3,z4,z5,z6,z7 )
+#else
+   SUBROUTINE sppt_ice( kstep, kopt)
+#endif
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE sppt_ice ***
+      !!
+      !! ** Purpose :   Apply collinear perturbation to ice fields
+      !!                For specific processes coded in LIM2/LIM3 
+      !!----------------------------------------------------------------------
+      !
+      INTEGER, INTENT(IN) :: kstep         ! Start (1) or End (2) of ice routine
+      INTEGER, INTENT(IN) :: kopt          ! Option for sea-ice routine
+#if defined key_lim2
+      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(INOUT) :: z1,z2,z3,z4,z5,z6,z7
+#endif
+#if defined key_lim3 || defined key_lim2
+      INTEGER :: jmt                       ! Number of sea-ice variables (depends on LIM version, process)
+      INTEGER :: jm, jl, jk
+
+#if defined key_lim3
+      jmt = 3*jpl+jpl*nlay_i
+      IF( kopt == 2 ) jmt=jmt+3*jpl+1
+#else
+      IF( kopt == 1 ) jmt=6
+      IF( kopt == 2 ) jmt=7
+      IF( kopt == 3 ) jmt=8
+#endif
+      IF( kopt == 4 ) jmt=2
+
+      IF( kstep == 1 ) THEN        ! Store the before values
+       IF ( ALLOCATED ( zicewrk) ) DEALLOCATE ( zicewrk )
+       ALLOCATE ( zicewrk(jpi,jpj,jmt) )
+       jm=1
+#if defined key_lim3
+       DO jl = 1, jpl
+          zicewrk(:,:,jm) = a_i(:,:,jl) ; jm=jm+1
+          zicewrk(:,:,jm) = v_i(:,:,jl) ; jm=jm+1
+          zicewrk(:,:,jm) =smv_i(:,:,jl); jm=jm+1
+          DO jk = 1, nlay_i
+             zicewrk(:,:,jm) = e_i(:,:,jk,jl) ; jm=jm+1
+          END DO
+       END DO
+       IF( kopt .EQ. 2 ) THEN
+         DO jl = 1, jpl
+           zicewrk(:,:,jm) =oa_i(:,:,jl) ; jm=jm+1
+           zicewrk(:,:,jm) = e_s(:,:,1,jl) ; jm=jm+1
+           zicewrk(:,:,jm) = v_s(:,:,jl) ; jm=jm+1
+         ENDDO
+         zicewrk(:,:,jm) =ato_i(:,:); jm=jm+1
+       ENDIF
+#else
+       IF( kopt .EQ. 1 ) THEN
+         zicewrk(:,:,jm) = frld   ; jm=jm+1
+         zicewrk(:,:,jm) = hsnif  ; jm=jm+1
+         zicewrk(:,:,jm) = hicif  ; jm=jm+1
+         zicewrk(:,:,jm) = tbif(:,:,1)   ; jm=jm+1
+         zicewrk(:,:,jm) = tbif(:,:,2)   ; jm=jm+1
+         zicewrk(:,:,jm) = tbif(:,:,3)   ; jm=jm+1
+       ENDIF
+       IF( kopt .EQ. 2 ) THEN
+         IF(.NOT. PRESENT(z1) ) CALL ctl_stop( 'sppt icehdf problem, step 1')
+         zicewrk(:,:,jm) = z1     ; jm=jm+1
+         zicewrk(:,:,jm) = z2     ; jm=jm+1
+         zicewrk(:,:,jm) = z3     ; jm=jm+1
+         zicewrk(:,:,jm) = z4     ; jm=jm+1
+         zicewrk(:,:,jm) = z5     ; jm=jm+1
+         zicewrk(:,:,jm) = z6     ; jm=jm+1
+         zicewrk(:,:,jm) = z7   
+       ENDIF
+       IF( kopt .EQ. 3 ) THEN
+         zicewrk(:,:,jm) = frld   ; jm=jm+1
+         zicewrk(:,:,jm) = hsnif  ; jm=jm+1
+         zicewrk(:,:,jm) = hicif  ; jm=jm+1
+         zicewrk(:,:,jm) = sist   ; jm=jm+1
+         zicewrk(:,:,jm) = tbif(:,:,1)   ; jm=jm+1
+         zicewrk(:,:,jm) = tbif(:,:,2)   ; jm=jm+1
+         zicewrk(:,:,jm) = tbif(:,:,3)   ; jm=jm+1
+       ENDIF
+#endif
+       IF ( kopt .EQ. 4 ) THEN
+         zicewrk(:,:,jm) = u_ice  ; jm=jm+1
+         zicewrk(:,:,jm) = v_ice  ; jm=jm+1
+       ENDIF
+      ELSEIF ( kstep == 2 ) THEN   ! Add collinear perturbation
+          jm=1
+#if defined key_lim3
+          DO jl = 1, jpl
+            a_i(:,:,jl) = a_i(:,:,jl) + (a_i(:,:,jl)-zicewrk(:,:,jm))*gauss_n_2d ; jm=jm+1 ; CALL lbc_lnk(a_i(:,:,jl), 'T', 1. )
+            v_i(:,:,jl) = v_i(:,:,jl) + (v_i(:,:,jl)-zicewrk(:,:,jm))*gauss_n_2d ; jm=jm+1 ; CALL lbc_lnk(v_i(:,:,jl), 'T', 1. )
+           smv_i(:,:,jl)=smv_i(:,:,jl)+ (smv_i(:,:,jl)-zicewrk(:,:,jm))*gauss_n_2d ; jm=jm+1 ; CALL lbc_lnk(smv_i(:,:,jl), 'T', 1. )
+            DO jk = 1, nlay_i
+               e_i(:,:,jk,jl) = e_i(:,:,jk,jl)+(e_i(:,:,jk,jl)-zicewrk(:,:,jm))*gauss_n_2d ; jm=jm+1 ; CALL lbc_lnk(e_i(:,:,jk,jl), 'T', 1. )
+            END DO
+          END DO
+          IF( kopt .EQ. 2 ) THEN
+            DO jl = 1, jpl
+              oa_i(:,:,jl)=oa_i(:,:,jl)+(oa_i(:,:,jl)-zicewrk(:,:,jm))*gauss_n_2d ; jm=jm+1 ; CALL lbc_lnk(oa_i(:,:,jl), 'T', 1. )
+               e_s(:,:,1,jl)= e_s(:,:,1,jl)+( e_s(:,:,1,jl)-zicewrk(:,:,jm))*gauss_n_2d ; jm=jm+1 ; CALL lbc_lnk(e_s(:,:,1,jl), 'T', 1. )
+               v_s(:,:,jl)= v_s(:,:,jl)+( v_s(:,:,jl)-zicewrk(:,:,jm))*gauss_n_2d ; jm=jm+1 ; CALL lbc_lnk(v_s(:,:,jl), 'T', 1. )
+            ENDDO
+            ato_i(:,:)=ato_i(:,:)+(ato_i(:,:)-zicewrk(:,:,jm))*gauss_n_2d ; jm=jm+1 ; CALL lbc_lnk(ato_i(:,:), 'T', 1. )
+          ENDIF
+          DEALLOCATE ( zicewrk )
+#else
+          IF( kopt .EQ. 1 ) THEN
+            frld   = frld   + gauss_n_2d * ( frld   - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk(frld, 'T', 1. )
+            hsnif  = hsnif  + gauss_n_2d * ( hsnif  - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk(hsnif, 'T', 1. )
+            hicif  = hicif  + gauss_n_2d * ( hicif  - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk(hicif, 'T', 1. )
+        tbif(:,:,1)=tbif(:,:,1) + gauss_n_2d * ( tbif(:,:,1) - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk(tbif(:,:,1), 'T', 1. )
+        tbif(:,:,2)=tbif(:,:,2) + gauss_n_2d * ( tbif(:,:,2) - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk(tbif(:,:,2), 'T', 1. )
+        tbif(:,:,3)=tbif(:,:,3) + gauss_n_2d * ( tbif(:,:,3) - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk(tbif(:,:,3), 'T', 1. )
+          ENDIF
+          IF( kopt .EQ. 2 ) THEN
+            IF(.NOT. PRESENT(z1) ) CALL ctl_stop( 'sppt icehdf problem, step 2')
+            z1 = z1 + gauss_n_2d * ( z1 - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk( z1, 'T', 1. )
+            z2 = z2 + gauss_n_2d * ( z2 - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk( z2, 'T', 1. )
+            z3 = z3 + gauss_n_2d * ( z3 - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk( z3, 'T', 1. )
+            z4 = z4 + gauss_n_2d * ( z4 - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk( z4, 'T', 1. )
+            z5 = z5 + gauss_n_2d * ( z5 - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk( z5, 'T', 1. )
+            z6 = z6 + gauss_n_2d * ( z6 - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk( z6, 'T', 1. )
+            z7 = z7 + gauss_n_2d * ( z7 - zicewrk(:,:,jm) )           ; CALL lbc_lnk( z7, 'T', 1. )
+          ENDIF
+          IF( kopt .EQ. 3 ) THEN
+            frld   = frld   + gauss_n_2d * ( frld   - zicewrk(:,:,jm) ) ; jm=jm+1  ; CALL lbc_lnk(frld, 'T', 1. )
+            hsnif  = hsnif  + gauss_n_2d * ( hsnif  - zicewrk(:,:,jm) ) ; jm=jm+1  ; CALL lbc_lnk(hsnif, 'T', 1. )
+            hicif  = hicif  + gauss_n_2d * ( hicif  - zicewrk(:,:,jm) ) ; jm=jm+1  ; CALL lbc_lnk(hicif, 'T', 1. )
+            sist   = sist   + gauss_n_2d * ( sist   - zicewrk(:,:,jm) ) ; jm=jm+1  ; CALL lbc_lnk(sist, 'T', 1. )
+        tbif(:,:,1)=tbif(:,:,1) + gauss_n_2d * ( tbif(:,:,1) - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk(tbif(:,:,1), 'T', 1. )
+        tbif(:,:,2)=tbif(:,:,2) + gauss_n_2d * ( tbif(:,:,2) - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk(tbif(:,:,2), 'T', 1. )
+        tbif(:,:,3)=tbif(:,:,3) + gauss_n_2d * ( tbif(:,:,3) - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk(tbif(:,:,3), 'T', 1. )
+          ENDIF
+#endif
+   ! EVP or VP rheology common to LIM2-EVP and LIM3
+          IF ( kopt .EQ. 4 ) THEN
+            u_ice = u_ice + (u_ice - zicewrk(:,:,jm) ) * gauss_n_2d ; jm=jm+1
+            v_ice = v_ice + (v_ice - zicewrk(:,:,jm) ) * gauss_n_2d ; jm=jm+1
+#if defined key_lim3 || (  defined key_lim2 && ! defined key_lim2_vp )
+            CALL lbc_lnk( u_ice, 'U', -1. )
+            CALL lbc_lnk( v_ice, 'V', -1. )
+#endif         
+#if defined key_lim2   &&   defined key_lim2_vp
+            CALL lbc_lnk( u_ice(:,1:jpj), 'I', -1. )
+            CALL lbc_lnk( v_ice(:,1:jpj), 'I', -1. )
+#endif         
+          ENDIF
+          DEALLOCATE ( zicewrk )
+      ENDIF
+#endif
+   END SUBROUTINE sppt_ice
+
+   SUBROUTINE spp_gen  ( kt, coeff, nn_type, rn_sd, kspp, klev )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE spp_gen ***
+      !!
+      !! ** Purpose :   Perturbing parameters (generic function)  
+      !!                Given a value of standard deviation, the 2D parameter
+      !!                coeff is perturbed following an additive Normal,
+      !!                a lognormal with mean the original coeff,
+      !!                a lognormal with median the original coeff,
+      !!		or the previous functions but with value bounded [0.1]
+      !!----------------------------------------------------------------------
+   INTEGER, INTENT( in ) :: kt
+   REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj) :: coeff
+   INTEGER, INTENT( in ) ::  nn_type
+   REAL(wp), INTENT( in ) :: rn_sd
+   INTEGER, INTENT( in ) ::  kspp
+   INTEGER, INTENT( in ), OPTIONAL ::  klev
+   REAL(wp), POINTER, DIMENSION(:,:) ::   gauss
+   REAL(wp) :: zsd,xme,mm
+   CHARACTER (LEN=99) :: cstrng
+   INTEGER :: jklev
+
+   CALL wrk_alloc(jpi,jpj,gauss)
+
+   IF( ln_spp_own_gauss ) THEN
+      gauss = gauss_n_2d_p
+   ELSE
+      gauss = rn_spp_stdev * gauss_n_2d / rn_sppt_stdev
+   ENDIF
+
+   IF( nn_type == 1 ) THEN
+       gauss = gauss * rn_sd
+       coeff = coeff * ( 1._wp + gauss )
+   ELSEIF ( nn_type  == 2 ) THEN
+       zsd = rn_sd
+       xme = -0.5_wp*zsd*zsd
+       gauss = gauss * zsd + xme
+       coeff = exp(gauss) * coeff
+   ELSEIF ( nn_type == 3 ) THEN
+       zsd = rn_sd
+       xme = 0._wp
+       gauss = gauss * zsd + xme
+       coeff = exp(gauss) * coeff
+   ELSEIF( nn_type == 4 ) THEN
+       gauss = gauss * rn_sd
+       coeff = coeff * ( 1._wp + gauss )
+       WHERE( coeff > 1._wp ) coeff = 1._wp
+       WHERE( coeff < 0._wp ) coeff = 0._wp
+   ELSEIF ( nn_type  == 5 ) THEN
+       zsd = rn_sd
+       xme = -0.5_wp*zsd*zsd
+       gauss = gauss * zsd + xme
+       coeff = exp(gauss) * coeff
+       WHERE( coeff > 1._wp ) coeff = 1._wp
+       WHERE( coeff < 0._wp ) coeff = 0._wp
+   ELSEIF ( nn_type == 6 ) THEN
+       zsd = rn_sd
+       xme = 0._wp
+       gauss = gauss * zsd + xme
+       coeff = exp(gauss) * coeff
+       WHERE( coeff > 1._wp ) coeff = 1._wp
+       WHERE( coeff < 0._wp ) coeff = 0._wp
+   ELSE
+       CALL ctl_stop( 'spp dqdt wrong option')
+   ENDIF
+
+#ifdef key_iomput
+   write(cstrng,'(A,I2.2)') 'spp_par',kspp
+   IF (iom_use(TRIM(cstrng)) ) CALL iom_put( TRIM(cstrng) , coeff )
+#endif
+
+   IF( ln_stopack_diags ) THEN
+     IF(PRESENT(klev)) THEN
+       jklev = klev
+     ELSE
+       jklev = 0 
+     ENDIF
+     CALL spp_stats(kt,kspp,jklev,coeff)
+   ENDIF
+
+   CALL wrk_dealloc(jpi,jpj,gauss)
+
+   END SUBROUTINE
+
+   SUBROUTINE spp_stats(mt,kp,kl,rcf)
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE spp_stats ***
+      !!
+      !! ** Purpose :   Compute and print basic SPP statistics
+      !!----------------------------------------------------------------------
+   IMPLICIT NONE
+   INTEGER, INTENT(IN) :: mt,kp,kl
+   REAL(wp), INTENT(IN) :: rcf(jpi,jpj)
+   REAL(wp) :: mi,ma
+   CHARACTER(LEN=16) :: cstr = '                '
+   SELECT CASE ( kp ) 
+     CASE( jk_spp_alb   ) 
+       cstr = 'ALBEDO      '
+     CASE( jk_spp_rhg   ) 
+       cstr = 'ICE RHEOLOGY' 
+     CASE( jk_spp_relw  ) 
+       cstr = 'RELATIVE WND' 
+     CASE( jk_spp_dqdt  ) 
+       cstr = 'SST RELAXAT.' 
+     CASE( jk_spp_deds  ) 
+       cstr = 'SSS RELAXAT.' 
+     CASE( jk_spp_arnf  ) 
+       cstr = 'RIVER MIXING' 
+     CASE( jk_spp_geot  ) 
+       cstr = 'GEOTHERM.FLX' 
+     CASE( jk_spp_qsi0  ) 
+       cstr = 'SOLAR EXTIN.' 
+     CASE( jk_spp_bfr   ) 
+       cstr = 'BOTTOM FRICT' 
+     CASE( jk_spp_aevd  ) 
+       cstr = 'EDDY VISCDIF' 
+     CASE( jk_spp_avt   ) 
+       cstr = 'VERT. DIFFUS'
+     CASE( jk_spp_avm   ) 
+       cstr = 'VERT. VISCOS'
+     CASE( jk_spp_tkelc ) 
+       cstr = 'TKE LANGMUIR'
+     CASE( jk_spp_tkedf ) 
+       cstr = 'TKE RN_EDIFF' 
+     CASE( jk_spp_tkeds ) 
+       cstr = 'TKE RN_EDISS' 
+     CASE( jk_spp_tkebb ) 
+       cstr = 'TKE RN_EBB  '
+     CASE( jk_spp_tkefr ) 
+       cstr = 'TKE RN_EFR  '
+     CASE( jk_spp_ahtu  ) 
+       cstr = 'TRALDF AHTU '
+     CASE( jk_spp_ahtv  ) 
+       cstr = 'TRALDF AHTV '
+     CASE( jk_spp_ahtw  ) 
+       cstr = 'TRALDF AHTW '
+     CASE( jk_spp_ahtt  ) 
+       cstr = 'TRALDF AHTT '
+     CASE( jk_spp_ahubbl )
+       cstr = 'BBL DIFFUSU '
+     CASE( jk_spp_ahvbbl )
+       cstr = 'BBL DIFFUSV '
+     CASE( jk_spp_ahm1 )
+       cstr = 'DYNLDF AHM1 '
+     CASE( jk_spp_ahm2 )
+       cstr = 'DYNLDF AHM2 '
+     CASE( jk_spp_ahm3 )
+       cstr = 'DYNLDF AHM3 '
+     CASE( jk_spp_ahm4 )
+       cstr = 'DYNLDF AHM4 '
+     CASE DEFAULT
+       CALL ctl_stop('Unrecognized SPP parameter: add it or turn off diagnostics')
+   END SELECT
+   mi = MINVAL(rcf)
+   ma = MAXVAL(rcf)
+   IF(lk_mpp) CALL mpp_min(mi)
+   IF(lk_mpp) CALL mpp_max(ma)
+   IF(lwp) THEN
+      IF ( kl > 0 ) write(cstr(14:16),'(I3.3)') kl
+      WRITE(numdiag,9300) lkt,cstr,mi,ma
+   ENDIF
+9300  FORMAT(' it :', i8, ' ', A16, ' Min: ',d10.3 , ' Max: ',d10.3)
+   rn_mmax ( kp ) =  MAX ( MAX( abs(mi), ma), rn_mmax ( kp ) )
+   END SUBROUTINE spp_stats
+
+   SUBROUTINE stopack_report
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE spp_stats ***
+      !!
+      !! ** Purpose :  Report at the end of the run
+      !!----------------------------------------------------------------------
+   IMPLICIT NONE
+   INTEGER :: jp, numrep
+   CHARACTER(LEN=16) :: cstr = '                '
+   REAL(wp) :: zmul
+
+   CALL ctl_opn(numrep, 'stopack.report', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+
+   WRITE(numrep,'(A)') ' REPORT from STOPACK stochastic physics package'
+   WRITE(numrep,'(A)')
+   WRITE(numrep,'(A)') ' SPP  part'
+
+   DO jp=1,jk_spp
+     SELECT CASE ( jp ) 
+       CASE( jk_spp_alb   ) 
+         cstr = 'ALBEDO      '
+       CASE( jk_spp_rhg   ) 
+         cstr = 'ICE RHEOLOGY' 
+       CASE( jk_spp_relw  ) 
+         cstr = 'RELATIVE WND' 
+       CASE( jk_spp_dqdt  ) 
+         cstr = 'SST RELAXAT.' 
+       CASE( jk_spp_deds  ) 
+         cstr = 'SSS RELAXAT.' 
+       CASE( jk_spp_arnf  ) 
+         cstr = 'RIVER MIXING' 
+       CASE( jk_spp_geot  ) 
+         cstr = 'GEOTHERM.FLX' 
+       CASE( jk_spp_qsi0  ) 
+         cstr = 'SOLAR EXTIN.' 
+       CASE( jk_spp_bfr   ) 
+         cstr = 'BOTTOM FRICT' 
+       CASE( jk_spp_aevd  ) 
+         cstr = 'EDDY VISCDIF' 
+       CASE( jk_spp_avt   ) 
+         cstr = 'VERT. DIFFUS'
+       CASE( jk_spp_avm   ) 
+         cstr = 'VERT. VISCOS'
+       CASE( jk_spp_tkelc ) 
+         cstr = 'TKE LANGMUIR'
+       CASE( jk_spp_tkedf ) 
+         cstr = 'TKE RN_EDIFF' 
+       CASE( jk_spp_tkeds ) 
+         cstr = 'TKE RN_EDISS' 
+       CASE( jk_spp_tkebb ) 
+         cstr = 'TKE RN_EBB  '
+       CASE( jk_spp_tkefr ) 
+         cstr = 'TKE RN_EFR  '
+       CASE( jk_spp_ahtu  ) 
+         cstr = 'TRALDF AHTU '
+       CASE( jk_spp_ahtv  ) 
+         cstr = 'TRALDF AHTV '
+       CASE( jk_spp_ahtw  ) 
+         cstr = 'TRALDF AHTW '
+       CASE( jk_spp_ahtt  ) 
+         cstr = 'TRALDF AHTT '
+       CASE( jk_spp_ahubbl )
+         cstr = 'BBL DIFFUSU '
+       CASE( jk_spp_ahvbbl )
+         cstr = 'BBL DIFFUSV '
+       CASE( jk_spp_ahm1 )
+         cstr = 'DYNLDF AHM1 '
+       CASE( jk_spp_ahm2 )
+         cstr = 'DYNLDF AHM2 '
+       CASE( jk_spp_ahm3 )
+         cstr = 'DYNLDF AHM3 '
+       CASE( jk_spp_ahm4 )
+         cstr = 'DYNLDF AHM4 '
+       CASE DEFAULT
+         CALL ctl_stop('Unrecognized SPP parameter: add it or turn off diagnostics')
+     END SELECT
+     IF ( rn_mmax(jp) .GT. 0._wp ) &
+     WRITE(numrep,'(A,A,A,D10.3)') ' Maximum of absolute values for parameter ', trim(cstr), ' = ', rn_mmax(jp)
+   ENDDO
+
+   IF(sppt_step .eq. 2 ) THEN
+        zmul = rdt
+      ELSEIF (sppt_step .eq. 0 .or. sppt_step .eq. 1 ) THEN
+        zmul = 1._wp
+   ENDIF
+   rn_mmax(jk_sppt_tem:jk_sppt_vvl) = rn_mmax(jk_sppt_tem:jk_sppt_vvl) * zmul
+
+   WRITE(numrep,'(A)')
+   WRITE(numrep,'(A)') ' SPPT part'
+   WRITE(numrep,'(A,D10.3)') ' Maximum of absolute values for TEM  ', rn_mmax(jk_sppt_tem)
+   IF( rn_mmax(jk_sppt_tem) > 0.5 ) WRITE(numrep,'(A)' ) 'Larger than 0.5, might be too big  '
+   WRITE(numrep,'(A,D10.3)') ' Maximum of absolute values for SAL  ', rn_mmax(jk_sppt_sal)
+   IF( rn_mmax(jk_sppt_sal) > 0.2 ) WRITE(numrep,'(A)' ) 'Larger than 0.2, might be too big  '
+   WRITE(numrep,'(A,D10.3)') ' Maximum of absolute values for U-VEL', rn_mmax(jk_sppt_uvl)
+   IF( rn_mmax(jk_sppt_uvl) > 0.1 ) WRITE(numrep,'(A)' ) 'Larger than 0.1, might be too big  '
+   WRITE(numrep,'(A,D10.3)') ' Maximum of absolute values for V-VEL', rn_mmax(jk_sppt_vvl)
+   IF( rn_mmax(jk_sppt_vvl) > 0.1 ) WRITE(numrep,'(A)' ) 'Larger than 0.1, might be too big  '
+
+   WRITE(numrep,'(A)')
+   WRITE(numrep,'(A)') ' SKEB part'
+   WRITE(numrep,'(A,A,A,D10.3)') ' Maximum of absolute values for NUM  ', trim(cstr), ' = ', rn_mmax(jk_skeb_dnum)
+   WRITE(numrep,'(A)'          ) ' (Perturbation from numerical dissipation)'
+   IF( rn_mmax(jk_skeb_dnum) < 0.5e-4 ) WRITE(numrep,'(A)' ) ' Smaller than 0.5e-4, might be too small'
+   IF( rn_mmax(jk_skeb_dnum) > 0.2e-2 ) WRITE(numrep,'(A)' ) ' Larger  than 0.2e-2, might be too big  '
+   WRITE(numrep,'(A)')
+   WRITE(numrep,'(A,A,A,D10.3)') ' Maximum of absolute values for CONV ', trim(cstr), ' = ', rn_mmax(jk_skeb_dcon)
+   WRITE(numrep,'(A)'          ) ' (Perturbation from convection dissipation)'
+   IF( rn_mmax(jk_skeb_dcon) < 0.5e-4 ) WRITE(numrep,'(A)' ) ' Smaller than 0.5e-4, might be too small'
+   IF( rn_mmax(jk_skeb_dcon) > 0.2e-2 ) WRITE(numrep,'(A)' ) ' Larger  than 0.2e-2, might be too big  '
+   WRITE(numrep,'(A)')
+   WRITE(numrep,'(A,A,A,D10.3)') ' Maximum of absolute values for TOTAL', trim(cstr), ' = ', rn_mmax(jk_skeb_tot)
+   WRITE(numrep,'(A)'          ) ' (Perturbation from total energy dissipation)'
+   IF( rn_mmax(jk_skeb_tot) < 0.5e-4 ) WRITE(numrep,'(A)' ) ' Smaller than 0.5e-4, might be too small'
+   IF( rn_mmax(jk_skeb_tot) > 0.2e-2 ) WRITE(numrep,'(A)' ) ' Larger  than 0.2e-2, might be too big  '
+
+   CLOSE(numrep)
+   CLOSE(numdiag)
+
+   END SUBROUTINE stopack_report
+
+   SUBROUTINE spp_aht ( kt, coeff,nn_type, rn_sd, kspp  )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE spp_aht ***
+      !!
+      !! ** Purpose :   Perturbing diffusivity parameters
+      !!                As spp_gen, but specifically designed for diffusivity
+      !!                where coeff can be 2D or 3D
+      !!----------------------------------------------------------------------
+   INTEGER, INTENT( in ) :: kt
+#if defined key_traldf_c3d
+   REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj,jpk) :: coeff
+#elif defined key_traldf_c2d
+   REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj) :: coeff
+#elif defined key_traldf_c1d
+   REAL(wp), INTENT( inout ), DIMENSION(jpk) :: coeff
+#else
+   REAL(wp), INTENT( inout ) :: coeff
+#endif
+   INTEGER, INTENT( in ) ::  kspp
+   INTEGER, INTENT( in ) ::  nn_type
+   REAL(wp), INTENT( in ) :: rn_sd
+   REAL(wp), POINTER, DIMENSION(:,:) ::   gauss
+   REAL(wp) :: zsd,xme
+   INTEGER :: jk
+
+   CALL wrk_alloc(jpi,jpj,gauss)
+
+   IF( ln_spp_own_gauss ) THEN
+      gauss = gauss_n_2d_p
+   ELSE
+      gauss = rn_spp_stdev * gauss_n_2d / rn_sppt_stdev
+   ENDIF
+
+   IF( nn_type == 1 ) THEN
+       gauss = gauss * rn_sd
+#if defined key_traldf_c3d
+       DO jk=1,jpk
+         coeff(:,:,jk) = coeff(:,:,jk) * ( 1._wp + gauss )
+       ENDDO
+#elif defined key_traldf_c2d
+       coeff = coeff * ( 1._wp + gauss )
+#endif
+   ELSEIF ( nn_type == 2 ) THEN
+       zsd = rn_sd
+       xme = -0.5_wp*zsd*zsd
+       gauss = gauss * zsd + xme
+#if defined key_traldf_c3d
+       DO jk=1,jpk
+         coeff(:,:,jk) = exp(gauss) * coeff(:,:,jk)
+       ENDDOF
+#elif defined key_traldf_c2d
+       coeff = exp(gauss) * coeff
+#endif
+   ELSEIF ( nn_type == 3 ) THEN
+       zsd = rn_sd
+       xme = 0._wp
+       gauss = gauss * zsd + xme
+#if defined key_traldf_c3d
+       DO jk=1,jpk
+         coeff(:,:,jk) = exp(gauss) * coeff(:,:,jk)
+       ENDDOF
+#elif defined key_traldf_c2d
+       coeff = exp(gauss) * coeff
+#endif
+   ELSE
+       CALL ctl_stop( 'spp aht wrong option')
+   ENDIF
+
+   IF( ln_stopack_diags ) THEN
+     CALL spp_stats(kt,kspp,0,coeff)
+   ENDIF
+
+   CALL wrk_dealloc(jpi,jpj,gauss)
+
+   END SUBROUTINE
+
+   SUBROUTINE spp_ahm ( kt, coeff,nn_type, rn_sd, kspp  )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE spp_ahm ***
+      !!
+      !! ** Purpose :   Perturbing viscosity parameters
+      !!                As spp_gen, but specifically designed for viscosity
+      !!                where coeff can be 2D or 3D
+      !!----------------------------------------------------------------------
+   INTEGER, INTENT( in ) :: kt
+#if defined key_dynldf_c3d
+   REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj,jpk) :: coeff
+#elif defined key_dynldf_c2d
+   REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj) :: coeff
+#elif defined key_dynldf_c1d
+   REAL(wp), INTENT( inout ), DIMENSION(jpk) :: coeff
+#else
+   REAL(wp), INTENT( inout ) :: coeff
+#endif
+   INTEGER, INTENT( in ) ::  kspp
+   INTEGER, INTENT( in ) ::  nn_type
+   REAL(wp), INTENT( in ) :: rn_sd
+   REAL(wp), POINTER, DIMENSION(:,:) ::   gauss
+   REAL(wp) :: zsd,xme
+   INTEGER :: jk
+
+   CALL wrk_alloc(jpi,jpj,gauss)
+
+   IF( ln_spp_own_gauss ) THEN
+      gauss = gauss_n_2d_p
+   ELSE
+      gauss = rn_spp_stdev * gauss_n_2d / rn_sppt_stdev
+   ENDIF
+
+   IF( nn_type == 1 ) THEN
+       gauss = gauss * rn_sd
+#if defined key_dynldf_c3d
+       DO jk=1,jpk
+         coeff(:,:,jk) = coeff(:,:,jk) * ( 1._wp + gauss )
+       ENDDO
+#elif defined key_dynldf_c2d
+       coeff = coeff * ( 1._wp + gauss )
+#endif
+   ELSEIF ( nn_type == 2 ) THEN
+       zsd = rn_sd
+       xme = -0.5_wp*zsd*zsd
+       gauss = gauss * zsd + xme
+#if defined key_dynldf_c3d
+       DO jk=1,jpk
+         coeff(:,:,jk) = exp(gauss) * coeff(:,:,jk)
+       ENDDOF
+#elif defined key_dynldf_c2d
+       coeff = exp(gauss) * coeff
+#endif
+   ELSEIF ( nn_type == 3 ) THEN
+       zsd = rn_sd
+       xme = 0._wp
+       gauss = gauss * zsd + xme
+#if defined key_dynldf_c3d
+       DO jk=1,jpk
+         coeff(:,:,jk) = exp(gauss) * coeff(:,:,jk)
+       ENDDOF
+#elif defined key_dynldf_c2d
+       coeff = exp(gauss) * coeff
+#endif
+   ELSE
+       CALL ctl_stop( 'spp ahm wrong option')
+   ENDIF
+
+   IF( ln_stopack_diags ) THEN
+     CALL spp_stats(kt,kspp,0,coeff)
+   ENDIF
+
+   CALL wrk_dealloc(jpi,jpj,gauss)
+
+   END SUBROUTINE
+
+   SUBROUTINE tra_sppt_apply ( kt , ks)
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_sppt_apply ***
+      !!
+      !! ** Purpose :   Apply perturbation to tracers
+      !! 		Step 0/1 for tendencies, step 2 for variables
+      !! 		after timestepping
+      !!----------------------------------------------------------------------
+
+      INTEGER, INTENT( in ) :: kt, ks     ! Main time step counter
+      REAL(wp) :: mi,ma
+
+      IF( ks .NE. sppt_step ) RETURN
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) ' Applying sppt on tracers at timestep ', kt
+         WRITE(numout,*)
+      ENDIF
+
+      ! Modify the tendencies
+      IF(sppt_step .eq. 2 ) THEN
+        ! At the end of the timestepping, after array swapping
+        tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) + gauss_n * spptt * rdt
+        tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) + gauss_n * sppts * rdt
+        tsb(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) + gauss_n * spptt * rdt
+        tsb(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) + gauss_n * sppts * rdt
+        IF( ln_sppt_glocon ) CALL sppt_glocon( tsn(:,:,:,jp_tem), 'T' )
+        IF( ln_sppt_glocon ) CALL sppt_glocon( tsn(:,:,:,jp_sal), 'T' )
+        IF( ln_sppt_glocon ) CALL sppt_glocon( tsb(:,:,:,jp_tem), 'T' )
+        IF( ln_sppt_glocon ) CALL sppt_glocon( tsb(:,:,:,jp_sal), 'T' )
+        CALL lbc_lnk( tsb(:,:,:,jp_tem) , 'T', 1._wp)
+        CALL lbc_lnk( tsb(:,:,:,jp_sal) , 'T', 1._wp)
+        CALL lbc_lnk( tsn(:,:,:,jp_tem) , 'T', 1._wp)
+        CALL lbc_lnk( tsn(:,:,:,jp_sal) , 'T', 1._wp)
+      ELSEIF (sppt_step .eq. 0 .or. sppt_step .eq. 1 ) THEN
+        ! At the beginning / before vertical diffusion
+        tsa(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) + gauss_n * spptt
+        tsa(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) + gauss_n * sppts
+        IF( ln_sppt_glocon ) CALL sppt_glocon( tsa(:,:,:,jp_tem), 'T' )
+        IF( ln_sppt_glocon ) CALL sppt_glocon( tsa(:,:,:,jp_sal), 'T' )
+        CALL lbc_lnk( tsa(:,:,:,jp_tem) , 'T', 1._wp)
+        CALL lbc_lnk( tsa(:,:,:,jp_sal) , 'T', 1._wp)
+      ENDIF
+
+      IF( ln_stopack_diags ) THEN
+        mi = MINVAL(spptt)
+        ma = MAXVAL(spptt)
+        IF(lk_mpp) CALL mpp_min(mi)
+        IF(lk_mpp) CALL mpp_max(ma)
+        IF(lwp) WRITE(numdiag,9301) lkt,'SPPT TEMPERATURE' ,mi,ma
+        rn_mmax ( jk_sppt_tem ) =  MAX ( MAX( abs(mi), ma), rn_mmax ( jk_sppt_tem ) )
+
+        mi = MINVAL(sppts)
+        ma = MAXVAL(sppts)
+        IF(lk_mpp) CALL mpp_min(mi)
+        IF(lk_mpp) CALL mpp_max(ma)
+        IF(lwp) WRITE(numdiag,9301) lkt,'SPPT SALINITY   ' ,mi,ma
+9301  FORMAT(' it :', i8, ' ', A16, ' Min: ',d10.3 , ' Max: ',d10.3)
+        rn_mmax ( jk_sppt_sal ) =  MAX ( MAX( abs(mi), ma), rn_mmax ( jk_sppt_sal ) )
+      ENDIF
+      ! Reset the tendencies
+      spptt = 0._wp  ; sppts = 0._wp
+
+   END SUBROUTINE tra_sppt_apply
+
+   SUBROUTINE dyn_sppt_apply ( kt , ks)
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dyn_sppt_apply ***
+      !!
+      !! ** Purpose :   Apply perturbation to momentum
+      !! 		Step 0/1 for tendencies, step 2 for variables
+      !! 		after timestepping
+      !!----------------------------------------------------------------------
+
+      INTEGER, INTENT( in ) :: kt, ks     ! Main time step counter
+      REAL(wp) :: mi,ma
+
+      IF( ks .NE. sppt_step ) RETURN
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) ' Applying sppt on currents at timestep ', kt
+         WRITE(numout,*)
+      ENDIF
+
+      ! Modify the tendencies
+      IF(sppt_step .eq. 2 ) THEN
+        ! At the end of the timestepping, after array swapping
+        un(:,:,:) = un(:,:,:) + gauss_nu* spptu * rdt
+        vn(:,:,:) = vn(:,:,:) + gauss_nv* spptv * rdt
+        ub(:,:,:) = ub(:,:,:) + gauss_nu* spptu * rdt
+        vb(:,:,:) = vb(:,:,:) + gauss_nv* spptv * rdt
+        IF( ln_sppt_glocon ) CALL sppt_glocon( ub, 'U' )
+        IF( ln_sppt_glocon ) CALL sppt_glocon( vb, 'V' )
+        IF( ln_sppt_glocon ) CALL sppt_glocon( un, 'U' )
+        IF( ln_sppt_glocon ) CALL sppt_glocon( vn, 'V' )
+        CALL lbc_lnk( un , 'U', -1._wp)
+        CALL lbc_lnk( vn , 'V', -1._wp)
+        CALL lbc_lnk( ub , 'U', -1._wp)
+        CALL lbc_lnk( vb , 'V', -1._wp)
+      ELSEIF (sppt_step .eq. 0 .or. sppt_step .eq. 1 ) THEN
+        ! At the beginning / before vertical diffusion
+        ua(:,:,:) = ua(:,:,:) + gauss_nu* spptu
+        va(:,:,:) = va(:,:,:) + gauss_nv* spptv
+        IF( ln_sppt_glocon ) CALL sppt_glocon( ua, 'U' )
+        IF( ln_sppt_glocon ) CALL sppt_glocon( va, 'V' )
+        CALL lbc_lnk( ua , 'U', -1._wp)
+        CALL lbc_lnk( va , 'V', -1._wp)
+      ENDIF
+
+      IF( ln_stopack_diags ) THEN
+        mi = MINVAL(spptu)
+        ma = MAXVAL(spptu)
+        IF(lk_mpp) CALL mpp_min(mi)
+        IF(lk_mpp) CALL mpp_max(ma)
+        IF(lwp) WRITE(numdiag,9301) lkt,'SPPT ZONAL CURR.' ,mi,ma
+        rn_mmax ( jk_sppt_uvl ) =  MAX ( MAX( abs(mi), ma), rn_mmax ( jk_sppt_uvl ) )
+
+        mi = MINVAL(spptv)
+        ma = MAXVAL(spptv)
+        IF(lk_mpp) CALL mpp_min(mi)
+        IF(lk_mpp) CALL mpp_max(ma)
+        IF(lwp) WRITE(numdiag,9301) lkt,'SPPT MERID CURR.' ,mi,ma
+9301  FORMAT(' it :', i8, ' ', A16, ' Min: ',d10.3 , ' Max: ',d10.3)
+        rn_mmax ( jk_sppt_vvl ) =  MAX ( MAX( abs(mi), ma), rn_mmax ( jk_sppt_vvl ) )
+      ENDIF
+
+      ! Reset the tendencies
+      spptu = 0._wp  ; spptv = 0._wp
+
+   END SUBROUTINE dyn_sppt_apply
+
+   SUBROUTINE sppt_glocon ( zts, cd_type )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE sppt_glocon ***
+      !!
+      !! ** Purpose :   Apply global conservation constraint
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(INOUT) :: zts(jpi,jpj,jpk)
+      CHARACTER(len=1), INTENT(IN) :: cd_type
+      INTEGER :: ji,jj,jk
+      REAL(wp) :: zv, vl
+      ! Calculate volume tendencies and renormalize
+      ! Note: sshn should be staggered before being used.
+      SELECT CASE ( cd_type )
+               CASE ( 'T' )  
+                jk=1
+                zv = SUM( tmask_i(:,:)*tmask(:,:,jk)*e1t(:,:)*e2t(:,:)*sshn(:,:)*zts(:,:,1) )
+                vl = SUM( tmask_i(:,:)*tmask(:,:,jk)*e1t(:,:)*e2t(:,:)*sshn(:,:) )
+                DO jk = 1, jpkm1
+                   DO jj=1,jpj
+                      DO ji=1,jpi
+                        zv = zv + SUM( tmask_i(:,:)*tmask(:,:,jk)*e1t(:,:)*e2t(:,:)*fse3t(:,:,jk)*zts(:,:,jk) )
+                        vl = vl + SUM( tmask_i(:,:)*tmask(:,:,jk)*e1t(:,:)*e2t(:,:)*fse3t(:,:,jk) )
+                      END DO
+                   END DO
+                END DO
+                IF(lk_mpp) CALL mpp_sum(zv)
+                IF(lk_mpp) CALL mpp_sum(vl)
+                zv = zv / vl
+                zts = zts - zv*tmask
+               CASE ( 'U' )
+                jk=1
+#if defined NEMO_V34
+                zv = SUM( umask(:,:,jk)*e1u(:,:)*e2u(:,:)*sshn(:,:)*zts(:,:,1) )
+                vl = SUM( umask(:,:,jk)*e1u(:,:)*e2u(:,:)*sshn(:,:) )
+#else
+                zv = SUM( umask_i(:,:)*umask(:,:,jk)*e1u(:,:)*e2u(:,:)*sshn(:,:)*zts(:,:,1) )
+                vl = SUM( umask_i(:,:)*umask(:,:,jk)*e1u(:,:)*e2u(:,:)*sshn(:,:) )
+#endif
+                DO jk = 1, jpkm1
+                   DO jj=1,jpj
+                      DO ji=1,jpi
+#if defined NEMO_V34
+                        zv = zv + SUM( umask(:,:,jk)*e1u(:,:)*e2u(:,:)*fse3u(:,:,jk)*zts(:,:,jk) )
+                        vl = vl + SUM( umask(:,:,jk)*e1u(:,:)*e2u(:,:)*fse3u(:,:,jk) )
+#else
+                        zv = zv + SUM( umask_i(:,:)*umask(:,:,jk)*e1u(:,:)*e2u(:,:)*fse3u(:,:,jk)*zts(:,:,jk) )
+                        vl = vl + SUM( umask_i(:,:)*umask(:,:,jk)*e1u(:,:)*e2u(:,:)*fse3u(:,:,jk) )
+#endif
+                      END DO
+                   END DO
+                END DO
+                IF(lk_mpp) CALL mpp_sum(zv)
+                IF(lk_mpp) CALL mpp_sum(vl)
+                zv = zv / vl
+                zts = zts - zv*umask
+               CASE ( 'V' )
+                jk=1
+#if defined NEMO_V34
+                zv = SUM( vmask(:,:,jk)*e1v(:,:)*e2v(:,:)*sshn(:,:)*zts(:,:,1) )
+                vl = SUM( vmask(:,:,jk)*e1v(:,:)*e2v(:,:)*sshn(:,:) )
+#else
+                zv = SUM( vmask_i(:,:)*vmask(:,:,jk)*e1v(:,:)*e2v(:,:)*sshn(:,:)*zts(:,:,1) )
+                vl = SUM( vmask_i(:,:)*vmask(:,:,jk)*e1v(:,:)*e2v(:,:)*sshn(:,:) )
+#endif
+                DO jk = 1, jpkm1
+                   DO jj=1,jpj
+                      DO ji=1,jpi
+#if defined NEMO_V34
+                        zv = zv + SUM( vmask(:,:,jk)*e1v(:,:)*e2v(:,:)*fse3v(:,:,jk)*zts(:,:,jk) )
+                        vl = vl + SUM( vmask(:,:,jk)*e1v(:,:)*e2v(:,:)*fse3v(:,:,jk) )
+#else
+                        zv = zv + SUM( vmask_i(:,:)*vmask(:,:,jk)*e1v(:,:)*e2v(:,:)*fse3v(:,:,jk)*zts(:,:,jk) )
+                        vl = vl + SUM( vmask_i(:,:)*vmask(:,:,jk)*e1v(:,:)*e2v(:,:)*fse3v(:,:,jk) )
+#endif
+                      END DO
+                   END DO
+                END DO
+                IF(lk_mpp) CALL mpp_sum(zv)
+                IF(lk_mpp) CALL mpp_sum(vl)
+                zv = zv / vl
+                zts = zts - zv*vmask
+              ENDSELECT
+
+   END SUBROUTINE sppt_glocon
+
+   SUBROUTINE gaussian_ar1_field ( gn, nk, wei, gb, a, b,gn0, fltf, istep)
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE gaussian_ar1_field ***
+      !!
+      !! ** Purpose :   Generate correlated perturbation field
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: a, b
+      REAL(wp), DIMENSION(jpk)    , INTENT(INOUT) :: wei
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: gb
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(INOUT) :: gn
+      REAL(wp), DIMENSION(jpi,jpj    ), INTENT(OUT) :: gn0
+      INTEGER, INTENT(IN) :: nk,istep
+      REAL(wp), INTENT(IN) :: fltf
+      REAL(wp) :: g2d(jpi,jpj)
+      INTEGER :: jf, jk, ji, jj, nbl
+
+      ! Random noise on 2d field
+      IF ( istep == 1 ) THEN
+         CALL sppt_rand2d( g2d ) 
+         CALL lbc_lnk( g2d , 'T', 1._wp)
+         g2d_save = g2d
+         IF( nn_sppt_step_bound .EQ. 1 ) THEN
+           WHERE ( g2d < -ABS(rn_sppt_bound) ) g2d = -ABS(rn_sppt_bound)
+           WHERE ( g2d >  ABS(rn_sppt_bound) ) g2d =  ABS(rn_sppt_bound)
+         ENDIF
+      ELSEIF ( istep == 2 ) THEN
+         g2d = rn_spp_stdev * g2d_save / rn_sppt_stdev
+      ELSEIF ( istep == 3 ) THEN
+         g2d = rn_skeb_stdev * g2d_save / rn_sppt_stdev
+      ENDIF
+   
+      ! Laplacian filter and re-normalization
+      DO jf = 1, nk
+         CALL sppt_flt( g2d )
+      END DO
+      g2d = g2d * fltf
+
+#ifdef key_iommput
+      ! Output the random field
+      IF(istep==1) THEN
+        CALL iom_put( "sppt_ran" , g2d )
+      ELSEIF (istep==2) THEN
+        CALL iom_put( "spp_ran" , g2d )
+      ELSEIF (istep==3) THEN
+        CALL iom_put( "skeb_ran" , g2d )
+      ENDIF
+#endif
+   
+      ! AR(1) process and array swap
+      g2d = a*gb + b*g2d
+      gb = g2d
+      gn0 = g2d * zdc(:,:,1)
+
+      IF (istep==2 .or. istep==3) RETURN
+
+      ! From 2- to 3-d and vertical weigth
+      IF(nn_vwei .eq. 2) THEN
+       DO jj=1,jpj
+        DO ji=1,jpi
+         nbl = mbathy(ji,jj)
+         wei(:) = 0._wp
+         IF(nbl>0) THEN
+           wei(1:nbl) = 1._wp
+           wei(1) = 0._wp
+           wei(2) = 0.5_wp
+           wei(nbl-1) = 0.5_wp
+           wei(nbl) = 0._wp
+           DO jk=1,jpk
+             gn(ji,jj,jk) = g2d(ji,jj) * wei(jk) * zdc(ji,jj,jk)
+           ENDDO
+         ENDIF
+        ENDDO
+       ENDDO
+      ELSEIF(nn_vwei .eq. 3) THEN
+       DO jj=1,jpj
+        DO ji=1,jpi
+         nbl = mbathy(ji,jj)
+         wei(:) = 0._wp
+         IF(nbl>0) THEN
+           wei(1:nbl) = 1._wp
+           wei(nbl-1) = 0.5_wp
+           wei(nbl) = 0._wp
+           DO jk=1,jpk
+             gn(ji,jj,jk) = g2d(ji,jj) * wei(jk) * zdc(ji,jj,jk)
+           ENDDO
+         ENDIF
+        ENDDO
+       ENDDO
+       ELSE
+        DO jk=1,jpk
+          gn(:,:,jk) = g2d * wei(jk) * zdc(:,:,jk)
+        ENDDO
+      ENDIF
+      
+      ! Bound
+      IF( nn_sppt_step_bound .EQ. 2 ) THEN
+        WHERE ( gn < -ABS(rn_sppt_bound) ) gn = -ABS(rn_sppt_bound)
+        WHERE ( gn >  ABS(rn_sppt_bound) ) gn =  ABS(rn_sppt_bound)
+        WHERE ( gn0< -ABS(rn_sppt_bound) ) gn0= -ABS(rn_sppt_bound)
+        WHERE ( gn0>  ABS(rn_sppt_bound) ) gn0=  ABS(rn_sppt_bound)
+      ENDIF
+
+   END SUBROUTINE gaussian_ar1_field
+   !
+   SUBROUTINE stopack_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE stopack_init ***
+      !!
+      !! ** Purpose :   Initialize stopack
+      !!----------------------------------------------------------------------
+      !!
+      INTEGER  :: ios, inum
+      INTEGER  :: nn_sppt_tra, nn_sppt_dyn, nn_spp_aht, nn_sppt_ice
+      INTEGER  :: ivals(8)
+      INTEGER(KIND=8)     ::   ziseed(4)           ! RNG seeds in integer type
+      REAL(KIND=8)        ::   zrseed(4)           ! RNG seeds in real type
+      REAL(KIND=8)        ::   rdt_sppt
+      !!
+      NAMELIST/namstopack/ ln_stopack, ln_sppt_taumap, rn_sppt_tau, &
+      ln_stopack_restart, sppt_filter_pass, rn_sppt_bound, sppt_step, nn_deftau,rn_sppt_stdev,&
+      ln_distcoast, rn_distcoast, nn_vwei, nn_stopack_seed, nn_sppt_step_bound, nn_rndm_freq, &
+      ln_sppt_glocon, ln_stopack_repr, &
+      rn_icestr_sd, rn_icealb_sd, &
+      ln_sppt_traxad, ln_sppt_trayad, ln_sppt_trazad, ln_sppt_trasad, ln_sppt_traldf, &
+      ln_sppt_trazdf, ln_sppt_trazdfp,ln_sppt_traevd, ln_sppt_trabbc, ln_sppt_trabbl, &
+      ln_sppt_tranpc, ln_sppt_tradmp, ln_sppt_traqsr, ln_sppt_transr, ln_sppt_traatf ,&
+      ln_sppt_dynhpg, ln_sppt_dynspg, ln_sppt_dynkeg, ln_sppt_dynrvo, ln_sppt_dynpvo, &
+      ln_sppt_dynzad, ln_sppt_dynldf, ln_sppt_dynzdf, ln_sppt_dynbfr, ln_sppt_dynatf, ln_sppt_dyntau,&
+      ln_sppt_icehdf, ln_sppt_icelat, ln_sppt_icezdf, & !ln_sppt_icealb, ln_sppt_icestr
+      nn_spp_icealb, nn_spp_icestr,&
+      spp_filter_pass,rn_spp_stdev,rn_spp_tau,&
+      nn_spp_bfr,nn_spp_dqdt,nn_spp_dedt,nn_spp_avt,nn_spp_avm,nn_spp_qsi0,nn_spp_ahtu,nn_spp_ahtv,nn_spp_ahm1,nn_spp_ahm2,&
+      nn_spp_ahtw,nn_spp_ahtt,nn_spp_relw,nn_spp_arnf,nn_spp_geot,nn_spp_aevd,nn_spp_ahubbl,nn_spp_ahvbbl,&
+      nn_spp_tkelc,nn_spp_tkedf,nn_spp_tkeds,nn_spp_tkebb,nn_spp_tkefr,&
+      rn_bfr_sd,rn_dqdt_sd,rn_dedt_sd,rn_avt_sd,rn_avm_sd,rn_qsi0_sd,rn_ahtu_sd,rn_ahtv_sd,rn_ahm1_sd,rn_ahm2_sd,&
+      rn_ahtw_sd,rn_ahtt_sd, rn_relw_sd, rn_arnf_sd,rn_geot_sd, rn_aevd_sd,rn_ahubbl_sd,rn_ahvbbl_sd,&
+      rn_tkelc_sd,rn_tkedf_sd,rn_tkeds_sd,rn_tkebb_sd,rn_tkefr_sd,&
+      ln_skeb,rn_skeb,nn_dcom_freq,rn_kh,rn_kc,ln_stopack_diags,skeb_filter_pass,rn_skeb_stdev,rn_skeb_tau,&
+      nn_dconv,rn_beta_num, rn_beta_con
+
+      ! Default values
+      rn_sppt_bound = 3._wp
+      ln_stopack = .false.
+      ln_sppt_taumap = .false.
+      rn_sppt_tau =  86400._wp * 5._wp
+      sppt_filter_pass = 30
+      nn_vwei = 0
+      ln_distcoast = .false.
+      ln_sppt_glocon = .false.
+      rn_distcoast = 10._wp
+      ln_stopack_restart = .true.
+      nn_stopack_seed = (/1,2,3,narea/)
+      nn_rndm_freq = 1
+      nn_sppt_step_bound = 2
+      nn_deftau = 1
+      ln_sppt_traxad       = .false.  ! Switch for x advection
+      ln_sppt_trayad       = .false.  ! Switch for y advection
+      ln_sppt_trazad       = .false.  ! Switch for z advection
+      ln_sppt_trasad       = .false.  ! Switch for z advection  (s- case)
+      ln_sppt_traldf       = .false.  ! Switch for lateral diffusion
+      ln_sppt_trazdf       = .false.  ! Switch for vertical diffusion
+      ln_sppt_trazdfp      = .false.  ! Switch for pure vertical diffusion
+      ln_sppt_traevd       = .false.  ! Switch for enhanced vertical diffusion
+      ln_sppt_trabbc       = .false.  ! Switch for bottom boundary condition
+      ln_sppt_trabbl       = .false.  ! Switch for bottom boundary layer
+      ln_sppt_tranpc       = .false.  ! Switch for non-penetrative convection
+      ln_sppt_tradmp       = .false.  ! Switch for tracer damping
+      ln_sppt_traqsr       = .false.  ! Switch for solar radiation
+      ln_sppt_transr       = .false.  ! Switch for non-solar radiation / freshwater flux
+      ln_sppt_traatf       = .false.  ! Switch for Asselin time-filter
+
+      ln_sppt_dynhpg       = .false.  ! Switch for hydrost. press. grad.
+      ln_sppt_dynspg       = .false.  ! Switch for surface  press. grad.
+      ln_sppt_dynkeg       = .false.  ! Switch for horiz. advcetion
+      ln_sppt_dynrvo       = .false.  ! Switch for Relative vorticity
+      ln_sppt_dynpvo       = .false.  ! Switch for planetary vortic.
+      ln_sppt_dynzad       = .false.  ! Switch for vertical advection
+      ln_sppt_dynldf       = .false.  ! Switch for lateral viscosity
+      ln_sppt_dynzdf       = .false.  ! Switch for vertical viscosity
+      ln_sppt_dynbfr       = .false.  ! Switch for bottom friction
+      ln_sppt_dynatf       = .false.  ! Switch for Asselin filter
+      ln_sppt_dyntau       = .false.  ! Switch for wind stress
+
+      ln_sppt_icehdf  = .false.
+      ln_sppt_icelat  = .false.
+      ln_sppt_icezdf  = .false.
+
+      nn_spp_icestr  = 0
+      nn_spp_icealb  = 0
+      rn_icestr_sd = 0.30_wp
+      rn_icealb_sd = 0.30_wp
+
+      nn_spp_bfr  =0
+      nn_spp_dqdt =0
+      nn_spp_arnf =0
+      nn_spp_aevd =0
+      nn_spp_geot =0
+      nn_spp_dedt =0
+      nn_spp_avt  =0
+      nn_spp_avm  =0
+      nn_spp_qsi0 =0
+      nn_spp_relw =0
+      nn_spp_tkelc =0
+      nn_spp_tkedf =0
+      nn_spp_tkeds =0
+      nn_spp_tkebb =0
+      nn_spp_tkefr =0
+
+      ln_skeb = .false.
+      ln_stopack_diags = .false.
+      rn_skeb_stdev = 1.0_wp
+      skeb_filter_pass = 50
+      rn_skeb_tau      = 50
+
+#ifdef NEMO_V34
+      REWIND( numnam )            
+      READ  ( numnam, namstopack )
+#else
+      REWIND( numnam_ref ) 
+      READ  ( numnam_ref, namstopack, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namstopack in reference namelist', lwp )
+
+      REWIND( numnam_cfg )  
+      READ  ( numnam_cfg, namstopack, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namstopack in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namstopack )
+#endif
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'init_stopack : Stochastic physics package'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namstopack: '
+         WRITE(numout,*)
+         WRITE(numout,*) '       Switch on STOPACK                                     ln_stopack : ',ln_stopack
+         WRITE(numout,*) '       Verbose diagnostics for STOPACK                  ln_stopack_diags:', ln_stopack_diags
+         WRITE(numout,*)
+         WRITE(numout,*) '  ***  Random generation'
+         WRITE(numout,*) '       Read from restart file previous perturbation  ln_stopack_restart :', ln_stopack_restart
+         WRITE(numout,*) '       Frequency of calls of the SPPT perturbation field   nn_rndm_freq :', nn_rndm_freq
+         WRITE(numout,*) '       Reproducibility of random generation             ln_stopack_repr :', ln_stopack_repr
+         WRITE(numout,*) '       Seed for random number generator (no restart)    nn_stopack_seed :', nn_stopack_seed(:)
+         WRITE(numout,*)
+         WRITE(numout,*) '  ***  SPPT scheme '
+         WRITE(numout,*) '       SPPT step (0: beginning; 1: before ZVD; 2: after ZVD)  sppt_step : ',sppt_step
+         WRITE(numout,*) '       Use map of decorr. time scale                     ln_sppt_taumap :', ln_sppt_taumap
+         WRITE(numout,*) '       If ln_sppt_taumap FALSE, use this constant (in days) rn_sppt_tau :', rn_sppt_tau
+         WRITE(numout,*) '       Number of filter passes to correlate in space   sppt_filter_pass :', sppt_filter_pass
+         WRITE(numout,*) '       Standard deviation of the white noise              rn_sppt_stdev :', rn_sppt_stdev
+         WRITE(numout,*) '       Apply global conservation constraints             ln_sppt_glocon :', ln_sppt_glocon
+         WRITE(numout,*)
+         WRITE(numout,*) '       Perturbation on tracers:'
+         WRITE(numout,*) '       Switch for x advection                            ln_sppt_traxad :', ln_sppt_traxad
+         WRITE(numout,*) '       Switch for y advection                            ln_sppt_trayad :', ln_sppt_trayad
+         WRITE(numout,*) '       Switch for z advection                            ln_sppt_trazad :', ln_sppt_trazad
+         WRITE(numout,*) '       Switch for z advection  (s- case)                 ln_sppt_trasad :', ln_sppt_trasad
+         WRITE(numout,*) '       Switch for lateral diffusion                      ln_sppt_traldf :', ln_sppt_traldf
+         WRITE(numout,*) '       Switch for vertical diffusion                     ln_sppt_trazdf :', ln_sppt_trazdf
+         WRITE(numout,*) '       Switch for pure vertical diffusion               ln_sppt_trazdfp :', ln_sppt_trazdfp
+         WRITE(numout,*) '       Switch for enhanced vertical diffusion            ln_sppt_traevd :', ln_sppt_traevd
+         WRITE(numout,*) '       Switch for bottom boundary condition              ln_sppt_trabbc :', ln_sppt_trabbc
+         WRITE(numout,*) '       Switch for bottom boundary layer                  ln_sppt_trabbl :', ln_sppt_trabbl
+         WRITE(numout,*) '       Switch for non-penetrative convection             ln_sppt_tranpc :', ln_sppt_tranpc
+         WRITE(numout,*) '       Switch for tracer damping                         ln_sppt_tradmp :', ln_sppt_tradmp
+         WRITE(numout,*) '       Switch for solar radiation                        ln_sppt_traqsr :', ln_sppt_traqsr
+         WRITE(numout,*) '       Switch for non-solar rad. / freshwater flx flux   ln_sppt_transr :', ln_sppt_transr
+         WRITE(numout,*) '       Switch for Asselin time-filter                    ln_sppt_traatf :', ln_sppt_traatf
+         WRITE(numout,*)
+         WRITE(numout,*) '       Perturbation on dynamics:'
+         WRITE(numout,*) '       Switch for horiz advection                        ln_sppt_dynkeg :', ln_sppt_dynkeg
+         WRITE(numout,*) '       Switch for z advection                            ln_sppt_dynzad :', ln_sppt_dynzad
+         WRITE(numout,*) '       Switch for wind stress                            ln_sppt_dyntau :', ln_sppt_dyntau
+         WRITE(numout,*) '       Switch for lateral diffusion                      ln_sppt_dynldf :', ln_sppt_dynldf
+         WRITE(numout,*) '       Switch for vertical diffusion                     ln_sppt_dynzdf :', ln_sppt_dynzdf
+         WRITE(numout,*) '       Switch for planet. vorticity                      ln_sppt_dynpvo :', ln_sppt_dynpvo
+         WRITE(numout,*) '       Switch for relative vorticity                     ln_sppt_dynrvo :', ln_sppt_dynrvo
+         WRITE(numout,*) '       Switch for hydrost. press. grad.                  ln_sppt_dynhpg :', ln_sppt_dynhpg
+         WRITE(numout,*) '       Switch for surface  press. grad.                  ln_sppt_dynspg :', ln_sppt_dynspg
+         WRITE(numout,*) '       Switch for bottom friction                        ln_sppt_dynbfr :', ln_sppt_dynbfr
+         WRITE(numout,*) '       Switch for Asselin time-filter                    ln_sppt_dynatf :', ln_sppt_dynatf
+         WRITE(numout,*)
+         WRITE(numout,*) '       Perturbation on sea-ice:'
+         WRITE(numout,*) '       Switch for sea-ice diffusivity                    ln_sppt_icehdf :', ln_sppt_icehdf
+         WRITE(numout,*) '       Switch for sea-ice lateral accretion              ln_sppt_icelat :', ln_sppt_icelat
+         WRITE(numout,*) '       Switch for sea-ice vertical thermodyn.            ln_sppt_icezdf :', ln_sppt_icezdf
+         WRITE(numout,*)
+         WRITE(numout,*) '       Bound for gaussian random numbers                  rn_sppt_bound :', rn_sppt_bound
+         WRITE(numout,*) '       Bound Step (0: nobound; 1: Gaussian; 2: Pert) nn_sppt_step_bound :', nn_sppt_step_bound
+         WRITE(numout,*) '       Definition of Tau (1: days; 2: timesteps)              nn_deftau :', nn_deftau
+         WRITE(numout,*) '       Smoothing of perturbation close to coast            ln_distcoast :', ln_distcoast
+         WRITE(numout,*) '       Spatial scale of the smoothing near coasts (m)      rn_distcoast :', rn_distcoast
+         WRITE(numout,*) '       Type of vertical weight:                                 nn_vwei :', nn_vwei
+         WRITE(numout,*) '                            0 : No weight '
+         WRITE(numout,*) '                            1 : First/Last level smoothing '
+         WRITE(numout,*) '                            2 : Top/Bottom smoothing'
+         WRITE(numout,*) '                            3 : Bottom smoothing'
+         WRITE(numout,*)
+         WRITE(numout,*) '  ***  SPP schemes (0: off ; 1: normal; 2:lognormal, mean as unpert.; 3: lognormal, median as unpert.'
+         WRITE(numout,*) '                   (the meaning of standard dev. depends on the distribution and parametr.)'
+         WRITE(numout,*)
+         WRITE(numout,*) '       Number of passes for spatial filter (AR1 field)     spp_filter_pass:', spp_filter_pass
+         WRITE(numout,*) '       Standard deviation of random generator (AR1 field)  rn_spp_stdev :', rn_spp_stdev      
+         WRITE(numout,*) '       Decorr. time scale                     (AR1 field)  rn_spp_tau   :', rn_spp_tau
+         WRITE(numout,*)
+         WRITE(numout,*) '       SPP for bottom friction coeff                       nn_spp_bfr   :', nn_spp_bfr  
+         WRITE(numout,*) '                                            STDEV          rn_bfr_sd    :', rn_bfr_sd   
+         WRITE(numout,*) '       SPP for SST relaxation  coeff                       nn_spp_dqdt  :', nn_spp_dqdt 
+         WRITE(numout,*) '                                            STDEV          rn_dqdt_sd   :', rn_dqdt_sd   
+         WRITE(numout,*) '       SPP for SSS relaxation  coeff                       nn_spp_dedt  :', nn_spp_dedt 
+         WRITE(numout,*) '                                            STDEV          rn_dedt_sd   :', rn_dedt_sd   
+         WRITE(numout,*) '       SPP for vertical tra mixing coeff (only TKE, GLS)   nn_spp_avt   :', nn_spp_avt  
+         WRITE(numout,*) '                                            STDEV          rn_avt_sd    :', rn_avt_sd   
+         WRITE(numout,*) '       SPP for vertical dyn mixing coeff (only TKE, GLS)   nn_spp_avm   :', nn_spp_avm  
+         WRITE(numout,*) '                                            STDEV          rn_avm_sd    :', rn_avm_sd   
+         WRITE(numout,*) '       SPP for solar penetration scheme  (only RGB)        nn_spp_qsi0  :', nn_spp_qsi0
+         WRITE(numout,*) '                                            STDEV          rn_qsi0_sd   :', rn_qsi0_sd   
+         WRITE(numout,*) '       SPP for horiz. diffusivity  U                       nn_spp_ahtu  :', nn_spp_ahtu
+         WRITE(numout,*) '                                            STDEV          rn_ahtu_sd   :', rn_ahtu_sd   
+         WRITE(numout,*) '       SPP for horiz. diffusivity  V                       nn_spp_ahtv  :', nn_spp_ahtv
+         WRITE(numout,*) '                                            STDEV          rn_ahtv_sd   :', rn_ahtv_sd   
+         WRITE(numout,*) '       SPP for horiz. diffusivity  W                       nn_spp_ahtw  :', nn_spp_ahtw
+         WRITE(numout,*) '                                            STDEV          rn_ahtw_sd   :', rn_ahtw_sd   
+         WRITE(numout,*) '       SPP for horiz. diffusivity  T                       nn_spp_ahtt  :', nn_spp_ahtt
+         WRITE(numout,*) '                                            STDEV          rn_ahtt_sd   :', rn_ahtt_sd   
+         WRITE(numout,*) '       SPP for horiz. viscosity (1/3)                      nn_spp_ahm1  :', nn_spp_ahm1
+         WRITE(numout,*) '                                            STDEV          rn_ahm1_sd   :', rn_ahm1_sd   
+         WRITE(numout,*) '       SPP for horiz. viscosity (2/4)                      nn_spp_ahm2  :', nn_spp_ahm2
+         WRITE(numout,*) '                                            STDEV          rn_ahm2_sd   :', rn_ahm2_sd   
+         WRITE(numout,*) '       SPP for relative wind factor                        nn_spp_relw  :', nn_spp_relw
+         WRITE(numout,*) '       (use 4, 5, 6 for nn_spp_relw to have options 1, 2, 3 with limits bounded to [0,1]'
+         WRITE(numout,*) '                                            STDEV          rn_relw_sd   :', rn_relw_sd   
+         WRITE(numout,*) '       SPP for mixing close to river mouth                 nn_spp_arnf  :', nn_spp_arnf
+         WRITE(numout,*) '                                            STDEV          rn_arnf_sd   :', rn_arnf_sd   
+         WRITE(numout,*) '       SPP for geothermal heating                          nn_spp_geot  :', nn_spp_geot
+         WRITE(numout,*) '                                            STDEV          rn_geot_sd   :', rn_geot_sd   
+         WRITE(numout,*) '       SPP for enhanced vertical diffusion                 nn_spp_aevd  :', nn_spp_aevd
+         WRITE(numout,*) '                                            STDEV          rn_aevd_sd   :', rn_aevd_sd   
+         WRITE(numout,*) '       SPP for TKE rn_lc    Langmuir cell coefficient      nn_spp_tkelc :', nn_spp_tkelc
+         WRITE(numout,*) '                                            STDEV          rn_tkelc_sd  :', rn_tkelc_sd   
+         WRITE(numout,*) '       SPP for TKE rn_ediff Eddy diff. coefficient         nn_spp_tkedf :', nn_spp_tkedf
+         WRITE(numout,*) '                                            STDEV          rn_tkedf_sd  :', rn_tkedf_sd   
+         WRITE(numout,*) '       SPP for TKE rn_ediss Kolmogoroff dissipation coeff. nn_spp_tkeds :', nn_spp_tkeds
+         WRITE(numout,*) '                                            STDEV          rn_tkeds_sd  :', rn_tkeds_sd   
+         WRITE(numout,*) '       SPP for TKE rn_ebb   Surface input of tke           nn_spp_tkebb :', nn_spp_tkebb
+         WRITE(numout,*) '                                            STDEV          rn_tkebb_sd  :', rn_tkebb_sd   
+         WRITE(numout,*) '       SPP for TKE rn_efr   Fraction of srf TKE below ML   nn_spp_tkefr :', nn_spp_tkefr
+         WRITE(numout,*) '                                            STDEV          rn_tkefr_sd  :', rn_tkefr_sd   
+         WRITE(numout,*) '       SPP for BBL U  diffusivity                          nn_spp_ahubbl:', nn_spp_ahubbl
+         WRITE(numout,*) '                                            STDEV          rn_ahubbl_sd :', rn_ahubbl_sd
+         WRITE(numout,*) '       SPP for BBL V  diffusivity                          nn_spp_ahvbbl:', nn_spp_ahvbbl
+         WRITE(numout,*) '                                            STDEV          rn_ahvbbl_sd :', rn_ahvbbl_sd
+         WRITE(numout,*)
+         WRITE(numout,*) '  ***  SPP schemes for sea-ice '
+         WRITE(numout,*) '       Albedo                                              nn_spp_icealb:', nn_spp_icealb
+         WRITE(numout,*) '            St. dev. for ice albedo                        rn_icealb_sd :', rn_icealb_sd
+         WRITE(numout,*) '       Ice Strength                                        nn_spp_icestr:', nn_spp_icestr
+         WRITE(numout,*) '            St. dev. for ice strength                      rn_icestr_sd :', rn_icestr_sd
+         WRITE(numout,*)
+         WRITE(numout,*) ' SKEB Perturbation scheme '
+         WRITE(numout,*) '       SKEB switch                                         ln_skeb      :', ln_skeb    
+         WRITE(numout,*) '       SKEB ratio of backscattered energy                  rn_skeb      :', rn_skeb    
+         WRITE(numout,*) '       Frequency update for dissipation mask               nn_dcom_freq :', nn_dcom_freq
+         WRITE(numout,*) '       Numerical dissipation factor (resolut. dependent)   rn_kh        :', rn_kh
+         WRITE(numout,*) '       Number of passes for spatial filter (AR1 field)     skeb_filter_pass:', skeb_filter_pass
+         WRITE(numout,*) '       Standard deviation of random generator (AR1 field)  rn_skeb_stdev:', rn_skeb_stdev      
+         WRITE(numout,*) '       Decorr. time scale                     (AR1 field)  rn_skeb_tau  :', rn_skeb_tau
+         WRITE(numout,*) '       Option of convection energy dissipation             nn_dconv     :', nn_dconv
+         WRITE(numout,*) '       Convection dissipation factor (resolut. dependent)  rn_kc        :', rn_kc
+         WRITE(numout,*) '       Multiplier for numerical dissipation                rn_beta_num  :', rn_beta_num
+         WRITE(numout,*) '       Multiplier for convection dissipation               rn_beta_con  :', rn_beta_con
+
+         WRITE(numout,*)
+      ENDIF
+
+      IF( .NOT. ln_stopack ) THEN
+         IF(lwp) WRITE(numout,*) '     STOPACK is switched off'
+         RETURN
+      ENDIF
+
+      IF( MOD( nitend - nit000 + 1, nn_rndm_freq) /= 0 .OR.   &
+          ( MOD( nstock, nn_rndm_freq) /= 0 .AND. nstock .GT. 0 ) ) THEN
+         WRITE(ctmp1,*) 'experiment length (', nitend - nit000 + 1, ') or nstock (', nstock,   &
+            &           ' is NOT a multiple of nn_rndm_freq (', nn_rndm_freq, ')'
+         CALL ctl_stop( ctmp1, 'Impossible to properly setup STOPACK restart' )
+      ENDIF
+
+      nn_sppt_tra = COUNT( (/ln_sppt_traxad, ln_sppt_trayad, ln_sppt_trazad, &
+                             ln_sppt_trasad, ln_sppt_traldf, ln_sppt_trazdf, &
+                            ln_sppt_trazdfp, ln_sppt_traevd, ln_sppt_trabbc, &
+                             ln_sppt_trabbl, ln_sppt_tranpc, ln_sppt_tradmp, &
+                             ln_sppt_traqsr, ln_sppt_transr, ln_sppt_traatf/) )
+      nn_sppt_dyn = COUNT( (/ln_sppt_dynhpg, ln_sppt_dynspg, ln_sppt_dynkeg, &
+                             ln_sppt_dynrvo, ln_sppt_dynpvo, ln_sppt_dynzad, &
+                             ln_sppt_dynldf, ln_sppt_dynzdf, ln_sppt_dynbfr, &
+                             ln_sppt_dynatf, ln_sppt_dyntau/) )
+      nn_sppt_ice = COUNT( (/ln_sppt_icehdf, ln_sppt_icelat, ln_sppt_icezdf/) )
+
+      ln_sppt_tra = ( nn_sppt_tra > 0 )
+      ln_sppt_dyn = ( nn_sppt_dyn > 0 )
+      ln_sppt_ice = ( nn_sppt_ice > 0 )
+
+      nn_spp = nn_spp_bfr+nn_spp_dqdt+nn_spp_dedt+nn_spp_avt+nn_spp_avm+nn_spp_qsi0+&
+      & nn_spp_ahtu+nn_spp_ahtv+nn_spp_ahtw+nn_spp_ahtt+nn_spp_relw+nn_spp_arnf+nn_spp_geot+nn_spp_aevd+&
+      & nn_spp_tkelc+nn_spp_tkedf+nn_spp_tkeds+nn_spp_tkebb+nn_spp_tkefr+nn_spp_ahubbl+nn_spp_ahvbbl+&
+      & nn_spp_ahm1+nn_spp_ahm2
+
+#ifdef NEMO_V34
+
+#ifndef key_trdtra
+      IF( ln_sppt_tra ) &
+         & CALL ctl_stop( ' SPPT on tracers on .AND. .NOT. key_trdtra ', &
+         &                ' SPPT cannot work without tracer tendency computation')
+#endif
+#ifndef key_trddyn
+      IF( ln_sppt_dyn ) &
+         & CALL ctl_stop( ' SPPT on momentum on .AND. .NOT. key_trddyn ', &
+         &                ' SPPT cannot work without dynamic tendency computation')
+#endif
+
+#else
+      IF( ln_sppt_tra .AND. .NOT.  ln_tra_trd ) &
+         & CALL ctl_stop( ' SPPT on tracers on .AND. .NOT. ln_tra_trd ', &
+         &                ' SPPT cannot work without tracer tendency computation')
+
+      IF( ln_sppt_dyn .AND. .NOT.  ln_dyn_trd ) &
+         & CALL ctl_stop( ' SPPT on momentum on .AND. .NOT. ln_dyn_trd ', &
+         &                ' SPPT cannot work without dynamic tendency computation')
+#endif
+
+      IF( ln_sppt_glocon .AND. lk_vvl ) &
+         & CALL ctl_stop( ' ln_sppt_glocal .AND. lk_vvl ', &
+         &                ' SPPT conservation not coded yet for VVL')
+
+      IF( nn_deftau .NE. 1 .AND. nn_deftau .NE. 2 ) &
+         & CALL ctl_stop( ' nn_deftau must be 1 or 2 ')
+
+      nn_spp_aht = nn_spp_ahtu + nn_spp_ahtv + nn_spp_ahtw + nn_spp_ahtt
+      IF(nn_spp_aht .GT. 0) THEN
+#if defined key_traldf_c3d
+        IF(lwp) WRITE(numout,*) 'SPP : diffusivity perturbation with 3D coefficients'
+#elif defined key_traldf_c2d
+        IF(lwp) WRITE(numout,*) 'SPP : diffusivity perturbation with 2D coefficients'
+#else
+        CALL ctl_stop( 'SPP : diffusivity perturbation requires key_traldf_c3d or key_traldf_c2d')
+#endif
+      ENDIF
+
+      IF(nn_spp_ahm1+nn_spp_ahm2 .GT. 0) THEN
+#if defined key_dynldf_c3d
+        IF(lwp) WRITE(numout,*) 'SPP : viscosity perturbation with 3D coefficients'
+#elif defined key_dynldf_c2d
+        IF(lwp) WRITE(numout,*) 'SPP : viscosity perturbation with 2D coefficients'
+#else
+        CALL ctl_stop( 'SPP : viscosity perturbation requires key_dynldf_c3d or key_dynldf_c2d')
+#endif
+      ENDIF
+
+      ln_skeb_own_gauss = .FALSE.
+      IF( ln_skeb ) THEN
+        IF( skeb_filter_pass > 0 .AND. skeb_filter_pass .NE. sppt_filter_pass ) ln_skeb_own_gauss = .TRUE.
+        IF( rn_skeb_tau > 0._wp .AND. rn_skeb_tau .NE. rn_sppt_tau ) ln_skeb_own_gauss = .TRUE.
+      ENDIF
+
+      ln_spp_own_gauss = .FALSE.
+      IF( nn_spp > 0 ) THEN
+        IF( spp_filter_pass > 0 .AND. spp_filter_pass .NE. sppt_filter_pass ) ln_spp_own_gauss = .TRUE.
+        IF( rn_spp_tau > 0._wp .AND. rn_spp_tau .NE. rn_sppt_tau ) ln_spp_own_gauss = .TRUE.
+      ENDIF
+
+      IF(lwp) THEN
+         WRITE(numout,*) '    SPPT on tracers     : ',ln_sppt_tra
+         WRITE(numout,*) '    SPPT on dynamics    : ',ln_sppt_dyn
+         WRITE(numout,*) '    SPPT on sea-ice     : ',ln_sppt_ice
+         WRITE(numout,*) '    SPP perturbed params: ',nn_spp
+         WRITE(numout,*) '    SKEB scheme         : ',ln_skeb
+         WRITE(numout,*) '    SPP with own scales : ',ln_spp_own_gauss
+         WRITE(numout,*) '    SKEB with own scales: ',ln_skeb_own_gauss
+      ENDIF
+
+      IF( sppt_tra_alloc() /= 0) CALL ctl_stop( 'STOP', 'sppt_alloc: unable to allocate arrays' )
+
+      spptt = 0._wp  ; sppts = 0._wp
+      spptu = 0._wp  ; spptv = 0._wp
+
+      ! Find filter attenuation factor
+   
+      flt_fac = sppt_fltfac( sppt_filter_pass )
+      rdt_sppt = nn_rndm_freq * rn_rdt
+   
+      IF( ln_sppt_taumap ) THEN
+         CALL iom_open ( 'sppt_tau_map', inum )
+         CALL iom_get  ( inum, jpdom_data, 'tau', sppt_tau )
+         CALL iom_close( inum )
+      ELSE
+         sppt_tau(:,:) = rn_sppt_tau
+      ENDIF
+
+      IF( nn_deftau .EQ. 1 ) THEN
+         ! Decorrelation time-scale defined in days
+         sppt_tau(:,:) = exp( -rdt_sppt / (sppt_tau(:,:)*86400._wp) )
+      ELSE
+         ! Decorrelation time-scale defined in timesteps
+         sppt_tau(:,:) = exp( -REAL( nn_rndm_freq, wp) / sppt_tau(:,:) )
+      ENDIF
+
+      IF( ln_distcoast ) THEN
+        CALL iom_open('dist.coast', inum )
+        CALL iom_get(inum,jpdom_data,'Tcoast',zdc)
+        CALL iom_close( inum )
+        zdc = 1._wp - exp(-zdc/rn_distcoast)
+        CALL lbc_lnk( zdc , 'T', 1._wp)
+      ELSE
+        zdc = 1._wp
+      ENDIF
+
+      ! Initialize
+      sppt_a = sppt_tau
+      sppt_b = sqrt ( 1._wp - sppt_tau*sppt_tau )
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '    **** SPPT SCHEME'
+         WRITE(numout,*) '    Definit. time-scale : ',nn_deftau
+         WRITE(numout,*) '     Decorr. time-scale : ',rn_sppt_tau
+         WRITE(numout,*) '        Mean AR1 a term : ',SUM(sppt_a)/REAL(jpi*jpj)
+         WRITE(numout,*) '        Mean AR1 b term : ',SUM(sppt_b)/REAL(jpi*jpj)
+         WRITE(numout,*)
+      ENDIF
+
+      gauss_b = 0._wp
+      ! Weigths
+      gauss_w(:)    = 1.0_wp 
+      IF( nn_vwei .eq. 1 ) THEN
+        gauss_w(1)    = 0.0_wp
+        gauss_w(2)    = 0.5_wp
+        gauss_w(jpk)  = 0.0_wp
+        gauss_w(jpkm1)= 0.5_wp
+      ENDIF
+
+      IF( ln_spp_own_gauss ) THEN
+        IF( spp_alloc() /= 0) CALL ctl_stop( 'STOP', 'spp_alloc: unable to allocate arrays' )
+        flt_fac_p = sppt_fltfac( spp_filter_pass )
+        spp_tau (:, :) = rn_spp_tau
+        IF( nn_deftau .EQ. 1 ) THEN
+          spp_tau(:,:) = spp_tau(:,:) * 86400._wp
+          spp_tau(:,:) = exp( -rdt_sppt / (spp_tau) )
+        ELSE
+          spp_tau(:,:) = exp( -1._wp / spp_tau )
+        ENDIF
+        spp_a = spp_tau
+        spp_b = sqrt ( 1._wp - spp_tau*spp_tau )
+        gauss_bp = 0._wp
+        IF(lwp) THEN
+          WRITE(numout,*)
+          WRITE(numout,*) '    **** SPP  SCHEME'
+          WRITE(numout,*) '    Definit. time-scale : ',nn_deftau
+          WRITE(numout,*) '     Decorr. time-scale : ',rn_spp_tau
+          WRITE(numout,*) '        Mean AR1 a term : ',SUM(spp_a)/REAL(jpi*jpj)
+          WRITE(numout,*) '        Mean AR1 b term : ',SUM(spp_b)/REAL(jpi*jpj)
+          WRITE(numout,*)
+        ENDIF
+      ENDIF
+
+      IF( ln_skeb_own_gauss ) THEN
+        IF( skeb_alloc() /= 0) CALL ctl_stop( 'STOP', 'skeb_alloc: unable to allocate arrays' )
+        flt_fac_k = sppt_fltfac( skeb_filter_pass )
+        skeb_tau (:, :) = rn_skeb_tau
+        IF( nn_deftau .EQ. 1 ) THEN
+          skeb_tau(:,:) = skeb_tau(:,:) * 86400._wp
+          skeb_tau(:,:) = exp( -rdt_sppt / (skeb_tau) )
+        ELSE
+          skeb_tau(:,:) = exp( -1._wp / skeb_tau )
+        ENDIF
+        skeb_a = skeb_tau
+        skeb_b = sqrt ( 1._wp - skeb_tau*skeb_tau )
+        gauss_bk = 0._wp
+        IF(lwp) THEN
+          WRITE(numout,*)
+          WRITE(numout,*) '    **** SEB  SCHEME'
+          WRITE(numout,*) '    Definit. time-scale : ',nn_deftau
+          WRITE(numout,*) '     Decorr. time-scale : ',rn_skeb_tau
+          WRITE(numout,*) '        Mean AR1 a term : ',SUM(skeb_a)/REAL(jpi*jpj)
+          WRITE(numout,*) '        Mean AR1 b term : ',SUM(skeb_b)/REAL(jpi*jpj)
+          WRITE(numout,*)
+        ENDIF
+      ENDIF
+
+      IF( ln_skeb_own_gauss .OR. ln_spp_own_gauss ) &
+      & ALLOCATE ( g2d_save(jpi,jpj) )
+
+      CALL stopack_rst( nit000, 'READ' )
+
+      IF(lwp .and. ln_stopack_diags) &
+      CALL ctl_opn(numdiag, 'stopack.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+   
+   END SUBROUTINE stopack_init
+   !
+   SUBROUTINE stopack_rst( kt, cdrw )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE stopack_rst ***
+      !!
+      !! ** Purpose :   Read/write from/to restarsts seed and perturbation field
+      !!----------------------------------------------------------------------
+      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
+      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
+      INTEGER :: id1, jseed
+      CHARACTER(LEN=10)   ::   clseed='spsd0_0000'
+      INTEGER(KIND=8)     ::   ziseed(4)           ! RNG seeds in integer type 
+      INTEGER(KIND=8)     ::   ivals(8) 
+      REAL(wp)            ::   zrseed4(4)           ! RNG seeds in integer type
+      REAL(wp)            ::   zrseed2d(jpi,jpj)
+      !
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
+         !                                   ! ---------------
+         IF(lwp) WRITE(numout,*) ' *** stopack-rst ***'
+         IF( ln_rstart ) THEN !* Read the restart file
+            id1 = iom_varid( numror, 'sppt_b', ldstop = .FALSE. )
+            !
+            IF( id1 > 0 ) THEN
+                CALL iom_get( numror, jpdom_autoglo, 'sppt_b', gauss_b )
+                DO jseed = 1 , 4
+                   WRITE(clseed(5:5) ,'(i1.1)') jseed
+                   CALL iom_get( numror, jpdom_autoglo, clseed(1:5) , zrseed2d )
+                   zrseed4(jseed) = zrseed2d(2,2)
+                   ziseed(jseed) = TRANSFER( zrseed4(jseed) , ziseed(jseed) )
+               END DO
+               IF (lwp) THEN
+                  WRITE(numout,*) 'Read ziseed4 = ',zrseed4(:)
+                  WRITE(numout,*) 'Read ziseed  = ',ziseed(:)
+               ENDIF
+               CALL kiss_seed( ziseed(1) , ziseed(2) , ziseed(3) , ziseed(4) )
+            ELSE
+               IF( lwp ) WRITE(numout,*) ' ===>>>> : previous run without sppt_b'
+               IF( ln_stopack_restart ) CALL ctl_stop( ' ln_stopack_restart TRUE :',&
+               & ' variable not found in restart file ')
+               IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of sppt_b to 0.'
+               gauss_b = 0._wp
+               IF ( .NOT. ln_stopack_repr ) THEN
+                  CALL date_and_time(VALUES=ivals)
+                  DO jseed=1,4
+                    nn_stopack_seed(jseed) = nn_stopack_seed(jseed) + INT(ivals(4+jseed))
+                  ENDDO
+               ENDIF
+               DO jseed=1,4
+                 ziseed(jseed) = TRANSFER(nn_stopack_seed(jseed)+narea, ziseed(jseed))
+               ENDDO
+               IF(lwp) WRITE(numout,*) ' ===>>>> Seed, nn_stopack_seed+narea',nn_stopack_seed+narea
+               IF(lwp) WRITE(numout,*) ' ===>>>> Seed, ziseed',ziseed
+               CALL kiss_seed( ziseed(1) , ziseed(2) , ziseed(3) , ziseed(4) )
+            ENDIF
+            IF( ln_spp_own_gauss ) THEN
+              id1 = iom_varid( numror, 'spp_b', ldstop = .FALSE. )
+              IF( id1 > 0 ) THEN
+                 CALL iom_get( numror, jpdom_autoglo, 'spp_b', gauss_bp )
+              ELSE
+                 IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of spp_b to 0.'
+                 gauss_bp = 0._wp
+              ENDIF
+            ENDIF
+            IF( ln_skeb_own_gauss ) THEN
+              id1 = iom_varid( numror, 'skeb_b', ldstop = .FALSE. )
+              IF( id1 > 0 ) THEN
+                 CALL iom_get( numror, jpdom_autoglo, 'skeb_b', gauss_bk )
+              ELSE
+                 IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of skeb_b to 0.'
+                 gauss_bk = 0._wp
+              ENDIF
+            ENDIF
+         ELSE
+            gauss_b = 0._wp
+            gauss_bp = 0._wp
+            gauss_bk = 0._wp
+            IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of STOPACK to 0.'
+            IF ( .NOT. ln_stopack_repr ) THEN
+                CALL date_and_time(VALUES=ivals)
+                DO jseed=1,4
+                   nn_stopack_seed(jseed) = nn_stopack_seed(jseed) + INT(ivals(4+jseed))
+                ENDDO
+            ENDIF
+            DO jseed=1,4
+               ziseed(jseed) = TRANSFER(nn_stopack_seed(jseed)+narea, ziseed(jseed))
+            ENDDO
+            IF(lwp) WRITE(numout,*) ' ===>>>> Seed, nn_stopack_seed+narea',nn_stopack_seed+narea
+            IF(lwp) WRITE(numout,*) ' ===>>>> Seed, ziseed',ziseed
+            CALL kiss_seed( ziseed(1) , ziseed(2) , ziseed(3) , ziseed(4) )
+         ENDIF
+         !
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+         !                                   ! -------------------
+         IF(lwp) WRITE(numout,*) '---- stopack-rst ----'
+         CALL iom_rstput( kt, nitrst, numrow, 'sppt_b', gauss_b )
+         IF(ln_spp_own_gauss) CALL iom_rstput( kt, nitrst, numrow, 'spp_b', gauss_bp )
+         IF(ln_skeb_own_gauss) CALL iom_rstput( kt, nitrst, numrow, 'skeb_b', gauss_bk )
+         CALL kiss_state( ziseed(1) , ziseed(2) , ziseed(3) , ziseed(4) )
+         DO jseed=1,4
+            zrseed4(jseed) = TRANSFER(ziseed(jseed), zrseed4(jseed))
+         ENDDO
+         IF (lwp) THEN
+            WRITE(numout,*) 'Write ziseed4 = ',zrseed4(:)
+            WRITE(numout,*) 'Write ziseed  = ',ziseed(:)
+         ENDIF
+         DO jseed = 1 , 4
+            WRITE(clseed(5:5) ,'(i1.1)') jseed
+            zrseed2d(:,:) = zrseed4(jseed)
+            CALL iom_rstput( kt, nitrst, numrow, clseed(1:5) , zrseed2d )
+         END DO
+         !
+      ENDIF
+      !
+   END SUBROUTINE stopack_rst
+   !
+   INTEGER FUNCTION sppt_tra_alloc()
+      !!---------------------------------------------------------------------
+      !!                  ***  FUNCTION sppt_tra_alloc  ***
+      !!---------------------------------------------------------------------
+      !
+      ALLOCATE( spptt(jpi,jpj,jpk) , sppts(jpi,jpj,jpk) , gauss_n(jpi,jpj,jpk) ,& 
+      gauss_nu(jpi,jpj,jpk) , gauss_nv(jpi,jpj,jpk) , & 
+      spptu(jpi,jpj,jpk) , spptv(jpi,jpj,jpk) , gauss_n_2d(jpi,jpj) ,&
+      gauss_b (jpi,jpj), sppt_tau(jpi,jpj), sppt_a(jpi,jpj), sppt_b(jpi,jpj), gauss_w(jpk),&
+      zdc(jpi,jpj,jpk), STAT= sppt_tra_alloc )
+      !
+      IF( lk_mpp             )   CALL mpp_sum ( sppt_tra_alloc )
+      IF( sppt_tra_alloc /= 0)   CALL ctl_warn('sppt_tra_alloc: failed to allocate arrays')
+
+      IF(nn_spp_bfr>0) THEN
+        ALLOCATE(coeff_bfr(jpi,jpj), STAT= sppt_tra_alloc )
+        IF( lk_mpp             )   CALL mpp_sum ( sppt_tra_alloc )
+        IF( sppt_tra_alloc /= 0)   CALL ctl_warn('sppt_tra_alloc: failed to allocate arrays')
+        coeff_bfr = 0._wp
+      ENDIF
+      gauss_b = 0._wp
+      gauss_n = 0._wp
+      gauss_nu= 0._wp
+      gauss_nv= 0._wp
+      gauss_n_2d = 0._wp
+
+   END FUNCTION sppt_tra_alloc
+
+   INTEGER FUNCTION skeb_alloc()
+      !!---------------------------------------------------------------------
+      !!                  ***  FUNCTION skeb_alloc  ***
+      !!---------------------------------------------------------------------
+      !
+      ALLOCATE( gauss_n_2d_k(jpi,jpj) , gauss_bk (jpi,jpj),skeb_tau(jpi,jpj),&
+      skeb_a(jpi,jpj), skeb_b(jpi,jpj), STAT= skeb_alloc )
+      !
+      IF( lk_mpp             )   CALL mpp_sum ( skeb_alloc )
+      IF( skeb_alloc /= 0)   CALL ctl_warn('sppt_tra_alloc: failed to allocate arrays')
+      gauss_n_2d_k = 0._wp
+      gauss_bk = 0._wp
+
+   END FUNCTION skeb_alloc
+
+   INTEGER FUNCTION spp_alloc()
+      !!---------------------------------------------------------------------
+      !!                  ***  FUNCTION spp_alloc  ***
+      !!---------------------------------------------------------------------
+      !
+      ALLOCATE( gauss_n_2d_p(jpi,jpj), gauss_bp (jpi,jpj),spp_tau(jpi,jpj), &
+      spp_a(jpi,jpj), spp_b(jpi,jpj), STAT= spp_alloc )
+      !
+      IF( lk_mpp             )   CALL mpp_sum ( spp_alloc )
+      IF( spp_alloc /= 0)   CALL ctl_warn('sppt_tra_alloc: failed to allocate arrays')
+      gauss_n_2d_p = 0._wp
+      gauss_bp = 0._wp
+
+   END FUNCTION spp_alloc
+
+   SUBROUTINE sppt_flt( psto )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE sppt_flt  ***
+      !!
+      !! ** Purpose :   apply horizontal Laplacian filter to input array
+      !!                Adapted from STO package
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psto
+      REAL(wp), DIMENSION(jpi,jpj)                ::   pstm
+      !!
+      INTEGER  :: ji, jj
+
+      pstm = psto
+      DO jj = 2, jpj-1
+         DO ji = 2, jpi-1
+            psto(ji,jj) = 0.5_wp * pstm(ji,jj) + 0.125_wp * &
+                              &  ( pstm(ji-1,jj) + pstm(ji+1,jj) +  &
+                              &    pstm(ji,jj-1) + pstm(ji,jj+1) )
+         END DO
+      END DO
+      CALL lbc_lnk( psto , 'T', 1._wp )
+
+   END SUBROUTINE sppt_flt
+
+   SUBROUTINE sppt_rand2d( psto )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE sppt_rand2d ***
+      !!
+      !! ** Purpose :   fill input array with white Gaussian noise
+      !!                Adapted from STO package
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(out)           ::   psto
+      !!
+      INTEGER  :: ji, jj
+      REAL(KIND=8) :: gran   ! Gaussian random number (forced KIND=8 as in kiss_gaussian)
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            CALL kiss_gaussian( gran )
+            psto(ji,jj) = gran*rn_sppt_stdev
+         END DO
+      END DO
+
+   END SUBROUTINE sppt_rand2d
+
+   FUNCTION sppt_fltfac( kpasses )
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION sppt_fltfac  ***
+      !!
+      !! ** Purpose :   compute factor to restore standard deviation
+      !!                as a function of the number of passes
+      !!                of the Laplacian filter
+      !!                From STO package
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) :: kpasses
+      REAL(wp) :: sppt_fltfac
+      !!
+      INTEGER :: jpasses, ji, jj, jflti, jfltj
+      INTEGER, DIMENSION(-1:1,-1:1) :: pflt0
+      ! 16-b reals to avoid overflows
+      INTEGER,  PARAMETER ::   qp = SELECTED_REAL_KIND(33,4931)
+      REAL(qp), DIMENSION(:,:), ALLOCATABLE :: pfltb
+      REAL(qp), DIMENSION(:,:), ALLOCATABLE :: pflta
+      REAL(qp) :: ratio
+      REAL(qp) :: aux0, aux1
+      pflt0(-1,-1) = 0 ; pflt0(-1,0) = 1 ; pflt0(-1,1) = 0
+      pflt0( 0,-1) = 1 ; pflt0( 0,0) = 4 ; pflt0( 0,1) = 1
+      pflt0( 1,-1) = 0 ; pflt0( 1,0) = 1 ; pflt0( 1,1) = 0
+      ALLOCATE(pfltb(-kpasses-1:kpasses+1,-kpasses-1:kpasses+1))
+      ALLOCATE(pflta(-kpasses-1:kpasses+1,-kpasses-1:kpasses+1))
+      pfltb(:,:) = 0
+      pfltb(0,0) = 1
+      DO jpasses = 1, kpasses
+        pflta(:,:) = 0._qp
+        DO jflti= -1, 1
+        DO jfltj= -1, 1
+          DO ji= -kpasses, kpasses
+          DO jj= -kpasses, kpasses
+            pflta(ji,jj) = pflta(ji,jj) + pfltb(ji+jflti,jj+jfltj) * pflt0(jflti,jfltj)
+          ENDDO
+          ENDDO
+        ENDDO
+        ENDDO
+        pfltb(:,:) = pflta(:,:)
+      ENDDO
+      ratio = SUM(pfltb(:,:))
+      aux0  = SUM(pfltb(:,:)*pfltb(:,:))
+      aux1  = ratio*ratio
+      ratio = aux1 / aux0
+      ratio = SQRT(ratio)
+      DEALLOCATE(pfltb,pflta)
+      sppt_fltfac = REAL(ratio, wp)
+   END FUNCTION sppt_fltfac
+
+   SUBROUTINE skeb_comp( lt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE skeb_comp  ***
+      !!
+      !! ** Purpose :   Computation of energy dissipation terms
+      !!                This is a wrapper to the enrgy terms computation
+      !!                routines
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(IN)  :: lt
+      CALL skeb_dnum ( lt )
+      CALL skeb_dcon ( lt )
+   END SUBROUTINE skeb_comp
+
+   SUBROUTINE bsfcomp( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE bsfcomp  ***
+      !!
+      !! ** Purpose :   Online computation of barotropic streamfunction
+      !!                For diagnostics purposes
+      !!                This routine is not explicitly called anywhere
+      !!                and utilizes low-level MPI routines
+      !!----------------------------------------------------------------------
+      USE MPI
+      !
+      INTEGER, INTENT(IN)  :: kt
+      !
+      REAL(wp) :: dtrpv(jpi,jpj)
+      INTEGER  :: jk,ji,jj
+      !
+#if defined key_mpp_mpi
+      INTEGER, DIMENSION(1) ::   ish
+      INTEGER, DIMENSION(2) ::   ish2
+      INTEGER               ::   ijp,tag,ierr,jp,isend,irecv
+      REAL(wp), DIMENSION(jpj) ::   zwrk    ! mask flux array at V-point
+      REAL(wp), DIMENSION(jpi) ::   zwr2    ! mask flux array at V-point
+      INTEGER :: istatus(MPI_STATUS_SIZE)
+#endif
+      IF(kt .EQ. nit000 ) THEN
+#ifdef NEMO_V34
+        ln_dpsiu = .FALSE.
+        ln_dpsiv = .FALSE.
+#else
+        IF (iom_use('bsft') ) ln_dpsiu = .TRUE.
+        IF (iom_use('bsftv') ) ln_dpsiv = .TRUE.
+#endif
+        IF( ln_dpsiv ) ALLOCATE ( dpsiv(jpi,jpj) )
+        IF( ln_dpsiu ) ALLOCATE ( dpsiu(jpi,jpj) )
+        IF( ln_dpsiv .OR. ln_dpsiu ) THEN
+          jpri = nint ( REAL(nimpp) / REAL(jpi) ) + 1
+          jprj = nint ( REAL(njmpp) / REAL(jpj) ) + 1
+        ENDIF
+      ENDIF
+
+      IF( ln_dpsiv ) THEN
+       dtrpv = 0._wp
+       DO jk = 1,jpkm1
+           dtrpv = dtrpv + fse3v(:,:,jk) * e1v(:,:) * vn(:,:,jk)
+       ENDDO
+       dpsiv(1,:)=0._wp
+       DO ji = 2,jpi
+          dpsiv(ji,:) = dpsiv(ji-1,:) + dtrpv(ji,:)
+       END DO
+      ENDIF
+
+      IF( ln_dpsiu ) THEN
+       dtrpv = 0._wp
+       DO jk = 1,jpkm1
+           dtrpv = dtrpv + fse3u(:,:,jk) * e2u(:,:) * un(:,:,jk)
+       ENDDO
+       dpsiu(:,1)=0._wp
+       DO jj = 2,jpj
+          dpsiu(:,jj) = dpsiu(:,jj-1) + dtrpv(:,jj)
+       END DO
+      ENDIF
+
+#if defined key_mpp_mpi
+      IF ( ln_dpsiv ) THEN
+       DO jp=1,jpni-1
+         IF( jpri == jp ) THEN ! SEND TO EAST 
+          zwrk(1:jpj) = dpsiv(jpi-1,:)
+          tag=2000+narea
+          CALL mpi_isend(zwrk, jpj, mpi_double_precision, noea, tag, mpi_comm_opa, isend,ierr)
+         ELSEIF ( jpri == jp+1 ) THEN ! RECEIVE FROM WEST
+          CALL mpi_irecv(zwrk, jpj, mpi_double_precision, nowe, mpi_any_tag, mpi_comm_opa, irecv,ierr)
+         ENDIF
+         IF(jpri == jp) CALL mpi_wait(isend, istatus, ierr)
+         IF(jpri == jp+1 ) THEN
+          CALL mpi_wait(irecv, istatus, ierr)
+          DO ji=1,jpi
+            dpsiv(ji,:) = dpsiv(ji,:) + zwrk(:)
+          ENDDO
+         ENDIF
+       ENDDO
+      ENDIF
+
+      IF ( ln_dpsiv ) THEN
+       DO jp=1,jpnj-1
+         IF( jprj == jp ) THEN ! SEND TO NORTH
+          zwr2(1:jpi) = dpsiu(:,jpj-1)
+          tag=3000+narea
+          CALL mpi_isend(zwr2, jpi, mpi_double_precision, nono, tag, mpi_comm_opa, isend,ierr)
+         ELSEIF ( jprj == jp+1 ) THEN ! RECEIVE FROM SOUTH
+          CALL mpi_irecv(zwr2, jpi, mpi_double_precision, noso, mpi_any_tag, mpi_comm_opa, irecv,ierr)
+         ENDIF
+         IF(jprj == jp) CALL mpi_wait(isend, istatus, ierr)
+         IF(jprj == jp+1 ) THEN
+          CALL mpi_wait(irecv, istatus, ierr)
+          DO ji=1,jpj
+            dpsiu(:,ji) = dpsiu(:,ji) + zwr2(:)
+          ENDDO
+         ENDIF
+       ENDDO
+      ENDIF
+#endif
+      IF (ln_dpsiu ) THEN
+          CALL lbc_lnk(dpsiu,'T',1._wp)
+#ifdef key_iomput
+          CALL iom_put( "bsft" , dpsiu )
+#endif
+      ENDIF
+      IF (ln_dpsiv ) THEN
+          CALL lbc_lnk(dpsiv,'T',1._wp)
+#ifdef key_iomput
+          CALL iom_put( "bsftv" , dpsiv )
+#endif
+      ENDIF
+
+      IF(kt .EQ. nitend ) THEN
+        IF (ln_dpsiv ) DEALLOCATE ( dpsiv )
+        IF (ln_dpsiu ) DEALLOCATE ( dpsiu )
+      ENDIF
+
+   END SUBROUTINE
+
+   SUBROUTINE skeb_dnum ( mt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE skeb_dnum  ***
+      !!
+      !! ** Purpose :   Computation of numerical energy dissipation
+      !!                For later use in the SKEB scheme
+      !!----------------------------------------------------------------------
+   INTEGER, INTENT(IN)  :: mt
+   REAL :: ds,dt,dtot,kh2
+   INTEGER :: ji,jj,jk
+   
+   IF ( mt .eq. nit000 ) THEN
+        ALLOCATE ( dnum(jpi,jpj,jpk) )
+        dnum (:,:,: ) = 0._wp
+   ENDIF
+
+   kh2 = rn_kh * rn_kh
+
+   IF( mt .eq. nit000 .OR. MOD( mt - 1, nn_dcom_freq ) == 0 ) THEN
+
+     DO jk=1,jpkm1
+      DO jj=2,jpj
+       DO ji=2,jpi
+          ! Shear
+          ds = (vn(ji,jj,jk)-vn(ji-1,jj,jk))*vmask(ji,jj,jk)*vmask(ji-1,jj,jk)/ e1v(ji,jj) + &
+               (un(ji,jj,jk)-un(ji,jj-1,jk))*umask(ji,jj,jk)*umask(ji,jj-1,jk)/ e2u(ji,jj)
+          ! Tension
+          dt = (vn(ji,jj,jk)-vn(ji-1,jj,jk))*vmask(ji,jj,jk)*vmask(ji-1,jj,jk)/ e2v(ji,jj) + &
+               (un(ji,jj,jk)-un(ji,jj-1,jk))*umask(ji,jj,jk)*umask(ji,jj-1,jk)/ e1u(ji,jj)
+  
+          dtot = sqrt( ds*ds + dt*dt ) * tmask(ji,jj,jk)
+          dnum(ji,jj,jk) = dtot*dtot*dtot*kh2*e1t(ji,jj)*e2t(ji,jj)
+       ENDDO
+      ENDDO
+     ENDDO
+   
+     CALL lbc_lnk(dnum,'T',1._wp)
+
+   ENDIF
+
+   END SUBROUTINE 
+
+   SUBROUTINE skeb_dcon ( mt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE skeb_dcon  ***
+      !!
+      !! ** Purpose :   Computation of convective energy dissipation
+      !!                For later use in the SKEB scheme
+      !!                Two formulation are implemented, based on buoyancy or
+      !!                convective mass flux formulations, respectively
+      !!----------------------------------------------------------------------
+   INTEGER, INTENT(IN)  :: mt
+   REAL(wp) :: zz, mf1,mf2,kc2
+   INTEGER  :: ji,jj,jk
+
+   IF ( mt .eq. nit000 ) THEN
+        ALLOCATE ( dcon(jpi,jpj,jpk) )
+        dcon (:,:,: ) = 0._wp
+   ENDIF
+
+   kc2 = rn_kc * rn_kc
+
+   IF( mt .eq. nit000 .OR. MOD( mt - 1, nn_dcom_freq ) == 0 ) THEN
+
+    IF(nn_dconv  .eq. 1 ) THEN
+
+     DO jk=2,jpkm1
+      DO jj=1,jpj
+       DO ji=1,jpi
+
+           zz = - grav*avt(ji,jj,jk) * ( rhd(ji,jj,jk)-rhd(ji,jj,jk-1) ) * wmask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) &
+              & / ( rau0 * fse3w(ji,jj,jk) ) 
+
+           dcon(ji,jj,jk) = kc2*zz*e1t(ji,jj)*e2t(ji,jj)*rau0 / fse3w(ji,jj,jk)
+
+       ENDDO
+      ENDDO
+     ENDDO
+
+    ELSEIF (nn_dconv .eq. 2 ) THEN
+
+     DO jk=2,jpkm1
+      DO jj=1,jpj
+       DO ji=1,jpi
+
+           mf1 = wn(ji,jj,jk+1)*e1t(ji,jj)*e2t(ji,jj)
+           mf2 = wn(ji,jj,jk-1)*e1t(ji,jj)*e2t(ji,jj)
+           dcon(ji,jj,jk) = rn_kc * (mf1-mf2)*(mf1-mf2) * tmask(ji,jj,jk+1) / (fse3w(ji,jj,jk)*rau0*rau0)
+
+       ENDDO
+      ENDDO
+     ENDDO
+
+    ENDIF
+
+     CALL lbc_lnk(dcon,'T',1._wp)
+
+   ENDIF
+
+   END SUBROUTINE
+
+   SUBROUTINE skeb_apply ( mt)
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE skeb_apply  ***
+      !!
+      !! ** Purpose :   Application of SKEB perturbation
+      !!                Convective and Numerical energy dissipation are
+      !!                multiplied by the beta terms
+      !!----------------------------------------------------------------------
+
+   INTEGER, INTENT(IN)  :: mt
+   INTEGER :: ji,jj,jk
+   REAL(wp), DIMENSION(jpi,jpj,jpk) :: ustar,vstar,psi
+   REAL(wp) :: mi,ma
+
+   IF( ln_stopack_diags ) THEN
+
+     psi = 0._wp
+     IF(ln_skeb_own_gauss) THEN
+       DO jk=1,jpkm1
+         psi(:,:,jk) = rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk) ) * gauss_n_2d_k(:,:) 
+       ENDDO
+     ELSE
+       DO jk=1,jpkm1
+         psi(:,:,jk) = rn_skeb_stdev * rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk) ) &
+         & * gauss_n_2d(:,:) / rn_sppt_stdev
+       ENDDO
+     ENDIF
+     ustar = 0._wp
+     vstar = 0._wp
+     DO jk=1,jpkm1
+      DO jj=2,jpj
+       DO ji=2,jpi
+          ustar(ji,jj,jk) =   ( psi(ji,jj,jk)-psi(ji,jj-1,jk) )*umask(ji,jj,jk)*umask(ji,jj-1,jk)/ e2u(ji,jj)
+          vstar(ji,jj,jk) = - ( psi(ji,jj,jk)-psi(ji-1,jj,jk) )*vmask(ji,jj,jk)*vmask(ji-1,jj,jk)/ e1v(ji,jj)
+       ENDDO
+      ENDDO
+     ENDDO
+     mi = MAXVAL(ABS(ustar))
+     ma = MAXVAL(ABS(vstar))
+     IF(lk_mpp) CALL mpp_max(mi)
+     IF(lk_mpp) CALL mpp_max(ma)
+     IF(lwp) WRITE(numdiag,9304) lkt,'DNUM ABS CURRENT' ,mi,ma
+
+     rn_mmax ( jk_skeb_dnum ) =  MAX ( MAX( mi, ma), rn_mmax ( jk_skeb_dnum ) )
+
+     psi = 0._wp
+     IF(ln_skeb_own_gauss) THEN
+       DO jk=1,jpkm1
+         psi(:,:,jk) = rn_skeb * sqrt( rn_beta_con * dcon(:,:,jk) ) * gauss_n_2d_k(:,:) 
+       ENDDO
+     ELSE
+       DO jk=1,jpkm1
+         psi(:,:,jk) = rn_skeb_stdev * rn_skeb * sqrt( rn_beta_con * dcon(:,:,jk) ) &
+         & * gauss_n_2d(:,:) / rn_sppt_stdev
+       ENDDO
+     ENDIF
+     ustar = 0._wp
+     vstar = 0._wp
+     DO jk=1,jpkm1
+      DO jj=2,jpj
+       DO ji=2,jpi
+          ustar(ji,jj,jk) =   ( psi(ji,jj,jk)-psi(ji,jj-1,jk) )*umask(ji,jj,jk)*umask(ji,jj-1,jk)/ e2u(ji,jj)
+          vstar(ji,jj,jk) = - ( psi(ji,jj,jk)-psi(ji-1,jj,jk) )*vmask(ji,jj,jk)*vmask(ji-1,jj,jk)/ e1v(ji,jj)
+       ENDDO
+      ENDDO
+     ENDDO
+     mi = MAXVAL(ABS(ustar))
+     ma = MAXVAL(ABS(vstar))
+     IF(lk_mpp) CALL mpp_max(mi)
+     IF(lk_mpp) CALL mpp_max(ma)
+     IF(lwp) WRITE(numdiag,9304) lkt,'DCON ABS CURRENT' ,mi,ma
+
+     rn_mmax ( jk_skeb_dcon ) =  MAX ( MAX( mi, ma), rn_mmax ( jk_skeb_dcon ) )
+
+9304 FORMAT(' it :', i8, ' ', A16, ' Min: ',d10.3 , ' Max: ',d10.3)
+
+   ENDIF
+
+   psi = 0._wp
+   IF(ln_skeb_own_gauss) THEN
+     DO jk=1,jpkm1
+       psi(:,:,jk) = rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk)+ rn_beta_con * dcon(:,:,jk) ) * gauss_n_2d_k(:,:) 
+     ENDDO
+   ELSE
+     DO jk=1,jpkm1
+       psi(:,:,jk) = rn_skeb_stdev * rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk)+ rn_beta_con * dcon(:,:,jk) ) &
+       & * gauss_n_2d(:,:) / rn_sppt_stdev
+     ENDDO
+   ENDIF
+
+   ustar = 0._wp
+   vstar = 0._wp
+
+   DO jk=1,jpkm1
+    DO jj=2,jpj
+     DO ji=2,jpi
+        ustar(ji,jj,jk) =   ( psi(ji,jj,jk)-psi(ji,jj-1,jk) )*umask(ji,jj,jk)*umask(ji,jj-1,jk)/ e2u(ji,jj)
+        vstar(ji,jj,jk) = - ( psi(ji,jj,jk)-psi(ji-1,jj,jk) )*vmask(ji,jj,jk)*vmask(ji-1,jj,jk)/ e1v(ji,jj)
+     ENDDO
+    ENDDO
+   ENDDO
+
+   IF( ln_stopack_diags ) THEN
+
+      mi = MAXVAL(ABS(dnum))
+      ma = MAXVAL(ABS(dcon))
+      IF(lk_mpp) CALL mpp_max(mi)
+      IF(lk_mpp) CALL mpp_max(ma)
+      IF(lwp) WRITE(numdiag,9302) lkt,'DNUM AND DCON MX' ,mi,ma
+
+      mi = MINVAL(ustar)
+      ma = MAXVAL(ustar)
+      IF(lk_mpp) CALL mpp_min(mi)
+      IF(lk_mpp) CALL mpp_max(ma)
+      IF(lwp) WRITE(numdiag,9302) lkt,'SKEB U-CURRENT  ' ,mi,ma
+
+      rn_mmax ( jk_skeb_tot ) =  MAX ( MAX( abs(mi), ma), rn_mmax ( jk_skeb_tot ) )
+
+      mi = MINVAL(vstar)
+      ma = MAXVAL(vstar)
+      IF(lk_mpp) CALL mpp_min(mi)
+      IF(lk_mpp) CALL mpp_max(ma)
+      IF(lwp) WRITE(numdiag,9302) lkt,'SKEB V-CURRENT  ' ,mi,ma
+9302  FORMAT(' it :', i8, ' ', A16, ' Min: ',d10.3 , ' Max: ',d10.3)
+
+      rn_mmax ( jk_skeb_tot ) =  MAX ( MAX( abs(mi), ma), rn_mmax ( jk_skeb_tot ) )
+
+   ENDIF
+
+   CALL lbc_lnk(ustar,'U',-1._wp)
+   CALL lbc_lnk(vstar,'V',-1._wp)
+
+   DO jk=1,jpkm1
+    DO jj=1,jpj
+     DO ji=1,jpi
+        un(ji,jj,jk)    = un(ji,jj,jk) + ustar(ji,jj,jk)
+        vn(ji,jj,jk)    = vn(ji,jj,jk) + vstar(ji,jj,jk)
+     ENDDO
+    ENDDO
+   ENDDO
+
+#ifdef key_iomput
+   IF (iom_use('ustar_skeb') ) CALL iom_put( "ustar_skeb" , ustar)
+   IF (iom_use('vstar_skeb') ) CALL iom_put( "vstar_skeb" , vstar)
+#endif
+
+   IF ( mt .eq. nitend ) THEN
+     IF(ALLOCATED(dnum)) DEALLOCATE ( dnum )
+     IF(ALLOCATED(dcon)) DEALLOCATE ( dcon )
+     IF (ln_stopack_diags .AND. lwp) CALL stopack_report
+   ENDIF
+
+   END SUBROUTINE
+
+   !!======================================================================
+   !! Random number generator, used in NEMO stochastic parameterization
+   !!
+   !!=====================================================================
+   !! History :  3.3  ! 2011-10 (J.-M. Brankart)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !! The module is based on (and includes) the
+   !! 64-bit KISS (Keep It Simple Stupid) random number generator
+   !! distributed by George Marsaglia :
+   !! http://groups.google.com/group/comp.lang.fortran/
+   !!        browse_thread/thread/a85bf5f2a97f5a55
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   kiss          : 64-bit KISS random number generator (period ~ 2^250)
+   !!   kiss_seed     : Define seeds for KISS random number generator
+   !!   kiss_state    : Get current state of KISS random number generator
+   !!   kiss_save     : Save current state of KISS (for future restart)
+   !!   kiss_load     : Load the saved state of KISS
+   !!   kiss_reset    : Reset the default seeds
+   !!   kiss_check    : Check the KISS pseudo-random sequence
+   !!   kiss_uniform  : Real random numbers with uniform distribution in [0,1]
+   !!   kiss_gaussian : Real random numbers with Gaussian distribution N(0,1)
+   !!   kiss_gamma    : Real random numbers with Gamma distribution Gamma(k,1)
+   !!   kiss_sample   : Select a random sample from a set of integers
+   !!
+   !!   ---CURRENTLY NOT USED IN NEMO :
+   !!   kiss_save, kiss_load, kiss_check, kiss_gamma, kiss_sample
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: dynhpg.F90 2528 2010-12-27 17:33:53Z rblod $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+   FUNCTION kiss()
+      !! --------------------------------------------------------------------
+      !!                  ***  FUNCTION kiss  ***
+      !!
+      !! ** Purpose :   64-bit KISS random number generator
+      !!
+      !! ** Method  :   combine several random number generators:
+      !!                (1) Xorshift (XSH), period 2^64-1,
+      !!                (2) Multiply-with-carry (MWC), period (2^121+2^63-1)
+      !!                (3) Congruential generator (CNG), period 2^64.
+      !!
+      !!                overall period:
+      !!                (2^250+2^192+2^64-2^186-2^129)/6
+      !!                            ~= 2^(247.42) or 10^(74.48)
+      !!
+      !!                set your own seeds with 'kiss_seed'
+      ! --------------------------------------------------------------------
+      IMPLICIT NONE
+      INTEGER(KIND=i8) :: kiss, t
+
+      t = ISHFT(x,58) + w
+      IF (s(x).eq.s(t)) THEN
+         w = ISHFT(x,-6) + s(x)
+      ELSE
+         w = ISHFT(x,-6) + 1 - s(x+t)
+      ENDIF
+      x = t + x
+      y = m( m( m(y,13_8), -17_8 ), 43_8 )
+      z = 6906969069_8 * z + 1234567_8
+
+      kiss = x + y + z
+
+      CONTAINS
+
+         FUNCTION s(k)
+            INTEGER(KIND=i8) :: s, k
+            s = ISHFT(k,-63)
+         END FUNCTION s
+
+         FUNCTION m(k, n)
+            INTEGER(KIND=i8) :: m, k, n
+            m =  IEOR(k, ISHFT(k, n) )
+         END FUNCTION m
+
+   END FUNCTION kiss
+
+
+   SUBROUTINE kiss_seed(ix, iy, iz, iw)
+      !! --------------------------------------------------------------------
+      !!                  ***  ROUTINE kiss_seed  ***
+      !!
+      !! ** Purpose :   Define seeds for KISS random number generator
+      !!
+      !! --------------------------------------------------------------------
+      IMPLICIT NONE
+      INTEGER(KIND=i8) :: ix, iy, iz, iw
+
+      IF(lwp) WRITE(numout,*) ' *** kiss_seed called with args:'
+      IF(lwp) WRITE(numout,*) ' *** ix = ',ix
+      IF(lwp) WRITE(numout,*) ' *** iy = ',iy
+      IF(lwp) WRITE(numout,*) ' *** iz = ',iz
+      IF(lwp) WRITE(numout,*) ' *** iw = ',iw
+
+      x = ix
+      y = iy
+      z = iz
+      w = iw
+
+   END SUBROUTINE kiss_seed
+
+
+   SUBROUTINE kiss_state(ix, iy, iz, iw)
+      !! --------------------------------------------------------------------
+      !!                  ***  ROUTINE kiss_state  ***
+      !!
+      !! ** Purpose :   Get current state of KISS random number generator
+      !!
+      !! --------------------------------------------------------------------
+      IMPLICIT NONE
+      INTEGER(KIND=i8) :: ix, iy, iz, iw
+
+      ix = x
+      iy = y
+      iz = z
+      iw = w
+
+   END SUBROUTINE kiss_state
+
+
+   SUBROUTINE kiss_reset()
+      !! --------------------------------------------------------------------
+      !!                  ***  ROUTINE kiss_reset  ***
+      !!
+      !! ** Purpose :   Reset the default seeds for KISS random number generator
+      !!
+      !! --------------------------------------------------------------------
+      IMPLICIT NONE
+
+      x=1234567890987654321_8
+      y=362436362436362436_8
+      z=1066149217761810_8
+      w=123456123456123456_8
+
+    END SUBROUTINE kiss_reset
+
+   SUBROUTINE kiss_uniform(uran)
+      !! --------------------------------------------------------------------
+      !!                  ***  ROUTINE kiss_uniform  ***
+      !!
+      !! ** Purpose :   Real random numbers with uniform distribution in [0,1]
+      !!
+      !! --------------------------------------------------------------------
+      IMPLICIT NONE
+      REAL(KIND=wp) :: uran
+
+      uran = half * ( one + REAL(kiss(),wp) / huge64 )
+
+   END SUBROUTINE kiss_uniform
+
+
+   SUBROUTINE kiss_gaussian(gran)
+      !! --------------------------------------------------------------------
+      !!                  ***  ROUTINE kiss_gaussian  ***
+      !!
+      !! ** Purpose :   Real random numbers with Gaussian distribution N(0,1)
+      !!
+      !! ** Method  :   Generate 2 new Gaussian draws (gran1 and gran2)
+      !!                from 2 uniform draws on [-1,1] (u1 and u2),
+      !!                using the Marsaglia polar method
+      !!                (see Devroye, Non-Uniform Random Variate Generation, p. 235-236)
+      !!
+      !! --------------------------------------------------------------------
+      IMPLICIT NONE
+      REAL(KIND=wp) :: gran, u1, u2, rsq, fac
+
+      IF (ig.EQ.1) THEN
+         rsq = two
+         DO WHILE ( (rsq.GE.one).OR. (rsq.EQ.zero) )
+            u1 = REAL(kiss(),wp) / huge64
+            u2 = REAL(kiss(),wp) / huge64
+            rsq = u1*u1 + u2*u2
+         ENDDO
+         fac = SQRT(-two*LOG(rsq)/rsq)
+         gran1 = u1 * fac
+         gran2 = u2 * fac
+      ENDIF
+
+      ! Output one of the 2 draws
+      IF (ig.EQ.1) THEN
+         gran = gran1 ; ig = 2
+      ELSE
+         gran = gran2 ; ig = 1
+      ENDIF
+
+   END SUBROUTINE kiss_gaussian
+
+
+   SUBROUTINE kiss_gamma(gamr,k)
+      !! --------------------------------------------------------------------
+      !!                  ***  ROUTINE kiss_gamma  ***
+      !!
+      !! ** Purpose :   Real random numbers with Gamma distribution Gamma(k,1)
+      !!
+      !! --------------------------------------------------------------------
+      IMPLICIT NONE
+      REAL(KIND=wp), PARAMETER :: p1 = 4.5_8
+      REAL(KIND=wp), PARAMETER :: p2 = 2.50407739677627_8  ! 1+LOG(9/2)
+      REAL(KIND=wp), PARAMETER :: p3 = 1.38629436111989_8  ! LOG(4)
+      REAL(KIND=wp) :: gamr, k, u1, u2, b, c, d, xx, yy, zz, rr, ee
+      LOGICAL :: accepted
+
+      IF (k.GT.one) THEN
+         ! Cheng's rejection algorithm
+         ! (see Devroye, Non-Uniform Random Variate Generation, p. 413)
+         b = k - p3 ; d = SQRT(two*k-one) ; c = k + d
+
+         accepted=.FALSE.
+         DO WHILE (.NOT.accepted)
+            CALL kiss_uniform(u1)
+            yy = LOG(u1/(one-u1)) / d  ! Mistake in Devroye: "* k" instead of "/ d"
+            xx = k * EXP(yy)
+            rr = b + c * yy - xx
+            CALL kiss_uniform(u2)
+            zz = u1 * u1 * u2
+
+            accepted = rr .GE. (zz*p1-p2)
+            IF (.NOT.accepted) accepted =  rr .GE. LOG(zz)
+         ENDDO
+
+         gamr = xx
+
+      ELSEIF (k.LT.one) THEN
+        ! Rejection from the Weibull density
+        ! (see Devroye, Non-Uniform Random Variate Generation, p. 415)
+        c = one/k ; d = (one-k) * EXP( (k/(one-k)) * LOG(k) )
+
+        accepted=.FALSE.
+        DO WHILE (.NOT.accepted)
+           CALL kiss_uniform(u1)
+           zz = -LOG(u1)
+           xx = EXP( c * LOG(zz) )
+           CALL kiss_uniform(u2)
+           ee = -LOG(u2)
+
+           accepted = (zz+ee) .GE. (d+xx)  ! Mistake in Devroye: "LE" instead of "GE"
+        ENDDO
+
+        gamr = xx
+
+      ELSE
+         ! Exponential distribution
+         CALL kiss_uniform(u1)
+         gamr = -LOG(u1)
+
+      ENDIF
+
+   END SUBROUTINE kiss_gamma
+
+
+   SUBROUTINE kiss_sample(a,n,k)
+      !! --------------------------------------------------------------------
+      !!                  ***  ROUTINE kiss_sample  ***
+      !!
+      !! ** Purpose :   Select a random sample of size k from a set of n integers
+      !!
+      !! ** Method  :   The sample is output in the first k elements of a
+      !!                Set k equal to n to obtain a random permutation
+      !!                  of the whole set of integers
+      !!
+      !! --------------------------------------------------------------------
+      IMPLICIT NONE
+      INTEGER(KIND=i8), DIMENSION(:) :: a
+      INTEGER(KIND=i8) :: n, k, i, j, atmp
+      REAL(KIND=wp) :: uran
+
+      ! Select the sample using the swapping method
+      ! (see Devroye, Non-Uniform Random Variate Generation, p. 612)
+      DO i=1,k
+         ! Randomly select the swapping element between i and n (inclusive)
+         CALL kiss_uniform(uran)
+         j = i - 1 + CEILING( REAL(n-i+1,8) * uran )
+         ! Swap elements i and j
+         atmp = a(i) ; a(i) = a(j) ; a(j) = atmp
+      ENDDO
+
+   END SUBROUTINE kiss_sample
+
+#ifdef NEMO_V34
+SUBROUTINE ctl_nam(ier,cstr,lout)
+INTEGER, INTENT(IN)            :: ier
+LOGICAL, INTENT(IN)            :: lout
+CHARACTER(LEN=*),INTENT(IN)    :: cstr
+IF(lout) WRITE(numout,'(A,I4,X,A)') 'Error ',ier,TRIM(cstr)
+END SUBROUTINE
+#endif
+
+END MODULE stopack
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90	(revision 11393)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90	(revision 11394)
@@ -28,4 +28,5 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE stopack
 
    IMPLICIT NONE
@@ -41,4 +42,5 @@
 
    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   qgh_trd0   ! geothermal heating trend
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   qgh_trd1   ! geothermal heating trend
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qgh              ! structure of input qgh (file informations, fields read)
  
@@ -89,8 +91,13 @@
       !
       !                             !  Add the geothermal heat flux trend on temperature
+
+      IF( nn_spp_geot .GT. 0) THEN
+          qgh_trd1 = qgh_trd0
+          CALL spp_gen(kt, qgh_trd0, nn_spp_geot, rn_geot_sd, jk_spp_geot)
+      ENDIF
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
             ik = mbkt(ji,jj)
-            zqgh_trd = qgh_trd0(ji,jj) / fse3t(ji,jj,ik)
+            zqgh_trd = qgh_trd1(ji,jj) / fse3t(ji,jj,ik)
             tsa(ji,jj,ik,jp_tem) = tsa(ji,jj,ik,jp_tem) + zqgh_trd
          END DO
@@ -163,4 +170,5 @@
          !
          ALLOCATE( qgh_trd0(jpi,jpj) )    ! allocation
+         ALLOCATE( qgh_trd1(jpi,jpj) )    ! allocation
          !
          SELECT CASE ( nn_geoflx )        ! geothermal heat flux / (rauO * Cp)
@@ -192,4 +200,5 @@
             !
          END SELECT
+         qgh_trd1 = qgh_trd0
          !
       ELSE
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 11393)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 11394)
@@ -39,4 +39,5 @@
    USE timing         ! Timing
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE stopack
 
    IMPLICIT NONE
@@ -67,4 +68,5 @@
    INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   mgrhu    , mgrhv       ! = +/-1, sign of grad(H) in u-(v-)direction (PUBLIC for TAM)
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ahu_bbl_0, ahv_bbl_0   ! diffusive bbl flux coefficients at u and v-points
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ahu_bbl_1, ahv_bbl_1   ! diffusive bbl flux coefficients at u and v-points
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   e3u_bbl_0, e3v_bbl_0   ! thichness of the bbl (e3) at u and v-points (PUBLIC for TAM)
 
@@ -86,4 +88,5 @@
          &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d  (jpi,jpj) , mgrhv(jpi,jpj) ,     &
          &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                          &
+         &      ahu_bbl_1(jpi,jpj) , ahv_bbl_1(jpi,jpj) ,                                          &
          &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) ,                                      STAT=tra_bbl_alloc )
          !
@@ -195,4 +198,13 @@
       ALLOCATE(zptb(1:jpi, 1:jpj))
       !
+      ahu_bbl_1 = ahu_bbl
+      IF( nn_spp_ahubbl .GT. 0 ) THEN
+          CALL spp_gen(1 , ahu_bbl_1, nn_spp_ahubbl, rn_ahubbl_sd, jk_spp_ahubbl )
+      ENDIF
+      ahv_bbl_1 = ahv_bbl
+      IF( nn_spp_ahvbbl .GT. 0 ) THEN
+          CALL spp_gen(1 , ahv_bbl_1, nn_spp_ahvbbl, rn_ahvbbl_sd, jk_spp_ahvbbl )
+      ENDIF
+      !
       DO jn = 1, kjpt                                     ! tracer loop
          !                                                ! ===========
@@ -209,8 +221,8 @@
                zbtr = r1_e12t(ji,jj)  / fse3t(ji,jj,ik)
                pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                         &
-                  &               + (   ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )   &
-                  &                   - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )   &
-                  &                   + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )   &
-                  &                   - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )   ) * zbtr
+                  &               + (   ahu_bbl_1(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )   &
+                  &                   - ahu_bbl_1(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )   &
+                  &                   + ahv_bbl_1(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )   &
+                  &                   - ahv_bbl_1(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )   ) * zbtr
             END DO
          END DO
@@ -549,10 +561,15 @@
       zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
 
-                                        !* sign of grad(H) at u- and v-points
-      mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0
+                                        !* sign of grad(H) at u- and v-points; zero if grad(H) = 0
+      mgrhu(:,:) = 0   ;   mgrhv(:,:) = 0
       DO jj = 1, jpjm1
          DO ji = 1, jpim1
-            mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
-            mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
+            IF( gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
+               mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
+            ENDIF
+            !
+            IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
+               mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
+            ENDIF
          END DO
       END DO
@@ -569,5 +586,4 @@
       ahu_bbl_0(:,:) = rn_ahtbbl * e2u(:,:) * e3u_bbl_0(:,:) / e1u(:,:)  * umask(:,:,1)
       ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:)  * vmask(:,:,1)
-
 
       IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : regional enhancement of ah_bbl
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90	(revision 11393)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90	(revision 11394)
@@ -32,4 +32,5 @@
    USE wrk_nemo        ! Memory allocation
    USE timing          ! Timing
+   USE stopack
 
    IMPLICIT NONE
@@ -43,4 +44,10 @@
    REAL, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   t0_ldf, s0_ldf   !: lateral diffusion trends of T & S for a cst profile
    !                                                               !  (key_traldf_ano only)
+#if defined key_traldf_c3d
+   REAL, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: ahtu0,ahtv0,ahtw0,ahtt0
+#endif
+#if defined key_traldf_c2d
+   REAL, SAVE, ALLOCATABLE, DIMENSION(:,:  ) :: ahtu0,ahtv0,ahtw0,ahtt0
+#endif
 
    !! * Substitutions
@@ -75,4 +82,43 @@
          ztrds(:,:,:) = tsa(:,:,:,jp_sal)
       ENDIF
+
+#if defined key_traldf_c3d
+         IF( kt .eq. nit000 .and. (nn_spp_ahtu+nn_spp_ahtv+nn_spp_ahtw+nn_spp_ahtt) .gt. 0 ) THEN
+             ALLOCATE ( ahtu0(jpi,jpj,jpk), ahtv0(jpi,jpj,jpk) )
+             ALLOCATE ( ahtt0(jpi,jpj,jpk), ahtw0(jpi,jpj,jpk) )
+             ahtu0 = ahtu
+             ahtv0 = ahtv
+             ahtw0 = ahtw
+             ahtt0 = ahtt
+         ENDIF
+#endif
+#if defined key_traldf_c2d
+         IF( kt .eq. nit000 .and. (nn_spp_ahtu+nn_spp_ahtv+nn_spp_ahtw+nn_spp_ahtt) .gt. 0 ) THEN
+             ALLOCATE ( ahtu0(jpi,jpj), ahtv0(jpi,jpj) )
+             ALLOCATE ( ahtt0(jpi,jpj), ahtw0(jpi,jpj) )
+             ahtu0 = ahtu
+             ahtv0 = ahtv
+             ahtw0 = ahtw
+             ahtt0 = ahtt
+         ENDIF
+#endif
+#if defined key_traldf_c3d || defined key_traldf_c2d
+         IF( nn_spp_ahtu .GT. 0) THEN
+             ahtu = ahtu0
+             CALL spp_aht(kt,ahtu,nn_spp_ahtu,rn_ahtu_sd,jk_spp_ahtu)
+         ENDIF
+         IF( nn_spp_ahtv .GT. 0) THEN
+             ahtv = ahtv0
+             CALL spp_aht(kt,ahtv,nn_spp_ahtv,rn_ahtv_sd,jk_spp_ahtv)
+         ENDIF
+         IF( nn_spp_ahtw .GT. 0) THEN
+             ahtw = ahtw0
+             CALL spp_aht(kt,ahtw,nn_spp_ahtw,rn_ahtw_sd,jk_spp_ahtw)
+         ENDIF
+         IF( nn_spp_ahtt .GT. 0) THEN
+             ahtt = ahtt0
+             CALL spp_aht(kt,ahtt,nn_spp_ahtt,rn_ahtt_sd,jk_spp_ahtt)
+         ENDIF
+#endif
 
       SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend
@@ -215,4 +261,6 @@
       IF( ierr == 1 )   CALL ctl_stop( ' iso-level in z-coordinate - partial step, not allowed' )
       IF( ierr == 2 )   CALL ctl_stop( ' isoneutral bilaplacian operator does not exist' )
+      IF( ln_traldf_hor .AND. ln_traldf_grif )    &
+            &   CALL ctl_stop( ' horizontal operator and Griffies triads not available; sitch to isoneutral operator' )
       IF( ln_traldf_grif .AND. ln_isfcav         )   &
            CALL ctl_stop( ' ice shelf and traldf_grif not tested')
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 11393)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 11394)
@@ -33,4 +33,5 @@
    USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
+   USE stopack
 
    IMPLICIT NONE
@@ -52,5 +53,5 @@
  
    ! Module variables
-   REAL(wp) ::   xsi0r                           !: inverse of rn_si0
+   REAL(wp), ALLOCATABLE ::   xsi0r(:,:)         !: inverse of rn_si0
    REAL(wp) ::   xsi1r                           !: inverse of rn_si1
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read)
@@ -182,5 +183,11 @@
          !                                        ! ============================================== !
          !
-         !                                                ! ------------------------- !
+         ! 
+         IF( nn_spp_qsi0 > 0 ) THEN
+             xsi0r = rn_si0
+             CALL spp_gen(kt, xsi0r, nn_spp_qsi0, rn_qsi0_sd, jk_spp_qsi0 )
+             xsi0r = 1.e0 / xsi0r
+         ENDIF
+         !                                               ! ------------------------- !
          IF( ln_qsr_rgb) THEN                             !  R-G-B  light penetration !
             !                                             ! ------------------------- !
@@ -251,5 +258,5 @@
 !CDIR NOVERRCHK   
                      DO ji = 1, jpi
-                        zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r     )
+                        zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r(ji,jj) )
                         zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekb(ji,jj) )
                         zc2 = ze2(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekg(ji,jj) )
@@ -263,4 +270,16 @@
                   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(ji,jj) )
+                        zzc1 = zcoef  * EXP( - fse3t(ji,jj,1) * zekb(ji,jj) )
+                        zzc2 = zcoef  * EXP( - fse3t(ji,jj,1) * zekg(ji,jj) )
+                        zzc3 = zcoef  * EXP( - fse3t(ji,jj,1) * zekr(ji,jj) )
+                        fraqsr_1lev(ji,jj) = 1.0 - ( zzc0 + zzc1 + zzc2  + zzc3  ) * tmask(ji,jj,2) 
+                     END DO
+                  END DO
+               ENDIF
                !
                DO jk = 1, nksr                                        ! compute and add qsr trend to ta
@@ -310,5 +329,6 @@
             !                                             ! ------------------------- !
             !
-            IF( lk_vvl ) THEN                                  !* variable volume
+            IF( lk_vvl .OR. nn_spp_qsi0 > 0 ) THEN        !* variable volume
+
                zz0   =        rn_abs   * r1_rau0_rcp
                zz1   = ( 1. - rn_abs ) * r1_rau0_rcp
@@ -316,6 +336,6 @@
                   DO jj = 1, jpj
                      DO ji = 1, jpi
-                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
-                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
+                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
+                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
                         qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) ) 
                      END DO
@@ -326,6 +346,6 @@
                   DO jj = 1, jpj
                      DO ji = 1, jpi
-                        zc0 = zz0 * EXP( -fsdepw(ji,jj,1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,1)*xsi1r )
-                        zc1 = zz0 * EXP( -fsdepw(ji,jj,2)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,2)*xsi1r )
+                        zc0 = zz0 * EXP( -fsdepw(ji,jj,1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,1)*xsi1r )
+                        zc1 = zz0 * EXP( -fsdepw(ji,jj,2)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,2)*xsi1r )
                         fraqsr_1lev(ji,jj) = ( zc0*tmask(ji,jj,1) - zc1*tmask(ji,jj,2) ) / r1_rau0_rcp
                      END DO
@@ -490,4 +510,5 @@
          !                       ! ===================================== !
          !
+         ALLOCATE( xsi0r(jpi,jpj) )
          xsi0r = 1.e0 / rn_si0
          xsi1r = 1.e0 / rn_si1
@@ -544,5 +565,5 @@
 !CDIR NOVERRCHK   
                         DO ji = 1, jpi
-                           zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r     )
+                           zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r(ji,jj) )
                            zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekb(ji,jj) )
                            zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekg(ji,jj) )
@@ -585,6 +606,6 @@
                   DO jj = 1, jpj                              ! top 400 meters
                      DO ji = 1, jpi
-                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
-                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
+                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
+                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
                         etot3(ji,jj,jk) = (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) * tmask(ji,jj,1) 
                      END DO
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRD/trddyn.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRD/trddyn.F90	(revision 11393)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRD/trddyn.F90	(revision 11394)
@@ -29,4 +29,5 @@
    USE lib_mpp        ! MPP library
    USE wrk_nemo       ! Memory allocation
+   USE stopack     
 
    IMPLICIT NONE
@@ -68,4 +69,5 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       IF( ln_dyn_trd )   CALL trd_dyn_iom( putrd, pvtrd, ktrd, kt )
+      IF( ln_dyn_trd .AND. ln_sppt_dyn ) CALL dyn_sppt_collect( putrd, pvtrd, ktrd, kt )
          
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90	(revision 11393)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90	(revision 11394)
@@ -29,4 +29,5 @@
    USE iom            ! I/O manager library
    USE lib_mpp        ! MPP library
+   USE stopack
    USE wrk_nemo       ! Memory allocation
 
@@ -249,6 +250,7 @@
       !                   ! 3D output of tracers trends using IOM interface
       IF( ln_tra_trd )   CALL trd_tra_iom ( ptrdx, ptrdy, ktrd, kt )
-
-      !                   ! Integral Constraints Properties for tracers trends                                       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( ln_tra_trd .AND. ln_sppt_tra )  CALL tra_sppt_collect( ptrdx, ptrdy, ktrd, kt )
+
+      !  Integral Constraints Properties for tracers trends                                       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       IF( ln_glo_trd )   CALL trd_glo( ptrdx, ptrdy, ktrd, 'TRA', kt )
 
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 11393)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 11394)
@@ -26,4 +26,5 @@
    USE wrk_nemo        ! Memory Allocation
    USE phycst, ONLY: vkarmn
+   USE stopack
 
    IMPLICIT NONE
@@ -52,5 +53,5 @@
    REAL(wp), PUBLIC ::   rn_tfrz0     ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM)
    LOGICAL , PUBLIC ::   ln_bfrimp    ! logical switch for implicit bottom friction
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::  bfrcoef2d, tfrcoef2d   ! 2D bottom/top drag coefficient (PUBLIC for TAM)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::  bfrcoef2d, tfrcoef2d,bfrcoef2d0   ! 2D bottom/top drag coefficient (PUBLIC for TAM)
 
    !! * Substitutions
@@ -68,5 +69,5 @@
       !!                ***  FUNCTION zdf_bfr_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( bfrcoef2d(jpi,jpj), tfrcoef2d(jpi,jpj), STAT=zdf_bfr_alloc )
+      ALLOCATE( bfrcoef2d(jpi,jpj), tfrcoef2d(jpi,jpj), bfrcoef2d0(jpi,jpj),STAT=zdf_bfr_alloc )
       !
       IF( lk_mpp             )   CALL mpp_sum ( zdf_bfr_alloc )
@@ -105,4 +106,9 @@
          WRITE(numout,*) 'zdf_bfr : Set bottom friction coefficient (non-linear case)'
          WRITE(numout,*) '~~~~~~~~'
+      ENDIF
+      ! 
+      IF( nn_spp_bfr > 0 ) THEN
+         bfrcoef2d = bfrcoef2d0
+         CALL spp_gen(kt, bfrcoef2d, nn_spp_bfr, rn_bfr_sd, jk_spp_bfr)
       ENDIF
       !
@@ -486,4 +492,5 @@
       ENDIF
       !
+      bfrcoef2d0(:,:) = bfrcoef2d(:,:)
       IF( nn_timing == 1 )  CALL timing_stop('zdf_bfr_init')
       !
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90	(revision 11393)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90	(revision 11394)
@@ -25,4 +25,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE timing          ! Timing
+   USE stopack         
 
    IMPLICIT NONE
@@ -30,4 +31,5 @@
 
    PUBLIC   zdf_evd    ! called by step.F90
+   REAL, SAVE, ALLOCATABLE, DIMENSION(:,:) :: rn_avevd0
 
    !! * Substitutions
@@ -69,7 +71,14 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~ '
          IF(lwp) WRITE(numout,*)
+         ALLOCATE ( rn_avevd0(jpi,jpj) )
+         rn_avevd0(:,:) = rn_avevd
       ENDIF
 
       zavt_evd(:,:,:) = avt(:,:,:)           ! set avt prior to evd application
+
+      IF(nn_spp_aevd.GT.0) THEN
+         rn_avevd0(:,:) = rn_avevd
+         CALL spp_gen(kt, rn_avevd0, nn_spp_aevd, rn_aevd_sd, jk_spp_aevd)
+      ENDIF
 
       SELECT CASE ( nn_evdm )
@@ -88,10 +97,10 @@
                   IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) THEN
 #endif
-                     avt (ji  ,jj  ,jk) = rn_avevd * tmask(ji  ,jj  ,jk)
-                     avm (ji  ,jj  ,jk) = rn_avevd * tmask(ji  ,jj  ,jk)
-                     avmu(ji  ,jj  ,jk) = rn_avevd * umask(ji  ,jj  ,jk)
-                     avmu(ji-1,jj  ,jk) = rn_avevd * umask(ji-1,jj  ,jk)
-                     avmv(ji  ,jj  ,jk) = rn_avevd * vmask(ji  ,jj  ,jk)
-                     avmv(ji  ,jj-1,jk) = rn_avevd * vmask(ji  ,jj-1,jk)
+                     avt (ji  ,jj  ,jk) = rn_avevd0(ji,jj) * tmask(ji  ,jj  ,jk)
+                     avm (ji  ,jj  ,jk) = rn_avevd0(ji,jj) * tmask(ji  ,jj  ,jk)
+                     avmu(ji  ,jj  ,jk) = rn_avevd0(ji,jj) * umask(ji  ,jj  ,jk)
+                     avmu(ji-1,jj  ,jk) = rn_avevd0(ji,jj) * umask(ji-1,jj  ,jk)
+                     avmv(ji  ,jj  ,jk) = rn_avevd0(ji,jj) * vmask(ji  ,jj  ,jk)
+                     avmv(ji  ,jj-1,jk) = rn_avevd0(ji,jj) * vmask(ji  ,jj-1,jk)
                   ENDIF
                END DO
@@ -115,5 +124,5 @@
                   IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 )   &
 #endif
-                     avt(ji,jj,jk) = rn_avevd * tmask(ji,jj,jk)
+                     avt(ji,jj,jk) = rn_avevd0(ji,jj) * tmask(ji,jj,jk)
                END DO
             END DO
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 11393)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 11394)
@@ -32,4 +32,5 @@
    USE timing         ! Timing
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE stopack
 
    IMPLICIT NONE
@@ -806,4 +807,6 @@
             END DO
          END DO
+         IF(nn_spp_avt > 0 ) CALL spp_gen(kt,avt(:,:,jk),nn_spp_avt,rn_avt_sd,jk)
+         IF(nn_spp_avm > 0 ) CALL spp_gen(kt,avm(:,:,jk),nn_spp_avm,rn_avm_sd,jk)
       END DO
       !
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 11393)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 11394)
@@ -57,6 +57,5 @@
    USE agrif_opa_update
 #endif
-
-
+   USE stopack
 
    IMPLICIT NONE
@@ -93,8 +92,13 @@
    REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3)
 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   en             !: now turbulent kinetic energy   [m2/s2]
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rn_lc0, rn_ediff0, rn_ediss0, rn_ebb0, rn_efr0
+
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_niw          !: TKE budget- near-inertial waves term
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   efr            ! surface boundary condition for nn_etau = 4
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avt_k , avm_k  ! not enhanced Kz
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avmu_k, avmv_k ! not enhanced Kz
 #if defined key_c1d
    !                                                                        !!** 1D cfg only  **   ('key_c1d')
@@ -184,5 +188,26 @@
          avmu(:,:,:) = avmu_k(:,:,:) 
          avmv(:,:,:) = avmv_k(:,:,:) 
-      ENDIF 
+      ENDIF
+      !
+      IF( nn_spp_tkelc > 0 ) THEN
+          rn_lc0 = rn_lc
+          CALL spp_gen(kt,rn_lc0,nn_spp_tkelc,rn_tkelc_sd,    jk_spp_tkelc )
+      ENDIF
+      IF( nn_spp_tkedf > 0 ) THEN
+          rn_ediff0 = rn_ediff
+          CALL spp_gen(kt,rn_ediff0,nn_spp_tkedf,rn_tkedf_sd, jk_spp_tkedf )
+      ENDIF
+      IF( nn_spp_tkeds > 0 ) THEN
+          rn_ediss0 = rn_ediss
+          CALL spp_gen(kt,rn_ediss0,nn_spp_tkeds,rn_tkeds_sd, jk_spp_tkeds )
+      ENDIF
+      IF( nn_spp_tkebb > 0 ) THEN
+          rn_ebb0 = rn_ebb
+          CALL spp_gen(kt,rn_ebb0,nn_spp_tkebb,rn_tkebb_sd,   jk_spp_tkebb )
+      ENDIF
+      IF( nn_spp_tkefr > 0 ) THEN
+          rn_efr0 = rn_efr
+          CALL spp_gen(kt,rn_efr0,nn_spp_tkefr,rn_tkefr_sd,   jk_spp_tkefr )
+      ENDIF
       !
       CALL tke_tke      ! now tke (en)
@@ -199,5 +224,9 @@
       IF( .NOT.Agrif_Root() )   CALL Agrif_Update_Tke( kt )      ! children only
 #endif      
-     ! 
+      IF (  kt .eq. nitend ) THEN
+        DEALLOCATE ( rn_lc0, rn_ediff0, rn_ediss0, rn_ebb0, rn_efr0 )
+      ENDIF
+      !
+
    END SUBROUTINE zdf_tke
 
@@ -225,6 +254,6 @@
       REAL(wp) ::   zrhoa  = 1.22                   ! Air density kg/m3
       REAL(wp) ::   zcdrag = 1.5e-3                 ! drag coefficient
-      REAL(wp) ::   zbbrau, zesh2                   ! temporary scalars
-      REAL(wp) ::   zfact1, zfact2, zfact3          !    -         -
+      REAL(wp) ::   zesh2                           ! temporary scalars
+      REAL(wp) ::   zfact1                          !    -         -
       REAL(wp) ::   ztx2  , zty2  , zcof            !    -         -
       REAL(wp) ::   ztau  , zdif                    !    -         -
@@ -233,5 +262,5 @@
 !!bfr      REAL(wp) ::   zebot                           !    -         -
       INTEGER , POINTER, DIMENSION(:,:  ) :: imlc
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zhlc
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zhlc, zbbrau,zfact2,zfact3
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw
       !!--------------------------------------------------------------------
@@ -240,11 +269,14 @@
       !
       CALL wrk_alloc( jpi,jpj, imlc )    ! integer
-      CALL wrk_alloc( jpi,jpj, zhlc ) 
+      CALL wrk_alloc( jpi,jpj, zhlc )
+      CALL wrk_alloc( jpi,jpj, zbbrau )
+      CALL wrk_alloc( jpi,jpj, zfact2 )
+      CALL wrk_alloc( jpi,jpj, zfact3 )
       CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 
       !
-      zbbrau = rn_ebb / rau0       ! Local constant initialisation
+      zbbrau = rn_ebb0 / rau0       ! Local constant initialisation
       zfact1 = -.5_wp * rdt 
-      zfact2 = 1.5_wp * rdt * rn_ediss
-      zfact3 = 0.5_wp       * rn_ediss
+      zfact2 = 1.5_wp * rdt * rn_ediss0
+      zfact3 = 0.5_wp       * rn_ediss0
       !
       !
@@ -261,5 +293,5 @@
       DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
+            en(ji,jj,1) = MAX( rn_emin0, zbbrau(ji,jj) * taum(ji,jj) ) * tmask(ji,jj,1)
          END DO
       END DO
@@ -326,8 +358,9 @@
                   !                                           ! vertical velocity due to LC
                   zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) )
-                  zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )
+                  zwlc = zind * rn_lc0(ji,jj) * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )
                   !                                           ! TKE Langmuir circulation source term
-                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( 1._wp - fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc ) /   &
+                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc ) /   &
                      &   zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1)
+
                END DO
             END DO
@@ -380,9 +413,9 @@
                zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
                zd_lw(ji,jj,jk) = zzd_lw
-               zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)
+               zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2(ji,jj) * dissl(ji,jj,jk) * tmask(ji,jj,jk)
                !
                !                                   ! right hand side in en
                en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    &
-                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) &
+                  &                                 + zfact3(ji,jj) * dissl(ji,jj,jk) * en (ji,jj,jk)  ) &
                   &                                 * wmask(ji,jj,jk)
             END DO
@@ -455,6 +488,6 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
-                     &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
+                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr0(ji,jj) * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
+                     &                                 * MAX(0.,1._wp - 2.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
                END DO
             END DO
@@ -464,6 +497,6 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                jk = nmln(ji,jj)
-               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
-                  &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
+               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr0(ji,jj) * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
+                  &                                 * MAX(0.,1._wp - 2.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
             END DO
          END DO
@@ -480,6 +513,6 @@
                   zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean 
                   zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
-                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
-                     &                        * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
+                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau(ji,jj) * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
+                     &                        * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
                END DO
             END DO
@@ -516,4 +549,7 @@
       CALL wrk_dealloc( jpi,jpj, imlc )    ! integer
       CALL wrk_dealloc( jpi,jpj, zhlc ) 
+      CALL wrk_dealloc( jpi,jpj, zbbrau ) 
+      CALL wrk_dealloc( jpi,jpj, zfact2 ) 
+      CALL wrk_dealloc( jpi,jpj, zfact3 ) 
       CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 
       !
@@ -698,5 +734,5 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zsqen = SQRT( en(ji,jj,jk) )
-               zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
+               zav   = rn_ediff0(ji,jj) * zmxlm(ji,jj,jk) * zsqen
                avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk)
                avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
@@ -704,4 +740,6 @@
             END DO
          END DO
+         IF(nn_spp_avt > 0 ) CALL spp_gen(1 ,avt(:,:,jk),nn_spp_avt,rn_avt_sd, jk_spp_avt, jk)
+         IF(nn_spp_avm > 0 ) CALL spp_gen(1 ,avm(:,:,jk),nn_spp_avm,rn_avm_sd, jk_spp_avm, jk)
       END DO
       CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
@@ -840,9 +878,11 @@
          rn_mxl0 = rmxl_min
       ENDIF
+
+      ALLOCATE( rn_lc0   (jpi,jpj) ) ; rn_lc0    = rn_lc
+      ALLOCATE( rn_ediff0(jpi,jpj) ) ; rn_ediff0 = rn_ediff
+      ALLOCATE( rn_ediss0(jpi,jpj) ) ; rn_ediss0 = rn_ediss
+      ALLOCATE( rn_ebb0  (jpi,jpj) ) ; rn_ebb0   = rn_ebb
+      ALLOCATE( rn_efr0  (jpi,jpj) ) ; rn_efr0   = rn_efr
       
-      IF( nn_etau == 2  ) THEN
-          ierr = zdf_mxl_alloc()
-          nmln(:,:) = nlb10           ! Initialization of nmln
-      ENDIF
 
       IF( nn_etau /= 0 .and. nn_htau == 2 ) THEN
@@ -950,21 +990,15 @@
               CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
               !
-              avt_k (:,:,:) = avt (:,:,:)
-              avm_k (:,:,:) = avm (:,:,:)
-              avmu_k(:,:,:) = avmu(:,:,:)
-              avmv_k(:,:,:) = avmv(:,:,:)
-              !
               DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO
            ENDIF
         ELSE                                   !* Start from rest
            en(:,:,:) = rn_emin * tmask(:,:,:)
-           DO jk = 1, jpk                           ! set the Kz to the background value
-              avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
-              avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)
-              avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)
-              avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)
-           END DO
         ENDIF
         !
+        avt_k (:,:,:) = avt (:,:,:)
+        avm_k (:,:,:) = avm (:,:,:)
+        avmu_k(:,:,:) = avmu(:,:,:)
+        avmv_k(:,:,:) = avmv(:,:,:)
+        
      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
         !                                   ! -------------------
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 11393)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 11394)
@@ -85,4 +85,5 @@
    USE lbcnfd, ONLY: isendto, nsndto, nfsloop, nfeloop ! Setup of north fold exchanges 
    USE sbc_oce, ONLY: lk_oasis
+   USE stopack
    USE stopar
    USE stopts
@@ -483,4 +484,5 @@
                             CALL dia_hsb_init   ! heat content, salt content and volume budgets
                             CALL     trd_init   ! Mixed-layer/Vorticity/Integral constraints trends
+                            CALL stopack_init   ! STOPACK scheme
       IF( lk_diaobs     ) THEN                  ! Observation & model comparison
                             CALL dia_obs_init            ! Initialize observational data
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 11393)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 11394)
@@ -118,12 +118,16 @@
          CALL lbc_lnk( tsb(:,:,:,tind), 'T', 1. )
       END DO
-
+      
+      IF( ln_stopack )   CALL stopack_pert( kstp )
                          CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
                                                       ! clem: moved here for bdy ice purpose
+
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Update stochastic parameters and random T/S fluctuations
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-       IF( ln_sto_eos ) CALL sto_par( kstp )          ! Stochastic parameters
-       IF( ln_sto_eos ) CALL sto_pts( tsn  )          ! Random T/S fluctuations
+
+      IF( ln_sto_eos )   CALL sto_par( kstp )          ! Stochastic parameters
+      IF( ln_sto_eos )   CALL sto_pts( tsn  )          ! Random T/S fluctuations
+      IF( ln_skeb  )     CALL skeb_comp( kstp )
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -149,5 +153,9 @@
       ENDIF
       IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths
-         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO
+         IF ( nn_spp_arnf .GT. 0 ) THEN
+              rn_avt_rnf0 = rn_avt_rnf
+              CALL spp_gen( kstp, rn_avt_rnf0,nn_spp_arnf,rn_arnf_sd,jk_spp_arnf )
+         ENDIF
+         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf0(:,:) * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO
       ENDIF
       IF( ln_zdfevd  )   CALL zdf_evd( kstp )         ! enhanced vertical eddy diffusivity
@@ -163,9 +171,11 @@
       IF( lrst_oce .AND. lk_zdftke )   CALL tke_rst( kstp, 'WRITE' )
       IF( lrst_oce .AND. lk_zdfgls )   CALL gls_rst( kstp, 'WRITE' )
+      IF( lrst_oce .AND. ln_stopack)   CALL stopack_rst( kstp, 'WRITE' )
       !
       !  LATERAL  PHYSICS
       !
       IF( lk_ldfslp ) THEN                            ! slope of lateral mixing
-         CALL eos( tsb, rhd, gdept_0(:,:,:) )               ! before in situ density
+         IF(ln_sto_eos ) CALL sto_pts( tsn )          ! Random T/S fluctuations
+                         CALL eos( tsb, rhd, gdept_0(:,:,:) )               ! before in situ density
          IF( ln_zps .AND. .NOT. ln_isfcav)                               &
             &            CALL zps_hde    ( kstp, jpts, tsb, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
@@ -202,5 +212,6 @@
           ! Note that the computation of vertical velocity above, hence "after" sea level
           ! is necessary to compute momentum advection for the rhs of barotropic loop:
-            CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation
+            IF(ln_sto_eos ) CALL sto_pts( tsn )                             ! Random T/S fluctuations
+                            CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation
             IF( ln_zps .AND. .NOT. ln_isfcav)                               &
                &            CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
@@ -214,11 +225,12 @@
                                   va(:,:,:) = 0.e0
           IF(  lk_asminc .AND. ln_asmiau .AND. &
-             & ln_dyninc       )  CALL dyn_asm_inc  ( kstp )   ! apply dynamics assimilation increment
-          IF( ln_neptsimp )       CALL dyn_nept_cor ( kstp )   ! subtract Neptune velocities (simplified)
-          IF( lk_bdy           )  CALL bdy_dyn3d_dmp( kstp )   ! bdy damping trends
-                                  CALL dyn_adv      ( kstp )   ! advection (vector or flux form)
-                                  CALL dyn_vor      ( kstp )   ! vorticity term including Coriolis
-                                  CALL dyn_ldf      ( kstp )   ! lateral mixing
-          IF( ln_neptsimp )       CALL dyn_nept_cor ( kstp )   ! add Neptune velocities (simplified)
+             & ln_dyninc       )  CALL dyn_asm_inc   ( kstp )   ! apply dynamics assimilation increment
+          IF( ln_sppt_dyn )       CALL dyn_sppt_apply( kstp, 0 )
+          IF( ln_neptsimp )       CALL dyn_nept_cor  ( kstp )   ! subtract Neptune velocities (simplified)
+          IF( lk_bdy           )  CALL bdy_dyn3d_dmp ( kstp )   ! bdy damping trends
+                                  CALL dyn_adv       ( kstp )   ! advection (vector or flux form)
+                                  CALL dyn_vor       ( kstp )   ! vorticity term including Coriolis
+                                  CALL dyn_ldf       ( kstp )   ! lateral mixing
+          IF( ln_neptsimp )       CALL dyn_nept_cor  ( kstp )   ! add Neptune velocities (simplified)
 #if defined key_agrif
           IF(.NOT. Agrif_Root())  CALL Agrif_Sponge_dyn        ! momentum sponge
@@ -262,6 +274,7 @@
                              tsa(:,:,:,:) = 0.e0            ! set tracer trends to zero
 
-      IF(  lk_asminc .AND. ln_asmiau .AND. &
-         & ln_trainc     )   CALL tra_asm_inc( kstp )       ! apply tracer assimilation increment
+      IF( ln_sppt_tra )      CALL tra_sppt_apply ( kstp, 0 )
+      IF( lk_asminc .AND. ln_asmiau .AND. &
+        & ln_trainc      )   CALL tra_asm_inc( kstp )       ! apply tracer assimilation increment
                              CALL tra_sbc    ( kstp )       ! surface boundary condition
       IF( ln_traqsr      )   CALL tra_qsr    ( kstp )       ! penetrative solar radiation qsr
@@ -280,4 +293,5 @@
       IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra          ! tracers sponge
 #endif
+      IF( ln_sppt_tra )      CALL tra_sppt_apply ( kstp, 1 )
                              CALL tra_zdf    ( kstp )       ! vertical mixing and after tracer fields
 
@@ -285,4 +299,5 @@
          IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
                              CALL tra_nxt( kstp )                ! tracer fields at next time step
+         IF( ln_sppt_tra )   CALL tra_sppt_apply ( kstp, 2 )
                              CALL eos    ( tsa, rhd, rhop, fsdept_n(:,:,:) )  ! Time-filtered in situ density for hpg computation
             IF( ln_zps .AND. .NOT. ln_isfcav)                                &
@@ -308,4 +323,15 @@
                              CALL tra_nxt( kstp )                ! tracer fields at next time step
          IF( ln_bias )       CALL dyn_bias( kstp )
+         IF( ln_sppt_tra )   THEN
+                                CALL tra_sppt_apply ( kstp, 2 )
+                                CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) )  ! now in situ density for hpg computation
+                                IF( ln_zps .AND. .NOT. ln_isfcav)                &
+                                & CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
+                                &                   rhd, gru , grv    )  ! of t, s, rd at the last ocean level
+                                IF( ln_zps .AND.       ln_isfcav)                                   & 
+                                & CALL zps_hde_isf( kstp, jpts, tsn, 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 last ocean level
+         ENDIF
       ENDIF
 
@@ -324,4 +350,5 @@
 
                                CALL dyn_bfr( kstp )         ! bottom friction
+            IF( ln_sppt_dyn )  CALL dyn_sppt_apply ( kstp, 1 )
                                CALL dyn_zdf( kstp )         ! vertical diffusion
       ELSE
@@ -331,4 +358,5 @@
         IF(  lk_asminc .AND. ln_asmiau .AND. &
            & ln_dyninc      )  CALL dyn_asm_inc( kstp )     ! apply dynamics assimilation increment
+        IF( ln_sppt_dyn )      CALL dyn_sppt_apply ( kstp, 0 )
         IF( ln_bkgwri )        CALL asm_bkg_wri( kstp )     ! output background fields
         IF( ln_neptsimp )      CALL dyn_nept_cor( kstp )    ! subtract Neptune velocities (simplified)
@@ -343,8 +371,11 @@
                                CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure
                                CALL dyn_bfr( kstp )         ! bottom friction
+            IF( ln_sppt_dyn )  CALL dyn_sppt_apply ( kstp, 1 )
                                CALL dyn_zdf( kstp )         ! vertical diffusion
                                CALL dyn_spg( kstp, indic )  ! surface pressure gradient
       ENDIF
                                CALL dyn_nxt( kstp )         ! lateral velocity at next time step
+            IF( ln_sppt_dyn )  CALL dyn_sppt_apply ( kstp, 2 )
+            IF( ln_skeb )      CALL skeb_apply ( kstp )
 
                                CALL ssh_swp( kstp )         ! swap of sea surface height
@@ -381,4 +412,5 @@
          IF( lwm.AND.numoni /= -1 ) CALL FLUSH    ( numoni )     ! flush output namelist ice
       ENDIF
+      IF( ln_sto_eos       )   CALL sto_rst_write( kstp )   ! write restart file for stochastic parameters
 
 
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 11393)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 11394)
@@ -96,4 +96,5 @@
    USE diahsb           ! heat, salt and volume budgets    (dia_hsb routine)
    USE diaharm
+   USE diaprod          ! ocean model: product diagnostics
    USE flo_oce          ! floats variables
    USE floats           ! floats computation               (flo_stp routine)
@@ -113,4 +114,6 @@
    USE timing           ! Timing
 
+   USE stopack          ! Stochastic physics
+
 #if defined key_agrif
    USE agrif_opa_sponge ! Momemtum and tracers sponges
