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_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA – NEMO

source: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/traqsr.F90 @ 14986

Last change on this file since 14986 was 14986, checked in by sparonuz, 3 years ago

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