Index: branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90
===================================================================
--- branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 13190)
+++ branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 13191)
@@ -21,5 +21,5 @@
    USE lib_mpp        ! MPP library
    USE wrk_nemo       ! work arrays
-   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
 
@@ -31,5 +31,5 @@
 
    INTEGER  ::   albd_init = 0      !: control flag for initialization
-  
+
    REAL(wp) ::   rmue     = 0.40    !  cosine of local solar altitude
    REAL(wp) ::   ralb_oce = 0.066   ! ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001)
@@ -37,5 +37,5 @@
    REAL(wp) ::   c2       = 0.10    !  "        "
    REAL(wp) ::   rcloud   = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0)
- 
+
    !                             !!* namelist namsbc_alb
    INTEGER  ::   nn_ice_alb
@@ -52,7 +52,7 @@
       !!----------------------------------------------------------------------
       !!               ***  ROUTINE albedo_ice  ***
-      !!          
-      !! ** Purpose :   Computation of the albedo of the snow/ice system 
-      !!       
+      !!
+      !! ** Purpose :   Computation of the albedo of the snow/ice system
+      !!
       !! ** Method  :   Two schemes are available (from namelist parameter nn_ice_alb)
       !!                  0: the scheme is that of Shine & Henderson-Sellers (JGR 1985) for clear-skies
@@ -73,12 +73,12 @@
       !! ** Note    :   The parameterization from Shine & Henderson-Sellers presents several misconstructions:
       !!                  1) ice albedo when ice thick. tends to 0 is different than ocean albedo
-      !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger 
+      !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger
       !!                     under melting conditions than under freezing conditions
-      !!                  3) the evolution of ice albedo as a function of ice thickness shows  
+      !!                  3) the evolution of ice albedo as a function of ice thickness shows
       !!                     3 sharp inflexion points (at 5cm, 100cm and 150cm) that look highly unrealistic
       !!
       !! References :   Shine & Henderson-Sellers 1985, JGR, 90(D1), 2243-2250.
       !!                Brandt et al. 2005, J. Climate, vol 18
-      !!                Grenfell & Perovich 2004, JGR, vol 109 
+      !!                Grenfell & Perovich 2004, JGR, vol 109
       !!----------------------------------------------------------------------
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice      !  ice surface temperature (Kelvin)
@@ -97,10 +97,10 @@
 
       ijpl = SIZE( pt_ice, 3 )                     ! number of ice categories
-      
+
       CALL wrk_alloc( jpi,jpj,ijpl, zalb, zalb_it )
 
-      IF( albd_init == 0 )   CALL albedo_init      ! initialization 
-
-      
+      IF( albd_init == 0 )   CALL albedo_init      ! initialization
+
+
       SELECT CASE ( nn_ice_alb )
 
@@ -109,15 +109,15 @@
       !------------------------------------------
       CASE( 0 )
-       
+
          ralb_sf = 0.80       ! dry snow
          ralb_sm = 0.65       ! melting snow
          ralb_if = 0.72       ! bare frozen ice
-         ralb_im = rn_albice  ! bare puddled ice 
-         
+         ralb_im = rn_albice  ! bare puddled ice
+
          !  Computation of ice albedo (free of snow)
          WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = ralb_im
          ELSE WHERE                                              ;   zalb(:,:,:) = ralb_if
          END  WHERE
-      
+
          WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
          ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 )
@@ -127,5 +127,5 @@
          ELSE WHERE                                       ;  zalb_it = 0.1    + 3.6    * ph_ice
          END WHERE
-     
+
          DO jl = 1, ijpl
             DO jj = 1, jpj
@@ -133,9 +133,9 @@
                   ! freezing snow
                   ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2
-                  !                                        !  freezing snow        
+                  !                                        !  freezing snow
                   zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) )
                   zalb_sf   = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  &
                      &                           + ph_snw(ji,jj,jl) * ( ralb_sf - zalb_it(ji,jj,jl) ) / c1  )   &
-                     &        +         zswitch   * ralb_sf  
+                     &        +         zswitch   * ralb_sf
 
                   ! melting snow
