Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90	(revision 11442)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90	(revision 13355)
@@ -12,5 +12,5 @@
 
    !!----------------------------------------------------------------------
-   !!   tra_bbc      : update the tracer trend at ocean bottom 
+   !!   tra_bbc      : update the tracer trend at ocean bottom
    !!   tra_bbc_init : initialization of geothermal heat flux trend
    !!----------------------------------------------------------------------
@@ -19,5 +19,5 @@
    USE phycst          ! physical constants
    USE trd_oce         ! trends: ocean variables
-   USE trdtra          ! trends manager: tracers 
+   USE trdtra          ! trends manager: tracers
    USE in_out_manager  ! I/O manager
    USE iom             ! I/O manager
@@ -44,5 +44,5 @@
    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   qgh_trd1   ! geothermal heating trend
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qgh              ! structure of input qgh (file informations, fields read)
- 
+
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -58,9 +58,9 @@
       !!                  ***  ROUTINE tra_bbc  ***
       !!
-      !! ** Purpose :   Compute the bottom boundary contition on temperature 
-      !!              associated with geothermal heating and add it to the 
+      !! ** Purpose :   Compute the bottom boundary contition on temperature
+      !!              associated with geothermal heating and add it to the
       !!              general trend of temperature equations.
       !!
-      !! ** Method  :   The geothermal heat flux set to its constant value of 
+      !! ** Method  :   The geothermal heat flux set to its constant value of
       !!              86.4 mW/m2 (Stein and Stein 1992, Huang 1999).
       !!       The temperature trend associated to this heat flux through the
@@ -92,8 +92,16 @@
       !                             !  Add the geothermal heat flux trend on temperature
 
+#if defined key_traldf_c2d || key_traldf_c3d
       IF( ln_stopack .AND. nn_spp_geot > 0) THEN
           qgh_trd1(:,:) = qgh_trd0(:,:)
           CALL spp_gen(kt, qgh_trd1, nn_spp_geot, rn_geot_sd, jk_spp_geot)
       ENDIF
+#else
+      IF ( ln_stopack .AND. nn_spp_geot > 0 ) &
+         & CALL ctl_stop( 'tra_bbc: parameter perturbation will only work with '// &
+                          'key_traldf_c2d or key_traldf_c3d')
+#endif
+
+
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
@@ -144,5 +152,5 @@
       CHARACTER(len=256) ::   cn_dir    ! Root directory for location of ssr files
       !
-      NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst, sn_qgh, cn_dir 
+      NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst, sn_qgh, cn_dir
       !!----------------------------------------------------------------------
 
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 11442)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 13355)
@@ -32,5 +32,5 @@
    USE trdtra         ! trends: active tracers
    !
-   USE iom            ! IOM library               
+   USE iom            ! IOM library
    USE in_out_manager ! I/O manager
    USE lbclnk         ! ocean lateral boundary conditions
@@ -38,5 +38,5 @@
    USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
    USE stopack
 
@@ -199,11 +199,26 @@
       !
       ahu_bbl_1(:,:) = ahu_bbl(:,:)
+#if defined key_traldf_c2d || key_traldf_c3d
       IF( ln_stopack .AND. nn_spp_ahubbl > 0 ) THEN
           CALL spp_gen(1, ahu_bbl_1, nn_spp_ahubbl, rn_ahubbl_sd, jk_spp_ahubbl )
       ENDIF
+#else
+      IF ( ln_stopack .AND. nn_spp_ahubbl > 0 ) &
+         & CALL ctl_stop( 'tra_bbl_dif: parameter perturbation will only work with '// &
+                          'key_traldf_c2d or key_traldf_c3d')
+#endif
+
+
       ahv_bbl_1(:,:) = ahv_bbl(:,:)
+#if defined key_traldf_c2d || key_traldf_c3d
       IF( ln_stopack .AND. nn_spp_ahvbbl > 0 ) THEN
           CALL spp_gen(1, ahv_bbl_1, nn_spp_ahvbbl, rn_ahvbbl_sd, jk_spp_ahvbbl )
       ENDIF
+#else
+      IF ( ln_stopack .AND. nn_spp_ahvbbl > 0 ) &
+         & CALL ctl_stop( 'tra_bbl_dif: parameter perturbation will only work with '// &
+                          'key_traldf_c2d or key_traldf_c3d')
+#endif
+
       !
       DO jn = 1, kjpt                                     ! tracer loop
