Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 4037)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 4038)
@@ -10,4 +10,5 @@
    !!             -   !  2005-11  (G. Madec) zco, zps, sco coordinate
    !!            3.2  !  2009-04  (G. Madec & NEMO team) 
+   !!            4.0  !  2012-05  (C. Rousset) store attenuation coef for use in ice model 
    !!----------------------------------------------------------------------
 
@@ -27,4 +28,5 @@
    USE iom             ! I/O manager
    USE fldread         ! read input fields
+   USE restart         ! ocean restart
    USE lib_mpp         ! MPP library
    USE wrk_nemo       ! Memory Allocation
@@ -43,4 +45,5 @@
    LOGICAL , PUBLIC ::   ln_qsr_2bd = .TRUE.    !: 2 band         light absorption flag
    LOGICAL , PUBLIC ::   ln_qsr_bio = .FALSE.   !: bio-model      light absorption flag
+   LOGICAL , PUBLIC ::   ln_qsr_ice = .FALSE.   !: light penetration for ice-model LIM3 (clem)
    INTEGER , PUBLIC ::   nn_chldta  = 0         !: use Chlorophyll data (=1) or not (=0)
    REAL(wp), PUBLIC ::   rn_abs     = 0.58_wp   !: fraction absorbed in the very near surface (RGB & 2 bands)
@@ -99,4 +102,5 @@
       REAL(wp) ::   zchl, zcoef, zfact   ! local scalars
       REAL(wp) ::   zc0, zc1, zc2, zc3   !    -         -
+      REAL(wp) ::   zzc0, zzc1, zzc2, zzc3   !    -         -
       REAL(wp) ::   zz0, zz1, z1_e3t     !    -         -
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zekb, zekg, zekr
@@ -158,4 +162,15 @@
          END DO
          CALL iom_put( 'qsr3d', etot3 )   ! Shortwave Radiation 3D distribution
+         ! clem: store attenuation coefficient of the first ocean level
+         IF ( ln_qsr_ice ) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF ( qsr(ji,jj) /= 0._wp ) THEN
+                     oatte(ji,jj) = ( qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) ) )
+                     iatte(ji,jj) = oatte(ji,jj)
+                  ENDIF
+               END DO
+            END DO
+         ENDIF
          !                                        ! ============================================== !
       ELSE                                        !  Ocean alone : 
@@ -216,4 +231,17 @@
                   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     )
+                        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) )
+                        oatte(ji,jj) = 1.0 - ( zzc0 + zzc1 + zzc2  + zzc3  ) * tmask(ji,jj,2) 
+                        iatte(ji,jj) = 1.0 - ( zzc0 + zzc1 + zcoef + zcoef ) * tmask(ji,jj,2)
+                     END DO
+                  END DO
+               ENDIF
                !
                DO jk = 1, nksr                                        ! compute and add qsr trend to ta
@@ -227,5 +255,10 @@
                   qsr_hc(:,:,jk) =  etot3(:,:,jk) * qsr(:,:)
                END DO
-            ENDIF
+               ! clem: store attenuation coefficient of the first ocean level
+               IF ( ln_qsr_ice ) THEN
+                  oatte(:,:) = etot3(:,:,1) / r1_rau0_rcp
+                  iatte(:,:) = oatte(:,:)
+               ENDIF
+           ENDIF
 
          ENDIF
@@ -246,4 +279,15 @@
                   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
+                        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 )
+                        oatte(ji,jj) = ( zc0*tmask(ji,jj,1) - zc1*tmask(ji,jj,2) ) / r1_rau0_rcp
+                        iatte(ji,jj) = oatte(ji,jj)
+                     END DO
+                  END DO
+               ENDIF
             ELSE                                               !* constant volume: coef. computed one for all
                DO jk = 1, nksr
@@ -254,4 +298,9 @@
                   END DO
                END DO
+               ! clem: store attenuation coefficient of the first ocean level
+               IF ( ln_qsr_ice ) THEN
+                  oatte(:,:) = etot3(:,:,1) / r1_rau0_rcp
+                  iatte(:,:) = oatte(:,:)
+               ENDIF
                !
             ENDIF
@@ -270,4 +319,15 @@
          !
       ENDIF
+      ! clem: store attenuation coefficient of the first ocean level
+      !IF (ln_traqsr) THEN
+      !   DO jj = 1, jpj
+      !      DO ji = 1, jpi
+      !         IF ( qsr(ji,jj) /= 0._wp ) THEN
+      !            oatte(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) )
+      !            iatte(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) )
+      !         ENDIF
+      !      END DO
+      !   END DO
+      !END IF
       !
       IF( lrst_oce ) THEN   !                  Write in the ocean restart file
@@ -325,5 +385,5 @@
       TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read
       !!
-      NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio,   &
+      NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice,  &
          &                  nn_chldta, rn_abs, rn_si0, rn_si1
       !!----------------------------------------------------------------------
@@ -331,4 +391,8 @@
       !
       IF( nn_timing == 1 )  CALL timing_start('tra_qsr_init')
+      !
+      ! clem init for oatte and iatte
+      oatte(:,:) = 1._wp
+      iatte(:,:) = 1._wp
       !
       CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        ) 
@@ -354,4 +418,5 @@
          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,*) '      light penetration for ice-model LIM3     ln_qsr_ice = ', ln_qsr_ice
          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
