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_r13383_HPC-02_Daley_Tiling/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/traqsr.F90 @ 13518

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

Tiling for modules before tra_adv

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