@@ -215,5 +230,5 @@
             END DO
          END DO
-         !               
+         !
          DO jj = 2, jpjm1                                    ! Compute the trend
             DO ji = 2, jpim1
@@ -426,5 +441,5 @@
                   za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
                   zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
-                  !                                                          ! 2*masked bottom density gradient 
+                  !                                                          ! 2*masked bottom density gradient
                   zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
                             - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
@@ -585,5 +600,5 @@
                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, &
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 11442)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 13355)
@@ -9,6 +9,6 @@
    !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
    !!             -   !  2005-11  (G. Madec) zco, zps, sco coordinate
-   !!            3.2  !  2009-04  (G. Madec & NEMO team) 
-   !!            3.4  !  2012-05  (C. Rousset) store attenuation coef for use in ice model 
+   !!            3.2  !  2009-04  (G. Madec & NEMO team)
+   !!            3.4  !  2012-05  (C. Rousset) store attenuation coef for use in ice model
    !!            3.6  !  2015-12  (O. Aumont, J. Jouanno, C. Ethe) use vertical profile of chlorophyll
    !!----------------------------------------------------------------------
@@ -43,5 +43,5 @@
    !                                 !!* Namelist namtra_qsr: penetrative solar radiation
    LOGICAL , PUBLIC ::   ln_traqsr    !: light absorption (qsr) flag
-   LOGICAL , PUBLIC ::   ln_qsr_rgb   !: Red-Green-Blue light absorption flag  
+   LOGICAL , PUBLIC ::   ln_qsr_rgb   !: Red-Green-Blue light absorption flag
    LOGICAL , PUBLIC ::   ln_qsr_2bd   !: 2 band         light absorption flag
    LOGICAL , PUBLIC ::   ln_qsr_bio   !: bio-model      light absorption flag
@@ -51,5 +51,5 @@
    REAL(wp), PUBLIC ::   rn_si0       !: very near surface depth of extinction      (RGB & 2 bands)
    REAL(wp), PUBLIC ::   rn_si1       !: deepest depth of extinction (water type I)       (2 bands)
- 
+
    ! Module variables
    REAL(wp), ALLOCATABLE ::   xsi0r(:,:)         !: inverse of rn_si0
@@ -80,5 +80,5 @@
       !!      Considering the 2 wavebands case:
       !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
-      !!         The temperature trend associated with the solar radiation penetration 
+      !!         The temperature trend associated with the solar radiation penetration
       !!         is given by : zta = 1/e3t dk[ I ] / (rau0*Cp)
       !!         At the bottom, boudary condition for the radiation is no flux :
@@ -86,5 +86,5 @@
       !!      in the last ocean level.
       !!         In z-coordinate case, the computation is only done down to the
-      !!      level where I(k) < 1.e-15 W/m2. In addition, the coefficients 
+      !!      level where I(k) < 1.e-15 W/m2. In addition, the coefficients
       !!      used for the computation are calculated one for once as they
       !!      depends on k only.
@@ -106,5 +106,5 @@
       REAL(wp) ::   zz0, zz1, z1_e3t     !    -         -
       REAL(wp) ::   zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze
-      REAL(wp) ::   zlogc, zlogc2, zlogc3 
+      REAL(wp) ::   zlogc, zlogc2, zlogc3
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zekb, zekg, zekr
       REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt, zchl3d
@@ -113,6 +113,6 @@
       IF( nn_timing == 1 )  CALL timing_start('tra_qsr')
       !
-      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        ) 
-      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 
+      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        )
+      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d )
       !
       IF( kt == nit000 ) THEN
@@ -124,5 +124,5 @@
 
       IF( l_trdtra ) THEN      ! Save ta and sa trends
-         CALL wrk_alloc( jpi, jpj, jpk, ztrdt ) 
+         CALL wrk_alloc( jpi, jpj, jpk, ztrdt )
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
       ENDIF
@@ -150,5 +150,5 @@
       !                                        Compute now qsr tracer content field
       !                                        ************************************
-      
+
       !                                           ! ============================================== !
       IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN      !  bio-model fluxes  : all vertical coordinates  !
