Index: NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zopt.F90
===================================================================
--- NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zopt.F90	(revision 11710)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zopt.F90	(revision 13233)
@@ -80,10 +80,15 @@
       ze2(:,:,:) = 0._wp
       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(:,:,:) = trb(:,:,:,jpnch) + trb(:,:,:,jpdch)
       IF( ln_p5z )   zchl3d(:,:,:) = zchl3d(:,:,:)    + trb(:,:,:,jppch)
       !
+      ! Computation of the light attenuation parameters based on a 
+      ! look-up table
       DO jk = 1, jpkm1   
          DO jj = 1, jpj
@@ -99,12 +104,26 @@
          END DO
       END DO
-      !                                        !* 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, 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, nksrp      
             etot_ndcy(:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk)
@@ -117,27 +136,40 @@
             END DO
          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, zqsr_corr, ze1, ze2, ze3 ) 
-         !
+
+         ! Total PAR computation at level jk that includes the diurnal cycle
          DO jk = 1, nksrp      
             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, 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, nksrp      
-            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, nksrp      
-              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
@@ -146,4 +178,5 @@
 
 
+      ! Biophysical feedback part (computation of vertical penetration of SW)
       IF( ln_qsr_bio ) THEN                    !* heat flux accros w-level (used in the dynamics)
          !                                     !  ------------------------
@@ -154,8 +187,12 @@
             etot3(:,:,jk) =  ( ze0(:,:,jk) + ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk) ) * tmask(:,:,jk)
          END DO
-         !                                     !  ------------------------
-      ENDIF
-      !                                        !* Euphotic depth and level
-      neln   (:,:) = 1                            !  ------------------------
+
+      ENDIF
+
+      ! 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.
+      ! -------------------------------------------------------------------
       heup   (:,:) = gdepw_n(:,:,2)
       heup_01(:,:) = gdepw_n(:,:,2)
@@ -176,8 +213,11 @@
       END DO
       !
+      ! 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
@@ -187,6 +227,6 @@
             DO ji = 1, jpi
                IF( gdepw_n(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
-                  zetmp1 (ji,jj) = zetmp1 (ji,jj) + etot     (ji,jj,jk) * e3t_n(ji,jj,jk) ! remineralisation
-                  zetmp2 (ji,jj) = zetmp2 (ji,jj) + etot_ndcy(ji,jj,jk) * e3t_n(ji,jj,jk) ! production
+                  zetmp1 (ji,jj) = zetmp1 (ji,jj) + etot     (ji,jj,jk) * e3t_n(ji,jj,jk) ! Actual PAR
+                  zetmp2 (ji,jj) = zetmp2 (ji,jj) + etot_ndcy(ji,jj,jk) * e3t_n(ji,jj,jk) ! Par averaged over 24h
                   zdepmoy(ji,jj) = zdepmoy(ji,jj) +                       e3t_n(ji,jj,jk)
                ENDIF
@@ -195,5 +235,5 @@
       END DO
       !
-      emoy(:,:,:) = etot(:,:,:)       ! remineralisation
+      emoy(:,:,:) = etot(:,:,:)       ! PAR
       zpar(:,:,:) = etot_ndcy(:,:,:)  ! diagnostic : PAR with no diurnal cycle 
       !
@@ -209,5 +249,7 @@
          END DO
       END DO
-      !
+
+      ! Computation of the mean usable light for the different phytoplankton
+      ! groups based on their absorption characteristics.
       zdepmoy(:,:)   = 0.e0
       zetmp3 (:,:)   = 0.e0
@@ -218,6 +260,6 @@
             DO ji = 1, jpi
                IF( gdepw_n(ji,jj,jk+1) <= MIN(hmld(ji,jj), heup_01(ji,jj)) ) THEN
-                  zetmp3 (ji,jj) = zetmp3 (ji,jj) + enano    (ji,jj,jk) * e3t_n(ji,jj,jk) ! production
-                  zetmp4 (ji,jj) = zetmp4 (ji,jj) + ediat    (ji,jj,jk) * e3t_n(ji,jj,jk) ! production
+                  zetmp3 (ji,jj) = zetmp3 (ji,jj) + enano    (ji,jj,jk) * e3t_n(ji,jj,jk) ! Nanophytoplankton
+                  zetmp4 (ji,jj) = zetmp4 (ji,jj) + ediat    (ji,jj,jk) * e3t_n(ji,jj,jk) ! Diatoms 
                   zdepmoy(ji,jj) = zdepmoy(ji,jj) +                       e3t_n(ji,jj,jk)
                ENDIF
@@ -241,4 +283,5 @@
       !
       IF( ln_p5z ) THEN
+         ! Picophytoplankton when using PISCES-QUOTA
          zetmp5 (:,:) = 0.e0
          DO jk = 1, nksrp
@@ -246,5 +289,5 @@
                DO ji = 1, jpi
                   IF( gdepw_n(ji,jj,jk+1) <= MIN(hmld(ji,jj), heup_01(ji,jj)) ) THEN
-                     zetmp5(ji,jj)  = zetmp5 (ji,jj) + epico(ji,jj,jk) * e3t_n(ji,jj,jk) ! production
+                     zetmp5(ji,jj)  = zetmp5 (ji,jj) + epico(ji,jj,jk) * e3t_n(ji,jj,jk)
                   ENDIF
                END DO
@@ -265,4 +308,6 @@
          END DO
       ENDIF
+
+      ! Output of the diagnostics
       IF( lk_iomput ) THEN
         IF( knt == nrdttrc ) THEN
