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

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90 @ 7525

Last change on this file since 7525 was 7525, checked in by mocavero, 7 years ago

changes after review
  • Property svn:keywords set to Id
File size: 24.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   !!            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 sbc_oce        ! surface boundary condition: ocean
25   USE trc_oce        ! share SMS/Ocean variables
26   USE trd_oce        ! trends: ocean variables
27   USE trdtra         ! trends manager: tracers
28   !
29   USE in_out_manager ! I/O manager
30   USE prtctl         ! Print control
31   USE iom            ! I/O manager
32   USE fldread        ! read input fields
33   USE restart        ! ocean restart
34   USE lib_mpp        ! MPP library
35   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
36   USE wrk_nemo       ! Memory Allocation
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   LOGICAL , PUBLIC ::   ln_qsr_ice   !: light penetration for ice-model LIM3 (clem)
51   INTEGER , PUBLIC ::   nn_chldta    !: use Chlorophyll data (=1) or not (=0)
52   REAL(wp), PUBLIC ::   rn_abs       !: fraction absorbed in the very near surface (RGB & 2 bands)
53   REAL(wp), PUBLIC ::   rn_si0       !: very near surface depth of extinction      (RGB & 2 bands)
54   REAL(wp), PUBLIC ::   rn_si1       !: deepest depth of extinction (water type I)       (2 bands)
55   !
56   INTEGER , PUBLIC ::   nksr         !: levels below which the light cannot penetrate (depth larger than 391 m)
57 
58   INTEGER, PARAMETER ::   np_RGB  = 1   ! R-G-B     light penetration with constant Chlorophyll
59   INTEGER, PARAMETER ::   np_RGBc = 2   ! R-G-B     light penetration with Chlorophyll data
60   INTEGER, PARAMETER ::   np_2BD  = 3   ! 2 bands   light penetration
61   INTEGER, PARAMETER ::   np_BIO  = 4   ! bio-model light penetration
62   !
63   INTEGER  ::   nqsr    ! user choice of the type of light penetration
64   REAL(wp) ::   xsi0r   ! inverse of rn_si0
65   REAL(wp) ::   xsi1r   ! inverse of rn_si1
66   !
67   REAL(wp) , DIMENSION(3,61)           ::   rkrgb    ! tabulated attenuation coefficients for RGB absorption
68   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read)
69
70   !! * Substitutions
71#  include "vectopt_loop_substitute.h90"
72   !!----------------------------------------------------------------------
73   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
74   !! $Id$
75   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
76   !!----------------------------------------------------------------------
77CONTAINS
78
79   SUBROUTINE tra_qsr( kt )
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 ] / (rau0*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      !
107      INTEGER  ::   ji, jj, jk               ! dummy loop indices
108      INTEGER  ::   irgb                     ! local integers
109      REAL(wp) ::   zchl, zcoef, z1_2        ! local scalars
110      REAL(wp) ::   zc0 , zc1 , zc2 , zc3    !    -         -
111      REAL(wp) ::   zzc0, zzc1, zzc2, zzc3   !    -         -
112      REAL(wp) ::   zz0 , zz1                !    -         -
113      REAL(wp) ::   zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze
114      REAL(wp) ::   zlogc, zlogc2, zlogc3 
115      REAL(wp), POINTER, DIMENSION(:,:)   :: zekb, zekg, zekr
116      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt
117      REAL(wp), POINTER, DIMENSION(:,:,:) :: zetot, zchl3d
118      !!----------------------------------------------------------------------
119      !
120      IF( nn_timing == 1 )  CALL timing_start('tra_qsr')
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      ENDIF
127      !
128      IF( l_trdtra ) THEN      ! trends diagnostic: save the input temperature trend
129         CALL wrk_alloc( jpi,jpj,jpk,   ztrdt ) 
130!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
131            DO jk = 1, jpk
132               DO jj = 1, jpj
133                  DO ji = 1, jpi
134                     ztrdt(ji,jj,jk) = tsa(ji,jj,jk,jp_tem)
135                  END DO
136               END DO
137            END DO
138      ENDIF
139      !
140      !                         !-----------------------------------!
141      !                         !  before qsr induced heat content  !
142      !                         !-----------------------------------!
143      IF( kt == nit000 ) THEN          !==  1st time step  ==!
144!!gm case neuler  not taken into account....
145         IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN    ! read in restart
146            IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field read in the restart file'
147            z1_2 = 0.5_wp
148            CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b )   ! before heat content trend due to Qsr flux
149         ELSE                                           ! No restart or restart not found: Euler forward time stepping
150            z1_2 = 1._wp
151!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
152            DO jk = 1, jpk
153               DO jj = 1, jpj
154                  DO ji = 1, jpi
155                     qsr_hc_b(ji,jj,jk) = 0._wp
156                  END DO
157               END DO
158            END DO
159         ENDIF
160      ELSE                             !==  Swap of qsr heat content  ==!
161         z1_2 = 0.5_wp
162!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
163            DO jk = 1, jpk
164               DO jj = 1, jpj
165                  DO ji = 1, jpi
166                     qsr_hc_b(ji,jj,jk) = qsr_hc(ji,jj,jk)
167                  END DO
168               END DO
169            END DO
170      ENDIF
171      !
172      !                         !--------------------------------!
173      SELECT CASE( nqsr )       !  now qsr induced heat content  !
174      !                         !--------------------------------!
175      !
176      CASE( np_BIO )                   !==  bio-model fluxes  ==!
177         !
178!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
179         DO jk = 1, nksr
180            DO jj = 1, jpj
181               DO ji = 1, jpi
182                  qsr_hc(ji,jj,jk) = r1_rau0_rcp * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) )
183               END DO
184             END DO
185         END DO
186         !
187      CASE( np_RGB , np_RGBc )         !==  R-G-B fluxes  ==!
188         !
189         CALL wrk_alloc( jpi,jpj,       zekb, zekg, zekr        ) 
190         CALL wrk_alloc( jpi,jpj,jpk,   ze0, ze1, ze2, ze3, zea, zchl3d ) 
191         !
192         IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll
193            CALL fld_read( kt, 1, sf_chl )         ! Read Chl data and provides it at the current time step
194!$OMP PARALLEL DO schedule(static) private(jk,jj,ji,zchl,zCtot,zze,zpsi,zlogc,zlogc2,zlogc3,zCb,zCmax,zpsimax,zdelpsi,zCze)
195            DO jk = 1, nksr + 1
196               DO jj = 2, jpjm1                       ! Separation in R-G-B depending of the surface Chl
197                  DO ji = fs_2, fs_jpim1
198                     zchl    = sf_chl(1)%fnow(ji,jj,1)
199                     zCtot   = 40.6  * zchl**0.459
200                     zze     = 568.2 * zCtot**(-0.746)
201                     IF( zze > 102. ) zze = 200.0 * zCtot**(-0.293)
202                     zpsi    = gdepw_n(ji,jj,jk) / zze
203                     !
204                     zlogc   = LOG( zchl )
205                     zlogc2  = zlogc * zlogc
206                     zlogc3  = zlogc * zlogc * zlogc
207                     zCb     = 0.768 + 0.087 * zlogc - 0.179 * zlogc2 - 0.025 * zlogc3
208                     zCmax   = 0.299 - 0.289 * zlogc + 0.579 * zlogc2
209                     zpsimax = 0.6   - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3
210                     zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2
211                     zCze    = 1.12  * (zchl)**0.803 
212                     !
213                     zchl3d(ji,jj,jk) = zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) / zdelpsi )**2 ) )
214                  END DO
215                  !
216               END DO
217            END DO
218         ELSE                                !* constant chrlorophyll
219!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
220           DO jk = 1, nksr + 1
221              DO jj = 1, jpj
222                 DO ji = 1, jpi
223                    zchl3d(ji,jj,jk) = 0.05 
224                 ENDDO
225              ENDDO
226            ENDDO
227         ENDIF
228         !
229         zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B
230!$OMP PARALLEL
231!$OMP DO schedule(static) private(jj,ji)
232         DO jj = 2, jpjm1
233            DO ji = fs_2, fs_jpim1
234               ze0(ji,jj,1) = rn_abs * qsr(ji,jj)
235               ze1(ji,jj,1) = zcoef  * qsr(ji,jj)
236               ze2(ji,jj,1) = zcoef  * qsr(ji,jj)
237               ze3(ji,jj,1) = zcoef  * qsr(ji,jj)
238               zea(ji,jj,1) =          qsr(ji,jj)
239            END DO
240         END DO
241!$OMP END DO NOWAIT
242         !
243         DO jk = 2, nksr+1                   !* interior equi-partition in R-G-B depending of vertical profile of Chl
244!$OMP DO schedule(static) private(jj,ji,zchl,irgb)
245            DO jj = 2, jpjm1
246               DO ji = fs_2, fs_jpim1
247                  zchl = MIN( 10. , MAX( 0.03, zchl3d(ji,jj,jk) ) )
248                  irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
249                  zekb(ji,jj) = rkrgb(1,irgb)
250                  zekg(ji,jj) = rkrgb(2,irgb)
251                  zekr(ji,jj) = rkrgb(3,irgb)
252               END DO
253            END DO
254
255!$OMP DO schedule(static) private(jj,ji,zc0,zc1,zc2,zc3)
256            DO jj = 2, jpjm1
257               DO ji = fs_2, fs_jpim1
258                  zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * xsi0r       )
259                  zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekb(ji,jj) )
260                  zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekg(ji,jj) )
261                  zc3 = ze3(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekr(ji,jj) )
262                  ze0(ji,jj,jk) = zc0
263                  ze1(ji,jj,jk) = zc1
264                  ze2(ji,jj,jk) = zc2
265                  ze3(ji,jj,jk) = zc3
266                  zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk)
267               END DO
268            END DO
269         END DO
270         !
271!$OMP DO schedule(static) private(jk,jj,ji)
272         DO jk = 1, nksr                     !* now qsr induced heat content
273            DO jj = 2, jpjm1
274               DO ji = fs_2, fs_jpim1
275                  qsr_hc(ji,jj,jk) = r1_rau0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) )
276               END DO
277            END DO
278         END DO
279!$OMP END DO NOWAIT
280!$OMP END PARALLEL
281         !
282         CALL wrk_dealloc( jpi,jpj,        zekb, zekg, zekr        ) 
283         CALL wrk_dealloc( jpi,jpj,jpk,   ze0, ze1, ze2, ze3, zea, zchl3d ) 
284         !
285      CASE( np_2BD  )            !==  2-bands fluxes  ==!
286         !
287         zz0 =        rn_abs   * r1_rau0_rcp      ! surface equi-partition in 2-bands
288         zz1 = ( 1. - rn_abs ) * r1_rau0_rcp
289!$OMP PARALLEL DO schedule(static) private(jk,jj,ji,zc0,zc1)
290         DO jk = 1, nksr                          ! solar heat absorbed at T-point in the top 400m
291            DO jj = 2, jpjm1
292               DO ji = fs_2, fs_jpim1
293                  zc0 = zz0 * EXP( -gdepw_n(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk  )*xsi1r )
294                  zc1 = zz0 * EXP( -gdepw_n(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk+1)*xsi1r )
295                  qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 
296               END DO
297            END DO
298         END DO
299         !
300      END SELECT
301      !
302      !                          !-----------------------------!
303!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
304      DO jk = 1, nksr            !  update to the temp. trend  !
305         DO jj = 2, jpjm1        !-----------------------------!
306            DO ji = fs_2, fs_jpim1   ! vector opt.
307               tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
308                  &                 + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t_n(ji,jj,jk)
309            END DO
310         END DO
311      END DO
312      !
313      IF( ln_qsr_ice ) THEN      ! sea-ice: store the 1st ocean level attenuation coefficient
314!$OMP PARALLEL DO schedule(static) private(jj,ji)
315         DO jj = 2, jpjm1 
316            DO ji = fs_2, fs_jpim1   ! vector opt.
317               IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) )
318               ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp
319               ENDIF
320            END DO
321         END DO
322         ! Update haloes since lim_thd needs fraqsr_1lev to be defined everywhere
323         CALL lbc_lnk( fraqsr_1lev(:,:), 'T', 1._wp )
324      ENDIF
325      !
326      IF( iom_use('qsr3d') ) THEN      ! output the shortwave Radiation distribution
327         CALL wrk_alloc( jpi,jpj,jpk,   zetot )
328         !
329!$OMP PARALLEL
330!$OMP DO schedule(static) private(jj,ji)
331         DO jj = 1, jpj 
332            DO ji = 1, jpi   ! vector opt.
333               zetot(ji,jj,nksr+1:jpk) = 0._wp     ! below ~400m set to zero
334            END DO
335         END DO
336         DO jk = nksr, 1, -1
337!$OMP DO schedule(static) private(jj,ji)
338            DO jj = 1, jpj 
339               DO ji = 1, jpi   ! vector opt.
340                  zetot(ji,jj,jk) = zetot(ji,jj,jk+1) + qsr_hc(ji,jj,jk) / r1_rau0_rcp
341               END DO
342            END DO
343!$OMP END DO NOWAIT
344         END DO         
345!$OMP END PARALLEL
346         CALL iom_put( 'qsr3d', zetot )   ! 3D distribution of shortwave Radiation
347         !
348         CALL wrk_dealloc( jpi,jpj,jpk,   zetot ) 
349      ENDIF
350      !
351      IF( lrst_oce ) THEN     ! write in the ocean restart file
352         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc      )
353         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev ) 
354      ENDIF
355      !
356      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
357!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
358         DO jk = 1, jpk
359            DO jj = 1, jpj
360               DO ji = 1, jpi
361                  ztrdt(ji,jj,jk) = tsa(ji,jj,jk,jp_tem) - ztrdt(ji,jj,jk)
362               END DO
363            END DO
364         END DO
365         CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt )
366         CALL wrk_dealloc( jpi,jpj,jpk,   ztrdt ) 
367      ENDIF
368      !                       ! print mean trends (used for debugging)
369      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
370      !
371      IF( nn_timing == 1 )  CALL timing_stop('tra_qsr')
372      !
373   END SUBROUTINE tra_qsr
374
375
376   SUBROUTINE tra_qsr_init
377      !!----------------------------------------------------------------------
378      !!                  ***  ROUTINE tra_qsr_init  ***
379      !!
380      !! ** Purpose :   Initialization for the penetrative solar radiation
381      !!
382      !! ** Method  :   The profile of solar radiation within the ocean is set
383      !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
384      !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
385      !!      default values correspond to clear water (type I in Jerlov'
386      !!      (1968) classification.
387      !!         called by tra_qsr at the first timestep (nit000)
388      !!
389      !! ** Action  : - initialize rn_si0, rn_si1 and rn_abs
390      !!
391      !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
392      !!----------------------------------------------------------------------
393      INTEGER  ::   ji, jj, jk                  ! dummy loop indices
394      INTEGER  ::   ios, irgb, ierror, ioptio   ! local integer
395      REAL(wp) ::   zz0, zc0 , zc1, zcoef      ! local scalars
396      REAL(wp) ::   zz1, zc2 , zc3, zchl       !   -      -
397      !
398      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files
399      TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read
400      !!
401      NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice,  &
402         &                  nn_chldta, rn_abs, rn_si0, rn_si1
403      !!----------------------------------------------------------------------
404      !
405      IF( nn_timing == 1 )   CALL timing_start('tra_qsr_init')
406      !
407      REWIND( numnam_ref )              ! Namelist namtra_qsr in reference     namelist
408      READ  ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901)
409901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_qsr in reference namelist', lwp )
410      !
411      REWIND( numnam_cfg )              ! Namelist namtra_qsr in configuration namelist
412      READ  ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 )
413902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist', lwp )
414      IF(lwm) WRITE ( numond, namtra_qsr )
415      !
416      IF(lwp) THEN                ! control print
417         WRITE(numout,*)
418         WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
419         WRITE(numout,*) '~~~~~~~~~~~~'
420         WRITE(numout,*) '   Namelist namtra_qsr : set the parameter of penetration'
421         WRITE(numout,*) '      RGB (Red-Green-Blue) light penetration       ln_qsr_rgb = ', ln_qsr_rgb
422         WRITE(numout,*) '      2 band               light penetration       ln_qsr_2bd = ', ln_qsr_2bd
423         WRITE(numout,*) '      bio-model            light penetration       ln_qsr_bio = ', ln_qsr_bio
424         WRITE(numout,*) '      light penetration for ice-model (LIM3)       ln_qsr_ice = ', ln_qsr_ice
425         WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)        nn_chldta  = ', nn_chldta
426         WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs     = ', rn_abs
427         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0     = ', rn_si0
428         WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1     = ', rn_si1
429         WRITE(numout,*)
430      ENDIF
431      !
432      ioptio = 0                    ! Parameter control
433      IF( ln_qsr_rgb  )   ioptio = ioptio + 1
434      IF( ln_qsr_2bd  )   ioptio = ioptio + 1
435      IF( ln_qsr_bio  )   ioptio = ioptio + 1
436      !
437      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE type of light penetration in namelist namtra_qsr',  &
438         &                               ' 2 bands, 3 RGB bands or bio-model light penetration' )
439      !
440      IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr = np_RGB 
441      IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr = np_RGBc
442      IF( ln_qsr_2bd                      )   nqsr = np_2BD
443      IF( ln_qsr_bio                      )   nqsr = np_BIO
444      !
445      !                             ! Initialisation
446      xsi0r = 1._wp / rn_si0
447      xsi1r = 1._wp / rn_si1
448      !
449      SELECT CASE( nqsr )
450      !                               
451      CASE( np_RGB , np_RGBc )         !==  Red-Green-Blue light penetration  ==!
452         !                             
453         IF(lwp)   WRITE(numout,*) '   R-G-B   light penetration '
454         !
455         CALL trc_oce_rgb( rkrgb )                 ! tabulated attenuation coef.
456         !                                   
457         nksr = trc_oce_ext_lev( r_si2, 33._wp )   ! level of light extinction
458         !
459         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
460         !
461         IF( nqsr == np_RGBc ) THEN                ! Chl data : set sf_chl structure
462            IF(lwp) WRITE(numout,*) '        Chlorophyll read in a file'
463            ALLOCATE( sf_chl(1), STAT=ierror )
464            IF( ierror > 0 ) THEN
465               CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' )   ;   RETURN
466            ENDIF
467            ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1)   )
468            IF( sn_chl%ln_tint )   ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) )
469            !                                        ! fill sf_chl with sn_chl and control print
470            CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',   &
471               &           'Solar penetration function of read chlorophyll', 'namtra_qsr' )
472         ENDIF
473         IF( nqsr == np_RGB ) THEN                 ! constant Chl
474            IF(lwp) WRITE(numout,*) '        Constant Chlorophyll concentration = 0.05'
475         ENDIF
476         !
477      CASE( np_2BD )                   !==  2 bands light penetration  ==!
478         !
479         IF(lwp)  WRITE(numout,*) '   2 bands light penetration'
480         !
481         nksr = trc_oce_ext_lev( rn_si1, 100._wp )    ! level of light extinction
482         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
483         !
484      CASE( np_BIO )                   !==  BIO light penetration  ==!
485         !
486         IF(lwp) WRITE(numout,*) '   bio-model light penetration'
487         IF( .NOT.lk_qsr_bio )   CALL ctl_stop( 'No bio model : ln_qsr_bio = true impossible ' )
488         !
489      END SELECT
490      !
491!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
492      DO jk = 1, jpk
493         DO jj = 1, jpj
494            DO ji = 1, jpi
495               qsr_hc(ji,jj,jk) = 0._wp     ! now qsr heat content set to zero where it will not be computed
496            END DO
497         END DO
498      END DO
499      !
500      ! 1st ocean level attenuation coefficient (used in sbcssm)
501      IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN
502         CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev'  , fraqsr_1lev  )
503      ELSE
504!$OMP PARALLEL DO schedule(static) private(jj,ji)
505         DO jj = 1, jpj
506            DO ji = 1, jpi
507               fraqsr_1lev(ji,jj) = 1._wp   ! default : no penetration
508            END DO
509         END DO
510      ENDIF
511      !
512      IF( nn_timing == 1 )   CALL timing_stop('tra_qsr_init')
513      !
514   END SUBROUTINE tra_qsr_init
515
516   !!======================================================================
517END MODULE traqsr
Note: See TracBrowser for help on using the repository browser.