@@ -159,5 +159,5 @@
          !                                        Add to the general trend
          DO jk = 1, jpkm1
-            DO jj = 2, jpjm1 
+            DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   z1_e3t = zfact / fse3t(ji,jj,jk)
@@ -180,13 +180,20 @@
          ENDIF
          !                                        ! ============================================== !
-      ELSE                                        !  Ocean alone : 
+      ELSE                                        !  Ocean alone :
          !                                        ! ============================================== !
          !
-         ! 
+         !
+#if defined key_traldf_c2d || key_traldf_c3d
          IF( ln_stopack .AND. ( 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
+            xsi0r = rn_si0
+            CALL spp_gen(kt, xsi0r, nn_spp_qsi0, rn_qsi0_sd, jk_spp_qsi0 )
+            xsi0r = 1.e0 / xsi0r
+         ENDIF
+#else
+         IF ( ln_stopack .AND. nn_spp_qsi0 > 0 ) &
+            & CALL ctl_stop( 'tra_qsr: parameter perturbation will only work with '// &
+                             'key_traldf_c2d or key_traldf_c3d')
+#endif
+
          !                                               ! ------------------------- !
          IF( ln_qsr_rgb) THEN                             !  R-G-B  light penetration !
@@ -199,5 +206,5 @@
                   CALL fld_read( kt, 1, sf_chl )            ! Read Chl data and provides it at the current time step
                   DO jk = 1, nksr + 1
-                     zchl3d(:,:,jk) = sf_chl(1)%fnow(:,:,1) 
+                     zchl3d(:,:,jk) = sf_chl(1)%fnow(:,:,1)
                   ENDDO
                   !
@@ -220,5 +227,5 @@
                         zpsimax = 0.6   - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3
                         zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2
-                        zCze    = 1.12  * (zchl)**0.803 
+                        zCze    = 1.12  * (zchl)**0.803
                         DO jk = 1, nksr + 1
                            zpsi = fsdept(ji,jj,jk) / zze
@@ -231,5 +238,5 @@
                ELSE                              !* Variable ocean volume but constant chrlorophyll
                   DO jk = 1, nksr + 1
-                     zchl3d(:,:,jk) = 0.05 
+                     zchl3d(:,:,jk) = 0.05
                   ENDDO
                ENDIF
@@ -256,5 +263,5 @@
 !CDIR NOVERRCHK
                   DO jj = 1, jpj
-!CDIR NOVERRCHK   
+!CDIR NOVERRCHK
                      DO ji = 1, jpi
                         zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r(ji,jj) )
@@ -289,5 +296,5 @@
                      END DO
                   END DO
-                  ! 
+                  !
                   DO jj = 1, jpj
                      DO ji = 1, jpi
@@ -296,5 +303,5 @@
                         zc2 = zcoef  * EXP( - fse3t(ji,jj,1) * zekg(ji,jj) )
                         zc3 = zcoef  * EXP( - fse3t(ji,jj,1) * zekr(ji,jj) )
-                        fraqsr_1lev(ji,jj) = 1.0 - ( zc0 + zc1 + zc2  + zc3  ) * tmask(ji,jj,2) 
+                        fraqsr_1lev(ji,jj) = 1.0 - ( zc0 + zc1 + zc2  + zc3  ) * tmask(ji,jj,2)
                      END DO
                   END DO
@@ -321,10 +328,10 @@
                zz0   =        rn_abs   * r1_rau0_rcp
                zz1   = ( 1. - rn_abs ) * r1_rau0_rcp
-               DO jk = 1, nksr                    ! solar heat absorbed at T-point in the top 400m 
+               DO jk = 1, nksr                    ! solar heat absorbed at T-point in the top 400m
                   DO jj = 1, jpj
                      DO ji = 1, jpi
                         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) ) 
+                        qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) )
                      END DO
                   END DO
@@ -344,5 +351,5 @@
                   DO jj = 2, jpjm1
                      DO ji = fs_2, fs_jpim1   ! vector opt.
-                        ! (ISF) no light penetration below the ice shelves         
+                        ! (ISF) no light penetration below the ice shelves
                         qsr_hc(ji,jj,jk) =  etot3(ji,jj,jk) * qsr(ji,jj) * tmask(ji,jj,1)
                      END DO
