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 branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z – NEMO

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90 @ 4147

Last change on this file since 4147 was 4147, checked in by cetlod, 10 years ago

merge in dev_LOCEAN_2013, the 1st development branch dev_r3853_CNRS9_Confsetting, from its starting point ( r3853 ) on the trunk: see ticket #1169
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: p4zopt.F90 3160 2011-11-20 14:27:18Z cetlod $
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      IF( ln_varpar ) THEN
114         ze1(:,:,1) = par_varsw(:,:) * qsr(:,:) * EXP( -0.5 * zekb(:,:,1) )
115         ze2(:,:,1) = par_varsw(:,:) * qsr(:,:) * EXP( -0.5 * zekg(:,:,1) )
116         ze3(:,:,1) = par_varsw(:,:) * qsr(:,:) * EXP( -0.5 * zekr(:,:,1) )
117      ELSE
118         ze1(:,:,1) = xparsw         * qsr(:,:) * EXP( -0.5 * zekb(:,:,1) )
119         ze2(:,:,1) = xparsw         * qsr(:,:) * EXP( -0.5 * zekg(:,:,1) )
120         ze3(:,:,1) = xparsw         * qsr(:,:) * EXP( -0.5 * zekr(:,:,1) )
121      ENDIF
122
123!CDIR NOVERRCHK
124      DO jj = 1, jpj
125!CDIR NOVERRCHK
126         DO ji = 1, jpi
127            zc1 = ze1(ji,jj,1)
128            zc2 = ze2(ji,jj,1) 
129            zc3 = ze3(ji,jj,1)
130            etot (ji,jj,1) = (       zc1 +        zc2 +       zc3 )
131            enano(ji,jj,1) = ( 2.1 * zc1 + 0.42 * zc2 + 0.4 * zc3 )
132            ediat(ji,jj,1) = ( 1.6 * zc1 + 0.69 * zc2 + 0.7 * zc3 )
133         END DO
134      END DO
135
136   
137      DO jk = 2, nksrp     
138!CDIR NOVERRCHK
139         DO jj = 1, jpj
140!CDIR NOVERRCHK
141            DO ji = 1, jpi
142               zc1 = ze1(ji,jj,jk-1) * EXP( -0.5 * ( zekb(ji,jj,jk-1) + zekb(ji,jj,jk) ) )
143               zc2 = ze2(ji,jj,jk-1) * EXP( -0.5 * ( zekg(ji,jj,jk-1) + zekg(ji,jj,jk) ) )
144               zc3 = ze3(ji,jj,jk-1) * EXP( -0.5 * ( zekr(ji,jj,jk-1) + zekr(ji,jj,jk) ) )
145               ze1  (ji,jj,jk) = zc1
146               ze2  (ji,jj,jk) = zc2
147               ze3  (ji,jj,jk) = zc3
148               etot (ji,jj,jk) = (       zc1 +        zc2 +       zc3 )
149               enano(ji,jj,jk) = ( 2.1 * zc1 + 0.42 * zc2 + 0.4 * zc3 )
150               ediat(ji,jj,jk) = ( 1.6 * zc1 + 0.69 * zc2 + 0.7 * zc3 )
151            END DO
152         END DO
153      END DO
154
155      IF( ln_qsr_bio ) THEN                    !* heat flux accros w-level (used in the dynamics)
156         !                                     !  ------------------------
157         zxsi0r = 1.e0 / rn_si0
158         !
159         ze0(:,:,1) = rn_abs * qsr(:,:)
160         !                                                    ! surface value : separation in R-G-B + near surface
161         IF( ln_varpar ) THEN
162            ze0(:,:,1) = ( 1. - 3. * par_varsw(:,:) ) * qsr(:,:)
163            ze1(:,:,1) = par_varsw(:,:)               * qsr(:,:)         
164            ze2(:,:,1) = par_varsw(:,:)               * qsr(:,:)
165            ze3(:,:,1) = par_varsw(:,:)               * qsr(:,:)
166         ELSE
167            ze0(:,:,1) = ( 1. - 3. * xparsw )  * qsr(:,:)
168            ze1(:,:,1) = xparsw                * qsr(:,:)         
169            ze2(:,:,1) = xparsw                * qsr(:,:)
170            ze3(:,:,1) = xparsw                * qsr(:,:)
171         ENDIF
172         etot3(:,:,1) =  qsr(:,:) * tmask(:,:,1)
173         !
174         !
175         DO jk = 2, nksrp + 1
176!CDIR NOVERRCHK
177            DO jj = 1, jpj
178!CDIR NOVERRCHK
179               DO ji = 1, jpi
180                  zc0 = ze0(ji,jj,jk-1) * EXP( -fse3t(ji,jj,jk-1) * zxsi0r )
181                  zc1 = ze1(ji,jj,jk-1) * EXP( -zekb(ji,jj,jk-1 ) )
182                  zc2 = ze2(ji,jj,jk-1) * EXP( -zekg(ji,jj,jk-1 ) )
183                  zc3 = ze3(ji,jj,jk-1) * EXP( -zekr(ji,jj,jk-1 ) )
184                  ze0(ji,jj,jk) = zc0
185                  ze1(ji,jj,jk) = zc1
186                  ze2(ji,jj,jk) = zc2
187                  ze3(ji,jj,jk) = zc3
188                  etot3(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk)
189              END DO
190              !
191            END DO
192            !
193        END DO
194        !
195      ENDIF
196
197      !                                        !* Euphotic depth and level
198      neln(:,:) = 1                            !  ------------------------
199      heup(:,:) = 300.
200
201      DO jk = 2, nksrp
202         DO jj = 1, jpj
203           DO ji = 1, jpi
204              IF( etot(ji,jj,jk) >= 0.0043 * qsr(ji,jj) )  THEN
205                 neln(ji,jj) = jk+1                    ! Euphotic level : 1rst T-level strictly below Euphotic layer
206                 !                                     ! nb: ensure the compatibility with nmld_trc definition in trd_mld_trc_zint
207                 heup(ji,jj) = fsdepw(ji,jj,jk+1)      ! Euphotic layer depth
208              ENDIF
209           END DO
210        END DO
211      END DO
212 
213      heup(:,:) = MIN( 300., heup(:,:) )
214
215      !                                        !* mean light over the mixed layer
216      zdepmoy(:,:)   = 0.e0                    !  -------------------------------
217      zetmp  (:,:)   = 0.e0
218      zetmp1 (:,:)   = 0.e0
219      zetmp2 (:,:)   = 0.e0
220
221      DO jk = 1, nksrp
222!CDIR NOVERRCHK
223         DO jj = 1, jpj
224!CDIR NOVERRCHK
225            DO ji = 1, jpi
226               IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
227                  zetmp  (ji,jj) = zetmp  (ji,jj) + etot (ji,jj,jk) * fse3t(ji,jj,jk)
228                  zetmp1 (ji,jj) = zetmp1 (ji,jj) + enano(ji,jj,jk) * fse3t(ji,jj,jk)
229                  zetmp2 (ji,jj) = zetmp2 (ji,jj) + ediat(ji,jj,jk) * fse3t(ji,jj,jk)
230                  zdepmoy(ji,jj) = zdepmoy(ji,jj) + fse3t(ji,jj,jk)
231               ENDIF
232            END DO
233         END DO
234      END DO
235      !
236      emoy(:,:,:) = etot(:,:,:)
237      !
238      DO jk = 1, nksrp
239!CDIR NOVERRCHK
240         DO jj = 1, jpj
241!CDIR NOVERRCHK
242            DO ji = 1, jpi
243               IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
244                  z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
245                  emoy (ji,jj,jk) = zetmp (ji,jj) * z1_dep
246                  enano(ji,jj,jk) = zetmp1(ji,jj) * z1_dep
247                  ediat(ji,jj,jk) = zetmp2(ji,jj) * z1_dep
248               ENDIF
249            END DO
250         END DO
251      END DO
252
253      IF( ln_diatrc ) THEN        ! save output diagnostics
254        !
255        IF( lk_iomput ) THEN
256           IF( jnt == nrdttrc ) THEN
257              CALL iom_put( "Heup", heup(:,:  ) * tmask(:,:,1) )  ! euphotic layer deptht
258              CALL iom_put( "PAR" , emoy(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation
259           ENDIF
260        ELSE
261           trc2d(:,:,  jp_pcs0_2d + 10) = heup(:,:  ) * tmask(:,:,1) 
262           trc3d(:,:,:,jp_pcs0_3d + 3)  = etot(:,:,:) * tmask(:,:,:)
263        ENDIF
264        !
265      ENDIF
266      !
267      CALL wrk_dealloc( jpi, jpj,      zdepmoy, zetmp, zetmp1, zetmp2 )
268      CALL wrk_dealloc( jpi, jpj, jpk, zekg, zekr, zekb, ze0, ze1, ze2, ze3 )
269      !
270      IF( nn_timing == 1 )  CALL timing_stop('p4z_opt')
271      !
272   END SUBROUTINE p4z_opt
273
274   SUBROUTINE p4z_optsbc( kt )
275      !!----------------------------------------------------------------------
276      !!                  ***  routine p4z_optsbc  ***
277      !!
278      !! ** purpose :   read and interpolate the variable PAR fraction
279      !!                of shortwave radiation
280      !!
281      !! ** method  :   read the files and interpolate the appropriate variables
282      !!
283      !! ** input   :   external netcdf files
284      !!
285      !!----------------------------------------------------------------------
286      !! * arguments
287      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
288
289      !! * local declarations
290      INTEGER  :: ji,jj
291      REAL(wp) :: zcoef
292      !!---------------------------------------------------------------------
293      !
294      IF( nn_timing == 1 )  CALL timing_start('p4z_optsbc')
295      !
296      ! Compute par_varsw at nit000 or only if there is more than 1 time record in par coefficient file
297      IF( ln_varpar ) THEN
298         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_par > 1 ) ) THEN
299            CALL fld_read( kt, 1, sf_par )
300            par_varsw(:,:) = ( sf_par(1)%fnow(:,:,1) )/3.0
301         ENDIF
302      ENDIF
303      !
304      IF( nn_timing == 1 )  CALL timing_stop('p4z_optsbc')
305      !
306   END SUBROUTINE p4z_optsbc
307
308   SUBROUTINE p4z_opt_init
309      !!----------------------------------------------------------------------
310      !!                  ***  ROUTINE p4z_opt_init  ***
311      !!
312      !! ** Purpose :   Initialization of tabulated attenuation coef
313      !!                and of the percentage of PAR in Shortwave
314      !!
315      !! ** Input   :   external ascii and netcdf files
316      !!----------------------------------------------------------------------
317      !
318      INTEGER :: numpar
319      INTEGER :: ierr
320      INTEGER :: ios                 ! Local integer output status for namelist read
321      REAL(wp), DIMENSION(nbtimes) :: zsteps                 ! times records
322      !
323      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
324      TYPE(FLD_N) ::   sn_par                ! informations about the fields to be read
325      !
326      NAMELIST/nampisopt/cn_dir, sn_par, ln_varpar, parlux
327
328      !!----------------------------------------------------------------------
329
330      IF( nn_timing == 1 )  CALL timing_start('p4z_opt_init')
331
332      REWIND( numnatp_ref )              ! Namelist nampisopt in reference namelist : Pisces attenuation coef. and PAR
333      READ  ( numnatp_ref, nampisopt, IOSTAT = ios, ERR = 901)
334901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisopt in reference namelist', lwp )
335
336      REWIND( numnatp_cfg )              ! Namelist nampisopt in configuration namelist : Pisces attenuation coef. and PAR
337      READ  ( numnatp_cfg, nampisopt, IOSTAT = ios, ERR = 902 )
338902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisopt in configuration namelist', lwp )
339      WRITE ( numonp, nampisopt )
340
341      IF(lwp) THEN
342         WRITE(numout,*) ' '
343         WRITE(numout,*) ' namelist : nampisopt '
344         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~ '
345         WRITE(numout,*) '    PAR as a variable fraction of SW     ln_varpar      = ', ln_varpar
346         WRITE(numout,*) '    Default value for the PAR fraction   parlux         = ', parlux
347      ENDIF
348      !
349      xparsw = parlux / 3.0
350      !
351      ! Variable PAR at the surface of the ocean
352      ! ----------------------------------------
353      IF( ln_varpar ) THEN
354         IF(lwp) WRITE(numout,*) '    initialize variable par fraction '
355         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
356         !
357         ALLOCATE( par_varsw(jpi,jpj) )
358         !
359         ALLOCATE( sf_par(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst
360         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_opt_init: unable to allocate sf_par structure' )
361         !
362         CALL fld_fill( sf_par, (/ sn_par /), cn_dir, 'p4z_opt_init', 'Variable PAR fraction ', 'nampisopt' )
363                                   ALLOCATE( sf_par(1)%fnow(jpi,jpj,1)   )
364         IF( sn_par%ln_tint )      ALLOCATE( sf_par(1)%fdta(jpi,jpj,1,2) )
365
366         CALL iom_open (  TRIM( sn_par%clname ) , numpar )
367         CALL iom_gettime( numpar, zsteps, kntime=ntimes_par)  ! get number of record in file
368      ENDIF
369      !
370      CALL trc_oce_rgb( xkrgb )                  ! tabulated attenuation coefficients
371      nksrp = trc_oce_ext_lev( r_si2, 0.33e2 )     ! max level of light extinction (Blue Chl=0.01)
372      !
373      IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksrp, ' ref depth = ', gdepw_0(nksrp+1), ' m'
374      !
375                         etot (:,:,:) = 0._wp
376                         enano(:,:,:) = 0._wp
377                         ediat(:,:,:) = 0._wp
378      IF( ln_qsr_bio )   etot3(:,:,:) = 0._wp
379      !
380      IF( nn_timing == 1 )  CALL timing_stop('p4z_opt_init')
381      !
382   END SUBROUTINE p4z_opt_init
383
384
385   INTEGER FUNCTION p4z_opt_alloc()
386      !!----------------------------------------------------------------------
387      !!                     ***  ROUTINE p4z_opt_alloc  ***
388      !!----------------------------------------------------------------------
389      ALLOCATE( enano(jpi,jpj,jpk), ediat(jpi,jpj,jpk), emoy (jpi,jpj,jpk), STAT=p4z_opt_alloc ) 
390         !
391      IF( p4z_opt_alloc /= 0 ) CALL ctl_warn('p4z_opt_alloc : failed to allocate arrays.')
392      !
393   END FUNCTION p4z_opt_alloc
394
395#else
396   !!----------------------------------------------------------------------
397   !!  Dummy module :                                   No PISCES bio-model
398   !!----------------------------------------------------------------------
399CONTAINS
400   SUBROUTINE p4z_opt                   ! Empty routine
401   END SUBROUTINE p4z_opt
402#endif 
403
404   !!======================================================================
405END MODULE  p4zopt
Note: See TracBrowser for help on using the repository browser.