@@ -143,10 +143,10 @@
                   zswitch   = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) )
                   zalb_sm = ( 1._wp - zswitch ) * ( ralb_im + ph_snw(ji,jj,jl) * ( ralb_sm - ralb_im ) / c2 )   &
-                      &     +         zswitch   *   ralb_sm 
+                      &     +         zswitch   *   ralb_sm
                   !
                   ! snow albedo
-                  zswitch  =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
+                  zswitch  =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )
                   zalb_st  =  zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
-               
+
                   ! Ice/snow albedo
                   zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
@@ -155,8 +155,13 @@
                END DO
             END DO
-            
+
+#if defined key_traldf_c2d || key_traldf_c3d
             IF ( ln_stopack .AND. nn_spp_icealb > 0 ) &
                & CALL spp_gen( 1, pa_ice_cs(:,:,jl), nn_spp_icealb, rn_icealb_sd, jk_spp_alb, jl )
-                        
+#else
+            IF ( ln_stopack .AND. nn_spp_icealb > 0 ) &
+               & CALL ctl_stop( 'albedo_ice: parameter perturbation will only work with '// &
+                                'key_traldf_c2d or key_traldf_c3d')
+#endif
          END DO
 
@@ -166,5 +171,5 @@
       !  New parameterization (2016)
       !------------------------------------------
-      CASE( 1 ) 
+      CASE( 1 )
 
          ralb_im = rn_albice  ! bare puddled ice
@@ -181,12 +186,12 @@
 !         ralb_sm = 0.82      ! melting snow
 !         ralb_if = 0.54      ! bare frozen ice
-! 
+!
          !  Computation of ice albedo (free of snow)
-         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
+         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) )
          z1_c2 = 1. / 0.05
          WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb = ralb_im
          ELSE WHERE                                              ;   zalb = ralb_if
          END  WHERE
-         
+
          WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
          ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = zalb     + ( 0.18 - zalb     ) * z1_c1 *  &
@@ -205,8 +210,8 @@
 
                    ! snow albedo
-                  zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
+                  zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )
                   zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
 
-                  ! Ice/snow albedo   
+                  ! Ice/snow albedo
                   zswitch             = MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
                   pa_ice_os(ji,jj,jl) = ( 1._wp - zswitch ) * zalb_st + zswitch *  zalb_it(ji,jj,jl)
@@ -214,14 +219,19 @@
               END DO
             END DO
-            
+
+#if defined key_traldf_c2d || key_traldf_c3d
             IF ( ln_stopack .AND. nn_spp_icealb > 0 ) &
                & CALL spp_gen( 1, pa_ice_os(:,:,jl), nn_spp_icealb, rn_icealb_sd, jk_spp_alb, jl )
-            
+#else
+            IF ( ln_stopack .AND. nn_spp_icealb > 0 ) &
+               & CALL ctl_stop( 'albedo_ice: parameter perturbation will only work with '// &
+                                'key_traldf_c2d or key_traldf_c3d')
+#endif
          END DO
          ! Effect of the clouds (2d order polynomial)
-         pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 ); 
+         pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 );
 
       END SELECT
-      
+
       CALL wrk_dealloc( jpi,jpj,ijpl, zalb, zalb_it )
       !
@@ -232,5 +242,5 @@
       !!----------------------------------------------------------------------
       !!               ***  ROUTINE albedo_oce  ***
-      !! 
+      !!
       !! ** Purpose :   Computation of the albedo of the ocean
       !!----------------------------------------------------------------------
@@ -238,9 +248,9 @@
       REAL(wp), DIMENSION(:,:), INTENT(out) ::   pa_oce_cs   !  albedo of ocean under clear sky
       !!
-      REAL(wp) :: zcoef 
+      REAL(wp) :: zcoef
       !!----------------------------------------------------------------------
       !
       zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )   ! Parameterization of Briegled and Ramanathan, 1982
-      pa_oce_cs(:,:) = zcoef 
+      pa_oce_cs(:,:) = zcoef
       pa_oce_os(:,:) = 0.06                       ! Parameterization of Kondratyev, 1969 and Payne, 1972
       !
@@ -257,5 +267,5 @@
       !!----------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice 
+      NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice
       !!----------------------------------------------------------------------
       !
