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 trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z – NEMO

source: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90 @ 4624

Last change on this file since 4624 was 4624, checked in by acc, 10 years ago

#1305. Fix slow start-up problems on some systems by introducing and using lwm logical to restrict output of merged namelists to the first (or only) processor. lwm is true only on the first processor regardless of ln_ctl. Small changes to all flavours of nemogcm.F90 are also required to write namctl and namcfg after the call to mynode which now opens output.namelist.dyn and writes nammpp.
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   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
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_mld_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( ln_diatrc ) THEN        ! save output diagnostics
255        !
256        IF( lk_iomput ) THEN
257           IF( jnt == nrdttrc ) THEN
258              CALL iom_put( "Heup", heup(:,:  ) * tmask(:,:,1) )  ! euphotic layer deptht
259              CALL iom_put( "PAR" , emoy(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation
260           ENDIF
261        ELSE
262           trc2d(:,:,  jp_pcs0_2d + 10) = heup(:,:  ) * tmask(:,:,1) 
263           trc3d(:,:,:,jp_pcs0_3d + 3)  = etot(:,:,:) * tmask(:,:,:)
264        ENDIF
265        !
266      ENDIF
267      !
268      CALL wrk_dealloc( jpi, jpj,      zdepmoy, zetmp, zetmp1, zetmp2 )
269      CALL wrk_dealloc( jpi, jpj, jpk, zekg, zekr, zekb, ze0, ze1, ze2, ze3 )
270      !
271      IF( nn_timing == 1 )  CALL timing_stop('p4z_opt')
272      !
273   END SUBROUTINE p4z_opt
274
275   SUBROUTINE p4z_optsbc( kt )
276      !!----------------------------------------------------------------------
277      !!                  ***  routine p4z_optsbc  ***
278      !!
279      !! ** purpose :   read and interpolate the variable PAR fraction
280      !!                of shortwave radiation
281      !!
282      !! ** method  :   read the files and interpolate the appropriate variables
283      !!
284      !! ** input   :   external netcdf files
285      !!
286      !!----------------------------------------------------------------------
287      !! * arguments
288      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
289
290      !! * local declarations
291      INTEGER  :: ji,jj
292      REAL(wp) :: zcoef
293      !!---------------------------------------------------------------------
294      !
295      IF( nn_timing == 1 )  CALL timing_start('p4z_optsbc')
296      !
297      ! Compute par_varsw at nit000 or only if there is more than 1 time record in par coefficient file
298      IF( ln_varpar ) THEN
299         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_par > 1 ) ) THEN
300            CALL fld_read( kt, 1, sf_par )
301            par_varsw(:,:) = ( sf_par(1)%fnow(:,:,1) )/3.0
302         ENDIF
303      ENDIF
304      !
305      IF( nn_timing == 1 )  CALL timing_stop('p4z_optsbc')
306      !
307   END SUBROUTINE p4z_optsbc
308
309   SUBROUTINE p4z_opt_init
310      !!----------------------------------------------------------------------
311      !!                  ***  ROUTINE p4z_opt_init  ***
312      !!
313      !! ** Purpose :   Initialization of tabulated attenuation coef
314      !!                and of the percentage of PAR in Shortwave
315      !!
316      !! ** Input   :   external ascii and netcdf files
317      !!----------------------------------------------------------------------
318      !
319      INTEGER :: numpar
320      INTEGER :: ierr
321      INTEGER :: ios                 ! Local integer output status for namelist read
322      REAL(wp), DIMENSION(nbtimes) :: zsteps                 ! times records
323      !
324      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
325      TYPE(FLD_N) ::   sn_par                ! informations about the fields to be read
326      !
327      NAMELIST/nampisopt/cn_dir, sn_par, ln_varpar, parlux
328
329      !!----------------------------------------------------------------------
330
331      IF( nn_timing == 1 )  CALL timing_start('p4z_opt_init')
332
333      REWIND( numnatp_ref )              ! Namelist nampisopt in reference namelist : Pisces attenuation coef. and PAR
334      READ  ( numnatp_ref, nampisopt, IOSTAT = ios, ERR = 901)
335901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisopt in reference namelist', lwp )
336
337      REWIND( numnatp_cfg )              ! Namelist nampisopt in configuration namelist : Pisces attenuation coef. and PAR
338      READ  ( numnatp_cfg, nampisopt, IOSTAT = ios, ERR = 902 )
339902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisopt in configuration namelist', lwp )
340      IF(lwm) WRITE ( numonp, nampisopt )
341
342      IF(lwp) THEN
343         WRITE(numout,*) ' '
344         WRITE(numout,*) ' namelist : nampisopt '
345         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~ '
346         WRITE(numout,*) '    PAR as a variable fraction of SW     ln_varpar      = ', ln_varpar
347         WRITE(numout,*) '    Default value for the PAR fraction   parlux         = ', parlux
348      ENDIF
349      !
350      xparsw = parlux / 3.0
351      !
352      ! Variable PAR at the surface of the ocean
353      ! ----------------------------------------
354      IF( ln_varpar ) THEN
355         IF(lwp) WRITE(numout,*) '    initialize variable par fraction '
356         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
357         !
358         ALLOCATE( par_varsw(jpi,jpj) )
359         !
360         ALLOCATE( sf_par(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst
361         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_opt_init: unable to allocate sf_par structure' )
362         !
363         CALL fld_fill( sf_par, (/ sn_par /), cn_dir, 'p4z_opt_init', 'Variable PAR fraction ', 'nampisopt' )
364                                   ALLOCATE( sf_par(1)%fnow(jpi,jpj,1)   )
365         IF( sn_par%ln_tint )      ALLOCATE( sf_par(1)%fdta(jpi,jpj,1,2) )
366
367         CALL iom_open (  TRIM( sn_par%clname ) , numpar )
368         CALL iom_gettime( numpar, zsteps, kntime=ntimes_par)  ! get number of record in file
369      ENDIF
370      !
371      CALL trc_oce_rgb( xkrgb )                  ! tabulated attenuation coefficients
372      nksrp = trc_oce_ext_lev( r_si2, 0.33e2 )     ! max level of light extinction (Blue Chl=0.01)
373      !
374      IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksrp, ' ref depth = ', gdepw_1d(nksrp+1), ' m'
375      !
376                         etot (:,:,:) = 0._wp
377                         enano(:,:,:) = 0._wp
378                         ediat(:,:,:) = 0._wp
379      IF( ln_qsr_bio )   etot3(:,:,:) = 0._wp
380      !
381      IF( nn_timing == 1 )  CALL timing_stop('p4z_opt_init')
382      !
383   END SUBROUTINE p4z_opt_init
384
385
386   INTEGER FUNCTION p4z_opt_alloc()
387      !!----------------------------------------------------------------------
388      !!                     ***  ROUTINE p4z_opt_alloc  ***
389      !!----------------------------------------------------------------------
390      ALLOCATE( enano(jpi,jpj,jpk), ediat(jpi,jpj,jpk), emoy (jpi,jpj,jpk), STAT=p4z_opt_alloc ) 
391         !
392      IF( p4z_opt_alloc /= 0 ) CALL ctl_warn('p4z_opt_alloc : failed to allocate arrays.')
393      !
394   END FUNCTION p4z_opt_alloc
395
396#else
397   !!----------------------------------------------------------------------
398   !!  Dummy module :                                   No PISCES bio-model
399   !!----------------------------------------------------------------------
400CONTAINS
401   SUBROUTINE p4z_opt                   ! Empty routine
402   END SUBROUTINE p4z_opt
403#endif 
404
405   !!======================================================================
406END MODULE  p4zopt
Note: See TracBrowser for help on using the repository browser.