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

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

Last change on this file was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

File size: 33.0 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
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 nemogcm.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      !!              Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562
97      !!----------------------------------------------------------------------
98      !
99      INTEGER, INTENT(in) ::   kt     ! ocean time-step
100      !
101      INTEGER  ::   ji, jj, jk           ! dummy loop indices
102      INTEGER  ::   irgb                 ! local integers
103      REAL(wp) ::   zchl, zcoef, zfact   ! local scalars
104      REAL(wp) ::   zc0, zc1, zc2, zc3   !    -         -
105      REAL(wp) ::   zz0, zz1, z1_e3t     !    -         -
106      REAL(wp) ::   zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze
107      REAL(wp) ::   zlogc, zlogc2, zlogc3 
108      REAL(wp), POINTER, DIMENSION(:,:  ) :: zekb, zekg, zekr
109      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt, zchl3d
110      !!--------------------------------------------------------------------------
111      !
112      IF( nn_timing == 1 )  CALL timing_start('tra_qsr')
113      !
114      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        ) 
115      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 
116      !
117      IF( kt == nit000 ) THEN
118         IF(lwp) WRITE(numout,*)
119         IF(lwp) WRITE(numout,*) 'tra_qsr : penetration of the surface solar radiation'
120         IF(lwp) WRITE(numout,*) '~~~~~~~'
121         IF(lwp .AND. lflush) CALL flush(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 .AND. nprint >0) THEN
139               WRITE(numout,*) '          nit000-1 qsr tracer content forcing field red in the restart file'
140               IF(lflush) CALL flush(numout)
141            ENDIF
142            zfact = 0.5e0
143            CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b )   ! before heat content trend due to Qsr flux
144         ELSE                                           ! No restart or restart not found: Euler forward time stepping
145            zfact = 1.e0
146            qsr_hc_b(:,:,:) = 0.e0
147         ENDIF
148      ELSE                                        ! Swap of forcing field
149         !                                        ! ---------------------
150         zfact = 0.5e0
151         qsr_hc_b(:,:,:) = qsr_hc(:,:,:)
152      ENDIF
153      !                                        Compute now qsr tracer content field
154      !                                        ************************************
155     
156      !                                           ! ============================================== !
157      IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN      !  bio-model fluxes  : all vertical coordinates  !
158         !                                        ! ============================================== !
159         DO jk = 1, jpkm1
160            qsr_hc(:,:,jk) = r1_rau0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) )
161         END DO
162         !                                        Add to the general trend
163         DO jk = 1, jpkm1
164            DO jj = 2, jpjm1 
165               DO ji = fs_2, fs_jpim1   ! vector opt.
166                  z1_e3t = zfact / fse3t(ji,jj,jk)
167                  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
168               END DO
169            END DO
170         END DO
171         CALL iom_put( 'qsr3d', etot3 )   ! Shortwave Radiation 3D distribution
172         ! clem: store attenuation coefficient of the first ocean level
173         IF ( ln_qsr_ice ) THEN
174            DO jj = 1, jpj
175               DO ji = 1, jpi
176                  IF ( qsr(ji,jj) /= 0._wp ) THEN
177                     fraqsr_1lev(ji,jj) = ( qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) ) )
178                  ELSE
179                     fraqsr_1lev(ji,jj) = 1.
180                  ENDIF
181               END DO
182            END DO
183         ENDIF
184         !                                        ! ============================================== !
185      ELSE                                        !  Ocean alone :
186         !                                        ! ============================================== !
187         !
188         !                                                ! ------------------------- !
189         IF( ln_qsr_rgb) THEN                             !  R-G-B  light penetration !
190            !                                             ! ------------------------- !
191            ! Set chlorophyl concentration
192            IF( nn_chldta == 1 .OR. nn_chldta == 2 .OR. lk_vvl ) THEN    !*  Variable Chlorophyll or ocean volume
193               !
194               IF( nn_chldta == 1 ) THEN        !*  2D Variable Chlorophyll
195                  !
196                  CALL fld_read( kt, 1, sf_chl )            ! Read Chl data and provides it at the current time step
197                  DO jk = 1, nksr + 1
198                     zchl3d(:,:,jk) = sf_chl(1)%fnow(:,:,1) 
199                  ENDDO
200                  !
201               ELSE IF( nn_chldta == 2 ) THEN    !*   -3-D Variable Chlorophyll
202                  !
203                  CALL fld_read( kt, 1, sf_chl )            ! Read Chl data and provides it at the current time step
204!CDIR NOVERRCHK   !
205                  DO jj = 1, jpj
206!CDIR NOVERRCHK
207                     DO ji = 1, jpi
208                        zchl    = sf_chl(1)%fnow(ji,jj,1)
209                        zCtot   = 40.6  * zchl**0.459
210                        zze     = 568.2 * zCtot**(-0.746)
211                        IF( zze > 102. ) zze = 200.0 * zCtot**(-0.293)
212                        zlogc   = LOG( zchl )
213                        zlogc2  = zlogc * zlogc
214                        zlogc3  = zlogc * zlogc * zlogc
215                        zCb     = 0.768 + 0.087 * zlogc - 0.179 * zlogc2 - 0.025 * zlogc3
216                        zCmax   = 0.299 - 0.289 * zlogc + 0.579 * zlogc2
217                        zpsimax = 0.6   - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3
218                        zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2
219                        zCze    = 1.12  * (zchl)**0.803 
220                        DO jk = 1, nksr + 1
221                           zpsi = fsdept(ji,jj,jk) / zze
222                           zchl3d(ji,jj,jk) = zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) / zdelpsi )**2 ) )
223                        END DO
224                        !
225                      END DO
226                   END DO
227                     !
228               ELSE                              !* Variable ocean volume but constant chrlorophyll
229                  DO jk = 1, nksr + 1
230                     zchl3d(:,:,jk) = 0.05 
231                  ENDDO
232               ENDIF
233               !
234               zcoef  = ( 1. - rn_abs ) / 3.e0                        !  equi-partition in R-G-B
235               ze0(:,:,1) = rn_abs  * qsr(:,:)
236               ze1(:,:,1) = zcoef * qsr(:,:)
237               ze2(:,:,1) = zcoef * qsr(:,:)
238               ze3(:,:,1) = zcoef * qsr(:,:)
239               zea(:,:,1) =         qsr(:,:)
240               !
241               DO jk = 2, nksr+1
242                  !
243                  DO jj = 1, jpj                                         ! Separation in R-G-B depending of vertical profile of Chl
244!CDIR NOVERRCHK
245                     DO ji = 1, jpi
246                        zchl = MIN( 10. , MAX( 0.03, zchl3d(ji,jj,jk) ) )
247                        irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
248                        zekb(ji,jj) = rkrgb(1,irgb)
249                        zekg(ji,jj) = rkrgb(2,irgb)
250                        zekr(ji,jj) = rkrgb(3,irgb)
251                     END DO
252                  END DO
253!CDIR NOVERRCHK
254                  DO jj = 1, jpj
255!CDIR NOVERRCHK   
256                     DO ji = 1, jpi
257                        zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r     )
258                        zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekb(ji,jj) )
259                        zc2 = ze2(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekg(ji,jj) )
260                        zc3 = ze3(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekr(ji,jj) )
261                        ze0(ji,jj,jk) = zc0
262                        ze1(ji,jj,jk) = zc1
263                        ze2(ji,jj,jk) = zc2
264                        ze3(ji,jj,jk) = zc3
265                        zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk)
266                     END DO
267                  END DO
268               END DO
269               !
270               DO jk = 1, nksr                                        ! compute and add qsr trend to ta
271                  qsr_hc(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) )
272               END DO
273               zea(:,:,nksr+1:jpk) = 0.e0     ! below 400m set to zero
274               CALL iom_put( 'qsr3d', zea )   ! Shortwave Radiation 3D distribution
275               !
276               IF ( ln_qsr_ice ) THEN    ! store attenuation coefficient of the first ocean level
277!CDIR NOVERRCHK
278                  DO jj = 1, jpj                                         ! Separation in R-G-B depending of the surface Chl
279!CDIR NOVERRCHK
280                     DO ji = 1, jpi
281                        zchl = MIN( 10. , MAX( 0.03, zchl3d(ji,jj,1) ) )
282                        irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
283                        zekb(ji,jj) = rkrgb(1,irgb)
284                        zekg(ji,jj) = rkrgb(2,irgb)
285                        zekr(ji,jj) = rkrgb(3,irgb)
286                     END DO
287                  END DO
288                  !
289                  DO jj = 1, jpj
290                     DO ji = 1, jpi
291                        zc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r     )
292                        zc1 = zcoef  * EXP( - fse3t(ji,jj,1) * zekb(ji,jj) )
293                        zc2 = zcoef  * EXP( - fse3t(ji,jj,1) * zekg(ji,jj) )
294                        zc3 = zcoef  * EXP( - fse3t(ji,jj,1) * zekr(ji,jj) )
295                        fraqsr_1lev(ji,jj) = 1.0 - ( zc0 + zc1 + zc2  + zc3  ) * tmask(ji,jj,2) 
296                     END DO
297                  END DO
298                  !
299               ENDIF
300               !
301            ELSE                                                 !*  Constant Chlorophyll
302               DO jk = 1, nksr
303                  qsr_hc(:,:,jk) =  etot3(:,:,jk) * qsr(:,:)
304               END DO
305               !  store attenuation coefficient of the first ocean level
306               IF( ln_qsr_ice ) THEN
307                  fraqsr_1lev(:,:) = etot3(:,:,1) / r1_rau0_rcp
308               ENDIF
309           ENDIF
310
311         ENDIF
312         !                                                ! ------------------------- !
313         IF( ln_qsr_2bd ) THEN                            !  2 band light penetration !
314            !                                             ! ------------------------- !
315            !
316            IF( lk_vvl ) THEN                                  !* variable volume
317               zz0   =        rn_abs   * r1_rau0_rcp
318               zz1   = ( 1. - rn_abs ) * r1_rau0_rcp
319               DO jk = 1, nksr                    ! solar heat absorbed at T-point in the top 400m
320                  DO jj = 1, jpj
321                     DO ji = 1, jpi
322                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
323                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
324                        qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) ) 
325                     END DO
326                  END DO
327               END DO
328               ! clem: store attenuation coefficient of the first ocean level
329               IF ( ln_qsr_ice ) THEN
330                  DO jj = 1, jpj
331                     DO ji = 1, jpi
332                        zc0 = zz0 * EXP( -fsdepw(ji,jj,1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,1)*xsi1r )
333                        zc1 = zz0 * EXP( -fsdepw(ji,jj,2)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,2)*xsi1r )
334                        fraqsr_1lev(ji,jj) = ( zc0*tmask(ji,jj,1) - zc1*tmask(ji,jj,2) ) / r1_rau0_rcp
335                     END DO
336                  END DO
337               ENDIF
338            ELSE                                               !* constant volume: coef. computed one for all
339               DO jk = 1, nksr
340                  DO jj = 2, jpjm1
341                     DO ji = fs_2, fs_jpim1   ! vector opt.
342                        ! (ISF) no light penetration below the ice shelves         
343                        qsr_hc(ji,jj,jk) =  etot3(ji,jj,jk) * qsr(ji,jj) * tmask(ji,jj,1)
344                     END DO
345                  END DO
346               END DO
347               ! clem: store attenuation coefficient of the first ocean level
348               IF ( ln_qsr_ice ) THEN
349                  fraqsr_1lev(:,:) = etot3(:,:,1) / r1_rau0_rcp
350               ENDIF
351               !
352            ENDIF
353            !
354         ENDIF
355         !
356         !                                        Add to the general trend
357         DO jk = 1, nksr
358            DO jj = 2, jpjm1 
359               DO ji = fs_2, fs_jpim1   ! vector opt.
360                  z1_e3t = zfact / fse3t(ji,jj,jk)
361                  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
362               END DO
363            END DO
364         END DO
365         !
366      ENDIF
367      !
368      IF( lrst_oce ) THEN   !                  Write in the ocean restart file
369         !                                     *******************************
370         IF(lwp .AND. nprint >0) THEN
371            WRITE(numout,*)
372            WRITE(numout,*) 'qsr tracer content forcing field written in ocean restart file ',   &
373            &                    'at it= ', kt,' date= ', ndastp
374            WRITE(numout,*) '~~~~'
375            IF(lflush) CALL flush(numout)
376         ENDIF
377         IF(nn_timing == 2)  CALL timing_start('iom_rstput') 
378         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc      )
379         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev )   ! default definition in sbcssm
380         IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
381         !
382      ENDIF
383
384      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
385         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
386         CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt )
387         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt ) 
388      ENDIF
389      !                       ! print mean trends (used for debugging)
390      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
391      !
392      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        ) 
393      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 
394      !
395      IF( nn_timing == 1 )  CALL timing_stop('tra_qsr')
396      !
397   END SUBROUTINE tra_qsr
398
399
400   SUBROUTINE tra_qsr_init
401      !!----------------------------------------------------------------------
402      !!                  ***  ROUTINE tra_qsr_init  ***
403      !!
404      !! ** Purpose :   Initialization for the penetrative solar radiation
405      !!
406      !! ** Method  :   The profile of solar radiation within the ocean is set
407      !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
408      !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
409      !!      default values correspond to clear water (type I in Jerlov'
410      !!      (1968) classification.
411      !!         called by tra_qsr at the first timestep (nit000)
412      !!
413      !! ** Action  : - initialize rn_si0, rn_si1 and rn_abs
414      !!
415      !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
416      !!----------------------------------------------------------------------
417      !
418      INTEGER  ::   ji, jj, jk                   ! dummy loop indices
419      INTEGER  ::   irgb, ierror, ioptio, nqsr   ! local integer
420      INTEGER  ::   ios                          ! Local integer output status for namelist read
421      REAL(wp) ::   zz0, zc0  , zc1, zcoef       ! local scalars
422      REAL(wp) ::   zz1, zc2  , zc3, zchl        !   -      -
423      REAL(wp), POINTER, DIMENSION(:,:  ) :: zekb, zekg, zekr
424      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea
425      !
426      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files
427      TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read
428      !!
429      NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice,  &
430         &                  nn_chldta, rn_abs, rn_si0, rn_si1
431      !!----------------------------------------------------------------------
432
433      !
434      IF( nn_timing == 1 )  CALL timing_start('tra_qsr_init')
435      !
436      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        ) 
437      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 
438      !
439
440      REWIND( numnam_ref )              ! Namelist namtra_qsr in reference namelist : Ratio and length of penetration
441      READ  ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901)
442901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in reference namelist', lwp )
443
444      REWIND( numnam_cfg )              !  Namelist namtra_qsr in configuration namelist : Ratio and length of penetration
445      READ  ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 )
446902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist', lwp )
447      IF(lwm .AND. nprint > 2) WRITE ( numond, namtra_qsr )
448      !
449      IF(lwp) THEN                ! control print
450         WRITE(numout,*)
451         WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
452         WRITE(numout,*) '~~~~~~~~~~~~'
453         WRITE(numout,*) '   Namelist namtra_qsr : set the parameter of penetration'
454         WRITE(numout,*) '      Light penetration (T) or not (F)         ln_traqsr  = ', ln_traqsr
455         WRITE(numout,*) '      RGB (Red-Green-Blue) light penetration   ln_qsr_rgb = ', ln_qsr_rgb
456         WRITE(numout,*) '      2 band               light penetration   ln_qsr_2bd = ', ln_qsr_2bd
457         WRITE(numout,*) '      bio-model            light penetration   ln_qsr_bio = ', ln_qsr_bio
458         WRITE(numout,*) '      light penetration for ice-model LIM3     ln_qsr_ice = ', ln_qsr_ice
459         WRITE(numout,*) '      RGB : Chl data (=1/2) or cst value (=0)  nn_chldta  = ', nn_chldta
460         WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs = ', rn_abs
461         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0 = ', rn_si0
462         WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1 = ', rn_si1
463         IF(lflush) CALL flush(numout)
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            IF(lflush) CALL flush(numout)
496         ENDIF
497         !
498      ENDIF
499      !                          ! ===================================== !
500      IF( ln_traqsr  ) THEN      !  Initialisation of Light Penetration  ! 
501         !                       ! ===================================== !
502         !
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     )
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            IF(lwp .AND. lflush) CALL flush(numout)
579            !
580         ENDIF
581            !                             ! ---------------------------------- !
582         IF( ln_qsr_2bd ) THEN            !    2 bands    light penetration    !
583            !                             ! ---------------------------------- !
584            !
585            !                                ! level of light extinction
586            nksr = trc_oce_ext_lev( rn_si1, 1.e2 )
587            IF(lwp) THEN
588               WRITE(numout,*)
589               WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
590               IF(lflush) CALL flush(numout)
591            ENDIF
592            !
593            IF( lk_vvl ) THEN                   ! variable volume
594               IF(lwp) WRITE(numout,*) '        key_vvl: light distribution will be computed at each time step'
595               IF(lwp .AND. lflush) CALL flush(numout)
596            ELSE                                ! constant volume: computes one for all
597               zz0 =        rn_abs   * r1_rau0_rcp
598               zz1 = ( 1. - rn_abs ) * r1_rau0_rcp
599               DO jk = 1, nksr                    !*  solar heat absorbed at T-point computed once for all
600                  DO jj = 1, jpj                              ! top 400 meters
601                     DO ji = 1, jpi
602                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
603                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
604                        etot3(ji,jj,jk) = (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) * tmask(ji,jj,1) 
605                     END DO
606                  END DO
607               END DO
608               etot3(:,:,nksr+1:jpk) = 0.e0                   ! below 400m set to zero
609               !
610            ENDIF
611         ENDIF
612         !                       ! ===================================== !
613      ELSE                       !        No light penetration           !                   
614         !                       ! ===================================== !
615         IF(lwp) THEN
616            WRITE(numout,*)
617            WRITE(numout,*) 'tra_qsr_init : NO solar flux penetration'
618            WRITE(numout,*) '~~~~~~~~~~~~'
619            IF(lflush) CALL flush(numout)
620         ENDIF
621      ENDIF
622      !
623      ! initialisation of fraqsr_1lev used in sbcssm
624      IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN
625         IF(nn_timing == 2)  CALL timing_start('iom_rstget')
626         CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev'  , fraqsr_1lev  )
627         IF(nn_timing == 2)  CALL timing_stop('iom_rstget')
628      ELSE
629         fraqsr_1lev(:,:) = 1._wp   ! default definition
630      ENDIF
631      !
632      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        ) 
633      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 
634      !
635      IF( nn_timing == 1 )  CALL timing_stop('tra_qsr_init')
636      !
637   END SUBROUTINE tra_qsr_init
638
639   !!======================================================================
640END MODULE traqsr
Note: See TracBrowser for help on using the repository browser.