New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
traqsr.F90 in NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/TRA – NEMO

source: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/TRA/traqsr.F90 @ 15574

Last change on this file since 15574 was 15574, checked in by techene, 3 years ago

#2605 #2715 trunk merged into dev_r14318_RK3_stage1
  • Property svn:keywords set to Id
File size: 59.5 KB
Line 
1MODULE traqsr
2   !!======================================================================
3   !!                       ***  MODULE  traqsr  ***
4   !! Ocean physics:   solar radiation penetration in the top ocean levels
5   !!======================================================================
6   !! History :  OPA  !  1990-10  (B. Blanke)  Original code
7   !!            7.0  !  1991-11  (G. Madec)
8   !!                 !  1996-01  (G. Madec)  s-coordinates
9   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
10   !!             -   !  2005-11  (G. Madec) zco, zps, sco coordinate
11   !!            3.2  !  2009-04  (G. Madec & NEMO team)
12   !!            3.6  !  2012-05  (C. Rousset) store attenuation coef for use in ice model
13   !!            3.6  !  2015-12  (O. Aumont, J. Jouanno, C. Ethe) use vertical profile of chlorophyll
14   !!            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
17   !!----------------------------------------------------------------------
18
19   !!----------------------------------------------------------------------
20   !!   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
25   !!   tra_qsr_init  : initialization of the qsr penetration
26   !!----------------------------------------------------------------------
27   USE oce            ! ocean dynamics and active tracers
28   USE phycst         ! physical constants
29   USE dom_oce        ! ocean space and time domain
30   USE domtile
31   USE sbc_oce        ! surface boundary condition: ocean
32   USE trc_oce        ! share SMS/Ocean variables
33   USE trd_oce        ! trends: ocean variables
34   USE trdtra         ! trends manager: tracers
35   !
36   USE in_out_manager ! I/O manager
37   USE prtctl         ! Print control
38   USE iom            ! I/O library
39   USE fldread        ! read input fields
40   USE restart        ! ocean restart
41   USE lib_mpp        ! MPP library
42   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
43   USE timing         ! Timing
44
45   IMPLICIT NONE
46   PRIVATE
47
48   PUBLIC   tra_qsr       ! routine called by step.F90 (ln_traqsr=T)
49   PUBLIC   tra_qsr_init  ! routine called by nemogcm.F90
50
51   !                                 !!* Namelist namtra_qsr: penetrative solar radiation
52   LOGICAL , PUBLIC ::   ln_traqsr    !: light absorption (qsr) flag
53   LOGICAL , PUBLIC ::   ln_qsr_rgb   !: Red-Green-Blue light absorption flag
54   LOGICAL , PUBLIC ::   ln_qsr_2bd   !: 2 band         light absorption flag
55   LOGICAL , PUBLIC ::   ln_qsr_bio   !: bio-model      light absorption flag
56   INTEGER , PUBLIC ::   nn_chldta    !: use Chlorophyll data (=1) or not (=0)
57   REAL(wp), PUBLIC ::   rn_abs       !: fraction absorbed in the very near surface (RGB & 2 bands)
58   REAL(wp), PUBLIC ::   rn_si0       !: very near surface depth of extinction      (RGB & 2 bands)
59   REAL(wp), PUBLIC ::   rn_si1       !: deepest depth of extinction (water type I)       (2 bands)
60   !
61   INTEGER, PARAMETER ::   np_RGB  = 1   ! R-G-B     light penetration with constant Chlorophyll
62   INTEGER, PARAMETER ::   np_RGBc = 2   ! R-G-B     light penetration with Chlorophyll data
63   INTEGER, PARAMETER ::   np_2BD  = 3   ! 2 bands   light penetration
64   INTEGER, PARAMETER ::   np_BIO  = 4   ! bio-model light penetration
65   !
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)
80   !
81   REAL(wp) , PUBLIC, DIMENSION(3,61)   ::   rkrgb    ! tabulated attenuation coefficients for RGB absorption
82   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read)
83
84   !! * Substitutions
85#  include "do_loop_substitute.h90"
86#  include "domzgr_substitute.h90"
87   !!----------------------------------------------------------------------
88   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
89   !! $Id$
90   !! Software governed by the CeCILL license (see ./LICENSE)
91   !!----------------------------------------------------------------------
92CONTAINS
93
94   SUBROUTINE tra_qsr( kt, Kmm, pts, Krhs )
95      !!----------------------------------------------------------------------
96      !!                  ***  ROUTINE tra_qsr  ***
97      !!
98      !! ** Purpose :   Compute the temperature trend due to the solar radiation
99      !!              penetration and add it to the general temperature trend.
100      !!
101      !! ** Method  : The profile of the solar radiation within the ocean is defined
102      !!      through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) or computed by
103      !!      the biogeochemical model
104      !!         The computation is only done down to the level where
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
108      !!              - send  trend for further diagnostics (l_trdtra=T)
109      !!----------------------------------------------------------------------
110      INTEGER,                                   INTENT(in   ) ::   kt, Kmm, Krhs   ! ocean time-step and time-level indices
111      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) ::   pts             ! active tracers and RHS of tracer equation
112      !
113      INTEGER  ::   ji, jj, jk               ! dummy loop indices
114      REAL(wp) ::   z1_2, ze3t                     ! local scalars
115      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, zetot
116      !!----------------------------------------------------------------------
117      !
118      IF( ln_timing )   CALL timing_start('tra_qsr')
119      !
120      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
121         IF( kt == nit000 ) THEN
122            IF(lwp) WRITE(numout,*)
123            IF(lwp) WRITE(numout,*) 'tra_qsr : penetration of the surface solar radiation'
124            IF(lwp) WRITE(numout,*) '~~~~~~~'
125         ENDIF
126      ENDIF
127      !
128      IF( l_trdtra ) THEN     ! trends diagnostic: save the input temperature trend
129         ALLOCATE( ztrdt(jpi,jpj,jpk) )
130         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
131      ENDIF
132      !
133#if ! defined key_RK3
134      !                        ! MLF only : heat content trend due to Qsr flux (qsr_hc)
135      !
136      !                         !-----------------------------------!
137      !                         !  before qsr induced heat content  !
138      !                         !-----------------------------------!
139      IF( kt == nit000 ) THEN          !==  1st time step  ==!
140         IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN    ! read in restart
141            z1_2 = 0.5_wp
142            IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                        ! Do only on the first tile
143               IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field read in the restart file'
144               CALL iom_get( numror, jpdom_auto, 'qsr_hc_b', qsr_hc_b )   ! before heat content trend due to Qsr flux
145            ENDIF
146         ELSE                                           ! No restart or Euler forward at 1st time step
147            z1_2 = 1._wp
148            DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
149               qsr_hc_b(ji,jj,jk) = 0._wp
150            END_3D
151         ENDIF
152      ELSE                             !==  Swap of qsr heat content  ==!
153         z1_2 = 0.5_wp
154         DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
155            qsr_hc_b(ji,jj,jk) = qsr_hc(ji,jj,jk)
156         END_3D
157      ENDIF
158#endif
159     
160      !                       !----------------------------!
161      SELECT CASE( nqsr )     !  qsr induced heat content  !
162      !                       !----------------------------!
163      !
164      CASE( np_RGBc )   ;   CALL qsr_RGBc( kt, Kmm, pts, Krhs )  !==  R-G-B fluxes using chlorophyll data     ==!    with Morel &Berthon (1989) vertical profile
165         !
166      CASE( np_RGB  )   ;   CALL qsr_RGB ( kt, Kmm, pts, Krhs )  !==  R-G-B fluxes with constant chlorophyll  ==!   
167         !
168      CASE( np_2BD  )   ;   CALL qsr_2BD (     Kmm, pts, Krhs )  !==  2-bands fluxes                          ==!
169         !
170      CASE( np_BIO )                                     !==  bio-model fluxes                        ==!
171         DO_3D( 0, 0, 0, 0, 1, nkV )
172#if defined key_RK3
173            !                                                  !- RK3 : temperature trend at jk t-level
174            ze3t   = e3t(ji,jj,jk,Kmm)
175            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
176#else
177            !                                                  !- MLF : heat content trend due to Qsr flux (qsr_hc)
178            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) )
179#endif
180         END_3D
181         !                                                     !- sea-ice : store the 1st level attenuation coefficient
182         WHERE( etot3(A2D(0),1) /= 0._wp )   ;   fraqsr_1lev(A2D(0)) = 1._wp - etot3(A2D(0),2) / etot3(A2D(0),1)
183         ELSEWHERE                           ;   fraqsr_1lev(A2D(0)) = 1._wp
184         END WHERE
185         !
186      END SELECT
187      !
188#if defined key_RK3
189      !                             ! RK3 : diagnostics/output
190      IF( l_trdtra .OR. iom_use('qsr3d') ) THEN     ! qsr diagnostics
191         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
192         !                                         ! qsr tracers trends saved for diagnostics
193         IF( l_trdtra )   CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt )
194         IF( iom_use('qsr3d') ) THEN               ! qsr distribution
195            DO jk = nkV, 1, -1
196               ztrdt(:,:,jk) = ztrdt(:,:,jk+1) + qsr_hc(:,:,jk) * rho0_rcp
197            END DO
198            CALL iom_put( 'qsr3d', ztrdt )   ! 3D distribution of shortwave Radiation
199         ENDIF
200         DEALLOCATE( ztrdt )
201      ENDIF
202#else
203      !                             ! MLF : add the temperature trend
204      DO_3D( 0, 0, 0, 0, 1, nksr )
205         pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
206            &                      + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) )   &
207            &                             / e3t(ji,jj,jk,Kmm)
208      END_3D
209      !
210!!st7-2
211      ! sea-ice: store the 1st ocean level attenuation coefficient
212      DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
213         IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) )
214         ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp
215         ENDIF
216      END_2D
217      !                                       !===>>> CAUTION: lbc_lnk is required on fraqsr_lev since sea ice computes on the full domain
218      !                                       !                otherwise restartability and reproducibility are broken
219      CALL lbc_lnk( 'tra_qsr', fraqsr_1lev(:,:), 'T', 1._wp )
220!!st      CALL lbc_lnk( 'tra_qsr', qsr_hc(:,:,:), 'T', 1._wp )
221      !
222      IF( iom_use('qsr3d') ) THEN      ! output the shortwave Radiation distribution
223         ALLOCATE( zetot(A2D(nn_hls),jpk) )
224         zetot(:,:,nksr+1:jpk) = 0._wp     ! below ~400m set to zero
225         DO_3DS(0, 0, 0, 0, nksr, 1, -1)
226            zetot(ji,jj,jk) = zetot(ji,jj,jk+1) + qsr_hc(ji,jj,jk) * rho0_rcp
227         END_3D
228         CALL iom_put( 'qsr3d', zetot )   ! 3D distribution of shortwave Radiation
229         DEALLOCATE( zetot )
230      ENDIF
231      !
232      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
233         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
234         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt )
235         DEALLOCATE( ztrdt )
236      ENDIF
237#endif
238      !
239      IF( .NOT. l_istiled .OR. ntile == nijtile )  THEN                ! Do only on the last tile
240         IF( lrst_oce ) THEN     ! write in the ocean restart file
241            CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc      )
242            CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev )
243         ENDIF
244      ENDIF
245      !
246      !                       ! print mean trends (used for debugging)
247      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
248      !
249      IF( ln_timing )   CALL timing_stop('tra_qsr')
250      !
251   END SUBROUTINE tra_qsr
252
253
254   SUBROUTINE qsr_RGBc( kt, Kmm, pts, Krhs )
255      !!----------------------------------------------------------------------
256      !!                  ***  ROUTINE qsr_RGBc  ***
257      !!
258      !! ** Purpose :   Red-Green-Blue solar radiation using chlorophyll data
259      !!
260      !! ** Method  : The profile of the solar radiation within the ocean is defined
261      !!      through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) and a ratio rn_abs
262      !!      Considering the 2 wavebands case:
263      !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
264      !!         The temperature trend associated with the solar radiation penetration
265      !!         is given by : zta = 1/e3t dk[ I ] / (rho0*Cp)
266      !!         At the bottom, boudary condition for the radiation is no flux :
267      !!      all heat which has not been absorbed in the above levels is put
268      !!      in the last ocean level.
269      !!         The computation is only done down to the level where
270      !!      I(k) < 1.e-15 W/m2 (i.e. over the top nk levels) .
271      !!
272      !! ** Action  : - update ta with the penetrative solar radiation trend
273      !!              - send  trend for further diagnostics (l_trdtra=T)
274      !!
275      !! Reference  : Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
276      !!              Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562
277      !!----------------------------------------------------------------------
278      INTEGER,                                   INTENT(in   ) ::   kt, Kmm, Krhs   ! ocean time-step and time-level indices
279      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) ::   pts         ! active tracers and RHS of tracer equation
280      !!
281      INTEGER  ::   ji, jj, jk               ! dummy loop indices
282      INTEGER  ::   irgb                     ! local integer
283      REAL(wp) ::   zc1 , zc2 , zc3, zchl    ! local scalars
284      REAL(wp) ::   zze0, zzeR, zzeG, zzeB, zzeT              !    -         -
285      REAL(wp) ::   zz0 , zz1 , ze3t                          !    -         -
286      REAL(wp) ::   zCb, zCmax, zpsi, zpsimax, zrdpsi, zCze   !    -         -
287      REAL(wp) ::   zlogc, zlogze, zlogCtot, zlogCze          !    -         -
288      REAL(wp), DIMENSION(A2D(0)    ) ::   ze0, zeR, zeG, zeB, zeT
289      REAL(wp), DIMENSION(A2D(0),0:3) ::   zc
290      !!----------------------------------------------------------------------
291      !
292      !
293      !                       !===========================================!
294      !                       !==  R-G-B fluxes using chlorophyll data  ==!    with Morel &Berthon (1989) vertical profile
295      !                       !===================================****====!
296      !
297      !                             !=  Chlorophyll data  =!
298      !
299      IF( ntile == 0 .OR. ntile == 1 )  THEN       ! Do only for the full domain
300         IF( ln_tile )   CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )   ! Use full domain
301         CALL fld_read( kt, 1, sf_chl )                                       ! Read Chl data and provides it at the current time step
302         IF( ln_tile )   CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 1 )   ! Revert to tile domain
303      ENDIF
304      !
305       DO_2D( 0, 0, 0, 0 )                          ! pre-calculated expensive coefficient
306         zlogc = LOG(  MAX( 0.03_wp, MIN( sf_chl(1)%fnow(ji,jj,1) ,10._wp ) )  ) ! zlogc = log(zchl)   with 0.03 <= Chl >= 10.
307         zc1   = 0.113328685307 + 0.803 * zlogc                               ! zc1 : log(zCze)  = log (1.12  * zchl**0.803)
308         zc2   = 3.703768066608 + 0.459 * zlogc                               ! zc2 : log(zCtot) = log(40.6  * zchl**0.459)
309         zc3   = 6.34247346942  - 0.746 * zc2                                 ! zc3 : log(zze)   = log(568.2 * zCtot**(-0.746))
310         IF( zc3 > 4.62497281328 )   zc3 = 5.298317366548 - 0.293 * zc2       ! IF(log(zze)>log(102)) log(zze) = log(200*zCtot**(-0.293))
311         !
312         zc(ji,jj,0) = zlogc                                                  ! ze(0) = log(zchl)
313         zc(ji,jj,1) = EXP( zc1 )                                             ! ze(1) = zCze
314         zc(ji,jj,2) = 1._wp / ( 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) )  ! ze(2) = 1/zdelpsi
315         zc(ji,jj,3) = EXP( - zc3 )                                           ! ze(3) = 1/zze
316      END_2D
317      !
318      !                             !=  surface light  =!
319      !
320      zz0 =           rn_abs              ! Infrared absorption
321      zz1 = ( 1._wp - rn_abs ) / 3._wp    ! R-G-B equi-partition
322      !
323      DO_2D( 0, 0, 0, 0 )                 ! surface light
324         ze0(ji,jj) = zz0 * qsr(ji,jj)   ;   zeR(ji,jj) = zz1 * qsr(ji,jj)    ! IR    ; Red
325         zeG(ji,jj) = zz1 * qsr(ji,jj)   ;   zeB(ji,jj) = zz1 * qsr(ji,jj)    ! Green ; Blue
326         zeT(ji,jj) =       qsr(ji,jj)                                        ! Total
327      END_2D
328      !             
329      !                             !=  interior light  =!
330      !
331      DO jk = 1, nk0                      !* near surface layers *!   (< ~12 meters : IR + RGB )
332         DO_2D( 0, 0, 0, 0 )
333            !                                      !- inverse of RGB attenuation lengths
334            zlogc     = zc(ji,jj,0)
335            zCb       = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) )
336            zCmax     = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 )
337            zpsimax   = 0.6   - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) )
338            ! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 )
339            zCze   = zc(ji,jj,1)
340            zrdpsi = zc(ji,jj,2)                                     ! 1/zdelpsi
341!!st05            zpsi   = zc(ji,jj,3) * gdepw(ji,jj,jk,Kmm)               ! gdepw/zze
342            zpsi   = zc(ji,jj,3) * gdepw(ji,jj,jk+1,Kmm)               ! gdepw/zze
343            !                                                        ! make sure zchl value is such that: 0.03 < zchl < 10.
344            zchl = MAX(  0.03_wp , MIN( zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) , 10._wp )  )
345            !                                                        ! Convert chlorophyll value to attenuation coefficient
346            irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )             ! look-up table index
347            !       Red             !         Green              !         Blue
348            r1_LR = rkrgb(3,irgb)   ;   r1_LG = rkrgb(2,irgb)    ;   r1_LB = rkrgb(1,irgb)
349            !
350            !                                      !- fluxes at jk+1 w-level
351            ze3t = e3t(ji,jj,jk,Kmm)
352            zze0 = ze0(ji,jj) * EXP( - ze3t*r1_si0 )   ;   zzeR = zeR(ji,jj) * EXP( - ze3t*r1_LR )   ! IR    ; Red  at jk+1 w-level
353            zzeG = zeG(ji,jj) * EXP( - ze3t*r1_LG  )   ;   zzeB = zeB(ji,jj) * EXP( - ze3t*r1_LB )   ! Green ; Blue      -      -
354            zzeT = ( zze0 + zzeB + zzeG + zzeR ) * wmask(ji,jj,jk+1)                                 ! Total             -      -
355!!st01            zzeT = ( zze0 + zzeR + zzeG + zzeB ) * wmask(ji,jj,jk+1)                                 ! Total             -      -
356            !
357#if defined key_RK3
358            !                                      !- RK3 : temperature trend at jk t-level
359            pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t
360#else
361            !                                      !- MLF : heat content trend due to Qsr flux (qsr_hc)
362            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT )
363#endif
364            ze0(ji,jj) = zze0   ;   zeR(ji,jj) = zzeR           ! IR    ; Red  store at jk+1 w-level
365            zeG(ji,jj) = zzeG   ;   zeB(ji,jj) = zzeB           ! Green ; Blue   -        -      -
366            zeT(ji,jj) = zzeT                                   ! total          -        -      -
367         END_2D
368         !
369      END DO
370      !
371      IF( kt == nit000 ) THEN
372         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'
373         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'
374         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'
375         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'
376         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'
377      ENDIF
378      !
379      DO jk = nk0+1, nkR                  !* down to Red extinction *!   (< ~71 meters : RGB , IR removed from calculation)
380          DO_2D( 0, 0, 0, 0 )
381            !                                      !- inverse of RGB attenuation lengths
382            zlogc     = zc(ji,jj,0)
383            zCb       = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) )
384            zCmax     = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 )
385            zpsimax   = 0.6   - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) )
386            ! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 )
387            zCze   = zc(ji,jj,1)
388            zrdpsi = zc(ji,jj,2)                               ! 1/zdelpsi
389            zpsi   = zc(ji,jj,3) * gdepw(ji,jj,jk+1,Kmm)         ! gdepw/zze
390!!st05            zpsi   = zc(ji,jj,3) * gdepw(ji,jj,jk,Kmm)         ! gdepw/zze
391            !                                                  ! make sure zchl value is such that: 0.03 < zchl < 10.
392            zchl = MAX(  0.03_wp , MIN( zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) , 10._wp )  )
393            !                                                  ! Convert chlorophyll value to attenuation coefficient
394            irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )       ! look-up table index
395            !       Red             !         Green              !         Blue
396            r1_LR = rkrgb(3,irgb)   ;   r1_LG = rkrgb(2,irgb)    ;   r1_LB = rkrgb(1,irgb)
397            !
398            !                                      !- fluxes at jk+1 w-level
399            ze3t = e3t(ji,jj,jk,Kmm)
400            zzeR = zeR(ji,jj) * EXP( - ze3t*r1_LR )                                                 ! Red          at jk+1 w-level
401            zzeG = zeG(ji,jj) * EXP( - ze3t*r1_LG )   ;   zzeB = zeB(ji,jj) * EXP( - ze3t*r1_LB )   ! Green ; Blue      -      -
402            zzeT = ( zzeR + zzeG + zzeB ) * wmask(ji,jj,jk+1)                                       ! Total             -      -
403            !
404#if defined key_RK3
405            !                                      !- RK3 : temperature trend at jk t-level
406            pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t
407#else
408            !                                      !- MLF : heat content trend due to Qsr flux (qsr_hc)
409            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT )
410#endif
411            zeR(ji,jj) = zzeR                                  ! Red          store at jk+1 w-level
412            zeG(ji,jj) = zzeG   ;   zeB(ji,jj) = zzeB          ! Green ; Blue   -        -      -
413            zeT(ji,jj) = zzeT                                  ! total          -        -      -
414         END_2D
415      END DO
416      !
417      IF( kt == nit000 ) THEN
418         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'
419         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'
420         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'
421         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'
422      ENDIF
423      !
424      DO jk = nkR+1, nkG                  !* down to Green extinction *!   (< ~350 m : GB , IR+R removed from calculation)
425          DO_2D( 0, 0, 0, 0 )
426            !                                      !- inverse of RGB attenuation lengths
427            zlogc     = zc(ji,jj,0)
428            zCb       = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) )
429            zCmax     = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 )
430            zpsimax   = 0.6   - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) )
431            ! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 )
432            zCze   = zc(ji,jj,1)
433            zrdpsi = zc(ji,jj,2)                               ! 1/zdelpsi
434            zpsi   = zc(ji,jj,3) * gdepw(ji,jj,jk+1,Kmm)         ! gdepw/zze
435!!st05            zpsi   = zc(ji,jj,3) * gdepw(ji,jj,jk,Kmm)         ! gdepw/zze
436            !                                                  ! make sure zchl value is such that: 0.03 < zchl < 10.
437            zchl = MAX(  0.03_wp , MIN( zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) , 10._wp )  )
438            !                                                  ! Convert chlorophyll value to attenuation coefficient
439            irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )       ! look-up table index
440            !     Green              !         Blue
441            r1_LG = rkrgb(2,irgb)    ;   r1_LB = rkrgb(1,irgb)
442            !
443            !                                      !- fluxes at jk+1 w-level
444            ze3t = e3t(ji,jj,jk,Kmm)
445            zzeG = zeG(ji,jj) * EXP( - ze3t * r1_LG )   ;   zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB ) ! Green ; Blue
446            zzeT = ( zzeG + zzeB ) * wmask(ji,jj,jk+1)                                                ! Total             -      -
447#if defined key_RK3
448            !                                      !- RK3 : temperature trend at jk t-level
449            pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t
450#else
451            !                                      !- MLF : heat content trend due to Qsr flux (qsr_hc)
452            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT )
453#endif
454            zeG(ji,jj) = zzeG   ;   zeB(ji,jj) = zzeB          ! Green ; Blue store at jk+1 w-level
455            zeT(ji,jj) = zzeT                                  ! total          -        -      -
456         END_2D
457      END DO
458      !
459      IF( kt == nit000 ) THEN
460         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'
461         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'
462         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'
463      ENDIF
464      !
465      DO jk = nkG+1, nkB                  !* down to Blue extinction *!   (< ~1300 m : B , IR+RG removed from calculation)
466          DO_2D( 0, 0, 0, 0 )
467            !                                      !- inverse of RGB attenuation lengths
468            zlogc     = zc(ji,jj,0)
469            zCb       = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) )
470            zCmax     = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 )
471            zpsimax   = 0.6   - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) )
472            ! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 )
473            zCze   = zc(ji,jj,1)
474            zrdpsi = zc(ji,jj,2)                               ! 1/zdelpsi
475            zpsi   = zc(ji,jj,3) * gdepw(ji,jj,jk+1,Kmm)         ! gdepw/zze
476!!st05            zpsi   = zc(ji,jj,3) * gdepw(ji,jj,jk,Kmm)         ! gdepw/zze
477            !                                                  ! make sure zchl value is such that: 0.03 < zchl < 10.
478            zchl = MAX(  0.03_wp , MIN( zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) , 10._wp )  )
479            !                                                  ! Convert chlorophyll value to attenuation coefficient
480            irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )       ! look-up table index
481            r1_LB = rkrgb(1,irgb)                              ! Blue
482            !
483            !                                      !- fluxes at jk+1 w-level
484            ze3t = e3t(ji,jj,jk,Kmm)
485            zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB )          ! Blue
486            zzeT = ( zzeB ) * wmask(ji,jj,jk+1)                ! Total             -      -
487#if defined key_RK3
488            !                                      !- RK3 : temperature trend at jk t-level
489            pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t
490#else
491            !                                      !- MLF : heat content trend due to Qsr flux (qsr_hc)
492            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT )
493#endif
494            zeB(ji,jj) = zzeB                                  ! Blue store at jk+1 w-level
495            zeT(ji,jj) = zzeT                                  ! total  -        -      -
496         END_2D
497      END DO
498      !
499      IF( kt == nit000 ) THEN
500         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'
501      ENDIF
502      !
503   END SUBROUTINE qsr_RGBc
504
505
506   SUBROUTINE qsr_RGB( kt, Kmm, pts, Krhs )
507      !!----------------------------------------------------------------------
508      !!                  ***  ROUTINE qsr_RGB  ***
509      !!
510      !! ** Purpose :   Red-Green-Blue solar radiation with constant chlorophyll
511      !!
512      !! ** Method  : The profile of the solar radiation within the ocean is defined
513      !!      through 2 wavebands (rn_si0,rn_si1) or 1 (rn_si0,rn_abs) + 3 wavebands (RGB)
514      !!         At the bottom, boudary condition for the radiation is no flux :
515      !!      all heat which has not been absorbed in the above levels is put
516      !!      in the last ocean level.
517      !!         For each band, the computation is only done down to the level where
518      !!      I(k) < 1.e-15 W/m2 (i.e. over the top nk levels) .
519      !!
520      !! ** Action  : - update ta with the penetrative solar radiation trend
521      !!              - send  trend for further diagnostics (l_trdtra=T)
522      !!
523      !! Reference  : Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
524      !!              Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562
525      !!----------------------------------------------------------------------
526      INTEGER,                                   INTENT(in   ) ::   kt, Kmm, Krhs   ! ocean time-step and time-level indices
527      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) ::   pts             ! active tracers and RHS of tracer equation
528      !!
529      INTEGER  ::   ji, jj, jk               ! dummy loop indices
530      REAL(wp) ::   zze0, zzeR, zzeG, zzeB, zzeT              !    -         -
531      REAL(wp) ::   zz0 , zz1 , ze3t                          !    -         -
532      REAL(wp), DIMENSION(A2D(0))   ::   ze0, zeR, zeG, zeB, zeT
533      !!----------------------------------------------------------------------
534      !     
535      !
536      !                       !==============================================!
537      !                       !==  R-G-B fluxes with constant chlorophyll  ==!   
538      !                       !======================********================!
539      !
540      !                             !=  surface light  =!
541      !
542      zz0 =           rn_abs              ! Infrared absorption
543      zz1 = ( 1._wp - rn_abs ) / 3._wp    ! surface equi-partition in R-G-B
544      !
545      DO_2D( 0, 0, 0, 0 )                 ! surface light
546         ze0(ji,jj) = zz0 * qsr(ji,jj)   ;   zeR(ji,jj) = zz1 * qsr(ji,jj)    ! IR    ; Red
547         zeG(ji,jj) = zz1 * qsr(ji,jj)   ;   zeB(ji,jj) = zz1 * qsr(ji,jj)    ! Green ; Blue
548         zeT(ji,jj) =       qsr(ji,jj)                                        ! Total
549      END_2D
550      !
551      !                             !=  interior light  =!
552      !
553      DO jk = 1, nk0                      !* near surface layers *!   (< ~12 meters : IR + RGB )
554          DO_2D( 0, 0, 0, 0 )
555            ze3t = e3t(ji,jj,jk,Kmm)
556            zze0 = ze0(ji,jj) * EXP( - ze3t * r1_si0 )   ;   zzeR = zeR(ji,jj) * EXP( - ze3t * r1_LR )   ! IR    ; Red  at jk+1 w-level
557            zzeG = zeG(ji,jj) * EXP( - ze3t * r1_LG  )   ;   zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB )   ! Green ; Blue      -      -
558            zzeT = ( zze0 + zzeB + zzeG + zzeR ) * wmask(ji,jj,jk+1)                                     ! Total             -      -
559!!st7-9            zzeT = ( zze0 + zzeR + zzeG + zzeB ) * wmask(ji,jj,jk+1)                                     ! Total             -      -
560#if defined key_RK3
561            !                                               ! RK3 : temperature trend at jk t-level
562            pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t
563#else
564            !                                               ! MLF : heat content trend due to Qsr flux (qsr_hc)
565            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT )
566#endif
567            ze0(ji,jj) = zze0   ;   zeR(ji,jj) = zzeR           ! IR    ; Red  store at jk+1 w-level
568            zeG(ji,jj) = zzeG   ;   zeB(ji,jj) = zzeB           ! Green ; Blue   -        -      -
569            zeT(ji,jj) = zzeT                                   ! total          -        -      -
570         END_2D
571!!stbug         IF( jk == 1 ) THEN               !* sea-ice *!   store the 1st level attenuation coeff.
572!!stbug            WHERE( qsr(A2D(0)) /= 0._wp )   ;   fraqsr_1lev(A2D(0)) = 1._wp - zeT(A2D(0)) / qsr(A2D(0))
573!!stbug            ELSEWHERE                       ;   fraqsr_1lev(A2D(0)) = 1._wp
574!!stbug            END WHERE
575!!stbug         ENDIF
576      END DO
577      !
578      IF( kt == nit000 ) THEN
579         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'
580         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'
581         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'
582         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'
583         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'
584      ENDIF
585      !
586      DO jk = nk0+1, nkR                  !* down to Red extinction *!   (< ~71 meters : RGB , IR removed from calculation)
587          DO_2D( 0, 0, 0, 0 )
588            ze3t = e3t(ji,jj,jk,Kmm)
589            zzeR = zeR(ji,jj) * EXP( - ze3t * r1_LR )                                                 ! Red          at jk+1 w-level
590            zzeG = zeG(ji,jj) * EXP( - ze3t * r1_LG )   ;   zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB ) ! Green ; Blue      -      -
591            zzeT = ( zzeB + zzeG + zzeR ) * wmask(ji,jj,jk+1)                                         ! Total             -      -
592!!st7-11            zzeT = ( zzeR + zzeG + zzeB ) * wmask(ji,jj,jk+1)                                         ! Total             -      -
593#if defined key_RK3
594            !                                               ! RK3 : temperature trend at jk t-level
595            pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t
596#else
597            !                                               ! MLF : heat content trend due to Qsr flux (qsr_hc)
598            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT )
599#endif
600            zeR(ji,jj) = zzeR                                   ! Red          store at jk+1 w-level
601            zeG(ji,jj) = zzeG   ;   zeB(ji,jj) = zzeB           ! Green ; Blue   -        -      -
602            zeT(ji,jj) = zzeT                                   ! total          -        -      -
603         END_2D
604      END DO
605      !
606      IF( kt == nit000 ) THEN
607         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'
608         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'
609         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'
610         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'
611      ENDIF
612      !
613      DO jk = nkR+1, nkG                  !* down to Green extinction *!   (< ~350 m : GB , IR+R removed from calculation)
614         DO_2D( 0, 0, 0, 0 )
615            ze3t = e3t(ji,jj,jk,Kmm)
616            zzeG = zeG(ji,jj) * EXP( - ze3t * r1_LG )   ;   zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB ) ! Green ; Blue at jk+1 w-level
617            zzeT = ( zzeG + zzeB ) * wmask(ji,jj,jk+1)                                                ! Total             -      -
618#if defined key_RK3
619            !                                               ! RK3 : temperature trend at jk t-level
620            pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t
621#else
622            !                                               ! MLF : heat content trend due to Qsr flux (qsr_hc)
623            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT )
624#endif
625            zeG(ji,jj) = zzeG   ;   zeB(ji,jj) = zzeB             ! Green ; Blue store at jk+1 w-level
626            zeT(ji,jj) = zzeT                                     ! total          -        -      -
627         END_2D
628      END DO
629      !
630      IF( kt == nit000 ) THEN
631         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'
632         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'
633         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'
634      ENDIF
635      !
636      DO jk = nkG+1, nkB                  !* down to Blue extinction *!   (< ~1300 m : B , IR+RG removed from calculation)
637         DO_2D( 0, 0, 0, 0 )
638            ze3t = e3t(ji,jj,jk,Kmm)
639            zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB )             ! Blue at jk+1 w-level
640            zzeT = ( zzeB ) * wmask(ji,jj,jk+1)                   ! Total     -      -
641#if defined key_RK3
642            !                                               ! RK3 : temperature trend at jk t-level
643            pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t
644#else
645            !                                               ! MLF : heat content trend due to Qsr flux (qsr_hc)
646            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT )
647#endif
648            zeB(ji,jj) = zzeB                                     ! Blue store at jk+1 w-level
649            zeT(ji,jj) = zzeT                                     ! total  -        -      -
650         END_2D
651      END DO
652      !
653      IF( kt == nit000 ) THEN
654         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'
655      ENDIF
656      !
657   END SUBROUTINE qsr_RGB
658
659
660   SUBROUTINE qsr_2BD( Kmm, pts, Krhs )
661      !!----------------------------------------------------------------------
662      !!                  ***  ROUTINE qsr_2BD  ***
663      !!
664      !! ** Purpose :   2 bands (IR+visible) solar radiation with constant chlorophyll
665      !!
666      !! ** Method  : The profile of the solar radiation within the ocean is defined
667      !!      through 2 wavebands (rn_si0,rn_si1) a ratio rn_abs for IR absorbtion.
668      !!      Considering the 2 wavebands case:
669      !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
670      !!         The temperature trend associated with the solar radiation penetration
671      !!         is given by : zta = 1/e3t dk[ I ] / (rho0*Cp)
672      !!         At the bottom, boudary condition for the radiation is no flux :
673      !!      all heat which has not been absorbed in the above levels is put
674      !!      in the last ocean level.
675      !!         The computation is only done down to the level where
676      !!      I(k) < 1.e-15 W/m2 (i.e. over the top nk levels) .
677      !!
678      !! ** Action  : - update ta with the penetrative solar radiation trend
679      !!              - send  trend for further diagnostics (l_trdtra=T)
680      !!
681      !! Reference  : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
682      !!              Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
683      !!              Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562
684      !!----------------------------------------------------------------------
685      INTEGER,                                   INTENT(in   ) ::   Kmm, Krhs   ! ocean time-step and time-level indices
686      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) ::   pts         ! active tracers and RHS of tracer equation
687      !!
688      INTEGER  ::   ji, jj, jk               ! dummy loop indices
689      REAL(wp) ::   zzatt                    !    -         -
690      REAL(wp) ::   zz0 , zz1 , ze3t         !    -         -
691      REAL(wp), DIMENSION(A2D(0)) ::   zatt
692      !!----------------------------------------------------------------------
693      !     
694      !                       !======================!
695      !                       !==  2-bands fluxes  ==!
696      !                       !======================!
697      !
698      zz0 =           rn_abs   * r1_rho0_rcp       ! surface equi-partition in 2-bands
699      zz1 = ( 1._wp - rn_abs ) * r1_rho0_rcp
700      !
701      zatt(A2D(0)) = r1_rho0_rcp                   !* surface value *!
702      !
703      DO_2D( 0, 0, 0, 0 )
704         zatt(ji,jj) = (  zz0 * EXP( -gdepw(ji,jj,1,Kmm)*r1_si0 ) + zz1 * EXP( -gdepw(ji,jj,1,Kmm)*r1_si1 )  )
705      END_2D
706      !
707!!st      IF(lwp) WRITE(numout,*) 'level = ', 1, ' qsr max = ' , MAXVAL(zatt)*rho0_rcp, ' W/m2', ' qsr min = ' , MINVAL(zatt)*rho0_rcp, ' W/m2'
708      !
709      DO jk = 1, nk0                               !* near surface layers *!   (< ~14 meters : IR + visible light )
710         DO_2D( 0, 0, 0, 0 )
711            ze3t  = e3t(ji,jj,jk,Kmm)                    ! light attenuation at jk+1 w-level (divided by rho0_rcp)
712            zzatt = (   zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*r1_si0 )     &
713               &      + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*r1_si1 )   ) * wmask(ji,jj,jk+1)
714#if defined key_RK3
715            !                                            ! RK3 : temperature trend at jk t-level
716            pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + qsr(ji,jj) * ( zatt(ji,jj) - zzatt ) / ze3t
717#else
718            !                                            ! MLF : heat content trend due to Qsr flux (qsr_hc)
719            qsr_hc(ji,jj,jk) =  qsr(ji,jj) * ( zatt(ji,jj) - zzatt )
720#endif
721            zatt(ji,jj) = zzatt                          ! save for the next level computation
722         END_2D
723!!stbug         !                                         !* sea-ice *!   store the 1st level attenuation coeff.
724!!stbug         IF( jk == 1 )   fraqsr_1lev(A2D(0)) = 1._wp - zatt(A2D(0)) * rho0_rcp
725      END DO
726!!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'
727      !
728      DO jk = nk0+1, nkV                           !* deeper layers *!   (visible light only)
729         DO_2D( 0, 0, 0, 0 )
730            ze3t  = e3t(ji,jj,jk,Kmm)                    ! light attenuation at jk+1 w-level (divided by rho0_rcp)
731            zzatt = (   zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*r1_si1 )   ) * wmask(ji,jj,jk+1)
732#if defined key_RK3
733            !                                            ! RK3 : temperature trend at jk t-level
734            pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + qsr(ji,jj) * ( zatt(ji,jj) - zzatt ) / ze3t
735#else
736            !                                            ! MLF : heat content trend due to Qsr flux (qsr_hc)
737            qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zatt(ji,jj) - zzatt )
738#endif
739            zatt(ji,jj) = zzatt                       ! save for the next level computation
740         END_2D
741      END DO     
742      !
743!!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'
744   END SUBROUTINE qsr_2bd
745
746
747   FUNCTION qsr_ext_lev( pL, pfr ) RESULT( klev )
748      !!----------------------------------------------------------------------
749      !!                 ***  ROUTINE trc_oce_ext_lev  ***
750      !!       
751      !! ** Purpose :   compute the maximum level of light penetration
752      !!         
753      !! ** Method  :   the function provides the level at which irradiance, I,
754      !!              has a negligible effect on temperature.
755      !!                T(n+1)-T(n) = ∆t dk[I] / ( rho0 Cp e3t_k ) 
756      !!              I(k) has a negligible effect on temperature at level k if:
757      !!                ∆t I(k) / ( rho0*Cp*e3t_k ) <= 1.e-15 °C
758      !!              with I(z) = Qsr*pfr*EXP(-z/L), therefore :
759      !!                z >= L * LOG( 1.e-15 * rho0*Cp*e3t_k / ( ∆t*Qsr*pfr ) )
760      !!              with Qsr being the maximum normal surface irradiance at sea
761      !!              level (~1000 W/m2).
762      !!                # pL is the longest depth of extinction:
763      !!                   - pL = 23.00 m (2 bands case)
764      !!                   - pL = 48.24 m (3 bands case: blue waveband & 0.03 mg/m2 for the chlorophyll)
765      !!                # pfr is the fraction of solar radiation which penetrates,
766      !!                considering Qsr=1000 W/m2 and rn_abs = 0.58:
767      !!                   - Qsr*pfr0 = Qsr *    rn_abs    = 580 W/m2   (top absorbtion)
768      !!                   - Qsr*pfr1 = Qsr * (1-rn_abs)   = 420 W/m2 (2 bands case)
769      !!                   - Qsr*pfr1 = Qsr * (1-rn_abs)/3 = 140 W/m2 (3 bands case & equi-partition)
770      !!
771      !!----------------------------------------------------------------------
772      INTEGER              ::   klev   ! result: maximum level of light penetration
773      REAL(wp), INTENT(in) ::   pL     ! depth of extinction
774      REAL(wp), INTENT(in) ::   pfr    ! frac. solar radiation which penetrates
775      !
776      INTEGER  ::   jk                 ! dummy loop index
777      REAL(wp) ::   zcoef              ! local scalar
778      REAL(wp) ::   zhext              ! deepest depth until which light penetrates
779      REAL(wp) ::   ze3t , zdw         ! max( e3t_k ) and min( w-depth_k+1 )
780      REAL(wp) ::   zprec = 10.e-15_wp ! required precision
781      REAL(wp) ::   zQmax= 1000._wp    ! maximum normal surface irradiance at sea level (W/m2)
782      !!----------------------------------------------------------------------
783      !
784      zQmax = 1000._wp     ! maximum normal surface irradiance at sea level (W/m2)
785      !
786      zcoef    =  zprec * rho0_rcp / ( rDt * zQmax * pfr)
787      !
788      IF( ln_zco .OR. ln_zps ) THEN      ! z- or zps coordinate (use 1D ref vertcial coordinate)
789         klev = jpkm1                              ! Level of light extinction zco / zps
790         DO jk = jpkm1, 1, -1
791            zdw  = gdepw_1d(jk+1)                  ! max w-depth at jk+1 level
792            ze3t =   e3t_1d(jk  )                  ! minimum e3t at jk   level
793            zhext =  - pL * LOG( zcoef * ze3t )    ! extinction depth
794            IF( zdw >= zhext )   klev = jk         ! last T-level reached by Qsr
795         END DO
796      ELSE                               ! s- or s-z- coordinate (use 3D vertical coordinate)
797         klev = jpkm1                              ! Level of light extinction
798         DO jk = jpkm1, 1, -1    !
799            IF( SUM( tmask(:,:,jk) ) > 0 ) THEN    ! ocean point at that level
800               zdw  = MAXVAL( gdepw_0(:,:,jk+1) *       wmask(:,:,jk)         )    ! max w-depth at jk+1 level
801               ze3t = MINVAL(   e3t_0(:,:,jk  ) , mask=(wmask(:,:,jk+1)==1)   )    ! minimum e3t at jk   level
802               zhext =  - pL * LOG( zcoef * ze3t )                                 ! extinction depth
803               IF( zdw >= zhext )   klev = jk                                      ! last T-level reached by Qsr
804            ELSE                                   ! only land point at level jk
805               klev = jk                                                           ! local domain sea-bed level
806            ENDIF
807         END DO
808         CALL mpp_max('tra_qsr', klev)             ! needed for reproducibility   !!st may be modified to avoid this comm.
809         !                                                                        !!st use ssmask to remove the comm ?
810      ENDIF
811      !
812!!st      IF(lwp) WRITE(numout,*) '                level of e3t light extinction = ', klev, ' ref depth = ', gdepw_1d(klev+1), ' m'
813   END FUNCTION qsr_ext_lev
814
815
816   SUBROUTINE tra_qsr_init
817      !!----------------------------------------------------------------------
818      !!                  ***  ROUTINE tra_qsr_init  ***
819      !!
820      !! ** Purpose :   Initialization for the penetrative solar radiation
821      !!
822      !! ** Method  :   The profile of solar radiation within the ocean is set
823      !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
824      !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
825      !!      default values correspond to clear water (type I in Jerlov'
826      !!      (1968) classification.
827      !!         called by tra_qsr at the first timestep (nit000)
828      !!
829      !! ** Action  : - initialize rn_si0, rn_si1 and rn_abs
830      !!
831      !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
832      !!----------------------------------------------------------------------
833      INTEGER  ::   ji, jj, jk                  ! dummy loop indices
834      INTEGER  ::   ios, ierror, ioptio   ! local integer
835      REAL(wp) ::   zLB, zLG, zLR         ! local scalar
836      REAL(wp) ::   zVlp, zchl            !   -      -
837      !
838      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files
839      TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read
840      !!
841      NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio,  &
842         &                  nn_chldta, rn_abs, rn_si0, rn_si1
843      !!----------------------------------------------------------------------
844      !
845      READ  ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901)
846901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_qsr in reference namelist' )
847      READ  ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902)
848902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist' )
849      IF(lwm) WRITE ( numond, namtra_qsr )
850      !
851      IF(lwp) THEN            !**  control print  **!
852         WRITE(numout,*)
853         WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
854         WRITE(numout,*) '~~~~~~~~~~~~'
855         WRITE(numout,*) '   Namelist namtra_qsr : set the parameter of penetration'
856         WRITE(numout,*) '      RGB (Red-Green-Blue) light penetration       ln_qsr_rgb = ', ln_qsr_rgb
857         WRITE(numout,*) '      2 band               light penetration       ln_qsr_2bd = ', ln_qsr_2bd
858         WRITE(numout,*) '      bio-model            light penetration       ln_qsr_bio = ', ln_qsr_bio
859         WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)        nn_chldta  = ', nn_chldta
860         WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs     = ', rn_abs
861         WRITE(numout,*) '      RGB & 2 bands: shortess attenuation depth    rn_si0     = ', rn_si0
862         WRITE(numout,*) '      2 bands: longest attenuation depth           rn_si1     = ', rn_si1
863         WRITE(numout,*)
864      ENDIF
865      !
866      ioptio = 0              !**  Parameter control  **!
867      IF( ln_qsr_rgb  )   ioptio = ioptio + 1
868      IF( ln_qsr_2bd  )   ioptio = ioptio + 1
869      IF( ln_qsr_bio  )   ioptio = ioptio + 1
870      !
871      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE type of light penetration in namelist namtra_qsr',  &
872         &                               ' 2 bands, 3 RGB bands or bio-model light penetration' )
873      !
874      IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr = np_RGB
875      IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr = np_RGBc
876      IF( ln_qsr_2bd                      )   nqsr = np_2BD
877      IF( ln_qsr_bio                      )   nqsr = np_BIO
878      !
879      !                       !**  Initialisation  **!
880      !
881      !                                !==  Infrared attenuation  ==!   (all schemes)
882      !                                !============================!
883      !
884      r1_si0 = 1._wp / rn_si0                ! inverse of infrared attenuation length
885      !
886      nk0 = qsr_ext_lev( rn_si0, rn_abs )    ! level of light extinction
887      !
888      IF(lwp) WRITE(numout,*) '   ==>>>   Infrared light attenuation'
889      IF(lwp) WRITE(numout,*) '              level of infrared extinction = ', nk0, ' ref depth = ', gdepw_1d(nk0+1), ' m'
890      IF(lwp) WRITE(numout,*)
891      !
892      SELECT CASE( nqsr )
893      !
894      CASE( np_RGBc, np_RGB )          !==  Red-Green-Blue light attenuation  ==!   (Chl data or constant)
895         !                             !========================================!
896         !
897         IF( nqsr == np_RGB ) THEN   ;   zchl   = 0.05        ! constant Chl value
898         ELSE                        ;   zchl   = 0.03        ! minimum  Chl value
899         ENDIF
900         zchl   = MAX( 0.03_wp , MIN( zchl , 10._wp) )     ! NB. make sure that chosen value verifies: 0.03 < zchl < 10
901         nc_rgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )    ! Convert Chl value to attenuation coefficient look-up table index
902         !
903         CALL trc_oce_rgb( rkrgb )                 ! tabulated attenuation coef.
904         !
905         zVlp =  ( 1._wp - rn_abs ) / 3._wp        ! visible light equi-partition
906         !
907         !     1 / length          !   attenuation  length   !         attenuation level
908         r1_LR = rkrgb(3,nc_rgb)   ;   zLR = 1._wp / r1_LR   ;   nkR = qsr_ext_lev( zLR, zVlp )   ! Red   
909         r1_LG = rkrgb(2,nc_rgb)   ;   zLG = 1._wp / r1_LG   ;   nkG = qsr_ext_lev( zLG, zVlp )   ! Green
910         r1_LB = rkrgb(1,nc_rgb)   ;   zLB = 1._wp / r1_LB   ;   nkB = qsr_ext_lev( zLB, zVlp )   ! Blue
911         !
912         nkV = nkB                                 ! maximum level of light penetration
913         !
914         IF( nqsr == np_RGB ) THEN
915            IF(lwp) WRITE(numout,*) '   ==>>>   RGB:  light attenuation with a constant Chlorophyll = ', zchl
916         ELSE
917            IF(lwp) WRITE(numout,*) '   ==>>>   RGB:  light attenuation using Chlorophyll data with min(Chl) = ', zchl
918         ENDIF           
919         IF(lwp) WRITE(numout,*) '                 level of Red   extinction = ', nkR, ' ref depth = ', gdepw_1d(nkR+1), ' m'
920         IF(lwp) WRITE(numout,*) '                 level of Green extinction = ', nkG, ' ref depth = ', gdepw_1d(nkG+1), ' m'
921         IF(lwp) WRITE(numout,*) '                 level of Blue  extinction = ', nkB, ' ref depth = ', gdepw_1d(nkB+1), ' m'
922         IF(lwp) WRITE(numout,*)
923         !
924         IF( nqsr == np_RGBc ) THEN                ! Chl data : set sf_chl structure
925            IF(lwp) WRITE(numout,*) '   ==>>>   Chlorophyll read in a file'
926            ALLOCATE( sf_chl(1), STAT=ierror )
927            IF( ierror > 0 ) THEN
928               CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' )   ;   RETURN
929            ENDIF
930            ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1) )
931            IF( sn_chl%ln_tint )   ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) )
932            !                                        ! fill sf_chl with sn_chl and control print
933            CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',                             &
934               &           'Solar penetration function of read chlorophyll', 'namtra_qsr' , no_print )
935         ENDIF
936         !
937      CASE( np_2BD )                   !==  2 bands light attenuation (IR+ visible light) ==!
938         !
939         !
940         r1_si1 = 1._wp / rn_si1                     ! inverse of visible light attenuation
941         zVlp =  ( 1._wp - rn_abs )                  ! visible light partition
942         nkV  = qsr_ext_lev( rn_si1, zVlp )          ! level of visible light extinction
943         !
944         IF(lwp) WRITE(numout,*) '   ==>>>   2 bands attenuation (Infrared + Visible light) '
945         IF(lwp) WRITE(numout,*) '                level of visible light extinction = ', nkV, ' ref depth = ', gdepw_1d(nkV+1), ' m'
946         IF(lwp) WRITE(numout,*)
947         !
948      CASE( np_BIO )                   !==  BIO light penetration  ==!
949         !
950         IF(lwp) WRITE(numout,*) '   ==>>>   bio-model light penetration'
951         IF( .NOT.lk_top )   CALL ctl_stop( 'No bio model : ln_qsr_bio = true impossible ' )
952         !
953         CALL trc_oce_rgb( rkrgb )                 ! tabulated attenuation coef.
954         !
955         nkV = trc_oce_ext_lev( r_si2, 33._wp )    ! maximum level of light extinction
956         !
957         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nkV, ' ref depth = ', gdepw_1d(nkV+1), ' m'
958         !
959      END SELECT
960      !
961      nksr = nkV       ! name of max level of light extinction used in traatf(_qco).F90
962      !
963#if ! defined key_RK3
964      qsr_hc(:,:,:) = 0._wp      ! MLF : now qsr heat content set to zero where it will not be computed
965#endif
966      !
967      !                          ! Sea-ice :   1st ocean level attenuation coefficient (used in sbcssm)
968      IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN
969         CALL iom_get( numror, jpdom_auto, 'fraqsr_1lev'  , fraqsr_1lev  )
970      ELSE
971         fraqsr_1lev(:,:) = 1._wp   ! default : no penetration
972      ENDIF
973      !
974   END SUBROUTINE tra_qsr_init
975
976   !!======================================================================
977END MODULE traqsr
Note: See TracBrowser for help on using the repository browser.