Changeset 14990


Ignore:
Timestamp:
2021-06-14T21:19:43+02:00 (5 months ago)
Author:
techene
Message:

#2605 : new optimisation of traqsr

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/TRA/traqsr.F90

    r14618 r14990  
    1313   !!            3.6  !  2015-12  (O. Aumont, J. Jouanno, C. Ethe) use vertical profile of chlorophyll 
    1414   !!            3.7  !  2015-11  (G. Madec, A. Coward)  remove optimisation for fix volume 
     15   !!            4.0  !  2020-11  (A. Coward)  optimisation 
     16   !!            4.5  !  2021-03  (G. Madec)  further optimisation + adaptation for RK3 
    1517   !!---------------------------------------------------------------------- 
    1618 
    1719   !!---------------------------------------------------------------------- 
    1820   !!   tra_qsr       : temperature trend due to the penetration of solar radiation 
     21   !!       qsr_RGBc  : IR + RGB light penetration with Chlorophyll data case 
     22   !!       qsr_RGB   : IR + RGB light penetration with constant Chlorophyll case 
     23   !!       qsr_2BD   : 2 bands (InfraRed + Visible light) case 
     24   !!       qsr_ext_lev : level of extinction for each bands 
    1925   !!   tra_qsr_init  : initialization of the qsr penetration 
    2026   !!---------------------------------------------------------------------- 
     
    5359   REAL(wp), PUBLIC ::   rn_si1       !: deepest depth of extinction (water type I)       (2 bands) 
    5460   ! 
    55    INTEGER , PUBLIC ::   nksr         !: levels below which the light cannot penetrate (depth larger than 391 m) 
    56  
    5761   INTEGER, PARAMETER ::   np_RGB  = 1   ! R-G-B     light penetration with constant Chlorophyll 
    5862   INTEGER, PARAMETER ::   np_RGBc = 2   ! R-G-B     light penetration with Chlorophyll data 
     
    6064   INTEGER, PARAMETER ::   np_BIO  = 4   ! bio-model light penetration 
    6165   ! 
    62    INTEGER  ::   nqsr    ! user choice of the type of light penetration 
    63    REAL(wp) ::   xsi0r   ! inverse of rn_si0 
    64    REAL(wp) ::   xsi1r   ! inverse of rn_si1 
     66   INTEGER  ::   nqsr     ! user choice of the type of light penetration 
     67   INTEGER  ::   nc_rgb   ! RGB with cst Chlorophyll: index associated with the chosen Chl value 
     68   ! 
     69   !                       ! extinction level  
     70   INTEGER  ::   nk0             !: IR (depth larger ~12 m) 
     71   INTEGER  ::   nkV             !: Visible light (depth larger than ~840 m)  
     72   INTEGER  ::   nkR, nkG, nkB   !: RGB (depth larger than ~100 m, ~470 m, ~1700 m, resp.) 
     73   ! 
     74   INTEGER, PUBLIC  ::   nksr    !: =nkV, i.e. maximum level of light extinction (used in traatf(_qco).F90) 
     75   ! 
     76   !                  ! inverse of attenuation length 
     77   REAL(wp) ::   r1_si0                 ! all schemes : infrared  = 1/rn_si0  
     78   REAL(wp) ::   r1_si1                 ! 2 band      : mean RGB  = 1/rn_si1    
     79   REAL(wp) ::   r1_LR, r1_LG, r1_LB    ! RGB with constant Chl (np_RGB) 
    6580   ! 
    6681   REAL(wp) , PUBLIC, DIMENSION(3,61)   ::   rkrgb    ! tabulated attenuation coefficients for RGB absorption 
     
    7792CONTAINS 
    7893 
    79    SUBROUTINE tra_qsr( kt, Kmm, pts, Krhs, kstg ) 
     94   SUBROUTINE tra_qsr( kt, Kmm, pts, Krhs ) 
    8095      !!---------------------------------------------------------------------- 
    8196      !!                  ***  ROUTINE tra_qsr  *** 
     
    85100      !! 
    86101      !! ** Method  : The profile of the solar radiation within the ocean is defined 
    87       !!      through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) and a ratio rn_abs 
    88       !!      Considering the 2 wavebands case: 
    89       !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) ) 
    90       !!         The temperature trend associated with the solar radiation penetration 
    91       !!         is given by : zta = 1/e3t dk[ I ] / (rho0*Cp) 
    92       !!         At the bottom, boudary condition for the radiation is no flux : 
    93       !!      all heat which has not been absorbed in the above levels is put 
    94       !!      in the last ocean level. 
     102      !!      through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) or computed by 
     103      !!      the biogeochemical model 
    95104      !!         The computation is only done down to the level where 
    96       !!      I(k) < 1.e-15 W/m2 (i.e. over the top nksr levels) . 
    97       !! 
    98       !! ** Action  : - update ta with the penetrative solar radiation trend 
     105      !!      I(k) < 1.e-15 W/m2 (i.e. over the top nk levels) . 
     106      !! 
     107      !! ** Action  : - update ts(jp_tem) with the penetrative solar radiation trend 
    99108      !!              - send  trend for further diagnostics (l_trdtra=T) 
    100       !! 
    101       !! Reference  : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp. 
    102       !!              Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516. 
    103       !!              Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562 
    104       !!---------------------------------------------------------------------- 
    105       INTEGER,                                   INTENT(in   ) ::   kt, Kmm, Krhs   ! ocean time-step and time level indices 
     109      !!---------------------------------------------------------------------- 
     110      INTEGER,                                   INTENT(in   ) ::   kt, Kmm, Krhs   ! ocean time-step and time-level indices 
    106111      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) ::   pts             ! active tracers and RHS of tracer equation 
    107       INTEGER , OPTIONAL                       , INTENT(in   ) ::   kstg            ! RK3 stage index 
    108112      ! 
    109113      INTEGER  ::   ji, jj, jk               ! dummy loop indices 
    110       INTEGER  ::   irgb, isi, iei, isj, iej ! local integers 
    111       INTEGER  ::   istg_1, istg_3           !   -       - 
    112       REAL(wp) ::   zchl, zcoef, z1_2        ! local scalars 
    113       REAL(wp) ::   zc0 , zc1 , zc2 , zc3    !   -      - 
    114       REAL(wp) ::   zzc0, zzc1, zzc2, zzc3   !   -      - 
    115       REAL(wp) ::   zz0 , zz1 , ze3t, zlui   !   -      - 
    116       REAL(wp) ::   zCb, zCmax, zpsi, zpsimax, zrdpsi, zCze 
    117       REAL(wp) ::   zlogc, zlogze, zlogCtot, zlogCze 
    118       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   ze0, ze1, ze2, ze3 
    119       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, zetot, ztmp3d 
     114      REAL(wp) ::   z1_2, ze3t                     ! local scalars 
     115      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, zetot 
    120116      !!---------------------------------------------------------------------- 
    121117      ! 
    122118      IF( ln_timing )   CALL timing_start('tra_qsr') 
    123       ! 
    124       IF( PRESENT( kstg ) ) THEN      ! RK3 : a few things have to be done at only a specific stage 
    125          istg_1 = kstg   ;   istg_3 = kstg 
    126       ELSE                            ! MLF : only one call by time step 
    127          istg_1 =   1    ;   istg_3 =   3 
    128       ENDIF 
    129119      ! 
    130120      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
     
    136126      ENDIF 
    137127      ! 
    138       IF( l_trdtra ) THEN      ! trends diagnostic: save the input temperature trend 
     128      IF( l_trdtra ) THEN     ! trends diagnostic: save the input temperature trend 
    139129         ALLOCATE( ztrdt(jpi,jpj,jpk) ) 
    140130         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) 
    141131      ENDIF 
    142       ! 
    143       !                         !-----------------------------------! 
    144       !                         !  before qsr induced heat content  ! 
    145       !                         !-----------------------------------! 
    146       IF( ntsi == Nis0 ) THEN ; isi = nn_hls ; ELSE ; isi = 0 ; ENDIF    ! Avoid double-counting when using tiling 
    147       IF( ntsj == Njs0 ) THEN ; isj = nn_hls ; ELSE ; isj = 0 ; ENDIF 
    148       IF( ntei == Nie0 ) THEN ; iei = nn_hls ; ELSE ; iei = 0 ; ENDIF 
    149       IF( ntej == Nje0 ) THEN ; iej = nn_hls ; ELSE ; iej = 0 ; ENDIF 
    150  
    151       IF( kt == nit000 .AND. istg_1 == 1 ) THEN   !==  1st time step  ==! (RK3: only at stage 1) 
     132      !  
     133#if ! defined key_RK3 
     134      !                        ! MLF only : heat content trend due to Qsr flux (qsr_hc) 
     135      ! 
     136      IF( kt == nit000 ) THEN          !==  1st time step  ==! 
    152137         IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN    ! read in restart 
    153138            z1_2 = 0.5_wp 
     
    158143         ELSE                                           ! No restart or Euler forward at 1st time step 
    159144            z1_2 = 1._wp 
    160             DO_3D( isi, iei, isj, iej, 1, jpk ) 
     145            DO_3D( 0,0, 0,0, 1, jpk ) 
    161146               qsr_hc_b(ji,jj,jk) = 0._wp 
    162147            END_3D 
    163148         ENDIF 
    164       ELSEIF( istg_3 == 3 ) THEN                  !==  Swap of qsr heat content  ==! 
     149      ELSE                             !==  Swap of qsr heat content  ==! 
    165150         z1_2 = 0.5_wp 
    166          DO_3D( isi, iei, isj, iej, 1, jpk ) 
     151         DO_3D( 0,0, 0,0, 1, jpk ) 
    167152            qsr_hc_b(ji,jj,jk) = qsr_hc(ji,jj,jk) 
    168153         END_3D 
    169154      ENDIF 
    170       ! 
    171       !                         !--------------------------------! 
    172       SELECT CASE( nqsr )       !  now qsr induced heat content  ! 
    173       !                         !--------------------------------! 
    174       ! 
    175       CASE( np_BIO )                   !==  bio-model fluxes  ==! 
    176          ! 
    177          DO_3D( isi, iei, isj, iej, 1, nksr ) 
     155#endif 
     156       
     157      !                       !----------------------------! 
     158      SELECT CASE( nqsr )     !  qsr induced heat content  ! 
     159      !                       !----------------------------! 
     160      ! 
     161      CASE( np_RGBc )   ;   CALL qsr_RGBc( kt, Kmm, pts, Krhs )  !==  R-G-B fluxes using chlorophyll data     ==!    with Morel &Berthon (1989) vertical profile 
     162         ! 
     163      CASE( np_RGB  )   ;   CALL qsr_RGB ( kt, Kmm, pts, Krhs )  !==  R-G-B fluxes with constant chlorophyll  ==!    
     164         ! 
     165      CASE( np_2BD  )   ;   CALL qsr_2BD (     Kmm, pts, Krhs )  !==  2-bands fluxes                          ==! 
     166         ! 
     167      CASE( np_BIO )                                     !==  bio-model fluxes                        ==! 
     168         DO_3D( 0, 0, 0, 0, 1, nkV ) 
     169#if defined key_RK3 
     170            !                                                  !- RK3 : temperature trend at jk t-level 
     171            ze3t   = e3t(ji,jj,jk,Kmm) 
     172            pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) / ze3t 
     173#else 
     174            !                                                  !- MLF : heat content trend due to Qsr flux (qsr_hc) 
    178175            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) 
     176#endif 
    179177         END_3D 
    180          ! 
    181       CASE( np_RGB , np_RGBc )         !==  R-G-B fluxes  ==! 
    182          ! 
    183          ALLOCATE( ze0 (A2D(nn_hls))           , ze1 (A2D(nn_hls)) ,   & 
    184             &      ze2 (A2D(nn_hls))           , ze3 (A2D(nn_hls)) ,   & 
    185             &      ztmp3d(A2D(nn_hls),nksr + 1)                     ) 
    186          ! 
    187          IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll 
    188             IF( ntile == 0 .OR. ntile == 1 )  THEN                                         ! Do only for the full domain 
    189                IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )            ! Use full domain 
    190                CALL fld_read( kt, 1, sf_chl )         ! Read Chl data and provides it at the current time step 
    191                IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 1 )            ! Revert to tile domain 
    192             ENDIF 
    193             ! 
    194             ! Separation in R-G-B depending on the surface Chl 
    195             ! perform and store as many of the 2D calculations as possible 
    196             ! before the 3D loop (use the temporary 2D arrays to replace the 
    197             ! most expensive calculations) 
    198             ! 
    199             DO_2D( isi, iei, isj, iej ) 
    200                        ! zlogc = log(zchl) 
    201                zlogc = LOG ( MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) ) 
    202                        ! zc1 : log(zCze)  = log (1.12  * zchl**0.803) 
    203                zc1   = 0.113328685307 + 0.803 * zlogc 
    204                        ! zc2 : log(zCtot) = log(40.6  * zchl**0.459) 
    205                zc2   = 3.703768066608 + 0.459 * zlogc 
    206                        ! zc3 : log(zze)   = log(568.2 * zCtot**(-0.746)) 
    207                zc3   = 6.34247346942  - 0.746 * zc2 
    208                        ! IF( log(zze) > log(102.) ) log(zze) = log(200.0 * zCtot**(-0.293)) 
    209                IF( zc3 > 4.62497281328 ) zc3 = 5.298317366548 - 0.293 * zc2 
    210                ! 
    211                ze0(ji,jj) = zlogc                                                 ! ze0 = log(zchl) 
    212                ze1(ji,jj) = EXP( zc1 )                                            ! ze1 = zCze 
    213                ze2(ji,jj) = 1._wp / ( 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) ) ! ze2 = 1/zdelpsi 
    214                ze3(ji,jj) = EXP( - zc3 )                                          ! ze3 = 1/zze 
    215             END_2D 
    216             ! 
    217             DO_3D( isi, iei, isj, iej, 1, nksr + 1 ) 
    218                ! zchl    = ALOG( ze0(ji,jj) ) 
    219                zlogc = ze0(ji,jj) 
    220                ! 
    221                zCb       = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) ) 
    222                zCmax     = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 ) 
    223                zpsimax   = 0.6   - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) ) 
    224                ! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) 
    225                ! 
    226                zCze   = ze1(ji,jj) 
    227                zrdpsi = ze2(ji,jj)                                                 ! 1/zdelpsi 
    228                zpsi   = ze3(ji,jj) * gdepw(ji,jj,jk,Kmm)                           ! gdepw/zze 
    229                ! 
    230                ! NB. make sure zchl value is such that: zchl = MIN( 10. , MAX( 0.03, zchl ) ) 
    231                zchl = MIN( 10. , MAX( 0.03, zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) ) ) 
    232                ! Convert chlorophyll value to attenuation coefficient look-up table index 
    233                ztmp3d(ji,jj,jk) = 41 + 20.*LOG10(zchl) + 1.e-15 
    234             END_3D 
    235          ELSE                                !* constant chlorophyll 
    236             zchl = 0.05 
    237             ! NB. make sure constant value is such that: 
    238             zchl = MIN( 10. , MAX( 0.03, zchl ) ) 
    239             ! Convert chlorophyll value to attenuation coefficient look-up table index 
    240             zlui = 41 + 20.*LOG10(zchl) + 1.e-15 
    241             DO jk = 1, nksr + 1 
    242                ztmp3d(:,:,jk) = zlui 
     178         !                                                     !- sea-ice : store the 1st level attenuation coefficient 
     179         WHERE( etot3(A2D(0),1) /= 0._wp )   ;   fraqsr_1lev(A2D(0)) = 1._wp - etot3(A2D(0),2) / etot3(A2D(0),1) 
     180         ELSEWHERE                           ;   fraqsr_1lev(A2D(0)) = 1._wp 
     181         END WHERE 
     182         ! 
     183      END SELECT 
     184      ! 
     185#if defined key_RK3 
     186      !                             ! RK3 : diagnostics/output 
     187      IF( l_trdtra .OR. iom_use('qsr3d') ) THEN     ! qsr diagnostics 
     188         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 
     189         !                                         ! qsr tracers trends saved for diagnostics 
     190         IF( l_trdtra )   CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt ) 
     191         IF( iom_use('qsr3d') ) THEN               ! qsr distribution 
     192            DO jk = nkV, 1, -1 
     193               ztrdt(:,:,jk) = ztrdt(:,:,jk+1) + qsr_hc(:,:,jk) * rho0_rcp 
    243194            END DO 
     195            CALL iom_put( 'qsr3d', ztrdt )   ! 3D distribution of shortwave Radiation 
    244196         ENDIF 
    245          ! 
    246          zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B 
    247          DO_2D( isi, iei, isj, iej ) 
    248             ze0(ji,jj) = rn_abs * qsr(ji,jj) 
    249             ze1(ji,jj) = zcoef  * qsr(ji,jj) 
    250             ze2(ji,jj) = zcoef  * qsr(ji,jj) 
    251             ze3(ji,jj) = zcoef  * qsr(ji,jj) 
    252             ! store the surface SW radiation; re-use the surface ztmp3d array 
    253             ! since the surface attenuation coefficient is not used 
    254             ztmp3d(ji,jj,1) =       qsr(ji,jj) 
    255          END_2D 
    256          ! 
    257          !                                    !* interior equi-partition in R-G-B depending on vertical profile of Chl 
    258          DO_3D( isi, iei, isj, iej, 2, nksr + 1 ) 
    259             ze3t = e3t(ji,jj,jk-1,Kmm) 
    260             irgb = NINT( ztmp3d(ji,jj,jk) ) 
    261             zc0 = ze0(ji,jj) * EXP( - ze3t * xsi0r ) 
    262             zc1 = ze1(ji,jj) * EXP( - ze3t * rkrgb(1,irgb) ) 
    263             zc2 = ze2(ji,jj) * EXP( - ze3t * rkrgb(2,irgb) ) 
    264             zc3 = ze3(ji,jj) * EXP( - ze3t * rkrgb(3,irgb) ) 
    265             ze0(ji,jj) = zc0 
    266             ze1(ji,jj) = zc1 
    267             ze2(ji,jj) = zc2 
    268             ze3(ji,jj) = zc3 
    269             ztmp3d(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk) 
    270          END_3D 
    271          ! 
    272          DO_3D( isi, iei, isj, iej, 1, nksr )          !* now qsr induced heat content 
    273             qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( ztmp3d(ji,jj,jk) - ztmp3d(ji,jj,jk+1) ) 
    274          END_3D 
    275          ! 
    276          DEALLOCATE( ze0 , ze1 , ze2 , ze3 , ztmp3d ) 
    277          ! 
    278       CASE( np_2BD  )            !==  2-bands fluxes  ==! 
    279          ! 
    280          zz0 =        rn_abs   * r1_rho0_rcp      ! surface equi-partition in 2-bands 
    281          zz1 = ( 1. - rn_abs ) * r1_rho0_rcp 
    282          DO_3D( isi, iei, isj, iej, 1, nksr )          !* now qsr induced heat content 
    283             zc0 = zz0 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi1r ) 
    284             zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r ) 
    285             qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 
    286          END_3D 
    287          ! 
    288       END SELECT 
    289       ! 
    290       !                          !-----------------------------! 
    291       !                          !  update to the temp. trend  ! 
    292       !                          !-----------------------------! 
     197         DEALLOCATE( ztrdt ) 
     198      ENDIF 
     199#else 
     200      !                             ! MLF : add the temperature trend 
    293201      DO_3D( 0, 0, 0, 0, 1, nksr ) 
    294202         pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   & 
     
    297205      END_3D 
    298206      ! 
     207!!st7-2 
    299208      ! sea-ice: store the 1st ocean level attenuation coefficient 
    300       DO_2D( isi, iei, isj, iej ) 
     209      DO_2D( 0, 0, 0, 0 ) 
    301210         IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) ) 
    302211         ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp 
    303212         ENDIF 
    304213      END_2D 
     214!!st7-2 end 
    305215      ! 
    306216      ! TEMP: [tiling] This change not necessary and working array can use A2D(nn_hls) if using XIOS (subdomain support) 
     
    317227      ENDIF 
    318228      ! 
     229      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics 
     230         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 
     231         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt ) 
     232         DEALLOCATE( ztrdt ) 
     233      ENDIF 
     234#endif 
     235      ! 
    319236      IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only on the last tile 
    320237         IF( lrst_oce ) THEN     ! write in the ocean restart file 
     
    324241      ENDIF 
    325242      ! 
    326       IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics 
    327          ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 
    328          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt ) 
    329          DEALLOCATE( ztrdt ) 
    330       ENDIF 
     243!!st      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics 
     244!!st         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 
     245!!st         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt ) 
     246!!st         DEALLOCATE( ztrdt ) 
     247!!st      ENDIF 
    331248      !                       ! print mean trends (used for debugging) 
    332249      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' ) 
     
    335252      ! 
    336253   END SUBROUTINE tra_qsr 
     254 
     255 
     256   SUBROUTINE qsr_RGBc( kt, Kmm, pts, Krhs ) 
     257      !!---------------------------------------------------------------------- 
     258      !!                  ***  ROUTINE qsr_RGBc  *** 
     259      !! 
     260      !! ** Purpose :   Red-Green-Blue solar radiation using chlorophyll data 
     261      !! 
     262      !! ** Method  : The profile of the solar radiation within the ocean is defined 
     263      !!      through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) and a ratio rn_abs 
     264      !!      Considering the 2 wavebands case: 
     265      !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) ) 
     266      !!         The temperature trend associated with the solar radiation penetration 
     267      !!         is given by : zta = 1/e3t dk[ I ] / (rho0*Cp) 
     268      !!         At the bottom, boudary condition for the radiation is no flux : 
     269      !!      all heat which has not been absorbed in the above levels is put 
     270      !!      in the last ocean level. 
     271      !!         The computation is only done down to the level where 
     272      !!      I(k) < 1.e-15 W/m2 (i.e. over the top nk levels) . 
     273      !! 
     274      !! ** Action  : - update ta with the penetrative solar radiation trend 
     275      !!              - send  trend for further diagnostics (l_trdtra=T) 
     276      !! 
     277      !! Reference  : Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516. 
     278      !!              Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562 
     279      !!---------------------------------------------------------------------- 
     280      INTEGER,                                   INTENT(in   ) ::   kt, Kmm, Krhs   ! ocean time-step and time-level indices 
     281      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) ::   pts         ! active tracers and RHS of tracer equation 
     282      !! 
     283      INTEGER  ::   ji, jj, jk               ! dummy loop indices 
     284      INTEGER  ::   irgb                     ! local integer 
     285      REAL(wp) ::   zc1 , zc2 , zc3, zchl    ! local scalars 
     286      REAL(wp) ::   zze0, zzeR, zzeG, zzeB, zzeT              !    -         - 
     287      REAL(wp) ::   zz0 , zz1 , ze3t                          !    -         - 
     288      REAL(wp) ::   zCb, zCmax, zpsi, zpsimax, zrdpsi, zCze   !    -         - 
     289      REAL(wp) ::   zlogc, zlogze, zlogCtot, zlogCze          !    -         - 
     290      REAL(wp), DIMENSION(A2D(0)    ) ::   ze0, zeR, zeG, zeB, zeT 
     291      REAL(wp), DIMENSION(A2D(0),0:3) ::   zc 
     292      !!---------------------------------------------------------------------- 
     293      ! 
     294      ! 
     295      !                       !===========================================! 
     296      !                       !==  R-G-B fluxes using chlorophyll data  ==!    with Morel &Berthon (1989) vertical profile 
     297      !                       !===================================****====! 
     298      ! 
     299      !                             !=  Chlorophyll data  =! 
     300      ! 
     301      IF( ntile == 0 .OR. ntile == 1 )  THEN       ! Do only for the full domain 
     302         IF( ln_tile )   CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )   ! Use full domain 
     303         CALL fld_read( kt, 1, sf_chl )                                       ! Read Chl data and provides it at the current time step 
     304         IF( ln_tile )   CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 1 )   ! Revert to tile domain 
     305      ENDIF 
     306      ! 
     307      DO_2D( 0, 0, 0, 0 )                          ! pre-calculated expensive coefficient 
     308         zlogc = LOG(  MAX( 0.03_wp, MIN( sf_chl(1)%fnow(ji,jj,1) ,10._wp ) )  ) ! zlogc = log(zchl)   with 0.03 <= Chl >= 10.  
     309         zc1   = 0.113328685307 + 0.803 * zlogc                               ! zc1 : log(zCze)  = log (1.12  * zchl**0.803) 
     310         zc2   = 3.703768066608 + 0.459 * zlogc                               ! zc2 : log(zCtot) = log(40.6  * zchl**0.459) 
     311         zc3   = 6.34247346942  - 0.746 * zc2                                 ! zc3 : log(zze)   = log(568.2 * zCtot**(-0.746)) 
     312         IF( zc3 > 4.62497281328 )   zc3 = 5.298317366548 - 0.293 * zc2       ! IF(log(zze)>log(102)) log(zze) = log(200*zCtot**(-0.293)) 
     313         ! 
     314         zc(ji,jj,0) = zlogc                                                  ! ze(0) = log(zchl) 
     315         zc(ji,jj,1) = EXP( zc1 )                                             ! ze(1) = zCze 
     316         zc(ji,jj,2) = 1._wp / ( 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) )  ! ze(2) = 1/zdelpsi 
     317         zc(ji,jj,3) = EXP( - zc3 )                                           ! ze(3) = 1/zze 
     318      END_2D 
     319      ! 
     320      !                             !=  surface light  =! 
     321      ! 
     322      zz0 =           rn_abs              ! Infrared absorption 
     323      zz1 = ( 1._wp - rn_abs ) / 3._wp    ! R-G-B equi-partition 
     324      ! 
     325      DO_2D( 0, 0, 0, 0 )                 ! surface light 
     326         ze0(ji,jj) = zz0 * qsr(ji,jj)   ;   zeR(ji,jj) = zz1 * qsr(ji,jj)    ! IR    ; Red 
     327         zeG(ji,jj) = zz1 * qsr(ji,jj)   ;   zeB(ji,jj) = zz1 * qsr(ji,jj)    ! Green ; Blue 
     328         zeT(ji,jj) =       qsr(ji,jj)                                        ! Total 
     329      END_2D 
     330      !               
     331      !                             !=  interior light  =! 
     332      ! 
     333      DO jk = 1, nk0                      !* near surface layers *!   (< ~12 meters : IR + RGB ) 
     334         DO_2D( 0, 0, 0, 0 ) 
     335            !                                      !- inverse of RGB attenuation lengths 
     336            zlogc     = zc(ji,jj,0) 
     337            zCb       = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) ) 
     338            zCmax     = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 ) 
     339            zpsimax   = 0.6   - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) ) 
     340            ! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) 
     341            zCze   = zc(ji,jj,1) 
     342            zrdpsi = zc(ji,jj,2)                                     ! 1/zdelpsi 
     343!!st05            zpsi   = zc(ji,jj,3) * gdepw(ji,jj,jk,Kmm)               ! gdepw/zze 
     344            zpsi   = zc(ji,jj,3) * gdepw(ji,jj,jk+1,Kmm)               ! gdepw/zze 
     345            !                                                        ! make sure zchl value is such that: 0.03 < zchl < 10.  
     346            zchl = MAX(  0.03_wp , MIN( zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) , 10._wp )  ) 
     347            !                                                        ! Convert chlorophyll value to attenuation coefficient 
     348            irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )             ! look-up table index 
     349            !       Red             !         Green              !         Blue 
     350            r1_LR = rkrgb(3,irgb)   ;   r1_LG = rkrgb(2,irgb)    ;   r1_LB = rkrgb(1,irgb) 
     351            ! 
     352            !                                      !- fluxes at jk+1 w-level 
     353            ze3t = e3t(ji,jj,jk,Kmm) 
     354            zze0 = ze0(ji,jj) * EXP( - ze3t*r1_si0 )   ;   zzeR = zeR(ji,jj) * EXP( - ze3t*r1_LR )   ! IR    ; Red  at jk+1 w-level 
     355            zzeG = zeG(ji,jj) * EXP( - ze3t*r1_LG  )   ;   zzeB = zeB(ji,jj) * EXP( - ze3t*r1_LB )   ! Green ; Blue      -      - 
     356            zzeT = ( zze0 + zzeB + zzeG + zzeR ) * wmask(ji,jj,jk+1)                                 ! Total             -      - 
     357!!st01            zzeT = ( zze0 + zzeR + zzeG + zzeB ) * wmask(ji,jj,jk+1)                                 ! Total             -      - 
     358            ! 
     359#if defined key_RK3 
     360            !                                      !- RK3 : temperature trend at jk t-level 
     361            pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t 
     362#else 
     363            !                                      !- MLF : heat content trend due to Qsr flux (qsr_hc) 
     364            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) 
     365#endif 
     366            ze0(ji,jj) = zze0   ;   zeR(ji,jj) = zzeR           ! IR    ; Red  store at jk+1 w-level 
     367            zeG(ji,jj) = zzeG   ;   zeB(ji,jj) = zzeB           ! Green ; Blue   -        -      - 
     368            zeT(ji,jj) = zzeT                                   ! total          -        -      - 
     369         END_2D 
     370!!stbug         IF( jk == 1 ) THEN                        !- sea-ice : store the 1st level attenuation coefficient 
     371!!stbug            WHERE( qsr(A2D(0)) /= 0._wp )   ;   fraqsr_1lev(A2D(0)) = 1._wp - zeT(A2D(0)) / qsr(A2D(0)) 
     372!!stbug            ELSEWHERE                       ;   fraqsr_1lev(A2D(0)) = 1._wp 
     373!!stbug            END WHERE 
     374!!stbug         ENDIF 
     375      END DO 
     376      ! 
     377      IF( kt == nit000 ) THEN 
     378         IF(lwp) WRITE(numout,*) 'nk0+1= ', nk0+1, ' qsr IR max = ' , MAXVAL(ze0(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(ze0(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     379         IF(lwp) WRITE(numout,*) '       ', nk0+1, ' qsr R max = ' , MAXVAL(zeR(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeR(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     380         IF(lwp) WRITE(numout,*) '       ', nk0+1, ' qsr G max = ' , MAXVAL(zeG(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeG(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     381         IF(lwp) WRITE(numout,*) '       ', nk0+1, ' qsr B max = ' , MAXVAL(zeB(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeB(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     382         IF(lwp) WRITE(numout,*) '       ', nk0+1, ' qsr T max = ' , MAXVAL(zeT(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeT(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     383      ENDIF 
     384      ! 
     385      DO jk = nk0+1, nkR                  !* down to Red extinction *!   (< ~71 meters : RGB , IR removed from calculation) 
     386         DO_2D( 0, 0, 0, 0 ) 
     387            !                                      !- inverse of RGB attenuation lengths 
     388            zlogc     = zc(ji,jj,0) 
     389            zCb       = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) ) 
     390            zCmax     = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 ) 
     391            zpsimax   = 0.6   - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) ) 
     392            ! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) 
     393            zCze   = zc(ji,jj,1) 
     394            zrdpsi = zc(ji,jj,2)                               ! 1/zdelpsi 
     395            zpsi   = zc(ji,jj,3) * gdepw(ji,jj,jk+1,Kmm)         ! gdepw/zze 
     396!!st05            zpsi   = zc(ji,jj,3) * gdepw(ji,jj,jk,Kmm)         ! gdepw/zze 
     397            !                                                  ! make sure zchl value is such that: 0.03 < zchl < 10.  
     398            zchl = MAX(  0.03_wp , MIN( zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) , 10._wp )  ) 
     399            !                                                  ! Convert chlorophyll value to attenuation coefficient 
     400            irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )       ! look-up table index 
     401            !       Red             !         Green              !         Blue 
     402            r1_LR = rkrgb(3,irgb)   ;   r1_LG = rkrgb(2,irgb)    ;   r1_LB = rkrgb(1,irgb) 
     403            ! 
     404            !                                      !- fluxes at jk+1 w-level 
     405            ze3t = e3t(ji,jj,jk,Kmm) 
     406            zzeR = zeR(ji,jj) * EXP( - ze3t*r1_LR )                                                 ! Red          at jk+1 w-level 
     407            zzeG = zeG(ji,jj) * EXP( - ze3t*r1_LG )   ;   zzeB = zeB(ji,jj) * EXP( - ze3t*r1_LB )   ! Green ; Blue      -      - 
     408            zzeT = ( zzeR + zzeG + zzeB ) * wmask(ji,jj,jk+1)                                       ! Total             -      - 
     409            ! 
     410#if defined key_RK3 
     411            !                                      !- RK3 : temperature trend at jk t-level 
     412            pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t 
     413#else 
     414            !                                      !- MLF : heat content trend due to Qsr flux (qsr_hc) 
     415            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) 
     416#endif 
     417            zeR(ji,jj) = zzeR                                  ! Red          store at jk+1 w-level 
     418            zeG(ji,jj) = zzeG   ;   zeB(ji,jj) = zzeB          ! Green ; Blue   -        -      - 
     419            zeT(ji,jj) = zzeT                                  ! total          -        -      - 
     420         END_2D 
     421      END DO 
     422      ! 
     423      IF( kt == nit000 ) THEN 
     424         IF(lwp) WRITE(numout,*) 'nkR+1= ', nkR+1, ' qsr R max = ' , MAXVAL(zeR(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeR(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     425         IF(lwp) WRITE(numout,*) '       ', nkR+1, ' qsr G max = ' , MAXVAL(zeG(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeG(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     426         IF(lwp) WRITE(numout,*) '       ', nkR+1, ' qsr B max = ' , MAXVAL(zeB(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeB(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     427         IF(lwp) WRITE(numout,*) '       ', nkR+1, ' qsr T max = ' , MAXVAL(zeT(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeT(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     428      ENDIF 
     429      ! 
     430      DO jk = nkR+1, nkG                  !* down to Green extinction *!   (< ~350 m : GB , IR+R removed from calculation) 
     431         DO_2D( 0, 0, 0, 0 ) 
     432            !                                      !- inverse of RGB attenuation lengths 
     433            zlogc     = zc(ji,jj,0) 
     434            zCb       = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) ) 
     435            zCmax     = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 ) 
     436            zpsimax   = 0.6   - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) ) 
     437            ! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) 
     438            zCze   = zc(ji,jj,1) 
     439            zrdpsi = zc(ji,jj,2)                               ! 1/zdelpsi 
     440            zpsi   = zc(ji,jj,3) * gdepw(ji,jj,jk+1,Kmm)         ! gdepw/zze 
     441!!st05            zpsi   = zc(ji,jj,3) * gdepw(ji,jj,jk,Kmm)         ! gdepw/zze 
     442            !                                                  ! make sure zchl value is such that: 0.03 < zchl < 10.  
     443            zchl = MAX(  0.03_wp , MIN( zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) , 10._wp )  ) 
     444            !                                                  ! Convert chlorophyll value to attenuation coefficient 
     445            irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )       ! look-up table index 
     446            !     Green              !         Blue 
     447            r1_LG = rkrgb(2,irgb)    ;   r1_LB = rkrgb(1,irgb) 
     448            ! 
     449            !                                      !- fluxes at jk+1 w-level 
     450            ze3t = e3t(ji,jj,jk,Kmm) 
     451            zzeG = zeG(ji,jj) * EXP( - ze3t * r1_LG )   ;   zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB ) ! Green ; Blue 
     452            zzeT = ( zzeG + zzeB ) * wmask(ji,jj,jk+1)                                                ! Total             -      - 
     453#if defined key_RK3 
     454            !                                      !- RK3 : temperature trend at jk t-level 
     455            pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t 
     456#else 
     457            !                                      !- MLF : heat content trend due to Qsr flux (qsr_hc) 
     458            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) 
     459#endif 
     460            zeG(ji,jj) = zzeG   ;   zeB(ji,jj) = zzeB          ! Green ; Blue store at jk+1 w-level 
     461            zeT(ji,jj) = zzeT                                  ! total          -        -      - 
     462         END_2D 
     463      END DO 
     464      ! 
     465      IF( kt == nit000 ) THEN 
     466         IF(lwp) WRITE(numout,*) 'nkG+1= ', nkG+1, ' qsr G max = ' , MAXVAL(zeG(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeG(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     467         IF(lwp) WRITE(numout,*) '       ', nkG+1, ' qsr B max = ' , MAXVAL(zeB(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeB(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     468         IF(lwp) WRITE(numout,*) '       ', nkG+1, ' qsr T max = ' , MAXVAL(zeT(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeT(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     469      ENDIF 
     470      ! 
     471      DO jk = nkG+1, nkB                  !* down to Blue extinction *!   (< ~1300 m : B , IR+RG removed from calculation) 
     472         DO_2D( 0, 0, 0, 0 ) 
     473            !                                      !- inverse of RGB attenuation lengths 
     474            zlogc     = zc(ji,jj,0) 
     475            zCb       = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) ) 
     476            zCmax     = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 ) 
     477            zpsimax   = 0.6   - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) ) 
     478            ! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) 
     479            zCze   = zc(ji,jj,1) 
     480            zrdpsi = zc(ji,jj,2)                               ! 1/zdelpsi 
     481            zpsi   = zc(ji,jj,3) * gdepw(ji,jj,jk+1,Kmm)         ! gdepw/zze 
     482!!st05            zpsi   = zc(ji,jj,3) * gdepw(ji,jj,jk,Kmm)         ! gdepw/zze 
     483            !                                                  ! make sure zchl value is such that: 0.03 < zchl < 10.  
     484            zchl = MAX(  0.03_wp , MIN( zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) , 10._wp )  ) 
     485            !                                                  ! Convert chlorophyll value to attenuation coefficient 
     486            irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )       ! look-up table index 
     487            r1_LB = rkrgb(1,irgb)                              ! Blue 
     488            ! 
     489            !                                      !- fluxes at jk+1 w-level 
     490            ze3t = e3t(ji,jj,jk,Kmm) 
     491            zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB )          ! Blue 
     492            zzeT = ( zzeB ) * wmask(ji,jj,jk+1)                ! Total             -      - 
     493#if defined key_RK3 
     494            !                                      !- RK3 : temperature trend at jk t-level 
     495            pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t 
     496#else 
     497            !                                      !- MLF : heat content trend due to Qsr flux (qsr_hc) 
     498            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) 
     499#endif 
     500            zeB(ji,jj) = zzeB                                  ! Blue store at jk+1 w-level 
     501            zeT(ji,jj) = zzeT                                  ! total  -        -      - 
     502         END_2D 
     503      END DO 
     504      ! 
     505      IF( kt == nit000 ) THEN 
     506         IF(lwp) WRITE(numout,*) 'nkB+1= ', nkB+1, ' qsr T max = ' , MAXVAL(zeT), ' W/m2' , MAXVAL(zeT(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)), ' K/s' 
     507      ENDIF 
     508      ! 
     509   END SUBROUTINE qsr_RGBc 
     510 
     511 
     512   SUBROUTINE qsr_RGB( kt, Kmm, pts, Krhs ) 
     513      !!---------------------------------------------------------------------- 
     514      !!                  ***  ROUTINE qsr_RGB  *** 
     515      !! 
     516      !! ** Purpose :   Red-Green-Blue solar radiation with constant chlorophyll 
     517      !! 
     518      !! ** Method  : The profile of the solar radiation within the ocean is defined 
     519      !!      through 2 wavebands (rn_si0,rn_si1) or 1 (rn_si0,rn_abs) + 3 wavebands (RGB)  
     520      !!         At the bottom, boudary condition for the radiation is no flux : 
     521      !!      all heat which has not been absorbed in the above levels is put 
     522      !!      in the last ocean level. 
     523      !!         For each band, the computation is only done down to the level where 
     524      !!      I(k) < 1.e-15 W/m2 (i.e. over the top nk levels) . 
     525      !! 
     526      !! ** Action  : - update ta with the penetrative solar radiation trend 
     527      !!              - send  trend for further diagnostics (l_trdtra=T) 
     528      !! 
     529      !! Reference  : Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516. 
     530      !!              Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562 
     531      !!---------------------------------------------------------------------- 
     532      INTEGER,                                   INTENT(in   ) ::   kt, Kmm, Krhs   ! ocean time-step and time-level indices 
     533      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) ::   pts             ! active tracers and RHS of tracer equation 
     534      !! 
     535      INTEGER  ::   ji, jj, jk               ! dummy loop indices 
     536      REAL(wp) ::   zze0, zzeR, zzeG, zzeB, zzeT              !    -         - 
     537      REAL(wp) ::   zz0 , zz1 , ze3t                          !    -         - 
     538      REAL(wp), DIMENSION(A2D(0))   ::   ze0, zeR, zeG, zeB, zeT 
     539      !!---------------------------------------------------------------------- 
     540      !       
     541      ! 
     542      !                       !==============================================! 
     543      !                       !==  R-G-B fluxes with constant chlorophyll  ==!    
     544      !                       !======================********================! 
     545      ! 
     546      !                             !=  surface light  =! 
     547      ! 
     548      zz0 =           rn_abs              ! Infrared absorption 
     549      zz1 = ( 1._wp - rn_abs ) / 3._wp    ! surface equi-partition in R-G-B 
     550      ! 
     551      DO_2D( 0, 0, 0, 0 )                 ! surface light 
     552         ze0(ji,jj) = zz0 * qsr(ji,jj)   ;   zeR(ji,jj) = zz1 * qsr(ji,jj)    ! IR    ; Red 
     553         zeG(ji,jj) = zz1 * qsr(ji,jj)   ;   zeB(ji,jj) = zz1 * qsr(ji,jj)    ! Green ; Blue 
     554         zeT(ji,jj) =       qsr(ji,jj)                                        ! Total 
     555      END_2D 
     556      ! 
     557      !                             !=  interior light  =! 
     558      ! 
     559      DO jk = 1, nk0                      !* near surface layers *!   (< ~12 meters : IR + RGB ) 
     560         DO_2D( 0, 0, 0, 0 ) 
     561            ze3t = e3t(ji,jj,jk,Kmm) 
     562            zze0 = ze0(ji,jj) * EXP( - ze3t * r1_si0 )   ;   zzeR = zeR(ji,jj) * EXP( - ze3t * r1_LR )   ! IR    ; Red  at jk+1 w-level 
     563            zzeG = zeG(ji,jj) * EXP( - ze3t * r1_LG  )   ;   zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB )   ! Green ; Blue      -      - 
     564            zzeT = ( zze0 + zzeB + zzeG + zzeR ) * wmask(ji,jj,jk+1)                                     ! Total             -      - 
     565!!st7-9            zzeT = ( zze0 + zzeR + zzeG + zzeB ) * wmask(ji,jj,jk+1)                                     ! Total             -      - 
     566#if defined key_RK3 
     567            !                                               ! RK3 : temperature trend at jk t-level 
     568            pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t 
     569#else 
     570            !                                               ! MLF : heat content trend due to Qsr flux (qsr_hc) 
     571            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) 
     572#endif 
     573            ze0(ji,jj) = zze0   ;   zeR(ji,jj) = zzeR           ! IR    ; Red  store at jk+1 w-level 
     574            zeG(ji,jj) = zzeG   ;   zeB(ji,jj) = zzeB           ! Green ; Blue   -        -      - 
     575            zeT(ji,jj) = zzeT                                   ! total          -        -      - 
     576         END_2D 
     577!!stbug         IF( jk == 1 ) THEN               !* sea-ice *!   store the 1st level attenuation coeff. 
     578!!stbug            WHERE( qsr(A2D(0)) /= 0._wp )   ;   fraqsr_1lev(A2D(0)) = 1._wp - zeT(A2D(0)) / qsr(A2D(0)) 
     579!!stbug            ELSEWHERE                       ;   fraqsr_1lev(A2D(0)) = 1._wp 
     580!!stbug            END WHERE 
     581!!stbug         ENDIF 
     582      END DO 
     583      ! 
     584      IF( kt == nit000 ) THEN 
     585         IF(lwp) WRITE(numout,*) 'nk0+1= ', nk0+1, ' qsr IR max = ' , MAXVAL(ze0(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(ze0(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     586         IF(lwp) WRITE(numout,*) '       ', nk0+1, ' qsr R max = ' , MAXVAL(zeR(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeR(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     587         IF(lwp) WRITE(numout,*) '       ', nk0+1, ' qsr G max = ' , MAXVAL(zeG(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeG(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     588         IF(lwp) WRITE(numout,*) '       ', nk0+1, ' qsr B max = ' , MAXVAL(zeB(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeB(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     589         IF(lwp) WRITE(numout,*) '       ', nk0+1, ' qsr T max = ' , MAXVAL(zeT(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeT(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     590      ENDIF 
     591      ! 
     592      DO jk = nk0+1, nkR                  !* down to Red extinction *!   (< ~71 meters : RGB , IR removed from calculation) 
     593         DO_2D( 0, 0, 0, 0 ) 
     594            ze3t = e3t(ji,jj,jk,Kmm) 
     595            zzeR = zeR(ji,jj) * EXP( - ze3t * r1_LR )                                                 ! Red          at jk+1 w-level 
     596            zzeG = zeG(ji,jj) * EXP( - ze3t * r1_LG )   ;   zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB ) ! Green ; Blue      -      - 
     597            zzeT = ( zzeB + zzeG + zzeR ) * wmask(ji,jj,jk+1)                                         ! Total             -      - 
     598!!st7-11            zzeT = ( zzeR + zzeG + zzeB ) * wmask(ji,jj,jk+1)                                         ! Total             -      - 
     599#if defined key_RK3 
     600            !                                               ! RK3 : temperature trend at jk t-level 
     601            pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t 
     602#else 
     603            !                                               ! MLF : heat content trend due to Qsr flux (qsr_hc) 
     604            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) 
     605#endif 
     606            zeR(ji,jj) = zzeR                                   ! Red          store at jk+1 w-level 
     607            zeG(ji,jj) = zzeG   ;   zeB(ji,jj) = zzeB           ! Green ; Blue   -        -      - 
     608            zeT(ji,jj) = zzeT                                   ! total          -        -      - 
     609         END_2D 
     610      END DO 
     611      ! 
     612      IF( kt == nit000 ) THEN 
     613         IF(lwp) WRITE(numout,*) 'nkR+1= ', nkR+1, ' qsr R max = ' , MAXVAL(zeR(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeR(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     614         IF(lwp) WRITE(numout,*) '       ', nkR+1, ' qsr G max = ' , MAXVAL(zeG(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeG(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     615         IF(lwp) WRITE(numout,*) '       ', nkR+1, ' qsr B max = ' , MAXVAL(zeB(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeB(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     616         IF(lwp) WRITE(numout,*) '       ', nkR+1, ' qsr T max = ' , MAXVAL(zeT(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeT(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     617      ENDIF 
     618      ! 
     619      DO jk = nkR+1, nkG                  !* down to Green extinction *!   (< ~350 m : GB , IR+R removed from calculation) 
     620         DO_2D( 0, 0, 0, 0 ) 
     621            ze3t = e3t(ji,jj,jk,Kmm) 
     622            zzeG = zeG(ji,jj) * EXP( - ze3t * r1_LG )   ;   zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB ) ! Green ; Blue at jk+1 w-level 
     623            zzeT = ( zzeG + zzeB ) * wmask(ji,jj,jk+1)                                                ! Total             -      - 
     624#if defined key_RK3 
     625            !                                               ! RK3 : temperature trend at jk t-level 
     626            pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t 
     627#else 
     628            !                                               ! MLF : heat content trend due to Qsr flux (qsr_hc) 
     629            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) 
     630#endif 
     631            zeG(ji,jj) = zzeG   ;   zeB(ji,jj) = zzeB             ! Green ; Blue store at jk+1 w-level 
     632            zeT(ji,jj) = zzeT                                     ! total          -        -      - 
     633         END_2D 
     634      END DO 
     635      ! 
     636      IF( kt == nit000 ) THEN 
     637         IF(lwp) WRITE(numout,*) 'nkG+1= ', nkG+1, ' qsr G max = ' , MAXVAL(zeG(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeG(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     638         IF(lwp) WRITE(numout,*) '       ', nkG+1, ' qsr B max = ' , MAXVAL(zeB(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeB(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     639         IF(lwp) WRITE(numout,*) '       ', nkG+1, ' qsr T max = ' , MAXVAL(zeT(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeT(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 
     640      ENDIF 
     641      ! 
     642      DO jk = nkG+1, nkB                  !* down to Blue extinction *!   (< ~1300 m : B , IR+RG removed from calculation) 
     643         DO_2D( 0, 0, 0, 0 ) 
     644            ze3t = e3t(ji,jj,jk,Kmm) 
     645            zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB )             ! Blue at jk+1 w-level 
     646            zzeT = ( zzeB ) * wmask(ji,jj,jk+1)                   ! Total     -      - 
     647#if defined key_RK3 
     648            !                                               ! RK3 : temperature trend at jk t-level 
     649            pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t 
     650#else 
     651            !                                               ! MLF : heat content trend due to Qsr flux (qsr_hc) 
     652            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) 
     653#endif 
     654            zeB(ji,jj) = zzeB                                     ! Blue store at jk+1 w-level 
     655            zeT(ji,jj) = zzeT                                     ! total  -        -      - 
     656         END_2D 
     657      END DO 
     658      ! 
     659      IF( kt == nit000 ) THEN 
     660         IF(lwp) WRITE(numout,*) 'nkB+1= ', nkB+1, ' qsr T max = ' , MAXVAL(zeT), ' W/m2' , MAXVAL(zeT(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)), ' K/s' 
     661      ENDIF 
     662      ! 
     663   END SUBROUTINE qsr_RGB 
     664 
     665 
     666   SUBROUTINE qsr_2BD( Kmm, pts, Krhs ) 
     667      !!---------------------------------------------------------------------- 
     668      !!                  ***  ROUTINE qsr_2BD  *** 
     669      !! 
     670      !! ** Purpose :   2 bands (IR+visible) solar radiation with constant chlorophyll 
     671      !! 
     672      !! ** Method  : The profile of the solar radiation within the ocean is defined 
     673      !!      through 2 wavebands (rn_si0,rn_si1) a ratio rn_abs for IR absorbtion. 
     674      !!      Considering the 2 wavebands case: 
     675      !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) ) 
     676      !!         The temperature trend associated with the solar radiation penetration 
     677      !!         is given by : zta = 1/e3t dk[ I ] / (rho0*Cp) 
     678      !!         At the bottom, boudary condition for the radiation is no flux : 
     679      !!      all heat which has not been absorbed in the above levels is put 
     680      !!      in the last ocean level. 
     681      !!         The computation is only done down to the level where 
     682      !!      I(k) < 1.e-15 W/m2 (i.e. over the top nk levels) . 
     683      !! 
     684      !! ** Action  : - update ta with the penetrative solar radiation trend 
     685      !!              - send  trend for further diagnostics (l_trdtra=T) 
     686      !! 
     687      !! Reference  : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp. 
     688      !!              Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516. 
     689      !!              Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562 
     690      !!---------------------------------------------------------------------- 
     691      INTEGER,                                   INTENT(in   ) ::   Kmm, Krhs   ! ocean time-step and time-level indices 
     692      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) ::   pts         ! active tracers and RHS of tracer equation 
     693      !! 
     694      INTEGER  ::   ji, jj, jk               ! dummy loop indices 
     695      REAL(wp) ::   zzatt                    !    -         - 
     696      REAL(wp) ::   zz0 , zz1 , ze3t         !    -         - 
     697      REAL(wp), DIMENSION(A2D(0)) ::   zatt 
     698      !!---------------------------------------------------------------------- 
     699      !       
     700      !                       !======================! 
     701      !                       !==  2-bands fluxes  ==! 
     702      !                       !======================! 
     703      ! 
     704      zz0 =           rn_abs   * r1_rho0_rcp       ! surface equi-partition in 2-bands 
     705      zz1 = ( 1._wp - rn_abs ) * r1_rho0_rcp 
     706      ! 
     707      zatt(A2D(0)) = r1_rho0_rcp                   !* surface value *! 
     708      ! 
     709      DO_2D( 0, 0, 0, 0 ) 
     710         zatt(ji,jj) = (  zz0 * EXP( -gdepw(ji,jj,1,Kmm)*r1_si0 ) + zz1 * EXP( -gdepw(ji,jj,1,Kmm)*r1_si1 )  ) 
     711      END_2D 
     712      ! 
     713!!st      IF(lwp) WRITE(numout,*) 'level = ', 1, ' qsr max = ' , MAXVAL(zatt)*rho0_rcp, ' W/m2', ' qsr min = ' , MINVAL(zatt)*rho0_rcp, ' W/m2'  
     714      ! 
     715      DO jk = 1, nk0                               !* near surface layers *!   (< ~14 meters : IR + visible light ) 
     716         DO_2D( 0, 0, 0, 0 ) 
     717            ze3t  = e3t(ji,jj,jk,Kmm)                    ! light attenuation at jk+1 w-level (divided by rho0_rcp) 
     718            zzatt = (   zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*r1_si0 )     & 
     719               &      + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*r1_si1 )   ) * wmask(ji,jj,jk+1) 
     720#if defined key_RK3 
     721            !                                            ! RK3 : temperature trend at jk t-level 
     722            pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + qsr(ji,jj) * ( zatt(ji,jj) - zzatt ) / ze3t 
     723#else 
     724            !                                            ! MLF : heat content trend due to Qsr flux (qsr_hc) 
     725            qsr_hc(ji,jj,jk) =  qsr(ji,jj) * ( zatt(ji,jj) - zzatt ) 
     726#endif 
     727            zatt(ji,jj) = zzatt                          ! save for the next level computation 
     728         END_2D 
     729!!stbug         !                                         !* sea-ice *!   store the 1st level attenuation coeff. 
     730!!stbug         IF( jk == 1 )   fraqsr_1lev(A2D(0)) = 1._wp - zatt(A2D(0)) * rho0_rcp 
     731      END DO 
     732!!st      IF(lwp) WRITE(numout,*) 'nk0+1= ', nk0+1, ' qsr max = ' , MAXVAL(zatt*qsr)*rho0_rcp, ' W/m2' , MAXVAL(zatt*qsr/e3t(:,:,nk0+1,Kmm)), ' K/s'  
     733      ! 
     734      DO jk = nk0+1, nkV                           !* deeper layers *!   (visible light only) 
     735         DO_2D( 0, 0, 0, 0 ) 
     736            ze3t  = e3t(ji,jj,jk,Kmm)                    ! light attenuation at jk+1 w-level (divided by rho0_rcp) 
     737            zzatt = (   zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*r1_si1 )   ) * wmask(ji,jj,jk+1) 
     738#if defined key_RK3 
     739            !                                            ! RK3 : temperature trend at jk t-level 
     740            pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + qsr(ji,jj) * ( zatt(ji,jj) - zzatt ) / ze3t 
     741#else 
     742            !                                            ! MLF : heat content trend due to Qsr flux (qsr_hc) 
     743            qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zatt(ji,jj) - zzatt ) 
     744#endif 
     745            zatt(ji,jj) = zzatt                       ! save for the next level computation 
     746         END_2D 
     747      END DO       
     748      ! 
     749!!st      IF(lwp) WRITE(numout,*) 'nkV+1= ', nkV+1, ' qsr max = ' , MAXVAL(zatt*qsr)*rho0_rcp, ' W/m2' , MAXVAL(zatt*qsr/e3t(:,:,nkV+1,Kmm)), ' K/s'  
     750   END SUBROUTINE qsr_2bd 
     751 
     752 
     753   FUNCTION qsr_ext_lev( pL, pfr ) RESULT( klev ) 
     754      !!---------------------------------------------------------------------- 
     755      !!                 ***  ROUTINE trc_oce_ext_lev  *** 
     756      !!        
     757      !! ** Purpose :   compute the maximum level of light penetration 
     758      !!           
     759      !! ** Method  :   the function provides the level at which irradiance, I,  
     760      !!              has a negligible effect on temperature. 
     761      !!                T(n+1)-T(n) = ∆t dk[I] / ( rho0 Cp e3t_k )   
     762      !!              I(k) has a negligible effect on temperature at level k if: 
     763      !!                ∆t I(k) / ( rho0*Cp*e3t_k ) <= 1.e-15 °C 
     764      !!              with I(z) = Qsr*pfr*EXP(-z/L), therefore : 
     765      !!                z >= L * LOG( 1.e-15 * rho0*Cp*e3t_k / ( ∆t*Qsr*pfr ) ) 
     766      !!              with Qsr being the maximum normal surface irradiance at sea  
     767      !!              level (~1000 W/m2). 
     768      !!                # pL is the longest depth of extinction: 
     769      !!                   - pL = 23.00 m (2 bands case) 
     770      !!                   - pL = 48.24 m (3 bands case: blue waveband & 0.03 mg/m2 for the chlorophyll) 
     771      !!                # pfr is the fraction of solar radiation which penetrates, 
     772      !!                considering Qsr=1000 W/m2 and rn_abs = 0.58: 
     773      !!                   - Qsr*pfr0 = Qsr *    rn_abs    = 580 W/m2   (top absorbtion) 
     774      !!                   - Qsr*pfr1 = Qsr * (1-rn_abs)   = 420 W/m2 (2 bands case) 
     775      !!                   - Qsr*pfr1 = Qsr * (1-rn_abs)/3 = 140 W/m2 (3 bands case & equi-partition) 
     776      !! 
     777      !!---------------------------------------------------------------------- 
     778      INTEGER              ::   klev   ! result: maximum level of light penetration 
     779      REAL(wp), INTENT(in) ::   pL     ! depth of extinction 
     780      REAL(wp), INTENT(in) ::   pfr    ! frac. solar radiation which penetrates  
     781      ! 
     782      INTEGER  ::   jk                 ! dummy loop index 
     783      REAL(wp) ::   zcoef              ! local scalar 
     784      REAL(wp) ::   zhext              ! deepest depth until which light penetrates 
     785      REAL(wp) ::   ze3t , zdw         ! max( e3t_k ) and min( w-depth_k+1 ) 
     786      REAL(wp) ::   zprec = 10.e-15_wp ! required precision  
     787      REAL(wp) ::   zQmax= 1000._wp    ! maximum normal surface irradiance at sea level (W/m2) 
     788      !!---------------------------------------------------------------------- 
     789      ! 
     790      zQmax = 1000._wp     ! maximum normal surface irradiance at sea level (W/m2) 
     791      ! 
     792      zcoef    =  zprec * rho0_rcp / ( rDt * zQmax * pfr) 
     793      ! 
     794      klev = jpkm1      ! Level of light extinction 
     795      DO jk = jpkm1, 1, -1 
     796         IF( SUM( tmask(:,:,jk) ) > 0 ) THEN    ! ocean point at that level 
     797            zdw  = MAXVAL( gdepw_0(:,:,jk+1) * wmask(:,:,jk) )    ! max w-depth at jk+1 level 
     798            ze3t = MINVAL(   e3t_0(:,:,jk  )                 )    ! minimum e3t at jk   level 
     799            zhext =  - pL * LOG( zcoef * ze3t )                   ! extinction depth 
     800            IF( zdw >= zhext )   klev = jk                        ! last T-level reached by Qsr 
     801         ELSE                                   ! only land point at level jk 
     802            klev = jk                                             ! local domain sea-bed level  
     803         ENDIF 
     804      END DO 
     805      ! 
     806!!st      IF(lwp) WRITE(numout,*) '                level of e3t light extinction = ', klev, ' ref depth = ', gdepw_1d(klev+1), ' m' 
     807!!st      ! 
     808!!st      klev = jpkm1      ! Level of light extinction 
     809!!st      DO jk = jpkm1, 1, -1 
     810!!st         IF( SUM( tmask(:,:,jk) ) > 0 ) THEN    ! ocean point at that level 
     811!!st            zdw  = MAXVAL( gdepw_0(:,:,jk+1) * wmask(:,:,jk) )    ! max w-depth at jk+1 level 
     812!!st            ze3t = MINVAL(   e3t_0(:,:,jk  )                 )    ! minimum e3t at jk   level 
     813!!st            zhext =  - pL * LOG( zcoef * pL )                   ! extinction depth 
     814!!st            IF( zdw >= zhext )   klev = jk                        ! last T-level reached by Qsr 
     815!!st         ELSE                                   ! only land point at level jk 
     816!!st            klev = jk                                             ! local domain sea-bed level  
     817!!st         ENDIF 
     818!!st      END DO 
     819!!st      ! 
     820!!st      IF(lwp) WRITE(numout,*) '                level of pL light extinction = ', klev, ' ref depth = ', gdepw_1d(klev+1), ' m' 
     821      ! 
     822   END FUNCTION qsr_ext_lev 
    337823 
    338824 
     
    355841      !!---------------------------------------------------------------------- 
    356842      INTEGER  ::   ji, jj, jk                  ! dummy loop indices 
    357       INTEGER  ::   ios, irgb, ierror, ioptio   ! local integer 
    358       REAL(wp) ::   zz0, zc0 , zc1, zcoef      ! local scalars 
    359       REAL(wp) ::   zz1, zc2 , zc3, zchl       !   -      - 
     843      INTEGER  ::   ios, ierror, ioptio   ! local integer 
     844      REAL(wp) ::   zLB, zLG, zLR         ! local scalar 
     845      REAL(wp) ::   zVlp, zchl            !   -      - 
    360846      ! 
    361847      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files 
     
    368854      READ  ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901) 
    369855901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_qsr in reference namelist' ) 
    370       ! 
    371       READ  ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 ) 
     856      READ  ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902) 
    372857902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist' ) 
    373858      IF(lwm) WRITE ( numond, namtra_qsr ) 
    374859      ! 
    375       IF(lwp) THEN                ! control print 
     860      IF(lwp) THEN            !**  control print  **! 
    376861         WRITE(numout,*) 
    377862         WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation' 
     
    383868         WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)        nn_chldta  = ', nn_chldta 
    384869         WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs     = ', rn_abs 
    385          WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0     = ', rn_si0 
    386          WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1     = ', rn_si1 
     870         WRITE(numout,*) '      RGB & 2 bands: shortess attenuation depth    rn_si0     = ', rn_si0 
     871         WRITE(numout,*) '      2 bands: longest attenuation depth           rn_si1     = ', rn_si1 
    387872         WRITE(numout,*) 
    388873      ENDIF 
    389874      ! 
    390       ioptio = 0                    ! Parameter control 
     875      ioptio = 0              !**  Parameter control  **! 
    391876      IF( ln_qsr_rgb  )   ioptio = ioptio + 1 
    392877      IF( ln_qsr_2bd  )   ioptio = ioptio + 1 
     
    401886      IF( ln_qsr_bio                      )   nqsr = np_BIO 
    402887      ! 
    403       !                             ! Initialisation 
    404       xsi0r = 1._wp / rn_si0 
    405       xsi1r = 1._wp / rn_si1 
     888      !                       !**  Initialisation  **! 
     889      ! 
     890      !                                !==  Infrared attenuation  ==!   (all schemes) 
     891      !                                !============================! 
     892      ! 
     893      r1_si0 = 1._wp / rn_si0                ! inverse of infrared attenuation length 
     894      ! 
     895      nk0 = qsr_ext_lev( rn_si0, rn_abs )    ! level of light extinction 
     896      ! 
     897      IF(lwp) WRITE(numout,*) '   ==>>>   Infrared light attenuation' 
     898      IF(lwp) WRITE(numout,*) '              level of infrared extinction = ', nk0, ' ref depth = ', gdepw_1d(nk0+1), ' m' 
     899      IF(lwp) WRITE(numout,*) 
    406900      ! 
    407901      SELECT CASE( nqsr ) 
    408902      ! 
    409       CASE( np_RGB , np_RGBc )         !==  Red-Green-Blue light penetration  ==! 
    410          ! 
    411          IF(lwp)   WRITE(numout,*) '   ==>>>   R-G-B   light penetration ' 
     903      CASE( np_RGBc, np_RGB )          !==  Red-Green-Blue light attenuation  ==!   (Chl data or constant) 
     904         !                             !========================================! 
     905         ! 
     906         IF( nqsr == np_RGB ) THEN   ;   zchl   = 0.05        ! constant Chl value 
     907         ELSE                        ;   zchl   = 0.03        ! minimum  Chl value 
     908         ENDIF 
     909         zchl   = MAX( 0.03_wp , MIN( zchl , 10._wp) )     ! NB. make sure that chosen value verifies: 0.03 < zchl < 10 
     910         nc_rgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )    ! Convert Chl value to attenuation coefficient look-up table index 
    412911         ! 
    413912         CALL trc_oce_rgb( rkrgb )                 ! tabulated attenuation coef. 
    414913         ! 
    415          nksr = trc_oce_ext_lev( r_si2, 33._wp )   ! level of light extinction 
    416          ! 
    417          IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m' 
     914         zVlp =  ( 1._wp - rn_abs ) / 3._wp        ! visible light equi-partition 
     915         ! 
     916         !     1 / length          !   attenuation  length   !         attenuation level 
     917         r1_LR = rkrgb(3,nc_rgb)   ;   zLR = 1._wp / r1_LR   ;   nkR = qsr_ext_lev( zLR, zVlp )   ! Red    
     918         r1_LG = rkrgb(2,nc_rgb)   ;   zLG = 1._wp / r1_LG   ;   nkG = qsr_ext_lev( zLG, zVlp )   ! Green 
     919         r1_LB = rkrgb(1,nc_rgb)   ;   zLB = 1._wp / r1_LB   ;   nkB = qsr_ext_lev( zLB, zVlp )   ! Blue 
     920         ! 
     921         nkV = nkB                                 ! maximum level of light penetration 
     922         ! 
     923         IF( nqsr == np_RGB ) THEN 
     924            IF(lwp) WRITE(numout,*) '   ==>>>   RGB:  light attenuation with a constant Chlorophyll = ', zchl 
     925         ELSE 
     926            IF(lwp) WRITE(numout,*) '   ==>>>   RGB:  light attenuation using Chlorophyll data with min(Chl) = ', zchl 
     927         ENDIF             
     928         IF(lwp) WRITE(numout,*) '                 level of Red   extinction = ', nkR, ' ref depth = ', gdepw_1d(nkR+1), ' m' 
     929         IF(lwp) WRITE(numout,*) '                 level of Green extinction = ', nkG, ' ref depth = ', gdepw_1d(nkG+1), ' m' 
     930         IF(lwp) WRITE(numout,*) '                 level of Blue  extinction = ', nkB, ' ref depth = ', gdepw_1d(nkB+1), ' m' 
     931         IF(lwp) WRITE(numout,*) 
    418932         ! 
    419933         IF( nqsr == np_RGBc ) THEN                ! Chl data : set sf_chl structure 
     
    423937               CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' )   ;   RETURN 
    424938            ENDIF 
    425             ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1)   ) 
     939            ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1) ) 
    426940            IF( sn_chl%ln_tint )   ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) ) 
    427941            !                                        ! fill sf_chl with sn_chl and control print 
    428             CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',   & 
     942            CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',                             & 
    429943               &           'Solar penetration function of read chlorophyll', 'namtra_qsr' , no_print ) 
    430944         ENDIF 
    431          IF( nqsr == np_RGB ) THEN                 ! constant Chl 
    432             IF(lwp) WRITE(numout,*) '   ==>>>   Constant Chlorophyll concentration = 0.05' 
    433          ENDIF 
    434          ! 
    435       CASE( np_2BD )                   !==  2 bands light penetration  ==! 
    436          ! 
    437          IF(lwp)  WRITE(numout,*) '   ==>>>   2 bands light penetration' 
    438          ! 
    439          nksr = trc_oce_ext_lev( rn_si1, 100._wp )    ! level of light extinction 
    440          IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m' 
     945         ! 
     946      CASE( np_2BD )                   !==  2 bands light attenuation (IR+ visible light) ==! 
     947         ! 
     948         ! 
     949         r1_si1 = 1._wp / rn_si1                     ! inverse of visible light attenuation 
     950         zVlp =  ( 1._wp - rn_abs )                  ! visible light partition 
     951         nkV  = qsr_ext_lev( rn_si1, zVlp )          ! level of visible light extinction 
     952         ! 
     953         IF(lwp) WRITE(numout,*) '   ==>>>   2 bands attenuation (Infrared + Visible light) ' 
     954         IF(lwp) WRITE(numout,*) '                level of visible light extinction = ', nkV, ' ref depth = ', gdepw_1d(nkV+1), ' m' 
     955         IF(lwp) WRITE(numout,*) 
    441956         ! 
    442957      CASE( np_BIO )                   !==  BIO light penetration  ==! 
     
    447962         CALL trc_oce_rgb( rkrgb )                 ! tabulated attenuation coef. 
    448963         ! 
    449          nksr = trc_oce_ext_lev( r_si2, 33._wp )   ! level of light extinction 
    450          ! 
    451          IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m' 
     964         nkV = trc_oce_ext_lev( r_si2, 33._wp )    ! maximum level of light extinction 
     965         ! 
     966         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nkV, ' ref depth = ', gdepw_1d(nkV+1), ' m' 
    452967         ! 
    453968      END SELECT 
    454969      ! 
    455       qsr_hc(:,:,:) = 0._wp     ! now qsr heat content set to zero where it will not be computed 
    456       ! 
    457       ! 1st ocean level attenuation coefficient (used in sbcssm) 
     970      nksr = nkV       ! name of max level of light extinction used in traatf(_qco).F90 
     971      ! 
     972#if ! defined key_RK3 
     973      qsr_hc(:,:,:) = 0._wp      ! MLF : now qsr heat content set to zero where it will not be computed 
     974#endif 
     975      ! 
     976      !                          ! Sea-ice :   1st ocean level attenuation coefficient (used in sbcssm) 
    458977      IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN 
    459978         CALL iom_get( numror, jpdom_auto, 'fraqsr_1lev'  , fraqsr_1lev  ) 
Note: See TracChangeset for help on using the changeset viewer.