Index: branches/CMIP5_IPSL/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- branches/CMIP5_IPSL/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 1803)
+++ branches/CMIP5_IPSL/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 1829)
@@ -27,4 +27,5 @@
    USE iom             ! I/O manager
    USE fldread         ! read input fields
+   USE dtachl
 
    IMPLICIT NONE
@@ -38,10 +39,10 @@
    LOGICAL , PUBLIC ::   ln_qsr_2bd = .TRUE.    !: 2 band         light absorption flag
    LOGICAL , PUBLIC ::   ln_qsr_bio = .FALSE.   !: bio-model      light absorption flag
-   INTEGER , PUBLIC ::   nn_chldta  = 0         !: use Chlorophyll data (=1) or not (=0)
+   INTEGER , PUBLIC ::   nn_chldta  = 0         !: use Chlorophyll 2D data (=1) 3D data (=2) or not (=0)
    REAL(wp), PUBLIC ::   rn_abs     = 0.58_wp   !: fraction absorbed in the very near surface (RGB & 2 bands)
    REAL(wp), PUBLIC ::   rn_si0     = 0.35_wp   !: very near surface depth of extinction      (RGB & 2 bands)
    REAL(wp), PUBLIC ::   rn_si1     = 23.0_wp   !: deepest depth of extinction (water type I)       (2 bands)
    REAL(wp), PUBLIC ::   rn_si2     = 61.8_wp   !: deepest depth of extinction (blue & 0.01 mg.m-3)     (RGB)
-   
+
    ! Module variables
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read)
@@ -96,5 +97,6 @@
       REAL(wp) ::   zchl, zcoef, zsi0r   ! temporary scalars
       REAL(wp) ::   zc0, zc1, zc2, zc3   !    -         -
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zekb, zekg, zekr            ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zekb2, zekg2, zekr2            ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zekb3, zekg3, zekr3            ! 3D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze0, ze1 , ze2, ze3, zea    ! 3D workspace
       !!----------------------------------------------------------------------
@@ -133,9 +135,9 @@
             !                                             ! ------------------------- !
             ! Set chlorophyl concentration
-            IF( nn_chldta ==1 ) THEN                             !*  Variable Chlorophyll
+            IF( nn_chldta == 1 ) THEN                             !*  2D Variable Chlorophyll
                !
                CALL fld_read( kt, 1, sf_chl )                         ! Read Chl data and provides it at the current time step
                !         
-!CDIR COLLAPSE
+!CDIR COLLAPSE 
 !CDIR NOVERRCHK
                DO jj = 1, jpj                                         ! Separation in R-G-B depending of the surface Chl
@@ -144,7 +146,7 @@
                      zchl = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj) ) )
                      irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
-                     zekb(ji,jj) = rkrgb(1,irgb)
-                     zekg(ji,jj) = rkrgb(2,irgb)
-                     zekr(ji,jj) = rkrgb(3,irgb)
+                     zekb2(ji,jj) = rkrgb(1,irgb)
+                     zekg2(ji,jj) = rkrgb(2,irgb)
+                     zekr2(ji,jj) = rkrgb(3,irgb)
                   END DO
                END DO
@@ -161,10 +163,61 @@
 !CDIR NOVERRCHK
                   DO jj = 1, jpj
+!CDIR NOVERRCHK
+                     DO ji = 1, jpi
+                        zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zsi0r     )
+                        zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekb2(ji,jj) )
+                        zc2 = ze2(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekg2(ji,jj) )
+                        zc3 = ze3(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekr2(ji,jj) )
+                        ze0(ji,jj,jk) = zc0
+                        ze1(ji,jj,jk) = zc1
+                        ze2(ji,jj,jk) = zc2
+                        ze3(ji,jj,jk) = zc3
+                        zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk)
+                     END DO
+                  END DO
+               END DO
+               !
+               DO jk = 1, nksr                                        ! compute and add qsr trend to ta
+                  ta(:,:,jk) = ta(:,:,jk) + ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk)
+               END DO
+               zea(:,:,nksr+1:jpk) = 0.e0     ! below 400m set to zero
+               CALL iom_put( 'qsr3d', zea )   ! Shortwave Radiation 3D distribution
+               !
+            ! Set chlorophyl concentration
+            ELSE IF( nn_chldta == 2) THEN                             !*  3D Variable Chlorophyll
+               !
+               CALL dta_chl( kt )
+               !         
+               DO jk = 1, jpkm1                                         ! Separation in R-G-B depending of the surface Chl
+!CDIR NOVERRCHK
+                 DO jj = 1, jpj                                         ! Separation in R-G-B depending of the surface Chl
+!CDIR NOVERRCHK
+                   DO ji = 1, jpi
+                     zchl = MIN( 10. , MAX( 0.03, chl_dta(ji,jj,jk) ) )
+                     irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
+                     zekb3(ji,jj,jk) = rkrgb(1,irgb) * fse3t(ji,jj,jk)
+                     zekg3(ji,jj,jk) = rkrgb(2,irgb) * fse3t(ji,jj,jk)
+                     zekr3(ji,jj,jk) = rkrgb(3,irgb) * fse3t(ji,jj,jk)
+                   END DO
+                 END DO
+               ENDDO
+               !
+               zsi0r = 1.e0 / rn_si0
+               zcoef  = ( 1. - rn_abs ) / 3.e0                        ! equi-partition in R-G-B
+               ze0(:,:,1) = rn_abs  * qsr(:,:)
+               ze1(:,:,1) = zcoef * qsr(:,:)
+               ze2(:,:,1) = zcoef * qsr(:,:)
+               ze3(:,:,1) = zcoef * qsr(:,:)
+               zea(:,:,1) =         qsr(:,:)
+               !
+               DO jk = 2, nksr+1
+!CDIR NOVERRCHK
+                  DO jj = 1, jpj
 !CDIR NOVERRCHK   
                      DO ji = 1, jpi
                         zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zsi0r     )
