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

source: branches/2017/dev_r8600_xios_read_write/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90 @ 8770

Last change on this file since 8770 was 8770, checked in by andmirek, 6 years ago

#1953 and #1962 merged (and some modifications) with write branch

  • Property svn:keywords set to Id
File size: 22.5 KB
RevLine 
[3]1MODULE traqsr
2   !!======================================================================
3   !!                       ***  MODULE  traqsr  ***
[6140]4   !! Ocean physics:   solar radiation penetration in the top ocean levels
[3]5   !!======================================================================
[1423]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)
[6140]12   !!            3.6  !  2012-05  (C. Rousset) store attenuation coef for use in ice model
[6403]13   !!            3.6  !  2015-12  (O. Aumont, J. Jouanno, C. Ethe) use vertical profile of chlorophyll
[6140]14   !!            3.7  !  2015-11  (G. Madec, A. Coward)  remove optimisation for fix volume
[3]15   !!----------------------------------------------------------------------
[503]16
17   !!----------------------------------------------------------------------
[6140]18   !!   tra_qsr       : temperature trend due to the penetration of solar radiation
19   !!   tra_qsr_init  : initialization of the qsr penetration
[3]20   !!----------------------------------------------------------------------
[6140]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
[4990]26   USE trd_oce        ! trends: ocean variables
27   USE trdtra         ! trends manager: tracers
[6140]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)
[3294]36   USE wrk_nemo       ! Memory Allocation
37   USE timing         ! Timing
[8770]38   USE iom_def, ONLY : lxios_read, lwxios
[3]39   IMPLICIT NONE
40   PRIVATE
41
[2528]42   PUBLIC   tra_qsr       ! routine called by step.F90 (ln_traqsr=T)
[5407]43   PUBLIC   tra_qsr_init  ! routine called by nemogcm.F90
[3]44
[4147]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
[4205]50   LOGICAL , PUBLIC ::   ln_qsr_ice   !: light penetration for ice-model LIM3 (clem)
[4147]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)
[6140]55   !
56   INTEGER , PUBLIC ::   nksr         !: levels below which the light cannot penetrate (depth larger than 391 m)
[5407]57 
[6140]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
[1423]68   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read)
[3]69
70   !! * Substitutions
71#  include "vectopt_loop_substitute.h90"
72   !!----------------------------------------------------------------------
[2528]73   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
[888]74   !! $Id$
[2715]75   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[3]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
[6140]84      !!              penetration and add it to the general temperature trend.
[3]85      !!
[1423]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)
[3]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.
[6140]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) .
[3]97      !!
98      !! ** Action  : - update ta with the penetrative solar radiation trend
[6140]99      !!              - send  trend for further diagnostics (l_trdtra=T)
[1423]100      !!
101      !! Reference  : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
102      !!              Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
[6403]103      !!              Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562
[503]104      !!----------------------------------------------------------------------
105      INTEGER, INTENT(in) ::   kt     ! ocean time-step
[2715]106      !
[6140]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    !    -         -
[4161]111      REAL(wp) ::   zzc0, zzc1, zzc2, zzc3   !    -         -
[6140]112      REAL(wp) ::   zz0 , zz1                !    -         -
[6403]113      REAL(wp) ::   zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze
114      REAL(wp) ::   zlogc, zlogc2, zlogc3 
[6140]115      REAL(wp), POINTER, DIMENSION(:,:)   :: zekb, zekg, zekr
[3294]116      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt
[6403]117      REAL(wp), POINTER, DIMENSION(:,:,:) :: zetot, zchl3d
[3]118      !!----------------------------------------------------------------------
[3294]119      !
120      IF( nn_timing == 1 )  CALL timing_start('tra_qsr')
121      !
[3]122      IF( kt == nit000 ) THEN
[503]123         IF(lwp) WRITE(numout,*)
124         IF(lwp) WRITE(numout,*) 'tra_qsr : penetration of the surface solar radiation'
125         IF(lwp) WRITE(numout,*) '~~~~~~~'
[3]126      ENDIF
[6140]127      !
128      IF( l_trdtra ) THEN      ! trends diagnostic: save the input temperature trend
129         CALL wrk_alloc( jpi,jpj,jpk,   ztrdt ) 
[7753]130         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
[216]131      ENDIF
[6140]132      !
133      !                         !-----------------------------------!
134      !                         !  before qsr induced heat content  !
135      !                         !-----------------------------------!
136      IF( kt == nit000 ) THEN          !==  1st time step  ==!
137!!gm case neuler  not taken into account....
138         IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN    ! read in restart
139            IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field read in the restart file'
140            z1_2 = 0.5_wp
[8668]141            CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b, ldxios = lxios_read )   ! before heat content trend due to Qsr flux
[2528]142         ELSE                                           ! No restart or restart not found: Euler forward time stepping
[6140]143            z1_2 = 1._wp
[7753]144            qsr_hc_b(:,:,:) = 0._wp
[2528]145         ENDIF
[6140]146      ELSE                             !==  Swap of qsr heat content  ==!
147         z1_2 = 0.5_wp
[7753]148         qsr_hc_b(:,:,:) = qsr_hc(:,:,:)
[2528]149      ENDIF
[6140]150      !
151      !                         !--------------------------------!
152      SELECT CASE( nqsr )       !  now qsr induced heat content  !
153      !                         !--------------------------------!
154      !
155      CASE( np_BIO )                   !==  bio-model fluxes  ==!
156         !
157         DO jk = 1, nksr
[7753]158            qsr_hc(:,:,jk) = r1_rau0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) )
[2528]159         END DO
[6140]160         !
161      CASE( np_RGB , np_RGBc )         !==  R-G-B fluxes  ==!
162         !
163         CALL wrk_alloc( jpi,jpj,       zekb, zekg, zekr        ) 
[6403]164         CALL wrk_alloc( jpi,jpj,jpk,   ze0, ze1, ze2, ze3, zea, zchl3d ) 
[6140]165         !
166         IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll
167            CALL fld_read( kt, 1, sf_chl )         ! Read Chl data and provides it at the current time step
[6403]168            DO jk = 1, nksr + 1
169               DO jj = 2, jpjm1                       ! Separation in R-G-B depending of the surface Chl
170                  DO ji = fs_2, fs_jpim1
171                     zchl    = sf_chl(1)%fnow(ji,jj,1)
172                     zCtot   = 40.6  * zchl**0.459
173                     zze     = 568.2 * zCtot**(-0.746)
174                     IF( zze > 102. ) zze = 200.0 * zCtot**(-0.293)
175                     zpsi    = gdepw_n(ji,jj,jk) / zze
176                     !
177                     zlogc   = LOG( zchl )
178                     zlogc2  = zlogc * zlogc
179                     zlogc3  = zlogc * zlogc * zlogc
180                     zCb     = 0.768 + 0.087 * zlogc - 0.179 * zlogc2 - 0.025 * zlogc3
181                     zCmax   = 0.299 - 0.289 * zlogc + 0.579 * zlogc2
182                     zpsimax = 0.6   - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3
183                     zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2
184                     zCze    = 1.12  * (zchl)**0.803 
185                     !
186                     zchl3d(ji,jj,jk) = zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) / zdelpsi )**2 ) )
187                  END DO
188                  !
[3]189               END DO
190            END DO
[6140]191         ELSE                                !* constant chrlorophyll
[6403]192           DO jk = 1, nksr + 1
[7753]193              zchl3d(:,:,jk) = 0.05 
[6403]194            ENDDO
[4161]195         ENDIF
[1423]196         !
[6140]197         zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B
198         DO jj = 2, jpjm1
199            DO ji = fs_2, fs_jpim1
200               ze0(ji,jj,1) = rn_abs * qsr(ji,jj)
201               ze1(ji,jj,1) = zcoef  * qsr(ji,jj)
202               ze2(ji,jj,1) = zcoef  * qsr(ji,jj)
203               ze3(ji,jj,1) = zcoef  * qsr(ji,jj)
204               zea(ji,jj,1) =          qsr(ji,jj)
205            END DO
206         END DO
207         !
[6403]208         DO jk = 2, nksr+1                   !* interior equi-partition in R-G-B depending of vertical profile of Chl
[6140]209            DO jj = 2, jpjm1
210               DO ji = fs_2, fs_jpim1
[6403]211                  zchl = MIN( 10. , MAX( 0.03, zchl3d(ji,jj,jk) ) )
212                  irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
213                  zekb(ji,jj) = rkrgb(1,irgb)
214                  zekg(ji,jj) = rkrgb(2,irgb)
215                  zekr(ji,jj) = rkrgb(3,irgb)
216               END DO
217            END DO
218
219            DO jj = 2, jpjm1
220               DO ji = fs_2, fs_jpim1
[6140]221                  zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * xsi0r       )
222                  zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekb(ji,jj) )
223                  zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekg(ji,jj) )
224                  zc3 = ze3(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekr(ji,jj) )
225                  ze0(ji,jj,jk) = zc0
226                  ze1(ji,jj,jk) = zc1
227                  ze2(ji,jj,jk) = zc2
228                  ze3(ji,jj,jk) = zc3
229                  zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk)
[1423]230               END DO
[6140]231            END DO
232         END DO
233         !
234         DO jk = 1, nksr                     !* now qsr induced heat content
235            DO jj = 2, jpjm1
236               DO ji = fs_2, fs_jpim1
237                  qsr_hc(ji,jj,jk) = r1_rau0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) )
[1423]238               END DO
[6140]239            END DO
240         END DO
[187]241         !
[6140]242         CALL wrk_dealloc( jpi,jpj,        zekb, zekg, zekr        ) 
[6403]243         CALL wrk_dealloc( jpi,jpj,jpk,   ze0, ze1, ze2, ze3, zea, zchl3d ) 
[6140]244         !
245      CASE( np_2BD  )            !==  2-bands fluxes  ==!
246         !
247         zz0 =        rn_abs   * r1_rau0_rcp      ! surface equi-partition in 2-bands
248         zz1 = ( 1. - rn_abs ) * r1_rau0_rcp
249         DO jk = 1, nksr                          ! solar heat absorbed at T-point in the top 400m
250            DO jj = 2, jpjm1
251               DO ji = fs_2, fs_jpim1
252                  zc0 = zz0 * EXP( -gdepw_n(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk  )*xsi1r )
253                  zc1 = zz0 * EXP( -gdepw_n(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk+1)*xsi1r )
254                  qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 
[2528]255               END DO
256            END DO
257         END DO
258         !
[6140]259      END SELECT
260      !
261      !                          !-----------------------------!
262      DO jk = 1, nksr            !  update to the temp. trend  !
263         DO jj = 2, jpjm1        !-----------------------------!
264            DO ji = fs_2, fs_jpim1   ! vector opt.
265               tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
266                  &                 + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t_n(ji,jj,jk)
267            END DO
268         END DO
269      END DO
270      !
271      IF( ln_qsr_ice ) THEN      ! sea-ice: store the 1st ocean level attenuation coefficient
272         DO jj = 2, jpjm1 
273            DO ji = fs_2, fs_jpim1   ! vector opt.
274               IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) )
275               ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp
276               ENDIF
277            END DO
278         END DO
279         ! Update haloes since lim_thd needs fraqsr_1lev to be defined everywhere
280         CALL lbc_lnk( fraqsr_1lev(:,:), 'T', 1._wp )
[3]281      ENDIF
[2528]282      !
[6140]283      IF( iom_use('qsr3d') ) THEN      ! output the shortwave Radiation distribution
284         CALL wrk_alloc( jpi,jpj,jpk,   zetot )
[2528]285         !
[7753]286         zetot(:,:,nksr+1:jpk) = 0._wp     ! below ~400m set to zero
[6140]287         DO jk = nksr, 1, -1
[7753]288            zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) / r1_rau0_rcp
[6140]289         END DO         
290         CALL iom_put( 'qsr3d', zetot )   ! 3D distribution of shortwave Radiation
291         !
292         CALL wrk_dealloc( jpi,jpj,jpk,   zetot ) 
[2528]293      ENDIF
[6140]294      !
295      IF( lrst_oce ) THEN     ! write in the ocean restart file
[8770]296         IF( lwxios ) CALL iom_swap(      cwxios_context          )
297         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc, ldxios = lwxios      )
298         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev, ldxios = lwxios ) 
299         IF( lwxios ) CALL iom_swap(      cxios_context          )
[6140]300      ENDIF
301      !
[503]302      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
[7753]303         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
[4990]304         CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt )
[6140]305         CALL wrk_dealloc( jpi,jpj,jpk,   ztrdt ) 
[3]306      ENDIF
[457]307      !                       ! print mean trends (used for debugging)
[2528]308      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
[503]309      !
[3294]310      IF( nn_timing == 1 )  CALL timing_stop('tra_qsr')
311      !
[3]312   END SUBROUTINE tra_qsr
313
314
315   SUBROUTINE tra_qsr_init
316      !!----------------------------------------------------------------------
317      !!                  ***  ROUTINE tra_qsr_init  ***
318      !!
319      !! ** Purpose :   Initialization for the penetrative solar radiation
320      !!
321      !! ** Method  :   The profile of solar radiation within the ocean is set
[1423]322      !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
[1601]323      !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
[3]324      !!      default values correspond to clear water (type I in Jerlov'
325      !!      (1968) classification.
326      !!         called by tra_qsr at the first timestep (nit000)
327      !!
[1423]328      !! ** Action  : - initialize rn_si0, rn_si1 and rn_abs
[3]329      !!
[503]330      !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
[3]331      !!----------------------------------------------------------------------
[6140]332      INTEGER  ::   ji, jj, jk                  ! dummy loop indices
333      INTEGER  ::   ios, irgb, ierror, ioptio   ! local integer
334      REAL(wp) ::   zz0, zc0 , zc1, zcoef      ! local scalars
335      REAL(wp) ::   zz1, zc2 , zc3, zchl       !   -      -
[2715]336      !
[1423]337      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files
338      TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read
[2715]339      !!
[6140]340      NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice,  &
[2528]341         &                  nn_chldta, rn_abs, rn_si0, rn_si1
[3]342      !!----------------------------------------------------------------------
[3294]343      !
[6140]344      IF( nn_timing == 1 )   CALL timing_start('tra_qsr_init')
[3294]345      !
[6140]346      REWIND( numnam_ref )              ! Namelist namtra_qsr in reference     namelist
347      READ  ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901)
348901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_qsr in reference namelist', lwp )
[3294]349      !
[6140]350      REWIND( numnam_cfg )              ! Namelist namtra_qsr in configuration namelist
[4147]351      READ  ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 )
[6140]352902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist', lwp )
[4624]353      IF(lwm) WRITE ( numond, namtra_qsr )
[1423]354      !
355      IF(lwp) THEN                ! control print
356         WRITE(numout,*)
357         WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
358         WRITE(numout,*) '~~~~~~~~~~~~'
[1601]359         WRITE(numout,*) '   Namelist namtra_qsr : set the parameter of penetration'
[6140]360         WRITE(numout,*) '      RGB (Red-Green-Blue) light penetration       ln_qsr_rgb = ', ln_qsr_rgb
361         WRITE(numout,*) '      2 band               light penetration       ln_qsr_2bd = ', ln_qsr_2bd
362         WRITE(numout,*) '      bio-model            light penetration       ln_qsr_bio = ', ln_qsr_bio
363         WRITE(numout,*) '      light penetration for ice-model (LIM3)       ln_qsr_ice = ', ln_qsr_ice
364         WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)        nn_chldta  = ', nn_chldta
365         WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs     = ', rn_abs
366         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0     = ', rn_si0
367         WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1     = ', rn_si1
368         WRITE(numout,*)
[1423]369      ENDIF
[6140]370      !
371      ioptio = 0                    ! Parameter control
372      IF( ln_qsr_rgb  )   ioptio = ioptio + 1
373      IF( ln_qsr_2bd  )   ioptio = ioptio + 1
374      IF( ln_qsr_bio  )   ioptio = ioptio + 1
375      !
376      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE type of light penetration in namelist namtra_qsr',  &
377         &                               ' 2 bands, 3 RGB bands or bio-model light penetration' )
378      !
379      IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr = np_RGB 
380      IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr = np_RGBc
381      IF( ln_qsr_2bd                      )   nqsr = np_2BD
382      IF( ln_qsr_bio                      )   nqsr = np_BIO
383      !
384      !                             ! Initialisation
385      xsi0r = 1._wp / rn_si0
386      xsi1r = 1._wp / rn_si1
387      !
388      SELECT CASE( nqsr )
389      !                               
390      CASE( np_RGB , np_RGBc )         !==  Red-Green-Blue light penetration  ==!
391         !                             
392         IF(lwp)   WRITE(numout,*) '   R-G-B   light penetration '
393         !
394         CALL trc_oce_rgb( rkrgb )                 ! tabulated attenuation coef.
395         !                                   
396         nksr = trc_oce_ext_lev( r_si2, 33._wp )   ! level of light extinction
397         !
398         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
399         !
400         IF( nqsr == np_RGBc ) THEN                ! Chl data : set sf_chl structure
401            IF(lwp) WRITE(numout,*) '        Chlorophyll read in a file'
402            ALLOCATE( sf_chl(1), STAT=ierror )
403            IF( ierror > 0 ) THEN
404               CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' )   ;   RETURN
405            ENDIF
406            ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1)   )
407            IF( sn_chl%ln_tint )   ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) )
408            !                                        ! fill sf_chl with sn_chl and control print
409            CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',   &
[7646]410               &           'Solar penetration function of read chlorophyll', 'namtra_qsr' , no_print )
[1448]411         ENDIF
[6140]412         IF( nqsr == np_RGB ) THEN                 ! constant Chl
413            IF(lwp) WRITE(numout,*) '        Constant Chlorophyll concentration = 0.05'
414         ENDIF
[1448]415         !
[6140]416      CASE( np_2BD )                   !==  2 bands light penetration  ==!
[1448]417         !
[6140]418         IF(lwp)  WRITE(numout,*) '   2 bands light penetration'
[1448]419         !
[6140]420         nksr = trc_oce_ext_lev( rn_si1, 100._wp )    ! level of light extinction
421         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
[1455]422         !
[6140]423      CASE( np_BIO )                   !==  BIO light penetration  ==!
[1448]424         !
[6140]425         IF(lwp) WRITE(numout,*) '   bio-model light penetration'
[7646]426         IF( .NOT.lk_top )   CALL ctl_stop( 'No bio model : ln_qsr_bio = true impossible ' )
[1423]427         !
[6140]428      END SELECT
[503]429      !
[7753]430      qsr_hc(:,:,:) = 0._wp     ! now qsr heat content set to zero where it will not be computed
[6140]431      !
432      ! 1st ocean level attenuation coefficient (used in sbcssm)
[5407]433      IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN
[8668]434         CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev'  , fraqsr_1lev, ldxios = lxios_read  )
[5407]435      ELSE
[7753]436         fraqsr_1lev(:,:) = 1._wp   ! default : no penetration
[5407]437      ENDIF
438      !
[8770]439      IF( lwxios ) THEN
440         CALL iom_set_rstw_var_active('qsr_hc_b')
441         CALL iom_set_rstw_var_active('fraqsr_1lev')
442      ENDIF
[6140]443      IF( nn_timing == 1 )   CALL timing_stop('tra_qsr_init')
[2715]444      !
[3]445   END SUBROUTINE tra_qsr_init
446
447   !!======================================================================
448END MODULE traqsr
Note: See TracBrowser for help on using the repository browser.