Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 5206)
@@ -262,4 +262,7 @@
          <field id="qns_ice"      long_name="non-solar heat flux at ice surface"                           unit="W/m2"     />
          <field id="qtr_ice"      long_name="solar heat flux transmitted thru the ice"                     unit="W/m2"     />
+         <field id="qsr_oce_mean" long_name="daily mean solar heat flux at ocean surface"                  unit="W/m2"     />
+         <field id="qsr_ice_mean" long_name="daily mean solar heat flux at ice surface"                    unit="W/m2"     />
+         <field id="qtr_ice_mean" long_name="daily meansolar heat flux transmitted thru the ice"           unit="W/m2"     />
          <field id="utau_ice"     long_name="Wind stress along i-axis over the ice at i-point"             unit="N/m2"     />
          <field id="vtau_ice"     long_name="Wind stress along j-axis over the ice at i-point"             unit="N/m2"     />
@@ -571,41 +574,70 @@
       <!-- ptrc on T grid -->
 
-      <field_group id="ptrc_T" grid_ref="grid_T_3D">
-       <field id="DIC"      long_name="Dissolved inorganic Concentration"        unit="mmol/m3" />
-       <field id="Alkalini" long_name="Total Alkalinity Concentration"           unit="mmol/m3" />
-       <field id="O2"       long_name="Oxygen Concentration"                     unit="mmol/m3" />
-       <field id="CaCO3"    long_name="Calcite Concentration"                    unit="mmol/m3" />
-       <field id="PO4"      long_name="Phosphate Concentration"                  unit="mmol/m3" />
-       <field id="POC"      long_name="Small organic carbon Concentration"       unit="mmol/m3" />
-       <field id="Si"       long_name="Silicate Concentration"                   unit="mmol/m3" />
-       <field id="PHY"      long_name="(Nano)Phytoplankton Concentration"        unit="mmol/m3" />
-       <field id="ZOO"      long_name="(Micro)Zooplankton Concentration"         unit="mmol/m3" />
-       <field id="DOC"      long_name="Dissolved organic Concentration"          unit="mmol/m3" />
-       <field id="PHY2"     long_name="Diatoms Concentration"                    unit="mmol/m3" />
-       <field id="ZOO2"     long_name="Mesozooplankton Concentration"            unit="mmol/m3" />
-       <field id="DSi"      long_name="Diatoms Silicate Concentration"           unit="mmol/m3" />
-       <field id="Fer"      long_name="Dissolved Iron Concentration"             unit="mmol/m3" />
-       <field id="BFe"      long_name="Big iron particles Concentration"         unit="mmol/m3" />
-       <field id="GOC"      long_name="Big organic carbon Concentration"         unit="mmol/m3" />
-       <field id="SFe"      long_name="Small iron particles Concentration"       unit="mmol/m3" />
-       <field id="DFe"      long_name="Diatoms iron  Concentration"              unit="mmol/m3" />
-       <field id="GSi"      long_name="Sinking biogenic Silicate Concentration"  unit="mmol/m3" />
-       <field id="NFe"      long_name="Nano iron Concentration"                  unit="mmol/m3" />
-       <field id="NCHL"     long_name="Nano chlorophyl Concentration"            unit="mg/m3"   />
-       <field id="DCHL"     long_name="Diatoms chlorophyl Concentration"         unit="mg/m3"   />
-       <field id="NO3"      long_name="Nitrate Concentration"                    unit="mmol/m3" />
-       <field id="NH4"      long_name="Ammonium Concentration"                   unit="mmol/m3" />
+     <field_group id="ptrc_T" grid_ref="grid_T_3D">
+       <field id="DIC"          long_name="Dissolved inorganic Concentration"        unit="mmol/m3" />
+       <field id="DIC_E3T"      long_name="DIC * E3T"                                unit="mmol/m2" > DIC * e3t </field >
+       <field id="Alkalini"     long_name="Total Alkalinity Concentration"           unit="mmol/m3" />
+       <field id="Alkalini_E3T" long_name="Alkalini * E3T"                           unit="mmol/m2"  > Alkalini * e3t </field >
+       <field id="O2"           long_name="Oxygen Concentration"                     unit="mmol/m3" />
+       <field id="O2_E3T"       long_name="O2 * E3T"                                 unit="mmol/m2"  > O2 * e3t </field >
+       <field id="CaCO3"        long_name="Calcite Concentration"                    unit="mmol/m3" />
+       <field id="CaCO3_E3T"    long_name="CaCO3 * E3T"                              unit="mmol/m2"  > CaCO3 * e3t </field >
+       <field id="PO4"          long_name="Phosphate Concentration"                  unit="mmol/m3" />
+       <field id="PO4_E3T"      long_name="PO4 * E3T"                                unit="mmol/m2"  > PO4 * e3t </field >
+       <field id="POC"          long_name="Small organic carbon Concentration"       unit="mmol/m3" />
+       <field id="POC_E3T"      long_name="POC * E3T"                                unit="mmol/m2"  > POC * e3t </field >
+       <field id="Si"           long_name="Silicate Concentration"                   unit="mmol/m3" />
+       <field id="Si_E3T"       long_name="Si * E3T"                                 unit="mmol/m2"  > Si * e3t </field >
+       <field id="PHY"          long_name="(Nano)Phytoplankton Concentration"        unit="mmol/m3" />
+       <field id="PHY_E3T"      long_name="PHY * E3T"                                unit="mmol/m2"  > PHY * e3t </field >
+       <field id="ZOO"          long_name="(Micro)Zooplankton Concentration"         unit="mmol/m3" />
+       <field id="ZOO_E3T"      long_name="ZOO2 * E3T"                               unit="mmol/m2"  > ZOO * e3t </field >
+       <field id="DOC"          long_name="Dissolved organic Concentration"          unit="mmol/m3" />
+       <field id="DOC_E3T"      long_name="DOC * E3T"                                unit="mmol/m2"  > DOC * e3t </field >
+       <field id="PHY2"         long_name="Diatoms Concentration"                    unit="mmol/m3" />
+       <field id="PHY2_E3T"     long_name="PHY2 * E3T"                               unit="mmol/m2"  > PHY2 * e3t </field >
+       <field id="ZOO2"         long_name="Mesozooplankton Concentration"            unit="mmol/m3" />
+       <field id="ZOO2_E3T"     long_name="ZOO2 * E3T"                               unit="mmol/m2"  > ZOO2 * e3t </field >
+       <field id="DSi"          long_name="Diatoms Silicate Concentration"           unit="mmol/m3" />
+       <field id="DSi_E3T"      long_name="Dsi * E3T"                                unit="mmol/m2"  > DSi * e3t </field >
+       <field id="Fer"          long_name="Dissolved Iron Concentration"             unit="mmol/m3" />
+       <field id="Fer_E3T"      long_name="Fer * E3T"                                unit="mmol/m2"  > Fer * e3t </field >
+       <field id="BFe"          long_name="Big iron particles Concentration"         unit="mmol/m3" />
+       <field id="BFe_E3T"      long_name="BFe * E3T"                                unit="mmol/m2"  > BFe * e3t </field >
+       <field id="GOC"          long_name="Big organic carbon Concentration"         unit="mmol/m3" />
+       <field id="GOC_E3T"      long_name="GOC * E3T"                                unit="mmol/m2"  > GOC * e3t </field >
+       <field id="SFe"          long_name="Small iron particles Concentration"       unit="mmol/m3" />
+       <field id="SFe_E3T"      long_name="SFe * E3T"                                unit="mmol/m2"  > SFe * e3t </field >
+       <field id="DFe"          long_name="Diatoms iron  Concentration"              unit="mmol/m3" />
+       <field id="DFe_E3T"      long_name="DFe * E3T"                                unit="mmol/m2"  > DFe * e3t </field >
+       <field id="GSi"          long_name="Sinking biogenic Silicate Concentration"  unit="mmol/m3" />
+       <field id="GSi_E3T"      long_name="GSi * E3T"                                unit="mmol/m2"  > GSi * e3t </field >
+       <field id="NFe"          long_name="Nano iron Concentration"                  unit="mmol/m3" />
+       <field id="NFe_E3T"      long_name="NFe * E3T"                                unit="mmol/m2"  > NFe * e3t </field >
+       <field id="NCHL"         long_name="Nano chlorophyl Concentration"            unit="mg/m3"   />
+       <field id="NCHL_E3T"     long_name="NCHL * E3T"                               unit="mmol/m2"  > NCHL * e3t </field >
+       <field id="DCHL"         long_name="Diatoms chlorophyl Concentration"         unit="mg/m3"   />
+       <field id="DCHL_E3T"     long_name="DCHL * E3T"                               unit="mmol/m2"  > DCHL * e3t </field >
+       <field id="NO3"          long_name="Nitrate Concentration"                    unit="mmol/m3" />
+       <field id="NO3_E3T"      long_name="NO3 * E3T"                                unit="mmol/m2"  > NO3 * e3t </field >
+       <field id="NH4"          long_name="Ammonium Concentration"                   unit="mmol/m3" />
+       <field id="NH4_E3T"      long_name="NH4 * E3T"                                unit="mmol/m2"  > NH4 * e3t </field >
 
        <!-- PISCES with Kriest parametisation : variables available with key_kriest -->
