Index: branches/UKMO/dev_r5518_GO6_stophys/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_stophys/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 15290)
+++ branches/UKMO/dev_r5518_GO6_stophys/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 15294)
@@ -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
    !!----------------------------------------------------------------------
@@ -33,4 +33,5 @@
    USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
+   USE stopack
 
    IMPLICIT NONE
@@ -42,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
@@ -50,7 +51,7 @@
    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) ::   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)
@@ -79,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 :
@@ -85,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.
@@ -105,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
@@ -112,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 +125,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
@@ -153,5 +154,5 @@
       !                                        Compute now qsr tracer content field
       !                                        ************************************
-      
+
       !                                           ! ============================================== !
       IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN      !  bio-model fluxes  : all vertical coordinates  !
@@ -162,5 +163,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)
@@ -183,8 +184,21 @@
          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
+#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 !
             !                                             ! ------------------------- !
@@ -196,5 +210,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
                   !
@@ -217,5 +231,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
@@ -228,5 +242,5 @@
                ELSE                              !* Variable ocean volume but constant chrlorophyll
                   DO jk = 1, nksr + 1
-                     zchl3d(:,:,jk) = 0.05 
+                     zchl3d(:,:,jk) = 0.05
                   ENDDO
                ENDIF
@@ -253,7 +267,7 @@
 !CDIR NOVERRCHK
                   DO jj = 1, jpj
-!CDIR NOVERRCHK   
-                     DO ji = 1, jpi
-                        zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r     )
+!CDIR NOVERRCHK
+                     DO ji = 1, jpi
+                        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) )
@@ -286,12 +300,12 @@
                      END DO
                   END DO
-                  ! 
+                  !
                   DO jj = 1, jpj
                      DO ji = 1, jpi
-                        zc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r     )
+                        zc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r(ji,jj) )
                         zc1 = zcoef  * EXP( - fse3t(ji,jj,1) * zekb(ji,jj) )
                         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
@@ -314,13 +328,14 @@
             !                                             ! ------------------------- !
             !
-            IF( lk_vvl ) THEN                                  !* variable volume
+            IF( lk_vvl .OR. ( ln_stopack .AND. ( nn_spp_qsi0 > 0 ) ) ) THEN        !* variable volume
+
                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 ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
-                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + 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) ) 
+                        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
                   END DO
@@ -330,6 +345,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
@@ -340,5 +355,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
@@ -356,5 +371,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)
@@ -375,7 +390,7 @@
             IF(lflush) CALL flush(numout)
          ENDIF
-         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')
          !
@@ -385,11 +400,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')
@@ -407,5 +422,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)
@@ -434,6 +449,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 +480,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 +495,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
@@ -498,7 +513,8 @@
       ENDIF
       !                          ! ===================================== !
-      IF( ln_traqsr  ) THEN      !  Initialisation of Light Penetration  !  
+      IF( ln_traqsr  ) THEN      !  Initialisation of Light Penetration  !
          !                       ! ===================================== !
          !
+         ALLOCATE( xsi0r(jpi,jpj) )
          xsi0r = 1.e0 / rn_si0
          xsi1r = 1.e0 / rn_si1
@@ -539,5 +555,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,14 +562,14 @@
                   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     )
+                           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) )
@@ -566,5 +582,5 @@
                         END DO
                      END DO
-                  END DO 
+                  END DO
                   !
                   DO jk = 1, nksr
@@ -600,7 +616,7 @@
                   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 )
-                        etot3(ji,jj,jk) = (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) * tmask(ji,jj,1) 
+                        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
                   END DO
@@ -611,5 +627,5 @@
          ENDIF
          !                       ! ===================================== !
-      ELSE                       !        No light penetration           !                   
+      ELSE                       !        No light penetration           !
          !                       ! ===================================== !
          IF(lwp) THEN
@@ -617,5 +633,5 @@
             WRITE(numout,*) 'tra_qsr_init : NO solar flux penetration'
             WRITE(numout,*) '~~~~~~~~~~~~'
-            IF(lflush) CALL flush(numout)
+            IF(lflush) CALL flush(numout)            
          ENDIF
       ENDIF
@@ -630,6 +646,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')