@@ -360,5 +367,5 @@
          !                                        Add to the general trend
          DO jk = 1, nksr
-            DO jj = 2, jpjm1 
+            DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   z1_e3t = zfact / fse3t(ji,jj,jk)
@@ -376,7 +383,7 @@
             &                    'at it= ', kt,' date= ', ndastp
          IF(lwp) WRITE(numout,*) '~~~~'
-         IF(nn_timing == 2)  CALL timing_start('iom_rstput') 
+         IF(nn_timing == 2)  CALL timing_start('iom_rstput')
          CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc      )
-         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev )   ! default definition in sbcssm 
+         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev )   ! default definition in sbcssm
          IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
          !
@@ -386,11 +393,11 @@
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt )
-         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt ) 
+         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt )
       ENDIF
       !                       ! print mean trends (used for debugging)
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
       !
-      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        ) 
-      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 
+      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        )
+      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d )
       !
       IF( nn_timing == 1 )  CALL timing_stop('tra_qsr')
@@ -408,5 +415,5 @@
       !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
       !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
-      !!      default values correspond to clear water (type I in Jerlov' 
+      !!      default values correspond to clear water (type I in Jerlov'
       !!      (1968) classification.
       !!         called by tra_qsr at the first timestep (nit000)
@@ -435,6 +442,6 @@
       IF( nn_timing == 1 )  CALL timing_start('tra_qsr_init')
       !
-      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        ) 
-      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 
+      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        )
+      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )
       !
 
@@ -465,5 +472,5 @@
 
       IF( ln_traqsr ) THEN     ! control consistency
-         !                      
+         !
          IF( .NOT.lk_qsr_bio .AND. ln_qsr_bio )   THEN
             CALL ctl_warn( 'No bio model : force ln_qsr_bio = FALSE ' )
@@ -480,5 +487,5 @@
             &              ' 2 bands, 3 RGB bands or bio-model light penetration' )
          !
-         IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr =  1 
+         IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr =  1
          IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr =  2
          IF( ln_qsr_rgb .AND. nn_chldta == 2 )   nqsr =  3
@@ -497,5 +504,5 @@
       ENDIF
       !                          ! ===================================== !
-      IF( ln_traqsr  ) THEN      !  Initialisation of Light Penetration  !  
+      IF( ln_traqsr  ) THEN      !  Initialisation of Light Penetration  !
          !                       ! ===================================== !
          !
@@ -539,5 +546,5 @@
                   zchl = 0.05                                 ! constant chlorophyll
                   irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
-                  zekb(:,:) = rkrgb(1,irgb)                   ! Separation in R-G-B depending of the chlorophyll 
+                  zekb(:,:) = rkrgb(1,irgb)                   ! Separation in R-G-B depending of the chlorophyll
                   zekg(:,:) = rkrgb(2,irgb)
                   zekr(:,:) = rkrgb(3,irgb)
@@ -546,12 +553,12 @@
                   ze0(:,:,1) = rn_abs
                   ze1(:,:,1) = zcoef
-                  ze2(:,:,1) = zcoef 
+                  ze2(:,:,1) = zcoef
                   ze3(:,:,1) = zcoef
                   zea(:,:,1) = tmask(:,:,1)                   ! = ( ze0+ze1+z2+ze3 ) * tmask
-               
+
                   DO jk = 2, nksr+1
 !CDIR NOVERRCHK
                      DO jj = 1, jpj
-!CDIR NOVERRCHK   
+!CDIR NOVERRCHK
                         DO ji = 1, jpi
                            zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r(ji,jj) )
@@ -566,5 +573,5 @@
                         END DO
                      END DO
-                  END DO 
+                  END DO
                   !
                   DO jk = 1, nksr
@@ -598,5 +605,5 @@
                         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) 
+                        etot3(ji,jj,jk) = (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) * tmask(ji,jj,1)
                      END DO
                   END DO
@@ -607,5 +614,5 @@
          ENDIF
          !                       ! ===================================== !
-      ELSE                       !        No light penetration           !                   
+      ELSE                       !        No light penetration           !
          !                       ! ===================================== !
          IF(lwp) THEN
@@ -625,6 +632,6 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        ) 
-      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 
+      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        )
+      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )
       !
       IF( nn_timing == 1 )  CALL timing_stop('tra_qsr_init')
