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 @ 14990

Last change on this file since 14990 was 14990, checked in by techene, 4 years ago

#2605 : new optimisation of traqsr

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