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/trunk/src/TOP/PISCES/P4Z – NEMO

source: NEMO/trunk/src/TOP/PISCES/P4Z/p4zopt.F90

Last change on this file was 15459, checked in by cetlod, 2 years ago

Various bug fixes and more comments in PISCES routines ; sette test OK in debug mode, nn_hls=1/2, with tiling ; run.stat unchanged ; of course tracer.stat different

  • Property svn:keywords set to Id
File size: 23.8 KB
RevLine 
[3443]1MODULE p4zopt
2   !!======================================================================
3   !!                         ***  MODULE p4zopt  ***
4   !! TOP - PISCES : Compute the light availability in the water column
5   !!======================================================================
[9169]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
[3443]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
[9169]17   USE fldread        !  time interpolation
[13286]18   USE prtctl         !  print control for debugging
[3443]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
[9169]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
[3443]33
34   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_par      ! structure of input par
[7970]35   INTEGER , PARAMETER :: nbtimes = 366  !: maximum number of times record in a file
[3443]36   INTEGER  :: ntimes_par                ! number of time steps in a file
[9169]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)
[3443]39   
[12377]40   !! * Substitutions
41#  include "do_loop_substitute.h90"
[13237]42#  include "domzgr_substitute.h90"
[3443]43   !!----------------------------------------------------------------------
[10067]44   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
[10069]45   !! $Id$
[10068]46   !! Software governed by the CeCILL license (see ./LICENSE)
[3443]47   !!----------------------------------------------------------------------
48CONTAINS
49
[12377]50   SUBROUTINE p4z_opt( kt, knt, Kbb, Kmm )
[3443]51      !!---------------------------------------------------------------------
52      !!                     ***  ROUTINE p4z_opt  ***
53      !!
54      !! ** Purpose :   Compute the light availability in the water column
55      !!              depending on the depth and the chlorophyll concentration
56      !!
57      !! ** Method  : - ???
58      !!---------------------------------------------------------------------
[5385]59      INTEGER, INTENT(in) ::   kt, knt   ! ocean time step
[12377]60      INTEGER, INTENT(in) ::   Kbb, Kmm  ! time level indices
[3443]61      !
62      INTEGER  ::   ji, jj, jk
63      INTEGER  ::   irgb
[5385]64      REAL(wp) ::   zchl
[3443]65      REAL(wp) ::   zc0 , zc1 , zc2, zc3, z1_dep
[9125]66      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zetmp5
67      REAL(wp), DIMENSION(jpi,jpj    ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4
68      REAL(wp), DIMENSION(jpi,jpj    ) :: zqsr100, zqsr_corr
[15459]69      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpar, ze0, ze1, ze2, ze3
[3443]70      !!---------------------------------------------------------------------
71      !
[9169]72      IF( ln_timing )   CALL timing_start('p4z_opt')
[3443]73
[9169]74      IF( knt == 1 .AND. ln_varpar )   CALL p4z_opt_sbc( kt )
[3443]75
76      !     Initialisation of variables used to compute PAR
77      !     -----------------------------------------------
[7753]78      ze1(:,:,:) = 0._wp
79      ze2(:,:,:) = 0._wp
80      ze3(:,:,:) = 0._wp
[15459]81
[7646]82      !
[15459]83      ! Attenuation coef. function of Chlorophyll and wavelength (Red-Green-Blue)
84      ! Thus the light penetration scheme is based on a decomposition of PAR
85      ! into three wave length domains. This was first officially published
86      ! in Lengaigne et al. (2007).
87      ! --------------------------------------------------------
[7646]88      !
[15459]89      ! Computation of the light attenuation parameters based on a
90      ! look-up table
[15090]91      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
[15459]92         zchl =  ( tr(ji,jj,jk,jpnch,Kbb) + tr(ji,jj,jk,jpdch,Kbb) + rtrn ) * 1.e6
93         IF( ln_p5z )   zchl = zchl + tr(ji,jj,jk,jppch,Kbb) * 1.e6
[12377]94         zchl = MIN(  10. , MAX( 0.05, zchl )  )
95         irgb = NINT( 41 + 20.* LOG10( zchl ) + rtrn )
96         !                                                         
[13333]97         ekb(ji,jj,jk) = rkrgb(1,irgb) * e3t(ji,jj,jk,Kmm)
98         ekg(ji,jj,jk) = rkrgb(2,irgb) * e3t(ji,jj,jk,Kmm)
99         ekr(ji,jj,jk) = rkrgb(3,irgb) * e3t(ji,jj,jk,Kmm)
[12377]100      END_3D
[15459]101
102
103      ! Photosynthetically Available Radiation (PAR)
104      ! Two cases are considered in the following :
105      ! (1) An explicit diunal cycle is activated. In that case, mean
106      ! QSR is used as PISCES in its current state has not been parameterized
107      ! for an explicit diurnal cycle
108      ! (2) no diurnal cycle of SW is active and in that case, QSR is used.
109      ! --------------------------------------------
[5385]110      IF( l_trcdm2dc ) THEN                     !  diurnal cycle
[15459]111         IF ( ln_p4z_dcyc ) THEN   ! Diurnal cycle in PISCES
112            !
113            !
114            ! SW over the ice free zone of the grid cell. This assumes that
115            ! SW is zero below sea ice which is a very crude assumption that is
116            ! not fully correct with LIM3 and SI3 but no information is
117            ! currently available to do a better job. SHould be improved in the
118            ! (near) future.
119            zqsr_corr(:,:) = qsr_mean(:,:) / ( 1.-fr_i(:,:) + rtrn )
120            !
121            CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 )
122            !
123            ! Used PAR is computed for each phytoplankton species
124            ! etot_ndcy is PAR at level jk averaged over 24h.
125            ! Due to their size, they have different light absorption characteristics
126            DO jk = 1, nksr
127               etot_ndcy(:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk)
128            END DO
129            !
130            ! SW over the ice free zone of the grid cell. This assumes that
131            ! SW is zero below sea ice which is a very crude assumption that is
132            ! not fully correct with LIM3 and SI3 but no information is
133            ! currently available to do a better job. SHould be improved in the
134            ! (near) future.
135            zqsr_corr(:,:) = qsr(:,:) / ( 1.-fr_i(:,:) + rtrn )
136            !
137            CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3 )
138            !
139            ! Total PAR computation at level jk that includes the diurnal cycle
140            DO jk = 1, nksr
141               etot (:,:,jk) =  ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk)
142               enano(:,:,jk) =  1.85 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.46 * ze3(:,:,jk)
143               ediat(:,:,jk) =  1.62 * ze1(:,:,jk) + 0.74 * ze2(:,:,jk) + 0.63 * ze3(:,:,jk)
144            END DO
145            IF( ln_p5z ) THEN
146               DO jk = 1, nksr
147                  epico  (:,:,jk) =  1.94 * ze1(:,:,jk) + 0.66 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)
148               END DO
149            ENDIF
150
151         ELSE ! No diurnal cycle in PISCES
152
153            !
154            !
155            ! SW over the ice free zone of the grid cell. This assumes that
156            ! SW is zero below sea ice which is a very crude assumption that is
157            ! not fully correct with LIM3 and SI3 but no information is
158            ! currently available to do a better job. SHould be improved in the
159            ! (near) future.
160            zqsr_corr(:,:) = qsr_mean(:,:) / ( 1.-fr_i(:,:) + rtrn )
161            !
162            CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 ) 
163            !
164            ! Used PAR is computed for each phytoplankton species
165            ! etot_ndcy is PAR at level jk averaged over 24h.
166            ! Due to their size, they have different light absorption characteristics
[13333]167            DO jk = 1, nksr     
[15459]168               etot_ndcy(:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk)
169               enano    (:,:,jk) =  1.85 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.46 * ze3(:,:,jk)
170               ediat    (:,:,jk) =  1.62 * ze1(:,:,jk) + 0.74 * ze2(:,:,jk) + 0.63 * ze3(:,:,jk)
[7646]171            END DO
[15459]172            IF( ln_p5z ) THEN
173               DO jk = 1, nksr     
174                  epico  (:,:,jk) =  1.94 * ze1(:,:,jk) + 0.66 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)
175               END DO
176            ENDIF
177            !
178            ! SW over the ice free zone of the grid cell. This assumes that
179            ! SW is zero below sea ice which is a very crude assumption that is
180            ! not fully correct with LIM3 and SI3 but no information is
181            ! currently available to do a better job. SHould be improved in the
182            ! (near) future.
183            zqsr_corr(:,:) = qsr(:,:) / ( 1.-fr_i(:,:) + rtrn )
184            !
185            CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3 ) 
186            !
187            ! Total PAR computation at level jk that includes the diurnal cycle
188            DO jk = 1, nksr     
189               etot(:,:,jk) =  ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk)
190            END DO
[7646]191         ENDIF
[5385]192         !
[15459]193      ELSE   ! no diurnal cycle
[5385]194         !
[6628]195         !
[15459]196         ! SW over the ice free zone of the grid cell. This assumes that
197         ! SW is zero below sea ice which is a very crude assumption that is
198         ! not fully correct with LIM3 and SI3 but no information is
199         ! currently available to do a better job. SHould be improved in the
200         ! (near) future.
[9169]201         zqsr_corr(:,:) = qsr(:,:) / ( 1.-fr_i(:,:) + rtrn )
[5385]202         !
[12377]203         CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100  ) 
[6628]204         !
[15459]205
206         ! Used PAR is computed for each phytoplankton species
207         ! Due to their size, they have different light absorption characteristics
[13333]208         DO jk = 1, nksr     
[15459]209            etot (:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk)    ! Total PAR
210            enano(:,:,jk) =  1.85 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.46 * ze3(:,:,jk)  ! Nanophytoplankton
211            ediat(:,:,jk) =  1.62 * ze1(:,:,jk) + 0.74 * ze2(:,:,jk) + 0.63 * ze3(:,:,jk)  ! Diatoms
[5385]212         END DO
[7646]213         IF( ln_p5z ) THEN
[13333]214            DO jk = 1, nksr     
[15459]215              epico(:,:,jk) =  1.94 * ze1(:,:,jk) + 0.66 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)  ! Picophytoplankton (PISCES-QUOTA)
[7646]216            END DO
217         ENDIF
[7753]218         etot_ndcy(:,:,:) =  etot(:,:,:) 
[3443]219      ENDIF
220
221
[15459]222      ! Biophysical feedback part (computation of vertical penetration of SW)
[3443]223      IF( ln_qsr_bio ) THEN                    !* heat flux accros w-level (used in the dynamics)
224         !                                     !  ------------------------
[12377]225         CALL p4z_opt_par( kt, Kmm, qsr, ze1, ze2, ze3, pe0=ze0 )
[3443]226         !
[7753]227         etot3(:,:,1) =  qsr(:,:) * tmask(:,:,1)
[13333]228         DO jk = 2, nksr + 1
[7753]229            etot3(:,:,jk) =  ( ze0(:,:,jk) + ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk) ) * tmask(:,:,jk)
[5385]230         END DO
231         !                                     !  ------------------------
[3443]232      ENDIF
[15459]233     
234      ! Euphotic depth and level
235      ! Two definitions of the euphotic zone are used here.
236      ! (1) The classical definition based on the relative threshold value
237      ! (2) An alternative definition based on a absolute threshold value.
238      ! -------------------------------------------------------------------
239      neln(:,:) = 1
[12377]240      heup   (:,:) = gdepw(:,:,2,Kmm)
241      heup_01(:,:) = gdepw(:,:,2,Kmm)
[3443]242
[15090]243      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, nksr)
[12377]244        IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >=  zqsr100(ji,jj) )  THEN
245           neln(ji,jj) = jk+1                    ! Euphotic level : 1rst T-level strictly below Euphotic layer
246           !                                     ! nb: ensure the compatibility with nmld_trc definition in trd_mld_trc_zint
247           heup(ji,jj) = gdepw(ji,jj,jk+1,Kmm)     ! Euphotic layer depth
248        ENDIF
[15459]249        IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.10 )  THEN
[12377]250           heup_01(ji,jj) = gdepw(ji,jj,jk+1,Kmm)  ! Euphotic layer depth (light level definition)
251        ENDIF
252      END_3D
[5385]253      !
[15459]254      ! The euphotic depth can not exceed 300 meters.
[7753]255      heup   (:,:) = MIN( 300., heup   (:,:) )
256      heup_01(:,:) = MIN( 300., heup_01(:,:) )
[15459]257     
258      ! Mean PAR over the mixed layer
259      ! -----------------------------
260      zdepmoy(:,:)   = 0.e0             
[7753]261      zetmp1 (:,:)   = 0.e0
262      zetmp2 (:,:)   = 0.e0
[3443]263
[15090]264      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr)
[12377]265         IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
[15459]266            zetmp1 (ji,jj) = zetmp1 (ji,jj) + etot     (ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! Actual PAR for remineralisation
267            zetmp2 (ji,jj) = zetmp2 (ji,jj) + etot_ndcy(ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! Par averaged over 24h for production
[12377]268            zdepmoy(ji,jj) = zdepmoy(ji,jj) +                       e3t(ji,jj,jk,Kmm)
269         ENDIF
270      END_3D
[3443]271      !
[7753]272      emoy(:,:,:) = etot(:,:,:)       ! remineralisation
273      zpar(:,:,:) = etot_ndcy(:,:,:)  ! diagnostic : PAR with no diurnal cycle
[3443]274      !
[15090]275      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr)
[12377]276         IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
277            z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
278            emoy (ji,jj,jk) = zetmp1(ji,jj) * z1_dep
279            zpar (ji,jj,jk) = zetmp2(ji,jj) * z1_dep
280         ENDIF
281      END_3D
[15459]282
283      ! Computation of the mean usable light for the different phytoplankton
284      ! groups based on their absorption characteristics.
[10362]285      zdepmoy(:,:)   = 0.e0
286      zetmp3 (:,:)   = 0.e0
287      zetmp4 (:,:)   = 0.e0
288      !
[15090]289      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr)
[12377]290         IF( gdepw(ji,jj,jk+1,Kmm) <= MIN(hmld(ji,jj), heup_01(ji,jj)) ) THEN
[15459]291            zetmp3 (ji,jj) = zetmp3 (ji,jj) + enano    (ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! Nanophytoplankton
292            zetmp4 (ji,jj) = zetmp4 (ji,jj) + ediat    (ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! Diatoms
[12377]293            zdepmoy(ji,jj) = zdepmoy(ji,jj) +                       e3t(ji,jj,jk,Kmm)
294         ENDIF
295      END_3D
[10362]296      enanom(:,:,:) = enano(:,:,:)
297      ediatm(:,:,:) = ediat(:,:,:)
298      !
[15090]299      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr)
[12377]300         IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
301            z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
302            enanom(ji,jj,jk) = zetmp3(ji,jj) * z1_dep
303            ediatm(ji,jj,jk) = zetmp4(ji,jj) * z1_dep
304         ENDIF
305      END_3D
[10362]306      !
[7646]307      IF( ln_p5z ) THEN
[15459]308         ! Picophytoplankton when using PISCES-QUOTA
[12276]309         ALLOCATE( zetmp5(jpi,jpj) )  ;   zetmp5 (:,:) = 0.e0
[15090]310         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr)
[12377]311            IF( gdepw(ji,jj,jk+1,Kmm) <= MIN(hmld(ji,jj), heup_01(ji,jj)) ) THEN
[15459]312               zetmp5(ji,jj)  = zetmp5 (ji,jj) + epico(ji,jj,jk) * e3t(ji,jj,jk,Kmm)
[12377]313            ENDIF
314         END_3D
[10362]315         !
316         epicom(:,:,:) = epico(:,:,:)
317         !
[15090]318         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr)
[12377]319            IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
320               z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
321               epicom(ji,jj,jk) = zetmp5(ji,jj) * z1_dep
322            ENDIF
323         END_3D
[12276]324         DEALLOCATE( zetmp5 )
[7646]325      ENDIF
[12276]326      !
327      IF( lk_iomput .AND.  knt == nrdttrc ) THEN
328         CALL iom_put( "Heup" , heup(:,:  ) * tmask(:,:,1) )  ! euphotic layer deptht
[14209]329         IF( iom_use( "PAR" ) ) THEN
[14213]330            zpar(:,:,1) = zpar(:,:,1) * ( 1._wp - fr_i(:,:) )
[14209]331            CALL iom_put( "PAR", zpar(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation
332         ENDIF
[3443]333      ENDIF
334      !
[9169]335      IF( ln_timing )   CALL timing_stop('p4z_opt')
[3443]336      !
337   END SUBROUTINE p4z_opt
338
[9124]339
[12377]340   SUBROUTINE p4z_opt_par( kt, Kmm, pqsr, pe1, pe2, pe3, pe0, pqsr100 ) 
[3443]341      !!----------------------------------------------------------------------
[5385]342      !!                  ***  routine p4z_opt_par  ***
[3443]343      !!
[5385]344      !! ** purpose :   compute PAR of each wavelength (Red-Green-Blue)
345      !!                for a given shortwave radiation
346      !!
347      !!----------------------------------------------------------------------
[9169]348      INTEGER                         , INTENT(in)              ::   kt                ! ocean time-step
[12377]349      INTEGER                         , INTENT(in)              ::   Kmm               ! ocean time-index
[9169]350      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   )           ::   pqsr              ! shortwave
351      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)           ::   pe1 , pe2 , pe3   ! PAR ( R-G-B)
352      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout), OPTIONAL ::   pe0               !
353      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(  out), OPTIONAL ::   pqsr100           !
354      !
[5385]355      INTEGER    ::   ji, jj, jk     ! dummy loop indices
[9169]356      REAL(wp), DIMENSION(jpi,jpj) ::  zqsr   ! shortwave
[5385]357      !!----------------------------------------------------------------------
358
359      !  Real shortwave
[7753]360      IF( ln_varpar ) THEN  ;  zqsr(:,:) = par_varsw(:,:) * pqsr(:,:)
361      ELSE                  ;  zqsr(:,:) = xparsw         * pqsr(:,:)
[5385]362      ENDIF
[6945]363     
364      !  Light at the euphotic depth
[9169]365      IF( PRESENT( pqsr100 ) )   pqsr100(:,:) = 0.01 * 3. * zqsr(:,:)
[6945]366
[5385]367      IF( PRESENT( pe0 ) ) THEN     !  W-level
368         !
[7753]369         pe0(:,:,1) = pqsr(:,:) - 3. * zqsr(:,:)    !   ( 1 - 3 * alpha ) * q
370         pe1(:,:,1) = zqsr(:,:)         
371         pe2(:,:,1) = zqsr(:,:)
372         pe3(:,:,1) = zqsr(:,:)
[5385]373         !
[15459]374         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, nksr + 1)
375            pe0(ji,jj,jk) = pe0(ji,jj,jk-1) * EXP( -e3t(ji,jj,jk-1,Kmm) * xsi0r )
376            pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -ekb  (ji,jj,jk-1 )        )
377            pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -ekg  (ji,jj,jk-1 )        )
378            pe3(ji,jj,jk) = pe3(ji,jj,jk-1) * EXP( -ekr  (ji,jj,jk-1 )        )
379        END_3D
[5385]380        !
381      ELSE   ! T- level
382        !
[7753]383        pe1(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekb(:,:,1) )
384        pe2(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekg(:,:,1) )
385        pe3(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekr(:,:,1) )
[5385]386        !
[15090]387        DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, nksr)
[12377]388           pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -0.5 * ( ekb(ji,jj,jk-1) + ekb(ji,jj,jk) ) )
389           pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -0.5 * ( ekg(ji,jj,jk-1) + ekg(ji,jj,jk) ) )
390           pe3(ji,jj,jk) = pe3(ji,jj,jk-1) * EXP( -0.5 * ( ekr(ji,jj,jk-1) + ekr(ji,jj,jk) ) )
391        END_3D
[5385]392        !
393      ENDIF
394      !
395   END SUBROUTINE p4z_opt_par
396
397
398   SUBROUTINE p4z_opt_sbc( kt )
399      !!----------------------------------------------------------------------
400      !!                  ***  routine p4z_opt_sbc  ***
401      !!
[3443]402      !! ** purpose :   read and interpolate the variable PAR fraction
403      !!                of shortwave radiation
404      !!
405      !! ** method  :   read the files and interpolate the appropriate variables
406      !!
407      !! ** input   :   external netcdf files
408      !!
409      !!----------------------------------------------------------------------
[9169]410      INTEGER, INTENT(in) ::   kt   ! ocean time step
[9124]411      !
[3443]412      INTEGER  :: ji,jj
413      REAL(wp) :: zcoef
414      !!---------------------------------------------------------------------
415      !
[9124]416      IF( ln_timing )  CALL timing_start('p4z_optsbc')
[3443]417      !
418      ! Compute par_varsw at nit000 or only if there is more than 1 time record in par coefficient file
419      IF( ln_varpar ) THEN
420         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_par > 1 ) ) THEN
421            CALL fld_read( kt, 1, sf_par )
[5385]422            par_varsw(:,:) = ( sf_par(1)%fnow(:,:,1) ) / 3.0
[3443]423         ENDIF
424      ENDIF
425      !
[9124]426      IF( ln_timing )  CALL timing_stop('p4z_optsbc')
[3443]427      !
[5385]428   END SUBROUTINE p4z_opt_sbc
[3443]429
[9124]430
[3443]431   SUBROUTINE p4z_opt_init
432      !!----------------------------------------------------------------------
433      !!                  ***  ROUTINE p4z_opt_init  ***
434      !!
435      !! ** Purpose :   Initialization of tabulated attenuation coef
436      !!                and of the percentage of PAR in Shortwave
437      !!
438      !! ** Input   :   external ascii and netcdf files
439      !!----------------------------------------------------------------------
[9169]440      INTEGER :: numpar, ierr, ios   ! Local integer
[3443]441      !
442      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
443      TYPE(FLD_N) ::   sn_par                ! informations about the fields to be read
444      !
[15459]445      NAMELIST/nampisopt/cn_dir, sn_par, ln_varpar, parlux, ln_p4z_dcyc
[3443]446      !!----------------------------------------------------------------------
[9169]447      IF(lwp) THEN
448         WRITE(numout,*)
449         WRITE(numout,*) 'p4z_opt_init : '
450         WRITE(numout,*) '~~~~~~~~~~~~ '
451      ENDIF
[4147]452      READ  ( numnatp_ref, nampisopt, IOSTAT = ios, ERR = 901)
[11536]453901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampisopt in reference namelist' )
[4147]454      READ  ( numnatp_cfg, nampisopt, IOSTAT = ios, ERR = 902 )
[11536]455902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisopt in configuration namelist' )
[4624]456      IF(lwm) WRITE ( numonp, nampisopt )
[3443]457
458      IF(lwp) THEN
[9169]459         WRITE(numout,*) '   Namelist : nampisopt '
[15459]460         WRITE(numout,*) '      PAR as a variable fraction of SW       ln_varpar      = ', ln_varpar
461         WRITE(numout,*) '      Default value for the PAR fraction     parlux         = ', parlux
462         WRITE(numout,*) '      Activate the diurnal cycle in PISCES   ln_p4z_dcyc    = ', ln_p4z_dcyc
[3443]463      ENDIF
464      !
465      xparsw = parlux / 3.0
[5385]466      xsi0r  = 1.e0 / rn_si0
[15459]467
468      ! Warning : activate the diurnal cycle with no diurnal cycle in the forcing fields makes no sense
469      ! That does not produce a bug because the model does not use the flag but a warning is necessary
470      ! ----------------------------------------------------------------------------------------------
471      IF ( ln_p4z_dcyc .AND. l_trcdm2dc ) THEN
472         IF (lwp) WRITE(numout,*) 'No diurnal cycle in the PAR forcing field '
473         IF (lwp) WRITE(numout,*) 'Activating the diurnal cycle in PISCES has no effect'
474      ENDIF
[3443]475      !
476      ! Variable PAR at the surface of the ocean
477      ! ----------------------------------------
478      IF( ln_varpar ) THEN
[9169]479         IF(lwp) WRITE(numout,*)
480         IF(lwp) WRITE(numout,*) '   ==>>>   initialize variable par fraction (ln_varpar=T)'
[3443]481         !
482         ALLOCATE( par_varsw(jpi,jpj) )
483         !
484         ALLOCATE( sf_par(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst
485         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_opt_init: unable to allocate sf_par structure' )
486         !
487         CALL fld_fill( sf_par, (/ sn_par /), cn_dir, 'p4z_opt_init', 'Variable PAR fraction ', 'nampisopt' )
488                                   ALLOCATE( sf_par(1)%fnow(jpi,jpj,1)   )
489         IF( sn_par%ln_tint )      ALLOCATE( sf_par(1)%fdta(jpi,jpj,1,2) )
490
491         CALL iom_open (  TRIM( sn_par%clname ) , numpar )
[10522]492         ntimes_par = iom_getszuld( numpar )   ! get number of record in file
[3443]493      ENDIF
494      !
[7753]495                         ekr      (:,:,:) = 0._wp
496                         ekb      (:,:,:) = 0._wp
497                         ekg      (:,:,:) = 0._wp
498                         etot     (:,:,:) = 0._wp
499                         etot_ndcy(:,:,:) = 0._wp
500                         enano    (:,:,:) = 0._wp
501                         ediat    (:,:,:) = 0._wp
502      IF( ln_p5z     )   epico    (:,:,:) = 0._wp
503      IF( ln_qsr_bio )   etot3    (:,:,:) = 0._wp
[3443]504      !
505   END SUBROUTINE p4z_opt_init
506
507
508   INTEGER FUNCTION p4z_opt_alloc()
509      !!----------------------------------------------------------------------
510      !!                     ***  ROUTINE p4z_opt_alloc  ***
511      !!----------------------------------------------------------------------
[7646]512      !
513      ALLOCATE( ekb(jpi,jpj,jpk), ekr(jpi,jpj,jpk),  &
514                ekg(jpi,jpj,jpk), STAT= p4z_opt_alloc  ) 
515      !
[10425]516      IF( p4z_opt_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p4z_opt_alloc : failed to allocate arrays.' )
[3443]517      !
518   END FUNCTION p4z_opt_alloc
519
520   !!======================================================================
[5656]521END MODULE p4zopt
Note: See TracBrowser for help on using the repository browser.