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

source: NEMO/branches/2020/dev_r13898_Tiling_Cleanup_MPI3/src/OCE/TRA/traqsr.F90 @ 13906

Last change on this file since 13906 was 13898, checked in by hadcv, 4 years ago

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