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.
p4zopt.F90 in NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z/p4zopt.F90 @ 11822

Last change on this file since 11822 was 11822, checked in by acc, 4 years ago

Branch 2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps. Sette tested updates to branch to align with trunk changes between 10721 and 11740. Sette tests are passing but results differ from branch before these changes (except for GYRE_PISCES and VORTEX) and branch results already differed from trunk because of algorithmic fixes. Will need more checks to confirm correctness.
  • Property svn:keywords set to Id
File size: 20.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   !!   p4z_opt       : light availability in the water column
12   !!----------------------------------------------------------------------
13   USE trc            ! tracer variables
14   USE oce_trc        ! tracer-ocean share variables
15   USE sms_pisces     ! Source Minus Sink of PISCES
16   USE iom            ! I/O manager
17   USE fldread        !  time interpolation
18   USE prtctl_trc     !  print control for debugging
19
20   IMPLICIT NONE
21   PRIVATE
22
23   PUBLIC   p4z_opt        ! called in p4zbio.F90 module
24   PUBLIC   p4z_opt_init   ! called in trcsms_pisces.F90 module
25   PUBLIC   p4z_opt_alloc
26
27   !! * Shared module variables
28
29   LOGICAL  ::   ln_varpar   ! boolean for variable PAR fraction
30   REAL(wp) ::   parlux      ! Fraction of shortwave as PAR
31   REAL(wp) ::   xparsw      ! parlux/3
32   REAL(wp) ::   xsi0r       ! 1. /rn_si0
33
34   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_par      ! structure of input par
35   INTEGER , PARAMETER :: nbtimes = 366  !: maximum number of times record in a file
36   INTEGER  :: ntimes_par                ! number of time steps in a file
37   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   par_varsw      ! PAR fraction of shortwave
38   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ekb, ekg, ekr  ! wavelength (Red-Green-Blue)
39
40   INTEGER  ::   nksrp   ! levels below which the light cannot penetrate ( depth larger than 391 m)
41
42   REAL(wp), DIMENSION(3,61) ::   xkrgb   ! tabulated attenuation coefficients for RGB absorption
43   
44   !!----------------------------------------------------------------------
45   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
46   !! $Id$
47   !! Software governed by the CeCILL license (see ./LICENSE)
48   !!----------------------------------------------------------------------
49CONTAINS
50
51   SUBROUTINE p4z_opt( kt, knt, Kbb, Kmm )
52      !!---------------------------------------------------------------------
53      !!                     ***  ROUTINE p4z_opt  ***
54      !!
55      !! ** Purpose :   Compute the light availability in the water column
56      !!              depending on the depth and the chlorophyll concentration
57      !!
58      !! ** Method  : - ???
59      !!---------------------------------------------------------------------
60      INTEGER, INTENT(in) ::   kt, knt   ! ocean time step
61      INTEGER, INTENT(in) ::   Kbb, Kmm  ! time level indices
62      !
63      INTEGER  ::   ji, jj, jk
64      INTEGER  ::   irgb
65      REAL(wp) ::   zchl
66      REAL(wp) ::   zc0 , zc1 , zc2, zc3, z1_dep
67      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zetmp5
68      REAL(wp), DIMENSION(jpi,jpj    ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4
69      REAL(wp), DIMENSION(jpi,jpj    ) :: zqsr100, zqsr_corr
70      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpar, ze0, ze1, ze2, ze3, zchl3d
71      !!---------------------------------------------------------------------
72      !
73      IF( ln_timing )   CALL timing_start('p4z_opt')
74      IF( ln_p5z    )   ALLOCATE( zetmp5(jpi,jpj) )
75
76      IF( knt == 1 .AND. ln_varpar )   CALL p4z_opt_sbc( kt )
77
78      !     Initialisation of variables used to compute PAR
79      !     -----------------------------------------------
80      ze1(:,:,:) = 0._wp
81      ze2(:,:,:) = 0._wp
82      ze3(:,:,:) = 0._wp
83      !
84      !                                        !* attenuation coef. function of Chlorophyll and wavelength (Red-Green-Blue)
85      !                                        !  --------------------------------------------------------
86                     zchl3d(:,:,:) = tr(:,:,:,jpnch,Kbb) + tr(:,:,:,jpdch,Kbb)
87      IF( ln_p5z )   zchl3d(:,:,:) = zchl3d(:,:,:)    + tr(:,:,:,jppch,Kbb)
88      !
89      DO jk = 1, jpkm1   
90         DO jj = 1, jpj
91            DO ji = 1, jpi
92               zchl = ( zchl3d(ji,jj,jk) + rtrn ) * 1.e6
93               zchl = MIN(  10. , MAX( 0.05, zchl )  )
94               irgb = NINT( 41 + 20.* LOG10( zchl ) + rtrn )
95               !                                                         
96               ekb(ji,jj,jk) = xkrgb(1,irgb) * e3t(ji,jj,jk,Kmm)
97               ekg(ji,jj,jk) = xkrgb(2,irgb) * e3t(ji,jj,jk,Kmm)
98               ekr(ji,jj,jk) = xkrgb(3,irgb) * e3t(ji,jj,jk,Kmm)
99            END DO
100         END DO
101      END DO
102      !                                        !* Photosynthetically Available Radiation (PAR)
103      !                                        !  --------------------------------------
104      IF( l_trcdm2dc ) THEN                     !  diurnal cycle
105         !
106         zqsr_corr(:,:) = qsr_mean(:,:) / ( 1.-fr_i(:,:) + rtrn )
107         !
108         CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 ) 
109         !
110         DO jk = 1, nksrp     
111            etot_ndcy(:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk)
112            enano    (:,:,jk) =  1.85 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.46 * ze3(:,:,jk)
113            ediat    (:,:,jk) =  1.62 * ze1(:,:,jk) + 0.74 * ze2(:,:,jk) + 0.63 * ze3(:,:,jk)
114         END DO
115         IF( ln_p5z ) THEN
116            DO jk = 1, nksrp     
117              epico  (:,:,jk) =  1.94 * ze1(:,:,jk) + 0.66 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)
118            END DO
119         ENDIF
120         !
121         zqsr_corr(:,:) = qsr(:,:) / ( 1.-fr_i(:,:) + rtrn )
122         !
123         CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3 ) 
124         !
125         DO jk = 1, nksrp     
126            etot(:,:,jk) =  ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk)
127         END DO
128         !
129      ELSE
130         !
131         zqsr_corr(:,:) = qsr(:,:) / ( 1.-fr_i(:,:) + rtrn )
132         !
133         CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100  ) 
134         !
135         DO jk = 1, nksrp     
136            etot (:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk)
137            enano(:,:,jk) =  1.85 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.46 * ze3(:,:,jk)
138            ediat(:,:,jk) =  1.62 * ze1(:,:,jk) + 0.74 * ze2(:,:,jk) + 0.63 * ze3(:,:,jk)
139         END DO
140         IF( ln_p5z ) THEN
141            DO jk = 1, nksrp     
142              epico(:,:,jk) =  1.94 * ze1(:,:,jk) + 0.66 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)
143            END DO
144         ENDIF
145         etot_ndcy(:,:,:) =  etot(:,:,:) 
146      ENDIF
147
148
149      IF( ln_qsr_bio ) THEN                    !* heat flux accros w-level (used in the dynamics)
150         !                                     !  ------------------------
151         CALL p4z_opt_par( kt, Kmm, qsr, ze1, ze2, ze3, pe0=ze0 )
152         !
153         etot3(:,:,1) =  qsr(:,:) * tmask(:,:,1)
154         DO jk = 2, nksrp + 1
155            etot3(:,:,jk) =  ( ze0(:,:,jk) + ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk) ) * tmask(:,:,jk)
156         END DO
157         !                                     !  ------------------------
158      ENDIF
159      !                                        !* Euphotic depth and level
160      neln   (:,:) = 1                            !  ------------------------
161      heup   (:,:) = gdepw(:,:,2,Kmm)
162      heup_01(:,:) = gdepw(:,:,2,Kmm)
163
164      DO jk = 2, nksrp
165         DO jj = 1, jpj
166           DO ji = 1, jpi
167              IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >=  zqsr100(ji,jj) )  THEN
168                 neln(ji,jj) = jk+1                    ! Euphotic level : 1rst T-level strictly below Euphotic layer
169                 !                                     ! nb: ensure the compatibility with nmld_trc definition in trd_mld_trc_zint
170                 heup(ji,jj) = gdepw(ji,jj,jk+1,Kmm)     ! Euphotic layer depth
171              ENDIF
172              IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.50 )  THEN
173                 heup_01(ji,jj) = gdepw(ji,jj,jk+1,Kmm)  ! Euphotic layer depth (light level definition)
174              ENDIF
175           END DO
176        END DO
177      END DO
178      !
179      heup   (:,:) = MIN( 300., heup   (:,:) )
180      heup_01(:,:) = MIN( 300., heup_01(:,:) )
181      !                                        !* mean light over the mixed layer
182      zdepmoy(:,:)   = 0.e0                    !  -------------------------------
183      zetmp1 (:,:)   = 0.e0
184      zetmp2 (:,:)   = 0.e0
185
186      DO jk = 1, nksrp
187         DO jj = 1, jpj
188            DO ji = 1, jpi
189               IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
190                  zetmp1 (ji,jj) = zetmp1 (ji,jj) + etot     (ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! remineralisation
191                  zetmp2 (ji,jj) = zetmp2 (ji,jj) + etot_ndcy(ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! production
192                  zdepmoy(ji,jj) = zdepmoy(ji,jj) +                       e3t(ji,jj,jk,Kmm)
193               ENDIF
194            END DO
195         END DO
196      END DO
197      !
198      emoy(:,:,:) = etot(:,:,:)       ! remineralisation
199      zpar(:,:,:) = etot_ndcy(:,:,:)  ! diagnostic : PAR with no diurnal cycle
200      !
201      DO jk = 1, nksrp
202         DO jj = 1, jpj
203            DO ji = 1, jpi
204               IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
205                  z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
206                  emoy (ji,jj,jk) = zetmp1(ji,jj) * z1_dep
207                  zpar (ji,jj,jk) = zetmp2(ji,jj) * z1_dep
208               ENDIF
209            END DO
210         END DO
211      END DO
212      !
213      zdepmoy(:,:)   = 0.e0
214      zetmp3 (:,:)   = 0.e0
215      zetmp4 (:,:)   = 0.e0
216      !
217      DO jk = 1, nksrp
218         DO jj = 1, jpj
219            DO ji = 1, jpi
220               IF( gdepw(ji,jj,jk+1,Kmm) <= MIN(hmld(ji,jj), heup_01(ji,jj)) ) THEN
221                  zetmp3 (ji,jj) = zetmp3 (ji,jj) + enano    (ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! production
222                  zetmp4 (ji,jj) = zetmp4 (ji,jj) + ediat    (ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! production
223                  zdepmoy(ji,jj) = zdepmoy(ji,jj) +                       e3t(ji,jj,jk,Kmm)
224               ENDIF
225            END DO
226         END DO
227      END DO
228      enanom(:,:,:) = enano(:,:,:)
229      ediatm(:,:,:) = ediat(:,:,:)
230      !
231      DO jk = 1, nksrp
232         DO jj = 1, jpj
233            DO ji = 1, jpi
234               IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
235                  z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
236                  enanom(ji,jj,jk) = zetmp3(ji,jj) * z1_dep
237                  ediatm(ji,jj,jk) = zetmp4(ji,jj) * z1_dep
238               ENDIF
239            END DO
240         END DO
241      END DO
242      !
243      IF( ln_p5z ) THEN
244         zetmp5 (:,:) = 0.e0
245         DO jk = 1, nksrp
246            DO jj = 1, jpj
247               DO ji = 1, jpi
248                  IF( gdepw(ji,jj,jk+1,Kmm) <= MIN(hmld(ji,jj), heup_01(ji,jj)) ) THEN
249                     zetmp5(ji,jj)  = zetmp5 (ji,jj) + epico(ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! production
250                  ENDIF
251               END DO
252            END DO
253         END DO
254         !
255         epicom(:,:,:) = epico(:,:,:)
256         !
257         DO jk = 1, nksrp
258            DO jj = 1, jpj
259               DO ji = 1, jpi
260                  IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
261                     z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
262                     epicom(ji,jj,jk) = zetmp5(ji,jj) * z1_dep
263                  ENDIF
264               END DO
265            END DO
266         END DO
267      ENDIF
268      IF( lk_iomput ) THEN
269        IF( knt == nrdttrc ) THEN
270           IF( iom_use( "Heup"  ) ) CALL iom_put( "Heup" , heup(:,:  ) * tmask(:,:,1) )  ! euphotic layer deptht
271           IF( iom_use( "PARDM" ) ) CALL iom_put( "PARDM", zpar(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation
272           IF( iom_use( "PAR"   ) ) CALL iom_put( "PAR"  , emoy(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation
273        ENDIF
274      ENDIF
275      !
276      IF( ln_p5z    )   DEALLOCATE( zetmp5 )
277      IF( ln_timing )   CALL timing_stop('p4z_opt')
278      !
279   END SUBROUTINE p4z_opt
280
281
282   SUBROUTINE p4z_opt_par( kt, Kmm, pqsr, pe1, pe2, pe3, pe0, pqsr100 ) 
283      !!----------------------------------------------------------------------
284      !!                  ***  routine p4z_opt_par  ***
285      !!
286      !! ** purpose :   compute PAR of each wavelength (Red-Green-Blue)
287      !!                for a given shortwave radiation
288      !!
289      !!----------------------------------------------------------------------
290      INTEGER                         , INTENT(in)              ::   kt                ! ocean time-step
291      INTEGER                         , INTENT(in)              ::   Kmm               ! ocean time-index
292      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   )           ::   pqsr              ! shortwave
293      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)           ::   pe1 , pe2 , pe3   ! PAR ( R-G-B)
294      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout), OPTIONAL ::   pe0               !
295      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(  out), OPTIONAL ::   pqsr100           !
296      !
297      INTEGER    ::   ji, jj, jk     ! dummy loop indices
298      REAL(wp), DIMENSION(jpi,jpj) ::  zqsr   ! shortwave
299      !!----------------------------------------------------------------------
300
301      !  Real shortwave
302      IF( ln_varpar ) THEN  ;  zqsr(:,:) = par_varsw(:,:) * pqsr(:,:)
303      ELSE                  ;  zqsr(:,:) = xparsw         * pqsr(:,:)
304      ENDIF
305     
306      !  Light at the euphotic depth
307      IF( PRESENT( pqsr100 ) )   pqsr100(:,:) = 0.01 * 3. * zqsr(:,:)
308
309      IF( PRESENT( pe0 ) ) THEN     !  W-level
310         !
311         pe0(:,:,1) = pqsr(:,:) - 3. * zqsr(:,:)    !   ( 1 - 3 * alpha ) * q
312         pe1(:,:,1) = zqsr(:,:)         
313         pe2(:,:,1) = zqsr(:,:)
314         pe3(:,:,1) = zqsr(:,:)
315         !
316         DO jk = 2, nksrp + 1
317            DO jj = 1, jpj
318               DO ji = 1, jpi
319                  pe0(ji,jj,jk) = pe0(ji,jj,jk-1) * EXP( -e3t(ji,jj,jk-1,Kmm) * xsi0r )
320                  pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -ekb  (ji,jj,jk-1 )        )
321                  pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -ekg  (ji,jj,jk-1 )        )
322                  pe3(ji,jj,jk) = pe3(ji,jj,jk-1) * EXP( -ekr  (ji,jj,jk-1 )        )
323               END DO
324              !
325            END DO
326            !
327         END DO
328        !
329      ELSE   ! T- level
330        !
331        pe1(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekb(:,:,1) )
332        pe2(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekg(:,:,1) )
333        pe3(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekr(:,:,1) )
334        !
335        DO jk = 2, nksrp     
336           DO jj = 1, jpj
337              DO ji = 1, jpi
338                 pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -0.5 * ( ekb(ji,jj,jk-1) + ekb(ji,jj,jk) ) )
339                 pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -0.5 * ( ekg(ji,jj,jk-1) + ekg(ji,jj,jk) ) )
340                 pe3(ji,jj,jk) = pe3(ji,jj,jk-1) * EXP( -0.5 * ( ekr(ji,jj,jk-1) + ekr(ji,jj,jk) ) )
341              END DO
342           END DO
343        END DO   
344        !
345      ENDIF
346      !
347   END SUBROUTINE p4z_opt_par
348
349
350   SUBROUTINE p4z_opt_sbc( kt )
351      !!----------------------------------------------------------------------
352      !!                  ***  routine p4z_opt_sbc  ***
353      !!
354      !! ** purpose :   read and interpolate the variable PAR fraction
355      !!                of shortwave radiation
356      !!
357      !! ** method  :   read the files and interpolate the appropriate variables
358      !!
359      !! ** input   :   external netcdf files
360      !!
361      !!----------------------------------------------------------------------
362      INTEGER, INTENT(in) ::   kt   ! ocean time step
363      !
364      INTEGER  :: ji,jj
365      REAL(wp) :: zcoef
366      !!---------------------------------------------------------------------
367      !
368      IF( ln_timing )  CALL timing_start('p4z_optsbc')
369      !
370      ! Compute par_varsw at nit000 or only if there is more than 1 time record in par coefficient file
371      IF( ln_varpar ) THEN
372         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_par > 1 ) ) THEN
373            CALL fld_read( kt, 1, sf_par )
374            par_varsw(:,:) = ( sf_par(1)%fnow(:,:,1) ) / 3.0
375         ENDIF
376      ENDIF
377      !
378      IF( ln_timing )  CALL timing_stop('p4z_optsbc')
379      !
380   END SUBROUTINE p4z_opt_sbc
381
382
383   SUBROUTINE p4z_opt_init
384      !!----------------------------------------------------------------------
385      !!                  ***  ROUTINE p4z_opt_init  ***
386      !!
387      !! ** Purpose :   Initialization of tabulated attenuation coef
388      !!                and of the percentage of PAR in Shortwave
389      !!
390      !! ** Input   :   external ascii and netcdf files
391      !!----------------------------------------------------------------------
392      INTEGER :: numpar, ierr, ios   ! Local integer
393      !
394      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
395      TYPE(FLD_N) ::   sn_par                ! informations about the fields to be read
396      !
397      NAMELIST/nampisopt/cn_dir, sn_par, ln_varpar, parlux
398      !!----------------------------------------------------------------------
399      IF(lwp) THEN
400         WRITE(numout,*)
401         WRITE(numout,*) 'p4z_opt_init : '
402         WRITE(numout,*) '~~~~~~~~~~~~ '
403      ENDIF
404      REWIND( numnatp_ref )              ! Namelist nampisopt in reference namelist : Pisces attenuation coef. and PAR
405      READ  ( numnatp_ref, nampisopt, IOSTAT = ios, ERR = 901)
406901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampisopt in reference namelist' )
407      REWIND( numnatp_cfg )              ! Namelist nampisopt in configuration namelist : Pisces attenuation coef. and PAR
408      READ  ( numnatp_cfg, nampisopt, IOSTAT = ios, ERR = 902 )
409902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisopt in configuration namelist' )
410      IF(lwm) WRITE ( numonp, nampisopt )
411
412      IF(lwp) THEN
413         WRITE(numout,*) '   Namelist : nampisopt '
414         WRITE(numout,*) '      PAR as a variable fraction of SW     ln_varpar      = ', ln_varpar
415         WRITE(numout,*) '      Default value for the PAR fraction   parlux         = ', parlux
416      ENDIF
417      !
418      xparsw = parlux / 3.0
419      xsi0r  = 1.e0 / rn_si0
420      !
421      ! Variable PAR at the surface of the ocean
422      ! ----------------------------------------
423      IF( ln_varpar ) THEN
424         IF(lwp) WRITE(numout,*)
425         IF(lwp) WRITE(numout,*) '   ==>>>   initialize variable par fraction (ln_varpar=T)'
426         !
427         ALLOCATE( par_varsw(jpi,jpj) )
428         !
429         ALLOCATE( sf_par(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst
430         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_opt_init: unable to allocate sf_par structure' )
431         !
432         CALL fld_fill( sf_par, (/ sn_par /), cn_dir, 'p4z_opt_init', 'Variable PAR fraction ', 'nampisopt' )
433                                   ALLOCATE( sf_par(1)%fnow(jpi,jpj,1)   )
434         IF( sn_par%ln_tint )      ALLOCATE( sf_par(1)%fdta(jpi,jpj,1,2) )
435
436         CALL iom_open (  TRIM( sn_par%clname ) , numpar )
437         ntimes_par = iom_getszuld( numpar )   ! get number of record in file
438      ENDIF
439      !
440      CALL trc_oce_rgb( xkrgb )                  ! tabulated attenuation coefficients
441      nksrp = trc_oce_ext_lev( r_si2, 0.33e2 )     ! max level of light extinction (Blue Chl=0.01)
442      !
443      IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksrp, ' ref depth = ', gdepw_1d(nksrp+1), ' m'
444      !
445                         ekr      (:,:,:) = 0._wp
446                         ekb      (:,:,:) = 0._wp
447                         ekg      (:,:,:) = 0._wp
448                         etot     (:,:,:) = 0._wp
449                         etot_ndcy(:,:,:) = 0._wp
450                         enano    (:,:,:) = 0._wp
451                         ediat    (:,:,:) = 0._wp
452      IF( ln_p5z     )   epico    (:,:,:) = 0._wp
453      IF( ln_qsr_bio )   etot3    (:,:,:) = 0._wp
454      !
455   END SUBROUTINE p4z_opt_init
456
457
458   INTEGER FUNCTION p4z_opt_alloc()
459      !!----------------------------------------------------------------------
460      !!                     ***  ROUTINE p4z_opt_alloc  ***
461      !!----------------------------------------------------------------------
462      !
463      ALLOCATE( ekb(jpi,jpj,jpk), ekr(jpi,jpj,jpk),  &
464                ekg(jpi,jpj,jpk), STAT= p4z_opt_alloc  ) 
465      !
466      IF( p4z_opt_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p4z_opt_alloc : failed to allocate arrays.' )
467      !
468   END FUNCTION p4z_opt_alloc
469
470   !!======================================================================
471END MODULE p4zopt
Note: See TracBrowser for help on using the repository browser.