source: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/TRA/traqsr.F90 @ 14618

Last change on this file since 14618 was 14618, checked in by techene, 9 months ago

#2506 RK3 work in progess : for 2 configurag

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