source: branches/CNRS/dev_r4826_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90 @ 5266

Last change on this file since 5266 was 5266, checked in by cetlod, 6 years ago

PISCES_QUOTA : First commits, see ticket #1516

File size: 18.3 KB
Line 
1MODULE p4zopt
2   !!======================================================================
3   !!                         ***  MODULE p4zopt  ***
4   !! TOP - PISCES : Compute the light availability in the water column
5   !!======================================================================
6   !! History :   1.0  !  2004     (O. Aumont) Original code
7   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
8   !!             3.2  !  2009-04  (C. Ethe, G. Madec)  optimisation
9   !!             3.4  !  2011-06  (O. Aumont, C. Ethe) Improve light availability of nano & diat
10   !!----------------------------------------------------------------------
11#if defined  key_pisces || defined key_pisces_quota
12   !!----------------------------------------------------------------------
13   !!   'key_pisces'                                       PISCES bio-model
14   !!----------------------------------------------------------------------
15   !!   p4z_opt       : light availability in the water column
16   !!----------------------------------------------------------------------
17   USE trc            ! tracer variables
18   USE oce_trc        ! tracer-ocean share variables
19   USE sms_pisces     ! Source Minus Sink of PISCES
20   USE iom            ! I/O manager
21   USE fldread         !  time interpolation
22   USE prtctl_trc      !  print control for debugging
23
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   p4z_opt        ! called in p4zbio.F90 module
29   PUBLIC   p4z_opt_init   ! called in trcsms_pisces.F90 module
30   PUBLIC   p4z_opt_alloc
31
32   !! * Shared module variables
33
34   LOGICAL  :: ln_varpar   !: boolean for variable PAR fraction
35   REAL(wp) :: parlux      !: Fraction of shortwave as PAR
36   REAL(wp) :: xparsw                 !: parlux/3
37
38   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_par      ! structure of input par
39   INTEGER , PARAMETER :: nbtimes = 365  !: maximum number of times record in a file
40   INTEGER  :: ntimes_par                ! number of time steps in a file
41   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:) :: par_varsw    !: PAR fraction of shortwave
42
43   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: enano, ediat   !: PAR for nano and diat
44#if defined key_pisces_quota
45   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: epico          !: PAR for pico
46#endif
47   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: emoy           !: averaged PAR in the mixed layer
48
49   INTEGER  ::   nksrp   ! levels below which the light cannot penetrate ( depth larger than 391 m)
50
51   REAL(wp), DIMENSION(3,61), PUBLIC ::   xkrgb   !: tabulated attenuation coefficients for RGB absorption
52   
53   !!* Substitution
54#  include "top_substitute.h90"
55   !!----------------------------------------------------------------------
56   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
57   !! $Id: p4zopt.F90 3160 2011-11-20 14:27:18Z cetlod $
58   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
59   !!----------------------------------------------------------------------
60CONTAINS
61
62   SUBROUTINE p4z_opt( kt, jnt )
63      !!---------------------------------------------------------------------
64      !!                     ***  ROUTINE p4z_opt  ***
65      !!
66      !! ** Purpose :   Compute the light availability in the water column
67      !!              depending on the depth and the chlorophyll concentration
68      !!
69      !! ** Method  : - ???
70      !!---------------------------------------------------------------------
71      !
72      INTEGER, INTENT(in) ::   kt, jnt   ! ocean time step
73      !
74      INTEGER  ::   ji, jj, jk
75      INTEGER  ::   irgb
76      REAL(wp) ::   zchl, zxsi0r
77      REAL(wp) ::   zc0 , zc1 , zc2, zc3, z1_dep
78      REAL(wp), POINTER, DIMENSION(:,:  ) :: zdepmoy, zetmp, zetmp1, zetmp2
79#if defined key_pisces_quota
80      REAL(wp), POINTER, DIMENSION(:,:  ) :: zetmp3
81#endif
82      REAL(wp), POINTER, DIMENSION(:,:,:) :: zekg, zekr, zekb, ze0, ze1, ze2, ze3
83      !!---------------------------------------------------------------------
84      !
85      IF( nn_timing == 1 )  CALL timing_start('p4z_opt')
86      !
87      ! Allocate temporary workspace
88      CALL wrk_alloc( jpi, jpj,      zdepmoy, zetmp, zetmp1, zetmp2 ) 
89      CALL wrk_alloc( jpi, jpj, jpk, zekg, zekr, zekb, ze0, ze1, ze2, ze3 )
90#if defined key_pisces_quota
91      CALL wrk_alloc( jpi, jpj,      zetmp3 )
92#endif
93
94      IF( jnt == 1 .AND. ln_varpar ) CALL p4z_optsbc( kt )
95
96      !     Initialisation of variables used to compute PAR
97      !     -----------------------------------------------
98      ze1(:,:,jpk) = 0._wp
99      ze2(:,:,jpk) = 0._wp
100      ze3(:,:,jpk) = 0._wp
101
102      !                                        !* attenuation coef. function of Chlorophyll and wavelength (Red-Green-Blue)
103      DO jk = 1, jpkm1                         !  --------------------------------------------------------
104!CDIR NOVERRCHK
105         DO jj = 1, jpj
106!CDIR NOVERRCHK
107            DO ji = 1, jpi
108#if defined key_pisces_quota
109               zchl = ( trn(ji,jj,jk,jpnch) + trn(ji,jj,jk,jpdch) + trn(ji,jj,jk,jppch) + rtrn ) * 1.e6
110#else
111               zchl = ( trn(ji,jj,jk,jpnch) + trn(ji,jj,jk,jpdch) + rtrn ) * 1.e6
112#endif
113               zchl = MIN(  10. , MAX( 0.05, zchl )  )
114               irgb = NINT( 41 + 20.* LOG10( zchl ) + rtrn )
115               !                                                         
116               zekb(ji,jj,jk) = xkrgb(1,irgb) * fse3t(ji,jj,jk)
117               zekg(ji,jj,jk) = xkrgb(2,irgb) * fse3t(ji,jj,jk)
118               zekr(ji,jj,jk) = xkrgb(3,irgb) * fse3t(ji,jj,jk)
119            END DO
120         END DO
121      END DO
122
123
124      !                                        !* Photosynthetically Available Radiation (PAR)
125      !                                        !  --------------------------------------
126
127      IF( ln_varpar ) THEN
128         ze1(:,:,1) = par_varsw(:,:) * qsr(:,:) * EXP( -0.5 * zekb(:,:,1) )
129         ze2(:,:,1) = par_varsw(:,:) * qsr(:,:) * EXP( -0.5 * zekg(:,:,1) )
130         ze3(:,:,1) = par_varsw(:,:) * qsr(:,:) * EXP( -0.5 * zekr(:,:,1) )
131      ELSE
132         ze1(:,:,1) = xparsw         * qsr(:,:) * EXP( -0.5 * zekb(:,:,1) )
133         ze2(:,:,1) = xparsw         * qsr(:,:) * EXP( -0.5 * zekg(:,:,1) )
134         ze3(:,:,1) = xparsw         * qsr(:,:) * EXP( -0.5 * zekr(:,:,1) )
135      ENDIF
136
137!CDIR NOVERRCHK
138      DO jj = 1, jpj
139!CDIR NOVERRCHK
140         DO ji = 1, jpi
141            zc1 = ze1(ji,jj,1)
142            zc2 = ze2(ji,jj,1) 
143            zc3 = ze3(ji,jj,1)
144            etot (ji,jj,1) = (       zc1 +        zc2 +       zc3 )
145            enano(ji,jj,1) = ( 2.1 * zc1 + 0.42 * zc2 + 0.4 * zc3 )
146            ediat(ji,jj,1) = ( 1.6 * zc1 + 0.69 * zc2 + 0.7 * zc3 )
147#if defined key_pisces_quota
148            epico(ji,jj,1) = ( 2.1 * zc1 + 0.42 * zc2 + 0.4 * zc3 )
149#endif
150         END DO
151      END DO
152
153   
154      DO jk = 2, nksrp     
155!CDIR NOVERRCHK
156         DO jj = 1, jpj
157!CDIR NOVERRCHK
158            DO ji = 1, jpi
159               zc1 = ze1(ji,jj,jk-1) * EXP( -0.5 * ( zekb(ji,jj,jk-1) + zekb(ji,jj,jk) ) )
160               zc2 = ze2(ji,jj,jk-1) * EXP( -0.5 * ( zekg(ji,jj,jk-1) + zekg(ji,jj,jk) ) )
161               zc3 = ze3(ji,jj,jk-1) * EXP( -0.5 * ( zekr(ji,jj,jk-1) + zekr(ji,jj,jk) ) )
162               ze1  (ji,jj,jk) = zc1
163               ze2  (ji,jj,jk) = zc2
164               ze3  (ji,jj,jk) = zc3
165               etot (ji,jj,jk) = (       zc1 +        zc2 +       zc3 )
166               enano(ji,jj,jk) = ( 2.1 * zc1 + 0.42 * zc2 + 0.4 * zc3 )
167               ediat(ji,jj,jk) = ( 1.6 * zc1 + 0.69 * zc2 + 0.7 * zc3 )
168#if defined key_pisces_quota
169               epico(ji,jj,jk) = ( 2.1 * zc1 + 0.42 * zc2 + 0.4 * zc3 )
170#endif
171            END DO
172         END DO
173      END DO
174 
175      IF( ln_qsr_bio ) THEN                    !* heat flux accros w-level (used in the dynamics)
176         !                                     !  ------------------------
177         zxsi0r = 1.e0 / rn_si0
178         !
179         ze0(:,:,1) = rn_abs * qsr(:,:)
180         !                                                    ! surface value : separation in R-G-B + near surface
181         IF( ln_varpar ) THEN
182            ze0(:,:,1) = ( 1. - 3. * par_varsw(:,:) ) * qsr(:,:)
183            ze1(:,:,1) = par_varsw(:,:)               * qsr(:,:)         
184            ze2(:,:,1) = par_varsw(:,:)               * qsr(:,:)
185            ze3(:,:,1) = par_varsw(:,:)               * qsr(:,:)
186         ELSE
187            ze0(:,:,1) = ( 1. - 3. * xparsw )  * qsr(:,:)
188            ze1(:,:,1) = xparsw                * qsr(:,:)         
189            ze2(:,:,1) = xparsw                * qsr(:,:)
190            ze3(:,:,1) = xparsw                * qsr(:,:)
191         ENDIF
192         etot3(:,:,1) =  qsr(:,:) * tmask(:,:,1)
193         !
194         !
195         DO jk = 2, nksrp + 1
196!CDIR NOVERRCHK
197            DO jj = 1, jpj
198!CDIR NOVERRCHK
199               DO ji = 1, jpi
200                  zc0 = ze0(ji,jj,jk-1) * EXP( -fse3t(ji,jj,jk-1) * zxsi0r )
201                  zc1 = ze1(ji,jj,jk-1) * EXP( -zekb(ji,jj,jk-1 ) )
202                  zc2 = ze2(ji,jj,jk-1) * EXP( -zekg(ji,jj,jk-1 ) )
203                  zc3 = ze3(ji,jj,jk-1) * EXP( -zekr(ji,jj,jk-1 ) )
204                  ze0(ji,jj,jk) = zc0
205                  ze1(ji,jj,jk) = zc1
206                  ze2(ji,jj,jk) = zc2
207                  ze3(ji,jj,jk) = zc3
208                  etot3(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk)
209              END DO
210              !
211            END DO
212            !
213        END DO
214        !
215      ENDIF
216
217      !                                        !* Euphotic depth and level
218      neln(:,:) = 1                            !  ------------------------
219      heup(:,:) = 300.
220
221      DO jk = 2, nksrp
222         DO jj = 1, jpj
223           DO ji = 1, jpi
224              IF( etot(ji,jj,jk) * tmask(ji,jj,jk) >= 0.0043 * qsr(ji,jj) )  THEN
225                 neln(ji,jj) = jk+1                    ! Euphotic level : 1rst T-level strictly below Euphotic layer
226                 !                                     ! nb: ensure the compatibility with nmld_trc definition in trd_mld_trc_zint
227                 heup(ji,jj) = fsdepw(ji,jj,jk+1)      ! Euphotic layer depth
228              ENDIF
229           END DO
230        END DO
231      END DO
232 
233      heup(:,:) = MIN( 300., heup(:,:) )
234
235      !                                        !* mean light over the mixed layer
236      zdepmoy(:,:)   = 0.e0                    !  -------------------------------
237      zetmp  (:,:)   = 0.e0
238      zetmp1 (:,:)   = 0.e0
239      zetmp2 (:,:)   = 0.e0
240#if defined key_pisces_quota
241      zetmp3 (:,:)   = 0.e0
242#endif
243
244      DO jk = 1, nksrp
245!CDIR NOVERRCHK
246         DO jj = 1, jpj
247!CDIR NOVERRCHK
248            DO ji = 1, jpi
249               IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
250                  zetmp  (ji,jj) = zetmp  (ji,jj) + etot (ji,jj,jk) * fse3t(ji,jj,jk)
251                  zetmp1 (ji,jj) = zetmp1 (ji,jj) + enano(ji,jj,jk) * fse3t(ji,jj,jk)
252                  zetmp2 (ji,jj) = zetmp2 (ji,jj) + ediat(ji,jj,jk) * fse3t(ji,jj,jk)
253#if defined key_pisces_quota
254                  zetmp3 (ji,jj) = zetmp3 (ji,jj) + epico(ji,jj,jk) * fse3t(ji,jj,jk)
255#endif
256                  zdepmoy(ji,jj) = zdepmoy(ji,jj) + fse3t(ji,jj,jk)
257               ENDIF
258            END DO
259         END DO
260      END DO
261      !
262      emoy(:,:,:) = etot(:,:,:)
263      !
264      DO jk = 1, nksrp
265!CDIR NOVERRCHK
266         DO jj = 1, jpj
267!CDIR NOVERRCHK
268            DO ji = 1, jpi
269               IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
270                  z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
271                  emoy (ji,jj,jk) = zetmp (ji,jj) * z1_dep
272                  enano(ji,jj,jk) = zetmp1(ji,jj) * z1_dep
273                  ediat(ji,jj,jk) = zetmp2(ji,jj) * z1_dep
274#if defined key_pisces_quota
275                  epico(ji,jj,jk) = zetmp3(ji,jj) * z1_dep
276#endif
277               ENDIF
278            END DO
279         END DO
280      END DO
281
282      IF( ln_diatrc ) THEN        ! save output diagnostics
283        !
284        IF( lk_iomput ) THEN
285           IF( jnt == nrdttrc ) THEN
286              CALL iom_put( "Heup", heup(:,:  ) * tmask(:,:,1) )  ! euphotic layer deptht
287              CALL iom_put( "PAR" , emoy(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation
288           ENDIF
289        ELSE
290           trc2d(:,:,  jp_pcs0_2d + 10) = heup(:,:  ) * tmask(:,:,1) 
291           trc3d(:,:,:,jp_pcs0_3d + 3)  = etot(:,:,:) * tmask(:,:,:)
292        ENDIF
293        !
294      ENDIF
295      !
296      CALL wrk_dealloc( jpi, jpj,      zdepmoy, zetmp, zetmp1, zetmp2 )
297      CALL wrk_dealloc( jpi, jpj, jpk, zekg, zekr, zekb, ze0, ze1, ze2, ze3 )
298#if defined key_pisces_quota
299      CALL wrk_dealloc( jpi, jpj,      zetmp3 )
300#endif
301      !
302      IF( nn_timing == 1 )  CALL timing_stop('p4z_opt')
303      !
304   END SUBROUTINE p4z_opt
305
306   SUBROUTINE p4z_optsbc( kt )
307      !!----------------------------------------------------------------------
308      !!                  ***  routine p4z_optsbc  ***
309      !!
310      !! ** purpose :   read and interpolate the variable PAR fraction
311      !!                of shortwave radiation
312      !!
313      !! ** method  :   read the files and interpolate the appropriate variables
314      !!
315      !! ** input   :   external netcdf files
316      !!
317      !!----------------------------------------------------------------------
318      !! * arguments
319      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
320
321      !! * local declarations
322      INTEGER  :: ji,jj
323      REAL(wp) :: zcoef
324      !!---------------------------------------------------------------------
325      !
326      IF( nn_timing == 1 )  CALL timing_start('p4z_optsbc')
327      !
328      ! Compute par_varsw at nit000 or only if there is more than 1 time record in par coefficient file
329      IF( ln_varpar ) THEN
330         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_par > 1 ) ) THEN
331            CALL fld_read( kt, 1, sf_par )
332            par_varsw(:,:) = ( sf_par(1)%fnow(:,:,1) )/3.0
333         ENDIF
334      ENDIF
335      !
336      IF( nn_timing == 1 )  CALL timing_stop('p4z_optsbc')
337      !
338   END SUBROUTINE p4z_optsbc
339
340   SUBROUTINE p4z_opt_init
341      !!----------------------------------------------------------------------
342      !!                  ***  ROUTINE p4z_opt_init  ***
343      !!
344      !! ** Purpose :   Initialization of tabulated attenuation coef
345      !!                and of the percentage of PAR in Shortwave
346      !!
347      !! ** Input   :   external ascii and netcdf files
348      !!----------------------------------------------------------------------
349      !
350      INTEGER :: numpar
351      INTEGER :: ierr
352      INTEGER :: ios                 ! Local integer output status for namelist read
353      REAL(wp), DIMENSION(nbtimes) :: zsteps                 ! times records
354      !
355      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
356      TYPE(FLD_N) ::   sn_par                ! informations about the fields to be read
357      !
358      NAMELIST/nampisopt/cn_dir, sn_par, ln_varpar, parlux
359
360      !!----------------------------------------------------------------------
361
362      IF( nn_timing == 1 )  CALL timing_start('p4z_opt_init')
363
364      REWIND( numnatp_ref )              ! Namelist nampisopt in reference namelist : Pisces attenuation coef. and PAR
365      READ  ( numnatp_ref, nampisopt, IOSTAT = ios, ERR = 901)
366901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisopt in reference namelist', lwp )
367
368      REWIND( numnatp_cfg )              ! Namelist nampisopt in configuration namelist : Pisces attenuation coef. and PAR
369      READ  ( numnatp_cfg, nampisopt, IOSTAT = ios, ERR = 902 )
370902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisopt in configuration namelist', lwp )
371      IF(lwm) WRITE ( numonp, nampisopt )
372
373      IF(lwp) THEN
374         WRITE(numout,*) ' '
375         WRITE(numout,*) ' namelist : nampisopt '
376         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~ '
377         WRITE(numout,*) '    PAR as a variable fraction of SW     ln_varpar      = ', ln_varpar
378         WRITE(numout,*) '    Default value for the PAR fraction   parlux         = ', parlux
379      ENDIF
380      !
381      xparsw = parlux / 3.0
382      !
383      ! Variable PAR at the surface of the ocean
384      ! ----------------------------------------
385      IF( ln_varpar ) THEN
386         IF(lwp) WRITE(numout,*) '    initialize variable par fraction '
387         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
388         !
389         ALLOCATE( par_varsw(jpi,jpj) )
390         !
391         ALLOCATE( sf_par(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst
392         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_opt_init: unable to allocate sf_par structure' )
393         !
394         CALL fld_fill( sf_par, (/ sn_par /), cn_dir, 'p4z_opt_init', 'Variable PAR fraction ', 'nampisopt' )
395                                   ALLOCATE( sf_par(1)%fnow(jpi,jpj,1)   )
396         IF( sn_par%ln_tint )      ALLOCATE( sf_par(1)%fdta(jpi,jpj,1,2) )
397
398         CALL iom_open (  TRIM( sn_par%clname ) , numpar )
399         CALL iom_gettime( numpar, zsteps, kntime=ntimes_par)  ! get number of record in file
400      ENDIF
401      !
402      CALL trc_oce_rgb( xkrgb )                  ! tabulated attenuation coefficients
403      nksrp = trc_oce_ext_lev( r_si2, 0.33e2 )     ! max level of light extinction (Blue Chl=0.01)
404      !
405      IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksrp, ' ref depth = ', gdepw_1d(nksrp+1), ' m'
406      !
407                         etot (:,:,:) = 0._wp
408                         enano(:,:,:) = 0._wp
409#if defined key_pisces_quota
410                         epico(:,:,:) = 0._wp
411#endif
412                         ediat(:,:,:) = 0._wp
413      IF( ln_qsr_bio )   etot3(:,:,:) = 0._wp
414      !
415      IF( nn_timing == 1 )  CALL timing_stop('p4z_opt_init')
416      !
417   END SUBROUTINE p4z_opt_init
418
419
420   INTEGER FUNCTION p4z_opt_alloc()
421      !!----------------------------------------------------------------------
422      !!                     ***  ROUTINE p4z_opt_alloc  ***
423      !!----------------------------------------------------------------------
424#if defined key_pisces_quota
425     ALLOCATE( enano(jpi,jpj,jpk), ediat(jpi,jpj,jpk), emoy (jpi,jpj,jpk), epico(jpi,jpj,jpk),   &
426      &         STAT=p4z_opt_alloc )
427#else
428      ALLOCATE( enano(jpi,jpj,jpk), ediat(jpi,jpj,jpk), emoy (jpi,jpj,jpk),   &
429      &         STAT=p4z_opt_alloc ) 
430#endif
431         !
432      IF( p4z_opt_alloc /= 0 ) CALL ctl_warn('p4z_opt_alloc : failed to allocate arrays.')
433      !
434   END FUNCTION p4z_opt_alloc
435
436#else
437   !!----------------------------------------------------------------------
438   !!  Dummy module :                                   No PISCES bio-model
439   !!----------------------------------------------------------------------
440CONTAINS
441   SUBROUTINE p4z_opt                   ! Empty routine
442   END SUBROUTINE p4z_opt
443#endif 
444
445   !!======================================================================
446END MODULE  p4zopt
Note: See TracBrowser for help on using the repository browser.