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 @ 13107

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

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