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 branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90 @ 6748

Last change on this file since 6748 was 6748, checked in by mocavero, 8 years ago

GYRE hybrid parallelization

  • Property svn:keywords set to Id
File size: 22.7 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 manager
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 wrk_nemo       ! Memory Allocation
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   LOGICAL , PUBLIC ::   ln_qsr_ice   !: light penetration for ice-model LIM3 (clem)
51   INTEGER , PUBLIC ::   nn_chldta    !: use Chlorophyll data (=1) or not (=0)
52   REAL(wp), PUBLIC ::   rn_abs       !: fraction absorbed in the very near surface (RGB & 2 bands)
53   REAL(wp), PUBLIC ::   rn_si0       !: very near surface depth of extinction      (RGB & 2 bands)
54   REAL(wp), PUBLIC ::   rn_si1       !: deepest depth of extinction (water type I)       (2 bands)
55   !
56   INTEGER , PUBLIC ::   nksr         !: levels below which the light cannot penetrate (depth larger than 391 m)
57 
58   INTEGER, PARAMETER ::   np_RGB  = 1   ! R-G-B     light penetration with constant Chlorophyll
59   INTEGER, PARAMETER ::   np_RGBc = 2   ! R-G-B     light penetration with Chlorophyll data
60   INTEGER, PARAMETER ::   np_2BD  = 3   ! 2 bands   light penetration
61   INTEGER, PARAMETER ::   np_BIO  = 4   ! bio-model light penetration
62   !
63   INTEGER  ::   nqsr    ! user choice of the type of light penetration
64   REAL(wp) ::   xsi0r   ! inverse of rn_si0
65   REAL(wp) ::   xsi1r   ! inverse of rn_si1
66   !
67   REAL(wp) , DIMENSION(3,61)           ::   rkrgb    ! tabulated attenuation coefficients for RGB absorption
68   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read)
69
70   !! * Substitutions
71#  include "vectopt_loop_substitute.h90"
72   !!----------------------------------------------------------------------
73   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
74   !! $Id$
75   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
76   !!----------------------------------------------------------------------
77CONTAINS
78
79   SUBROUTINE tra_qsr( kt )
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 ] / (rau0*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      !
107      INTEGER  ::   ji, jj, jk               ! dummy loop indices
108      INTEGER  ::   irgb                     ! local integers
109      REAL(wp) ::   zchl, zcoef, z1_2        ! local scalars
110      REAL(wp) ::   zc0 , zc1 , zc2 , zc3    !    -         -
111      REAL(wp) ::   zzc0, zzc1, zzc2, zzc3   !    -         -
112      REAL(wp) ::   zz0 , zz1                !    -         -
113      REAL(wp) ::   zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze
114      REAL(wp) ::   zlogc, zlogc2, zlogc3 
115      REAL(wp), POINTER, DIMENSION(:,:)   :: zekb, zekg, zekr
116      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt
117      REAL(wp), POINTER, DIMENSION(:,:,:) :: zetot, zchl3d
118      !!----------------------------------------------------------------------
119      !
120      IF( nn_timing == 1 )  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         CALL wrk_alloc( jpi,jpj,jpk,   ztrdt ) 
130!$OMP PARALLEL WORKSHARE
131         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
132!$OMP END PARALLEL WORKSHARE
133      ENDIF
134      !
135      !                         !-----------------------------------!
136      !                         !  before qsr induced heat content  !
137      !                         !-----------------------------------!
138      IF( kt == nit000 ) THEN          !==  1st time step  ==!
139!!gm case neuler  not taken into account....
140         IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN    ! read in restart
141            IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field read in the restart file'
142            z1_2 = 0.5_wp
143            CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b )   ! before heat content trend due to Qsr flux
144         ELSE                                           ! No restart or restart not found: Euler forward time stepping
145            z1_2 = 1._wp
146            qsr_hc_b(:,:,:) = 0._wp
147         ENDIF
148      ELSE                             !==  Swap of qsr heat content  ==!
149         z1_2 = 0.5_wp
150!$OMP PARALLEL WORKSHARE
151         qsr_hc_b(:,:,:) = qsr_hc(:,:,:)
152!$OMP END PARALLEL WORKSHARE
153      ENDIF
154      !
155      !                         !--------------------------------!
156      SELECT CASE( nqsr )       !  now qsr induced heat content  !
157      !                         !--------------------------------!
158      !
159      CASE( np_BIO )                   !==  bio-model fluxes  ==!
160         !
161!$OMP PARALLEL DO schedule(static) private(jk)
162         DO jk = 1, nksr
163            qsr_hc(:,:,jk) = r1_rau0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) )
164         END DO
165         !
166      CASE( np_RGB , np_RGBc )         !==  R-G-B fluxes  ==!
167         !
168         CALL wrk_alloc( jpi,jpj,       zekb, zekg, zekr        ) 
169         CALL wrk_alloc( jpi,jpj,jpk,   ze0, ze1, ze2, ze3, zea, zchl3d ) 
170         !
171         IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll
172            CALL fld_read( kt, 1, sf_chl )         ! Read Chl data and provides it at the current time step
173!$OMP PARALLEL DO schedule(static) private(jk,jj,ji,zchl,zCtot,zze,zpsi,zlogc,zlogc2,zlogc3,zCb,zCmax,zpsimax,zdelpsi,zCze)
174            DO jk = 1, nksr + 1
175               DO jj = 2, jpjm1                       ! Separation in R-G-B depending of the surface Chl
176                  DO ji = fs_2, fs_jpim1
177                     zchl    = sf_chl(1)%fnow(ji,jj,1)
178                     zCtot   = 40.6  * zchl**0.459
179                     zze     = 568.2 * zCtot**(-0.746)
180                     IF( zze > 102. ) zze = 200.0 * zCtot**(-0.293)
181                     zpsi    = gdepw_n(ji,jj,jk) / zze
182                     !
183                     zlogc   = LOG( zchl )
184                     zlogc2  = zlogc * zlogc
185                     zlogc3  = zlogc * zlogc * zlogc
186                     zCb     = 0.768 + 0.087 * zlogc - 0.179 * zlogc2 - 0.025 * zlogc3
187                     zCmax   = 0.299 - 0.289 * zlogc + 0.579 * zlogc2
188                     zpsimax = 0.6   - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3
189                     zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2
190                     zCze    = 1.12  * (zchl)**0.803 
191                     !
192                     zchl3d(ji,jj,jk) = zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) / zdelpsi )**2 ) )
193                  END DO
194                  !
195               END DO
196            END DO
197         ELSE                                !* constant chrlorophyll
198!$OMP PARALLEL DO schedule(static) private(jk)
199           DO jk = 1, nksr + 1
200              zchl3d(:,:,jk) = 0.05 
201            ENDDO
202         ENDIF
203         !
204         zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B
205         DO jj = 2, jpjm1
206            DO ji = fs_2, fs_jpim1
207               ze0(ji,jj,1) = rn_abs * qsr(ji,jj)
208               ze1(ji,jj,1) = zcoef  * qsr(ji,jj)
209               ze2(ji,jj,1) = zcoef  * qsr(ji,jj)
210               ze3(ji,jj,1) = zcoef  * qsr(ji,jj)
211               zea(ji,jj,1) =          qsr(ji,jj)
212            END DO
213         END DO
214         !
215!$OMP PARALLEL DO schedule(static) private(jk,jj,ji,zchl,irgb,zc0,zc1,zc2,zc3)
216         DO jk = 2, nksr+1                   !* interior equi-partition in R-G-B depending of vertical profile of Chl
217            DO jj = 2, jpjm1
218               DO ji = fs_2, fs_jpim1
219                  zchl = MIN( 10. , MAX( 0.03, zchl3d(ji,jj,jk) ) )
220                  irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
221                  zekb(ji,jj) = rkrgb(1,irgb)
222                  zekg(ji,jj) = rkrgb(2,irgb)
223                  zekr(ji,jj) = rkrgb(3,irgb)
224               END DO
225            END DO
226
227            DO jj = 2, jpjm1
228               DO ji = fs_2, fs_jpim1
229                  zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * xsi0r       )
230                  zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekb(ji,jj) )
231                  zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekg(ji,jj) )
232                  zc3 = ze3(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekr(ji,jj) )
233                  ze0(ji,jj,jk) = zc0
234                  ze1(ji,jj,jk) = zc1
235                  ze2(ji,jj,jk) = zc2
236                  ze3(ji,jj,jk) = zc3
237                  zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk)
238               END DO
239            END DO
240         END DO
241         !
242!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
243         DO jk = 1, nksr                     !* now qsr induced heat content
244            DO jj = 2, jpjm1
245               DO ji = fs_2, fs_jpim1
246                  qsr_hc(ji,jj,jk) = r1_rau0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) )
247               END DO
248            END DO
249         END DO
250         !
251         CALL wrk_dealloc( jpi,jpj,        zekb, zekg, zekr        ) 
252         CALL wrk_dealloc( jpi,jpj,jpk,   ze0, ze1, ze2, ze3, zea, zchl3d ) 
253         !
254      CASE( np_2BD  )            !==  2-bands fluxes  ==!
255         !
256         zz0 =        rn_abs   * r1_rau0_rcp      ! surface equi-partition in 2-bands
257         zz1 = ( 1. - rn_abs ) * r1_rau0_rcp
258!$OMP PARALLEL DO schedule(static) private(jk,jj,ji,zc0,zc1)
259         DO jk = 1, nksr                          ! solar heat absorbed at T-point in the top 400m
260            DO jj = 2, jpjm1
261               DO ji = fs_2, fs_jpim1
262                  zc0 = zz0 * EXP( -gdepw_n(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk  )*xsi1r )
263                  zc1 = zz0 * EXP( -gdepw_n(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk+1)*xsi1r )
264                  qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 
265               END DO
266            END DO
267         END DO
268         !
269      END SELECT
270      !
271      !                          !-----------------------------!
272!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
273      DO jk = 1, nksr            !  update to the temp. trend  !
274         DO jj = 2, jpjm1        !-----------------------------!
275            DO ji = fs_2, fs_jpim1   ! vector opt.
276               tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
277                  &                 + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t_n(ji,jj,jk)
278            END DO
279         END DO
280      END DO
281      !
282      IF( ln_qsr_ice ) THEN      ! sea-ice: store the 1st ocean level attenuation coefficient
283         DO jj = 2, jpjm1 
284            DO ji = fs_2, fs_jpim1   ! vector opt.
285               IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) )
286               ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp
287               ENDIF
288            END DO
289         END DO
290         ! Update haloes since lim_thd needs fraqsr_1lev to be defined everywhere
291         CALL lbc_lnk( fraqsr_1lev(:,:), 'T', 1._wp )
292      ENDIF
293      !
294      IF( iom_use('qsr3d') ) THEN      ! output the shortwave Radiation distribution
295         CALL wrk_alloc( jpi,jpj,jpk,   zetot )
296         !
297         zetot(:,:,nksr+1:jpk) = 0._wp     ! below ~400m set to zero
298         DO jk = nksr, 1, -1
299            zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) / r1_rau0_rcp
300         END DO         
301         CALL iom_put( 'qsr3d', zetot )   ! 3D distribution of shortwave Radiation
302         !
303         CALL wrk_dealloc( jpi,jpj,jpk,   zetot ) 
304      ENDIF
305      !
306      IF( lrst_oce ) THEN     ! write in the ocean restart file
307         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc      )
308         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev ) 
309      ENDIF
310      !
311      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
312         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
313         CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt )
314         CALL wrk_dealloc( jpi,jpj,jpk,   ztrdt ) 
315      ENDIF
316      !                       ! print mean trends (used for debugging)
317      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
318      !
319      IF( nn_timing == 1 )  CALL timing_stop('tra_qsr')
320      !
321   END SUBROUTINE tra_qsr
322
323
324   SUBROUTINE tra_qsr_init
325      !!----------------------------------------------------------------------
326      !!                  ***  ROUTINE tra_qsr_init  ***
327      !!
328      !! ** Purpose :   Initialization for the penetrative solar radiation
329      !!
330      !! ** Method  :   The profile of solar radiation within the ocean is set
331      !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
332      !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
333      !!      default values correspond to clear water (type I in Jerlov'
334      !!      (1968) classification.
335      !!         called by tra_qsr at the first timestep (nit000)
336      !!
337      !! ** Action  : - initialize rn_si0, rn_si1 and rn_abs
338      !!
339      !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
340      !!----------------------------------------------------------------------
341      INTEGER  ::   ji, jj, jk                  ! dummy loop indices
342      INTEGER  ::   ios, irgb, ierror, ioptio   ! local integer
343      REAL(wp) ::   zz0, zc0 , zc1, zcoef      ! local scalars
344      REAL(wp) ::   zz1, zc2 , zc3, zchl       !   -      -
345      !
346      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files
347      TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read
348      !!
349      NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice,  &
350         &                  nn_chldta, rn_abs, rn_si0, rn_si1
351      !!----------------------------------------------------------------------
352      !
353      IF( nn_timing == 1 )   CALL timing_start('tra_qsr_init')
354      !
355      REWIND( numnam_ref )              ! Namelist namtra_qsr in reference     namelist
356      READ  ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901)
357901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_qsr in reference namelist', lwp )
358      !
359      REWIND( numnam_cfg )              ! Namelist namtra_qsr in configuration namelist
360      READ  ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 )
361902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist', lwp )
362      IF(lwm) WRITE ( numond, namtra_qsr )
363      !
364      IF(lwp) THEN                ! control print
365         WRITE(numout,*)
366         WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
367         WRITE(numout,*) '~~~~~~~~~~~~'
368         WRITE(numout,*) '   Namelist namtra_qsr : set the parameter of penetration'
369         WRITE(numout,*) '      RGB (Red-Green-Blue) light penetration       ln_qsr_rgb = ', ln_qsr_rgb
370         WRITE(numout,*) '      2 band               light penetration       ln_qsr_2bd = ', ln_qsr_2bd
371         WRITE(numout,*) '      bio-model            light penetration       ln_qsr_bio = ', ln_qsr_bio
372         WRITE(numout,*) '      light penetration for ice-model (LIM3)       ln_qsr_ice = ', ln_qsr_ice
373         WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)        nn_chldta  = ', nn_chldta
374         WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs     = ', rn_abs
375         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0     = ', rn_si0
376         WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1     = ', rn_si1
377         WRITE(numout,*)
378      ENDIF
379      !
380      ioptio = 0                    ! Parameter control
381      IF( ln_qsr_rgb  )   ioptio = ioptio + 1
382      IF( ln_qsr_2bd  )   ioptio = ioptio + 1
383      IF( ln_qsr_bio  )   ioptio = ioptio + 1
384      !
385      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE type of light penetration in namelist namtra_qsr',  &
386         &                               ' 2 bands, 3 RGB bands or bio-model light penetration' )
387      !
388      IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr = np_RGB 
389      IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr = np_RGBc
390      IF( ln_qsr_2bd                      )   nqsr = np_2BD
391      IF( ln_qsr_bio                      )   nqsr = np_BIO
392      !
393      !                             ! Initialisation
394      xsi0r = 1._wp / rn_si0
395      xsi1r = 1._wp / rn_si1
396      !
397      SELECT CASE( nqsr )
398      !                               
399      CASE( np_RGB , np_RGBc )         !==  Red-Green-Blue light penetration  ==!
400         !                             
401         IF(lwp)   WRITE(numout,*) '   R-G-B   light penetration '
402         !
403         CALL trc_oce_rgb( rkrgb )                 ! tabulated attenuation coef.
404         !                                   
405         nksr = trc_oce_ext_lev( r_si2, 33._wp )   ! level of light extinction
406         !
407         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
408         !
409         IF( nqsr == np_RGBc ) THEN                ! Chl data : set sf_chl structure
410            IF(lwp) WRITE(numout,*) '        Chlorophyll read in a file'
411            ALLOCATE( sf_chl(1), STAT=ierror )
412            IF( ierror > 0 ) THEN
413               CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' )   ;   RETURN
414            ENDIF
415            ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1)   )
416            IF( sn_chl%ln_tint )   ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) )
417            !                                        ! fill sf_chl with sn_chl and control print
418            CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',   &
419               &           'Solar penetration function of read chlorophyll', 'namtra_qsr' )
420         ENDIF
421         IF( nqsr == np_RGB ) THEN                 ! constant Chl
422            IF(lwp) WRITE(numout,*) '        Constant Chlorophyll concentration = 0.05'
423         ENDIF
424         !
425      CASE( np_2BD )                   !==  2 bands light penetration  ==!
426         !
427         IF(lwp)  WRITE(numout,*) '   2 bands light penetration'
428         !
429         nksr = trc_oce_ext_lev( rn_si1, 100._wp )    ! level of light extinction
430         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
431         !
432      CASE( np_BIO )                   !==  BIO light penetration  ==!
433         !
434         IF(lwp) WRITE(numout,*) '   bio-model light penetration'
435         IF( .NOT.lk_qsr_bio )   CALL ctl_stop( 'No bio model : ln_qsr_bio = true impossible ' )
436         !
437      END SELECT
438      !
439!$OMP PARALLEL WORKSHARE
440      qsr_hc(:,:,:) = 0._wp     ! now qsr heat content set to zero where it will not be computed
441!$OMP END PARALLEL WORKSHARE
442      !
443      ! 1st ocean level attenuation coefficient (used in sbcssm)
444      IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN
445         CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev'  , fraqsr_1lev  )
446      ELSE
447         fraqsr_1lev(:,:) = 1._wp   ! default : no penetration
448      ENDIF
449      !
450      IF( nn_timing == 1 )   CALL timing_stop('tra_qsr_init')
451      !
452   END SUBROUTINE tra_qsr_init
453
454   !!======================================================================
455END MODULE traqsr
Note: See TracBrowser for help on using the repository browser.