-                        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) )
-                        zc3 = ze3(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekr(ji,jj) )
+                        zc1 = ze1(ji,jj,jk-1) * EXP( - zekb3(ji,jj,jk-1) )
+                        zc2 = ze2(ji,jj,jk-1) * EXP( - zekg3(ji,jj,jk-1) )
+                        zc3 = ze3(ji,jj,jk-1) * EXP( - zekr3(ji,jj,jk-1) )
                         ze0(ji,jj,jk) = zc0
                         ze1(ji,jj,jk) = zc1
@@ -261,13 +314,13 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namtra_qsr : set the parameter of penetration'
-         WRITE(numout,*) '      Light penetration (T) or not (F)         ln_traqsr  = ', ln_traqsr
-         WRITE(numout,*) '      RGB (Red-Green-Blue) light penetration   ln_qsr_rgb = ', ln_qsr_rgb
-         WRITE(numout,*) '      2 band               light penetration   ln_qsr_2bd = ', ln_qsr_2bd
-         WRITE(numout,*) '      bio-model            light penetration   ln_qsr_bio = ', ln_qsr_bio
-         WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)    nn_chldta  = ', nn_chldta
-         WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs = ', rn_abs
-         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0 = ', rn_si0
-         WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1 = ', rn_si1
-         WRITE(numout,*) '      3 bands: longest depth of extinction         rn_si2 = ', rn_si2
+         WRITE(numout,*) '      Light penetration (T) or not (F)                 ln_traqsr  = ', ln_traqsr
+         WRITE(numout,*) '      RGB (Red-Green-Blue) light penetration           ln_qsr_rgb = ', ln_qsr_rgb
+         WRITE(numout,*) '      2 band               light penetration           ln_qsr_2bd = ', ln_qsr_2bd
+         WRITE(numout,*) '      bio-model            light penetration           ln_qsr_bio = ', ln_qsr_bio
+         WRITE(numout,*) '      RGB : Chl 2D/3D data (=1/2) or cst value (=0)    nn_chldta  = ', nn_chldta
+         WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)        rn_abs = ', rn_abs
+         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction      rn_si0 = ', rn_si0
+         WRITE(numout,*) '      2 bands: longest depth of extinction             rn_si1 = ', rn_si1
+         WRITE(numout,*) '      3 bands: longest depth of extinction             rn_si2 = ', rn_si2
       ENDIF
 
@@ -278,4 +331,8 @@
             ln_qsr_bio = .FALSE.
          ENDIF
+         IF( .NOT.lk_dtachl .AND. ln_qsr_rgb .AND. nn_chldta == 2 ) THEN 
+            CALL ctl_stop( 'You want to use a Chl 3D data to force your light penetration', &
+            &              'key_dtachl is required in compilation '      )
+         ENDIF
          !
          ioptio = 0                      ! Parameter control
@@ -284,24 +341,20 @@
          IF( ln_qsr_bio  )   ioptio = ioptio + 1
          !
-         IF( ioptio /= 1 ) THEN
-            ln_qsr_rgb = .TRUE.
-            nn_chldta  = 0
-            ln_qsr_2bd = .FALSE.
-            ln_qsr_bio = .FALSE.
-            CALL ctl_warn( '          Choose ONE type of light penetration in namelist namtra_qsr',   &
-           &               ' otherwise, we force the model to run with RGB light penetration' )
-         ENDIF
+         IF( ioptio /= 1 )  &
+            CALL ctl_stop( '          Choose ONE type of light penetration in namelist namtra_qsr' )
          !
          IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr =  1 
          IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr =  2
-         IF( ln_qsr_2bd                      )   nqsr =  3
-         IF( ln_qsr_bio                      )   nqsr =  4
+         IF( ln_qsr_rgb .AND. nn_chldta == 2 )   nqsr =  3
+         IF( ln_qsr_2bd                      )   nqsr =  4 
+         IF( ln_qsr_bio                      )   nqsr =  5
          !
          IF(lwp) THEN                   ! Print the choice
             WRITE(numout,*)
             IF( nqsr ==  1 )   WRITE(numout,*) '         R-G-B  light penetration - Constant Chlorophyll'
-            IF( nqsr ==  2 )   WRITE(numout,*) '         R-G-B  light penetration - Chl data '
-            IF( nqsr ==  3 )   WRITE(numout,*) '         2 band light penetration'
-            IF( nqsr ==  4 )   WRITE(numout,*) '         bio-model light penetration'
+            IF( nqsr ==  2 )   WRITE(numout,*) '         R-G-B  light penetration - 2D Chl data '
+            IF( nqsr ==  3 )   WRITE(numout,*) '         R-G-B  light penetration - 3D Chl data '
+            IF( nqsr ==  4 )   WRITE(numout,*) '         2 band light penetration'
+            IF( nqsr ==  5 )   WRITE(numout,*) '         bio-model light penetration'
          ENDIF
          !
@@ -328,5 +381,5 @@
             !
             IF( nn_chldta == 1 ) THEN           !* Chl data : set sf_chl structure
-               IF(lwp) WRITE(numout,*)
+             IF(lwp) WRITE(numout,*)
                IF(lwp) WRITE(numout,*) '        Chlorophyll read in a file'
                ALLOCATE( sf_chl(1), STAT=ierror )
@@ -339,5 +392,5 @@
                CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',   &
                   &                                         'Solar penetration function of read chlorophyll', 'namtra_qsr' )
-               !
+               ! 
             ELSE                                !* constant Chl : compute once for all the distribution of light (etot3)
                IF(lwp) WRITE(numout,*)
