source: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90 @ 11442

Last change on this file since 11442 was 11442, checked in by mattmartin, 21 months ago

Introduction of stochastic physics in NEMO, based on Andrea Storto's code.
For details, see ticket https://code.metoffice.gov.uk/trac/utils/ticket/251.

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