-       <field id="Num"      long_name="Number of organic particles"              unit="nbr" />
+       <field id="Num"         long_name="Number of organic particles"              unit="nbr" />
+       <field id="Num_E3T"     long_name="Num * E3T"                                unit="nbr*m"  > Num * e3t </field >
 
        <!-- PISCES light : variables available with key_pisces_reduced -->
-       <field id="DET"      long_name="Detritus"                                 unit="mmol-N/m3" />
-       <field id="DOM"      long_name="Dissolved Organic Matter"                 unit="mmol-N/m3" />
+       <field id="DET"         long_name="Detritus"                                 unit="mmol-N/m3" />
+       <field id="DET_E3T"     long_name="DET * E3T"                                unit="mmol-N/m2"  > DET * e3t </field >
+       <field id="DOM"         long_name="Dissolved Organic Matter"                 unit="mmol-N/m3" />
+       <field id="DOM_E3T"     long_name="DOM * E3T"                                unit="mmol-N/m2"  > DOM * e3t </field >
 
        <!-- CFC11 : variables available with key_cfc -->
-       <field id="CFC11"    long_name="CFC-11 Concentration"                     unit="umol/L" />
+       <field id="CFC11"       long_name="CFC-11 Concentration"                     unit="umol/L" />
+       <field id="CFC11_E3T"   long_name="CFC11 * E3T"                              unit="umol/m2"  > CFC11 * e3t </field >
        <!-- Bomb C14 : variables available with key_c14b -->
-       <field id="C14B"     long_name="Bomb C14 Concentration"                   unit="ration" />
+       <field id="C14B"        long_name="Bomb C14 Concentration"                   unit="ration" />
+       <field id="C14B_E3T"    long_name="C14B * E3T"                               unit="ration*m"  > C14B * e3t </field >
      </field_group>
 
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90	(revision 5206)
@@ -65,7 +65,4 @@
 
    !!* Ice Rheology
-
-   LOGICAL , PUBLIC::  ltrcdm2dc_ice = .FALSE.              !: In case of Diurnal Cycle short wave, compute a Daily Mean short waves flux
-
 # if defined key_lim2_vp
    !                                                      !!* VP rheology *
@@ -113,5 +110,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qdtcn         !: energy from the ice to the ocean point (at a factor 2)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   thcm          !: part of the solar energy used in the lead heat budget
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fstric_daymean!: Solar flux transmitted trough the ice, for day mean of qsr ( PISCES, with dirunal cycle )
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fstric_mean   !: Solar flux transmitted trough the ice, for day mean of qsr ( PISCES, with dirunal cycle )
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fstric        !: Solar flux transmitted trough the ice
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ffltbif       !: linked with the max heat contained in brine pockets (?)
@@ -173,6 +170,4 @@
          &     tbif  (jpi,jpj,jplayersp1)                           , STAT=ierr(5))
 
-      IF( ltrcdm2dc_ice ) ALLOCATE(fstric_daymean(jpi,jpj), STAT=ierr(6) )
-
       !* moment used in the advection scheme
       ALLOCATE(sxice (jpi,jpj) , syice (jpi,jpj) , sxxice(jpi,jpj) ,     &
@@ -201,5 +196,4 @@
    !!   Default option         Empty module        NO LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-   LOGICAL , PUBLIC::  ltrcdm2dc_ice = .FALSE.              !: In case of Diurnal Cycle short wave, compute a Daily Mean short waves flux
 #endif
    !!-----------------------------------------------------------------
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_2/iceini_2.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_2/iceini_2.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_2/iceini_2.F90	(revision 5206)
@@ -60,9 +60,4 @@
       ENDIF
       !                                
-      ! When Diurnal cycle, core bulk and LIM2  are activated, 
-      ! a daily mean qsr is computed for tracer/biogeochemistery model                                 !
-      IF( ltrcdm2dc )THEN  ; ltrcdm2dc_ice = .TRUE.
-      ELSE                 ; ltrcdm2dc_ice = .FALSE.
-      ENDIF
       !                                ! Allocate the ice arrays
       ierr =        ice_alloc_2    ()       ! ice variables
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 5206)
@@ -46,5 +46,4 @@
    PUBLIC   lim_sbc_flx_2      ! called by sbc_ice_lim_2
    PUBLIC   lim_sbc_tau_2      ! called by sbc_ice_lim_2
-   PUBLIC   lim_bio_meanqsr_2  ! called by sbc_ice_lim_2
 
    REAL(wp)  ::   r1_rdtice            ! = 1. / rdt_ice 
@@ -264,4 +263,8 @@
       ENDIF
 
+      !   daily mean qsr when diurnal cycle is applied on physics - for BGC models
+      IF( l_trcdm2dc ) qsr_mean(:,:) =  pfrld(:,:) * qsr_mean(:,:) + ( 1. - pfrld(:,:) ) * fstric_mean(:,:)
+
+
       IF(ln_ctl) THEN            ! control print
          CALL prt_ctl(tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns   , clinfo2=' qns     : ')
@@ -430,18 +433,4 @@
    END SUBROUTINE lim_sbc_tau_2
 
-   SUBROUTINE lim_bio_meanqsr_2
-      !!---------------------------------------------------------------------
-      !!                     ***  ROUTINE lim_bio_meanqsr
-      !!
-      !! ** Purpose :   provide daily qsr_mean for PISCES when 
-      !!                analytic diurnal cycle is applied in physic
-      !!
-      !! ** Method  :   add part under ice
-      !!
-      !!---------------------------------------------------------------------
-
-      qsr_mean(:,:) =  pfrld(:,:) * qsr_mean(:,:) + ( 1. - pfrld(:,:) ) * fstric_daymean(:,:)
-
-   END SUBROUTINE lim_bio_meanqsr_2
 
    SUBROUTINE lim_sbc_init_2
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90	(revision 5206)
@@ -114,5 +114,8 @@
       CALL wrk_alloc( jpi, jpj, jpk, zmsk )
 
-      IF( kt == nit000 )   CALL lim_thd_init_2  ! Initialization (first time-step only)
+      IF( kt == nit000 ) THEN
+         CALL lim_thd_init_2  ! Initialization (first time-step only)
+         IF( l_trcdm2dc ) ALLOCATE( fstric_mean(jpi,jpj), fstbif_mean_1d(jpij), qsr_ice_mean_1d(jpij) )
+      ENDIF
    
       !-------------------------------------------!
@@ -137,5 +140,6 @@
       rdq_ice(:,:) = 0.e0   ! heat content associated with rdm_ice
       zmsk (:,:,:) = 0.e0
-      IF( ltrcdm2dc_ice  ) fstric_daymean (:,:) = 0.e0   ! part of solar radiation absorbing inside the ice
+      !
+      IF( l_trcdm2dc  ) fstric_mean(:,:) = 0.e0   ! part of solar radiation absorbing inside the ice
 
       ! set to zero snow thickness smaller than epsi04
