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/trunk/src/OCE/TRA – NEMO

source: NEMO/trunk/src/OCE/TRA/traqsr.F90 @ 14587

Last change on this file since 14587 was 14215, checked in by acc, 3 years ago

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