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