@@ -285,5 +289,5 @@
          CALL tab_2d_1d_2( nbpb, tbif_1d    (1:nbpb , 3 ), tbif(:,:,3)    , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d_2( nbpb, qsr_ice_1d (1:nbpb)     , qsr_ice(:,:,1) , jpi, jpj, npb(1:nbpb) )
-         IF( ltrcdm2dc_ice ) &
+         IF( l_trcdm2dc ) &
          & CALL tab_2d_1d_2( nbpb, qsr_ice_mean_1d (1:nbpb), qsr_ice_mean(:,:,1), jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d_2( nbpb, fr1_i0_1d  (1:nbpb)     , fr1_i0         , jpi, jpj, npb(1:nbpb) )
@@ -336,9 +340,9 @@
          CALL tab_1d_2d_2( nbpb, qsr_ice(:,:,1), npb, qsr_ice_1d(1:nbpb)  , jpi, jpj )
          CALL tab_1d_2d_2( nbpb, qns_ice(:,:,1), npb, qns_ice_1d(1:nbpb)  , jpi, jpj )
-         IF( ltrcdm2dc_ice )THEN
-            CALL tab_1d_2d_2( nbpb, fstric_daymean     , npb, fstbif_daymean_1d (1:nbpb)  , jpi, jpj )
-            CALL tab_1d_2d_2( nbpb, qsr_ice_mean(:,:,1), npb,    qsr_ice_mean_1d(1:nbpb)  , jpi, jpj )
+         IF( l_trcdm2dc ) &
+            CALL tab_1d_2d_2( nbpb, fstric_mean        , npb, fstbif_mean_1d (1:nbpb), jpi, jpj )
+            CALL tab_1d_2d_2( nbpb, qsr_ice_mean(:,:,1), npb, qsr_ice_mean_1d(1:nbpb), jpi, jpj )
          ENDIF
-         IF( .NOT. lk_cpl )   CALL tab_1d_2d_2( nbpb, qla_ice(:,:,1), npb, qla_ice_1d(1:nbpb)             , jpi, jpj )
+         IF( .NOT. lk_cpl )   CALL tab_1d_2d_2( nbpb, qla_ice(:,:,1), npb, qla_ice_1d(1:nbpb), jpi, jpj )
          !
       ENDIF
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_2/limthd_zdf_2.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_2/limthd_zdf_2.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_2/limthd_zdf_2.F90	(revision 5206)
@@ -273,14 +273,14 @@
        END DO
 
-       IF( ltrcdm2dc_ice )THEN
-
+       IF( l_topdcy )THEN
+          !
           DO ji = kideb , kiut
-             zihsn  = MAX( zzero , SIGN (zone , -h_snow_1d(ji) ) )
-             zihic  = MAX( zzero , 1.0 - ( h_ice_1d(ji) / zhsu ) )
+             zihsn    = MAX( zzero , SIGN (zone , -h_snow_1d(ji) ) )
+             zihic    = MAX( zzero , 1.0 - ( h_ice_1d(ji) / zhsu ) )
              zi0(ji)  = zihsn * ( fr1_i0_1d(ji) + zihic * fr2_i0_1d(ji) )
              zexp     = MIN( zone , EXP( -1.5 * ( h_ice_1d(ji) - zhsu ) ) )
-             fstbif_daymean_1d(ji) = zi0(ji) * qsr_ice_mean_1d(ji) * zexp
+             fstbif_mean_1d(ji) = zi0(ji) * qsr_ice_mean_1d(ji) * zexp
           END DO
-
+          !
        ENDIF
 
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_2/thd_ice_2.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_2/thd_ice_2.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_2/thd_ice_2.F90	(revision 5206)
@@ -55,5 +55,5 @@
       fstbif_1d   ,     &  !:    "                  "      fstric
       fltbif_1d   ,     &  !:    "                  "      ffltbif
-      fstbif_daymean_1d, & !:    "                  " fstric_daymean
+      fstbif_mean_1d,   & !:    "                   "      fstric_mean
       fscbq_1d    ,     &  !:    "                  "      fscmcbq
       qsr_ice_1d  ,     &  !:    "                  "      qsr_ice
@@ -122,6 +122,4 @@
          &      tbif_1d(jpij, jplayersp1), Stat=ierr(4))
          !
-      IF( ltrcdm2dc_ice )ALLOCATE(fstbif_daymean_1d(jpij),qsr_ice_mean_1d(jpij),Stat=ierr(5))
-         !
       thd_ice_alloc_2 = MAXVAL(ierr)
       IF( thd_ice_alloc_2 /= 0 )   CALL ctl_warn('thd_ice_alloc_2: failed to allocate arrays')
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 5206)
@@ -301,4 +301,6 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ftr_ice   !: transmitted solar radiation under ice
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ftr_ice_mean  !: transmitted daily mean solar radiation under ice (diurnal cycle) 
+
 
    !!--------------------------------------------------------------------------
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 5206)
@@ -118,4 +118,9 @@
       IF( iom_use('qt_oce' ) )   CALL iom_put( "qt_oce"  , ( qsr(:,:) + qns(:,:) ) * pfrld(:,:) )  
       IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) ) * a_i_b(:,:,:), dim=3 ) )
+      IF( l_trcdm2dc ) THEN
+        IF( iom_use('qsr_oce_mean') )   CALL iom_put( "qsr_oce_mean" , qsr_mean(:,:) * pfrld(:,:) )   !   daily mean solar flux at ocean surface
+        IF( iom_use('qsr_ice_mean') )   CALL iom_put( "qsr_ice_mean" , SUM( qsr_ice_mean(:,:,:) * a_i_b(:,:,:), dim=3 ) ) !  daily mean solar flux at ice surface
+        IF( iom_use('qtr_ice_mean') )   CALL iom_put( "qtr_ice_mean" , SUM( ftr_ice_mean(:,:,:) * a_i_b(:,:,:), dim=3 ) ) !  daily mean solar flux transmitted thru ice
+      ENDIF
 
       ! pfrld is the lead fraction at the previous time step (actually between TRP and THD)
@@ -212,4 +217,8 @@
       fr_i  (:,:)   = at_i(:,:)             ! Sea-ice fraction            
       tn_ice(:,:,:) = t_su(:,:,:)           ! Ice surface temperature                      
+
+
+      !   daily mean qsr when diurnal cycle is applied on physics - for BGC models
+      IF( l_trcdm2dc ) qsr_mean(:,:) =  pfrld(:,:) * qsr_mean(:,:) +  ftr_ice_mean(:,:)
 
       !------------------------------------------------!
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 5206)
@@ -98,4 +98,6 @@
       IF( nn_timing == 1 )  CALL timing_start('limthd')
 
+      IF( kt == nit000 .AND. l_trcdm2dc )  ALLOCATE( ftr_ice_mean(jpi,jpj,jpl), ftr_ice_mean_1d(jpij), qsr_ice_mean_1d(jpij) )
+
       ! conservation test
       IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
@@ -106,4 +108,6 @@
       !------------------------------------------------------------------------!
       ftr_ice(:,:,:) = 0._wp  ! part of solar radiation transmitted through the ice
+      IF( l_trcdm2dc ) ftr_ice_mean(:,:,:) = 0._wp   ! part of daily mean solar radiation absorbing inside the ice
+
 
       !--------------------
@@ -576,4 +580,8 @@
          CALL tab_2d_1d( nbpb, qns_ice_1d (1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, ftr_ice_1d (1:nbpb), ftr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
+         IF( l_trcdm2dc )  THEN
+            CALL tab_2d_1d( nbpb, qsr_ice_mean_1d (1:nbpb), qsr_ice_mean(:,:,jl) , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nbpb, ftr_ice_mean_1d (1:nbpb), ftr_ice_mean(:,:,jl) , jpi, jpj, npb(1:nbpb) )
+         ENDIF
          IF( .NOT. lk_cpl ) THEN
             CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb), qla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
@@ -670,5 +678,10 @@
          CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj)
          CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj )
-                  
+         !         
+         IF( l_trcdm2dc )  THEN
+            CALL tab_1d_2d( nbpb, qns_ice_mean(:,:,jl), npb, qns_ice_mean_1d(1:nbpb) , jpi, jpj)
+            CALL tab_1d_2d( nbpb, ftr_ice_mean(:,:,jl), npb, ftr_ice_mean_1d(1:nbpb) , jpi, jpj )
+         ENDIF
+         !
       END SELECT
 
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 5206)
@@ -126,4 +126,5 @@
       REAL(wp), POINTER, DIMENSION(:)     ::   zdifcase    ! case of the equation resolution (1->4)
       REAL(wp), POINTER, DIMENSION(:)     ::   zftrice     ! solar radiation transmitted through the ice
+      REAL(wp), POINTER, DIMENSION(:)     ::   zftrice_mean ! daily mean solar radiation transmitted through the ice
       REAL(wp), POINTER, DIMENSION(:)     ::   zihic
       
@@ -174,4 +175,5 @@
       CALL wrk_alloc( jpij,nlay_i+3, zindterm, zindtbis, zdiagbis  )
       CALL wrk_alloc( jpij,nlay_i+3,3, ztrid )
+      IF( l_trcdm2dc )  CALL wrk_alloc( jpij, zftrice_mean )
 
       CALL wrk_alloc( jpij, zdq, zq_ini, zhfx_err )
@@ -249,4 +251,10 @@
       END DO
 
