New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
traqsr.F90 in NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traqsr.F90 @ 10946

Last change on this file since 10946 was 10946, checked in by acc, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert STO, TRD and USR modules and all knock on effects of these conversions. Note change to USR module may have implications for the TEST CASES (not tested yet). Standard SETTE tested only

  • Property svn:keywords set to Id
File size: 22.0 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
[9019]31   USE iom            ! I/O library
[6140]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 timing         ! Timing
[3]37
38   IMPLICIT NONE
39   PRIVATE
40
[2528]41   PUBLIC   tra_qsr       ! routine called by step.F90 (ln_traqsr=T)
[5407]42   PUBLIC   tra_qsr_init  ! routine called by nemogcm.F90
[3]43
[4147]44   !                                 !!* Namelist namtra_qsr: penetrative solar radiation
45   LOGICAL , PUBLIC ::   ln_traqsr    !: light absorption (qsr) flag
46   LOGICAL , PUBLIC ::   ln_qsr_rgb   !: Red-Green-Blue light absorption flag 
47   LOGICAL , PUBLIC ::   ln_qsr_2bd   !: 2 band         light absorption flag
48   LOGICAL , PUBLIC ::   ln_qsr_bio   !: bio-model      light absorption flag
49   INTEGER , PUBLIC ::   nn_chldta    !: use Chlorophyll data (=1) or not (=0)
50   REAL(wp), PUBLIC ::   rn_abs       !: fraction absorbed in the very near surface (RGB & 2 bands)
51   REAL(wp), PUBLIC ::   rn_si0       !: very near surface depth of extinction      (RGB & 2 bands)
52   REAL(wp), PUBLIC ::   rn_si1       !: deepest depth of extinction (water type I)       (2 bands)
[6140]53   !
54   INTEGER , PUBLIC ::   nksr         !: levels below which the light cannot penetrate (depth larger than 391 m)
[5407]55 
[6140]56   INTEGER, PARAMETER ::   np_RGB  = 1   ! R-G-B     light penetration with constant Chlorophyll
57   INTEGER, PARAMETER ::   np_RGBc = 2   ! R-G-B     light penetration with Chlorophyll data
58   INTEGER, PARAMETER ::   np_2BD  = 3   ! 2 bands   light penetration
59   INTEGER, PARAMETER ::   np_BIO  = 4   ! bio-model light penetration
60   !
61   INTEGER  ::   nqsr    ! user choice of the type of light penetration
62   REAL(wp) ::   xsi0r   ! inverse of rn_si0
63   REAL(wp) ::   xsi1r   ! inverse of rn_si1
64   !
65   REAL(wp) , DIMENSION(3,61)           ::   rkrgb    ! tabulated attenuation coefficients for RGB absorption
[1423]66   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read)
[3]67
68   !! * Substitutions
69#  include "vectopt_loop_substitute.h90"
70   !!----------------------------------------------------------------------
[9598]71   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[888]72   !! $Id$
[10068]73   !! Software governed by the CeCILL license (see ./LICENSE)
[3]74   !!----------------------------------------------------------------------
75CONTAINS
76
[10946]77   SUBROUTINE tra_qsr( kt, Kmm, Krhs )
[3]78      !!----------------------------------------------------------------------
79      !!                  ***  ROUTINE tra_qsr  ***
80      !!
81      !! ** Purpose :   Compute the temperature trend due to the solar radiation
[6140]82      !!              penetration and add it to the general temperature trend.
[3]83      !!
[1423]84      !! ** Method  : The profile of the solar radiation within the ocean is defined
85      !!      through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) and a ratio rn_abs
86      !!      Considering the 2 wavebands case:
87      !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
88      !!         The temperature trend associated with the solar radiation penetration
89      !!         is given by : zta = 1/e3t dk[ I ] / (rau0*Cp)
[3]90      !!         At the bottom, boudary condition for the radiation is no flux :
91      !!      all heat which has not been absorbed in the above levels is put
92      !!      in the last ocean level.
[6140]93      !!         The computation is only done down to the level where
94      !!      I(k) < 1.e-15 W/m2 (i.e. over the top nksr levels) .
[3]95      !!
96      !! ** Action  : - update ta with the penetrative solar radiation trend
[6140]97      !!              - send  trend for further diagnostics (l_trdtra=T)
[1423]98      !!
99      !! Reference  : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
100      !!              Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
[6403]101      !!              Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562
[503]102      !!----------------------------------------------------------------------
103      INTEGER, INTENT(in) ::   kt     ! ocean time-step
[10946]104      INTEGER, INTENT(in) ::   Kmm, Krhs     ! time level indices
[2715]105      !
[6140]106      INTEGER  ::   ji, jj, jk               ! dummy loop indices
107      INTEGER  ::   irgb                     ! local integers
108      REAL(wp) ::   zchl, zcoef, z1_2        ! local scalars
109      REAL(wp) ::   zc0 , zc1 , zc2 , zc3    !    -         -
[4161]110      REAL(wp) ::   zzc0, zzc1, zzc2, zzc3   !    -         -
[6140]111      REAL(wp) ::   zz0 , zz1                !    -         -
[6403]112      REAL(wp) ::   zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze
113      REAL(wp) ::   zlogc, zlogc2, zlogc3 
[9019]114      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   :: zekb, zekg, zekr
115      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt
116      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zetot, zchl3d
[3]117      !!----------------------------------------------------------------------
[3294]118      !
[9019]119      IF( ln_timing )   CALL timing_start('tra_qsr')
[3294]120      !
[3]121      IF( kt == nit000 ) THEN
[503]122         IF(lwp) WRITE(numout,*)
123         IF(lwp) WRITE(numout,*) 'tra_qsr : penetration of the surface solar radiation'
124         IF(lwp) WRITE(numout,*) '~~~~~~~'
[3]125      ENDIF
[6140]126      !
127      IF( l_trdtra ) THEN      ! trends diagnostic: save the input temperature trend
[9019]128         ALLOCATE( ztrdt(jpi,jpj,jpk) ) 
[10874]129         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
[216]130      ENDIF
[6140]131      !
132      !                         !-----------------------------------!
133      !                         !  before qsr induced heat content  !
134      !                         !-----------------------------------!
135      IF( kt == nit000 ) THEN          !==  1st time step  ==!
136!!gm case neuler  not taken into account....
137         IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN    ! read in restart
138            IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field read in the restart file'
139            z1_2 = 0.5_wp
[9367]140            CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b, ldxios = lrxios )   ! before heat content trend due to Qsr flux
[2528]141         ELSE                                           ! No restart or restart not found: Euler forward time stepping
[6140]142            z1_2 = 1._wp
[7753]143            qsr_hc_b(:,:,:) = 0._wp
[2528]144         ENDIF
[6140]145      ELSE                             !==  Swap of qsr heat content  ==!
146         z1_2 = 0.5_wp
[7753]147         qsr_hc_b(:,:,:) = qsr_hc(:,:,:)
[2528]148      ENDIF
[6140]149      !
150      !                         !--------------------------------!
151      SELECT CASE( nqsr )       !  now qsr induced heat content  !
152      !                         !--------------------------------!
153      !
154      CASE( np_BIO )                   !==  bio-model fluxes  ==!
155         !
156         DO jk = 1, nksr
[7753]157            qsr_hc(:,:,jk) = r1_rau0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) )
[2528]158         END DO
[6140]159         !
160      CASE( np_RGB , np_RGBc )         !==  R-G-B fluxes  ==!
161         !
[9019]162         ALLOCATE( zekb(jpi,jpj)     , zekg(jpi,jpj)     , zekr  (jpi,jpj)     , &
163            &      ze0 (jpi,jpj,jpk) , ze1 (jpi,jpj,jpk) , ze2   (jpi,jpj,jpk) , &
164            &      ze3 (jpi,jpj,jpk) , zea (jpi,jpj,jpk) , zchl3d(jpi,jpj,jpk)   ) 
[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)
[10874]175                     zpsi    = gdepw_n(ji,jj,jk) / zze
[6403]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
[10874]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) )
[6140]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         !
[9019]242         DEALLOCATE( zekb , zekg , zekr , ze0 , ze1 , ze2 , ze3 , zea , zchl3d ) 
[6140]243         !
244      CASE( np_2BD  )            !==  2-bands fluxes  ==!
245         !
246         zz0 =        rn_abs   * r1_rau0_rcp      ! surface equi-partition in 2-bands
247         zz1 = ( 1. - rn_abs ) * r1_rau0_rcp
248         DO jk = 1, nksr                          ! solar heat absorbed at T-point in the top 400m
249            DO jj = 2, jpjm1
250               DO ji = fs_2, fs_jpim1
[10874]251                  zc0 = zz0 * EXP( -gdepw_n(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk  )*xsi1r )
252                  zc1 = zz0 * EXP( -gdepw_n(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk+1)*xsi1r )
[6140]253                  qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 
[2528]254               END DO
255            END DO
256         END DO
257         !
[6140]258      END SELECT
259      !
260      !                          !-----------------------------!
261      DO jk = 1, nksr            !  update to the temp. trend  !
262         DO jj = 2, jpjm1        !-----------------------------!
263            DO ji = fs_2, fs_jpim1   ! vector opt.
[10874]264               tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
265                  &                 + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t_n(ji,jj,jk)
[6140]266            END DO
267         END DO
268      END DO
269      !
[9019]270      ! sea-ice: store the 1st ocean level attenuation coefficient
271      DO jj = 2, jpjm1 
272         DO ji = fs_2, fs_jpim1   ! vector opt.
273            IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) )
274            ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp
275            ENDIF
[6140]276         END DO
[9019]277      END DO
[10425]278      CALL lbc_lnk( 'traqsr', fraqsr_1lev(:,:), 'T', 1._wp )
[2528]279      !
[6140]280      IF( iom_use('qsr3d') ) THEN      ! output the shortwave Radiation distribution
[9019]281         ALLOCATE( zetot(jpi,jpj,jpk) )
[7753]282         zetot(:,:,nksr+1:jpk) = 0._wp     ! below ~400m set to zero
[6140]283         DO jk = nksr, 1, -1
[9019]284            zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rau0_rcp
[6140]285         END DO         
286         CALL iom_put( 'qsr3d', zetot )   ! 3D distribution of shortwave Radiation
[9019]287         DEALLOCATE( zetot ) 
[2528]288      ENDIF
[6140]289      !
290      IF( lrst_oce ) THEN     ! write in the ocean restart file
[9367]291         IF( lwxios ) CALL iom_swap(      cwxios_context          )
292         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc     , ldxios = lwxios )
293         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev, ldxios = lwxios ) 
294         IF( lwxios ) CALL iom_swap(      cxios_context          )
[6140]295      ENDIF
296      !
[503]297      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
[10874]298         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
[10946]299         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt )
[9019]300         DEALLOCATE( ztrdt ) 
[3]301      ENDIF
[457]302      !                       ! print mean trends (used for debugging)
[2528]303      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
[503]304      !
[9019]305      IF( ln_timing )   CALL timing_stop('tra_qsr')
[3294]306      !
[3]307   END SUBROUTINE tra_qsr
308
309
310   SUBROUTINE tra_qsr_init
311      !!----------------------------------------------------------------------
312      !!                  ***  ROUTINE tra_qsr_init  ***
313      !!
314      !! ** Purpose :   Initialization for the penetrative solar radiation
315      !!
316      !! ** Method  :   The profile of solar radiation within the ocean is set
[1423]317      !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
[1601]318      !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
[3]319      !!      default values correspond to clear water (type I in Jerlov'
320      !!      (1968) classification.
321      !!         called by tra_qsr at the first timestep (nit000)
322      !!
[1423]323      !! ** Action  : - initialize rn_si0, rn_si1 and rn_abs
[3]324      !!
[503]325      !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
[3]326      !!----------------------------------------------------------------------
[6140]327      INTEGER  ::   ji, jj, jk                  ! dummy loop indices
328      INTEGER  ::   ios, irgb, ierror, ioptio   ! local integer
329      REAL(wp) ::   zz0, zc0 , zc1, zcoef      ! local scalars
330      REAL(wp) ::   zz1, zc2 , zc3, zchl       !   -      -
[2715]331      !
[1423]332      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files
333      TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read
[2715]334      !!
[9019]335      NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio,  &
[2528]336         &                  nn_chldta, rn_abs, rn_si0, rn_si1
[3]337      !!----------------------------------------------------------------------
[3294]338      !
[6140]339      REWIND( numnam_ref )              ! Namelist namtra_qsr in reference     namelist
340      READ  ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901)
341901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_qsr in reference namelist', lwp )
[3294]342      !
[6140]343      REWIND( numnam_cfg )              ! Namelist namtra_qsr in configuration namelist
[4147]344      READ  ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 )
[9168]345902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist', lwp )
[4624]346      IF(lwm) WRITE ( numond, namtra_qsr )
[1423]347      !
348      IF(lwp) THEN                ! control print
349         WRITE(numout,*)
350         WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
351         WRITE(numout,*) '~~~~~~~~~~~~'
[1601]352         WRITE(numout,*) '   Namelist namtra_qsr : set the parameter of penetration'
[6140]353         WRITE(numout,*) '      RGB (Red-Green-Blue) light penetration       ln_qsr_rgb = ', ln_qsr_rgb
354         WRITE(numout,*) '      2 band               light penetration       ln_qsr_2bd = ', ln_qsr_2bd
355         WRITE(numout,*) '      bio-model            light penetration       ln_qsr_bio = ', ln_qsr_bio
356         WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)        nn_chldta  = ', nn_chldta
357         WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs     = ', rn_abs
358         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0     = ', rn_si0
359         WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1     = ', rn_si1
360         WRITE(numout,*)
[1423]361      ENDIF
[6140]362      !
363      ioptio = 0                    ! Parameter control
364      IF( ln_qsr_rgb  )   ioptio = ioptio + 1
365      IF( ln_qsr_2bd  )   ioptio = ioptio + 1
366      IF( ln_qsr_bio  )   ioptio = ioptio + 1
367      !
368      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE type of light penetration in namelist namtra_qsr',  &
369         &                               ' 2 bands, 3 RGB bands or bio-model light penetration' )
370      !
371      IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr = np_RGB 
372      IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr = np_RGBc
373      IF( ln_qsr_2bd                      )   nqsr = np_2BD
374      IF( ln_qsr_bio                      )   nqsr = np_BIO
375      !
376      !                             ! Initialisation
377      xsi0r = 1._wp / rn_si0
378      xsi1r = 1._wp / rn_si1
379      !
380      SELECT CASE( nqsr )
381      !                               
382      CASE( np_RGB , np_RGBc )         !==  Red-Green-Blue light penetration  ==!
383         !                             
[9169]384         IF(lwp)   WRITE(numout,*) '   ==>>>   R-G-B   light penetration '
[6140]385         !
386         CALL trc_oce_rgb( rkrgb )                 ! tabulated attenuation coef.
387         !                                   
388         nksr = trc_oce_ext_lev( r_si2, 33._wp )   ! level of light extinction
389         !
390         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
391         !
392         IF( nqsr == np_RGBc ) THEN                ! Chl data : set sf_chl structure
[9169]393            IF(lwp) WRITE(numout,*) '   ==>>>   Chlorophyll read in a file'
[6140]394            ALLOCATE( sf_chl(1), STAT=ierror )
395            IF( ierror > 0 ) THEN
396               CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' )   ;   RETURN
397            ENDIF
398            ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1)   )
399            IF( sn_chl%ln_tint )   ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) )
400            !                                        ! fill sf_chl with sn_chl and control print
401            CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',   &
[7646]402               &           'Solar penetration function of read chlorophyll', 'namtra_qsr' , no_print )
[1448]403         ENDIF
[6140]404         IF( nqsr == np_RGB ) THEN                 ! constant Chl
[9169]405            IF(lwp) WRITE(numout,*) '   ==>>>   Constant Chlorophyll concentration = 0.05'
[6140]406         ENDIF
[1448]407         !
[6140]408      CASE( np_2BD )                   !==  2 bands light penetration  ==!
[1448]409         !
[9169]410         IF(lwp)  WRITE(numout,*) '   ==>>>   2 bands light penetration'
[1448]411         !
[6140]412         nksr = trc_oce_ext_lev( rn_si1, 100._wp )    ! level of light extinction
413         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
[1455]414         !
[6140]415      CASE( np_BIO )                   !==  BIO light penetration  ==!
[1448]416         !
[9169]417         IF(lwp) WRITE(numout,*) '   ==>>>   bio-model light penetration'
[7646]418         IF( .NOT.lk_top )   CALL ctl_stop( 'No bio model : ln_qsr_bio = true impossible ' )
[1423]419         !
[6140]420      END SELECT
[503]421      !
[7753]422      qsr_hc(:,:,:) = 0._wp     ! now qsr heat content set to zero where it will not be computed
[6140]423      !
424      ! 1st ocean level attenuation coefficient (used in sbcssm)
[5407]425      IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN
[9367]426         CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev'  , fraqsr_1lev, ldxios = lrxios  )
[5407]427      ELSE
[7753]428         fraqsr_1lev(:,:) = 1._wp   ! default : no penetration
[5407]429      ENDIF
430      !
[9367]431      IF( lwxios ) THEN
432         CALL iom_set_rstw_var_active('qsr_hc_b')
433         CALL iom_set_rstw_var_active('fraqsr_1lev')
434      ENDIF
435      !
[3]436   END SUBROUTINE tra_qsr_init
437
438   !!======================================================================
439END MODULE traqsr
Note: See TracBrowser for help on using the repository browser.