Index: branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 13190)
+++ branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 13191)
@@ -161,5 +161,5 @@
          !                                      ! ====================== !
          !
-         rn_sfac = 1._wp       ! Default to one if missing from namelist 
+         rn_sfac = 1._wp       ! Default to one if missing from namelist
          REWIND( numnam_ref )              ! Namelist namsbc_core in reference namelist : CORE bulk parameters
          READ  ( numnam_ref, namsbc_core, IOSTAT = ios, ERR = 901)
@@ -205,8 +205,14 @@
       ENDIF
 
+#if defined key_traldf_c2d || key_traldf_c3d
       IF( ln_stopack .AND. 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
+#else
+      IF ( ln_stopack .AND. nn_spp_relw > 0 ) &
+         & CALL ctl_stop( 'sbc_blk_core: parameter perturbation will only work with '// &
+                          'key_traldf_c2d or key_traldf_c3d')
+#endif
 
       CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step
@@ -217,9 +223,9 @@
 #if defined key_cice
       IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN
-         qlw_ice(:,:,1)   = sf(jp_qlw)%fnow(:,:,1) 
+         qlw_ice(:,:,1)   = sf(jp_qlw)%fnow(:,:,1)
          IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
          ELSE                ; qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1)
          ENDIF
-         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)         
+         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)
          qatm_ice(:,:)    = sf(jp_humi)%fnow(:,:,1)
          tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac
@@ -231,6 +237,6 @@
       !
    END SUBROUTINE sbc_blk_core
-   
-   
+
+
    SUBROUTINE blk_oce_core( kt, sf, pst, pu, pv )
       !!---------------------------------------------------------------------
@@ -242,5 +248,5 @@
       !! ** Method  :   CORE bulk formulea for the ocean using atmospheric
       !!      fields read in sbc_read
-      !! 
+      !!
       !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2)
       !!              - vtau    : j-component of the stress at V-point  (N/m2)
@@ -280,5 +286,5 @@
       ! local scalars ( place there for vector optimisation purposes)
       zcoef_qsatw = 0.98 * 640380. / rhoa
-      
+
       zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
 
@@ -288,5 +294,5 @@
 
       ! ... components ( U10m - U_oce ) at T-point (unmasked)
-      zwnd_i(:,:) = 0.e0  
+      zwnd_i(:,:) = 0.e0
       zwnd_j(:,:) = 0.e0
 #if defined key_cyclone
@@ -332,5 +338,5 @@
       CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, sf(jp_humi)%fnow, wndm,   &
          &               Cd, Ch, Ce, zt_zu, zq_zu )
-    
+
       ! ... tau module, i and j component
       DO jj = 1, jpj
@@ -363,5 +369,5 @@
       CALL lbc_lnk( vtau(:,:), 'V', -1. )
 
-    
+
       !  Turbulent fluxes over ocean
       ! -----------------------------
@@ -388,5 +394,5 @@
          CALL prt_ctl( tab2d_1=zst   , clinfo1=' blk_oce_core: zst    : ')
       ENDIF
-       
+
       ! ----------------------------------------------------------------------------- !
       !     III    Total FLUXES                                                       !
@@ -396,5 +402,5 @@
          &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1)
       !
-      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar 
+      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar
          &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus                         &   ! remove latent melting heat for solid precip
          &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST
@@ -437,6 +443,6 @@
       !
    END SUBROUTINE blk_oce_core
- 
-   
+
+
 #if defined key_lim2 || defined key_lim3
    SUBROUTINE blk_ice_core_tau
@@ -531,5 +537,5 @@
 
       IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_tau')
-      
+
    END SUBROUTINE blk_ice_core_tau
 
@@ -544,5 +550,5 @@
       !!                between atmosphere and sea-ice using CORE bulk
       !!                formulea, ice variables and read atmmospheric fields.
-      !! 
+      !!
       !! caution : the net upward water flux has with mm/day unit
       !!---------------------------------------------------------------------
@@ -564,5 +570,5 @@
       IF( nn_timing == 1 )  CALL timing_start('blk_ice_core_flx')
       !
-      CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 
+      CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb )
 
       ! local scalars ( place there for vector optimisation purposes)