+      IF( l_trcdm2dc ) THEN
+         DO ji = kideb , kiut
+            zftrice_mean(ji) =  qsr_ice_mean_1d(ji) * i0(ji) ! Solar radiation transmitted below the surface layer
+         END DO
+      ENDIF
+
       !---------------------------------------------------------
       ! Transmission - absorption of solar radiation in the ice
@@ -282,4 +290,14 @@
          ftr_ice_1d(ji) = zradtr_i(ji,nlay_i) 
       END DO
+
+
+      IF( l_trcdm2dc ) THEN
+         DO ji = kideb , kiut
+             ftr_ice_mean_1d(ji) =  ftr_ice_mean_1d(ji) & 
+                 &                + a_i_1d(ji) * zftrice_mean(ji) &
+                 &                             * EXP( - rn_kappa_i * ( MAX ( 0._wp , h_i_1d(ji) ) ) ) &
+                                               * EXP( - zraext_s   * ( MAX ( 0._wp , h_s_1d(ji) ) ) )
+         END DO
+      ENDIF
 
       !
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 5206)
@@ -41,4 +41,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ftr_ice_1d   
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qsr_ice_1d  
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qsr_ice_mean_1d  
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ftr_ice_mean_1d   
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fr1_i0_1d   
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fr2_i0_1d   
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90	(revision 5206)
@@ -152,6 +152,6 @@
 #endif
          !
-#if defined key_lim2
-      IF( ltrcdm2dc_ice )   ALLOCATE( qsr_ice_mean (jpi,jpj,jpl), STAT=ierr(3) )
+#if defined key_lim3 || defined key_lim2 
+      IF( l_trcdm2dc )   ALLOCATE( qsr_ice_mean (jpi,jpj,jpl), STAT=ierr(3) )
 #endif
          !
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 5206)
@@ -81,6 +81,7 @@
    !!              Ocean Surface Boundary Condition fields
    !!----------------------------------------------------------------------
+   LOGICAL , PUBLIC ::   l_trcdm2dc               !: In case of Diurnal Cycle short wave, compute a Daily Mean short waves flux
+   !
    LOGICAL , PUBLIC ::   lhftau = .FALSE.        !: HF tau used in TKE: mean(stress module) - module(mean stress)
-   LOGICAL , PUBLIC ::   ltrcdm2dc               !: In case of Diurnal Cycle short wave, compute a Daily Mean short waves flux
    !!                                   !!   now    ! before   !!
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau   , utau_b   !: sea surface i-stress (ocean referential)     [N/m2]
@@ -136,5 +137,5 @@
       !!                  ***  FUNCTION sbc_oce_alloc  ***
       !!---------------------------------------------------------------------
-      INTEGER :: ierr(5)
+      INTEGER :: ierr(6)
       !!---------------------------------------------------------------------
       ierr(:) = 0
@@ -158,9 +159,9 @@
          &      ssv_m  (jpi,jpj) , sss_m  (jpi,jpj), ssh_m(jpi,jpj) , STAT=ierr(4) )
          !
+         IF( l_trcdm2dc )   ALLOCATE( qsr_mean(jpi,jpj) , STAT=ierr(5) )
+         !
 #if defined key_vvl
-      ALLOCATE( e3t_m(jpi,jpj) , STAT=ierr(5) )
-#endif
-         !
-      IF( ltrcdm2dc ) ALLOCATE( qsr_mean(jpi,jpj) , STAT=ierr(5) )
+      ALLOCATE( e3t_m(jpi,jpj) , STAT=ierr(6) )
+#endif
          !
       sbc_oce_alloc = MAXVAL( ierr )
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 5206)
@@ -196,6 +196,5 @@
       IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce_core( kt, sf, sst_m, ssu_m, ssv_m )
 
-      ! If diurnal cycle is activated, compute a daily mean short waves flux for biogeochemistery
-      IF( ltrcdm2dc )   CALL blk_bio_meanqsr
+      IF( l_trcdm2dc )   CALL blk_bio_meanqsr  ! diurnal cycle : daily mean short waves flux for biogeochemistery
 
 #if defined key_cice
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 5206)
@@ -172,4 +172,6 @@
             IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   &
                &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx )
+
+            IF( l_trcdm2dc )       CALL blk_ice_meanqsr( zalb_ice, qsr_ice_mean, jpl )
             !
          CASE ( jp_cpl )
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90	(revision 5206)
@@ -194,4 +194,5 @@
                &                      tprecip    , sprecip    ,                         &
                &                      fr1_i0     , fr2_i0     , cp_ice_msh , jpl  )
+
             IF( ltrcdm2dc_ice )   CALL blk_ice_meanqsr( zalb_ice, qsr_ice_mean, jpl )
 
@@ -237,7 +238,4 @@
                            CALL lim_thd_2      ( kt )      ! Ice thermodynamics 
                            CALL lim_sbc_flx_2  ( kt )      ! update surface ocean mass, heat & salt fluxes 
-#if defined key_top
-        IF( ltrcdm2dc_ice )CALL lim_bio_meanqsr_2
-#endif
 
          IF(  .NOT. lk_mpp )THEN
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 5206)
@@ -24,4 +24,5 @@
    USE phycst           ! physical constants
    USE sbc_oce          ! Surface boundary condition: ocean fields
+   USE trc_oce          ! shared ocean-passive tracers variables
    USE sbc_ice          ! Surface boundary condition: ice fields
    USE sbcdcy           ! surface boundary condition: diurnal cycle
@@ -151,17 +152,14 @@
       END SELECT
       !
-#if defined key_top && ! defined key_offline
-      ltrcdm2dc = (ln_dm2dc .AND. ln_blk_core .AND. nn_ice==2)
-      IF( ltrcdm2dc )THEN
-         IF(lwp)THEN
-            WRITE(numout,*)"analytical diurnal cycle, core bulk formulation and LIM2 use: "
-            WRITE(numout,*)"Diurnal cycle on physics but not in passive tracers"
+      IF( ln_dm2dc .AND. lk_top  )  THEN   
+         l_trcdm2dc = .TRUE.
+         IF( lwp ) THEN
+            WRITE(numout,*) "In case of diurnal cycle coupled with passive tracers"
+            WRITE(numout,*) "Computation of a daily mean shortwave for some biogeochemical models"
          ENDIF
-      ENDIF
-#else 
-      ltrcdm2dc =  .FALSE.
-#endif
-
-      !
+      ELSE             
+         l_trcdm2dc = .FALSE.
+      ENDIF
+
       !                              ! allocate sbc arrays
       IF( sbc_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_oce arrays' )
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/trc_oce.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/trc_oce.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/trc_oce.F90	(revision 5206)
@@ -32,4 +32,5 @@
    !!   'key_top'                                                 bio-model          
    !!----------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER ::   lk_top     = .TRUE.   !: TOP model
    LOGICAL, PUBLIC, PARAMETER ::   lk_qsr_bio = .TRUE.   !: bio-model light absorption flag
 #else
@@ -37,4 +38,5 @@
    !! Default option                          No bio-model light absorption      
    !!----------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER ::   lk_top     = .FALSE.   !: TOP model
    LOGICAL, PUBLIC, PARAMETER ::   lk_qsr_bio = .FALSE.   !: bio-model light absorption flag
 #endif
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zopt.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zopt.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zopt.F90	(revision 5206)
@@ -89,5 +89,8 @@
 
       !                                          ! surface irradiance
-      zpar0m (:,:)   = qsr   (:,:) * 0.43        ! ------------------
+      !                                          ! ------------------
+      IF( ln_dm2dc ) THEN   ;   zpar0m(:,:) = qsr_mean(:,:) * 0.43
+      ELSE                  ;   zpar0m(:,:) = qsr     (:,:) * 0.43
+      ENDIF
       zpar100(:,:)   = zpar0m(:,:) * 0.01
       zparr  (:,:,1) = zpar0m(:,:) * 0.5
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 5206)
@@ -172,6 +172,6 @@
                ztem1  = MAX( 0., tsn(ji,jj,jk,jp_tem) )
                ztem2  = tsn(ji,jj,jk,jp_tem) - 10.
-               zetot1 = MAX( 0., etot(ji,jj,jk) - 1.) / ( 4. + etot(ji,jj,jk) ) 
-               zetot2 = 30. / ( 30. + etot(ji,jj,jk) ) 
+               zetot1 = MAX( 0., etot_ndcy(ji,jj,jk) - 1.) / ( 4. + etot_ndcy(ji,jj,jk) ) 
+               zetot2 = 30. / ( 30. + etot_ndcy(ji,jj,jk) ) 
 
                xfracal(ji,jj,jk) = caco3r * MIN( zlim1, zlim2, zlim3 )                  &
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 5206)
@@ -35,4 +35,5 @@
    REAL(wp) :: parlux      !: Fraction of shortwave as PAR
    REAL(wp) :: xparsw                 !: parlux/3
