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.
trcopt.F90 in NEMO/branches/2020/dev_r12985_TOP-04_IMMERSE_BGC_interface/src/TOP – NEMO

source: NEMO/branches/2020/dev_r12985_TOP-04_IMMERSE_BGC_interface/src/TOP/trcopt.F90 @ 14362

Last change on this file since 14362 was 14362, checked in by lovato, 4 years ago

trcopt reference version for reference light at top level (#2489)
File size: 14.6 KB
Line 
1MODULE trcopt
2   !!======================================================================
3   !!                         ***  MODULE trcopt  ***
4   !! TOP :   Compute the light in the water column for RGB wavelengths
5   !!======================================================================
6   !! History :  1.0  !  2020     (T. Lovato) Initial code
7   !!----------------------------------------------------------------------
8   !!   trc_opt       : light availability in the water column
9   !!----------------------------------------------------------------------
10   USE trc            ! tracer variables
11   USE oce_trc        ! tracer-ocean share variables
12   USE iom            ! I/O manager
13   USE fldread        ! time interpolation
14
15   IMPLICIT NONE
16   PRIVATE
17
18   PUBLIC   trc_opt       ! called in spefici BGC model routines
19   PUBLIC   trc_opt_ini   ! called in trcini.F90
20   PUBLIC   trc_opt_alloc
21
22   !! * Shared module variables
23
24   LOGICAL  ::   ln_varpar   ! boolean for variable PAR fraction
25   REAL(wp) ::   parlux      ! Fraction of shortwave as PAR
26   CHARACTER (len=25) :: light_loc
27
28   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_par        ! structure of input par
29   INTEGER                              ::   ntimes_par    ! number of time steps in par file
30   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   par_varsw      ! PAR fraction of shortwave
31   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ekb, ekg, ekr  ! wavelength (Red-Green-Blue)
32   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xeps  ! weighted diffusion coefficient
33
34   INTEGER  ::   nksrp   ! levels below which the light cannot penetrate ( depth larger than 391 m)
35
36   ! TL: This array should come directly from traqsr module
37   REAL(wp), DIMENSION(3,61) ::   xkrgb   ! tabulated attenuation coefficients for RGB absorption
38   
39   !! * Substitutions
40#  include "do_loop_substitute.h90"
41#  include "domzgr_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/TOP 4.0 , NEMO Consortium (2020)
44   !! $Id: trcopt.F90 12377 2020-02-12 14:39:06Z acc $
45   !! Software governed by the CeCILL license (see ./LICENSE)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   SUBROUTINE trc_opt( kt, knt, Kbb, Kmm, zchl, ze1, ze2, ze3)
50      !!---------------------------------------------------------------------
51      !!                     ***  ROUTINE trc_opt  ***
52      !!
53      !! ** Purpose : Compute the light availability in the water column
54      !!              depending on depth and chlorophyll concentration
55      !!
56      !! ** Method  : Morel
57      !!---------------------------------------------------------------------
58      INTEGER, INTENT(in) ::   kt, knt   ! ocean time step
59      INTEGER, INTENT(in) ::   Kbb, Kmm  ! time level indices
60      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   zchl  ! chlorophyll field
61      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(out),OPTIONAL :: ze1, ze2, ze3 ! PAR for individual wavelength
62      !
63      INTEGER  ::   ji, jj, jk, irgb
64      REAL(wp) ::   ztmp
65      REAL(wp), DIMENSION(jpi,jpj    ) :: parsw, zqsr100, zqsr_corr
66      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze0
67      !!---------------------------------------------------------------------
68      !
69      IF( ln_timing )   CALL timing_start('trc_opt')
70
71      !     Initialisation of variables used to compute PAR
72      !     -----------------------------------------------
73      ze0(:,:,:) = 0._wp
74      ze1(:,:,:) = 0._wp
75      ze2(:,:,:) = 0._wp
76      ze3(:,:,:) = 0._wp
77
78      !     PAR conversion factor
79      !     --------------------
80      IF( knt == 1 .AND. ln_varpar )   CALL trc_opt_sbc( kt )
81      !
82      IF( ln_varpar ) THEN  ;  parsw(:,:) = par_varsw(:,:)
83      ELSE                  ;  parsw(:,:) = parlux / 3.0
84      ENDIF
85
86      !     Attenuation coef. function of Chlorophyll and wavelength (RGB)
87      !     --------------------------------------------------------------
88      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
89         ztmp = ( zchl(ji,jj,jk) + rtrn ) * 1.e6
90         ztmp = MIN(  10. , MAX( 0.05, ztmp )  )
91         irgb = NINT( 41 + 20.* LOG10( ztmp ) + rtrn )
92         !                                                         
93         ekb(ji,jj,jk) = xkrgb(1,irgb) * e3t(ji,jj,jk,Kmm)
94         ekg(ji,jj,jk) = xkrgb(2,irgb) * e3t(ji,jj,jk,Kmm)
95         ekr(ji,jj,jk) = xkrgb(3,irgb) * e3t(ji,jj,jk,Kmm)
96      END_3D
97
98      !     Heat flux across w-level (used in the dynamics)
99      !     -----------------------------------------------
100      IF( ln_qsr_bio ) THEN
101         !
102         zqsr_corr(:,:) = parsw(:,:) * qsr(:,:)
103         !
104         ze0(:,:,1) = (1._wp - 3._wp * parsw(:,:)) * qsr(:,:)  !  ( 1 - 3 * alpha ) * q
105         ze1(:,:,1) = zqsr_corr(:,:)
106         ze2(:,:,1) = zqsr_corr(:,:)
107         ze3(:,:,1) = zqsr_corr(:,:)
108         !
109         DO jk = 2, nksrp + 1
110            DO_2D(1, 1, 1, 1)
111                  ze0(ji,jj,jk) = ze0(ji,jj,jk-1) * EXP( -e3t(ji,jj,jk-1,Kmm) * (1. / rn_si0) )
112                  ze1(ji,jj,jk) = ze1(ji,jj,jk-1) * EXP( -ekb  (ji,jj,jk-1 )        )
113                  ze2(ji,jj,jk) = ze2(ji,jj,jk-1) * EXP( -ekg  (ji,jj,jk-1 )        )
114                  ze3(ji,jj,jk) = ze3(ji,jj,jk-1) * EXP( -ekr  (ji,jj,jk-1 )        )
115            END_2D
116         END DO
117         !
118         etot3(:,:,1) = qsr(:,:) * tmask(:,:,1)
119         DO jk = 2, nksrp + 1
120            etot3(:,:,jk) =  ( ze0(:,:,jk) + ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk) ) * tmask(:,:,jk)
121         END DO
122         !                                     !  ------------------------
123      ENDIF
124
125      !     Photosynthetically Available Radiation (PAR)
126      !     --------------------------------------------
127      zqsr_corr(:,:) = parsw(:,:) * qsr(:,:) / ( 1.-fr_i(:,:) + rtrn )
128      !
129      CALL trc_opt_par( kt, zqsr_corr, ze1, ze2, ze3 )
130      !
131      DO jk = 1, nksrp
132         etot (:,:,jk) = ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk)
133      ENDDO
134
135      ! weighted broadband diffuse attenuation coefficient
136      WHERE(etot .ne. 0.) &
137         xeps = ( ze1 * ekr + ze2 * ekb + ze3 * ekg ) / e3t(:,:,:,Kmm) / etot
138      ! Compute PAR at cell center (T-level) or integrate over cell depth
139      IF ( TRIM(light_loc) == 'center' ) THEN
140          IF(lwp) WRITE(numout,*) 'trcopt : center'
141          etot = etot * EXP( -xeps * 0.5 * e3t(:,:,:,Kmm))
142      ELSE IF ( TRIM(light_loc) == 'integral' ) THEN
143          IF(lwp) WRITE(numout,*) 'trcopt : integral'
144          WHERE(etot == 0.) &
145          etot = etot / xeps * (1. - EXP(-xeps*e3t(:,:,:,Kmm)))
146      ENDIF
147
148      !
149      ! No Diurnal cycle PAR
150      IF( l_trcdm2dc ) THEN
151         zqsr_corr(:,:) = parsw(:,:) * qsr_mean(:,:) / ( 1.-fr_i(:,:) + rtrn )
152         !
153         CALL trc_opt_par( kt, zqsr_corr, ze1, ze2, ze3 )
154         DO jk = 1, nksrp
155            etot_ndcy(:,:,jk) = ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk)
156         END DO
157      ELSE
158         etot_ndcy(:,:,:) = etot(:,:,:)
159      ENDIF
160
161      !     Light at the euphotic depth
162      !     ---------------------------
163      zqsr100 = 0.01 * 3. * zqsr_corr(:,:)
164
165      !     Euphotic depth and level
166      !     ------------------------
167      neln   (:,:) = 1
168      heup   (:,:) = gdepw(:,:,2,Kmm)
169      heup_01(:,:) = gdepw(:,:,2,Kmm)
170      !
171      DO_3D( 1, 1, 1, 1, 2, nksrp )
172        IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >=  zqsr100(ji,jj) )  THEN
173           ! Euphotic level (1st T-level strictly below Euphotic layer)
174           ! NOTE: ensure compatibility with nmld_trc definition in trdmxl_trc
175           neln(ji,jj) = jk+1
176           !
177           ! Euphotic layer depth
178           heup(ji,jj) = gdepw(ji,jj,jk+1,Kmm)
179        ENDIF
180        ! Euphotic layer depth (light level definition)
181        IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.50 )  THEN
182           heup_01(ji,jj) = gdepw(ji,jj,jk+1,Kmm)
183        ENDIF
184      END_3D
185      !
186      heup   (:,:) = MIN( 300., heup   (:,:) )
187      heup_01(:,:) = MIN( 300., heup_01(:,:) )
188      !
189      IF( lk_iomput ) THEN
190         CALL iom_put( "xeps" , xeps(:,:,:) * tmask(:,:,:) )
191         CALL iom_put( "Heup" , heup(:,:  ) * tmask(:,:,1) )
192      ENDIF
193      !
194      IF( ln_timing )   CALL timing_stop('trc_opt')
195      !
196   END SUBROUTINE trc_opt
197
198
199   SUBROUTINE trc_opt_par( kt, zqsr, pe1, pe2, pe3) 
200      !!----------------------------------------------------------------------
201      !!                  ***  routine trc_opt_par  ***
202      !!
203      !! ** purpose :   compute PAR of each wavelength (Red-Green-Blue)
204      !!                for a given shortwave radiation at w-level
205      !!
206      !!----------------------------------------------------------------------
207      INTEGER                         , INTENT(in)      ::   kt                ! ocean time-step
208      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in)      ::   zqsr              ! real shortwave
209      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(out)     ::   pe1 , pe2 , pe3   ! PAR (R-G-B)
210      !
211      INTEGER    ::   ji, jj, jk                   ! dummy loop indices
212      !!----------------------------------------------------------------------
213      pe1(:,:,:) = 0. ; pe2(:,:,:) = 0. ; pe3(:,:,:) = 0.
214      !
215      pe1(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekb(:,:,1) )
216      pe2(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekg(:,:,1) )
217      pe3(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekr(:,:,1) )
218      !
219      DO_3D( 1, 1, 1, 1, 2, nksrp )
220         pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -0.5 * ( ekb(ji,jj,jk-1) + ekb(ji,jj,jk) ) )
221         pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -0.5 * ( ekg(ji,jj,jk-1) + ekg(ji,jj,jk) ) )
222         pe3(ji,jj,jk) = pe3(ji,jj,jk-1) * EXP( -0.5 * ( ekr(ji,jj,jk-1) + ekr(ji,jj,jk) ) )
223      END_3D
224      !
225      !
226   END SUBROUTINE trc_opt_par
227
228
229   SUBROUTINE trc_opt_sbc( kt )
230      !!----------------------------------------------------------------------
231      !!                  ***  routine trc_opt_sbc  ***
232      !!
233      !! ** purpose :   read and interpolate the variable PAR fraction
234      !!                of shortwave radiation
235      !!
236      !! ** method  :   read the files and interpolate the appropriate variables
237      !!
238      !! ** input   :   external netcdf files
239      !!
240      !!----------------------------------------------------------------------
241      INTEGER, INTENT(in) ::   kt   ! ocean time step
242      !
243      INTEGER  :: ji,jj
244      REAL(wp) :: zcoef
245      !!---------------------------------------------------------------------
246      !
247      IF( ln_timing )  CALL timing_start('trc_opt_sbc')
248      !
249      ! Compute par_varsw at nit000 or only if there is more than 1 time record in par coefficient file
250      IF( ln_varpar ) THEN
251         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_par > 1 ) ) THEN
252            CALL fld_read( kt, 1, sf_par )
253            par_varsw(:,:) = ( sf_par(1)%fnow(:,:,1) ) / 3.0
254         ENDIF
255      ENDIF
256      !
257      IF( ln_timing )  CALL timing_stop('trc_opt_sbc')
258      !
259   END SUBROUTINE trc_opt_sbc
260
261
262   SUBROUTINE trc_opt_ini
263      !!----------------------------------------------------------------------
264      !!                  ***  ROUTINE trc_opt_ini  ***
265      !!
266      !! ** Purpose :   Initialization of tabulated attenuation coefficients
267      !!                and percentage of PAR in Shortwave
268      !!
269      !! ** Input   :   external ascii and netcdf files
270      !!----------------------------------------------------------------------
271      INTEGER :: numpar, ierr, ios   ! Local integer
272      !
273      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
274      TYPE(FLD_N) ::   sn_par                ! informations about the fields to be read
275      !
276      NAMELIST/namtrc_opt/cn_dir, sn_par, ln_varpar, parlux, light_loc
277      !!----------------------------------------------------------------------
278      IF(lwp) THEN
279         WRITE(numout,*)
280         WRITE(numout,*) 'trc_opt_ini : Initialize light module'
281         WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
282      ENDIF
283      READ  ( numnat_ref, namtrc_opt, IOSTAT = ios, ERR = 901)
284901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtrc_opt in reference namelist' )
285      READ  ( numnat_cfg, namtrc_opt, IOSTAT = ios, ERR = 902 )
286902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtrc_opt in configuration namelist' )
287      IF(lwm) WRITE ( numont, namtrc_opt )
288
289      IF(lwp) THEN
290         WRITE(numout,*) '   Namelist : namtrc_opt '
291         WRITE(numout,*) '      PAR as a variable fraction of SW     ln_varpar      = ', ln_varpar
292         WRITE(numout,*) '      Fraction of shortwave as PAR         parlux         = ', parlux
293         WRITE(numout,*) '      Light location in the water cell     light_loc      = ', light_loc
294      ENDIF
295      !
296      ! Variable PAR at the surface of the ocean
297      ! ----------------------------------------
298      IF( ln_varpar ) THEN
299         IF(lwp) WRITE(numout,*)
300         IF(lwp) WRITE(numout,*) '   ==>>>   initialize variable par fraction (ln_varpar=T)'
301         !
302         ALLOCATE( par_varsw(jpi,jpj) )
303         !
304         ALLOCATE( sf_par(1), STAT=ierr )           !* allocate and fill sf_par (forcing structure) with sn_par
305         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'trc_opt_ini: unable to allocate sf_par structure' )
306         !
307         CALL fld_fill( sf_par, (/ sn_par /), cn_dir, 'trc_opt_ini', 'Initialize prescribed PAR forcing ', 'namtrc_opt' )
308                                   ALLOCATE( sf_par(1)%fnow(jpi,jpj,1)   )
309         IF( sn_par%ln_tint )      ALLOCATE( sf_par(1)%fdta(jpi,jpj,1,2) )
310
311         CALL iom_open (  TRIM( sn_par%clname ) , numpar )
312         ntimes_par = iom_getszuld( numpar )   ! get number of record in file
313      ENDIF
314      !
315      CALL trc_oce_rgb( xkrgb )                  ! tabulated attenuation coefficients
316      nksrp = trc_oce_ext_lev( r_si2, 0.33e2 )   ! max level of light extinction (Blue Chl=0.01)
317      !
318      IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksrp, ' ref depth = ', gdepw_1d(nksrp+1), ' m'
319      !
320                         ekr      (:,:,:) = 0._wp
321                         ekb      (:,:,:) = 0._wp
322                         ekg      (:,:,:) = 0._wp
323                         etot     (:,:,:) = 0._wp
324                         etot_ndcy(:,:,:) = 0._wp
325      IF( ln_qsr_bio )   etot3    (:,:,:) = 0._wp
326      !
327   END SUBROUTINE trc_opt_ini
328
329
330   INTEGER FUNCTION trc_opt_alloc()
331      !!----------------------------------------------------------------------
332      !!                     ***  ROUTINE trc_opt_alloc  ***
333      !!----------------------------------------------------------------------
334      !
335      ALLOCATE( ekb(jpi,jpj,jpk), ekr(jpi,jpj,jpk),  &
336                ekg(jpi,jpj,jpk), xeps(jpi,jpj,jpk), STAT= trc_opt_alloc  ) 
337      !
338      IF( trc_opt_alloc /= 0 ) CALL ctl_stop( 'STOP', 'trc_opt_alloc : failed to allocate arrays.' )
339      !
340   END FUNCTION trc_opt_alloc
341
342   !!======================================================================
343END MODULE trcopt
Note: See TracBrowser for help on using the repository browser.