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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90 @ 4431

Last change on this file since 4431 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
File size: 30.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   !!----------------------------------------------------------------------
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   USE lib_mpp         ! MPP library
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   tra_qsr       ! routine called by step.F90 (ln_traqsr=T)
36   PUBLIC   tra_qsr_init  ! routine called by opa.F90
37
38   !                                           !!* Namelist namtra_qsr: penetrative solar radiation
39   LOGICAL , PUBLIC ::   ln_traqsr  = .TRUE.    !: light absorption (qsr) flag
40   LOGICAL , PUBLIC ::   ln_qsr_rgb = .FALSE.   !: Red-Green-Blue light absorption flag 
41   LOGICAL , PUBLIC ::   ln_qsr_2bd = .TRUE.    !: 2 band         light absorption flag
42   LOGICAL , PUBLIC ::   ln_qsr_bio = .FALSE.   !: bio-model      light absorption flag
43   INTEGER , PUBLIC ::   nn_chldta  = 0         !: use Chlorophyll data (=1) or not (=0)
44   REAL(wp), PUBLIC ::   rn_abs     = 0.58_wp   !: fraction absorbed in the very near surface (RGB & 2 bands)
45   REAL(wp), PUBLIC ::   rn_si0     = 0.35_wp   !: very near surface depth of extinction      (RGB & 2 bands)
46   REAL(wp), PUBLIC ::   rn_si1     = 23.0_wp   !: deepest depth of extinction (water type I)       (2 bands)
47   
48   ! Module variables
49   REAL(wp) ::   xsi0r                           !: inverse of rn_si0
50   REAL(wp) ::   xsi1r                           !: inverse of rn_si1
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
55   !! * Control permutation of array indices
56#  include "oce_ftrans.h90"
57#  include "dom_oce_ftrans.h90"
58#  include "sbc_oce_ftrans.h90"
59#  include "trc_oce_ftrans.h90"
60
61   !! * Substitutions
62#  include "domzgr_substitute.h90"
63#  include "vectopt_loop_substitute.h90"
64   !!----------------------------------------------------------------------
65   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
66   !! $Id$
67   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
68   !!----------------------------------------------------------------------
69CONTAINS
70
71   SUBROUTINE tra_qsr( kt )
72      !!----------------------------------------------------------------------
73      !!                  ***  ROUTINE tra_qsr  ***
74      !!
75      !! ** Purpose :   Compute the temperature trend due to the solar radiation
76      !!      penetration and add it to the general temperature trend.
77      !!
78      !! ** Method  : The profile of the solar radiation within the ocean is defined
79      !!      through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) and a ratio rn_abs
80      !!      Considering the 2 wavebands case:
81      !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
82      !!         The temperature trend associated with the solar radiation penetration
83      !!         is given by : zta = 1/e3t dk[ I ] / (rau0*Cp)
84      !!         At the bottom, boudary condition for the radiation is no flux :
85      !!      all heat which has not been absorbed in the above levels is put
86      !!      in the last ocean level.
87      !!         In z-coordinate case, the computation is only done down to the
88      !!      level where I(k) < 1.e-15 W/m2. In addition, the coefficients
89      !!      used for the computation are calculated one for once as they
90      !!      depends on k only.
91      !!
92      !! ** Action  : - update ta with the penetrative solar radiation trend
93      !!              - save the trend in ttrd ('key_trdtra')
94      !!
95      !! Reference  : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
96      !!              Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
97      !!----------------------------------------------------------------------
98      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
99      USE wrk_nemo, ONLY:   zekb => wrk_2d_1 , zekg => wrk_2d_2 , zekr => wrk_2d_3
100      USE wrk_nemo, ONLY:   ze0  => wrk_3d_1 , ze1  => wrk_3d_2 , ze2  => wrk_3d_3
101      USE wrk_nemo, ONLY:   ze3  => wrk_3d_4 , zea  => wrk_3d_5
102
103      !! DCSE_NEMO: need additional directives for renamed module variables
104!FTRANS ze0 ze1 ze2 ze3 zea :I :I :z
105      !
106      INTEGER, INTENT(in) ::   kt     ! ocean time-step
107      !
108      INTEGER  ::   ji, jj, jk           ! dummy loop indices
109      INTEGER  ::   irgb                 ! local integers
110      REAL(wp) ::   zchl, zcoef, zfact   ! local scalars
111      REAL(wp) ::   zc0, zc1, zc2, zc3   !    -         -
112      REAL(wp) ::   zz0, zz1, z1_e3t     !    -         -
113
114!FTRANS ztrdt :I :I :z
115      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ztrdt
116      !!----------------------------------------------------------------------
117
118      IF( wrk_in_use(3, 1,2,3,4,5) .OR. wrk_in_use(2, 1,2,3) )THEN
119         CALL ctl_stop('tra_qsr: requested workspace arrays unavailable')   ;   RETURN
120      ENDIF
121
122      IF( kt == nit000 ) THEN
123         IF(lwp) WRITE(numout,*)
124         IF(lwp) WRITE(numout,*) 'tra_qsr : penetration of the surface solar radiation'
125         IF(lwp) WRITE(numout,*) '~~~~~~~'
126         IF( .NOT.ln_traqsr )   RETURN
127      ENDIF
128
129      IF( l_trdtra ) THEN      ! Save ta and sa trends
130         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
131      ENDIF
132
133      !                                        Set before qsr tracer content field
134      !                                        ***********************************
135      IF( kt == nit000 ) THEN                     ! Set the forcing field at nit000 - 1
136         !                                        ! -----------------------------------
137         IF( ln_rstart .AND.    &                    ! Restart: read in restart file
138              & iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN
139            IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field red in the restart file'
140            zfact = 0.5e0
141            CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b )   ! before heat content trend due to Qsr flux
142         ELSE                                           ! No restart or restart not found: Euler forward time stepping
143            zfact = 1.e0
144            qsr_hc_b(:,:,:) = 0.e0
145         ENDIF
146      ELSE                                        ! Swap of forcing field
147         !                                        ! ---------------------
148         zfact = 0.5e0
149         qsr_hc_b(:,:,:) = qsr_hc(:,:,:)
150      ENDIF
151      !                                        Compute now qsr tracer content field
152      !                                        ************************************
153     
154      !                                           ! ============================================== !
155      IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN      !  bio-model fluxes  : all vertical coordinates  !
156         !                                        ! ============================================== !
157#if defined key_z_first
158         DO jj = 1, jpj
159            DO ji = 1, jpi
160               DO jk = 1, jpkm1
161                  qsr_hc(:,:,jk) = ro0cpr * ( etot3(:,:,jk) - etot3(:,:,jk+1) )
162               END DO
163            END DO
164         END DO
165#else
166         DO jk = 1, jpkm1
167            qsr_hc(:,:,jk) = ro0cpr * ( etot3(:,:,jk) - etot3(:,:,jk+1) )
168         END DO
169#endif
170         !                                        Add to the general trend
171#if defined key_z_first
172         DO jj = 2, jpjm1 
173            DO ji = 2, jpim1
174               DO jk = 1, jpkm1
175#else
176         DO jk = 1, jpkm1
177            DO jj = 2, jpjm1 
178               DO ji = fs_2, fs_jpim1   ! vector opt.
179#endif
180                  z1_e3t = zfact / fse3t(ji,jj,jk)
181                  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
182               END DO
183            END DO
184         END DO
185         CALL iom_put( 'qsr3d', etot3 )   ! Shortwave Radiation 3D distribution
186         !                                        ! ============================================== !
187      ELSE                                        !  Ocean alone :
188         !                                        ! ============================================== !
189         !
190         !                                                ! ------------------------- !
191         IF( ln_qsr_rgb) THEN                             !  R-G-B  light penetration !
192            !                                             ! ------------------------- !
193            ! Set chlorophyl concentration
194            IF( nn_chldta == 1 .OR. lk_vvl ) THEN            !*  Variable Chlorophyll or ocean volume
195               !
196               IF( nn_chldta == 1 ) THEN                             !*  Variable Chlorophyll
197                  !
198                  CALL fld_read( kt, 1, sf_chl )                         ! Read Chl data and provides it at the current time step
199                  !         
200!CDIR COLLAPSE
201!CDIR NOVERRCHK
202                  DO jj = 1, jpj                                         ! Separation in R-G-B depending of the surface Chl
203!CDIR NOVERRCHK
204                     DO ji = 1, jpi
205                        zchl = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) )
206                        irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
207                        zekb(ji,jj) = rkrgb(1,irgb)
208                        zekg(ji,jj) = rkrgb(2,irgb)
209                        zekr(ji,jj) = rkrgb(3,irgb)
210                     END DO
211                  END DO
212               ELSE                                            ! Variable ocean volume but constant chrlorophyll
213                  zchl = 0.05                                     ! constant chlorophyll
214                  irgb = NINT( 41 + 20.*LOG10( zchl ) + 1.e-15 )
215                  zekb(:,:) = rkrgb(1,irgb)                       ! Separation in R-G-B depending of the chlorophyll
216                  zekg(:,:) = rkrgb(2,irgb)
217                  zekr(:,:) = rkrgb(3,irgb)
218               ENDIF
219               !
220               zcoef  = ( 1. - rn_abs ) / 3.e0                        ! equi-partition in R-G-B
221               ze0(:,:,1) = rn_abs  * qsr(:,:)
222               ze1(:,:,1) = zcoef * qsr(:,:)
223               ze2(:,:,1) = zcoef * qsr(:,:)
224               ze3(:,:,1) = zcoef * qsr(:,:)
225               zea(:,:,1) =         qsr(:,:)
226               !
227#if defined key_z_first
228               DO jj = 1, jpj
229                  DO ji = 1, jpi
230                     DO jk = 2, nksr+1
231#else
232               DO jk = 2, nksr+1
233!CDIR NOVERRCHK
234                  DO jj = 1, jpj
235!CDIR NOVERRCHK   
236                     DO ji = 1, jpi
237#endif
238                        zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r     )
239                        zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekb(ji,jj) )
240                        zc2 = ze2(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekg(ji,jj) )
241                        zc3 = ze3(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekr(ji,jj) )
242                        ze0(ji,jj,jk) = zc0
243                        ze1(ji,jj,jk) = zc1
244                        ze2(ji,jj,jk) = zc2
245                        ze3(ji,jj,jk) = zc3
246                        zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk)
247                     END DO
248                  END DO
249               END DO
250               !
251#if defined key_z_first
252               DO jj = 1, jpj
253                  DO ji = 1, jpi
254                     DO jk = 1, nksr                                  ! compute and add qsr trend to ta
255                        qsr_hc(ji,jj,jk) = ro0cpr * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) )
256                     END DO
257                  END DO
258               END DO
259#else
260               DO jk = 1, nksr                                        ! compute and add qsr trend to ta
261                  qsr_hc(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) )
262               END DO
263#endif
264               zea(:,:,nksr+1:jpk) = 0.e0     ! below 400m set to zero
265               CALL iom_put( 'qsr3d', zea )   ! Shortwave Radiation 3D distribution
266               !
267            ELSE                                                 !*  Constant Chlorophyll
268#if defined key_z_first
269               DO jj = 1, jpj
270                  DO ji = 1, jpi
271                     DO jk = 1, nksr
272                        qsr_hc(ji,jj,jk) =  etot3(ji,jj,jk) * qsr(ji,jj)
273                     END DO
274                  END DO
275               END DO
276#else
277               DO jk = 1, nksr
278                  qsr_hc(:,:,jk) =  etot3(:,:,jk) * qsr(:,:)
279               END DO
280#endif
281            ENDIF
282
283         ENDIF
284         !                                                ! ------------------------- !
285         IF( ln_qsr_2bd ) THEN                            !  2 band light penetration !
286            !                                             ! ------------------------- !
287            !
288            IF( lk_vvl ) THEN                                  !* variable volume
289               zz0   =        rn_abs   * ro0cpr
290               zz1   = ( 1. - rn_abs ) * ro0cpr
291#if defined key_z_first
292               DO jj = 2, jpjm1
293                  DO ji = 2, jpim1
294                     DO jk = 1, nksr              ! solar heat absorbed at T-point in the top 400m
295#else
296               DO jk = 1, nksr                    ! solar heat absorbed at T-point in the top 400m
297                  DO jj = 2, jpjm1
298                     DO ji = 2, jpim1
299#endif
300                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
301                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
302                        qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) ) 
303                     END DO
304                  END DO
305               END DO
306            ELSE                                               !* constant volume: coef. computed one for all
307#if defined key_z_first
308               DO jj = 2, jpjm1
309                  DO ji = 2, jpim1
310                     DO jk = 1, nksr
311#else
312               DO jk = 1, nksr
313                  DO jj = 2, jpjm1
314                     DO ji = fs_2, fs_jpim1   ! vector opt.
315#endif
316                        qsr_hc(ji,jj,jk) =  etot3(ji,jj,jk) * qsr(ji,jj)
317                     END DO
318                  END DO
319               END DO
320               !
321            ENDIF
322            !
323         ENDIF
324         !
325         !                                        Add to the general trend
326#if defined key_z_first
327         DO jj = 2, jpjm1 
328            DO ji = 2, jpim1
329               DO jk = 1, nksr
330#else
331         DO jk = 1, nksr
332            DO jj = 2, jpjm1 
333               DO ji = fs_2, fs_jpim1   ! vector opt.
334#endif
335                  z1_e3t = zfact / fse3t(ji,jj,jk)
336                  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
337               END DO
338            END DO
339         END DO
340         !
341      ENDIF
342      !
343      IF( lrst_oce ) THEN   !                  Write in the ocean restart file
344         !                                     *******************************
345         IF(lwp) WRITE(numout,*)
346         IF(lwp) WRITE(numout,*) 'qsr tracer content forcing field written in ocean restart file ',   &
347            &                    'at it= ', kt,' date= ', ndastp
348         IF(lwp) WRITE(numout,*) '~~~~'
349         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b', qsr_hc )
350         !
351      ENDIF
352
353      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
354         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
355         CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_qsr, ztrdt )
356         DEALLOCATE( ztrdt )
357      ENDIF
358      !                       ! print mean trends (used for debugging)
359      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
360      !
361      IF( wrk_not_released(3, 1,2,3,4,5) .OR.   &
362          wrk_not_released(2, 1,2,3)     )   CALL ctl_stop('tra_qsr: failed to release workspace arrays')
363      !
364   END SUBROUTINE tra_qsr
365
366   !! * Reset control of array index permutation
367!FTRANS CLEAR
368#  include "oce_ftrans.h90"
369#  include "dom_oce_ftrans.h90"
370#  include "sbc_oce_ftrans.h90"
371#  include "trc_oce_ftrans.h90"
372
373   SUBROUTINE tra_qsr_init
374      !!----------------------------------------------------------------------
375      !!                  ***  ROUTINE tra_qsr_init  ***
376      !!
377      !! ** Purpose :   Initialization for the penetrative solar radiation
378      !!
379      !! ** Method  :   The profile of solar radiation within the ocean is set
380      !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
381      !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
382      !!      default values correspond to clear water (type I in Jerlov'
383      !!      (1968) classification.
384      !!         called by tra_qsr at the first timestep (nit000)
385      !!
386      !! ** Action  : - initialize rn_si0, rn_si1 and rn_abs
387      !!
388      !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
389      !!----------------------------------------------------------------------
390      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
391      USE wrk_nemo, ONLY:   zekb => wrk_2d_1 , zekg => wrk_2d_2 , zekr => wrk_2d_3
392      USE wrk_nemo, ONLY:   ze0  => wrk_3d_1 , ze1  => wrk_3d_2 , ze2 => wrk_3d_3
393      USE wrk_nemo, ONLY:   ze3  => wrk_3d_4 , zea  => wrk_3d_5
394
395      !! DCSE_NEMO: Need additional directives for renamed module variables
396!FTRANS ze0 ze1 ze2 ze3 zea :I :I :z
397
398      !
399      INTEGER  ::   ji, jj, jk     ! dummy loop indices
400      INTEGER  ::   irgb, ierror, ioptio, nqsr   ! local integer
401      REAL(wp) ::   zz0, zc0  , zc1, zcoef       ! local scalars
402      REAL(wp) ::   zz1, zc2  , zc3, zchl        !   -      -
403      !
404      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files
405      TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read
406      !!
407      NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio,   &
408         &                  nn_chldta, rn_abs, rn_si0, rn_si1
409      !!----------------------------------------------------------------------
410
411      IF( wrk_in_use(2, 1,2,3) .OR. wrk_in_use(3, 1,2,3,4,5) )THEN
412         CALL ctl_stop('tra_qsr_init: requested workspace arrays unavailable')   ;   RETURN
413      ENDIF
414
415      cn_dir = './'       ! directory in which the model is executed
416      ! ... default values (NB: frequency positive => hours, negative => months)
417      !            !     file       ! frequency !  variable  ! time interp !  clim   ! 'yearly' or ! weights  ! rotation   !
418      !            !     name       !  (hours)  !    name    !    (T/F)    !  (T/F)  ! 'monthly'   ! filename ! pairs      !
419      sn_chl = FLD_N( 'chlorophyll' ,    -1     ,  'CHLA'    ,  .true.     , .true.  ,   'yearly'  , ''       , ''         )
420      !
421      REWIND( numnam )            ! Read Namelist namtra_qsr : ratio and length of penetration
422      READ  ( numnam, namtra_qsr )
423      !
424      IF(lwp) THEN                ! control print
425         WRITE(numout,*)
426         WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
427         WRITE(numout,*) '~~~~~~~~~~~~'
428         WRITE(numout,*) '   Namelist namtra_qsr : set the parameter of penetration'
429         WRITE(numout,*) '      Light penetration (T) or not (F)         ln_traqsr  = ', ln_traqsr
430         WRITE(numout,*) '      RGB (Red-Green-Blue) light penetration   ln_qsr_rgb = ', ln_qsr_rgb
431         WRITE(numout,*) '      2 band               light penetration   ln_qsr_2bd = ', ln_qsr_2bd
432         WRITE(numout,*) '      bio-model            light penetration   ln_qsr_bio = ', ln_qsr_bio
433         WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)    nn_chldta  = ', nn_chldta
434         WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs = ', rn_abs
435         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0 = ', rn_si0
436         WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1 = ', rn_si1
437      ENDIF
438
439      IF( ln_traqsr ) THEN     ! control consistency
440         !                     
441         IF( .NOT.lk_qsr_bio .AND. ln_qsr_bio )   THEN
442            CALL ctl_warn( 'No bio model : force ln_qsr_bio = FALSE ' )
443            ln_qsr_bio = .FALSE.
444         ENDIF
445         !
446         ioptio = 0                      ! Parameter control
447         IF( ln_qsr_rgb  )   ioptio = ioptio + 1
448         IF( ln_qsr_2bd  )   ioptio = ioptio + 1
449         IF( ln_qsr_bio  )   ioptio = ioptio + 1
450         !
451         IF( ioptio /= 1 ) &
452            CALL ctl_stop( '          Choose ONE type of light penetration in namelist namtra_qsr',  &
453            &              ' 2 bands, 3 RGB bands or bio-model light penetration' )
454         !
455         IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr =  1 
456         IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr =  2
457         IF( ln_qsr_2bd                      )   nqsr =  3
458         IF( ln_qsr_bio                      )   nqsr =  4
459         !
460         IF(lwp) THEN                   ! Print the choice
461            WRITE(numout,*)
462            IF( nqsr ==  1 )   WRITE(numout,*) '         R-G-B   light penetration - Constant Chlorophyll'
463            IF( nqsr ==  2 )   WRITE(numout,*) '         R-G-B   light penetration - Chl data '
464            IF( nqsr ==  3 )   WRITE(numout,*) '         2 bands light penetration'
465            IF( nqsr ==  4 )   WRITE(numout,*) '         bio-model light penetration'
466         ENDIF
467         !
468      ENDIF
469      !                          ! ===================================== !
470      IF( ln_traqsr  ) THEN      !  Initialisation of Light Penetration  ! 
471         !                       ! ===================================== !
472         !
473         xsi0r = 1.e0 / rn_si0
474         xsi1r = 1.e0 / rn_si1
475         !                                ! ---------------------------------- !
476         IF( ln_qsr_rgb ) THEN            !  Red-Green-Blue light penetration  !
477            !                             ! ---------------------------------- !
478            !
479            CALL trc_oce_rgb( rkrgb )           !* tabulated attenuation coef.
480            !
481            !                                   !* level of light extinction
482            IF(  ln_sco ) THEN   ;   nksr = jpkm1
483            ELSE                 ;   nksr = trc_oce_ext_lev( r_si2, 0.33e2 )
484            ENDIF
485
486            IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_0(nksr+1), ' m'
487            !
488            IF( nn_chldta == 1 ) THEN           !* Chl data : set sf_chl structure
489               IF(lwp) WRITE(numout,*)
490               IF(lwp) WRITE(numout,*) '        Chlorophyll read in a file'
491               ALLOCATE( sf_chl(1), STAT=ierror )
492               IF( ierror > 0 ) THEN
493                  CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' )   ;   RETURN
494               ENDIF
495               ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1)   )
496               IF( sn_chl%ln_tint )ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) )
497               !                                        ! fill sf_chl with sn_chl and control print
498               CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',   &
499                  &                                         'Solar penetration function of read chlorophyll', 'namtra_qsr' )
500               !
501            ELSE                                !* constant Chl : compute once for all the distribution of light (etot3)
502               IF(lwp) WRITE(numout,*)
503               IF(lwp) WRITE(numout,*) '        Constant Chlorophyll concentration = 0.05'
504               IF( lk_vvl ) THEN                   ! variable volume
505                  IF(lwp) WRITE(numout,*) '        key_vvl: light distribution will be computed at each time step'
506               ELSE                                ! constant volume: computes one for all
507                  IF(lwp) WRITE(numout,*) '        fixed volume: light distribution computed one for all'
508                  !
509                  zchl = 0.05                                 ! constant chlorophyll
510                  irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
511                  zekb(:,:) = rkrgb(1,irgb)                   ! Separation in R-G-B depending of the chlorophyll
512                  zekg(:,:) = rkrgb(2,irgb)
513                  zekr(:,:) = rkrgb(3,irgb)
514                  !
515                  zcoef = ( 1. - rn_abs ) / 3.e0              ! equi-partition in R-G-B
516               
517#if defined key_z_first
518                  DO jj = 1, jpj
519                     DO ji = 1, jpi
520                        ze0(ji,jj,1) = rn_abs
521                        ze1(ji,jj,1) = zcoef
522                        ze2(ji,jj,1) = zcoef 
523                        ze3(ji,jj,1) = zcoef
524                        zea(ji,jj,1) = tmask(ji,jj,1)         ! = ( ze0+ze1+z2+ze3 ) * tmask
525                        DO jk = 2, nksr+1
526#else
527                  ze0(:,:,1) = rn_abs
528                  ze1(:,:,1) = zcoef
529                  ze2(:,:,1) = zcoef 
530                  ze3(:,:,1) = zcoef
531                  zea(:,:,1) = tmask(:,:,1)                   ! = ( ze0+ze1+z2+ze3 ) * tmask
532                  DO jk = 2, nksr+1
533!CDIR NOVERRCHK
534                     DO jj = 1, jpj
535!CDIR NOVERRCHK   
536                        DO ji = 1, jpi
537#endif
538                           zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t_0(ji,jj,jk-1) * xsi0r     )
539                           zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t_0(ji,jj,jk-1) * zekb(ji,jj) )
540                           zc2 = ze2(ji,jj,jk-1) * EXP( - fse3t_0(ji,jj,jk-1) * zekg(ji,jj) )
541                           zc3 = ze3(ji,jj,jk-1) * EXP( - fse3t_0(ji,jj,jk-1) * zekr(ji,jj) )
542                           ze0(ji,jj,jk) = zc0
543                           ze1(ji,jj,jk) = zc1
544                           ze2(ji,jj,jk) = zc2
545                           ze3(ji,jj,jk) = zc3
546                           zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk)
547                        END DO
548                     END DO
549                  END DO 
550                  !
551#if defined key_z_first
552                  DO jj = 1, jpj
553                     DO ji = 1, jpi
554                        DO jk = 1, nksr
555                           etot3(ji,jj,jk) = ro0cpr * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) ) 
556                        END DO
557                     END DO
558                  END DO
559#else
560                  DO jk = 1, nksr
561                     etot3(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) 
562                  END DO
563#endif
564                  etot3(:,:,nksr+1:jpk) = 0.e0                ! below 400m set to zero
565               ENDIF
566            ENDIF
567            !
568         ENDIF
569            !                             ! ---------------------------------- !
570         IF( ln_qsr_2bd ) THEN            !    2 bands    light penetration    !
571            !                             ! ---------------------------------- !
572            !
573            !                                ! level of light extinction
574            nksr = trc_oce_ext_lev( rn_si1, 1.e2 )
575            IF(lwp) THEN
576               WRITE(numout,*)
577            IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_0(nksr+1), ' m'
578            ENDIF
579            !
580            IF( lk_vvl ) THEN                   ! variable volume
581               IF(lwp) WRITE(numout,*) '        key_vvl: light distribution will be computed at each time step'
582            ELSE                                ! constant volume: computes one for all
583               zz0 =        rn_abs   * ro0cpr
584               zz1 = ( 1. - rn_abs ) * ro0cpr
585#if defined key_z_first
586               DO jj = 1, jpj                     !*  solar heat absorbed at T-point computed once for all
587                  DO ji = 1, jpi
588                     DO jk = 1, nksr                         ! top 400 meters
589                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
590                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
591                        etot3(ji,jj,jk) = (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) 
592                     END DO
593                     DO jk = nksr+1, jpk
594                        etot3(ji,jj,jk) = 0.e0       ! below 400m set to zero
595                     END DO
596                  END DO
597               END DO
598#else
599               DO jk = 1, nksr                    !*  solar heat absorbed at T-point computed once for all
600                  DO jj = 1, jpj                              ! top 400 meters
601                     DO ji = 1, jpi
602                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
603                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
604                        etot3(ji,jj,jk) = (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) 
605                     END DO
606                  END DO
607               END DO
608               etot3(:,:,nksr+1:jpk) = 0.e0                   ! below 400m set to zero
609#endif
610               !
611            ENDIF
612         ENDIF
613         !                       ! ===================================== !
614      ELSE                       !        No light penetration           !                   
615         !                       ! ===================================== !
616         IF(lwp) THEN
617            WRITE(numout,*)
618            WRITE(numout,*) 'tra_qsr_init : NO solar flux penetration'
619            WRITE(numout,*) '~~~~~~~~~~~~'
620         ENDIF
621      ENDIF
622      !
623      IF( wrk_not_released(2, 1,2,3)     .OR.   &
624          wrk_not_released(3, 1,2,3,4,5) )   CALL ctl_stop('tra_qsr_init: failed to release workspace arrays')
625      !
626   END SUBROUTINE tra_qsr_init
627
628   !!======================================================================
629END MODULE traqsr
Note: See TracBrowser for help on using the repository browser.