Index: NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/P4Z/p4zopt.F90
===================================================================
--- NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/P4Z/p4zopt.F90	(revision 14384)
+++ NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/P4Z/p4zopt.F90	(revision 14385)
@@ -80,9 +80,14 @@
       ze3(:,:,:) = 0._wp
       !
-      !                                        !* attenuation coef. function of Chlorophyll and wavelength (Red-Green-Blue)
-      !                                        !  --------------------------------------------------------
+      ! Attenuation coef. function of Chlorophyll and wavelength (Red-Green-Blue)
+      ! Thus the light penetration scheme is based on a decomposition of PAR
+      ! into three wave length domains. This was first officially published
+      ! in Lengaigne et al. (2007).
+      ! --------------------------------------------------------
                      zchl3d(:,:,:) = tr(:,:,:,jpnch,Kbb) + tr(:,:,:,jpdch,Kbb)
-      IF( ln_p5z )   zchl3d(:,:,:) = zchl3d(:,:,:)    + tr(:,:,:,jppch,Kbb)
-      !
+      IF( ln_p5z )   zchl3d(:,:,:) = zchl3d(:,:,:)       + tr(:,:,:,jppch,Kbb)
+      !
+      ! Computation of the light attenuation parameters based on a 
+      ! look-up table
       DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          zchl = ( zchl3d(ji,jj,jk) + rtrn ) * 1.e6
@@ -94,12 +99,27 @@
          ekr(ji,jj,jk) = rkrgb(3,irgb) * e3t(ji,jj,jk,Kmm)
       END_3D
-      !                                        !* Photosynthetically Available Radiation (PAR)
-      !                                        !  --------------------------------------
+
+      ! Photosynthetically Available Radiation (PAR)
+      ! Two cases are considered in the following : 
+      ! (1) An explicit diunal cycle is activated. In that case, mean 
+      ! QSR is used as PISCES in its current state has not been parameterized
+      ! for an explicit diurnal cycle
+      ! (2) no diurnal cycle of SW is active and in that case, QSR is used.
+      ! --------------------------------------------
       IF( l_trcdm2dc ) THEN                     !  diurnal cycle
          !
+          !
+         ! SW over the ice free zone of the grid cell. This assumes that
+         ! SW is zero below sea ice which is a very crude assumption that is 
+         ! not fully correct with LIM3 and SI3 but no information is 
+         ! currently available to do a better job. SHould be improved in the 
+         ! (near) future.
          zqsr_corr(:,:) = qsr_mean(:,:) / ( 1.-fr_i(:,:) + rtrn )
          !
          CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 ) 
          !
+         ! Used PAR is computed for each phytoplankton species
+         ! etot_ndcy is PAR at level jk averaged over 24h.
+         ! Due to their size, they have different light absorption characteristics
          DO jk = 1, nksr      
             etot_ndcy(:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk)
@@ -113,26 +133,42 @@
          ENDIF
          !
+         ! SW over the ice free zone of the grid cell. This assumes that
+         ! SW is zero below sea ice which is a very crude assumption that is 
+         ! not fully correct with LIM3 and SI3 but no information is 
+         ! currently available to do a better job. SHould be improved in the 
+         ! (near) future.
+
          zqsr_corr(:,:) = qsr(:,:) / ( 1.-fr_i(:,:) + rtrn )
          !
          CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3 ) 
          !
+         ! Total PAR computation at level jk that includes the diurnal cycle
          DO jk = 1, nksr      
             etot(:,:,jk) =  ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk)
          END DO
          !
-      ELSE
-         !
+      ELSE   ! no diurnal cycle
+         !
+         !
+         ! SW over the ice free zone of the grid cell. This assumes that
+         ! SW is zero below sea ice which is a very crude assumption that is 
+         ! not fully correct with LIM3 and SI3 but no information is 
+         ! currently available to do a better job. SHould be improved in the 
+         ! (near) future.
          zqsr_corr(:,:) = qsr(:,:) / ( 1.-fr_i(:,:) + rtrn )
          !
          CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100  ) 
          !
+
+         ! Used PAR is computed for each phytoplankton species
+         ! Due to their size, they have different light absorption characteristics
          DO jk = 1, nksr      
-            etot (:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk)
-            enano(:,:,jk) =  1.85 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.46 * ze3(:,:,jk)
-            ediat(:,:,jk) =  1.62 * ze1(:,:,jk) + 0.74 * ze2(:,:,jk) + 0.63 * ze3(:,:,jk)
+            etot (:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk)    ! Total PAR
+            enano(:,:,jk) =  1.85 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.46 * ze3(:,:,jk)  ! Nanophytoplankton
+            ediat(:,:,jk) =  1.62 * ze1(:,:,jk) + 0.74 * ze2(:,:,jk) + 0.63 * ze3(:,:,jk)  ! Diatoms
          END DO
          IF( ln_p5z ) THEN
             DO jk = 1, nksr      
