source: branches/NERC/dev_r5107_NOC_MEDUSA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90 @ 5712

Last change on this file since 5712 was 5712, checked in by acc, 5 years ago

Branch NERC/dev_r5107_NOC_MEDUSA. Remove extra blanks from the ends of the SVN keyword lines that cause unnecessary conflicts with fcm_make

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