+   REAL(wp) :: xsi0r                 !:  1. /rn_si0
 
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_par      ! structure of input par
@@ -42,4 +43,5 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: enano, ediat   !: PAR for phyto, nano and diat 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: etot_ndcy      !: PAR over 24h in case of diurnal cycle
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: emoy           !: averaged PAR in the mixed layer
 
@@ -71,8 +73,8 @@
       INTEGER  ::   ji, jj, jk
       INTEGER  ::   irgb
-      REAL(wp) ::   zchl, zxsi0r
+      REAL(wp) ::   zchl
       REAL(wp) ::   zc0 , zc1 , zc2, zc3, z1_dep
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zdepmoy, zetmp, zetmp1, zetmp2
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zekg, zekr, zekb, ze0, ze1, ze2, ze3
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4, zqsr100
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpar, zekg, zekr, zekb, ze0, ze1, ze2, ze3
       !!---------------------------------------------------------------------
       !
@@ -80,15 +82,14 @@
       !
       ! Allocate temporary workspace
-      CALL wrk_alloc( jpi, jpj,      zdepmoy, zetmp, zetmp1, zetmp2 ) 
-      CALL wrk_alloc( jpi, jpj, jpk, zekg, zekr, zekb, ze0, ze1, ze2, ze3 )
-
-      IF( jnt == 1 .AND. ln_varpar ) CALL p4z_optsbc( kt )
+      CALL wrk_alloc( jpi, jpj,      zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 )
+      CALL wrk_alloc( jpi, jpj, jpk, zpar, zekg, zekr, zekb, ze0, ze1, ze2, ze3 )
+
+      IF( jnt == 1 .AND. ln_varpar ) CALL p4z_opt_sbc( kt )
 
       !     Initialisation of variables used to compute PAR
       !     -----------------------------------------------
-      ze1(:,:,jpk) = 0._wp
-      ze2(:,:,jpk) = 0._wp
-      ze3(:,:,jpk) = 0._wp
-
+      ze1(:,:,:) = 0._wp
+      ze2(:,:,:) = 0._wp
+      ze3(:,:,:) = 0._wp
       !                                        !* attenuation coef. function of Chlorophyll and wavelength (Red-Green-Blue)
       DO jk = 1, jpkm1                         !  --------------------------------------------------------
@@ -107,93 +108,49 @@
          END DO
       END DO
-
-
       !                                        !* Photosynthetically Available Radiation (PAR)
       !                                        !  --------------------------------------
-
-      IF( ln_varpar ) THEN
-         ze1(:,:,1) = par_varsw(:,:) * qsr(:,:) * EXP( -0.5 * zekb(:,:,1) )
-         ze2(:,:,1) = par_varsw(:,:) * qsr(:,:) * EXP( -0.5 * zekg(:,:,1) )
-         ze3(:,:,1) = par_varsw(:,:) * qsr(:,:) * EXP( -0.5 * zekr(:,:,1) )
+      IF( ln_dm2dc ) THEN                     !  diurnal cycle
+         ! 1% of qsr to compute euphotic layer
+         zqsr100(:,:) = 0.01 * qsr_mean(:,:)     !  daily mean qsr
+         !
+         CALL p4z_opt_par( kt, qsr_mean, zekb, zekg, zekr, ze1, ze2, ze3 ) 
+         !
+         DO jk = 1, nksrp      
+            etot_ndcy(:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk)
+            enano    (:,:,jk) =  2.1 * ze1(:,:,jk) + 0.42 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)
+            ediat    (:,:,jk) =  1.6 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.7 * ze3(:,:,jk)
+         END DO
+         !
+         CALL p4z_opt_par( kt, qsr, zekb, zekg, zekr, ze1, ze2, ze3 ) 
+         !
+         DO jk = 1, nksrp      
+            etot(:,:,jk) =  ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk)
+         END DO
+         !
       ELSE
-         ze1(:,:,1) = xparsw         * qsr(:,:) * EXP( -0.5 * zekb(:,:,1) )
-         ze2(:,:,1) = xparsw         * qsr(:,:) * EXP( -0.5 * zekg(:,:,1) )
-         ze3(:,:,1) = xparsw         * qsr(:,:) * EXP( -0.5 * zekr(:,:,1) )
-      ENDIF
-
-!CDIR NOVERRCHK
-      DO jj = 1, jpj
-!CDIR NOVERRCHK
-         DO ji = 1, jpi
-            zc1 = ze1(ji,jj,1)
-            zc2 = ze2(ji,jj,1) 
-            zc3 = ze3(ji,jj,1)
-            etot (ji,jj,1) = (       zc1 +        zc2 +       zc3 )
-            enano(ji,jj,1) = ( 2.1 * zc1 + 0.42 * zc2 + 0.4 * zc3 )
-            ediat(ji,jj,1) = ( 1.6 * zc1 + 0.69 * zc2 + 0.7 * zc3 )
-         END DO
-      END DO
-
-    
-      DO jk = 2, nksrp      
-!CDIR NOVERRCHK
-         DO jj = 1, jpj
-!CDIR NOVERRCHK
-            DO ji = 1, jpi
-               zc1 = ze1(ji,jj,jk-1) * EXP( -0.5 * ( zekb(ji,jj,jk-1) + zekb(ji,jj,jk) ) )
-               zc2 = ze2(ji,jj,jk-1) * EXP( -0.5 * ( zekg(ji,jj,jk-1) + zekg(ji,jj,jk) ) )
-               zc3 = ze3(ji,jj,jk-1) * EXP( -0.5 * ( zekr(ji,jj,jk-1) + zekr(ji,jj,jk) ) )
-               ze1  (ji,jj,jk) = zc1
-               ze2  (ji,jj,jk) = zc2
-               ze3  (ji,jj,jk) = zc3
-               etot (ji,jj,jk) = (       zc1 +        zc2 +       zc3 )
-               enano(ji,jj,jk) = ( 2.1 * zc1 + 0.42 * zc2 + 0.4 * zc3 )
-               ediat(ji,jj,jk) = ( 1.6 * zc1 + 0.69 * zc2 + 0.7 * zc3 )
-            END DO
-         END DO
-      END DO
+         ! 1% of qsr to compute euphotic layer
+         zqsr100(:,:) = 0.01 * qsr(:,:)
+         !
+         CALL p4z_opt_par( kt, qsr, zekb, zekg, zekr, ze1, ze2, ze3 ) 
+         !
+         DO jk = 1, nksrp      
+            etot (:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk)
+            enano(:,:,jk) =  2.1 * ze1(:,:,jk) + 0.42 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)
+            ediat(:,:,jk) =  1.6 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.7 * ze3(:,:,jk)
+         END DO
+         etot_ndcy(:,:,:) =  etot(:,:,:) 
+      ENDIF
+
 
       IF( ln_qsr_bio ) THEN                    !* heat flux accros w-level (used in the dynamics)
          !                                     !  ------------------------
-         zxsi0r = 1.e0 / rn_si0
-         !
-         ze0(:,:,1) = rn_abs * qsr(:,:)
-         !                                                    ! surface value : separation in R-G-B + near surface
-         IF( ln_varpar ) THEN
-            ze0(:,:,1) = ( 1. - 3. * par_varsw(:,:) ) * qsr(:,:)
-            ze1(:,:,1) = par_varsw(:,:)               * qsr(:,:)         
-            ze2(:,:,1) = par_varsw(:,:)               * qsr(:,:)
-            ze3(:,:,1) = par_varsw(:,:)               * qsr(:,:)
-         ELSE
-            ze0(:,:,1) = ( 1. - 3. * xparsw )  * qsr(:,:)
-            ze1(:,:,1) = xparsw                * qsr(:,:)         
-            ze2(:,:,1) = xparsw                * qsr(:,:)
-            ze3(:,:,1) = xparsw                * qsr(:,:)
-         ENDIF
+         CALL p4z_opt_par( kt, qsr, zekb, zekg, zekr, ze1, ze2, ze3, pe0=ze0 )
+         !
          etot3(:,:,1) =  qsr(:,:) * tmask(:,:,1)
-         !
-         !
          DO jk = 2, nksrp + 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) * zxsi0r )
