2020WP/ENHANCE-10_acc_fix_traqsr: traqsr_minmem.F90

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