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 trunk/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90 @ 4624

Last change on this file since 4624 was 4624, checked in by acc, 10 years ago

#1305. Fix slow start-up problems on some systems by introducing and using lwm logical to restrict output of merged namelists to the first (or only) processor. lwm is true only on the first processor regardless of ln_ctl. Small changes to all flavours of nemogcm.F90 are also required to write namctl and namcfg after the call to mynode which now opens output.namelist.dyn and writes nammpp.
  • 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 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
32   USE wrk_nemo       ! Memory Allocation
33   USE timing         ! Timing
34   USE sbc_ice, ONLY : lk_lim3
35
36   IMPLICIT NONE
37   PRIVATE
38
39   PUBLIC   tra_qsr       ! routine called by step.F90 (ln_traqsr=T)
40   PUBLIC   tra_qsr_init  ! routine called by opa.F90
41
42   !                                 !!* Namelist namtra_qsr: penetrative solar radiation
43   LOGICAL , PUBLIC ::   ln_traqsr    !: light absorption (qsr) flag
44   LOGICAL , PUBLIC ::   ln_qsr_rgb   !: Red-Green-Blue light absorption flag 
45   LOGICAL , PUBLIC ::   ln_qsr_2bd   !: 2 band         light absorption flag
46   LOGICAL , PUBLIC ::   ln_qsr_bio   !: bio-model      light absorption flag
47   LOGICAL , PUBLIC ::   ln_qsr_ice   !: light penetration for ice-model LIM3 (clem)
48   INTEGER , PUBLIC ::   nn_chldta    !: use Chlorophyll data (=1) or not (=0)
49   REAL(wp), PUBLIC ::   rn_abs       !: fraction absorbed in the very near surface (RGB & 2 bands)
50   REAL(wp), PUBLIC ::   rn_si0       !: very near surface depth of extinction      (RGB & 2 bands)
51   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
56   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
59
60   !! * Substitutions
61#  include "domzgr_substitute.h90"
62#  include "vectopt_loop_substitute.h90"
63   !!----------------------------------------------------------------------
64   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
65   !! $Id$
66   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
67   !!----------------------------------------------------------------------
68CONTAINS
69
70   SUBROUTINE tra_qsr( kt )
71      !!----------------------------------------------------------------------
72      !!                  ***  ROUTINE tra_qsr  ***
73      !!
74      !! ** Purpose :   Compute the temperature trend due to the solar radiation
75      !!      penetration and add it to the general temperature trend.
76      !!
77      !! ** Method  : The profile of the solar radiation within the ocean is defined
78      !!      through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) and a ratio rn_abs
79      !!      Considering the 2 wavebands case:
80      !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
81      !!         The temperature trend associated with the solar radiation penetration
82      !!         is given by : zta = 1/e3t dk[ I ] / (rau0*Cp)
83      !!         At the bottom, boudary condition for the radiation is no flux :
84      !!      all heat which has not been absorbed in the above levels is put
85      !!      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.
90      !!
91      !! ** Action  : - update ta with the penetrative solar radiation trend
92      !!              - save the trend in ttrd ('key_trdtra')
93      !!
94      !! Reference  : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
95      !!              Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
96      !!----------------------------------------------------------------------
97      !
98      INTEGER, INTENT(in) ::   kt     ! ocean time-step
99      !
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   !    -         -
104      REAL(wp) ::   zzc0, zzc1, zzc2, zzc3   !    -         -
105      REAL(wp) ::   zz0, zz1, z1_e3t     !    -         -
106      REAL(wp), POINTER, DIMENSION(:,:  ) :: zekb, zekg, zekr
107      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt
108      !!----------------------------------------------------------------------
109      !
110      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 ) 
114      !
115      IF( kt == nit000 ) THEN
116         IF(lwp) WRITE(numout,*)
117         IF(lwp) WRITE(numout,*) 'tra_qsr : penetration of the surface solar radiation'
118         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 ) 
124         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
125      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
135            CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b )   ! before heat content trend due to Qsr flux
136         ELSE                                           ! No restart or restart not found: Euler forward time stepping
137            zfact = 1.e0
138            qsr_hc_b(:,:,:) = 0.e0
139         ENDIF
140      ELSE                                        ! Swap of forcing field
141         !                                        ! ---------------------
142         zfact = 0.5e0
143         qsr_hc_b(:,:,:) = qsr_hc(:,:,:)
144      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
152            qsr_hc(:,:,jk) = r1_rau0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) )
153         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
160               END DO
161            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
172               END DO
173            END DO
174         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)
207               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
332      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
333         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 ) 
336      ENDIF
337      !                       ! print mean trends (used for debugging)
338      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 ) 
342      !
343      IF( nn_timing == 1 )  CALL timing_stop('tra_qsr')
344      !
345   END SUBROUTINE tra_qsr
346
347
348   SUBROUTINE tra_qsr_init
349      !!----------------------------------------------------------------------
350      !!                  ***  ROUTINE tra_qsr_init  ***
351      !!
352      !! ** Purpose :   Initialization for the penetrative solar radiation
353      !!
354      !! ** Method  :   The profile of solar radiation within the ocean is set
355      !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
356      !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
357      !!      default values correspond to clear water (type I in Jerlov'
358      !!      (1968) classification.
359      !!         called by tra_qsr at the first timestep (nit000)
360      !!
361      !! ** Action  : - initialize rn_si0, rn_si1 and rn_abs
362      !!
363      !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
364      !!----------------------------------------------------------------------
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
373      !
374      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files
375      TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read
376      !!
377      NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice,  &
378         &                  nn_chldta, rn_abs, rn_si0, rn_si1
379      !!----------------------------------------------------------------------
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
395      READ  ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901)
396901   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
399      READ  ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 )
400902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist', lwp )
401      IF(lwm) WRITE ( numond, namtra_qsr )
402      !
403      IF(lwp) THEN                ! control print
404         WRITE(numout,*)
405         WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
406         WRITE(numout,*) '~~~~~~~~~~~~'
407         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.
425         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'
447         ENDIF
448         !
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')
572      !
573   END SUBROUTINE tra_qsr_init
574
575   !!======================================================================
576END MODULE traqsr
Note: See TracBrowser for help on using the repository browser.