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/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90 @ 2317

Last change on this file since 2317 was 2317, checked in by cetlod, 14 years ago

Improve the computation of the divergence of the downward solar irradiance in traqsr.F90, see ticket #726

  • change fsdepw & fse3t into fsdepw_0 & fse3t_0 in tra_qsr_init
  • modify tra_qsr to recompute systematically etot3 in vvl case
  • suppress the namelist parameter rn_si2 which is computed as the RGB longest depth of extinction in trc_oce.F90
  • Property svn:keywords set to Id
File size: 26.5 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   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   tra_qsr      : trend due to the solar radiation penetration
16   !!   tra_qsr_init : solar radiation penetration initialization
17   !!----------------------------------------------------------------------
18   USE oce             ! ocean dynamics and active tracers
19   USE dom_oce         ! ocean space and time domain
20   USE sbc_oce         ! surface boundary condition: ocean
21   USE trc_oce         ! share SMS/Ocean variables
22   USE trdmod_oce      ! ocean variables trends
23   USE trdtra          ! ocean active tracers trends
24   USE in_out_manager  ! I/O manager
25   USE phycst          ! physical constants
26   USE prtctl          ! Print control
27   USE iom             ! I/O manager
28   USE fldread         ! read input fields
29   USE restart         ! ocean restart
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   tra_qsr       ! routine called by step.F90 (ln_traqsr=T)
35   PUBLIC   tra_qsr_init  ! routine called by opa.F90
36
37   !                                           !!* Namelist namtra_qsr: penetrative solar radiation
38   LOGICAL , PUBLIC ::   ln_traqsr  = .TRUE.    !: light absorption (qsr) flag
39   LOGICAL , PUBLIC ::   ln_qsr_rgb = .FALSE.   !: Red-Green-Blue light absorption flag 
40   LOGICAL , PUBLIC ::   ln_qsr_2bd = .TRUE.    !: 2 band         light absorption flag
41   LOGICAL , PUBLIC ::   ln_qsr_bio = .FALSE.   !: bio-model      light absorption flag
42   INTEGER , PUBLIC ::   nn_chldta  = 0         !: use Chlorophyll data (=1) or not (=0)
43   REAL(wp), PUBLIC ::   rn_abs     = 0.58_wp   !: fraction absorbed in the very near surface (RGB & 2 bands)
44   REAL(wp), PUBLIC ::   rn_si0     = 0.35_wp   !: very near surface depth of extinction      (RGB & 2 bands)
45   REAL(wp), PUBLIC ::   rn_si1     = 23.0_wp   !: deepest depth of extinction (water type I)       (2 bands)
46   
47   ! Module variables
48   REAL(wp) ::   xsi0r                           !: inverse of rn_si0
49   REAL(wp) ::   xsi1r                           !: inverse of rn_si1
50!$AGRIF_DO_NOT_TREAT
51   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read)
52   INTEGER, PUBLIC ::   nksr              ! levels below which the light cannot penetrate ( depth larger than 391 m)
53   REAL(wp), DIMENSION(3,61) ::   rkrgb   !: tabulated attenuation coefficients for RGB absorption
54!$AGRIF_END_DO_NOT_TREAT
55
56   !! * Substitutions
57#  include "domzgr_substitute.h90"
58#  include "vectopt_loop_substitute.h90"
59   !!----------------------------------------------------------------------
60   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
61   !! $Id$
62   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
63   !!----------------------------------------------------------------------
64
65CONTAINS
66
67   SUBROUTINE tra_qsr( kt )
68      !!----------------------------------------------------------------------
69      !!                  ***  ROUTINE tra_qsr  ***
70      !!
71      !! ** Purpose :   Compute the temperature trend due to the solar radiation
72      !!      penetration and add it to the general temperature trend.
73      !!
74      !! ** Method  : The profile of the solar radiation within the ocean is defined
75      !!      through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) and a ratio rn_abs
76      !!      Considering the 2 wavebands case:
77      !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
78      !!         The temperature trend associated with the solar radiation penetration
79      !!         is given by : zta = 1/e3t dk[ I ] / (rau0*Cp)
80      !!         At the bottom, boudary condition for the radiation is no flux :
81      !!      all heat which has not been absorbed in the above levels is put
82      !!      in the last ocean level.
83      !!         In z-coordinate case, the computation is only done down to the
84      !!      level where I(k) < 1.e-15 W/m2. In addition, the coefficients
85      !!      used for the computation are calculated one for once as they
86      !!      depends on k only.
87      !!
88      !! ** Action  : - update ta with the penetrative solar radiation trend
89      !!              - save the trend in ttrd ('key_trdtra')
90      !!
91      !! Reference  : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
92      !!              Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
93      !!----------------------------------------------------------------------
94      !!
95      INTEGER, INTENT(in) ::   kt     ! ocean time-step
96      !!
97      INTEGER  ::   ji, jj, jk           ! dummy loop indices
98      INTEGER  ::   irgb                 ! temporary integers
99      REAL(wp) ::   zchl, zcoef          ! temporary scalars
100      REAL(wp) ::   zc0, zc1, zc2, zc3   !    -         -
101      REAL(wp) ::   zz0, zz1             !    -         -
102      REAL(wp) ::   z1_e3t, zfact        !    -         -
103      REAL(wp), DIMENSION(jpi,jpj)     ::   zekb, zekg, zekr            ! 2D workspace
104      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze0, ze1 , ze2, ze3, zea    ! 3D workspace
105      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ztrdt,  ztrds
106      !!----------------------------------------------------------------------
107
108      IF( kt == nit000 ) THEN
109         IF(lwp) WRITE(numout,*)
110         IF(lwp) WRITE(numout,*) 'tra_qsr : penetration of the surface solar radiation'
111         IF(lwp) WRITE(numout,*) '~~~~~~~'
112         IF( .NOT.ln_traqsr )   RETURN
113      ENDIF
114
115      IF( l_trdtra ) THEN      ! Save ta and sa trends
116         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
117         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = 0.
118      ENDIF
119
120      !                                        Set before qsr tracer content field
121      !                                        ***********************************
122      IF( kt == nit000 ) THEN                     ! Set the forcing field at nit000 - 1
123         !                                        ! -----------------------------------
124         IF( ln_rstart .AND.    &                    ! Restart: read in restart file
125              & iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN
126            IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field red in the restart file'
127            zfact = 0.5e0
128            CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b )   ! before heat content trend due to Qsr flux
129         ELSE                                           ! No restart or restart not found: Euler forward time stepping
130            zfact = 1.e0
131            qsr_hc_b(:,:,:) = 0.e0
132         ENDIF
133      ELSE                                        ! Swap of forcing field
134         !                                        ! ---------------------
135         zfact = 0.5e0
136         qsr_hc_b(:,:,:) = qsr_hc(:,:,:)
137      ENDIF
138      !                                        Compute now qsr tracer content field
139      !                                        ************************************
140     
141      !                                           ! ============================================== !
142      IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN      !  bio-model fluxes  : all vertical coordinates  !
143         !                                        ! ============================================== !
144         DO jk = 1, jpkm1
145            DO jj = 2, jpjm1
146               DO ji = fs_2, fs_jpim1   ! vector opt.
147                  qsr_hc(ji,jj,jk) =  ro0cpr * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) / fse3t(ji,jj,jk) 
148               END DO
149            END DO
150         END DO
151         CALL iom_put( 'qsr3d', etot3 )   ! Shortwave Radiation 3D distribution
152         !                                        ! ============================================== !
153      ELSE                                        !  Ocean alone :
154         !                                        ! ============================================== !
155         !
156         !                                                ! ------------------------- !
157         IF( ln_qsr_rgb) THEN                             !  R-G-B  light penetration !
158            !                                             ! ------------------------- !
159            ! Set chlorophyl concentration
160            IF( nn_chldta == 1 .OR. lk_vvl ) THEN            !*  Variable Chlorophyll or ocean volume
161               !
162               IF( nn_chldta == 1 ) THEN                             !*  Variable Chlorophyll
163                  !
164                  CALL fld_read( kt, 1, sf_chl )                         ! Read Chl data and provides it at the current time step
165                  !         
166!CDIR COLLAPSE
167!CDIR NOVERRCHK
168                  DO jj = 1, jpj                                         ! Separation in R-G-B depending of the surface Chl
169!CDIR NOVERRCHK
170                     DO ji = 1, jpi
171                        zchl = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) )
172                        irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
173                        zekb(ji,jj) = rkrgb(1,irgb)
174                        zekg(ji,jj) = rkrgb(2,irgb)
175                        zekr(ji,jj) = rkrgb(3,irgb)
176                     END DO
177                  END DO
178               ELSE                                            ! Variable ocean volume but constant chrlorophyll
179                  zchl = 0.05                                     ! constant chlorophyll
180                  irgb = NINT( 41 + 20.*LOG10( zchl ) + 1.e-15 )
181                  zekb(:,:) = rkrgb(1,irgb)                       ! Separation in R-G-B depending of the chlorophyll
182                  zekg(:,:) = rkrgb(2,irgb)
183                  zekr(:,:) = rkrgb(3,irgb)
184               ENDIF
185               !
186               zcoef  = ( 1. - rn_abs ) / 3.e0                        ! equi-partition in R-G-B
187               ze0(:,:,1) = rn_abs  * qsr(:,:)
188               ze1(:,:,1) = zcoef * qsr(:,:)
189               ze2(:,:,1) = zcoef * qsr(:,:)
190               ze3(:,:,1) = zcoef * qsr(:,:)
191               zea(:,:,1) =         qsr(:,:)
192               !
193               DO jk = 2, nksr+1
194!CDIR NOVERRCHK
195                  DO jj = 1, jpj
196!CDIR NOVERRCHK   
197                     DO ji = 1, jpi
198                        zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r     )
199                        zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekb(ji,jj) )
200                        zc2 = ze2(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekg(ji,jj) )
201                        zc3 = ze3(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekr(ji,jj) )
202                        ze0(ji,jj,jk) = zc0
203                        ze1(ji,jj,jk) = zc1
204                        ze2(ji,jj,jk) = zc2
205                        ze3(ji,jj,jk) = zc3
206                        zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk)
207                     END DO
208                  END DO
209               END DO
210               !
211               DO jk = 1, nksr                                        ! compute and add qsr trend to ta
212                  qsr_hc(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) )
213               END DO
214               zea(:,:,nksr+1:jpk) = 0.e0     ! below 400m set to zero
215               CALL iom_put( 'qsr3d', zea )   ! Shortwave Radiation 3D distribution
216               !
217            ELSE                                                 !*  Constant Chlorophyll
218               DO jk = 1, nksr
219                  qsr_hc(:,:,jk) =  etot3(:,:,jk) * qsr(:,:)
220               END DO
221            ENDIF
222
223         ENDIF
224         !                                                ! ------------------------- !
225         IF( ln_qsr_2bd ) THEN                            !  2 band light penetration !
226            !                                             ! ------------------------- !
227            !
228            IF( lk_vvl ) THEN                                  !* variable volume
229               zz0   =        rn_abs   * ro0cpr
230               zz1   = ( 1. - rn_abs ) * ro0cpr
231               DO jk = 1, nksr                    ! solar heat absorbed at T-point in the top 400m
232                  DO jj = 2, jpjm1
233                     DO ji = 2, jpim1
234                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
235                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
236                        qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
237                     END DO
238                  END DO
239               END DO
240            ELSE                                               !* constant volume: coef. computed one for all
241               DO jk = 1, nksr
242                  DO jj = 2, jpjm1
243                     DO ji = fs_2, fs_jpim1   ! vector opt.
244                        qsr_hc(ji,jj,jk) =  etot3(ji,jj,jk) * qsr(ji,jj)
245                     END DO
246                  END DO
247               END DO
248               !
249            ENDIF
250            !
251         ENDIF
252         !
253      ENDIF
254      !                                        Add to the general trend
255      !                                        ************************
256      DO jk = 1, nksr
257         DO jj = 2, jpjm1 
258            DO ji = fs_2, fs_jpim1   ! vector opt.
259               z1_e3t = zfact / fse3t(ji,jj,jk)
260               tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) * z1_e3t
261            END DO
262         END DO
263      END DO
264      !
265      IF( lrst_oce ) THEN   !                  Write in the ocean restart file
266         !                                     *******************************
267         IF(lwp) WRITE(numout,*)
268         IF(lwp) WRITE(numout,*) 'qsr tracer content forcing field written in ocean restart file ',   &
269            &                    'at it= ', kt,' date= ', ndastp
270         IF(lwp) WRITE(numout,*) '~~~~'
271         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b', qsr_hc )
272         !
273      ENDIF
274
275      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
276         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
277         CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_qsr, ztrdt )
278         CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_qsr, ztrds )
279         DEALLOCATE( ztrdt )    ;        DEALLOCATE( ztrds )
280      ENDIF
281      !                       ! print mean trends (used for debugging)
282      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
283      !
284   END SUBROUTINE tra_qsr
285
286
287   SUBROUTINE tra_qsr_init
288      !!----------------------------------------------------------------------
289      !!                  ***  ROUTINE tra_qsr_init  ***
290      !!
291      !! ** Purpose :   Initialization for the penetrative solar radiation
292      !!
293      !! ** Method  :   The profile of solar radiation within the ocean is set
294      !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
295      !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
296      !!      default values correspond to clear water (type I in Jerlov'
297      !!      (1968) classification.
298      !!         called by tra_qsr at the first timestep (nit000)
299      !!
300      !! ** Action  : - initialize rn_si0, rn_si1 and rn_abs
301      !!
302      !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
303      !!----------------------------------------------------------------------
304      INTEGER  ::   ji, jj, jk            ! dummy loop indices
305      INTEGER  ::   irgb, ierror          ! temporary integer
306      INTEGER  ::   ioptio, nqsr          ! temporary integer
307      REAL(wp) ::   zc0  , zc1, zcoef     ! temporary scalars
308      REAL(wp) ::   zc2  , zc3  , zchl    !    -         -
309      REAL(wp) ::   zz0  , zz1            !    -         -
310      REAL(wp), DIMENSION(jpi,jpj)     ::   zekb, zekg, zekr              ! 2D workspace
311      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze0 , ze1 , ze2 , ze3 , zea   ! 3D workspace
312      !!
313      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files
314      TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read
315      NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio,   &
316         &                  nn_chldta, rn_abs, rn_si0, rn_si1
317      !!----------------------------------------------------------------------
318
319      cn_dir = './'       ! directory in which the model is executed
320      ! ... default values (NB: frequency positive => hours, negative => months)
321      !            !     file       ! frequency !  variable  ! time interp !  clim   ! 'yearly' or ! weights  ! rotation   !
322      !            !     name       !  (hours)  !    name    !    (T/F)    !  (T/F)  ! 'monthly'   ! filename ! pairs      !
323      sn_chl = FLD_N( 'chlorophyll' ,    -1     ,  'CHLA'    ,  .true.     , .true.  ,   'yearly'  , ''       , ''         )
324      !
325      REWIND( numnam )            ! Read Namelist namtra_qsr : ratio and length of penetration
326      READ  ( numnam, namtra_qsr )
327      !
328      IF(lwp) THEN                ! control print
329         WRITE(numout,*)
330         WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
331         WRITE(numout,*) '~~~~~~~~~~~~'
332         WRITE(numout,*) '   Namelist namtra_qsr : set the parameter of penetration'
333         WRITE(numout,*) '      Light penetration (T) or not (F)         ln_traqsr  = ', ln_traqsr
334         WRITE(numout,*) '      RGB (Red-Green-Blue) light penetration   ln_qsr_rgb = ', ln_qsr_rgb
335         WRITE(numout,*) '      2 band               light penetration   ln_qsr_2bd = ', ln_qsr_2bd
336         WRITE(numout,*) '      bio-model            light penetration   ln_qsr_bio = ', ln_qsr_bio
337         WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)    nn_chldta  = ', nn_chldta
338         WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs = ', rn_abs
339         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0 = ', rn_si0
340         WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1 = ', rn_si1
341      ENDIF
342
343      IF( ln_traqsr ) THEN     ! control consistency
344         !                     
345         IF( .NOT.lk_qsr_bio .AND. ln_qsr_bio )   THEN
346            CALL ctl_warn( 'No bio model : force ln_qsr_bio = FALSE ' )
347            ln_qsr_bio = .FALSE.
348         ENDIF
349         !
350         ioptio = 0                      ! Parameter control
351         IF( ln_qsr_rgb  )   ioptio = ioptio + 1
352         IF( ln_qsr_2bd  )   ioptio = ioptio + 1
353         IF( ln_qsr_bio  )   ioptio = ioptio + 1
354         !
355         IF( ioptio /= 1 ) &
356            CALL ctl_stop( '          Choose ONE type of light penetration in namelist namtra_qsr',  &
357            &              ' 2 bands, 3 RGB bands or bio-model light penetration' )
358         !
359         IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr =  1 
360         IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr =  2
361         IF( ln_qsr_2bd                      )   nqsr =  3
362         IF( ln_qsr_bio                      )   nqsr =  4
363         !
364         IF(lwp) THEN                   ! Print the choice
365            WRITE(numout,*)
366            IF( nqsr ==  1 )   WRITE(numout,*) '         R-G-B   light penetration - Constant Chlorophyll'
367            IF( nqsr ==  2 )   WRITE(numout,*) '         R-G-B   light penetration - Chl data '
368            IF( nqsr ==  3 )   WRITE(numout,*) '         2 bands light penetration'
369            IF( nqsr ==  4 )   WRITE(numout,*) '         bio-model light penetration'
370         ENDIF
371         !
372      ENDIF
373      !                          ! ===================================== !
374      IF( ln_traqsr  ) THEN      !  Initialisation of Light Penetration  ! 
375         !                       ! ===================================== !
376         !
377         xsi0r = 1.e0 / rn_si0
378         xsi1r = 1.e0 / rn_si1
379         !                                ! ---------------------------------- !
380         IF( ln_qsr_rgb ) THEN            !  Red-Green-Blue light penetration  !
381            !                             ! ---------------------------------- !
382            !
383            CALL trc_oce_rgb( rkrgb )           !* tabulated attenuation coef.
384            !
385            !                                   !* level of light extinction
386            IF(  ln_sco ) THEN   ;   nksr = jpkm1
387            ELSE                 ;   nksr = trc_oce_ext_lev( r_si2, 0.33e2 )
388            ENDIF
389
390            IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_0(nksr+1), ' m'
391            !
392            IF( nn_chldta == 1 ) THEN           !* Chl data : set sf_chl structure
393               IF(lwp) WRITE(numout,*)
394               IF(lwp) WRITE(numout,*) '        Chlorophyll read in a file'
395               ALLOCATE( sf_chl(1), STAT=ierror )
396               IF( ierror > 0 ) THEN
397                  CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' )   ;   RETURN
398               ENDIF
399               ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1)   )
400               IF( sn_chl%ln_tint )ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) )
401               !                                        ! fill sf_chl with sn_chl and control print
402               CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',   &
403                  &                                         'Solar penetration function of read chlorophyll', 'namtra_qsr' )
404               !
405            ELSE                                !* constant Chl : compute once for all the distribution of light (etot3)
406               IF(lwp) WRITE(numout,*)
407               IF(lwp) WRITE(numout,*) '        Constant Chlorophyll concentration = 0.05'
408               IF( lk_vvl ) THEN                   ! variable volume
409                  IF(lwp) WRITE(numout,*) '        key_vvl: light distribution will be computed at each time step'
410               ELSE                                ! constant volume: computes one for all
411                  IF(lwp) WRITE(numout,*) '        fixed volume: light distribution computed one for all'
412                  !
413                  zchl = 0.05                                 ! constant chlorophyll
414                  irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
415                  zekb(:,:) = rkrgb(1,irgb)                   ! Separation in R-G-B depending of the chlorophyll
416                  zekg(:,:) = rkrgb(2,irgb)
417                  zekr(:,:) = rkrgb(3,irgb)
418                  !
419                  zcoef = ( 1. - rn_abs ) / 3.e0              ! equi-partition in R-G-B
420                  ze0(:,:,1) = rn_abs
421                  ze1(:,:,1) = zcoef
422                  ze2(:,:,1) = zcoef 
423                  ze3(:,:,1) = zcoef
424                  zea(:,:,1) = tmask(:,:,1)                   ! = ( ze0+ze1+z2+ze3 ) * tmask
425               
426                  DO jk = 2, nksr+1
427!CDIR NOVERRCHK
428                     DO jj = 1, jpj
429!CDIR NOVERRCHK   
430                        DO ji = 1, jpi
431                           zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t_0(ji,jj,jk-1) * xsi0r     )
432                           zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t_0(ji,jj,jk-1) * zekb(ji,jj) )
433                           zc2 = ze2(ji,jj,jk-1) * EXP( - fse3t_0(ji,jj,jk-1) * zekg(ji,jj) )
434                           zc3 = ze3(ji,jj,jk-1) * EXP( - fse3t_0(ji,jj,jk-1) * zekr(ji,jj) )
435                           ze0(ji,jj,jk) = zc0
436                           ze1(ji,jj,jk) = zc1
437                           ze2(ji,jj,jk) = zc2
438                           ze3(ji,jj,jk) = zc3
439                           zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk)
440                        END DO
441                     END DO
442                  END DO 
443                  !
444                  DO jk = 1, nksr
445                     etot3(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t_0(:,:,jk)
446                  END DO
447                  etot3(:,:,nksr+1:jpk) = 0.e0                ! below 400m set to zero
448               ENDIF
449            ENDIF
450            !
451         ENDIF
452            !                             ! ---------------------------------- !
453         IF( ln_qsr_2bd ) THEN            !    2 bands    light penetration    !
454            !                             ! ---------------------------------- !
455            !
456            !                                ! level of light extinction
457            nksr = trc_oce_ext_lev( rn_si1, 1.e2 )
458            IF(lwp) THEN
459               WRITE(numout,*)
460            IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_0(nksr+1), ' m'
461            ENDIF
462            !
463            IF( lk_vvl ) THEN                   ! variable volume
464               IF(lwp) WRITE(numout,*) '        key_vvl: light distribution will be computed at each time step'
465            ELSE                                ! constant volume: computes one for all
466               zz0 =        rn_abs   * ro0cpr
467               zz1 = ( 1. - rn_abs ) * ro0cpr
468               DO jk = 1, nksr                    !*  solar heat absorbed at T-point computed once for all
469                  DO jj = 1, jpj                              ! top 400 meters
470                     DO ji = 1, jpi
471                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
472                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
473                        etot3(ji,jj,jk) = (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) / fse3t_0(ji,jj,jk)
474                     END DO
475                  END DO
476               END DO
477               etot3(:,:,nksr+1:jpk) = 0.e0                   ! below 400m set to zero
478               !
479            ENDIF
480         ENDIF
481         !                       ! ===================================== !
482      ELSE                       !        No light penetration           !                   
483         !                       ! ===================================== !
484         IF(lwp) THEN
485            WRITE(numout,*)
486            WRITE(numout,*) 'tra_qsr_init : NO solar flux penetration'
487            WRITE(numout,*) '~~~~~~~~~~~~'
488         ENDIF
489      ENDIF
490      !
491   END SUBROUTINE tra_qsr_init
492
493   !!======================================================================
494END MODULE traqsr
Note: See TracBrowser for help on using the repository browser.