source: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/p4zopt.F90 @ 2268

Last change on this file since 2268 was 2268, checked in by omamce, 10 years ago

O.M. : temp correction of a bu in p4zopt. Ticket 1355 on NEMO Trac

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