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 @ 13551

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

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