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/trunk/src/OCE/TRA – NEMO

source: NEMO/trunk/src/OCE/TRA/traqsr.F90 @ 13333

Last change on this file since 13333 was 13333, checked in by cetlod, 4 years ago

trunk : bugfix in the heat content trend when activating bio-model light penetration, see ticket #2499

  • Property svn:keywords set to Id
File size: 22.9 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   !
[13333]65   REAL(wp) , PUBLIC, 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
[12377]69#  include "do_loop_substitute.h90"
[13237]70#  include "domzgr_substitute.h90"
[3]71   !!----------------------------------------------------------------------
[9598]72   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[888]73   !! $Id$
[10068]74   !! Software governed by the CeCILL license (see ./LICENSE)
[3]75   !!----------------------------------------------------------------------
76CONTAINS
77
[12377]78   SUBROUTINE tra_qsr( kt, Kmm, pts, Krhs )
[3]79      !!----------------------------------------------------------------------
80      !!                  ***  ROUTINE tra_qsr  ***
81      !!
82      !! ** Purpose :   Compute the temperature trend due to the solar radiation
[6140]83      !!              penetration and add it to the general temperature trend.
[3]84      !!
[1423]85      !! ** Method  : The profile of the solar radiation within the ocean is defined
86      !!      through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) and a ratio rn_abs
87      !!      Considering the 2 wavebands case:
88      !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
89      !!         The temperature trend associated with the solar radiation penetration
[12489]90      !!         is given by : zta = 1/e3t dk[ I ] / (rho0*Cp)
[3]91      !!         At the bottom, boudary condition for the radiation is no flux :
92      !!      all heat which has not been absorbed in the above levels is put
93      !!      in the last ocean level.
[6140]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) .
[3]96      !!
97      !! ** Action  : - update ta with the penetrative solar radiation trend
[6140]98      !!              - send  trend for further diagnostics (l_trdtra=T)
[1423]99      !!
100      !! Reference  : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
101      !!              Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
[6403]102      !!              Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562
[503]103      !!----------------------------------------------------------------------
[12377]104      INTEGER,                                   INTENT(in   ) :: kt            ! ocean time-step
105      INTEGER,                                   INTENT(in   ) :: Kmm, Krhs     ! time level indices
106      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts           ! active tracers and RHS of tracer equation
[2715]107      !
[6140]108      INTEGER  ::   ji, jj, jk               ! dummy loop indices
109      INTEGER  ::   irgb                     ! local integers
110      REAL(wp) ::   zchl, zcoef, z1_2        ! local scalars
111      REAL(wp) ::   zc0 , zc1 , zc2 , zc3    !    -         -
[4161]112      REAL(wp) ::   zzc0, zzc1, zzc2, zzc3   !    -         -
[13205]113      REAL(wp) ::   zz0 , zz1 , ze3t, zlui   !    -         -
114      REAL(wp) ::   zCb, zCmax, zpsi, zpsimax, zrdpsi, zCze
115      REAL(wp) ::   zlogc, zlogze, zlogCtot, zlogCze
116      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   :: ze0, ze1, ze2, ze3
117      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, zetot, ztmp3d
[3]118      !!----------------------------------------------------------------------
[3294]119      !
[9019]120      IF( ln_timing )   CALL timing_start('tra_qsr')
[3294]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
[9019]129         ALLOCATE( ztrdt(jpi,jpj,jpk) ) 
[12377]130         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
[216]131      ENDIF
[6140]132      !
133      !                         !-----------------------------------!
134      !                         !  before qsr induced heat content  !
135      !                         !-----------------------------------!
136      IF( kt == nit000 ) THEN          !==  1st time step  ==!
[12489]137         IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0  .AND. .NOT.l_1st_euler ) THEN    ! read in restart
[6140]138            IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field read in the restart file'
139            z1_2 = 0.5_wp
[13286]140            CALL iom_get( numror, jpdom_auto, '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
[12489]157            qsr_hc(:,:,jk) = r1_rho0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) )
[2528]158         END DO
[6140]159         !
160      CASE( np_RGB , np_RGBc )         !==  R-G-B fluxes  ==!
161         !
[13205]162         ALLOCATE( ze0 (jpi,jpj)           , ze1 (jpi,jpj) ,   &
163            &      ze2 (jpi,jpj)           , ze3 (jpi,jpj) ,   &
164            &      ztmp3d(jpi,jpj,nksr + 1)                     )
[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
[13205]168            !
169            ! Separation in R-G-B depending on the surface Chl
170            ! perform and store as many of the 2D calculations as possible
171            ! before the 3D loop (use the temporary 2D arrays to replace the
172            ! most expensive calculations)
173            !
[13295]174            DO_2D( 0, 0, 0, 0 )
[13205]175                       ! zlogc = log(zchl)
176               zlogc = LOG ( MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) )     
177                       ! zc1 : log(zCze)  = log (1.12  * zchl**0.803)
178               zc1   = 0.113328685307 + 0.803 * zlogc                         
179                       ! zc2 : log(zCtot) = log(40.6  * zchl**0.459)
180               zc2   = 3.703768066608 + 0.459 * zlogc                           
181                       ! zc3 : log(zze)   = log(568.2 * zCtot**(-0.746))
182               zc3   = 6.34247346942  - 0.746 * zc2                           
183                       ! IF( log(zze) > log(102.) ) log(zze) = log(200.0 * zCtot**(-0.293))
184               IF( zc3 > 4.62497281328 ) zc3 = 5.298317366548 - 0.293 * zc2       
185               !   
186               ze0(ji,jj) = zlogc                                                 ! ze0 = log(zchl)
187               ze1(ji,jj) = EXP( zc1 )                                            ! ze1 = zCze
188               ze2(ji,jj) = 1._wp / ( 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) ) ! ze2 = 1/zdelpsi
189               ze3(ji,jj) = EXP( - zc3 )                                          ! ze3 = 1/zze
190            END_2D
191           
192!
[13295]193            DO_3D( 0, 0, 0, 0, 1, nksr + 1 )
[13205]194               ! zchl    = ALOG( ze0(ji,jj) )
195               zlogc = ze0(ji,jj)
196               !
197               zCb       = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) )
198               zCmax     = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 )
199               zpsimax   = 0.6   - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) )
200               ! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 )
201               !
202               zCze   = ze1(ji,jj)
203               zrdpsi = ze2(ji,jj)                                                 ! 1/zdelpsi
204               zpsi   = ze3(ji,jj) * gdepw(ji,jj,jk,Kmm)                           ! gdepw/zze
205               !
206               ! NB. make sure zchl value is such that: zchl = MIN( 10. , MAX( 0.03, zchl ) )
207               zchl = MIN( 10. , MAX( 0.03, zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) ) )
208               ! Convert chlorophyll value to attenuation coefficient look-up table index
209               ztmp3d(ji,jj,jk) = 41 + 20.*LOG10(zchl) + 1.e-15
210            END_3D
211         ELSE                                !* constant chlorophyll
212            zchl = 0.05
213            ! NB. make sure constant value is such that:
214            zchl = MIN( 10. , MAX( 0.03, zchl ) )
215            ! Convert chlorophyll value to attenuation coefficient look-up table index
216            zlui = 41 + 20.*LOG10(zchl) + 1.e-15
[6403]217            DO jk = 1, nksr + 1
[13205]218               ztmp3d(:,:,jk) = zlui 
[3]219            END DO
[4161]220         ENDIF
[1423]221         !
[6140]222         zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B
[13295]223         DO_2D( 0, 0, 0, 0 )
[13205]224            ze0(ji,jj) = rn_abs * qsr(ji,jj)
225            ze1(ji,jj) = zcoef  * qsr(ji,jj)
226            ze2(ji,jj) = zcoef  * qsr(ji,jj)
227            ze3(ji,jj) = zcoef  * qsr(ji,jj)
228            ! store the surface SW radiation; re-use the surface ztmp3d array
229            ! since the surface attenuation coefficient is not used
230            ztmp3d(ji,jj,1) =       qsr(ji,jj)
[12377]231         END_2D
[6140]232         !
[13205]233         !* interior equi-partition in R-G-B depending on vertical profile of Chl
[13295]234         DO_3D( 0, 0, 0, 0, 2, nksr + 1 )
[13205]235            ze3t = e3t(ji,jj,jk-1,Kmm)
236            irgb = NINT( ztmp3d(ji,jj,jk) )
237            zc0 = ze0(ji,jj) * EXP( - ze3t * xsi0r )
238            zc1 = ze1(ji,jj) * EXP( - ze3t * rkrgb(1,irgb) )
239            zc2 = ze2(ji,jj) * EXP( - ze3t * rkrgb(2,irgb) )
240            zc3 = ze3(ji,jj) * EXP( - ze3t * rkrgb(3,irgb) )
241            ze0(ji,jj) = zc0
242            ze1(ji,jj) = zc1
243            ze2(ji,jj) = zc2
244            ze3(ji,jj) = zc3
245            ztmp3d(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk)
246         END_3D
[6140]247         !
[13295]248         DO_3D( 0, 0, 0, 0, 1, nksr )
[13205]249            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( ztmp3d(ji,jj,jk) - ztmp3d(ji,jj,jk+1) )
[12377]250         END_3D
[187]251         !
[13205]252         DEALLOCATE( ze0 , ze1 , ze2 , ze3 , ztmp3d ) 
[6140]253         !
254      CASE( np_2BD  )            !==  2-bands fluxes  ==!
255         !
[12489]256         zz0 =        rn_abs   * r1_rho0_rcp      ! surface equi-partition in 2-bands
257         zz1 = ( 1. - rn_abs ) * r1_rho0_rcp
[13295]258         DO_3D( 0, 0, 0, 0, 1, nksr )
[12377]259            zc0 = zz0 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi1r )
260            zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r )
261            qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 
262         END_3D
[2528]263         !
[6140]264      END SELECT
265      !
266      !                          !-----------------------------!
[13295]267      DO_3D( 0, 0, 0, 0, 1, nksr )
[12377]268         pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
[13237]269            &                      + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) )   &
270            &                             / e3t(ji,jj,jk,Kmm)
[12377]271      END_3D
[6140]272      !
[9019]273      ! sea-ice: store the 1st ocean level attenuation coefficient
[13295]274      DO_2D( 0, 0, 0, 0 )
[12489]275         IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) )
[12377]276         ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp
277         ENDIF
278      END_2D
[10425]279      CALL lbc_lnk( 'traqsr', fraqsr_1lev(:,:), 'T', 1._wp )
[2528]280      !
[6140]281      IF( iom_use('qsr3d') ) THEN      ! output the shortwave Radiation distribution
[9019]282         ALLOCATE( zetot(jpi,jpj,jpk) )
[7753]283         zetot(:,:,nksr+1:jpk) = 0._wp     ! below ~400m set to zero
[6140]284         DO jk = nksr, 1, -1
[12489]285            zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rho0_rcp
[6140]286         END DO         
287         CALL iom_put( 'qsr3d', zetot )   ! 3D distribution of shortwave Radiation
[9019]288         DEALLOCATE( zetot ) 
[2528]289      ENDIF
[6140]290      !
291      IF( lrst_oce ) THEN     ! write in the ocean restart file
[9367]292         IF( lwxios ) CALL iom_swap(      cwxios_context          )
293         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc     , ldxios = lwxios )
294         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev, ldxios = lwxios ) 
295         IF( lwxios ) CALL iom_swap(      cxios_context          )
[6140]296      ENDIF
297      !
[503]298      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
[12377]299         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
300         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt )
[9019]301         DEALLOCATE( ztrdt ) 
[3]302      ENDIF
[457]303      !                       ! print mean trends (used for debugging)
[12377]304      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
[503]305      !
[9019]306      IF( ln_timing )   CALL timing_stop('tra_qsr')
[3294]307      !
[3]308   END SUBROUTINE tra_qsr
309
310
311   SUBROUTINE tra_qsr_init
312      !!----------------------------------------------------------------------
313      !!                  ***  ROUTINE tra_qsr_init  ***
314      !!
315      !! ** Purpose :   Initialization for the penetrative solar radiation
316      !!
317      !! ** Method  :   The profile of solar radiation within the ocean is set
[1423]318      !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
[1601]319      !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
[3]320      !!      default values correspond to clear water (type I in Jerlov'
321      !!      (1968) classification.
322      !!         called by tra_qsr at the first timestep (nit000)
323      !!
[1423]324      !! ** Action  : - initialize rn_si0, rn_si1 and rn_abs
[3]325      !!
[503]326      !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
[3]327      !!----------------------------------------------------------------------
[6140]328      INTEGER  ::   ji, jj, jk                  ! dummy loop indices
329      INTEGER  ::   ios, irgb, ierror, ioptio   ! local integer
330      REAL(wp) ::   zz0, zc0 , zc1, zcoef      ! local scalars
331      REAL(wp) ::   zz1, zc2 , zc3, zchl       !   -      -
[2715]332      !
[1423]333      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files
334      TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read
[2715]335      !!
[9019]336      NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio,  &
[2528]337         &                  nn_chldta, rn_abs, rn_si0, rn_si1
[3]338      !!----------------------------------------------------------------------
[3294]339      !
[6140]340      READ  ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901)
[11536]341901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_qsr in reference namelist' )
[3294]342      !
[4147]343      READ  ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 )
[11536]344902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist' )
[4624]345      IF(lwm) WRITE ( numond, namtra_qsr )
[1423]346      !
347      IF(lwp) THEN                ! control print
348         WRITE(numout,*)
349         WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
350         WRITE(numout,*) '~~~~~~~~~~~~'
[1601]351         WRITE(numout,*) '   Namelist namtra_qsr : set the parameter of penetration'
[6140]352         WRITE(numout,*) '      RGB (Red-Green-Blue) light penetration       ln_qsr_rgb = ', ln_qsr_rgb
353         WRITE(numout,*) '      2 band               light penetration       ln_qsr_2bd = ', ln_qsr_2bd
354         WRITE(numout,*) '      bio-model            light penetration       ln_qsr_bio = ', ln_qsr_bio
355         WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)        nn_chldta  = ', nn_chldta
356         WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs     = ', rn_abs
357         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0     = ', rn_si0
358         WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1     = ', rn_si1
359         WRITE(numout,*)
[1423]360      ENDIF
[6140]361      !
362      ioptio = 0                    ! Parameter control
363      IF( ln_qsr_rgb  )   ioptio = ioptio + 1
364      IF( ln_qsr_2bd  )   ioptio = ioptio + 1
365      IF( ln_qsr_bio  )   ioptio = ioptio + 1
366      !
367      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE type of light penetration in namelist namtra_qsr',  &
368         &                               ' 2 bands, 3 RGB bands or bio-model light penetration' )
369      !
370      IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr = np_RGB 
371      IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr = np_RGBc
372      IF( ln_qsr_2bd                      )   nqsr = np_2BD
373      IF( ln_qsr_bio                      )   nqsr = np_BIO
374      !
375      !                             ! Initialisation
376      xsi0r = 1._wp / rn_si0
377      xsi1r = 1._wp / rn_si1
378      !
379      SELECT CASE( nqsr )
380      !                               
381      CASE( np_RGB , np_RGBc )         !==  Red-Green-Blue light penetration  ==!
382         !                             
[9169]383         IF(lwp)   WRITE(numout,*) '   ==>>>   R-G-B   light penetration '
[6140]384         !
385         CALL trc_oce_rgb( rkrgb )                 ! tabulated attenuation coef.
386         !                                   
387         nksr = trc_oce_ext_lev( r_si2, 33._wp )   ! level of light extinction
388         !
389         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
390         !
391         IF( nqsr == np_RGBc ) THEN                ! Chl data : set sf_chl structure
[9169]392            IF(lwp) WRITE(numout,*) '   ==>>>   Chlorophyll read in a file'
[6140]393            ALLOCATE( sf_chl(1), STAT=ierror )
394            IF( ierror > 0 ) THEN
395               CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' )   ;   RETURN
396            ENDIF
397            ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1)   )
398            IF( sn_chl%ln_tint )   ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) )
399            !                                        ! fill sf_chl with sn_chl and control print
400            CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',   &
[7646]401               &           'Solar penetration function of read chlorophyll', 'namtra_qsr' , no_print )
[1448]402         ENDIF
[6140]403         IF( nqsr == np_RGB ) THEN                 ! constant Chl
[9169]404            IF(lwp) WRITE(numout,*) '   ==>>>   Constant Chlorophyll concentration = 0.05'
[6140]405         ENDIF
[1448]406         !
[6140]407      CASE( np_2BD )                   !==  2 bands light penetration  ==!
[1448]408         !
[9169]409         IF(lwp)  WRITE(numout,*) '   ==>>>   2 bands light penetration'
[1448]410         !
[6140]411         nksr = trc_oce_ext_lev( rn_si1, 100._wp )    ! level of light extinction
412         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
[1455]413         !
[6140]414      CASE( np_BIO )                   !==  BIO light penetration  ==!
[1448]415         !
[9169]416         IF(lwp) WRITE(numout,*) '   ==>>>   bio-model light penetration'
[7646]417         IF( .NOT.lk_top )   CALL ctl_stop( 'No bio model : ln_qsr_bio = true impossible ' )
[1423]418         !
[13333]419         CALL trc_oce_rgb( rkrgb )                 ! tabulated attenuation coef.
420         !                                   
421         nksr = trc_oce_ext_lev( r_si2, 33._wp )   ! level of light extinction
422         !
423         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
424         !
[6140]425      END SELECT
[503]426      !
[7753]427      qsr_hc(:,:,:) = 0._wp     ! now qsr heat content set to zero where it will not be computed
[6140]428      !
429      ! 1st ocean level attenuation coefficient (used in sbcssm)
[5407]430      IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN
[13286]431         CALL iom_get( numror, jpdom_auto, 'fraqsr_1lev'  , fraqsr_1lev, ldxios = lrxios  )
[5407]432      ELSE
[7753]433         fraqsr_1lev(:,:) = 1._wp   ! default : no penetration
[5407]434      ENDIF
435      !
[9367]436      IF( lwxios ) THEN
437         CALL iom_set_rstw_var_active('qsr_hc_b')
438         CALL iom_set_rstw_var_active('fraqsr_1lev')
439      ENDIF
440      !
[3]441   END SUBROUTINE tra_qsr_init
442
443   !!======================================================================
444END MODULE traqsr
Note: See TracBrowser for help on using the repository browser.