-                  zc1 = ze1(ji,jj,jk-1) * EXP( -zekb(ji,jj,jk-1 ) )
-                  zc2 = ze2(ji,jj,jk-1) * EXP( -zekg(ji,jj,jk-1 ) )
-                  zc3 = ze3(ji,jj,jk-1) * EXP( -zekr(ji,jj,jk-1 ) )
-                  ze0(ji,jj,jk) = zc0
-                  ze1(ji,jj,jk) = zc1
-                  ze2(ji,jj,jk) = zc2
-                  ze3(ji,jj,jk) = zc3
-                  etot3(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk)
-              END DO
-              !
-            END DO
-            !
-        END DO
-        !
-      ENDIF
-
+            etot3(:,:,jk) =  ( ze0(:,:,jk) + ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk) ) * tmask(:,:,jk)
+         END DO
+         !                                     !  ------------------------
+      ENDIF
       !                                        !* Euphotic depth and level
       neln(:,:) = 1                            !  ------------------------
@@ -203,7 +160,7 @@
          DO jj = 1, jpj
            DO ji = 1, jpi
-              IF( etot(ji,jj,jk) * tmask(ji,jj,jk) >= 0.0043 * qsr(ji,jj) )  THEN
+              IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.43 * zqsr100(ji,jj) )  THEN
                  neln(ji,jj) = jk+1                    ! Euphotic level : 1rst T-level strictly below Euphotic layer
-                 !                                     ! nb: ensure the compatibility with nmld_trc definition in trd_mxl_trc_zint
+                 !                                     ! nb: ensure the compatibility with nmld_trc definition in trd_mld_trc_zint
                  heup(ji,jj) = fsdepw(ji,jj,jk+1)      ! Euphotic layer depth
               ENDIF
@@ -211,12 +168,12 @@
         END DO
       END DO
- 
+      !
       heup(:,:) = MIN( 300., heup(:,:) )
-
       !                                        !* mean light over the mixed layer
       zdepmoy(:,:)   = 0.e0                    !  -------------------------------
-      zetmp  (:,:)   = 0.e0
       zetmp1 (:,:)   = 0.e0
       zetmp2 (:,:)   = 0.e0
+      zetmp3 (:,:)   = 0.e0
+      zetmp4 (:,:)   = 0.e0
 
       DO jk = 1, nksrp
@@ -226,7 +183,8 @@
             DO ji = 1, jpi
                IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
-                  zetmp  (ji,jj) = zetmp  (ji,jj) + etot (ji,jj,jk) * fse3t(ji,jj,jk)
-                  zetmp1 (ji,jj) = zetmp1 (ji,jj) + enano(ji,jj,jk) * fse3t(ji,jj,jk)
-                  zetmp2 (ji,jj) = zetmp2 (ji,jj) + ediat(ji,jj,jk) * fse3t(ji,jj,jk)
+                  zetmp1 (ji,jj) = zetmp1 (ji,jj) + etot     (ji,jj,jk) * fse3t(ji,jj,jk) ! remineralisation
+                  zetmp2 (ji,jj) = zetmp2 (ji,jj) + etot_ndcy(ji,jj,jk) * fse3t(ji,jj,jk) ! production
+                  zetmp3 (ji,jj) = zetmp3 (ji,jj) + enano    (ji,jj,jk) * fse3t(ji,jj,jk) ! production
+                  zetmp4 (ji,jj) = zetmp4 (ji,jj) + ediat    (ji,jj,jk) * fse3t(ji,jj,jk) ! production
                   zdepmoy(ji,jj) = zdepmoy(ji,jj) + fse3t(ji,jj,jk)
                ENDIF
@@ -235,5 +193,6 @@
       END DO
       !
-      emoy(:,:,:) = etot(:,:,:)
+      emoy(:,:,:) = etot(:,:,:)       ! remineralisation
+      zpar(:,:,:) = etot_ndcy(:,:,:)  ! diagnostic : PAR with no diurnal cycle 
       !
       DO jk = 1, nksrp
@@ -244,12 +203,13 @@
                IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
                   z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
-                  emoy (ji,jj,jk) = zetmp (ji,jj) * z1_dep
-                  enano(ji,jj,jk) = zetmp1(ji,jj) * z1_dep
-                  ediat(ji,jj,jk) = zetmp2(ji,jj) * z1_dep
+                  emoy (ji,jj,jk) = zetmp1(ji,jj) * z1_dep
+                  zpar (ji,jj,jk) = zetmp2(ji,jj) * z1_dep
+                  enano(ji,jj,jk) = zetmp3(ji,jj) * z1_dep
+                  ediat(ji,jj,jk) = zetmp4(ji,jj) * z1_dep
                ENDIF
             END DO
          END DO
       END DO
-
+      !
       IF( lk_iomput ) THEN
         IF( jnt == nrdttrc  ) THEN
@@ -259,11 +219,11 @@
       ELSE
          IF( ln_diatrc ) THEN        ! save output diagnostics
-            trc2d(:,:,  jp_pcs0_2d + 10) = heup(:,:  ) * tmask(:,:,1)  
+            trc2d(:,:,  jp_pcs0_2d + 10) = heup(:,:  ) * tmask(:,:,1)
             trc3d(:,:,:,jp_pcs0_3d + 3)  = etot(:,:,:) * tmask(:,:,:)
          ENDIF
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj,      zdepmoy, zetmp, zetmp1, zetmp2 )
-      CALL wrk_dealloc( jpi, jpj, jpk, zekg, zekr, zekb, ze0, ze1, ze2, ze3 )
+      CALL wrk_dealloc( jpi, jpj,      zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 )
+      CALL wrk_dealloc( jpi, jpj, jpk, zpar, zekg, zekr, zekb, ze0, ze1, ze2, ze3 )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_opt')
@@ -271,7 +231,76 @@
    END SUBROUTINE p4z_opt
 
-   SUBROUTINE p4z_optsbc( kt )
-      !!----------------------------------------------------------------------
-      !!                  ***  routine p4z_optsbc  ***
+   SUBROUTINE p4z_opt_par( kt, pqsr, pekb, pekg, pekr, pe1, pe2, pe3, pe0 ) 
+      !!----------------------------------------------------------------------
+      !!                  ***  routine p4z_opt_par  ***
+      !!
+      !! ** purpose :   compute PAR of each wavelength (Red-Green-Blue)
+      !!                for a given shortwave radiation
+      !!
+      !!----------------------------------------------------------------------
+      !! * arguments
+      INTEGER, INTENT(in)                                       ::  kt            !   ocean time-step
+      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in)              ::  pqsr          !   shortwave
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)              ::  pekb, pekg, pekr   ! wavelength (Red-Green-Blue)
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)           ::  pe1 , pe2 , pe3   !  PAR ( R-G-B)
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout), OPTIONAL ::  pe0  
+      !! * local variables
+      INTEGER    ::   ji, jj, jk     ! dummy loop indices
+      REAL(wp), DIMENSION(jpi,jpj)     ::  zqsr          !   shortwave
+      !!----------------------------------------------------------------------
+
+      !  Real shortwave
+      IF( ln_varpar ) THEN  ;  zqsr(:,:) = par_varsw(:,:) * pqsr(:,:)
+      ELSE                  ;  zqsr(:,:) = xparsw         * pqsr(:,:)
+      ENDIF
+      !
+      IF( PRESENT( pe0 ) ) THEN     !  W-level
+         !
+         pe0(:,:,1) = pqsr(:,:) - 3. * zqsr(:,:)    !   ( 1 - 3 * alpha ) * q
+         pe1(:,:,1) = zqsr(:,:)         
+         pe2(:,:,1) = zqsr(:,:)
+         pe3(:,:,1) = zqsr(:,:)
+         !
+         DO jk = 2, nksrp + 1
+!CDIR NOVERRCHK
+            DO jj = 1, jpj
+!CDIR NOVERRCHK
+               DO ji = 1, jpi
+                  pe0(ji,jj,jk) = pe0(ji,jj,jk-1) * EXP( -fse3t(ji,jj,jk-1) * xsi0r )
+                  pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -pekb(ji,jj,jk-1 ) )
+                  pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -pekg(ji,jj,jk-1 ) )
+                  pe3(ji,jj,jk) = pe3(ji,jj,jk-1) * EXP( -pekr(ji,jj,jk-1 ) )
+               END DO
+              !
+            END DO
+            !
+         END DO
+        !
+      ELSE   ! T- level
+        !
+        pe1(:,:,1) = zqsr(:,:) * EXP( -0.5 * pekb(:,:,1) )
+        pe2(:,:,1) = zqsr(:,:) * EXP( -0.5 * pekg(:,:,1) )
+        pe3(:,:,1) = zqsr(:,:) * EXP( -0.5 * pekr(:,:,1) )
+        !
+        DO jk = 2, nksrp      
+!CDIR NOVERRCHK
+           DO jj = 1, jpj
+!CDIR NOVERRCHK
+              DO ji = 1, jpi
+                 pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -0.5 * ( pekb(ji,jj,jk-1) + pekb(ji,jj,jk) ) )
+                 pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -0.5 * ( pekg(ji,jj,jk-1) + pekg(ji,jj,jk) ) )
+                 pe3(ji,jj,jk) = pe3(ji,jj,jk-1) * EXP( -0.5 * ( pekr(ji,jj,jk-1) + pekr(ji,jj,jk) ) )
+              END DO
+           END DO
+        END DO    
+        !
+      ENDIF
+      ! 
+   END SUBROUTINE p4z_opt_par
+
+
+   SUBROUTINE p4z_opt_sbc( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  routine p4z_opt_sbc  ***
       !!
       !! ** purpose :   read and interpolate the variable PAR fraction
@@ -284,5 +313,5 @@
       !!----------------------------------------------------------------------
       !! * arguments
-      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
+      INTEGER ,                INTENT(in) ::   kt     ! ocean time step
 
       !! * local declarations
@@ -297,5 +326,5 @@
          IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_par > 1 ) ) THEN
             CALL fld_read( kt, 1, sf_par )
