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

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

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