2020WP/ENHANCE-10_acc_fix_traqsr: traqsr_trunk.F90

File traqsr_trunk.F90, 21.3 KB (added by acc, 7 weeks 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                !    -         -
113      REAL(wp) ::   zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze
114      REAL(wp) ::   zlogc, zlogc2, zlogc3 
115      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   :: zekb, zekg, zekr
116      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt
117      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zetot, zchl3d
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( zekb(jpi,jpj)     , zekg(jpi,jpj)     , zekr  (jpi,jpj)     , &
163            &      ze0 (jpi,jpj,jpk) , ze1 (jpi,jpj,jpk) , ze2   (jpi,jpj,jpk) , &
164            &      ze3 (jpi,jpj,jpk) , zea (jpi,jpj,jpk) , zchl3d(jpi,jpj,jpk)   ) 
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            DO jk = 1, nksr + 1
169               DO jj = 2, jpjm1                       ! Separation in R-G-B depending of the surface Chl
170                  DO ji = 2, jpim1
171                     zchl    = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) )
172                     zCtot   = 40.6  * zchl**0.459
173                     zze     = 568.2 * zCtot**(-0.746)
174                     IF( zze > 102. ) zze = 200.0 * zCtot**(-0.293)
175                     zpsi    = gdepw(ji,jj,jk,Kmm) / zze
176                     !
177                     zlogc   = LOG( zchl )
178                     zlogc2  = zlogc * zlogc
179                     zlogc3  = zlogc * zlogc * zlogc
180                     zCb     = 0.768 + 0.087 * zlogc - 0.179 * zlogc2 - 0.025 * zlogc3
181                     zCmax   = 0.299 - 0.289 * zlogc + 0.579 * zlogc2
182                     zpsimax = 0.6   - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3
183                     zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2
184                     zCze    = 1.12  * (zchl)**0.803 
185                     !
186                     zchl3d(ji,jj,jk) = zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) / zdelpsi )**2 ) )
187                  END DO
188                  !
189               END DO
190            END DO
191         ELSE                                !* constant chrlorophyll
192           DO jk = 1, nksr + 1
193              zchl3d(:,:,jk) = 0.05 
194            ENDDO
195         ENDIF
196         !
197         zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B
198         DO_2D_00_00
199            ze0(ji,jj,1) = rn_abs * qsr(ji,jj)
200            ze1(ji,jj,1) = zcoef  * qsr(ji,jj)
201            ze2(ji,jj,1) = zcoef  * qsr(ji,jj)
202            ze3(ji,jj,1) = zcoef  * qsr(ji,jj)
203            zea(ji,jj,1) =          qsr(ji,jj)
204         END_2D
205         !
206         DO jk = 2, nksr+1                   !* interior equi-partition in R-G-B depending of vertical profile of Chl
207            DO_2D_00_00
208               zchl = MIN( 10. , MAX( 0.03, zchl3d(ji,jj,jk) ) )
209               irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
210               zekb(ji,jj) = rkrgb(1,irgb)
211               zekg(ji,jj) = rkrgb(2,irgb)
212               zekr(ji,jj) = rkrgb(3,irgb)
213            END_2D
214
215            DO_2D_00_00
216               zc0 = ze0(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * xsi0r       )
217               zc1 = ze1(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekb(ji,jj) )
218               zc2 = ze2(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekg(ji,jj) )
219               zc3 = ze3(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekr(ji,jj) )
220               ze0(ji,jj,jk) = zc0
221               ze1(ji,jj,jk) = zc1
222               ze2(ji,jj,jk) = zc2
223               ze3(ji,jj,jk) = zc3
224               zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk)
225            END_2D
226         END DO
227         !
228         DO_3D_00_00( 1, nksr )
229            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) )
230         END_3D
231         !
232         DEALLOCATE( zekb , zekg , zekr , ze0 , ze1 , ze2 , ze3 , zea , zchl3d ) 
233         !
234      CASE( np_2BD  )            !==  2-bands fluxes  ==!
235         !
236         zz0 =        rn_abs   * r1_rho0_rcp      ! surface equi-partition in 2-bands
237         zz1 = ( 1. - rn_abs ) * r1_rho0_rcp
238         DO_3D_00_00( 1, nksr )
239            zc0 = zz0 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi1r )
240            zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r )
241            qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 
242         END_3D
243         !
244      END SELECT
245      !
246      !                          !-----------------------------!
247      DO_3D_00_00( 1, nksr )
248         pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
249            &                      + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t(ji,jj,jk,Kmm)
250      END_3D
251      !
252      ! sea-ice: store the 1st ocean level attenuation coefficient
253      DO_2D_00_00
254         IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) )
255         ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp
256         ENDIF
257      END_2D
258      CALL lbc_lnk( 'traqsr', fraqsr_1lev(:,:), 'T', 1._wp )
259      !
260      IF( iom_use('qsr3d') ) THEN      ! output the shortwave Radiation distribution
261         ALLOCATE( zetot(jpi,jpj,jpk) )
262         zetot(:,:,nksr+1:jpk) = 0._wp     ! below ~400m set to zero
263         DO jk = nksr, 1, -1
264            zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rho0_rcp
265         END DO         
266         CALL iom_put( 'qsr3d', zetot )   ! 3D distribution of shortwave Radiation
267         DEALLOCATE( zetot ) 
268      ENDIF
269      !
270      IF( lrst_oce ) THEN     ! write in the ocean restart file
271         IF( lwxios ) CALL iom_swap(      cwxios_context          )
272         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc     , ldxios = lwxios )
273         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev, ldxios = lwxios ) 
274         IF( lwxios ) CALL iom_swap(      cxios_context          )
275      ENDIF
276      !
277      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
278         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
279         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt )
280         DEALLOCATE( ztrdt ) 
281      ENDIF
282      !                       ! print mean trends (used for debugging)
283      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
284      !
285      IF( ln_timing )   CALL timing_stop('tra_qsr')
286      !
287   END SUBROUTINE tra_qsr
288
289
290   SUBROUTINE tra_qsr_init
291      !!----------------------------------------------------------------------
292      !!                  ***  ROUTINE tra_qsr_init  ***
293      !!
294      !! ** Purpose :   Initialization for the penetrative solar radiation
295      !!
296      !! ** Method  :   The profile of solar radiation within the ocean is set
297      !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
298      !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
299      !!      default values correspond to clear water (type I in Jerlov'
300      !!      (1968) classification.
301      !!         called by tra_qsr at the first timestep (nit000)
302      !!
303      !! ** Action  : - initialize rn_si0, rn_si1 and rn_abs
304      !!
305      !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
306      !!----------------------------------------------------------------------
307      INTEGER  ::   ji, jj, jk                  ! dummy loop indices
308      INTEGER  ::   ios, irgb, ierror, ioptio   ! local integer
309      REAL(wp) ::   zz0, zc0 , zc1, zcoef      ! local scalars
310      REAL(wp) ::   zz1, zc2 , zc3, zchl       !   -      -
311      !
312      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files
313      TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read
314      !!
315      NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio,  &
316         &                  nn_chldta, rn_abs, rn_si0, rn_si1
317      !!----------------------------------------------------------------------
318      !
319      READ  ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901)
320901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_qsr in reference namelist' )
321      !
322      READ  ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 )
323902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist' )
324      IF(lwm) WRITE ( numond, namtra_qsr )
325      !
326      IF(lwp) THEN                ! control print
327         WRITE(numout,*)
328         WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
329         WRITE(numout,*) '~~~~~~~~~~~~'
330         WRITE(numout,*) '   Namelist namtra_qsr : set the parameter of penetration'
331         WRITE(numout,*) '      RGB (Red-Green-Blue) light penetration       ln_qsr_rgb = ', ln_qsr_rgb
332         WRITE(numout,*) '      2 band               light penetration       ln_qsr_2bd = ', ln_qsr_2bd
333         WRITE(numout,*) '      bio-model            light penetration       ln_qsr_bio = ', ln_qsr_bio
334         WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)        nn_chldta  = ', nn_chldta
335         WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs     = ', rn_abs
336         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0     = ', rn_si0
337         WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1     = ', rn_si1
338         WRITE(numout,*)
339      ENDIF
340      !
341      ioptio = 0                    ! Parameter control
342      IF( ln_qsr_rgb  )   ioptio = ioptio + 1
343      IF( ln_qsr_2bd  )   ioptio = ioptio + 1
344      IF( ln_qsr_bio  )   ioptio = ioptio + 1
345      !
346      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE type of light penetration in namelist namtra_qsr',  &
347         &                               ' 2 bands, 3 RGB bands or bio-model light penetration' )
348      !
349      IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr = np_RGB 
350      IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr = np_RGBc
351      IF( ln_qsr_2bd                      )   nqsr = np_2BD
352      IF( ln_qsr_bio                      )   nqsr = np_BIO
353      !
354      !                             ! Initialisation
355      xsi0r = 1._wp / rn_si0
356      xsi1r = 1._wp / rn_si1
357      !
358      SELECT CASE( nqsr )
359      !                               
360      CASE( np_RGB , np_RGBc )         !==  Red-Green-Blue light penetration  ==!
361         !                             
362         IF(lwp)   WRITE(numout,*) '   ==>>>   R-G-B   light penetration '
363         !
364         CALL trc_oce_rgb( rkrgb )                 ! tabulated attenuation coef.
365         !                                   
366         nksr = trc_oce_ext_lev( r_si2, 33._wp )   ! level of light extinction
367         !
368         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
369         !
370         IF( nqsr == np_RGBc ) THEN                ! Chl data : set sf_chl structure
371            IF(lwp) WRITE(numout,*) '   ==>>>   Chlorophyll read in a file'
372            ALLOCATE( sf_chl(1), STAT=ierror )
373            IF( ierror > 0 ) THEN
374               CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' )   ;   RETURN
375            ENDIF
376            ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1)   )
377            IF( sn_chl%ln_tint )   ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) )
378            !                                        ! fill sf_chl with sn_chl and control print
379            CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',   &
380               &           'Solar penetration function of read chlorophyll', 'namtra_qsr' , no_print )
381         ENDIF
382         IF( nqsr == np_RGB ) THEN                 ! constant Chl
383            IF(lwp) WRITE(numout,*) '   ==>>>   Constant Chlorophyll concentration = 0.05'
384         ENDIF
385         !
386      CASE( np_2BD )                   !==  2 bands light penetration  ==!
387         !
388         IF(lwp)  WRITE(numout,*) '   ==>>>   2 bands light penetration'
389         !
390         nksr = trc_oce_ext_lev( rn_si1, 100._wp )    ! level of light extinction
391         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
392         !
393      CASE( np_BIO )                   !==  BIO light penetration  ==!
394         !
395         IF(lwp) WRITE(numout,*) '   ==>>>   bio-model light penetration'
396         IF( .NOT.lk_top )   CALL ctl_stop( 'No bio model : ln_qsr_bio = true impossible ' )
397         !
398      END SELECT
399      !
400      qsr_hc(:,:,:) = 0._wp     ! now qsr heat content set to zero where it will not be computed
401      !
402      ! 1st ocean level attenuation coefficient (used in sbcssm)
403      IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN
404         CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev'  , fraqsr_1lev, ldxios = lrxios  )
405      ELSE
406         fraqsr_1lev(:,:) = 1._wp   ! default : no penetration
407      ENDIF
408      !
409      IF( lwxios ) THEN
410         CALL iom_set_rstw_var_active('qsr_hc_b')
411         CALL iom_set_rstw_var_active('fraqsr_1lev')
412      ENDIF
413      !
414   END SUBROUTINE tra_qsr_init
415
416   !!======================================================================
417END MODULE traqsr