@@ -587,5 +593,5 @@
                z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
                ! lw sensitivity
-               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3                                               
+               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
 
                ! ----------------------------!
@@ -597,5 +603,5 @@
                z_qsb(ji,jj,jl) = rhoa * cpa * Cice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) )
                ! Latent Heat
-               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * wndm_ice(ji,jj)   &                           
+               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * wndm_ice(ji,jj)   &
                   &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) )
               ! Latent heat sensitivity for ice (Dqla/Dt)
@@ -628,5 +634,5 @@
 
 #if defined  key_lim3
-      CALL wrk_alloc( jpi,jpj, zevap, zsnw ) 
+      CALL wrk_alloc( jpi,jpj, zevap, zsnw )
 
       ! --- evaporation --- !
@@ -638,5 +644,5 @@
       ! --- evaporation minus precipitation --- !
       zsnw(:,:) = 0._wp
-      CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing 
+      CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing
       emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
       emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
@@ -661,8 +667,8 @@
       DO jl = 1, jpl
          qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) )
-                                   ! But we do not have Tice => consider it at 0 degC => evap=0 
+                                   ! But we do not have Tice => consider it at 0 degC => evap=0
       END DO
 
-      CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 
+      CALL wrk_dealloc( jpi,jpj, zevap, zsnw )
 #endif
 
@@ -688,5 +694,5 @@
       !
       IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_flx')
-      
+
    END SUBROUTINE blk_ice_core_flx
 #endif
@@ -701,5 +707,5 @@
       !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
       !!
-      !! ** Method : Monin Obukhov Similarity Theory 
+      !! ** Method : Monin Obukhov Similarity Theory
       !!             + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10)
       !!
@@ -711,5 +717,5 @@
       !!    - better first guess of stability by checking air-sea difference of virtual temperature
       !!       rather than temperature difference only...
-      !!    - added function "cd_neutral_10m" that uses the improved parametrization of 
+      !!    - added function "cd_neutral_10m" that uses the improved parametrization of
       !!      Large & Yeager 2008. Drag-coefficient reduction for Cyclone conditions!
       !!    - using code-wide physical constants defined into "phycst.mod" rather than redifining them
@@ -746,5 +752,5 @@
 
       IF( nn_timing == 1 )  CALL timing_start('turb_core_2z')
-    
+
       CALL wrk_alloc( jpi,jpj, U_zu, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd )
       CALL wrk_alloc( jpi,jpj, zeta_u, stab )
@@ -758,5 +764,5 @@
       U_zu = MAX( 0.5 , dU )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s
 
-      !! First guess of stability: 
+      !! First guess of stability:
       ztmp0 = T_zt*(1. + 0.608*q_zt) - sst*(1. + 0.608*q_sat) ! air-sea difference of virtual pot. temp. at zt
       stab  = 0.5 + sign(0.5,ztmp0)                           ! stab = 1 if dTv > 0  => STABLE, 0 if unstable
@@ -772,5 +778,5 @@
       Ce_n10  = 1.e-3*( 34.6 * sqrt_Cd_n10 )
       Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))
-    
+
       !! Initializing transf. coeff. with their first guess neutral equivalents :
       Cd = ztmp0   ;   Ce = Ce_n10   ;   Ch = Ch_n10   ;   sqrt_Cd = sqrt_Cd_n10
@@ -783,14 +789,14 @@
          !
          ztmp1 = T_zu - sst   ! Updating air/sea differences
-         ztmp2 = q_zu - q_sat 
+         ztmp2 = q_zu - q_sat
 
          ! Updating turbulent scales :   (L&Y 2004 eq. (7))
          ztmp1  = Ch/sqrt_Cd*ztmp1    ! theta*
          ztmp2  = Ce/sqrt_Cd*ztmp2    ! q*
-       
+
          ztmp0 = T_zu*(1. + 0.608*q_zu) ! virtual potential temperature at zu
 
          ! Estimate the inverse of Monin-Obukov length (1/L) at height zu:
