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/2011/dev_r2787_PISCES_improvment/NEMOGCM/NEMO/TOP_SRC/PISCES – NEMO

source: branches/2011/dev_r2787_PISCES_improvment/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zopt.F90 @ 2823

Last change on this file since 2823 was 2823, checked in by cetlod, 13 years ago

Add new parameterisation in PISCES, see ticket #854
  • Property svn:keywords set to Id
File size: 11.8 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
22   IMPLICIT NONE
23   PRIVATE
24
25   PUBLIC   p4z_opt        ! called in p4zbio.F90 module
26   PUBLIC   p4z_opt_init   ! called in trcsms_pisces.F90 module
27   PUBLIC   p4z_opt_alloc
28
29   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: etot, enano, ediat   !: PAR for phyto, nano and diat
30   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: emoy                 !: averaged PAR in the mixed layer
31
32   INTEGER  ::   nksrp   ! levels below which the light cannot penetrate ( depth larger than 391 m)
33   REAL(wp) ::   parlux = 0.43_wp / 3._wp
34
35   REAL(wp), DIMENSION(3,61), PUBLIC ::   xkrgb   !: tabulated attenuation coefficients for RGB absorption
36   
37   !!* Substitution
38#  include "top_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
41   !! $Id$
42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE p4z_opt( kt, jnt )
47      !!---------------------------------------------------------------------
48      !!                     ***  ROUTINE p4z_opt  ***
49      !!
50      !! ** Purpose :   Compute the light availability in the water column
51      !!              depending on the depth and the chlorophyll concentration
52      !!
53      !! ** Method  : - ???
54      !!---------------------------------------------------------------------
55      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
56      USE wrk_nemo, ONLY:   zdepmoy => wrk_2d_1 , zetmp  => wrk_2d_2
57      USE wrk_nemo, ONLY:   zetmp1  => wrk_2d_3 , zetmp2 => wrk_2d_4
58      USE wrk_nemo, ONLY:   zekg    => wrk_3d_2 , zekr   => wrk_3d_3 , zekb => wrk_3d_4
59      USE wrk_nemo, ONLY:   ze0     => wrk_3d_5 , ze1    => wrk_3d_6
60      USE wrk_nemo, ONLY:   ze2     => wrk_3d_7 , ze3    => wrk_3d_8
61      !
62      INTEGER, INTENT(in) ::   kt, jnt   ! ocean time step
63      !
64      INTEGER  ::   ji, jj, jk
65      INTEGER  ::   irgb
66      REAL(wp) ::   zchl, zxsi0r
67      REAL(wp) ::   zc0 , zc1 , zc2, zc3, z1_dep
68      !!---------------------------------------------------------------------
69
70      IF(  wrk_in_use(2, 1,2,3,4)   .OR.   wrk_in_use(3, 2,3,4,5,6,7,8)   ) THEN
71         CALL ctl_stop('p4z_opt: requested workspace arrays unavailable')   ;   RETURN
72      ENDIF
73
74      !     Initialisation of variables used to compute PAR
75      !     -----------------------------------------------
76      ze1 (:,:,jpk) = 0._wp
77      ze2 (:,:,jpk) = 0._wp
78      ze3 (:,:,jpk) = 0._wp
79
80      !                                        !* attenuation coef. function of Chlorophyll and wavelength (Red-Green-Blue)
81      DO jk = 1, jpkm1                         !  --------------------------------------------------------
82!CDIR NOVERRCHK
83         DO jj = 1, jpj
84!CDIR NOVERRCHK
85            DO ji = 1, jpi
86               zchl = ( trn(ji,jj,jk,jpnch) + trn(ji,jj,jk,jpdch) + rtrn ) * 1.e6
87               zchl = MIN(  10. , MAX( 0.05, zchl )  )
88               irgb = NINT( 41 + 20.* LOG10( zchl ) + rtrn )
89               !                                                         
90               zekb(ji,jj,jk) = xkrgb(1,irgb) * fse3t(ji,jj,jk)
91               zekg(ji,jj,jk) = xkrgb(2,irgb) * fse3t(ji,jj,jk)
92               zekr(ji,jj,jk) = xkrgb(3,irgb) * fse3t(ji,jj,jk)
93            END DO
94         END DO
95      END DO
96
97
98      !                                        !* Photosynthetically Available Radiation (PAR)
99      !                                        !  --------------------------------------
100!CDIR NOVERRCHK
101      DO jj = 1, jpj
102!CDIR NOVERRCHK
103         DO ji = 1, jpi
104            zc1 = parlux * qsr(ji,jj) * EXP( -0.5 * zekb(ji,jj,1) )
105            zc2 = parlux * qsr(ji,jj) * EXP( -0.5 * zekg(ji,jj,1) )
106            zc3 = parlux * qsr(ji,jj) * EXP( -0.5 * zekr(ji,jj,1) )
107            ze1  (ji,jj,1) = zc1
108            ze2  (ji,jj,1) = zc2
109            ze3  (ji,jj,1) = zc3
110            etot (ji,jj,1) = (       zc1 +        zc2 +       zc3 )
111            enano(ji,jj,1) = ( 2.1 * zc1 + 0.42 * zc2 + 0.4 * zc3 )
112            ediat(ji,jj,1) = ( 1.6 * zc1 + 0.69 * zc2 + 0.7 * zc3 )
113         END DO
114      END DO
115
116   
117      DO jk = 2, nksrp     
118!CDIR NOVERRCHK
119         DO jj = 1, jpj
120!CDIR NOVERRCHK
121            DO ji = 1, jpi
122               zc1 = ze1(ji,jj,jk-1) * EXP( -0.5 * ( zekb(ji,jj,jk-1) + zekb(ji,jj,jk) ) )
123               zc2 = ze2(ji,jj,jk-1) * EXP( -0.5 * ( zekg(ji,jj,jk-1) + zekg(ji,jj,jk) ) )
124               zc3 = ze3(ji,jj,jk-1) * EXP( -0.5 * ( zekr(ji,jj,jk-1) + zekr(ji,jj,jk) ) )
125               ze1  (ji,jj,jk) = zc1
126               ze2  (ji,jj,jk) = zc2
127               ze3  (ji,jj,jk) = zc3
128               etot (ji,jj,jk) = (       zc1 +        zc2 +       zc3 )
129               enano(ji,jj,jk) = ( 2.1 * zc1 + 0.42 * zc2 + 0.4 * zc3 )
130               ediat(ji,jj,jk) = ( 1.6 * zc1 + 0.69 * zc2 + 0.7 * zc3 )
131            END DO
132         END DO
133      END DO
134
135      IF( ln_qsr_bio ) THEN                    !* heat flux accros w-level (used in the dynamics)
136         !                                     !  ------------------------
137         zxsi0r = 1.e0 / rn_si0
138         !
139         ze0  (:,:,1) = rn_abs * qsr(:,:)
140         ze1  (:,:,1) = parlux * qsr(:,:)             ! surface value : separation in R-G-B + near surface
141         ze2  (:,:,1) = parlux * qsr(:,:)
142         ze3  (:,:,1) = parlux * qsr(:,:)
143         etot3(:,:,1) =          qsr(:,:) * tmask(:,:,1)
144         !
145         DO jk = 2, nksrp + 1
146!CDIR NOVERRCHK
147            DO jj = 1, jpj
148!CDIR NOVERRCHK
149               DO ji = 1, jpi
150                  zc0 = ze0(ji,jj,jk-1) * EXP( -fse3t(ji,jj,jk-1) * zxsi0r )
151                  zc1 = ze1(ji,jj,jk-1) * EXP( -zekb(ji,jj,jk-1 ) )
152                  zc2 = ze2(ji,jj,jk-1) * EXP( -zekg(ji,jj,jk-1 ) )
153                  zc3 = ze3(ji,jj,jk-1) * EXP( -zekr(ji,jj,jk-1 ) )
154                  ze0(ji,jj,jk) = zc0
155                  ze1(ji,jj,jk) = zc1
156                  ze2(ji,jj,jk) = zc2
157                  ze3(ji,jj,jk) = zc3
158                  etot3(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk)
159              END DO
160              !
161            END DO
162            !
163        END DO
164        !
165      ENDIF
166
167      !                                        !* Euphotic depth and level
168      neln(:,:) = 1                            !  ------------------------
169      heup(:,:) = 300.
170
171      DO jk = 2, nksrp
172         DO jj = 1, jpj
173           DO ji = 1, jpi
174              IF( etot(ji,jj,jk) >= 0.0043 * qsr(ji,jj) )  THEN
175                 neln(ji,jj) = jk+1                    ! Euphotic level : 1rst T-level strictly below Euphotic layer
176                 !                                     ! nb: ensure the compatibility with nmld_trc definition in trd_mld_trc_zint
177                 heup(ji,jj) = fsdepw(ji,jj,jk+1)      ! Euphotic layer depth
178              ENDIF
179           END DO
180        END DO
181      END DO
182 
183      heup(:,:) = MIN( 300., heup(:,:) )
184
185      !                                        !* mean light over the mixed layer
186      zdepmoy(:,:)   = 0.e0                    !  -------------------------------
187      zetmp  (:,:)   = 0.e0
188      zetmp1 (:,:)   = 0.e0
189      zetmp2 (:,:)   = 0.e0
190
191      DO jk = 1, nksrp
192!CDIR NOVERRCHK
193         DO jj = 1, jpj
194!CDIR NOVERRCHK
195            DO ji = 1, jpi
196               IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
197                  zetmp  (ji,jj) = zetmp  (ji,jj) + etot (ji,jj,jk) * fse3t(ji,jj,jk)
198                  zetmp1 (ji,jj) = zetmp1 (ji,jj) + enano(ji,jj,jk) * fse3t(ji,jj,jk)
199                  zetmp2 (ji,jj) = zetmp2 (ji,jj) + ediat(ji,jj,jk) * fse3t(ji,jj,jk)
200                  zdepmoy(ji,jj) = zdepmoy(ji,jj) + fse3t(ji,jj,jk)
201               ENDIF
202            END DO
203         END DO
204      END DO
205      !
206      emoy(:,:,:) = etot(:,:,:)
207      !
208      DO jk = 1, nksrp
209!CDIR NOVERRCHK
210         DO jj = 1, jpj
211!CDIR NOVERRCHK
212            DO ji = 1, jpi
213               IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
214                  z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
215                  emoy (ji,jj,jk) = zetmp (ji,jj) * z1_dep
216                  enano(ji,jj,jk) = zetmp1(ji,jj) * z1_dep
217                  ediat(ji,jj,jk) = zetmp2(ji,jj) * z1_dep
218               ENDIF
219            END DO
220         END DO
221      END DO
222
223#if defined key_diatrc
224# if ! defined key_iomput
225      ! save for outputs
226      trc2d(:,:,  jp_pcs0_2d + 10) = heup(:,:  ) * tmask(:,:,1)
227      trc3d(:,:,:,jp_pcs0_3d + 3)  = etot(:,:,:) * tmask(:,:,:)
228# else
229      ! write diagnostics
230      IF( jnt == nrdttrc ) then
231         CALL iom_put( "Heup", heup(:,:  ) * tmask(:,:,1) )  ! euphotic layer deptht
232         CALL iom_put( "PAR" , etot(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation
233      ENDIF
234# endif
235#endif
236      !
237      IF( wrk_not_released(2, 1,2,3,4)           .OR.   &
238          wrk_not_released(3, 2,3,4,5,6,7,8)   )   CALL ctl_stop('p4z_opt: failed to release workspace arrays')
239      !
240   END SUBROUTINE p4z_opt
241
242
243   SUBROUTINE p4z_opt_init
244      !!----------------------------------------------------------------------
245      !!                  ***  ROUTINE p4z_opt_init  ***
246      !!
247      !! ** Purpose :   Initialization of tabulated attenuation coef
248      !!----------------------------------------------------------------------
249      !
250      CALL trc_oce_rgb( xkrgb )                  ! tabulated attenuation coefficients
251      nksrp = trc_oce_ext_lev( r_si2, 0.33e2 )     ! max level of light extinction (Blue Chl=0.01)
252      !
253      IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksrp, ' ref depth = ', gdepw_0(nksrp+1), ' m'
254      !
255                         etot (:,:,:) = 0._wp
256                         enano(:,:,:) = 0._wp
257                         ediat(:,:,:) = 0._wp
258      IF( ln_qsr_bio )   etot3(:,:,:) = 0._wp
259      !
260   END SUBROUTINE p4z_opt_init
261
262
263   INTEGER FUNCTION p4z_opt_alloc()
264      !!----------------------------------------------------------------------
265      !!                     ***  ROUTINE p4z_opt_alloc  ***
266      !!----------------------------------------------------------------------
267      ALLOCATE( etot (jpi,jpj,jpk) , enano(jpi,jpj,jpk) ,     &
268         &      ediat(jpi,jpj,jpk) , emoy (jpi,jpj,jpk) , STAT=p4z_opt_alloc )
269         !
270      IF( p4z_opt_alloc /= 0 ) CALL ctl_warn('p4z_opt_alloc : failed to allocate arrays.')
271      !
272   END FUNCTION p4z_opt_alloc
273
274#else
275   !!----------------------------------------------------------------------
276   !!  Dummy module :                                   No PISCES bio-model
277   !!----------------------------------------------------------------------
278CONTAINS
279   SUBROUTINE p4z_opt                   ! Empty routine
280   END SUBROUTINE p4z_opt
281#endif 
282
283   !!======================================================================
284END MODULE  p4zopt
Note: See TracBrowser for help on using the repository browser.