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/trunk/src/TOP – NEMO

source: NEMO/trunk/src/TOP/trcopt.F90 @ 14558

Last change on this file since 14558 was 14558, checked in by lovato, 3 years ago

include TOP optical module (#2489)
File size: 15.5 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 ! Light location in the water cell ('center', 'integral')
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      ! No Diurnal cycle PAR
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      !     Weighted broadband attenuation coefficient
148      !     ------------------------------------------
149      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
150         ztmp = ze1(ji,jj,jk)* ekb(ji,jj,jk) + ze2(ji,jj,jk) * ekg(ji,jj,jk) + ze3(ji,jj,jk) * ekr(ji,jj,jk)
151         xeps(ji,jj,jk) = ztmp / e3t(ji,jj,jk,Kmm) / (etot(ji,jj,jk) + rtrn)
152      END_3D
153      !xeps = (ze1(:,:,:) * ekb(:,:,:) + ze2(:,:,:) * ekg(:,:,:) + ze3(:,:,:) * ekr(:,:,:)) / e3t(:,:,:,Kmm) / (etot(:,:,:) + rtrn)
154
155      !     Light at the euphotic depth
156      !     ---------------------------
157      zqsr100 = 0.01 * 3. * zqsr_corr(:,:)
158
159      !     Euphotic depth and level
160      !     ------------------------
161      neln   (:,:) = 1
162      heup   (:,:) = gdepw(:,:,2,Kmm)
163      heup_01(:,:) = gdepw(:,:,2,Kmm)
164      !
165      DO_3D( 1, 1, 1, 1, 2, nksrp )
166        IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >=  zqsr100(ji,jj) )  THEN
167           ! Euphotic level (1st T-level strictly below Euphotic layer)
168           ! NOTE: ensure compatibility with nmld_trc definition in trdmxl_trc
169           neln(ji,jj) = jk+1
170           !
171           ! Euphotic layer depth
172           heup(ji,jj) = gdepw(ji,jj,jk+1,Kmm)
173        ENDIF
174        ! Euphotic layer depth (light level definition)
175        IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.50 )  THEN
176           heup_01(ji,jj) = gdepw(ji,jj,jk+1,Kmm)
177        ENDIF
178      END_3D
179      !
180      heup   (:,:) = MIN( 300., heup   (:,:) )
181      heup_01(:,:) = MIN( 300., heup_01(:,:) )
182      !
183      IF( lk_iomput ) THEN
184         CALL iom_put( "xeps" , xeps(:,:,:) * tmask(:,:,:) )
185         CALL iom_put( "Heup" , heup(:,:  ) * tmask(:,:,1) )
186      ENDIF
187      !
188      IF( ln_timing )   CALL timing_stop('trc_opt')
189      !
190   END SUBROUTINE trc_opt
191
192
193   SUBROUTINE trc_opt_par( kt, zqsr, pe1, pe2, pe3) 
194      !!----------------------------------------------------------------------
195      !!                  ***  routine trc_opt_par  ***
196      !!
197      !! ** purpose :   compute PAR of each wavelength (Red-Green-Blue)
198      !!                from given surface shortwave radiation
199      !!
200      !!----------------------------------------------------------------------
201      INTEGER                         , INTENT(in)      ::   kt                ! ocean time-step
202      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in)      ::   zqsr              ! real shortwave
203      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(out)     ::   pe1 , pe2 , pe3   ! PAR (R-G-B)
204      !
205      INTEGER                       ::   ji, jj, jk        ! dummy loop indices
206      REAL(wp), DIMENSION(jpi,jpj)  ::   we1, we2, we3     ! PAR (R-G-B) at w-level
207      !!----------------------------------------------------------------------
208      pe1(:,:,:) = 0. ; pe2(:,:,:) = 0. ; pe3(:,:,:) = 0.
209      !
210      IF ( TRIM(light_loc) == 'center' ) THEN
211         ! cell-center (t-pivot)
212         pe1(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekb(:,:,1) )
213         pe2(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekg(:,:,1) )
214         pe3(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekr(:,:,1) )
215         !
216         DO_3D( 1, 1, 1, 1, 2, nksrp )
217            pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -0.5 * ( ekb(ji,jj,jk-1) + ekb(ji,jj,jk) ) )
218            pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -0.5 * ( ekg(ji,jj,jk-1) + ekg(ji,jj,jk) ) )
219            pe3(ji,jj,jk) = pe3(ji,jj,jk-1) * EXP( -0.5 * ( ekr(ji,jj,jk-1) + ekr(ji,jj,jk) ) )
220         END_3D
221         !
222      ELSE IF ( TRIM(light_loc) == 'integral' ) THEN
223         ! integrate over cell thickness
224         we1(:,:) = zqsr(:,:)
225         we2(:,:) = zqsr(:,:)
226         we3(:,:) = zqsr(:,:)
227         !
228         DO_3D( 1, 1, 1, 1, 1, nksrp )
229            ! integrate PAR over current t-level
230            pe1(ji,jj,jk) = we1(ji,jj) / (ekb(ji,jj,jk) + rtrn) * (1. - EXP( -ekb(ji,jj,jk) ))
231            pe2(ji,jj,jk) = we2(ji,jj) / (ekg(ji,jj,jk) + rtrn) * (1. - EXP( -ekg(ji,jj,jk) ))
232            pe3(ji,jj,jk) = we3(ji,jj) / (ekr(ji,jj,jk) + rtrn) * (1. - EXP( -ekr(ji,jj,jk) ))
233            ! PAR at next w-level
234            we1(ji,jj) = we1(ji,jj) * EXP( -ekb(ji,jj,jk) )
235            we2(ji,jj) = we2(ji,jj) * EXP( -ekg(ji,jj,jk) )
236            we3(ji,jj) = we3(ji,jj) * EXP( -ekr(ji,jj,jk) )
237         END_3D
238         !
239      ENDIF 
240      !
241      !
242   END SUBROUTINE trc_opt_par
243
244
245   SUBROUTINE trc_opt_sbc( kt )
246      !!----------------------------------------------------------------------
247      !!                  ***  routine trc_opt_sbc  ***
248      !!
249      !! ** purpose :   read and interpolate the variable PAR fraction
250      !!                of shortwave radiation
251      !!
252      !! ** method  :   read the files and interpolate the appropriate variables
253      !!
254      !! ** input   :   external netcdf files
255      !!
256      !!----------------------------------------------------------------------
257      INTEGER, INTENT(in) ::   kt   ! ocean time step
258      !
259      INTEGER  :: ji,jj
260      REAL(wp) :: zcoef
261      !!---------------------------------------------------------------------
262      !
263      IF( ln_timing )  CALL timing_start('trc_opt_sbc')
264      !
265      ! Compute par_varsw at nit000 or only if there is more than 1 time record in par coefficient file
266      IF( ln_varpar ) THEN
267         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_par > 1 ) ) THEN
268            CALL fld_read( kt, 1, sf_par )
269            par_varsw(:,:) = ( sf_par(1)%fnow(:,:,1) ) / 3.0
270         ENDIF
271      ENDIF
272      !
273      IF( ln_timing )  CALL timing_stop('trc_opt_sbc')
274      !
275   END SUBROUTINE trc_opt_sbc
276
277
278   SUBROUTINE trc_opt_ini
279      !!----------------------------------------------------------------------
280      !!                  ***  ROUTINE trc_opt_ini  ***
281      !!
282      !! ** Purpose :   Initialization of tabulated attenuation coefficients
283      !!                and percentage of PAR in Shortwave
284      !!
285      !! ** Input   :   external ascii and netcdf files
286      !!----------------------------------------------------------------------
287      INTEGER :: numpar, ierr, ios   ! Local integer
288      !
289      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
290      TYPE(FLD_N) ::   sn_par                ! informations about the fields to be read
291      !
292      NAMELIST/namtrc_opt/cn_dir, sn_par, ln_varpar, parlux, light_loc
293      !!----------------------------------------------------------------------
294      IF(lwp) THEN
295         WRITE(numout,*)
296         WRITE(numout,*) 'trc_opt_ini : Initialize light module'
297         WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
298      ENDIF
299      READ  ( numnat_ref, namtrc_opt, IOSTAT = ios, ERR = 901)
300901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtrc_opt in reference namelist' )
301      READ  ( numnat_cfg, namtrc_opt, IOSTAT = ios, ERR = 902 )
302902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtrc_opt in configuration namelist' )
303      IF(lwm) WRITE ( numont, namtrc_opt )
304
305      IF(lwp) THEN
306         WRITE(numout,*) '   Namelist : namtrc_opt '
307         WRITE(numout,*) '      PAR as a variable fraction of SW     ln_varpar      = ', ln_varpar
308         WRITE(numout,*) '      Fraction of shortwave as PAR         parlux         = ', parlux
309         WRITE(numout,*) '      Light location in the water cell     light_loc      = ', light_loc
310      ENDIF
311      !
312      ! Variable PAR at the surface of the ocean
313      ! ----------------------------------------
314      IF( ln_varpar ) THEN
315         IF(lwp) WRITE(numout,*)
316         IF(lwp) WRITE(numout,*) '   ==>>>   initialize variable par fraction (ln_varpar=T)'
317         !
318         ALLOCATE( par_varsw(jpi,jpj) )
319         !
320         ALLOCATE( sf_par(1), STAT=ierr )           !* allocate and fill sf_par (forcing structure) with sn_par
321         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'trc_opt_ini: unable to allocate sf_par structure' )
322         !
323         CALL fld_fill( sf_par, (/ sn_par /), cn_dir, 'trc_opt_ini', 'Initialize prescribed PAR forcing ', 'namtrc_opt' )
324                                   ALLOCATE( sf_par(1)%fnow(jpi,jpj,1)   )
325         IF( sn_par%ln_tint )      ALLOCATE( sf_par(1)%fdta(jpi,jpj,1,2) )
326
327         CALL iom_open (  TRIM( sn_par%clname ) , numpar )
328         ntimes_par = iom_getszuld( numpar )   ! get number of record in file
329      ENDIF
330      !
331      CALL trc_oce_rgb( xkrgb )                  ! tabulated attenuation coefficients
332      nksrp = trc_oce_ext_lev( r_si2, 0.33e2 )   ! max level of light extinction (Blue Chl=0.01)
333      !
334      IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksrp, ' ref depth = ', gdepw_1d(nksrp+1), ' m'
335      !
336                         ekr      (:,:,:) = 0._wp
337                         ekb      (:,:,:) = 0._wp
338                         ekg      (:,:,:) = 0._wp
339                         etot     (:,:,:) = 0._wp
340                         etot_ndcy(:,:,:) = 0._wp
341      IF( ln_qsr_bio )   etot3    (:,:,:) = 0._wp
342      !
343   END SUBROUTINE trc_opt_ini
344
345
346   INTEGER FUNCTION trc_opt_alloc()
347      !!----------------------------------------------------------------------
348      !!                     ***  ROUTINE trc_opt_alloc  ***
349      !!----------------------------------------------------------------------
350      !
351      ALLOCATE( ekb(jpi,jpj,jpk), ekr(jpi,jpj,jpk),  &
352                ekg(jpi,jpj,jpk), xeps(jpi,jpj,jpk), STAT= trc_opt_alloc  ) 
353      !
354      IF( trc_opt_alloc /= 0 ) CALL ctl_stop( 'STOP', 'trc_opt_alloc : failed to allocate arrays.' )
355      !
356   END FUNCTION trc_opt_alloc
357
358   !!======================================================================
359END MODULE trcopt
Note: See TracBrowser for help on using the repository browser.