-         ztmp0 =  (vkarmn*grav/ztmp0*(ztmp1*(1.+0.608*q_zu) + 0.608*T_zu*ztmp2)) / (Cd*U_zu*U_zu) 
+         ztmp0 =  (vkarmn*grav/ztmp0*(ztmp1*(1.+0.608*q_zu) + 0.608*T_zu*ztmp2)) / (Cd*U_zu*U_zu)
          !                                                                     ( Cd*U_zu*U_zu is U*^2 at zu)
 
@@ -799,5 +805,5 @@
          zpsi_h_u = psi_h( zeta_u )
          zpsi_m_u = psi_m( zeta_u )
-       
+
          !! Shifting temperature and humidity at zu (L&Y 2004 eq. (9b-9c))
          IF ( .NOT. l_zt_equal_zu ) THEN
@@ -808,7 +814,7 @@
             q_zu = max(0., q_zu)
          END IF
-       
+
          IF( ln_cdgw ) THEN      ! surface wave case
-            sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u ) 
+            sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u )
             Cd      = sqrt_Cd * sqrt_Cd
          ELSE
@@ -820,5 +826,5 @@
            ztmp0 = cd_neutral_10m(ztmp0)                                                 ! Cd_n10
            sqrt_Cd_n10 = sqrt(ztmp0)
-       
+
            Ce_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)                     ! L&Y 2004 eq. (6b)
            stab    = 0.5 + sign(0.5,zeta_u)                           ! update stability
@@ -827,5 +833,5 @@
            !! Update of transfer coefficients:
            ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - zpsi_m_u)   ! L&Y 2004 eq. (10a)
-           Cd      = ztmp0 / ( ztmp1*ztmp1 )   
+           Cd      = ztmp0 / ( ztmp1*ztmp1 )
            sqrt_Cd = SQRT( Cd )
          ENDIF
@@ -833,8 +839,8 @@
          ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10
          ztmp2 = sqrt_Cd / sqrt_Cd_n10
-         ztmp1 = 1. + Ch_n10*ztmp0               
+         ztmp1 = 1. + Ch_n10*ztmp0
          Ch  = Ch_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10b)
          !
-         ztmp1 = 1. + Ce_n10*ztmp0               
+         ztmp1 = 1. + Ce_n10*ztmp0
          Ce  = Ce_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c)
          !
@@ -854,5 +860,5 @@
    FUNCTION cd_neutral_10m( zw10 )
       !!----------------------------------------------------------------------
-      !! Estimate of the neutral drag coefficient at 10m as a function 
+      !! Estimate of the neutral drag coefficient at 10m as a function
       !! of neutral wind  speed at 10m
       !!
@@ -860,15 +866,15 @@
       !!
       !! Author: L. Brodeau, june 2014
-      !!----------------------------------------------------------------------    
+      !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   zw10           ! scalar wind speed at 10m (m/s)
       REAL(wp), DIMENSION(jpi,jpj)             ::   cd_neutral_10m
       !
       REAL(wp), DIMENSION(:,:), POINTER ::   rgt33
-      !!----------------------------------------------------------------------    
+      !!----------------------------------------------------------------------
       !
       CALL wrk_alloc( jpi,jpj, rgt33 )
       !
       !! When wind speed > 33 m/s => Cyclone conditions => special treatment
-      rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) )   ! If zw10 < 33. => 0, else => 1  
+      rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) )   ! If zw10 < 33. => 0, else => 1
       cd_neutral_10m = 1.e-3 * ( &
          &       (1._wp - rgt33)*( 2.7_wp/zw10 + 0.142_wp + zw10/13.09_wp - 3.14807E-10*zw10**6) & ! zw10< 33.
Index: branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90
===================================================================
--- branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 13190)
+++ branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 13191)
@@ -24,5 +24,5 @@
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    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
    USE wrk_nemo       ! Memory Allocation
@@ -42,5 +42,5 @@
    REAL(wp)        ::   rn_dqdt         ! restoring factor on SST and SSS
    REAL(wp)        ::   rn_deds         ! restoring factor on SST and SSS
-   LOGICAL         ::   ln_sssr_bnd     ! flag to bound erp term 
+   LOGICAL         ::   ln_sssr_bnd     ! flag to bound erp term
    REAL(wp)        ::   rn_sssr_bnd     ! ABS(Max./Min.) value of erp term [mm/day]
 
