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_r13296_HPC-07_mocavero_mpi3/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/TRA/traqsr.F90 @ 13630

Last change on this file since 13630 was 13630, checked in by mocavero, 4 years ago

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