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 @ 12385

Last change on this file since 12385 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

  • Property svn:keywords set to Id
File size: 21.2 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
[12377]69#  include "do_loop_substitute.h90"
[3]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
[12377]77   SUBROUTINE tra_qsr( kt, Kmm, pts, 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      !!----------------------------------------------------------------------
[12377]103      INTEGER,                                   INTENT(in   ) :: kt            ! ocean time-step
104      INTEGER,                                   INTENT(in   ) :: Kmm, Krhs     ! time level indices
105      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts           ! active tracers and RHS of tracer equation
[2715]106      !
[6140]107      INTEGER  ::   ji, jj, jk               ! dummy loop indices
108      INTEGER  ::   irgb                     ! local integers
109      REAL(wp) ::   zchl, zcoef, z1_2        ! local scalars
110      REAL(wp) ::   zc0 , zc1 , zc2 , zc3    !    -         -
[4161]111      REAL(wp) ::   zzc0, zzc1, zzc2, zzc3   !    -         -
[6140]112      REAL(wp) ::   zz0 , zz1                !    -         -
[6403]113      REAL(wp) ::   zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze
114      REAL(wp) ::   zlogc, zlogc2, zlogc3 
[9019]115      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   :: zekb, zekg, zekr
116      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt
117      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zetot, zchl3d
[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  ==!
137!!gm case neuler  not taken into account....
138         IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN    ! read in restart
139            IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field read in the restart file'
140            z1_2 = 0.5_wp
[9367]141            CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b, ldxios = lrxios )   ! before heat content trend due to Qsr flux
[2528]142         ELSE                                           ! No restart or restart not found: Euler forward time stepping
[6140]143            z1_2 = 1._wp
[7753]144            qsr_hc_b(:,:,:) = 0._wp
[2528]145         ENDIF
[6140]146      ELSE                             !==  Swap of qsr heat content  ==!
147         z1_2 = 0.5_wp
[7753]148         qsr_hc_b(:,:,:) = qsr_hc(:,:,:)
[2528]149      ENDIF
[6140]150      !
151      !                         !--------------------------------!
152      SELECT CASE( nqsr )       !  now qsr induced heat content  !
153      !                         !--------------------------------!
154      !
155      CASE( np_BIO )                   !==  bio-model fluxes  ==!
156         !
157         DO jk = 1, nksr
[7753]158            qsr_hc(:,:,jk) = r1_rau0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) )
[2528]159         END DO
[6140]160         !
161      CASE( np_RGB , np_RGBc )         !==  R-G-B fluxes  ==!
162         !
[9019]163         ALLOCATE( zekb(jpi,jpj)     , zekg(jpi,jpj)     , zekr  (jpi,jpj)     , &
164            &      ze0 (jpi,jpj,jpk) , ze1 (jpi,jpj,jpk) , ze2   (jpi,jpj,jpk) , &
165            &      ze3 (jpi,jpj,jpk) , zea (jpi,jpj,jpk) , zchl3d(jpi,jpj,jpk)   ) 
[6140]166         !
167         IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll
168            CALL fld_read( kt, 1, sf_chl )         ! Read Chl data and provides it at the current time step
[6403]169            DO jk = 1, nksr + 1
170               DO jj = 2, jpjm1                       ! Separation in R-G-B depending of the surface Chl
[12377]171                  DO ji = 2, jpim1
[11410]172                     zchl    = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) )
[6403]173                     zCtot   = 40.6  * zchl**0.459
174                     zze     = 568.2 * zCtot**(-0.746)
175                     IF( zze > 102. ) zze = 200.0 * zCtot**(-0.293)
[12377]176                     zpsi    = gdepw(ji,jj,jk,Kmm) / zze
[6403]177                     !
178                     zlogc   = LOG( zchl )
179                     zlogc2  = zlogc * zlogc
180                     zlogc3  = zlogc * zlogc * zlogc
181                     zCb     = 0.768 + 0.087 * zlogc - 0.179 * zlogc2 - 0.025 * zlogc3
182                     zCmax   = 0.299 - 0.289 * zlogc + 0.579 * zlogc2
183                     zpsimax = 0.6   - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3
184                     zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2
185                     zCze    = 1.12  * (zchl)**0.803 
186                     !
187                     zchl3d(ji,jj,jk) = zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) / zdelpsi )**2 ) )
188                  END DO
189                  !
[3]190               END DO
191            END DO
[6140]192         ELSE                                !* constant chrlorophyll
[6403]193           DO jk = 1, nksr + 1
[7753]194              zchl3d(:,:,jk) = 0.05 
[6403]195            ENDDO
[4161]196         ENDIF
[1423]197         !
[6140]198         zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B
[12377]199         DO_2D_00_00
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_2D
[6140]206         !
[6403]207         DO jk = 2, nksr+1                   !* interior equi-partition in R-G-B depending of vertical profile of Chl
[12377]208            DO_2D_00_00
209               zchl = MIN( 10. , MAX( 0.03, zchl3d(ji,jj,jk) ) )
210               irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
211               zekb(ji,jj) = rkrgb(1,irgb)
212               zekg(ji,jj) = rkrgb(2,irgb)
213               zekr(ji,jj) = rkrgb(3,irgb)
214            END_2D
[6403]215
[12377]216            DO_2D_00_00
217               zc0 = ze0(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * xsi0r       )
218               zc1 = ze1(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekb(ji,jj) )
219               zc2 = ze2(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekg(ji,jj) )
220               zc3 = ze3(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekr(ji,jj) )
221               ze0(ji,jj,jk) = zc0
222               ze1(ji,jj,jk) = zc1
223               ze2(ji,jj,jk) = zc2
224               ze3(ji,jj,jk) = zc3
225               zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk)
226            END_2D
[6140]227         END DO
228         !
[12377]229         DO_3D_00_00( 1, nksr )
230            qsr_hc(ji,jj,jk) = r1_rau0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) )
231         END_3D
[187]232         !
[9019]233         DEALLOCATE( zekb , zekg , zekr , ze0 , ze1 , ze2 , ze3 , zea , zchl3d ) 
[6140]234         !
235      CASE( np_2BD  )            !==  2-bands fluxes  ==!
236         !
237         zz0 =        rn_abs   * r1_rau0_rcp      ! surface equi-partition in 2-bands
238         zz1 = ( 1. - rn_abs ) * r1_rau0_rcp
[12377]239         DO_3D_00_00( 1, nksr )
240            zc0 = zz0 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi1r )
241            zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r )
242            qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 
243         END_3D
[2528]244         !
[6140]245      END SELECT
246      !
247      !                          !-----------------------------!
[12377]248      DO_3D_00_00( 1, nksr )
249         pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
250            &                      + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t(ji,jj,jk,Kmm)
251      END_3D
[6140]252      !
[9019]253      ! sea-ice: store the 1st ocean level attenuation coefficient
[12377]254      DO_2D_00_00
255         IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) )
256         ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp
257         ENDIF
258      END_2D
[10425]259      CALL lbc_lnk( 'traqsr', fraqsr_1lev(:,:), 'T', 1._wp )
[2528]260      !
[6140]261      IF( iom_use('qsr3d') ) THEN      ! output the shortwave Radiation distribution
[9019]262         ALLOCATE( zetot(jpi,jpj,jpk) )
[7753]263         zetot(:,:,nksr+1:jpk) = 0._wp     ! below ~400m set to zero
[6140]264         DO jk = nksr, 1, -1
[9019]265            zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rau0_rcp
[6140]266         END DO         
267         CALL iom_put( 'qsr3d', zetot )   ! 3D distribution of shortwave Radiation
[9019]268         DEALLOCATE( zetot ) 
[2528]269      ENDIF
[6140]270      !
271      IF( lrst_oce ) THEN     ! write in the ocean restart file
[9367]272         IF( lwxios ) CALL iom_swap(      cwxios_context          )
273         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc     , ldxios = lwxios )
274         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev, ldxios = lwxios ) 
275         IF( lwxios ) CALL iom_swap(      cxios_context          )
[6140]276      ENDIF
277      !
[503]278      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
[12377]279         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
280         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt )
[9019]281         DEALLOCATE( ztrdt ) 
[3]282      ENDIF
[457]283      !                       ! print mean trends (used for debugging)
[12377]284      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
[503]285      !
[9019]286      IF( ln_timing )   CALL timing_stop('tra_qsr')
[3294]287      !
[3]288   END SUBROUTINE tra_qsr
289
290
291   SUBROUTINE tra_qsr_init
292      !!----------------------------------------------------------------------
293      !!                  ***  ROUTINE tra_qsr_init  ***
294      !!
295      !! ** Purpose :   Initialization for the penetrative solar radiation
296      !!
297      !! ** Method  :   The profile of solar radiation within the ocean is set
[1423]298      !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
[1601]299      !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
[3]300      !!      default values correspond to clear water (type I in Jerlov'
301      !!      (1968) classification.
302      !!         called by tra_qsr at the first timestep (nit000)
303      !!
[1423]304      !! ** Action  : - initialize rn_si0, rn_si1 and rn_abs
[3]305      !!
[503]306      !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
[3]307      !!----------------------------------------------------------------------
[6140]308      INTEGER  ::   ji, jj, jk                  ! dummy loop indices
309      INTEGER  ::   ios, irgb, ierror, ioptio   ! local integer
310      REAL(wp) ::   zz0, zc0 , zc1, zcoef      ! local scalars
311      REAL(wp) ::   zz1, zc2 , zc3, zchl       !   -      -
[2715]312      !
[1423]313      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files
314      TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read
[2715]315      !!
[9019]316      NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio,  &
[2528]317         &                  nn_chldta, rn_abs, rn_si0, rn_si1
[3]318      !!----------------------------------------------------------------------
[3294]319      !
[6140]320      READ  ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901)
[11536]321901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_qsr in reference namelist' )
[3294]322      !
[4147]323      READ  ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 )
[11536]324902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist' )
[4624]325      IF(lwm) WRITE ( numond, namtra_qsr )
[1423]326      !
327      IF(lwp) THEN                ! control print
328         WRITE(numout,*)
329         WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
330         WRITE(numout,*) '~~~~~~~~~~~~'
[1601]331         WRITE(numout,*) '   Namelist namtra_qsr : set the parameter of penetration'
[6140]332         WRITE(numout,*) '      RGB (Red-Green-Blue) light penetration       ln_qsr_rgb = ', ln_qsr_rgb
333         WRITE(numout,*) '      2 band               light penetration       ln_qsr_2bd = ', ln_qsr_2bd
334         WRITE(numout,*) '      bio-model            light penetration       ln_qsr_bio = ', ln_qsr_bio
335         WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)        nn_chldta  = ', nn_chldta
336         WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs     = ', rn_abs
337         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0     = ', rn_si0
338         WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1     = ', rn_si1
339         WRITE(numout,*)
[1423]340      ENDIF
[6140]341      !
342      ioptio = 0                    ! Parameter control
343      IF( ln_qsr_rgb  )   ioptio = ioptio + 1
344      IF( ln_qsr_2bd  )   ioptio = ioptio + 1
345      IF( ln_qsr_bio  )   ioptio = ioptio + 1
346      !
347      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE type of light penetration in namelist namtra_qsr',  &
348         &                               ' 2 bands, 3 RGB bands or bio-model light penetration' )
349      !
350      IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr = np_RGB 
351      IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr = np_RGBc
352      IF( ln_qsr_2bd                      )   nqsr = np_2BD
353      IF( ln_qsr_bio                      )   nqsr = np_BIO
354      !
355      !                             ! Initialisation
356      xsi0r = 1._wp / rn_si0
357      xsi1r = 1._wp / rn_si1
358      !
359      SELECT CASE( nqsr )
360      !                               
361      CASE( np_RGB , np_RGBc )         !==  Red-Green-Blue light penetration  ==!
362         !                             
[9169]363         IF(lwp)   WRITE(numout,*) '   ==>>>   R-G-B   light penetration '
[6140]364         !
365         CALL trc_oce_rgb( rkrgb )                 ! tabulated attenuation coef.
366         !                                   
367         nksr = trc_oce_ext_lev( r_si2, 33._wp )   ! level of light extinction
368         !
369         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
370         !
371         IF( nqsr == np_RGBc ) THEN                ! Chl data : set sf_chl structure
[9169]372            IF(lwp) WRITE(numout,*) '   ==>>>   Chlorophyll read in a file'
[6140]373            ALLOCATE( sf_chl(1), STAT=ierror )
374            IF( ierror > 0 ) THEN
375               CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' )   ;   RETURN
376            ENDIF
377            ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1)   )
378            IF( sn_chl%ln_tint )   ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) )
379            !                                        ! fill sf_chl with sn_chl and control print
380            CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',   &
[7646]381               &           'Solar penetration function of read chlorophyll', 'namtra_qsr' , no_print )
[1448]382         ENDIF
[6140]383         IF( nqsr == np_RGB ) THEN                 ! constant Chl
[9169]384            IF(lwp) WRITE(numout,*) '   ==>>>   Constant Chlorophyll concentration = 0.05'
[6140]385         ENDIF
[1448]386         !
[6140]387      CASE( np_2BD )                   !==  2 bands light penetration  ==!
[1448]388         !
[9169]389         IF(lwp)  WRITE(numout,*) '   ==>>>   2 bands light penetration'
[1448]390         !
[6140]391         nksr = trc_oce_ext_lev( rn_si1, 100._wp )    ! level of light extinction
392         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
[1455]393         !
[6140]394      CASE( np_BIO )                   !==  BIO light penetration  ==!
[1448]395         !
[9169]396         IF(lwp) WRITE(numout,*) '   ==>>>   bio-model light penetration'
[7646]397         IF( .NOT.lk_top )   CALL ctl_stop( 'No bio model : ln_qsr_bio = true impossible ' )
[1423]398         !
[6140]399      END SELECT
[503]400      !
[7753]401      qsr_hc(:,:,:) = 0._wp     ! now qsr heat content set to zero where it will not be computed
[6140]402      !
403      ! 1st ocean level attenuation coefficient (used in sbcssm)
[5407]404      IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN
[9367]405         CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev'  , fraqsr_1lev, ldxios = lrxios  )
[5407]406      ELSE
[7753]407         fraqsr_1lev(:,:) = 1._wp   ! default : no penetration
[5407]408      ENDIF
409      !
[9367]410      IF( lwxios ) THEN
411         CALL iom_set_rstw_var_active('qsr_hc_b')
412         CALL iom_set_rstw_var_active('fraqsr_1lev')
413      ENDIF
414      !
[3]415   END SUBROUTINE tra_qsr_init
416
417   !!======================================================================
418END MODULE traqsr
Note: See TracBrowser for help on using the repository browser.