2020WP/ENHANCE-10_acc_fix_traqsr: traqsr_lomem.F90

File traqsr_lomem.F90, 21.4 KB (added by acc, 5 months ago)
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   !!----------------------------------------------------------------------
16
17   !!----------------------------------------------------------------------
18   !!   tra_qsr       : temperature trend due to the penetration of solar radiation
19   !!   tra_qsr_init  : initialization of the qsr penetration
20   !!----------------------------------------------------------------------
21   USE oce            ! ocean dynamics and active tracers
22   USE phycst         ! physical constants
23   USE dom_oce        ! ocean space and time domain
24   USE sbc_oce        ! surface boundary condition: ocean
25   USE trc_oce        ! share SMS/Ocean variables
26   USE trd_oce        ! trends: ocean variables
27   USE trdtra         ! trends manager: tracers
28   !
29   USE in_out_manager ! I/O manager
30   USE prtctl         ! Print control
31   USE iom            ! I/O library
32   USE fldread        ! read input fields
33   USE restart        ! ocean restart
34   USE lib_mpp        ! MPP library
35   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
36   USE timing         ! Timing
37
38   IMPLICIT NONE
39   PRIVATE
40
41   PUBLIC   tra_qsr       ! routine called by step.F90 (ln_traqsr=T)
42   PUBLIC   tra_qsr_init  ! routine called by nemogcm.F90
43
44   !                                 !!* Namelist namtra_qsr: penetrative solar radiation
45   LOGICAL , PUBLIC ::   ln_traqsr    !: light absorption (qsr) flag
46   LOGICAL , PUBLIC ::   ln_qsr_rgb   !: Red-Green-Blue light absorption flag 
47   LOGICAL , PUBLIC ::   ln_qsr_2bd   !: 2 band         light absorption flag
48   LOGICAL , PUBLIC ::   ln_qsr_bio   !: bio-model      light absorption flag
49   INTEGER , PUBLIC ::   nn_chldta    !: use Chlorophyll data (=1) or not (=0)
50   REAL(wp), PUBLIC ::   rn_abs       !: fraction absorbed in the very near surface (RGB & 2 bands)
51   REAL(wp), PUBLIC ::   rn_si0       !: very near surface depth of extinction      (RGB & 2 bands)
52   REAL(wp), PUBLIC ::   rn_si1       !: deepest depth of extinction (water type I)       (2 bands)
53   !
54   INTEGER , PUBLIC ::   nksr         !: levels below which the light cannot penetrate (depth larger than 391 m)
55 
56   INTEGER, PARAMETER ::   np_RGB  = 1   ! R-G-B     light penetration with constant Chlorophyll
57   INTEGER, PARAMETER ::   np_RGBc = 2   ! R-G-B     light penetration with Chlorophyll data
58   INTEGER, PARAMETER ::   np_2BD  = 3   ! 2 bands   light penetration
59   INTEGER, PARAMETER ::   np_BIO  = 4   ! bio-model light penetration
60   !
61   INTEGER  ::   nqsr    ! user choice of the type of light penetration
62   REAL(wp) ::   xsi0r   ! inverse of rn_si0
63   REAL(wp) ::   xsi1r   ! inverse of rn_si1
64   !
65   REAL(wp) , DIMENSION(3,61)           ::   rkrgb    ! tabulated attenuation coefficients for RGB absorption
66   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read)
67
68   !! * Substitutions
69#  include "do_loop_substitute.h90"
70   !!----------------------------------------------------------------------
71   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
72   !! $Id: traqsr.F90 12489 2020-02-28 15:55:11Z davestorkey $
73   !! Software governed by the CeCILL license (see ./LICENSE)
74   !!----------------------------------------------------------------------
75CONTAINS
76
77   SUBROUTINE tra_qsr( kt, Kmm, pts, Krhs )
78      !!----------------------------------------------------------------------
79      !!                  ***  ROUTINE tra_qsr  ***
80      !!
81      !! ** Purpose :   Compute the temperature trend due to the solar radiation
82      !!              penetration and add it to the general temperature trend.
83      !!
84      !! ** Method  : The profile of the solar radiation within the ocean is defined
85      !!      through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) and a ratio rn_abs
86      !!      Considering the 2 wavebands case:
87      !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
88      !!         The temperature trend associated with the solar radiation penetration
89      !!         is given by : zta = 1/e3t dk[ I ] / (rho0*Cp)
90      !!         At the bottom, boudary condition for the radiation is no flux :
91      !!      all heat which has not been absorbed in the above levels is put
92      !!      in the last ocean level.
93      !!         The computation is only done down to the level where
94      !!      I(k) < 1.e-15 W/m2 (i.e. over the top nksr levels) .
95      !!
96      !! ** Action  : - update ta with the penetrative solar radiation trend
97      !!              - send  trend for further diagnostics (l_trdtra=T)
98      !!
99      !! Reference  : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
100      !!              Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
101      !!              Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562
102      !!----------------------------------------------------------------------
103      INTEGER,                                   INTENT(in   ) :: kt            ! ocean time-step
104      INTEGER,                                   INTENT(in   ) :: Kmm, Krhs     ! time level indices
105      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts           ! active tracers and RHS of tracer equation
106      !
107      INTEGER  ::   ji, jj, jk               ! dummy loop indices
108      INTEGER  ::   irgb                     ! local integers
109      REAL(wp) ::   zchl, zcoef, z1_2        ! local scalars
110      REAL(wp) ::   zc0 , zc1 , zc2 , zc3    !    -         -
111      REAL(wp) ::   zzc0, zzc1, zzc2, zzc3   !    -         -
112      REAL(wp) ::   zz0 , zz1 , ze3t, zlui   !    -         -
113      REAL(wp) ::   zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze
114      REAL(wp) ::   zlogc, zlogc2, zlogc3 
115      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   :: zeka, zekb, zekg, zekr
116      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   :: ze0, ze1, ze2, ze3
117      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, zetot, ztmp3d
118      !!----------------------------------------------------------------------
119      !
120      IF( ln_timing )   CALL timing_start('tra_qsr')
121      !
122      IF( kt == nit000 ) THEN
123         IF(lwp) WRITE(numout,*)
124         IF(lwp) WRITE(numout,*) 'tra_qsr : penetration of the surface solar radiation'
125         IF(lwp) WRITE(numout,*) '~~~~~~~'
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      !                         !-----------------------------------!
134      !                         !  before qsr induced heat content  !
135      !                         !-----------------------------------!
136      IF( kt == nit000 ) THEN          !==  1st time step  ==!
137         IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0  .AND. .NOT.l_1st_euler ) THEN    ! read in restart
138            IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field read in the restart file'
139            z1_2 = 0.5_wp
140            CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b, ldxios = lrxios )   ! before heat content trend due to Qsr flux
141         ELSE                                           ! No restart or restart not found: Euler forward time stepping
142            z1_2 = 1._wp
143            qsr_hc_b(:,:,:) = 0._wp
144         ENDIF
145      ELSE                             !==  Swap of qsr heat content  ==!
146         z1_2 = 0.5_wp
147         qsr_hc_b(:,:,:) = qsr_hc(:,:,:)
148      ENDIF
149      !
150      !                         !--------------------------------!
151      SELECT CASE( nqsr )       !  now qsr induced heat content  !
152      !                         !--------------------------------!
153      !
154      CASE( np_BIO )                   !==  bio-model fluxes  ==!
155         !
156         DO jk = 1, nksr
157            qsr_hc(:,:,jk) = r1_rho0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) )
158         END DO
159         !
160      CASE( np_RGB , np_RGBc )         !==  R-G-B fluxes  ==!
161         !
162         ALLOCATE( ze0 (jpi,jpj)           , ze1 (jpi,jpj) ,   &
163            &      ze2 (jpi,jpj)           , ze3 (jpi,jpj) ,   &
164            &      ztmp3d(jpi,jpj,nksr + 1)                     )
165         !
166         IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll
167            CALL fld_read( kt, 1, sf_chl )         ! Read Chl data and provides it at the current time step
168            ! Separation in R-G-B depending of the surface Chl
169            DO_3D_00_00 ( 1, nksr + 1 )
170               zchl    = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) )
171               zCtot   = 40.6  * zchl**0.459
172               zze     = 568.2 * zCtot**(-0.746)
173               IF( zze > 102. ) zze = 200.0 * zCtot**(-0.293)
174               zpsi    = gdepw(ji,jj,jk,Kmm) / zze
175               !
176               zlogc   = LOG( zchl )
177               zlogc2  = zlogc * zlogc
178               zlogc3  = zlogc * zlogc * zlogc
179               zCb     = 0.768 + 0.087 * zlogc - 0.179 * zlogc2 - 0.025 * zlogc3
180               zCmax   = 0.299 - 0.289 * zlogc + 0.579 * zlogc2
181               zpsimax = 0.6   - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3
182               zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2
183               zCze    = 1.12  * (zchl)**0.803 
184               !
185               ! NB. make sure zchl value is such that: zchl = MIN( 10. , MAX( 0.03, zchl ) )
186               zchl = MIN( 10. , MAX( 0.03, zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) / zdelpsi )**2 ) ) ) )
187               ! Convert chlorophyll value to attenuation coefficient look-up table index
188               ztmp3d(ji,jj,jk) = 41 + 20.*LOG10(zchl) + 1.e-15
189            END_3D
190         ELSE                                !* constant chlorophyll
191            zchl = 0.05
192            ! NB. make sure constant value is such that:
193            zchl = MIN( 10. , MAX( 0.03, zchl ) )
194            ! Convert chlorophyll value to attenuation coefficient look-up table index
195            zlui = 41 + 20.*LOG10(zchl) + 1.e-15
196            DO jk = 1, nksr + 1
197               ztmp3d(:,:,jk) = zlui 
198            END DO
199         ENDIF
200         !
201         zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B
202         DO_2D_00_00
203            ze0(ji,jj) = rn_abs * qsr(ji,jj)
204            ze1(ji,jj) = zcoef  * qsr(ji,jj)
205            ze2(ji,jj) = zcoef  * qsr(ji,jj)
206            ze3(ji,jj) = zcoef  * qsr(ji,jj)
207            ! store the surface SW radiation; re-use the surface ztmp3d array
208            ! since the surface attenuation coefficient is not used
209            ztmp3d(ji,jj,1) =       qsr(ji,jj)
210         END_2D
211         !
212         !* interior equi-partition in R-G-B depending of vertical profile of Chl
213         DO_3D_00_00( 2, nksr+1 )
214            irgb = NINT( ztmp3d(ji,jj,jk) )
215            ze3t = e3t(ji,jj,jk-1,Kmm)
216            ze0(ji,jj) = ze0(ji,jj) * EXP( - ze3t * xsi0r )
217            ze1(ji,jj) = ze1(ji,jj) * EXP( - ze3t * rkrgb(1,irgb) )
218            ze2(ji,jj) = ze2(ji,jj) * EXP( - ze3t * rkrgb(2,irgb) )
219            ze3(ji,jj) = ze3(ji,jj) * EXP( - ze3t * rkrgb(3,irgb) )
220            ! store the SW radiation penetrating to this location
221            ! re-use the ztmp3d array since the attenuation coefficient
222            ! at this point will not be needed again
223            ztmp3d(ji,jj,jk) = ( ze0(ji,jj) + ze1(ji,jj) + ze2(ji,jj) + ze3(ji,jj) ) * wmask(ji,jj,jk)
224         END_3D
225         !
226         DO_3D_00_00( 1, nksr )
227            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( ztmp3d(ji,jj,jk) - ztmp3d(ji,jj,jk+1) )
228         END_3D
229         !
230         DEALLOCATE( ze0 , ze1 , ze2 , ze3 , ztmp3d ) 
231         !
232      CASE( np_2BD  )            !==  2-bands fluxes  ==!
233         !
234         zz0 =        rn_abs   * r1_rho0_rcp      ! surface equi-partition in 2-bands
235         zz1 = ( 1. - rn_abs ) * r1_rho0_rcp
236         DO_3D_00_00( 1, nksr )
237            zc0 = zz0 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi1r )
238            zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r )
239            qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 
240         END_3D
241         !
242      END SELECT
243      !
244      !                          !-----------------------------!
245      DO_3D_00_00( 1, nksr )
246         pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
247            &                      + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t(ji,jj,jk,Kmm)
248      END_3D
249      !
250      ! sea-ice: store the 1st ocean level attenuation coefficient
251      DO_2D_00_00
252         IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) )
253         ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp
254         ENDIF
255      END_2D
256      CALL lbc_lnk( 'traqsr', fraqsr_1lev(:,:), 'T', 1._wp )
257      !
258      IF( iom_use('qsr3d') ) THEN      ! output the shortwave Radiation distribution
259         ALLOCATE( zetot(jpi,jpj,jpk) )
260         zetot(:,:,nksr+1:jpk) = 0._wp     ! below ~400m set to zero
261         DO jk = nksr, 1, -1
262            zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rho0_rcp
263         END DO         
264         CALL iom_put( 'qsr3d', zetot )   ! 3D distribution of shortwave Radiation
265         DEALLOCATE( zetot ) 
266      ENDIF
267      !
268      IF( lrst_oce ) THEN     ! write in the ocean restart file
269         IF( lwxios ) CALL iom_swap(      cwxios_context          )
270         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc     , ldxios = lwxios )
271         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev, ldxios = lwxios ) 
272         IF( lwxios ) CALL iom_swap(      cxios_context          )
273      ENDIF
274      !
275      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
276         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
277         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt )
278         DEALLOCATE( ztrdt ) 
279      ENDIF
280      !                       ! print mean trends (used for debugging)
281      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
282      !
283      IF( ln_timing )   CALL timing_stop('tra_qsr')
284      !
285   END SUBROUTINE tra_qsr
286
287
288   SUBROUTINE tra_qsr_init
289      !!----------------------------------------------------------------------
290      !!                  ***  ROUTINE tra_qsr_init  ***
291      !!
292      !! ** Purpose :   Initialization for the penetrative solar radiation
293      !!
294      !! ** Method  :   The profile of solar radiation within the ocean is set
295      !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
296      !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
297      !!      default values correspond to clear water (type I in Jerlov'
298      !!      (1968) classification.
299      !!         called by tra_qsr at the first timestep (nit000)
300      !!
301      !! ** Action  : - initialize rn_si0, rn_si1 and rn_abs
302      !!
303      !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
304      !!----------------------------------------------------------------------
305      INTEGER  ::   ji, jj, jk                  ! dummy loop indices
306      INTEGER  ::   ios, irgb, ierror, ioptio   ! local integer
307      REAL(wp) ::   zz0, zc0 , zc1, zcoef      ! local scalars
308      REAL(wp) ::   zz1, zc2 , zc3, zchl       !   -      -
309      !
310      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files
311      TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read
312      !!
313      NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio,  &
314         &                  nn_chldta, rn_abs, rn_si0, rn_si1
315      !!----------------------------------------------------------------------
316      !
317      READ  ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901)
318901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_qsr in reference namelist' )
319      !
320      READ  ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 )
321902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist' )
322      IF(lwm) WRITE ( numond, namtra_qsr )
323      !
324      IF(lwp) THEN                ! control print
325         WRITE(numout,*)
326         WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
327         WRITE(numout,*) '~~~~~~~~~~~~'
328         WRITE(numout,*) '   Namelist namtra_qsr : set the parameter of penetration'
329         WRITE(numout,*) '      RGB (Red-Green-Blue) light penetration       ln_qsr_rgb = ', ln_qsr_rgb
330         WRITE(numout,*) '      2 band               light penetration       ln_qsr_2bd = ', ln_qsr_2bd
331         WRITE(numout,*) '      bio-model            light penetration       ln_qsr_bio = ', ln_qsr_bio
332         WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)        nn_chldta  = ', nn_chldta
333         WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs     = ', rn_abs
334         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0     = ', rn_si0
335         WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1     = ', rn_si1
336         WRITE(numout,*)
337      ENDIF
338      !
339      ioptio = 0                    ! Parameter control
340      IF( ln_qsr_rgb  )   ioptio = ioptio + 1
341      IF( ln_qsr_2bd  )   ioptio = ioptio + 1
342      IF( ln_qsr_bio  )   ioptio = ioptio + 1
343      !
344      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE type of light penetration in namelist namtra_qsr',  &
345         &                               ' 2 bands, 3 RGB bands or bio-model light penetration' )
346      !
347      IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr = np_RGB 
348      IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr = np_RGBc
349      IF( ln_qsr_2bd                      )   nqsr = np_2BD
350      IF( ln_qsr_bio                      )   nqsr = np_BIO
351      !
352      !                             ! Initialisation
353      xsi0r = 1._wp / rn_si0
354      xsi1r = 1._wp / rn_si1
355      !
356      SELECT CASE( nqsr )
357      !                               
358      CASE( np_RGB , np_RGBc )         !==  Red-Green-Blue light penetration  ==!
359         !                             
360         IF(lwp)   WRITE(numout,*) '   ==>>>   R-G-B   light penetration '
361         !
362         CALL trc_oce_rgb( rkrgb )                 ! tabulated attenuation coef.
363         !                                   
364         nksr = trc_oce_ext_lev( r_si2, 33._wp )   ! level of light extinction
365         !
366         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
367         !
368         IF( nqsr == np_RGBc ) THEN                ! Chl data : set sf_chl structure
369            IF(lwp) WRITE(numout,*) '   ==>>>   Chlorophyll read in a file'
370            ALLOCATE( sf_chl(1), STAT=ierror )
371            IF( ierror > 0 ) THEN
372               CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' )   ;   RETURN
373            ENDIF
374            ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1)   )
375            IF( sn_chl%ln_tint )   ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) )
376            !                                        ! fill sf_chl with sn_chl and control print
377            CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',   &
378               &           'Solar penetration function of read chlorophyll', 'namtra_qsr' , no_print )
379         ENDIF
380         IF( nqsr == np_RGB ) THEN                 ! constant Chl
381            IF(lwp) WRITE(numout,*) '   ==>>>   Constant Chlorophyll concentration = 0.05'
382         ENDIF
383         !
384      CASE( np_2BD )                   !==  2 bands light penetration  ==!
385         !
386         IF(lwp)  WRITE(numout,*) '   ==>>>   2 bands light penetration'
387         !
388         nksr = trc_oce_ext_lev( rn_si1, 100._wp )    ! level of light extinction
389         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
390         !
391      CASE( np_BIO )                   !==  BIO light penetration  ==!
392         !
393         IF(lwp) WRITE(numout,*) '   ==>>>   bio-model light penetration'
394         IF( .NOT.lk_top )   CALL ctl_stop( 'No bio model : ln_qsr_bio = true impossible ' )
395         !
396      END SELECT
397      !
398      qsr_hc(:,:,:) = 0._wp     ! now qsr heat content set to zero where it will not be computed
399      !
400      ! 1st ocean level attenuation coefficient (used in sbcssm)
401      IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN
402         CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev'  , fraqsr_1lev, ldxios = lrxios  )
403      ELSE
404         fraqsr_1lev(:,:) = 1._wp   ! default : no penetration
405      ENDIF
406      !
407      IF( lwxios ) THEN
408         CALL iom_set_rstw_var_active('qsr_hc_b')
409         CALL iom_set_rstw_var_active('fraqsr_1lev')
410      ENDIF
411      !
412   END SUBROUTINE tra_qsr_init
413
414   !!======================================================================
415END MODULE traqsr