-            par_varsw(:,:) = ( sf_par(1)%fnow(:,:,1) )/3.0
+            par_varsw(:,:) = ( sf_par(1)%fnow(:,:,1) ) / 3.0
          ENDIF
       ENDIF
@@ -303,5 +332,5 @@
       IF( nn_timing == 1 )  CALL timing_stop('p4z_optsbc')
       !
-   END SUBROUTINE p4z_optsbc
+   END SUBROUTINE p4z_opt_sbc
 
    SUBROUTINE p4z_opt_init
@@ -347,4 +376,5 @@
       !
       xparsw = parlux / 3.0
+      xsi0r  = 1.e0 / rn_si0
       !
       ! Variable PAR at the surface of the ocean
@@ -372,8 +402,9 @@
       IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksrp, ' ref depth = ', gdepw_1d(nksrp+1), ' m'
       !
-                         etot (:,:,:) = 0._wp
-                         enano(:,:,:) = 0._wp
-                         ediat(:,:,:) = 0._wp
-      IF( ln_qsr_bio )   etot3(:,:,:) = 0._wp
+                         etot     (:,:,:) = 0._wp
+                         etot_ndcy(:,:,:) = 0._wp
+                         enano    (:,:,:) = 0._wp
+                         ediat    (:,:,:) = 0._wp
+      IF( ln_qsr_bio )   etot3    (:,:,:) = 0._wp
       ! 
       IF( nn_timing == 1 )  CALL timing_stop('p4z_opt_init')
@@ -386,5 +417,6 @@
       !!                     ***  ROUTINE p4z_opt_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( enano(jpi,jpj,jpk), ediat(jpi,jpj,jpk), emoy (jpi,jpj,jpk), STAT=p4z_opt_alloc ) 
+      ALLOCATE( enano    (jpi,jpj,jpk), ediat(jpi,jpj,jpk), &
+        &       etot_ndcy(jpi,jpj,jpk), emoy (jpi,jpj,jpk), STAT=p4z_opt_alloc ) 
          !
       IF( p4z_opt_alloc /= 0 ) CALL ctl_warn('p4z_opt_alloc : failed to allocate arrays.')
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 5206)
@@ -134,5 +134,5 @@
             DO jj = 1 ,jpj
                DO ji = 1, jpi
-                  IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+                  IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                      zval = MAX( 1., zstrn(ji,jj) )
                      zval = 1.5 * zval / ( 12. + zval )
@@ -157,5 +157,5 @@
                DO ji = 1, jpi
                   ! Computation of the P-I slope for nanos and diatoms
-                  IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+                  IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                       ztn         = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. )
                       zadap       = xadap * ztn / ( 2.+ ztn )
@@ -196,5 +196,5 @@
 
                   ! Computation of the P-I slope for nanos and diatoms
-                  IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+                  IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                       ztn         = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. )
                       zadap       = ztn / ( 2.+ ztn )
@@ -252,5 +252,5 @@
             DO ji = 1, jpi
 
-                IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+                IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                    !    Si/C of diatoms
                    !    ------------------------
@@ -302,5 +302,5 @@
 !CDIR NOVERRCHK
             DO ji = 1, jpi
-               IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                   !  production terms for nanophyto.
                   zprorca(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * trn(ji,jj,jk,jpphy) * rfact2
@@ -341,5 +341,5 @@
                      zprdch(ji,jj,jk) = zprdch(ji,jj,jk) * zmixdiat(ji,jj)
                   ENDIF
-                  IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+                  IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                      !  production terms for nanophyto. ( chlorophyll )
                      znanotot = enano(ji,jj,jk) * zstrn(ji,jj)
@@ -365,5 +365,5 @@
 !CDIR NOVERRCHK
                DO ji = 1, jpi
-                  IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+                  IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                      !  production terms for nanophyto. ( chlorophyll )
                      znanotot = enano(ji,jj,jk) * zstrn(ji,jj)
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 5206)
@@ -136,5 +136,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zcoef = ryyss * cvol(ji,jj,1) 
+                  zcoef = ryyss * e1e2t(ji,jj) * h_rnf(ji,jj) 
                   rivalk(ji,jj) =   sf_river(jr_dic)%fnow(ji,jj,1)                                    &
                      &              * 1.E3        / ( 12. * zcoef + rtrn )
@@ -187,4 +187,5 @@
       INTEGER  :: ierr, ierr1, ierr2, ierr3
       INTEGER  :: ios                 ! Local integer output status for namelist read
+      INTEGER  :: ik50               !  last level where depth less than 50 m
       REAL(wp) :: zexpide, zdenitide, zmaskt
       REAL(wp) :: ztimes_dust, ztimes_riv, ztimes_ndep 
@@ -247,4 +248,11 @@
       ENDIF
 
+      ! set the number of level over which river runoffs are applied 
+      ! online configuration : computed in sbcrnf
+      IF( lk_offline ) THEN
+        nk_rnf(:,:) = 1
+        h_rnf (:,:) = fsdept(:,:,1)
+      ENDIF
+
       ! dust input from the atmosphere
       ! ------------------------------
@@ -410,5 +418,12 @@
          CALL iom_close( numiron )
          !
-         DO jk = 1, 5
+         ik50 = 5        !  last level where depth less than 50 m
+         DO jk = jpkm1, 1, -1
+            IF( gdept_1d(jk) > 50. )  ik50 = jk - 1
+         END DO
+         IF (lwp) WRITE(numout,*)
+         IF (lwp) WRITE(numout,*) ' Level corresponding to 50m depth ',  ik50,' ', gdept_1d(ik50+1)
+         IF (lwp) WRITE(numout,*)
+         DO jk = 1, ik50
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
@@ -421,28 +436,66 @@
             END DO
          END DO