@@ -101,6 +101,13 @@
                rn_dqdt_s(:,:) = rn_dqdt
 
-               IF( ln_stopack .AND. nn_spp_dqdt > 0 ) &
-                  & CALL spp_gen( kt, rn_dqdt_s, nn_spp_dqdt, rn_dqdt_sd, jk_spp_dqdt )
+#if defined key_traldf_c2d || key_traldf_c3d
+            IF( ln_stopack .AND. nn_spp_dqdt > 0 ) &
+               & CALL spp_gen( kt, rn_dqdt_s, nn_spp_dqdt, rn_dqdt_sd, jk_spp_dqdt )
+#else
+            IF ( ln_stopack .AND. nn_spp_dqdt > 0 ) &
+               & CALL ctl_stop( 'sbc_ssr: parameter perturbation will only work with '// &
+                                'key_traldf_c2d or key_traldf_c3d')
+#endif
+
                DO jj = 1, jpj
                   DO ji = 1, jpi
@@ -117,11 +124,19 @@
                CALL wrk_alloc( jpi, jpj, zsrp)
                zsrp(:,:) = rn_deds
+#if defined key_traldf_c2d || key_traldf_c3d
                IF( ln_stopack .AND. nn_spp_dedt > 0 ) &
                   & CALL spp_gen(kt, zsrp, nn_spp_dedt, rn_dedt_sd, jk_spp_deds )
+#else
+            IF ( ln_stopack .AND. nn_spp_icealb > 0 ) &
+               & CALL ctl_stop( 'sbc_ssr: parameter perturbation will only work with '// &
+                                'key_traldf_c2d or key_traldf_c3d')
+#endif
+
+
 !CDIR COLLAPSE
                DO jj = 1, jpj
                   DO ji = 1, jpi
                      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) ) 
+                        &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) )
                      sfx(ji,jj) = sfx(ji,jj) + zerp                 ! salt flux
                      erp(ji,jj) = zerp / MAX( sss_m(ji,jj), 1.e-20 ) ! converted into an equivalent volume flux (diagnostic only)
@@ -134,10 +149,16 @@
                CALL wrk_alloc( jpi, jpj, zsrp)
                zsrp(:,:) = rn_deds
+#if defined key_traldf_c2d || key_traldf_c3d
                IF( ln_stopack .AND. nn_spp_dedt > 0 ) &
                   & CALL spp_gen( kt, zsrp, nn_spp_dedt, rn_dedt_sd, jk_spp_deds )
-               zerp_bnd = rn_sssr_bnd / rday                          !       -              -    
+#else
+               IF ( ln_stopack .AND. nn_spp_dedt > 0 ) &
+                  & CALL ctl_stop( 'sbc_ssr: parameter perturbation will only work with '// &
+                                   'key_traldf_c2d or key_traldf_c3d')
+#endif
+               zerp_bnd = rn_sssr_bnd / rday                          !       -              -
 !CDIR COLLAPSE
                DO jj = 1, jpj
-                  DO ji = 1, jpi                            
+                  DO ji = 1, jpi
                      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) )   &
@@ -161,5 +182,5 @@
    END SUBROUTINE sbc_ssr
 
- 
+
    SUBROUTINE sbc_ssr_init
       !!---------------------------------------------------------------------
@@ -184,6 +205,6 @@
       !!----------------------------------------------------------------------
       !
- 
-      REWIND( numnam_ref )              ! Namelist namsbc_ssr in reference namelist : 
+
+      REWIND( numnam_ref )              ! Namelist namsbc_ssr in reference namelist :
       READ  ( numnam_ref, namsbc_ssr, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_ssr in reference namelist', lwp )
@@ -240,10 +261,10 @@
       ENDIF
       !
-      !                            !* Initialize qrp and erp if no restoring 
+      !                            !* Initialize qrp and erp if no restoring
       IF( nn_sstr /= 1                   )   qrp(:,:) = 0._wp
       IF( nn_sssr /= 1 .OR. nn_sssr /= 2 )   erp(:,:) = 0._wp
       !
    END SUBROUTINE sbc_ssr_init
-      
+
    !!======================================================================
 END MODULE sbcssr
