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/2020/dev_r12953_ENHANCE-10_acc_fix_traqsr/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r12953_ENHANCE-10_acc_fix_traqsr/src/OCE/TRA/traqsr.F90 @ 12956

Last change on this file since 12956 was 12956, checked in by acc, 4 years ago

Branch: 2020/dev_r12953_ENHANCE-10_acc_fix_traqsr. Optimised version of traqsr.F90 as discussed in wiki:2020WP/ENHANCE-10_acc_fix_traqsr (#2469). This version uses less memory than the original and achieves the same results faster. Successfully SETTE tested.
  • Property svn:keywords set to Id
File size: 22.3 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   !!----------------------------------------------------------------------
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$
73   !! Software governed by the CeCILL license (see ./LICENSE)
74   !!----------------------------------------------------------------------
75CONTAINS
76
77   SUBROUTINE tra_qsr( kt, Kmm, pts, Krhs )
78      !!----------------------------------------------------------------------
79      !!                  ***  ROUTINE tra_qsr  ***
80      !!
81      !! ** Purpose :   Compute the temperature trend due to the solar radiation
82      !!              penetration and add it to the general temperature trend.
83      !!
84      !! ** Method  : The profile of the solar radiation within the ocean is defined
85      !!      through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) and a ratio rn_abs
86      !!      Considering the 2 wavebands case:
87      !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
88      !!         The temperature trend associated with the solar radiation penetration
89      !!         is given by : zta = 1/e3t dk[ I ] / (rho0*Cp)
90      !!         At the bottom, boudary condition for the radiation is no flux :
91      !!      all heat which has not been absorbed in the above levels is put
92      !!      in the last ocean level.
93      !!         The computation is only done down to the level where
94      !!      I(k) < 1.e-15 W/m2 (i.e. over the top nksr levels) .
95      !!
96      !! ** Action  : - update ta with the penetrative solar radiation trend
97      !!              - send  trend for further diagnostics (l_trdtra=T)
98      !!
99      !! Reference  : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
100      !!              Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
101      !!              Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562
102      !!----------------------------------------------------------------------
103      INTEGER,                                   INTENT(in   ) :: kt            ! ocean time-step
104      INTEGER,                                   INTENT(in   ) :: Kmm, Krhs     ! time level indices
105      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts           ! active tracers and RHS of tracer equation
106      !
107      INTEGER  ::   ji, jj, jk               ! dummy loop indices
108      INTEGER  ::   irgb                     ! local integers
109      REAL(wp) ::   zchl, zcoef, z1_2        ! local scalars
110      REAL(wp) ::   zc0 , zc1 , zc2 , zc3    !    -         -
111      REAL(wp) ::   zzc0, zzc1, zzc2, zzc3   !    -         -
112      REAL(wp) ::   zz0 , zz1 , ze3t, zlui   !    -         -
113      REAL(wp) ::   zCb, zCmax, zpsi, zpsimax, zrdpsi, zCze
114      REAL(wp) ::   zlogc, zlogze, zlogCtot, zlogCze
115      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   :: ze0, ze1, ze2, ze3
116      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, zetot, 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( ze0 (jpi,jpj)           , ze1 (jpi,jpj) ,   &
162            &      ze2 (jpi,jpj)           , ze3 (jpi,jpj) ,   &
163            &      ztmp3d(jpi,jpj,nksr + 1)                     )
164         !
165         IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll
166            CALL fld_read( kt, 1, sf_chl )         ! Read Chl data and provides it at the current time step
167            !
168            ! Separation in R-G-B depending on the surface Chl
169            ! perform and store as many of the 2D calculations as possible
170            ! before the 3D loop (use the temporary 2D arrays to replace the
171            ! most expensive calculations)
172            !
173            ze0(:,:) = LOG ( MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(:,:,1) ) ) )     ! ze0 = log(zchl)
174            ze1(:,:) = 0.113328685307 + 0.803 * ze0(:,:)                           ! ze1 : log(zCze)  = log (1.12  * zchl**0.803)
175            ze2(:,:) = 3.703768066608 + 0.459 * ze0(:,:)                           ! ze2 : log(zCtot) = log(40.6  * zchl**0.459)
176            ze3(:,:) = 6.34247346942  - 0.746 * ze2(:,:)                           ! ze3 : log(zze)   = log(568.2 * zCtot**(-0.746))
177            WHERE( ze3 > 4.62497281328 ) ze3 = 5.298317366548 - 0.293 * ze2        !       IF( log(zze) > log(102.) ) log(zze) =
178            !                                                                      !                    log(200.0 * zCtot**(-0.293))
179            ze1(:,:) = EXP(ze1(:,:))                                               ! ze1 = zCze
180            ze2(:,:) = 1._wp / ( 0.710 + ze0(:,:) * ( 0.159 + ze0(:,:) * 0.021 ) ) ! ze2 = 1/zdelpsi
181            ze3(:,:) = 1._wp / EXP(ze3(:,:))                                       ! ze3 = 1/zze
182            !
183            DO_3D_00_00 ( 1, nksr + 1 )
184               ! zchl    = ALOG( ze0(ji,jj) )
185               zlogc = ze0(ji,jj)
186               !
187               zCb       = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) )
188               zCmax     = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 )
189               zpsimax   = 0.6   - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) )
190               ! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 )
191               !
192               zCze   = ze1(ji,jj)
193               zrdpsi = ze2(ji,jj)                                                 ! 1/zdelpsi
194               zpsi   = ze3(ji,jj) * gdepw(ji,jj,jk,Kmm)                           ! gdepw/zze
195               !
196               ! NB. make sure zchl value is such that: zchl = MIN( 10. , MAX( 0.03, zchl ) )
197               zchl = MIN( 10. , MAX( 0.03, zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) ) )
198               ! Convert chlorophyll value to attenuation coefficient look-up table index
199               ztmp3d(ji,jj,jk) = 41 + 20.*LOG10(zchl) + 1.e-15
200            END_3D
201         ELSE                                !* constant chlorophyll
202            zchl = 0.05
203            ! NB. make sure constant value is such that:
204            zchl = MIN( 10. , MAX( 0.03, zchl ) )
205            ! Convert chlorophyll value to attenuation coefficient look-up table index
206            zlui = 41 + 20.*LOG10(zchl) + 1.e-15
207            DO jk = 1, nksr + 1
208               ztmp3d(:,:,jk) = zlui 
209            END DO
210         ENDIF
211         !
212         zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B
213         DO_2D_00_00
214            ze0(ji,jj) = rn_abs * qsr(ji,jj)
215            ze1(ji,jj) = zcoef  * qsr(ji,jj)
216            ze2(ji,jj) = zcoef  * qsr(ji,jj)
217            ze3(ji,jj) = zcoef  * qsr(ji,jj)
218            ! store the surface SW radiation; re-use the surface ztmp3d array
219            ! since the surface attenuation coefficient is not used
220            ztmp3d(ji,jj,1) =       qsr(ji,jj)
221         END_2D
222         !
223         !* interior equi-partition in R-G-B depending on vertical profile of Chl
224         DO_3D_00_00 ( 2, nksr + 1 )
225            ze3t = e3t(ji,jj,jk-1,Kmm)
226            irgb = NINT( ztmp3d(ji,jj,jk) )
227            zc0 = ze0(ji,jj) * EXP( - ze3t * xsi0r )
228            zc1 = ze1(ji,jj) * EXP( - ze3t * rkrgb(1,irgb) )
229            zc2 = ze2(ji,jj) * EXP( - ze3t * rkrgb(2,irgb) )
230            zc3 = ze3(ji,jj) * EXP( - ze3t * rkrgb(3,irgb) )
231            ze0(ji,jj) = zc0
232            ze1(ji,jj) = zc1
233            ze2(ji,jj) = zc2
234            ze3(ji,jj) = zc3
235            ztmp3d(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk)
236         END_3D
237         !
238         DO_3D_00_00( 1, nksr )
239            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( ztmp3d(ji,jj,jk) - ztmp3d(ji,jj,jk+1) )
240         END_3D
241         !
242         DEALLOCATE( ze0 , ze1 , ze2 , ze3 , ztmp3d ) 
243         !
244      CASE( np_2BD  )            !==  2-bands fluxes  ==!
245         !
246         zz0 =        rn_abs   * r1_rho0_rcp      ! surface equi-partition in 2-bands
247         zz1 = ( 1. - rn_abs ) * r1_rho0_rcp
248         DO_3D_00_00( 1, nksr )
249            zc0 = zz0 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi1r )
250            zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r )
251            qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 
252         END_3D
253         !
254      END SELECT
255      !
256      !                          !-----------------------------!
257      DO_3D_00_00( 1, nksr )
258         pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
259            &                      + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t(ji,jj,jk,Kmm)
260      END_3D
261      !
262      ! sea-ice: store the 1st ocean level attenuation coefficient
263      DO_2D_00_00
264         IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) )
265         ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp
266         ENDIF
267      END_2D
268      CALL lbc_lnk( 'traqsr', fraqsr_1lev(:,:), 'T', 1._wp )
269      !
270      IF( iom_use('qsr3d') ) THEN      ! output the shortwave Radiation distribution
271         ALLOCATE( zetot(jpi,jpj,jpk) )
272         zetot(:,:,nksr+1:jpk) = 0._wp     ! below ~400m set to zero
273         DO jk = nksr, 1, -1
274            zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rho0_rcp
275         END DO         
276         CALL iom_put( 'qsr3d', zetot )   ! 3D distribution of shortwave Radiation
277         DEALLOCATE( zetot ) 
278      ENDIF
279      !
280      IF( lrst_oce ) THEN     ! write in the ocean restart file
281         IF( lwxios ) CALL iom_swap(      cwxios_context          )
282         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc     , ldxios = lwxios )
283         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev, ldxios = lwxios ) 
284         IF( lwxios ) CALL iom_swap(      cxios_context          )
285      ENDIF
286      !
287      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
288         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
289         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt )
290         DEALLOCATE( ztrdt ) 
291      ENDIF
292      !                       ! print mean trends (used for debugging)
293      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
294      !
295      IF( ln_timing )   CALL timing_stop('tra_qsr')
296      !
297   END SUBROUTINE tra_qsr
298
299
300   SUBROUTINE tra_qsr_init
301      !!----------------------------------------------------------------------
302      !!                  ***  ROUTINE tra_qsr_init  ***
303      !!
304      !! ** Purpose :   Initialization for the penetrative solar radiation
305      !!
306      !! ** Method  :   The profile of solar radiation within the ocean is set
307      !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
308      !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
309      !!      default values correspond to clear water (type I in Jerlov'
310      !!      (1968) classification.
311      !!         called by tra_qsr at the first timestep (nit000)
312      !!
313      !! ** Action  : - initialize rn_si0, rn_si1 and rn_abs
314      !!
315      !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
316      !!----------------------------------------------------------------------
317      INTEGER  ::   ji, jj, jk                  ! dummy loop indices
318      INTEGER  ::   ios, irgb, ierror, ioptio   ! local integer
319      REAL(wp) ::   zz0, zc0 , zc1, zcoef      ! local scalars
320      REAL(wp) ::   zz1, zc2 , zc3, zchl       !   -      -
321      !
322      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files
323      TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read
324      !!
325      NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio,  &
326         &                  nn_chldta, rn_abs, rn_si0, rn_si1
327      !!----------------------------------------------------------------------
328      !
329      READ  ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901)
330901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_qsr in reference namelist' )
331      !
332      READ  ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 )
333902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist' )
334      IF(lwm) WRITE ( numond, namtra_qsr )
335      !
336      IF(lwp) THEN                ! control print
337         WRITE(numout,*)
338         WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
339         WRITE(numout,*) '~~~~~~~~~~~~'
340         WRITE(numout,*) '   Namelist namtra_qsr : set the parameter of penetration'
341         WRITE(numout,*) '      RGB (Red-Green-Blue) light penetration       ln_qsr_rgb = ', ln_qsr_rgb
342         WRITE(numout,*) '      2 band               light penetration       ln_qsr_2bd = ', ln_qsr_2bd
343         WRITE(numout,*) '      bio-model            light penetration       ln_qsr_bio = ', ln_qsr_bio
344         WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)        nn_chldta  = ', nn_chldta
345         WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs     = ', rn_abs
346         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0     = ', rn_si0
347         WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1     = ', rn_si1
348         WRITE(numout,*)
349      ENDIF
350      !
351      ioptio = 0                    ! Parameter control
352      IF( ln_qsr_rgb  )   ioptio = ioptio + 1
353      IF( ln_qsr_2bd  )   ioptio = ioptio + 1
354      IF( ln_qsr_bio  )   ioptio = ioptio + 1
355      !
356      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE type of light penetration in namelist namtra_qsr',  &
357         &                               ' 2 bands, 3 RGB bands or bio-model light penetration' )
358      !
359      IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr = np_RGB 
360      IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr = np_RGBc
361      IF( ln_qsr_2bd                      )   nqsr = np_2BD
362      IF( ln_qsr_bio                      )   nqsr = np_BIO
363      !
364      !                             ! Initialisation
365      xsi0r = 1._wp / rn_si0
366      xsi1r = 1._wp / rn_si1
367      !
368      SELECT CASE( nqsr )
369      !                               
370      CASE( np_RGB , np_RGBc )         !==  Red-Green-Blue light penetration  ==!
371         !                             
372         IF(lwp)   WRITE(numout,*) '   ==>>>   R-G-B   light penetration '
373         !
374         CALL trc_oce_rgb( rkrgb )                 ! tabulated attenuation coef.
375         !                                   
376         nksr = trc_oce_ext_lev( r_si2, 33._wp )   ! level of light extinction
377         !
378         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
379         !
380         IF( nqsr == np_RGBc ) THEN                ! Chl data : set sf_chl structure
381            IF(lwp) WRITE(numout,*) '   ==>>>   Chlorophyll read in a file'
382            ALLOCATE( sf_chl(1), STAT=ierror )
383            IF( ierror > 0 ) THEN
384               CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' )   ;   RETURN
385            ENDIF
386            ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1)   )
387            IF( sn_chl%ln_tint )   ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) )
388            !                                        ! fill sf_chl with sn_chl and control print
389            CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',   &
390               &           'Solar penetration function of read chlorophyll', 'namtra_qsr' , no_print )
391         ENDIF
392         IF( nqsr == np_RGB ) THEN                 ! constant Chl
393            IF(lwp) WRITE(numout,*) '   ==>>>   Constant Chlorophyll concentration = 0.05'
394         ENDIF
395         !
396      CASE( np_2BD )                   !==  2 bands light penetration  ==!
397         !
398         IF(lwp)  WRITE(numout,*) '   ==>>>   2 bands light penetration'
399         !
400         nksr = trc_oce_ext_lev( rn_si1, 100._wp )    ! level of light extinction
401         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
402         !
403      CASE( np_BIO )                   !==  BIO light penetration  ==!
404         !
405         IF(lwp) WRITE(numout,*) '   ==>>>   bio-model light penetration'
406         IF( .NOT.lk_top )   CALL ctl_stop( 'No bio model : ln_qsr_bio = true impossible ' )
407         !
408      END SELECT
409      !
410      qsr_hc(:,:,:) = 0._wp     ! now qsr heat content set to zero where it will not be computed
411      !
412      ! 1st ocean level attenuation coefficient (used in sbcssm)
413      IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN
414         CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev'  , fraqsr_1lev, ldxios = lrxios  )
415      ELSE
416         fraqsr_1lev(:,:) = 1._wp   ! default : no penetration
417      ENDIF
418      !
419      IF( lwxios ) THEN
420         CALL iom_set_rstw_var_active('qsr_hc_b')
421         CALL iom_set_rstw_var_active('fraqsr_1lev')
422      ENDIF
423      !
424   END SUBROUTINE tra_qsr_init
425
426   !!======================================================================
427END MODULE traqsr
Note: See TracBrowser for help on using the repository browser.