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.
Changeset 6225 for branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90 – NEMO

Ignore:
Timestamp:
2016-01-08T10:35:19+01:00 (8 years ago)
Author:
jamesharle
Message:

Update MPP_BDY_UPDATE branch to be consistent with head of trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

    r4624 r6225  
    22   !!====================================================================== 
    33   !!                       ***  MODULE  traqsr  *** 
    4    !! Ocean physics: solar radiation penetration in the top ocean levels 
     4   !! Ocean physics:   solar radiation penetration in the top ocean levels 
    55   !!====================================================================== 
    66   !! History :  OPA  !  1990-10  (B. Blanke)  Original code 
     
    1010   !!             -   !  2005-11  (G. Madec) zco, zps, sco coordinate 
    1111   !!            3.2  !  2009-04  (G. Madec & NEMO team)  
    12    !!            4.0  !  2012-05  (C. Rousset) store attenuation coef for use in ice model  
     12   !!            3.6  !  2012-05  (C. Rousset) store attenuation coef for use in ice model  
     13   !!            3.7  !  2015-11  (G. Madec, A. Coward)  remove optimisation for fix volume  
    1314   !!---------------------------------------------------------------------- 
    1415 
    1516   !!---------------------------------------------------------------------- 
    16    !!   tra_qsr      : trend due to the solar radiation penetration 
    17    !!   tra_qsr_init : solar radiation penetration initialization 
     17   !!   tra_qsr       : temperature trend due to the penetration of solar radiation  
     18   !!   tra_qsr_init  : initialization of the qsr penetration  
    1819   !!---------------------------------------------------------------------- 
    19    USE oce             ! ocean dynamics and active tracers 
    20    USE dom_oce         ! ocean space and time domain 
    21    USE sbc_oce         ! surface boundary condition: ocean 
    22    USE trc_oce         ! share SMS/Ocean variables 
    23    USE trdmod_oce      ! ocean variables trends 
    24    USE trdtra          ! ocean active tracers trends  
    25    USE in_out_manager  ! I/O manager 
    26    USE phycst          ! physical constants 
    27    USE prtctl          ! Print control 
    28    USE iom             ! I/O manager 
    29    USE fldread         ! read input fields 
    30    USE restart         ! ocean restart 
    31    USE lib_mpp         ! MPP library 
     20   USE oce            ! ocean dynamics and active tracers 
     21   USE phycst         ! physical constants 
     22   USE dom_oce        ! ocean space and time domain 
     23   USE sbc_oce        ! surface boundary condition: ocean 
     24   USE trc_oce        ! share SMS/Ocean variables 
     25   USE trd_oce        ! trends: ocean variables 
     26   USE trdtra         ! trends manager: tracers 
     27   ! 
     28   USE in_out_manager ! I/O manager 
     29   USE prtctl         ! Print control 
     30   USE iom            ! I/O manager 
     31   USE fldread        ! read input fields 
     32   USE restart        ! ocean restart 
     33   USE lib_mpp        ! MPP library 
     34   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    3235   USE wrk_nemo       ! Memory Allocation 
    3336   USE timing         ! Timing 
    34    USE sbc_ice, ONLY : lk_lim3 
    3537 
    3638   IMPLICIT NONE 
     
    3840 
    3941   PUBLIC   tra_qsr       ! routine called by step.F90 (ln_traqsr=T) 
    40    PUBLIC   tra_qsr_init  ! routine called by opa.F90 
     42   PUBLIC   tra_qsr_init  ! routine called by nemogcm.F90 
    4143 
    4244   !                                 !!* Namelist namtra_qsr: penetrative solar radiation 
     
    5052   REAL(wp), PUBLIC ::   rn_si0       !: very near surface depth of extinction      (RGB & 2 bands) 
    5153   REAL(wp), PUBLIC ::   rn_si1       !: deepest depth of extinction (water type I)       (2 bands) 
    52     
    53    ! Module variables 
    54    REAL(wp) ::   xsi0r                           !: inverse of rn_si0 
    55    REAL(wp) ::   xsi1r                           !: inverse of rn_si1 
     54   ! 
     55   INTEGER , PUBLIC ::   nksr         !: levels below which the light cannot penetrate (depth larger than 391 m) 
     56  
     57   INTEGER, PARAMETER ::   np_RGB  = 1   ! R-G-B     light penetration with constant Chlorophyll 
     58   INTEGER, PARAMETER ::   np_RGBc = 2   ! R-G-B     light penetration with Chlorophyll data 
     59   INTEGER, PARAMETER ::   np_2BD  = 3   ! 2 bands   light penetration 
     60   INTEGER, PARAMETER ::   np_BIO  = 4   ! bio-model light penetration 
     61   ! 
     62   INTEGER  ::   nqsr    ! user choice of the type of light penetration 
     63   REAL(wp) ::   xsi0r   ! inverse of rn_si0 
     64   REAL(wp) ::   xsi1r   ! inverse of rn_si1 
     65   ! 
     66   REAL(wp) , DIMENSION(3,61)           ::   rkrgb    ! tabulated attenuation coefficients for RGB absorption 
    5667   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read) 
    57    INTEGER, PUBLIC ::   nksr              ! levels below which the light cannot penetrate ( depth larger than 391 m) 
    58    REAL(wp), DIMENSION(3,61) ::   rkrgb   !: tabulated attenuation coefficients for RGB absorption 
    5968 
    6069   !! * Substitutions 
    61 #  include "domzgr_substitute.h90" 
    6270#  include "vectopt_loop_substitute.h90" 
    6371   !!---------------------------------------------------------------------- 
     
    7381      !! 
    7482      !! ** Purpose :   Compute the temperature trend due to the solar radiation 
    75       !!      penetration and add it to the general temperature trend. 
     83      !!              penetration and add it to the general temperature trend. 
    7684      !! 
    7785      !! ** Method  : The profile of the solar radiation within the ocean is defined 
     
    8492      !!      all heat which has not been absorbed in the above levels is put 
    8593      !!      in the last ocean level. 
    86       !!         In z-coordinate case, the computation is only done down to the 
    87       !!      level where I(k) < 1.e-15 W/m2. In addition, the coefficients  
    88       !!      used for the computation are calculated one for once as they 
    89       !!      depends on k only. 
     94      !!         The computation is only done down to the level where  
     95      !!      I(k) < 1.e-15 W/m2 (i.e. over the top nksr levels) .  
    9096      !! 
    9197      !! ** Action  : - update ta with the penetrative solar radiation trend 
    92       !!              - save the trend in ttrd ('key_trdtra') 
     98      !!              - send  trend for further diagnostics (l_trdtra=T) 
    9399      !! 
    94100      !! Reference  : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp. 
    95101      !!              Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516. 
    96102      !!---------------------------------------------------------------------- 
    97       ! 
    98103      INTEGER, INTENT(in) ::   kt     ! ocean time-step 
    99104      ! 
    100       INTEGER  ::   ji, jj, jk           ! dummy loop indices 
    101       INTEGER  ::   irgb                 ! local integers 
    102       REAL(wp) ::   zchl, zcoef, zfact   ! local scalars 
    103       REAL(wp) ::   zc0, zc1, zc2, zc3   !    -         - 
     105      INTEGER  ::   ji, jj, jk               ! dummy loop indices 
     106      INTEGER  ::   irgb                     ! local integers 
     107      REAL(wp) ::   zchl, zcoef, z1_2        ! local scalars 
     108      REAL(wp) ::   zc0 , zc1 , zc2 , zc3    !    -         - 
    104109      REAL(wp) ::   zzc0, zzc1, zzc2, zzc3   !    -         - 
    105       REAL(wp) ::   zz0, zz1, z1_e3t     !    -         - 
    106       REAL(wp), POINTER, DIMENSION(:,:  ) :: zekb, zekg, zekr 
     110      REAL(wp) ::   zz0 , zz1                !    -         - 
     111      REAL(wp), POINTER, DIMENSION(:,: :: zekb, zekg, zekr 
    107112      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt 
     113      REAL(wp), POINTER, DIMENSION(:,:,:) :: zetot 
    108114      !!---------------------------------------------------------------------- 
    109115      ! 
    110116      IF( nn_timing == 1 )  CALL timing_start('tra_qsr') 
    111       ! 
    112       CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        )  
    113       CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )  
    114117      ! 
    115118      IF( kt == nit000 ) THEN 
     
    117120         IF(lwp) WRITE(numout,*) 'tra_qsr : penetration of the surface solar radiation' 
    118121         IF(lwp) WRITE(numout,*) '~~~~~~~' 
    119          IF( .NOT.ln_traqsr )   RETURN 
    120       ENDIF 
    121  
    122       IF( l_trdtra ) THEN      ! Save ta and sa trends 
    123          CALL wrk_alloc( jpi, jpj, jpk, ztrdt )  
     122      ENDIF 
     123      ! 
     124      IF( l_trdtra ) THEN      ! trends diagnostic: save the input temperature trend 
     125         CALL wrk_alloc( jpi,jpj,jpk,   ztrdt )  
    124126         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
    125127      ENDIF 
    126  
    127       !                                        Set before qsr tracer content field 
    128       !                                        *********************************** 
    129       IF( kt == nit000 ) THEN                     ! Set the forcing field at nit000 - 1 
    130          !                                        ! ----------------------------------- 
    131          IF( ln_rstart .AND.    &                    ! Restart: read in restart file 
    132               & iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN 
    133             IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field red in the restart file' 
    134             zfact = 0.5e0 
     128      ! 
     129      !                         !-----------------------------------! 
     130      !                         !  before qsr induced heat content  ! 
     131      !                         !-----------------------------------! 
     132      IF( kt == nit000 ) THEN          !==  1st time step  ==! 
     133!!gm case neuler  not taken into account.... 
     134         IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN    ! read in restart 
     135            IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field read in the restart file' 
     136            z1_2 = 0.5_wp 
    135137            CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b )   ! before heat content trend due to Qsr flux 
    136138         ELSE                                           ! No restart or restart not found: Euler forward time stepping 
    137             zfact = 1.e0 
    138             qsr_hc_b(:,:,:) = 0.e0 
     139            z1_2 = 1._wp 
     140            qsr_hc_b(:,:,:) = 0._wp 
    139141         ENDIF 
    140       ELSE                                        ! Swap of forcing field 
    141          !                                        ! --------------------- 
    142          zfact = 0.5e0 
     142      ELSE                             !==  Swap of qsr heat content  ==! 
     143         z1_2 = 0.5_wp 
    143144         qsr_hc_b(:,:,:) = qsr_hc(:,:,:) 
    144145      ENDIF 
    145       !                                        Compute now qsr tracer content field 
    146       !                                        ************************************ 
    147        
    148       !                                           ! ============================================== ! 
    149       IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN      !  bio-model fluxes  : all vertical coordinates  ! 
    150          !                                        ! ============================================== ! 
    151          DO jk = 1, jpkm1 
     146      ! 
     147      !                         !--------------------------------! 
     148      SELECT CASE( nqsr )       !  now qsr induced heat content  ! 
     149      !                         !--------------------------------! 
     150      ! 
     151      CASE( np_BIO )                   !==  bio-model fluxes  ==! 
     152         ! 
     153         DO jk = 1, nksr 
    152154            qsr_hc(:,:,jk) = r1_rau0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) ) 
    153155         END DO 
    154          !                                        Add to the general trend 
    155          DO jk = 1, jpkm1 
    156             DO jj = 2, jpjm1  
    157                DO ji = fs_2, fs_jpim1   ! vector opt. 
    158                   z1_e3t = zfact / fse3t(ji,jj,jk) 
    159                   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 
     156         ! 
     157      CASE( np_RGB , np_RGBc )         !==  R-G-B fluxes  ==! 
     158         ! 
     159         CALL wrk_alloc( jpi,jpj,       zekb, zekg, zekr        )  
     160         CALL wrk_alloc( jpi,jpj,jpk,   ze0, ze1, ze2, ze3, zea )  
     161         ! 
     162         IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll 
     163            CALL fld_read( kt, 1, sf_chl )         ! Read Chl data and provides it at the current time step 
     164            DO jj = 2, jpjm1                       ! Separation in R-G-B depending of the surface Chl 
     165               DO ji = fs_2, fs_jpim1 
     166                  zchl = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) 
     167                  irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 
     168                  zekb(ji,jj) = rkrgb(1,irgb) 
     169                  zekg(ji,jj) = rkrgb(2,irgb) 
     170                  zekr(ji,jj) = rkrgb(3,irgb) 
    160171               END DO 
    161172            END DO 
    162          END DO 
    163          CALL iom_put( 'qsr3d', etot3 )   ! Shortwave Radiation 3D distribution 
    164          ! clem: store attenuation coefficient of the first ocean level 
    165          IF ( lk_lim3 .AND. ln_qsr_ice ) THEN 
    166             DO jj = 1, jpj 
    167                DO ji = 1, jpi 
    168                   IF ( qsr(ji,jj) /= 0._wp ) THEN 
    169                      oatte(ji,jj) = ( qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) ) ) 
    170                      iatte(ji,jj) = oatte(ji,jj) 
    171                   ENDIF 
     173         ELSE                                !* constant chrlorophyll 
     174            zchl = 0.05                            ! constant chlorophyll 
     175            !                                      ! Separation in R-G-B depending of the chlorophyll 
     176            irgb = NINT( 41 + 20.*LOG10( zchl ) + 1.e-15 ) 
     177            DO jj = 2, jpjm1 
     178               DO ji = fs_2, fs_jpim1 
     179                  zekb(ji,jj) = rkrgb(1,irgb)                       
     180                  zekg(ji,jj) = rkrgb(2,irgb) 
     181                  zekr(ji,jj) = rkrgb(3,irgb) 
    172182               END DO 
    173183            END DO 
    174184         ENDIF 
    175          !                                        ! ============================================== ! 
    176       ELSE                                        !  Ocean alone :  
    177          !                                        ! ============================================== ! 
    178          ! 
    179          !                                                ! ------------------------- ! 
    180          IF( ln_qsr_rgb) THEN                             !  R-G-B  light penetration ! 
    181             !                                             ! ------------------------- ! 
    182             ! Set chlorophyl concentration 
    183             IF( nn_chldta == 1 .OR. lk_vvl ) THEN            !*  Variable Chlorophyll or ocean volume 
    184                ! 
    185                IF( nn_chldta == 1 ) THEN                             !*  Variable Chlorophyll 
    186                   ! 
    187                   CALL fld_read( kt, 1, sf_chl )                         ! Read Chl data and provides it at the current time step 
    188                   !          
    189 !CDIR COLLAPSE 
    190 !CDIR NOVERRCHK 
    191                   DO jj = 1, jpj                                         ! Separation in R-G-B depending of the surface Chl 
    192 !CDIR NOVERRCHK 
    193                      DO ji = 1, jpi 
    194                         zchl = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) 
    195                         irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 
    196                         zekb(ji,jj) = rkrgb(1,irgb) 
    197                         zekg(ji,jj) = rkrgb(2,irgb) 
    198                         zekr(ji,jj) = rkrgb(3,irgb) 
    199                      END DO 
    200                   END DO 
    201                ELSE                                            ! Variable ocean volume but constant chrlorophyll 
    202                   zchl = 0.05                                     ! constant chlorophyll 
    203                   irgb = NINT( 41 + 20.*LOG10( zchl ) + 1.e-15 ) 
    204                   zekb(:,:) = rkrgb(1,irgb)                       ! Separation in R-G-B depending of the chlorophyll  
    205                   zekg(:,:) = rkrgb(2,irgb) 
    206                   zekr(:,:) = rkrgb(3,irgb) 
     185         ! 
     186         zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B 
     187         DO jj = 2, jpjm1 
     188            DO ji = fs_2, fs_jpim1 
     189               ze0(ji,jj,1) = rn_abs * qsr(ji,jj) 
     190               ze1(ji,jj,1) = zcoef  * qsr(ji,jj) 
     191               ze2(ji,jj,1) = zcoef  * qsr(ji,jj) 
     192               ze3(ji,jj,1) = zcoef  * qsr(ji,jj) 
     193               zea(ji,jj,1) =          qsr(ji,jj) 
     194            END DO 
     195         END DO 
     196         ! 
     197         DO jk = 2, nksr+1                   !* interior equi-partition in R-G-B 
     198            DO jj = 2, jpjm1 
     199               DO ji = fs_2, fs_jpim1 
     200                  zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * xsi0r       ) 
     201                  zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekb(ji,jj) ) 
     202                  zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekg(ji,jj) ) 
     203                  zc3 = ze3(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekr(ji,jj) ) 
     204                  ze0(ji,jj,jk) = zc0 
     205                  ze1(ji,jj,jk) = zc1 
     206                  ze2(ji,jj,jk) = zc2 
     207                  ze3(ji,jj,jk) = zc3 
     208                  zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk) 
     209               END DO 
     210            END DO 
     211         END DO 
     212         ! 
     213         DO jk = 1, nksr                     !* now qsr induced heat content 
     214            DO jj = 2, jpjm1 
     215               DO ji = fs_2, fs_jpim1 
     216                  qsr_hc(ji,jj,jk) = r1_rau0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) ) 
     217               END DO 
     218            END DO 
     219         END DO 
     220         ! 
     221         CALL wrk_dealloc( jpi,jpj,        zekb, zekg, zekr        )  
     222         CALL wrk_dealloc( jpi,jpj,jpk,   ze0, ze1, ze2, ze3, zea )  
     223         ! 
     224      CASE( np_2BD  )            !==  2-bands fluxes  ==! 
     225         ! 
     226         zz0 =        rn_abs   * r1_rau0_rcp      ! surface equi-partition in 2-bands 
     227         zz1 = ( 1. - rn_abs ) * r1_rau0_rcp 
     228         DO jk = 1, nksr                          ! solar heat absorbed at T-point in the top 400m  
     229            DO jj = 2, jpjm1 
     230               DO ji = fs_2, fs_jpim1 
     231                  zc0 = zz0 * EXP( -gdepw_n(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk  )*xsi1r ) 
     232                  zc1 = zz0 * EXP( -gdepw_n(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk+1)*xsi1r ) 
     233                  qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) )  
     234               END DO 
     235            END DO 
     236         END DO 
     237         ! 
     238      END SELECT 
     239      ! 
     240      !                          !-----------------------------! 
     241      DO jk = 1, nksr            !  update to the temp. trend  ! 
     242         DO jj = 2, jpjm1        !-----------------------------! 
     243            DO ji = fs_2, fs_jpim1   ! vector opt. 
     244               tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   & 
     245                  &                 + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t_n(ji,jj,jk) 
     246            END DO 
     247         END DO 
     248      END DO 
     249      ! 
     250      IF( ln_qsr_ice ) THEN      ! sea-ice: store the 1st ocean level attenuation coefficient 
     251         DO jj = 2, jpjm1  
     252            DO ji = fs_2, fs_jpim1   ! vector opt. 
     253               IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) ) 
     254               ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp 
    207255               ENDIF 
    208                ! 
    209                zcoef  = ( 1. - rn_abs ) / 3.e0                        ! equi-partition in R-G-B 
    210                ze0(:,:,1) = rn_abs  * qsr(:,:) 
    211                ze1(:,:,1) = zcoef * qsr(:,:) 
    212                ze2(:,:,1) = zcoef * qsr(:,:) 
    213                ze3(:,:,1) = zcoef * qsr(:,:) 
    214                zea(:,:,1) =         qsr(:,:) 
    215                ! 
    216                DO jk = 2, nksr+1 
    217 !CDIR NOVERRCHK 
    218                   DO jj = 1, jpj 
    219 !CDIR NOVERRCHK    
    220                      DO ji = 1, jpi 
    221                         zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r     ) 
    222                         zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekb(ji,jj) ) 
    223                         zc2 = ze2(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekg(ji,jj) ) 
    224                         zc3 = ze3(ji,jj,jk-1) * EXP( - fse3t(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 ) * tmask(ji,jj,jk) 
    230                      END DO 
    231                   END DO 
    232                END DO 
    233                ! clem: store attenuation coefficient of the first ocean level 
    234                IF ( lk_lim3 .AND. ln_qsr_ice ) THEN 
    235                   DO jj = 1, jpj 
    236                      DO ji = 1, jpi 
    237                         zzc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r     ) 
    238                         zzc1 = zcoef  * EXP( - fse3t(ji,jj,1) * zekb(ji,jj) ) 
    239                         zzc2 = zcoef  * EXP( - fse3t(ji,jj,1) * zekg(ji,jj) ) 
    240                         zzc3 = zcoef  * EXP( - fse3t(ji,jj,1) * zekr(ji,jj) ) 
    241                         oatte(ji,jj) = 1.0 - ( zzc0 + zzc1 + zzc2  + zzc3  ) * tmask(ji,jj,2)  
    242                         iatte(ji,jj) = 1.0 - ( zzc0 + zzc1 + zcoef + zcoef ) * tmask(ji,jj,2) 
    243                      END DO 
    244                   END DO 
    245                ENDIF 
    246                ! 
    247                DO jk = 1, nksr                                        ! compute and add qsr trend to ta 
    248                   qsr_hc(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) ) 
    249                END DO 
    250                zea(:,:,nksr+1:jpk) = 0.e0     ! below 400m set to zero 
    251                CALL iom_put( 'qsr3d', zea )   ! Shortwave Radiation 3D distribution 
    252                ! 
    253             ELSE                                                 !*  Constant Chlorophyll 
    254                DO jk = 1, nksr 
    255                   qsr_hc(:,:,jk) =  etot3(:,:,jk) * qsr(:,:) 
    256                END DO 
    257                ! clem: store attenuation coefficient of the first ocean level 
    258                IF ( lk_lim3 .AND. ln_qsr_ice ) THEN 
    259                   oatte(:,:) = etot3(:,:,1) / r1_rau0_rcp 
    260                   iatte(:,:) = oatte(:,:) 
    261                ENDIF 
    262            ENDIF 
    263  
    264          ENDIF 
    265          !                                                ! ------------------------- ! 
    266          IF( ln_qsr_2bd ) THEN                            !  2 band light penetration ! 
    267             !                                             ! ------------------------- ! 
    268             ! 
    269             IF( lk_vvl ) THEN                                  !* variable volume 
    270                zz0   =        rn_abs   * r1_rau0_rcp 
    271                zz1   = ( 1. - rn_abs ) * r1_rau0_rcp 
    272                DO jk = 1, nksr                    ! solar heat absorbed at T-point in the top 400m  
    273                   DO jj = 1, jpj 
    274                      DO ji = 1, jpi 
    275                         zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r ) 
    276                         zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 
    277                         qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) )  
    278                      END DO 
    279                   END DO 
    280                END DO 
    281                ! clem: store attenuation coefficient of the first ocean level 
    282                IF ( lk_lim3 .AND. ln_qsr_ice ) THEN 
    283                   DO jj = 1, jpj 
    284                      DO ji = 1, jpi 
    285                         zc0 = zz0 * EXP( -fsdepw(ji,jj,1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,1)*xsi1r ) 
    286                         zc1 = zz0 * EXP( -fsdepw(ji,jj,2)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,2)*xsi1r ) 
    287                         oatte(ji,jj) = ( zc0*tmask(ji,jj,1) - zc1*tmask(ji,jj,2) ) / r1_rau0_rcp 
    288                         iatte(ji,jj) = oatte(ji,jj) 
    289                      END DO 
    290                   END DO 
    291                ENDIF 
    292             ELSE                                               !* constant volume: coef. computed one for all 
    293                DO jk = 1, nksr 
    294                   DO jj = 2, jpjm1 
    295                      DO ji = fs_2, fs_jpim1   ! vector opt. 
    296                         qsr_hc(ji,jj,jk) =  etot3(ji,jj,jk) * qsr(ji,jj) 
    297                      END DO 
    298                   END DO 
    299                END DO 
    300                ! clem: store attenuation coefficient of the first ocean level 
    301                IF ( lk_lim3 .AND. ln_qsr_ice ) THEN 
    302                   oatte(:,:) = etot3(:,:,1) / r1_rau0_rcp 
    303                   iatte(:,:) = oatte(:,:) 
    304                ENDIF 
    305                ! 
    306             ENDIF 
    307             ! 
    308          ENDIF 
    309          ! 
    310          !                                        Add to the general trend 
    311          DO jk = 1, nksr 
    312             DO jj = 2, jpjm1  
    313                DO ji = fs_2, fs_jpim1   ! vector opt. 
    314                   z1_e3t = zfact / fse3t(ji,jj,jk) 
    315                   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 
    316                END DO 
    317             END DO 
    318          END DO 
    319          ! 
    320       ENDIF 
    321       ! 
    322       IF( lrst_oce ) THEN   !                  Write in the ocean restart file 
    323          !                                     ******************************* 
    324          IF(lwp) WRITE(numout,*) 
    325          IF(lwp) WRITE(numout,*) 'qsr tracer content forcing field written in ocean restart file ',   & 
    326             &                    'at it= ', kt,' date= ', ndastp 
    327          IF(lwp) WRITE(numout,*) '~~~~' 
    328          CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b', qsr_hc ) 
    329          ! 
    330       ENDIF 
    331  
     256            END DO 
     257         END DO 
     258         ! Update haloes since lim_thd needs fraqsr_1lev to be defined everywhere 
     259         CALL lbc_lnk( fraqsr_1lev(:,:), 'T', 1._wp ) 
     260      ENDIF 
     261      ! 
     262      IF( iom_use('qsr3d') ) THEN      ! output the shortwave Radiation distribution 
     263         CALL wrk_alloc( jpi,jpj,jpk,   zetot ) 
     264         ! 
     265         zetot(:,:,nksr+1:jpk) = 0._wp     ! below ~400m set to zero 
     266         DO jk = nksr, 1, -1 
     267            zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) / r1_rau0_rcp 
     268         END DO          
     269         CALL iom_put( 'qsr3d', zetot )   ! 3D distribution of shortwave Radiation 
     270         ! 
     271         CALL wrk_dealloc( jpi,jpj,jpk,   zetot )  
     272      ENDIF 
     273      ! 
     274      IF( lrst_oce ) THEN     ! write in the ocean restart file 
     275         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc      ) 
     276         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev )  
     277      ENDIF 
     278      ! 
    332279      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics 
    333280         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    334          CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_qsr, ztrdt ) 
    335          CALL wrk_dealloc( jpi, jpj, jpk, ztrdt )  
     281         CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt ) 
     282         CALL wrk_dealloc( jpi,jpj,jpk,  ztrdt )  
    336283      ENDIF 
    337284      !                       ! print mean trends (used for debugging) 
    338285      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' ) 
    339       ! 
    340       CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        )  
    341       CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )  
    342286      ! 
    343287      IF( nn_timing == 1 )  CALL timing_stop('tra_qsr') 
     
    363307      !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp. 
    364308      !!---------------------------------------------------------------------- 
    365       ! 
    366       INTEGER  ::   ji, jj, jk                   ! dummy loop indices 
    367       INTEGER  ::   irgb, ierror, ioptio, nqsr   ! local integer 
    368       INTEGER  ::   ios                          ! Local integer output status for namelist read 
    369       REAL(wp) ::   zz0, zc0  , zc1, zcoef       ! local scalars 
    370       REAL(wp) ::   zz1, zc2  , zc3, zchl        !   -      - 
    371       REAL(wp), POINTER, DIMENSION(:,:  ) :: zekb, zekg, zekr 
    372       REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea 
     309      INTEGER  ::   ji, jj, jk                  ! dummy loop indices 
     310      INTEGER  ::   ios, irgb, ierror, ioptio   ! local integer 
     311      REAL(wp) ::   zz0, zc0 , zc1, zcoef      ! local scalars 
     312      REAL(wp) ::   zz1, zc2 , zc3, zchl       !   -      - 
    373313      ! 
    374314      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files 
    375315      TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read 
    376316      !! 
    377       NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice,  & 
     317      NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice,  & 
    378318         &                  nn_chldta, rn_abs, rn_si0, rn_si1 
    379319      !!---------------------------------------------------------------------- 
    380  
    381       ! 
    382       IF( nn_timing == 1 )  CALL timing_start('tra_qsr_init') 
    383       ! 
    384       ! clem init for oatte and iatte 
    385       IF( .NOT. ln_rstart ) THEN 
    386          oatte(:,:) = 1._wp 
    387          iatte(:,:) = 1._wp 
    388       ENDIF 
    389       ! 
    390       CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        )  
    391       CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )  
    392       ! 
    393  
    394       REWIND( numnam_ref )              ! Namelist namtra_qsr in reference namelist : Ratio and length of penetration 
     320      ! 
     321      IF( nn_timing == 1 )   CALL timing_start('tra_qsr_init') 
     322      ! 
     323      REWIND( numnam_ref )              ! Namelist namtra_qsr in reference     namelist 
    395324      READ  ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901) 
    396 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in reference namelist', lwp ) 
    397  
    398       REWIND( numnam_cfg )              !  Namelist namtra_qsr in configuration namelist : Ratio and length of penetration 
     325901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_qsr in reference namelist', lwp ) 
     326      ! 
     327      REWIND( numnam_cfg )              ! Namelist namtra_qsr in configuration namelist 
    399328      READ  ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 ) 
    400 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist', lwp ) 
     329902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist', lwp ) 
    401330      IF(lwm) WRITE ( numond, namtra_qsr ) 
    402331      ! 
     
    406335         WRITE(numout,*) '~~~~~~~~~~~~' 
    407336         WRITE(numout,*) '   Namelist namtra_qsr : set the parameter of penetration' 
    408          WRITE(numout,*) '      Light penetration (T) or not (F)         ln_traqsr  = ', ln_traqsr 
    409          WRITE(numout,*) '      RGB (Red-Green-Blue) light penetration   ln_qsr_rgb = ', ln_qsr_rgb 
    410          WRITE(numout,*) '      2 band               light penetration   ln_qsr_2bd = ', ln_qsr_2bd 
    411          WRITE(numout,*) '      bio-model            light penetration   ln_qsr_bio = ', ln_qsr_bio 
    412          WRITE(numout,*) '      light penetration for ice-model LIM3     ln_qsr_ice = ', ln_qsr_ice 
    413          WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)    nn_chldta  = ', nn_chldta 
    414          WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs = ', rn_abs 
    415          WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0 = ', rn_si0 
    416          WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1 = ', rn_si1 
    417          WRITE(numout,*) '      light penetration for ice-model LIM3     ln_qsr_ice = ', ln_qsr_ice     
    418       ENDIF 
    419  
    420       IF( ln_traqsr ) THEN     ! control consistency 
    421          !                       
    422          IF( .NOT.lk_qsr_bio .AND. ln_qsr_bio )   THEN 
    423             CALL ctl_warn( 'No bio model : force ln_qsr_bio = FALSE ' ) 
    424             ln_qsr_bio = .FALSE. 
     337         WRITE(numout,*) '      RGB (Red-Green-Blue) light penetration       ln_qsr_rgb = ', ln_qsr_rgb 
     338         WRITE(numout,*) '      2 band               light penetration       ln_qsr_2bd = ', ln_qsr_2bd 
     339         WRITE(numout,*) '      bio-model            light penetration       ln_qsr_bio = ', ln_qsr_bio 
     340         WRITE(numout,*) '      light penetration for ice-model (LIM3)       ln_qsr_ice = ', ln_qsr_ice 
     341         WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)        nn_chldta  = ', nn_chldta 
     342         WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs     = ', rn_abs 
     343         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0     = ', rn_si0 
     344         WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1     = ', rn_si1 
     345         WRITE(numout,*) 
     346      ENDIF 
     347      ! 
     348      ioptio = 0                    ! Parameter control 
     349      IF( ln_qsr_rgb  )   ioptio = ioptio + 1 
     350      IF( ln_qsr_2bd  )   ioptio = ioptio + 1 
     351      IF( ln_qsr_bio  )   ioptio = ioptio + 1 
     352      ! 
     353      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE type of light penetration in namelist namtra_qsr',  & 
     354         &                               ' 2 bands, 3 RGB bands or bio-model light penetration' ) 
     355      ! 
     356      IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr = np_RGB  
     357      IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr = np_RGBc 
     358      IF( ln_qsr_2bd                      )   nqsr = np_2BD 
     359      IF( ln_qsr_bio                      )   nqsr = np_BIO 
     360      ! 
     361      !                             ! Initialisation 
     362      xsi0r = 1._wp / rn_si0 
     363      xsi1r = 1._wp / rn_si1 
     364      ! 
     365      SELECT CASE( nqsr ) 
     366      !                                
     367      CASE( np_RGB , np_RGBc )         !==  Red-Green-Blue light penetration  ==! 
     368         !                              
     369         IF(lwp)   WRITE(numout,*) '   R-G-B   light penetration ' 
     370         ! 
     371         CALL trc_oce_rgb( rkrgb )                 ! tabulated attenuation coef. 
     372         !                                    
     373         nksr = trc_oce_ext_lev( r_si2, 33._wp )   ! level of light extinction 
     374         ! 
     375         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m' 
     376         ! 
     377         IF( nqsr == np_RGBc ) THEN                ! Chl data : set sf_chl structure 
     378            IF(lwp) WRITE(numout,*) '        Chlorophyll read in a file' 
     379            ALLOCATE( sf_chl(1), STAT=ierror ) 
     380            IF( ierror > 0 ) THEN 
     381               CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' )   ;   RETURN 
     382            ENDIF 
     383            ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1)   ) 
     384            IF( sn_chl%ln_tint )   ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) ) 
     385            !                                        ! fill sf_chl with sn_chl and control print 
     386            CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',   & 
     387               &           'Solar penetration function of read chlorophyll', 'namtra_qsr' ) 
    425388         ENDIF 
    426          ! 
    427          ioptio = 0                      ! Parameter control 
    428          IF( ln_qsr_rgb  )   ioptio = ioptio + 1 
    429          IF( ln_qsr_2bd  )   ioptio = ioptio + 1 
    430          IF( ln_qsr_bio  )   ioptio = ioptio + 1 
    431          ! 
    432          IF( ioptio /= 1 ) & 
    433             CALL ctl_stop( '          Choose ONE type of light penetration in namelist namtra_qsr',  & 
    434             &              ' 2 bands, 3 RGB bands or bio-model light penetration' ) 
    435          ! 
    436          IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr =  1  
    437          IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr =  2 
    438          IF( ln_qsr_2bd                      )   nqsr =  3 
    439          IF( ln_qsr_bio                      )   nqsr =  4 
    440          ! 
    441          IF(lwp) THEN                   ! Print the choice 
    442             WRITE(numout,*) 
    443             IF( nqsr ==  1 )   WRITE(numout,*) '         R-G-B   light penetration - Constant Chlorophyll' 
    444             IF( nqsr ==  2 )   WRITE(numout,*) '         R-G-B   light penetration - Chl data ' 
    445             IF( nqsr ==  3 )   WRITE(numout,*) '         2 bands light penetration' 
    446             IF( nqsr ==  4 )   WRITE(numout,*) '         bio-model light penetration' 
     389         IF( nqsr == np_RGB ) THEN                 ! constant Chl 
     390            IF(lwp) WRITE(numout,*) '        Constant Chlorophyll concentration = 0.05' 
    447391         ENDIF 
    448392         ! 
    449       ENDIF 
    450       !                          ! ===================================== ! 
    451       IF( ln_traqsr  ) THEN      !  Initialisation of Light Penetration  !   
    452          !                       ! ===================================== ! 
    453          ! 
    454          xsi0r = 1.e0 / rn_si0 
    455          xsi1r = 1.e0 / rn_si1 
    456          !                                ! ---------------------------------- ! 
    457          IF( ln_qsr_rgb ) THEN            !  Red-Green-Blue light penetration  ! 
    458             !                             ! ---------------------------------- ! 
    459             ! 
    460             CALL trc_oce_rgb( rkrgb )           !* tabulated attenuation coef. 
    461             ! 
    462             !                                   !* level of light extinction 
    463             IF(  ln_sco ) THEN   ;   nksr = jpkm1 
    464             ELSE                 ;   nksr = trc_oce_ext_lev( r_si2, 0.33e2 ) 
    465             ENDIF 
    466  
    467             IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m' 
    468             ! 
    469             IF( nn_chldta == 1 ) THEN           !* Chl data : set sf_chl structure 
    470                IF(lwp) WRITE(numout,*) 
    471                IF(lwp) WRITE(numout,*) '        Chlorophyll read in a file' 
    472                ALLOCATE( sf_chl(1), STAT=ierror ) 
    473                IF( ierror > 0 ) THEN 
    474                   CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' )   ;   RETURN 
    475                ENDIF 
    476                ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1)   ) 
    477                IF( sn_chl%ln_tint )ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) ) 
    478                !                                        ! fill sf_chl with sn_chl and control print 
    479                CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',   & 
    480                   &                                         'Solar penetration function of read chlorophyll', 'namtra_qsr' ) 
    481                ! 
    482             ELSE                                !* constant Chl : compute once for all the distribution of light (etot3) 
    483                IF(lwp) WRITE(numout,*) 
    484                IF(lwp) WRITE(numout,*) '        Constant Chlorophyll concentration = 0.05' 
    485                IF( lk_vvl ) THEN                   ! variable volume 
    486                   IF(lwp) WRITE(numout,*) '        key_vvl: light distribution will be computed at each time step' 
    487                ELSE                                ! constant volume: computes one for all 
    488                   IF(lwp) WRITE(numout,*) '        fixed volume: light distribution computed one for all' 
    489                   ! 
    490                   zchl = 0.05                                 ! constant chlorophyll 
    491                   irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 
    492                   zekb(:,:) = rkrgb(1,irgb)                   ! Separation in R-G-B depending of the chlorophyll  
    493                   zekg(:,:) = rkrgb(2,irgb) 
    494                   zekr(:,:) = rkrgb(3,irgb) 
    495                   ! 
    496                   zcoef = ( 1. - rn_abs ) / 3.e0              ! equi-partition in R-G-B 
    497                   ze0(:,:,1) = rn_abs 
    498                   ze1(:,:,1) = zcoef 
    499                   ze2(:,:,1) = zcoef  
    500                   ze3(:,:,1) = zcoef 
    501                   zea(:,:,1) = tmask(:,:,1)                   ! = ( ze0+ze1+z2+ze3 ) * tmask 
    502                 
    503                   DO jk = 2, nksr+1 
    504 !CDIR NOVERRCHK 
    505                      DO jj = 1, jpj 
    506 !CDIR NOVERRCHK    
    507                         DO ji = 1, jpi 
    508                            zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r     ) 
    509                            zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekb(ji,jj) ) 
    510                            zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekg(ji,jj) ) 
    511                            zc3 = ze3(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekr(ji,jj) ) 
    512                            ze0(ji,jj,jk) = zc0 
    513                            ze1(ji,jj,jk) = zc1 
    514                            ze2(ji,jj,jk) = zc2 
    515                            ze3(ji,jj,jk) = zc3 
    516                            zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk) 
    517                         END DO 
    518                      END DO 
    519                   END DO  
    520                   ! 
    521                   DO jk = 1, nksr 
    522                      etot3(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) )  
    523                   END DO 
    524                   etot3(:,:,nksr+1:jpk) = 0.e0                ! below 400m set to zero 
    525                ENDIF 
    526             ENDIF 
    527             ! 
    528          ENDIF 
    529             !                             ! ---------------------------------- ! 
    530          IF( ln_qsr_2bd ) THEN            !    2 bands    light penetration    ! 
    531             !                             ! ---------------------------------- ! 
    532             ! 
    533             !                                ! level of light extinction 
    534             nksr = trc_oce_ext_lev( rn_si1, 1.e2 ) 
    535             IF(lwp) THEN 
    536                WRITE(numout,*) 
    537             IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m' 
    538             ENDIF 
    539             ! 
    540             IF( lk_vvl ) THEN                   ! variable volume 
    541                IF(lwp) WRITE(numout,*) '        key_vvl: light distribution will be computed at each time step' 
    542             ELSE                                ! constant volume: computes one for all 
    543                zz0 =        rn_abs   * r1_rau0_rcp 
    544                zz1 = ( 1. - rn_abs ) * r1_rau0_rcp 
    545                DO jk = 1, nksr                    !*  solar heat absorbed at T-point computed once for all 
    546                   DO jj = 1, jpj                              ! top 400 meters 
    547                      DO ji = 1, jpi 
    548                         zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r ) 
    549                         zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 
    550                         etot3(ji,jj,jk) = (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  )  
    551                      END DO 
    552                   END DO 
    553                END DO 
    554                etot3(:,:,nksr+1:jpk) = 0.e0                   ! below 400m set to zero 
    555                ! 
    556             ENDIF 
    557          ENDIF 
    558          !                       ! ===================================== ! 
    559       ELSE                       !        No light penetration           !                    
    560          !                       ! ===================================== ! 
    561          IF(lwp) THEN 
    562             WRITE(numout,*) 
    563             WRITE(numout,*) 'tra_qsr_init : NO solar flux penetration' 
    564             WRITE(numout,*) '~~~~~~~~~~~~' 
    565          ENDIF 
    566       ENDIF 
    567       ! 
    568       CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        )  
    569       CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )  
    570       ! 
    571       IF( nn_timing == 1 )  CALL timing_stop('tra_qsr_init') 
     393      CASE( np_2BD )                   !==  2 bands light penetration  ==! 
     394         ! 
     395         IF(lwp)  WRITE(numout,*) '   2 bands light penetration' 
     396         ! 
     397         nksr = trc_oce_ext_lev( rn_si1, 100._wp )    ! level of light extinction 
     398         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m' 
     399         ! 
     400      CASE( np_BIO )                   !==  BIO light penetration  ==! 
     401         ! 
     402         IF(lwp) WRITE(numout,*) '   bio-model light penetration' 
     403         IF( .NOT.lk_qsr_bio )   CALL ctl_stop( 'No bio model : ln_qsr_bio = true impossible ' ) 
     404         ! 
     405      END SELECT 
     406      ! 
     407      qsr_hc(:,:,:) = 0._wp     ! now qsr heat content set to zero where it will not be computed 
     408      ! 
     409      ! 1st ocean level attenuation coefficient (used in sbcssm) 
     410      IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN 
     411         CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev'  , fraqsr_1lev  ) 
     412      ELSE 
     413         fraqsr_1lev(:,:) = 1._wp   ! default : no penetration 
     414      ENDIF 
     415      ! 
     416      IF( nn_timing == 1 )   CALL timing_stop('tra_qsr_init') 
    572417      ! 
    573418   END SUBROUTINE tra_qsr_init 
Note: See TracChangeset for help on using the changeset viewer.