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.
Changeset 5989 for branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90 – NEMO

Ignore:
Timestamp:
2015-12-03T09:10:32+01:00 (8 years ago)
Author:
deazer
Message:

Merging TMB and 25h diagnostics to head of trunk
added brief documentation

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90

    r5260 r5989  
    3535   REAL(wp) :: parlux      !: Fraction of shortwave as PAR 
    3636   REAL(wp) :: xparsw                 !: parlux/3 
     37   REAL(wp) :: xsi0r                 !:  1. /rn_si0 
    3738 
    3839   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_par      ! structure of input par 
     
    4243 
    4344   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: enano, ediat   !: PAR for phyto, nano and diat  
     45   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: etot_ndcy      !: PAR over 24h in case of diurnal cycle 
    4446   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: emoy           !: averaged PAR in the mixed layer 
     47   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ekb, ekg, ekr  !: wavelength (Red-Green-Blue) 
    4548 
    4649   INTEGER  ::   nksrp   ! levels below which the light cannot penetrate ( depth larger than 391 m) 
     
    4851   REAL(wp), DIMENSION(3,61), PUBLIC ::   xkrgb   !: tabulated attenuation coefficients for RGB absorption 
    4952    
    50    !!* Substitution 
    51 #  include "top_substitute.h90" 
     53   !! * Substitutions 
     54#  include "domzgr_substitute.h90" 
    5255   !!---------------------------------------------------------------------- 
    5356   !! NEMO/TOP 3.3 , NEMO Consortium (2010) 
     
    5760CONTAINS 
    5861 
    59    SUBROUTINE p4z_opt( kt, jnt ) 
     62   SUBROUTINE p4z_opt( kt, knt ) 
    6063      !!--------------------------------------------------------------------- 
    6164      !!                     ***  ROUTINE p4z_opt  *** 
     
    6770      !!--------------------------------------------------------------------- 
    6871      ! 
    69       INTEGER, INTENT(in) ::   kt, jnt   ! ocean time step 
     72      INTEGER, INTENT(in) ::   kt, knt   ! ocean time step 
    7073      ! 
    7174      INTEGER  ::   ji, jj, jk 
    7275      INTEGER  ::   irgb 
    73       REAL(wp) ::   zchl, zxsi0r 
     76      REAL(wp) ::   zchl 
    7477      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 
     78      REAL(wp), POINTER, DIMENSION(:,:  ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4, zqsr100 
     79      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpar, ze0, ze1, ze2, ze3 
    7780      !!--------------------------------------------------------------------- 
    7881      ! 
     
    8083      ! 
    8184      ! 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 ) 
     85      CALL wrk_alloc( jpi, jpj,      zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 ) 
     86      CALL wrk_alloc( jpi, jpj, jpk, zpar, ze0, ze1, ze2, ze3 ) 
     87 
     88      IF( knt == 1 .AND. ln_varpar ) CALL p4z_opt_sbc( kt ) 
    8689 
    8790      !     Initialisation of variables used to compute PAR 
    8891      !     ----------------------------------------------- 
    89       ze1(:,:,jpk) = 0._wp 
    90       ze2(:,:,jpk) = 0._wp 
    91       ze3(:,:,jpk) = 0._wp 
    92  
     92      ze1(:,:,:) = 0._wp 
     93      ze2(:,:,:) = 0._wp 
     94      ze3(:,:,:) = 0._wp 
    9395      !                                        !* attenuation coef. function of Chlorophyll and wavelength (Red-Green-Blue) 
    9496      DO jk = 1, jpkm1                         !  -------------------------------------------------------- 
    95 !CDIR NOVERRCHK 
    9697         DO jj = 1, jpj 
    97 !CDIR NOVERRCHK 
    9898            DO ji = 1, jpi 
    99                zchl = ( trn(ji,jj,jk,jpnch) + trn(ji,jj,jk,jpdch) + rtrn ) * 1.e6 
     99               zchl = ( trb(ji,jj,jk,jpnch) + trb(ji,jj,jk,jpdch) + rtrn ) * 1.e6 
    100100               zchl = MIN(  10. , MAX( 0.05, zchl )  ) 
    101101               irgb = NINT( 41 + 20.* LOG10( zchl ) + rtrn ) 
    102102               !                                                          
    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) 
     103               ekb(ji,jj,jk) = xkrgb(1,irgb) * fse3t(ji,jj,jk) 
     104               ekg(ji,jj,jk) = xkrgb(2,irgb) * fse3t(ji,jj,jk) 
     105               ekr(ji,jj,jk) = xkrgb(3,irgb) * fse3t(ji,jj,jk) 
    106106            END DO 
    107107         END DO 
    108108      END DO 
    109  
    110  
    111109      !                                        !* Photosynthetically Available Radiation (PAR) 
    112110      !                                        !  -------------------------------------- 
    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) ) 
     111      IF( l_trcdm2dc ) THEN                     !  diurnal cycle 
     112         ! 1% of qsr to compute euphotic layer 
     113         zqsr100(:,:) = 0.01 * qsr_mean(:,:)     !  daily mean qsr 
     114         ! 
     115         CALL p4z_opt_par( kt, qsr_mean, ze1, ze2, ze3 )  
     116         ! 
     117         DO jk = 1, nksrp       
     118            etot_ndcy(:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk) 
     119            enano    (:,:,jk) =  2.1 * ze1(:,:,jk) + 0.42 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk) 
     120            ediat    (:,:,jk) =  1.6 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.7 * ze3(:,:,jk) 
     121         END DO 
     122         ! 
     123         CALL p4z_opt_par( kt, qsr, ze1, ze2, ze3 )  
     124         ! 
     125         DO jk = 1, nksrp       
     126            etot(:,:,jk) =  ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk) 
     127         END DO 
     128         ! 
    118129      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 
     130         ! 1% of qsr to compute euphotic layer 
     131         zqsr100(:,:) = 0.01 * qsr(:,:) 
     132         ! 
     133         CALL p4z_opt_par( kt, qsr, ze1, ze2, ze3 )  
     134         ! 
     135         DO jk = 1, nksrp       
     136            etot (:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk) 
     137            enano(:,:,jk) =  2.1 * ze1(:,:,jk) + 0.42 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk) 
     138            ediat(:,:,jk) =  1.6 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.7 * ze3(:,:,jk) 
     139         END DO 
     140         etot_ndcy(:,:,:) =  etot(:,:,:)  
     141      ENDIF 
     142 
    155143 
    156144      IF( ln_qsr_bio ) THEN                    !* heat flux accros w-level (used in the dynamics) 
    157145         !                                     !  ------------------------ 
    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 
     146         CALL p4z_opt_par( kt, qsr, ze1, ze2, ze3, pe0=ze0 ) 
     147         ! 
    173148         etot3(:,:,1) =  qsr(:,:) * tmask(:,:,1) 
    174          ! 
    175          ! 
    176149         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  
     150            etot3(:,:,jk) =  ( ze0(:,:,jk) + ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk) ) * tmask(:,:,jk) 
     151         END DO 
     152         !                                     !  ------------------------ 
     153      ENDIF 
    198154      !                                        !* Euphotic depth and level 
    199155      neln(:,:) = 1                            !  ------------------------ 
     
    203159         DO jj = 1, jpj 
    204160           DO ji = 1, jpi 
    205               IF( etot(ji,jj,jk) * tmask(ji,jj,jk) >= 0.0043 * qsr(ji,jj) )  THEN 
     161              IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.43 * zqsr100(ji,jj) )  THEN 
    206162                 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 
     163                 !                                     ! nb: ensure the compatibility with nmld_trc definition in trd_mld_trc_zint 
    208164                 heup(ji,jj) = fsdepw(ji,jj,jk+1)      ! Euphotic layer depth 
    209165              ENDIF 
     
    211167        END DO 
    212168      END DO 
    213   
     169      ! 
    214170      heup(:,:) = MIN( 300., heup(:,:) ) 
    215  
    216171      !                                        !* mean light over the mixed layer 
    217172      zdepmoy(:,:)   = 0.e0                    !  ------------------------------- 
    218       zetmp  (:,:)   = 0.e0 
    219173      zetmp1 (:,:)   = 0.e0 
    220174      zetmp2 (:,:)   = 0.e0 
     175      zetmp3 (:,:)   = 0.e0 
     176      zetmp4 (:,:)   = 0.e0 
    221177 
    222178      DO jk = 1, nksrp 
    223 !CDIR NOVERRCHK 
    224179         DO jj = 1, jpj 
    225 !CDIR NOVERRCHK 
    226180            DO ji = 1, jpi 
    227181               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) 
     182                  zetmp1 (ji,jj) = zetmp1 (ji,jj) + etot     (ji,jj,jk) * fse3t(ji,jj,jk) ! remineralisation 
     183                  zetmp2 (ji,jj) = zetmp2 (ji,jj) + etot_ndcy(ji,jj,jk) * fse3t(ji,jj,jk) ! production 
     184                  zetmp3 (ji,jj) = zetmp3 (ji,jj) + enano    (ji,jj,jk) * fse3t(ji,jj,jk) ! production 
     185                  zetmp4 (ji,jj) = zetmp4 (ji,jj) + ediat    (ji,jj,jk) * fse3t(ji,jj,jk) ! production 
    231186                  zdepmoy(ji,jj) = zdepmoy(ji,jj) + fse3t(ji,jj,jk) 
    232187               ENDIF 
     
    235190      END DO 
    236191      ! 
    237       emoy(:,:,:) = etot(:,:,:) 
     192      emoy(:,:,:) = etot(:,:,:)       ! remineralisation 
     193      zpar(:,:,:) = etot_ndcy(:,:,:)  ! diagnostic : PAR with no diurnal cycle  
    238194      ! 
    239195      DO jk = 1, nksrp 
    240 !CDIR NOVERRCHK 
    241196         DO jj = 1, jpj 
    242 !CDIR NOVERRCHK 
    243197            DO ji = 1, jpi 
    244198               IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN 
    245199                  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 
     200                  emoy (ji,jj,jk) = zetmp1(ji,jj) * z1_dep 
     201                  zpar (ji,jj,jk) = zetmp2(ji,jj) * z1_dep 
     202                  enano(ji,jj,jk) = zetmp3(ji,jj) * z1_dep 
     203                  ediat(ji,jj,jk) = zetmp4(ji,jj) * z1_dep 
    249204               ENDIF 
    250205            END DO 
    251206         END DO 
    252207      END DO 
    253  
     208      ! 
    254209      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 
     210        IF( knt == nrdttrc ) THEN 
     211           IF( iom_use( "Heup"  ) ) CALL iom_put( "Heup" , heup(:,:  ) * tmask(:,:,1) )  ! euphotic layer deptht 
     212           IF( iom_use( "PARDM" ) ) CALL iom_put( "PARDM", zpar(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation 
     213           IF( iom_use( "PAR"   ) ) CALL iom_put( "PAR"  , emoy(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation 
    258214        ENDIF 
    259215      ELSE 
    260216         IF( ln_diatrc ) THEN        ! save output diagnostics 
    261             trc2d(:,:,  jp_pcs0_2d + 10) = heup(:,:  ) * tmask(:,:,1)   
     217            trc2d(:,:,  jp_pcs0_2d + 10) = heup(:,:  ) * tmask(:,:,1) 
    262218            trc3d(:,:,:,jp_pcs0_3d + 3)  = etot(:,:,:) * tmask(:,:,:) 
    263219         ENDIF 
    264220      ENDIF 
    265221      ! 
    266       CALL wrk_dealloc( jpi, jpj,      zdepmoy, zetmp, zetmp1, zetmp2 ) 
    267       CALL wrk_dealloc( jpi, jpj, jpk, zekg, zekr, zekb, ze0, ze1, ze2, ze3 ) 
     222      CALL wrk_dealloc( jpi, jpj,      zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 ) 
     223      CALL wrk_dealloc( jpi, jpj, jpk, zpar, ze0, ze1, ze2, ze3 ) 
    268224      ! 
    269225      IF( nn_timing == 1 )  CALL timing_stop('p4z_opt') 
     
    271227   END SUBROUTINE p4z_opt 
    272228 
    273    SUBROUTINE p4z_optsbc( kt ) 
    274       !!---------------------------------------------------------------------- 
    275       !!                  ***  routine p4z_optsbc  *** 
     229   SUBROUTINE p4z_opt_par( kt, pqsr, pe1, pe2, pe3, pe0 )  
     230      !!---------------------------------------------------------------------- 
     231      !!                  ***  routine p4z_opt_par  *** 
     232      !! 
     233      !! ** purpose :   compute PAR of each wavelength (Red-Green-Blue) 
     234      !!                for a given shortwave radiation 
     235      !! 
     236      !!---------------------------------------------------------------------- 
     237      !! * arguments 
     238      INTEGER, INTENT(in)                                       ::  kt            !   ocean time-step 
     239      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in)              ::  pqsr          !   shortwave 
     240      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)           ::  pe1 , pe2 , pe3   !  PAR ( R-G-B) 
     241      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout), OPTIONAL ::  pe0   
     242      !! * local variables 
     243      INTEGER    ::   ji, jj, jk     ! dummy loop indices 
     244      REAL(wp), DIMENSION(jpi,jpj)     ::  zqsr          !   shortwave 
     245      !!---------------------------------------------------------------------- 
     246 
     247      !  Real shortwave 
     248      IF( ln_varpar ) THEN  ;  zqsr(:,:) = par_varsw(:,:) * pqsr(:,:) 
     249      ELSE                  ;  zqsr(:,:) = xparsw         * pqsr(:,:) 
     250      ENDIF 
     251      ! 
     252      IF( PRESENT( pe0 ) ) THEN     !  W-level 
     253         ! 
     254         pe0(:,:,1) = pqsr(:,:) - 3. * zqsr(:,:)    !   ( 1 - 3 * alpha ) * q 
     255         pe1(:,:,1) = zqsr(:,:)          
     256         pe2(:,:,1) = zqsr(:,:) 
     257         pe3(:,:,1) = zqsr(:,:) 
     258         ! 
     259         DO jk = 2, nksrp + 1 
     260            DO jj = 1, jpj 
     261               DO ji = 1, jpi 
     262                  pe0(ji,jj,jk) = pe0(ji,jj,jk-1) * EXP( -fse3t(ji,jj,jk-1) * xsi0r ) 
     263                  pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -ekb(ji,jj,jk-1 ) ) 
     264                  pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -ekg(ji,jj,jk-1 ) ) 
     265                  pe3(ji,jj,jk) = pe3(ji,jj,jk-1) * EXP( -ekr(ji,jj,jk-1 ) ) 
     266               END DO 
     267              ! 
     268            END DO 
     269            ! 
     270         END DO 
     271        ! 
     272      ELSE   ! T- level 
     273        ! 
     274        pe1(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekb(:,:,1) ) 
     275        pe2(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekg(:,:,1) ) 
     276        pe3(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekr(:,:,1) ) 
     277        ! 
     278        DO jk = 2, nksrp       
     279           DO jj = 1, jpj 
     280              DO ji = 1, jpi 
     281                 pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -0.5 * ( ekb(ji,jj,jk-1) + ekb(ji,jj,jk) ) ) 
     282                 pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -0.5 * ( ekg(ji,jj,jk-1) + ekg(ji,jj,jk) ) ) 
     283                 pe3(ji,jj,jk) = pe3(ji,jj,jk-1) * EXP( -0.5 * ( ekr(ji,jj,jk-1) + ekr(ji,jj,jk) ) ) 
     284              END DO 
     285           END DO 
     286        END DO     
     287        ! 
     288      ENDIF 
     289      !  
     290   END SUBROUTINE p4z_opt_par 
     291 
     292 
     293   SUBROUTINE p4z_opt_sbc( kt ) 
     294      !!---------------------------------------------------------------------- 
     295      !!                  ***  routine p4z_opt_sbc  *** 
    276296      !! 
    277297      !! ** purpose :   read and interpolate the variable PAR fraction 
     
    284304      !!---------------------------------------------------------------------- 
    285305      !! * arguments 
    286       INTEGER, INTENT( in  ) ::   kt   ! ocean time step 
     306      INTEGER ,                INTENT(in) ::   kt     ! ocean time step 
    287307 
    288308      !! * local declarations 
     
    297317         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_par > 1 ) ) THEN 
    298318            CALL fld_read( kt, 1, sf_par ) 
    299             par_varsw(:,:) = ( sf_par(1)%fnow(:,:,1) )/3.0 
     319            par_varsw(:,:) = ( sf_par(1)%fnow(:,:,1) ) / 3.0 
    300320         ENDIF 
    301321      ENDIF 
     
    303323      IF( nn_timing == 1 )  CALL timing_stop('p4z_optsbc') 
    304324      ! 
    305    END SUBROUTINE p4z_optsbc 
     325   END SUBROUTINE p4z_opt_sbc 
    306326 
    307327   SUBROUTINE p4z_opt_init 
     
    347367      ! 
    348368      xparsw = parlux / 3.0 
     369      xsi0r  = 1.e0 / rn_si0 
    349370      ! 
    350371      ! Variable PAR at the surface of the ocean 
     
    372393      IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksrp, ' ref depth = ', gdepw_1d(nksrp+1), ' m' 
    373394      ! 
    374                          etot (:,:,:) = 0._wp 
    375                          enano(:,:,:) = 0._wp 
    376                          ediat(:,:,:) = 0._wp 
    377       IF( ln_qsr_bio )   etot3(:,:,:) = 0._wp 
     395                         ekr      (:,:,:) = 0._wp 
     396                         ekb      (:,:,:) = 0._wp 
     397                         ekg      (:,:,:) = 0._wp 
     398                         etot     (:,:,:) = 0._wp 
     399                         etot_ndcy(:,:,:) = 0._wp 
     400                         enano    (:,:,:) = 0._wp 
     401                         ediat    (:,:,:) = 0._wp 
     402      IF( ln_qsr_bio )   etot3    (:,:,:) = 0._wp 
    378403      !  
    379404      IF( nn_timing == 1 )  CALL timing_stop('p4z_opt_init') 
     
    386411      !!                     ***  ROUTINE p4z_opt_alloc  *** 
    387412      !!---------------------------------------------------------------------- 
    388       ALLOCATE( enano(jpi,jpj,jpk), ediat(jpi,jpj,jpk), emoy (jpi,jpj,jpk), STAT=p4z_opt_alloc )  
     413      ALLOCATE( ekb(jpi,jpj,jpk)      , ekr(jpi,jpj,jpk), ekg(jpi,jpj,jpk),   & 
     414        &       enano(jpi,jpj,jpk)    , ediat(jpi,jpj,jpk), & 
     415        &       etot_ndcy(jpi,jpj,jpk), emoy (jpi,jpj,jpk), STAT=p4z_opt_alloc )  
    389416         ! 
    390417      IF( p4z_opt_alloc /= 0 ) CALL ctl_warn('p4z_opt_alloc : failed to allocate arrays.') 
     
    402429 
    403430   !!====================================================================== 
    404 END MODULE  p4zopt 
     431END MODULE p4zopt 
Note: See TracChangeset for help on using the changeset viewer.