-              epico(:,:,jk) =  1.94 * ze1(:,:,jk) + 0.66 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)
+              epico(:,:,jk) =  1.94 * ze1(:,:,jk) + 0.66 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)  ! Picophytoplankton (PISCES-QUOTA)
             END DO
          ENDIF
@@ -141,4 +177,5 @@
 
 
+      ! Biophysical feedback part (computation of vertical penetration of SW)
       IF( ln_qsr_bio ) THEN                    !* heat flux accros w-level (used in the dynamics)
          !                                     !  ------------------------
@@ -151,6 +188,11 @@
          !                                     !  ------------------------
       ENDIF
-      !                                        !* Euphotic depth and level
-      neln   (:,:) = 1                            !  ------------------------
+      
+      ! Euphotic depth and level
+      ! Two definitions of the euphotic zone are used here. 
+      ! (1) The classical definition based on the relative threshold value
+      ! (2) An alternative definition based on a absolute threshold value.
+      ! -------------------------------------------------------------------
+      neln(:,:) = 1
       heup   (:,:) = gdepw(:,:,2,Kmm)
       heup_01(:,:) = gdepw(:,:,2,Kmm)
@@ -162,13 +204,16 @@
            heup(ji,jj) = gdepw(ji,jj,jk+1,Kmm)     ! Euphotic layer depth
         ENDIF
-        IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.50 )  THEN
+        IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.10 )  THEN
            heup_01(ji,jj) = gdepw(ji,jj,jk+1,Kmm)  ! Euphotic layer depth (light level definition)
         ENDIF
       END_3D
       !
+      ! The euphotic depth can not exceed 300 meters.
       heup   (:,:) = MIN( 300., heup   (:,:) )
       heup_01(:,:) = MIN( 300., heup_01(:,:) )
-      !                                        !* mean light over the mixed layer
-      zdepmoy(:,:)   = 0.e0                    !  -------------------------------
+      
+      ! Mean PAR over the mixed layer
+      ! -----------------------------
+      zdepmoy(:,:)   = 0.e0             
       zetmp1 (:,:)   = 0.e0
       zetmp2 (:,:)   = 0.e0
@@ -176,6 +221,6 @@
       DO_3D( 1, 1, 1, 1, 1, nksr )
          IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
-            zetmp1 (ji,jj) = zetmp1 (ji,jj) + etot     (ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! remineralisation
-            zetmp2 (ji,jj) = zetmp2 (ji,jj) + etot_ndcy(ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! production
+            zetmp1 (ji,jj) = zetmp1 (ji,jj) + etot     (ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! Actual PAR for remineralisation
+            zetmp2 (ji,jj) = zetmp2 (ji,jj) + etot_ndcy(ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! Par averaged over 24h for production
             zdepmoy(ji,jj) = zdepmoy(ji,jj) +                       e3t(ji,jj,jk,Kmm)
          ENDIF
@@ -192,5 +237,7 @@
          ENDIF
       END_3D
-      !
+
+      ! Computation of the mean usable light for the different phytoplankton
+      ! groups based on their absorption characteristics.
       zdepmoy(:,:)   = 0.e0
       zetmp3 (:,:)   = 0.e0
@@ -199,6 +246,6 @@
       DO_3D( 1, 1, 1, 1, 1, nksr )
          IF( gdepw(ji,jj,jk+1,Kmm) <= MIN(hmld(ji,jj), heup_01(ji,jj)) ) THEN
-            zetmp3 (ji,jj) = zetmp3 (ji,jj) + enano    (ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! production
-            zetmp4 (ji,jj) = zetmp4 (ji,jj) + ediat    (ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! production
+            zetmp3 (ji,jj) = zetmp3 (ji,jj) + enano    (ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! Nanophytoplankton
+            zetmp4 (ji,jj) = zetmp4 (ji,jj) + ediat    (ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! Diatoms
             zdepmoy(ji,jj) = zdepmoy(ji,jj) +                       e3t(ji,jj,jk,Kmm)
          ENDIF
@@ -216,8 +263,9 @@
       !
       IF( ln_p5z ) THEN
+         ! Picophytoplankton when using PISCES-QUOTA
          ALLOCATE( zetmp5(jpi,jpj) )  ;   zetmp5 (:,:) = 0.e0
          DO_3D( 1, 1, 1, 1, 1, nksr )
             IF( gdepw(ji,jj,jk+1,Kmm) <= MIN(hmld(ji,jj), heup_01(ji,jj)) ) THEN
-               zetmp5(ji,jj)  = zetmp5 (ji,jj) + epico(ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! production
+               zetmp5(ji,jj)  = zetmp5 (ji,jj) + epico(ji,jj,jk) * e3t(ji,jj,jk,Kmm)
             ENDIF
          END_3D