-         IF( cp_cfg == 'orca' .AND. jp_cfg == 2 ) THEN
-            ii0 = 176   ;   ii1 =  176        ! Southern Island : Kerguelen
-            ij0 =  37   ;   ij1 =   37  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+         IF( cp_cfg == 'orca' ) THEN
+            IF( jp_cfg == 2 ) THEN
+               ii0 = 176   ;   ii1 =  176        ! Southern Island : Kerguelen
+               ij0 =  37   ;   ij1 =   37  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+               !
+               ii0 = 119   ;   ii1 =  119        ! South Georgia
+               ij0 =  29   ;   ij1 =   29  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+               !
+               ii0 = 111   ;   ii1 =  111        ! Falklands
+               ij0 =  35   ;   ij1 =   35  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+               !
+               ii0 = 168   ;   ii1 =  168        ! Crozet
+               ij0 =  40   ;   ij1 =   40  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+               !
+               ii0 = 119   ;   ii1 =  119        ! South Orkney
+               ij0 =  28   ;   ij1 =   28  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+               !
+               ii0 = 140   ;   ii1 =  140        ! Bouvet Island
+               ij0 =  33   ;   ij1 =   33  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+               !
+               ii0 = 178   ;   ii1 =  178        ! Prince edwards
+               ij0 =  34   ;   ij1 =   34  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+               !
+               ii0 =  43   ;   ii1 =   43        ! Balleny islands
+               ij0 =  21   ;   ij1 =   21  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+               !
+            ELSE IF( jp_cfg == 1 ) THEN
+               !
+               ii0 = 357   ;   ii1 =  357        ! Southern Island : Kerguelen
+               ij0 =  75   ;   ij1 =   76  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp
+               !
+               ii0 = 243   ;   ii1 =  243        ! South Georgia
+               ij0 =  57   ;   ij1 =   59  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp
+               !
+               ii0 = 227   ;   ii1 =  227        ! Falklands
+               ij0 =  71   ;   ij1 =   73  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp
+               !
+               ii0 = 326   ;   ii1 =  327        ! Crozet
+               ij0 =  79   ;   ij1 =   79  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp
+               !
+               ii0 = 243   ;   ii1 =  243        ! South Orkney
+               ij0 =  56   ;   ij1 =   56  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp
+               !
+               ii0 = 283   ;   ii1 =  286        ! Bouvet Island
+               ij0 =  66   ;   ij1 =   66  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp
+               !
+               ii0 = 361   ;   ii1 =  361        ! Prince edwards
+               ij0 =  67   ;   ij1 =   68  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp
+               !
+               ii0 =  91   ;   ii1 =   92        ! Balleny islands
+               ij0 =  42   ;   ij1 =   43  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp
+               !
+            ELSE IF( jp_cfg == 05 ) THEN
+               !
+               CALL ctl_warn( ' Coastal supply of iron modifications for southern Islands in ORCA_R05' )
+               !
+            ELSE IF( jp_cfg == 025 ) THEN
+               !
+               CALL ctl_warn( ' Coastal supply of iron modifications for southern Islands in ORCA_R025' )
+               !
+            ENDIF
             !
-            ii0 = 119   ;   ii1 =  119        ! South Georgia
-            ij0 =  29   ;   ij1 =   29  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
-            !
-            ii0 = 111   ;   ii1 =  111        ! Falklands
-            ij0 =  35   ;   ij1 =   35  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
-            !
-            ii0 = 168   ;   ii1 =  168        ! Crozet
-            ij0 =  40   ;   ij1 =   40  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
-            !
-            ii0 = 119   ;   ii1 =  119        ! South Orkney
-            ij0 =  28   ;   ij1 =   28  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
-            !
-            ii0 = 140   ;   ii1 =  140        ! Bouvet Island
-            ij0 =  33   ;   ij1 =   33  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
-            !
-            ii0 = 178   ;   ii1 =  178        ! Prince edwards
-            ij0 =  34   ;   ij1 =   34  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
-            !
-            ii0 =  43   ;   ii1 =   43        ! Balleny islands
-            ij0 =  21   ;   ij1 =   21  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
          ENDIF
          CALL lbc_lnk( zcmask , 'T', 1. )      ! lateral boundary conditions on cmask   (sign unchanged)
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90	(revision 5206)
@@ -167,10 +167,16 @@
       ! ----------------------------------------------------------
       IF( ln_river ) THEN
-         trn(:,:,1,jppo4) = trn(:,:,1,jppo4) + rivdip(:,:) * rfact2
-         trn(:,:,1,jpno3) = trn(:,:,1,jpno3) + rivdin(:,:) * rfact2
-         trn(:,:,1,jpfer) = trn(:,:,1,jpfer) + rivdic(:,:) * 5.e-5 * rfact2
-         trn(:,:,1,jpsil) = trn(:,:,1,jpsil) + rivdsi(:,:) * rfact2
-         trn(:,:,1,jpdic) = trn(:,:,1,jpdic) + rivdic(:,:) * rfact2
-         trn(:,:,1,jptal) = trn(:,:,1,jptal) + ( rivalk(:,:) - rno3 * rivdin(:,:) ) * rfact2
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               DO jk = 1, nk_rnf(ji,jj)
+                  trn(ji,jj,jk,jppo4) = trn(ji,jj,jk,jppo4) +  rivdip(ji,jj) * rfact2
+                  trn(ji,jj,jk,jpno3) = trn(ji,jj,jk,jpno3) +  rivdin(ji,jj) * rfact2
+                  trn(ji,jj,jk,jpfer) = trn(ji,jj,jk,jpfer) +  rivdic(ji,jj) * 5.e-5 * rfact2
+                  trn(ji,jj,jk,jpsil) = trn(ji,jj,jk,jpsil) +  rivdsi(ji,jj) * rfact2
+                  trn(ji,jj,jk,jpdic) = trn(ji,jj,jk,jpdic) +  rivdic(ji,jj) * rfact2
+                  trn(ji,jj,jk,jptal) = trn(ji,jj,jk,jptal) +  ( rivalk(ji,jj) - rno3 * rivdin(ji,jj) ) * rfact2
+               ENDDO
+            ENDDO
+         ENDDO
       ENDIF
       
@@ -357,5 +363,5 @@
                ztrfer = biron(ji,jj,jk)       / ( concfediaz + biron(ji,jj,jk)       )
                ztrpo4 = trn  (ji,jj,jk,jppo4) / ( concnnh4   + trn  (ji,jj,jk,jppo4) ) 
-               zlight =  ( 1.- EXP( -etot(ji,jj,jk) / diazolight ) ) 
+               zlight =  ( 1.- EXP( -etot_ndcy(ji,jj,jk) / diazolight ) ) 
                znitrpot(ji,jj,jk) =  MAX( 0.e0, ( 0.6 * tgfunc(ji,jj,jk) - 2.15 ) * r1_rday )   &
                  &         *  zfact * MIN( ztrfer, ztrpo4 ) * zlight
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90	(revision 5206)
@@ -83,4 +83,5 @@
    USE sbc_oce , ONLY :   wndm       =>    wndm       !: 10m wind speed 
    USE sbc_oce , ONLY :   qsr        =>    qsr        !: penetrative solar radiation (w m-2)
+   USE sbc_oce , ONLY :   qsr_mean   =>    qsr_mean   !: daily mean solar heat flux
    USE sbc_oce , ONLY :   emp        =>    emp        !: freshwater budget: volume flux               [Kg/m2/s]
    USE sbc_oce , ONLY :   emp_b      =>    emp_b      !: freshwater budget: volume flux               [Kg/m2/s]
@@ -96,4 +97,5 @@
    USE sbcrnf  , ONLY :   rnfmsk_z   =>    rnfmsk_z   !: mixed adv scheme in runoffs vicinity (vert.)
    USE sbcrnf  , ONLY :   h_rnf      =>    h_rnf      !: river runoff   [Kg/m2/s]
+   USE sbcrnf  , ONLY :   nk_rnf     =>    nk_rnf     !: depth of runoff in model level
 
    USE trc_oce
@@ -135,4 +137,5 @@
 # endif
 
+   USE diaar5 , ONLY :   lk_diaar5  =>   lk_diaar5
 #else
    !!----------------------------------------------------------------------
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/trcini.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 5206)
@@ -71,9 +71,5 @@
       CALL top_alloc()              ! allocate TOP arrays
 
-#if defined key_offline
-      ltrcdm2dc = .FALSE.
-#endif
-
-      IF( ltrcdm2dc )CALL ctl_warn( ' Diurnal cycle on physics but not in PISCES or LOBSTER ' )
+      IF( lk_offline )  ln_dm2dc = .FALSE.
 
       IF( nn_cla == 1 )   &
Index: /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/trcstp.F90
===================================================================
--- /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/trcstp.F90	(revision 5205)
+++ /branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/TOP_SRC/trcstp.F90	(revision 5206)
@@ -69,12 +69,4 @@
       ENDIF
       !    
-      IF( ltrcdm2dc ) THEN
-         ! When Diurnal cycle, core bulk and LIM2  are activated, put daily mean qsr in qsr for TOP/biogeochemistery time-step
-         ! and save qsr with diurnal cycle in qsr_tmp
-         CALL wrk_alloc( jpi,jpj, zqsr_tmp )
-         zqsr_tmp(:,:) = qsr     (:,:)
-         qsr     (:,:) = qsr_mean(:,:)    
-      ENDIF
-      !    
       IF( nn_dttrc /= 1 )   CALL trc_sub_stp( kt )  ! averaging physical variables for sub-stepping
       !    
@@ -106,10 +98,4 @@
       ENDIF
       !
-      IF( ltrcdm2dc ) THEN
-         ! put back qsr with diurnal cycle in qsr
-         qsr(:,:) = zqsr_tmp(:,:)
-         CALL wrk_dealloc( jpi,jpj, zqsr_tmp )
-      ENDIF
-      !
       ztrai = 0._wp                                                   !  content of all tracers
       DO jn = 1, jptra
