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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
traqsr.F90 in branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago

Branch 2015/dev_r5803_NOC_WAD. Merge in trunk changes from 5803 to 5869 in preparation for merge. Also tidied and reorganised some wetting and drying code. Renamed wadlmt.F90 to wetdry.F90. Wetting drying code changes restricted to domzgr.F90, domvvl.F90 nemogcm.F90 sshwzv.F90, dynspg_ts.F90, wetdry.F90 and dynhpg.F90. Code passes full SETTE tests with ln_wd=.false.. Still awaiting test case for checking with ln_wd=.false.
  • Property svn:keywords set to Id
File size: 29.9 KB
Line 
1MODULE traqsr
2   !!======================================================================
3   !!                       ***  MODULE  traqsr  ***
4   !! Ocean physics: solar radiation penetration in the top ocean levels
5   !!======================================================================
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)
12   !!            4.0  !  2012-05  (C. Rousset) store attenuation coef for use in ice model
13   !!----------------------------------------------------------------------
14
15   !!----------------------------------------------------------------------
16   !!   tra_qsr      : trend due to the solar radiation penetration
17   !!   tra_qsr_init : solar radiation penetration initialization
18   !!----------------------------------------------------------------------
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 trd_oce        ! trends: ocean variables
24   USE trdtra         ! trends manager: tracers
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
32   USE wrk_nemo       ! Memory Allocation
33   USE timing         ! Timing
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   tra_qsr       ! routine called by step.F90 (ln_traqsr=T)
39   PUBLIC   tra_qsr_init  ! routine called by nemogcm.F90
40
41   !                                 !!* Namelist namtra_qsr: penetrative solar radiation
42   LOGICAL , PUBLIC ::   ln_traqsr    !: light absorption (qsr) flag
43   LOGICAL , PUBLIC ::   ln_qsr_rgb   !: Red-Green-Blue light absorption flag 
44   LOGICAL , PUBLIC ::   ln_qsr_2bd   !: 2 band         light absorption flag
45   LOGICAL , PUBLIC ::   ln_qsr_bio   !: bio-model      light absorption flag
46   LOGICAL , PUBLIC ::   ln_qsr_ice   !: light penetration for ice-model LIM3 (clem)
47   INTEGER , PUBLIC ::   nn_chldta    !: use Chlorophyll data (=1) or not (=0)
48   REAL(wp), PUBLIC ::   rn_abs       !: fraction absorbed in the very near surface (RGB & 2 bands)
49   REAL(wp), PUBLIC ::   rn_si0       !: very near surface depth of extinction      (RGB & 2 bands)
50   REAL(wp), PUBLIC ::   rn_si1       !: deepest depth of extinction (water type I)       (2 bands)
51 
52   ! Module variables
53   REAL(wp) ::   xsi0r                           !: inverse of rn_si0
54   REAL(wp) ::   xsi1r                           !: inverse of rn_si1
55   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read)
56   INTEGER, PUBLIC ::   nksr              ! levels below which the light cannot penetrate ( depth larger than 391 m)
57   REAL(wp), DIMENSION(3,61) ::   rkrgb   !: tabulated attenuation coefficients for RGB absorption
58
59   !! * Substitutions
60#  include "domzgr_substitute.h90"
61#  include "vectopt_loop_substitute.h90"
62   !!----------------------------------------------------------------------
63   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
64   !! $Id$
65   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
66   !!----------------------------------------------------------------------
67CONTAINS
68
69   SUBROUTINE tra_qsr( kt )
70      !!----------------------------------------------------------------------
71      !!                  ***  ROUTINE tra_qsr  ***
72      !!
73      !! ** Purpose :   Compute the temperature trend due to the solar radiation
74      !!      penetration and add it to the general temperature trend.
75      !!
76      !! ** Method  : The profile of the solar radiation within the ocean is defined
77      !!      through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) and a ratio rn_abs
78      !!      Considering the 2 wavebands case:
79      !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
80      !!         The temperature trend associated with the solar radiation penetration
81      !!         is given by : zta = 1/e3t dk[ I ] / (rau0*Cp)
82      !!         At the bottom, boudary condition for the radiation is no flux :
83      !!      all heat which has not been absorbed in the above levels is put
84      !!      in the last ocean level.
85      !!         In z-coordinate case, the computation is only done down to the
86      !!      level where I(k) < 1.e-15 W/m2. In addition, the coefficients
87      !!      used for the computation are calculated one for once as they
88      !!      depends on k only.
89      !!
90      !! ** Action  : - update ta with the penetrative solar radiation trend
91      !!              - save the trend in ttrd ('key_trdtra')
92      !!
93      !! Reference  : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
94      !!              Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
95      !!----------------------------------------------------------------------
96      !
97      INTEGER, INTENT(in) ::   kt     ! ocean time-step
98      !
99      INTEGER  ::   ji, jj, jk           ! dummy loop indices
100      INTEGER  ::   irgb                 ! local integers
101      REAL(wp) ::   zchl, zcoef, zfact   ! local scalars
102      REAL(wp) ::   zc0, zc1, zc2, zc3   !    -         -
103      REAL(wp) ::   zzc0, zzc1, zzc2, zzc3   !    -         -
104      REAL(wp) ::   zz0, zz1, z1_e3t     !    -         -
105      REAL(wp), POINTER, DIMENSION(:,:  ) :: zekb, zekg, zekr
106      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt
107      !!----------------------------------------------------------------------
108      !
109      IF( nn_timing == 1 )  CALL timing_start('tra_qsr')
110      !
111      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        ) 
112      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 
113      !
114      IF( kt == nit000 ) THEN
115         IF(lwp) WRITE(numout,*)
116         IF(lwp) WRITE(numout,*) 'tra_qsr : penetration of the surface solar radiation'
117         IF(lwp) WRITE(numout,*) '~~~~~~~'
118         IF( .NOT.ln_traqsr )   RETURN
119      ENDIF
120
121      IF( l_trdtra ) THEN      ! Save ta and sa trends
122         CALL wrk_alloc( jpi, jpj, jpk, ztrdt ) 
123         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
124      ENDIF
125
126      !                                        Set before qsr tracer content field
127      !                                        ***********************************
128      IF( kt == nit000 ) THEN                     ! Set the forcing field at nit000 - 1
129         !                                        ! -----------------------------------
130         qsr_hc(:,:,:) = 0.e0
131         !
132         IF( ln_rstart .AND.    &                    ! Restart: read in restart file
133              & iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN
134            IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field red in the restart file'
135            zfact = 0.5e0
136            CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b )   ! before heat content trend due to Qsr flux
137         ELSE                                           ! No restart or restart not found: Euler forward time stepping
138            zfact = 1.e0
139            qsr_hc_b(:,:,:) = 0.e0
140         ENDIF
141      ELSE                                        ! Swap of forcing field
142         !                                        ! ---------------------
143         zfact = 0.5e0
144         qsr_hc_b(:,:,:) = qsr_hc(:,:,:)
145      ENDIF
146      !                                        Compute now qsr tracer content field
147      !                                        ************************************
148     
149      !                                           ! ============================================== !
150      IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN      !  bio-model fluxes  : all vertical coordinates  !
151         !                                        ! ============================================== !
152         DO jk = 1, jpkm1
153            qsr_hc(:,:,jk) = r1_rau0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) )
154         END DO
155         !                                        Add to the general trend
156         DO jk = 1, jpkm1
157            DO jj = 2, jpjm1 
158               DO ji = fs_2, fs_jpim1   ! vector opt.
159                  z1_e3t = zfact / fse3t(ji,jj,jk)
160                  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
161               END DO
162            END DO
163         END DO
164         CALL iom_put( 'qsr3d', etot3 )   ! Shortwave Radiation 3D distribution
165         ! clem: store attenuation coefficient of the first ocean level
166         IF ( ln_qsr_ice ) THEN
167            DO jj = 1, jpj
168               DO ji = 1, jpi
169                  IF ( qsr(ji,jj) /= 0._wp ) THEN
170                     fraqsr_1lev(ji,jj) = ( qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) ) )
171                  ELSE
172                     fraqsr_1lev(ji,jj) = 1.
173                  ENDIF
174               END DO
175            END DO
176         ENDIF
177         !                                        ! ============================================== !
178      ELSE                                        !  Ocean alone :
179         !                                        ! ============================================== !
180         !
181         !                                                ! ------------------------- !
182         IF( ln_qsr_rgb) THEN                             !  R-G-B  light penetration !
183            !                                             ! ------------------------- !
184            ! Set chlorophyl concentration
185            IF( nn_chldta == 1 .OR. lk_vvl ) THEN            !*  Variable Chlorophyll or ocean volume
186               !
187               IF( nn_chldta == 1 ) THEN                             !*  Variable Chlorophyll
188                  !
189                  CALL fld_read( kt, 1, sf_chl )                         ! Read Chl data and provides it at the current time step
190                  !         
191                  DO jj = 1, jpj                                         ! Separation in R-G-B depending of the surface Chl
192                     DO ji = 1, jpi
193                        zchl = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) )
194                        irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
195                        zekb(ji,jj) = rkrgb(1,irgb)
196                        zekg(ji,jj) = rkrgb(2,irgb)
197                        zekr(ji,jj) = rkrgb(3,irgb)
198                     END DO
199                  END DO
200               ELSE                                            ! Variable ocean volume but constant chrlorophyll
201                  zchl = 0.05                                     ! constant chlorophyll
202                  irgb = NINT( 41 + 20.*LOG10( zchl ) + 1.e-15 )
203                  zekb(:,:) = rkrgb(1,irgb)                       ! Separation in R-G-B depending of the chlorophyll
204                  zekg(:,:) = rkrgb(2,irgb)
205                  zekr(:,:) = rkrgb(3,irgb)
206               ENDIF
207               !
208               zcoef  = ( 1. - rn_abs ) / 3.e0                        ! equi-partition in R-G-B
209               ze0(:,:,1) = rn_abs  * qsr(:,:)
210               ze1(:,:,1) = zcoef * qsr(:,:)
211               ze2(:,:,1) = zcoef * qsr(:,:)
212               ze3(:,:,1) = zcoef * qsr(:,:)
213               zea(:,:,1) =         qsr(:,:)
214               !
215               DO jk = 2, nksr+1
216                  DO jj = 1, jpj
217                     DO ji = 1, jpi
218                        zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r     )
219                        zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekb(ji,jj) )
220                        zc2 = ze2(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekg(ji,jj) )
221                        zc3 = ze3(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekr(ji,jj) )
222                        ze0(ji,jj,jk) = zc0
223                        ze1(ji,jj,jk) = zc1
224                        ze2(ji,jj,jk) = zc2
225                        ze3(ji,jj,jk) = zc3
226                        zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk)
227                     END DO
228                  END DO
229               END DO
230               ! clem: store attenuation coefficient of the first ocean level
231               IF ( ln_qsr_ice ) THEN
232                  DO jj = 1, jpj
233                     DO ji = 1, jpi
234                        zzc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r     )
235                        zzc1 = zcoef  * EXP( - fse3t(ji,jj,1) * zekb(ji,jj) )
236                        zzc2 = zcoef  * EXP( - fse3t(ji,jj,1) * zekg(ji,jj) )
237                        zzc3 = zcoef  * EXP( - fse3t(ji,jj,1) * zekr(ji,jj) )
238                        fraqsr_1lev(ji,jj) = 1.0 - ( zzc0 + zzc1 + zzc2  + zzc3  ) * tmask(ji,jj,2) 
239                     END DO
240                  END DO
241               ENDIF
242               !
243               DO jk = 1, nksr                                        ! compute and add qsr trend to ta
244                  qsr_hc(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) )
245               END DO
246               zea(:,:,nksr+1:jpk) = 0.e0     ! below 400m set to zero
247               CALL iom_put( 'qsr3d', zea )   ! Shortwave Radiation 3D distribution
248               !
249            ELSE                                                 !*  Constant Chlorophyll
250               DO jk = 1, nksr
251                  qsr_hc(:,:,jk) =  etot3(:,:,jk) * qsr(:,:)
252               END DO
253               ! clem: store attenuation coefficient of the first ocean level
254               IF ( ln_qsr_ice ) THEN
255                  fraqsr_1lev(:,:) = etot3(:,:,1) / r1_rau0_rcp
256               ENDIF
257           ENDIF
258
259         ENDIF
260         !                                                ! ------------------------- !
261         IF( ln_qsr_2bd ) THEN                            !  2 band light penetration !
262            !                                             ! ------------------------- !
263            !
264            IF( lk_vvl ) THEN                                  !* variable volume
265               zz0   =        rn_abs   * r1_rau0_rcp
266               zz1   = ( 1. - rn_abs ) * r1_rau0_rcp
267               DO jk = 1, nksr                    ! solar heat absorbed at T-point in the top 400m
268                  DO jj = 1, jpj
269                     DO ji = 1, jpi
270                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
271                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
272                        qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) ) 
273                     END DO
274                  END DO
275               END DO
276               ! clem: store attenuation coefficient of the first ocean level
277               IF ( ln_qsr_ice ) THEN
278                  DO jj = 1, jpj
279                     DO ji = 1, jpi
280                        zc0 = zz0 * EXP( -fsdepw(ji,jj,1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,1)*xsi1r )
281                        zc1 = zz0 * EXP( -fsdepw(ji,jj,2)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,2)*xsi1r )
282                        fraqsr_1lev(ji,jj) = ( zc0*tmask(ji,jj,1) - zc1*tmask(ji,jj,2) ) / r1_rau0_rcp
283                     END DO
284                  END DO
285               ENDIF
286            ELSE                                               !* constant volume: coef. computed one for all
287               DO jk = 1, nksr
288                  DO jj = 2, jpjm1
289                     DO ji = fs_2, fs_jpim1   ! vector opt.
290                        ! (ISF) no light penetration below the ice shelves         
291                        qsr_hc(ji,jj,jk) =  etot3(ji,jj,jk) * qsr(ji,jj) * tmask(ji,jj,1)
292                     END DO
293                  END DO
294               END DO
295               ! clem: store attenuation coefficient of the first ocean level
296               IF ( ln_qsr_ice ) THEN
297                  fraqsr_1lev(:,:) = etot3(:,:,1) / r1_rau0_rcp
298               ENDIF
299               !
300            ENDIF
301            !
302         ENDIF
303         !
304         !                                        Add to the general trend
305         DO jk = 1, nksr
306            DO jj = 2, jpjm1 
307               DO ji = fs_2, fs_jpim1   ! vector opt.
308                  z1_e3t = zfact / fse3t(ji,jj,jk)
309                  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
310               END DO
311            END DO
312         END DO
313         !
314      ENDIF
315      !
316      IF( lrst_oce ) THEN   !                  Write in the ocean restart file
317         !                                     *******************************
318         IF(lwp) WRITE(numout,*)
319         IF(lwp) WRITE(numout,*) 'qsr tracer content forcing field written in ocean restart file ',   &
320            &                    'at it= ', kt,' date= ', ndastp
321         IF(lwp) WRITE(numout,*) '~~~~'
322         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc      )
323         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev )   ! default definition in sbcssm
324         !
325      ENDIF
326
327      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
328         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
329         CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt )
330         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt ) 
331      ENDIF
332      !                       ! print mean trends (used for debugging)
333      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
334      !
335      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        ) 
336      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 
337      !
338      IF( nn_timing == 1 )  CALL timing_stop('tra_qsr')
339      !
340   END SUBROUTINE tra_qsr
341
342
343   SUBROUTINE tra_qsr_init
344      !!----------------------------------------------------------------------
345      !!                  ***  ROUTINE tra_qsr_init  ***
346      !!
347      !! ** Purpose :   Initialization for the penetrative solar radiation
348      !!
349      !! ** Method  :   The profile of solar radiation within the ocean is set
350      !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
351      !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
352      !!      default values correspond to clear water (type I in Jerlov'
353      !!      (1968) classification.
354      !!         called by tra_qsr at the first timestep (nit000)
355      !!
356      !! ** Action  : - initialize rn_si0, rn_si1 and rn_abs
357      !!
358      !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
359      !!----------------------------------------------------------------------
360      !
361      INTEGER  ::   ji, jj, jk                   ! dummy loop indices
362      INTEGER  ::   irgb, ierror, ioptio, nqsr   ! local integer
363      INTEGER  ::   ios                          ! Local integer output status for namelist read
364      REAL(wp) ::   zz0, zc0  , zc1, zcoef       ! local scalars
365      REAL(wp) ::   zz1, zc2  , zc3, zchl        !   -      -
366      REAL(wp), POINTER, DIMENSION(:,:  ) :: zekb, zekg, zekr
367      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea
368      !
369      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files
370      TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read
371      !!
372      NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice,  &
373         &                  nn_chldta, rn_abs, rn_si0, rn_si1
374      !!----------------------------------------------------------------------
375
376      !
377      IF( nn_timing == 1 )  CALL timing_start('tra_qsr_init')
378      !
379      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        ) 
380      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 
381      !
382
383      REWIND( numnam_ref )              ! Namelist namtra_qsr in reference namelist : Ratio and length of penetration
384      READ  ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901)
385901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in reference namelist', lwp )
386
387      REWIND( numnam_cfg )              !  Namelist namtra_qsr in configuration namelist : Ratio and length of penetration
388      READ  ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 )
389902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist', lwp )
390      IF(lwm) WRITE ( numond, namtra_qsr )
391      !
392      IF(lwp) THEN                ! control print
393         WRITE(numout,*)
394         WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
395         WRITE(numout,*) '~~~~~~~~~~~~'
396         WRITE(numout,*) '   Namelist namtra_qsr : set the parameter of penetration'
397         WRITE(numout,*) '      Light penetration (T) or not (F)         ln_traqsr  = ', ln_traqsr
398         WRITE(numout,*) '      RGB (Red-Green-Blue) light penetration   ln_qsr_rgb = ', ln_qsr_rgb
399         WRITE(numout,*) '      2 band               light penetration   ln_qsr_2bd = ', ln_qsr_2bd
400         WRITE(numout,*) '      bio-model            light penetration   ln_qsr_bio = ', ln_qsr_bio
401         WRITE(numout,*) '      light penetration for ice-model LIM3     ln_qsr_ice = ', ln_qsr_ice
402         WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)    nn_chldta  = ', nn_chldta
403         WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs = ', rn_abs
404         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0 = ', rn_si0
405         WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1 = ', rn_si1
406      ENDIF
407
408      IF( ln_traqsr ) THEN     ! control consistency
409         !                     
410         IF( .NOT.lk_qsr_bio .AND. ln_qsr_bio )   THEN
411            CALL ctl_warn( 'No bio model : force ln_qsr_bio = FALSE ' )
412            ln_qsr_bio = .FALSE.
413         ENDIF
414         !
415         ioptio = 0                      ! Parameter control
416         IF( ln_qsr_rgb  )   ioptio = ioptio + 1
417         IF( ln_qsr_2bd  )   ioptio = ioptio + 1
418         IF( ln_qsr_bio  )   ioptio = ioptio + 1
419         !
420         IF( ioptio /= 1 ) &
421            CALL ctl_stop( '          Choose ONE type of light penetration in namelist namtra_qsr',  &
422            &              ' 2 bands, 3 RGB bands or bio-model light penetration' )
423         !
424         IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr =  1 
425         IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr =  2
426         IF( ln_qsr_2bd                      )   nqsr =  3
427         IF( ln_qsr_bio                      )   nqsr =  4
428         !
429         IF(lwp) THEN                   ! Print the choice
430            WRITE(numout,*)
431            IF( nqsr ==  1 )   WRITE(numout,*) '         R-G-B   light penetration - Constant Chlorophyll'
432            IF( nqsr ==  2 )   WRITE(numout,*) '         R-G-B   light penetration - Chl data '
433            IF( nqsr ==  3 )   WRITE(numout,*) '         2 bands light penetration'
434            IF( nqsr ==  4 )   WRITE(numout,*) '         bio-model light penetration'
435         ENDIF
436         !
437      ENDIF
438      !                          ! ===================================== !
439      IF( ln_traqsr  ) THEN      !  Initialisation of Light Penetration  ! 
440         !                       ! ===================================== !
441         !
442         xsi0r = 1.e0 / rn_si0
443         xsi1r = 1.e0 / rn_si1
444         !                                ! ---------------------------------- !
445         IF( ln_qsr_rgb ) THEN            !  Red-Green-Blue light penetration  !
446            !                             ! ---------------------------------- !
447            !
448            CALL trc_oce_rgb( rkrgb )           !* tabulated attenuation coef.
449            !
450            !                                   !* level of light extinction
451            IF(  ln_sco ) THEN   ;   nksr = jpkm1
452            ELSE                 ;   nksr = trc_oce_ext_lev( r_si2, 0.33e2 )
453            ENDIF
454
455            IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
456            !
457            IF( nn_chldta == 1 ) THEN           !* Chl data : set sf_chl structure
458               IF(lwp) WRITE(numout,*)
459               IF(lwp) WRITE(numout,*) '        Chlorophyll read in a file'
460               ALLOCATE( sf_chl(1), STAT=ierror )
461               IF( ierror > 0 ) THEN
462                  CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' )   ;   RETURN
463               ENDIF
464               ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1)   )
465               IF( sn_chl%ln_tint )ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) )
466               !                                        ! fill sf_chl with sn_chl and control print
467               CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',   &
468                  &                                         'Solar penetration function of read chlorophyll', 'namtra_qsr' )
469               !
470            ELSE                                !* constant Chl : compute once for all the distribution of light (etot3)
471               IF(lwp) WRITE(numout,*)
472               IF(lwp) WRITE(numout,*) '        Constant Chlorophyll concentration = 0.05'
473               IF( lk_vvl ) THEN                   ! variable volume
474                  IF(lwp) WRITE(numout,*) '        key_vvl: light distribution will be computed at each time step'
475               ELSE                                ! constant volume: computes one for all
476                  IF(lwp) WRITE(numout,*) '        fixed volume: light distribution computed one for all'
477                  !
478                  zchl = 0.05                                 ! constant chlorophyll
479                  irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
480                  zekb(:,:) = rkrgb(1,irgb)                   ! Separation in R-G-B depending of the chlorophyll
481                  zekg(:,:) = rkrgb(2,irgb)
482                  zekr(:,:) = rkrgb(3,irgb)
483                  !
484                  zcoef = ( 1. - rn_abs ) / 3.e0              ! equi-partition in R-G-B
485                  ze0(:,:,1) = rn_abs
486                  ze1(:,:,1) = zcoef
487                  ze2(:,:,1) = zcoef 
488                  ze3(:,:,1) = zcoef
489                  zea(:,:,1) = tmask(:,:,1)                   ! = ( ze0+ze1+z2+ze3 ) * tmask
490               
491                  DO jk = 2, nksr+1
492                     DO jj = 1, jpj
493                        DO ji = 1, jpi
494                           zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r     )
495                           zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekb(ji,jj) )
496                           zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekg(ji,jj) )
497                           zc3 = ze3(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekr(ji,jj) )
498                           ze0(ji,jj,jk) = zc0
499                           ze1(ji,jj,jk) = zc1
500                           ze2(ji,jj,jk) = zc2
501                           ze3(ji,jj,jk) = zc3
502                           zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk)
503                        END DO
504                     END DO
505                  END DO 
506                  !
507                  DO jk = 1, nksr
508                     ! (ISF) no light penetration below the ice shelves
509                     etot3(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) ) * tmask(:,:,1)
510                  END DO
511                  etot3(:,:,nksr+1:jpk) = 0.e0                ! below 400m set to zero
512               ENDIF
513            ENDIF
514            !
515         ENDIF
516            !                             ! ---------------------------------- !
517         IF( ln_qsr_2bd ) THEN            !    2 bands    light penetration    !
518            !                             ! ---------------------------------- !
519            !
520            !                                ! level of light extinction
521            nksr = trc_oce_ext_lev( rn_si1, 1.e2 )
522            IF(lwp) THEN
523               WRITE(numout,*)
524            IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
525            ENDIF
526            !
527            IF( lk_vvl ) THEN                   ! variable volume
528               IF(lwp) WRITE(numout,*) '        key_vvl: light distribution will be computed at each time step'
529            ELSE                                ! constant volume: computes one for all
530               zz0 =        rn_abs   * r1_rau0_rcp
531               zz1 = ( 1. - rn_abs ) * r1_rau0_rcp
532               DO jk = 1, nksr                    !*  solar heat absorbed at T-point computed once for all
533                  DO jj = 1, jpj                              ! top 400 meters
534                     DO ji = 1, jpi
535                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
536                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
537                        etot3(ji,jj,jk) = (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) * tmask(ji,jj,1) 
538                     END DO
539                  END DO
540               END DO
541               etot3(:,:,nksr+1:jpk) = 0.e0                   ! below 400m set to zero
542               !
543            ENDIF
544         ENDIF
545         !                       ! ===================================== !
546      ELSE                       !        No light penetration           !                   
547         !                       ! ===================================== !
548         IF(lwp) THEN
549            WRITE(numout,*)
550            WRITE(numout,*) 'tra_qsr_init : NO solar flux penetration'
551            WRITE(numout,*) '~~~~~~~~~~~~'
552         ENDIF
553      ENDIF
554      !
555      ! initialisation of fraqsr_1lev used in sbcssm
556      IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN
557         CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev'  , fraqsr_1lev  )
558      ELSE
559         fraqsr_1lev(:,:) = 1._wp   ! default definition
560      ENDIF
561      !
562      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        ) 
563      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 
564      !
565      IF( nn_timing == 1 )  CALL timing_stop('tra_qsr_init')
566      !
567   END SUBROUTINE tra_qsr_init
568
569   !!======================================================================
570END MODULE traqsr
Note: See TracBrowser for help on using the repository browser.