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

Ignore:
Timestamp:
2016-01-08T10:35:19+01:00 (8 years ago)
Author:
jamesharle
Message:

Update MPP_BDY_UPDATE branch to be consistent with head of trunk

File:
1 edited

Legend:

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

    r4624 r6225  
    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" 
    5253   !!---------------------------------------------------------------------- 
    5354   !! NEMO/TOP 3.3 , NEMO Consortium (2010) 
     
    5758CONTAINS 
    5859 
    59    SUBROUTINE p4z_opt( kt, jnt ) 
     60   SUBROUTINE p4z_opt( kt, knt ) 
    6061      !!--------------------------------------------------------------------- 
    6162      !!                     ***  ROUTINE p4z_opt  *** 
     
    6768      !!--------------------------------------------------------------------- 
    6869      ! 
    69       INTEGER, INTENT(in) ::   kt, jnt   ! ocean time step 
     70      INTEGER, INTENT(in) ::   kt, knt   ! ocean time step 
    7071      ! 
    7172      INTEGER  ::   ji, jj, jk 
    7273      INTEGER  ::   irgb 
    73       REAL(wp) ::   zchl, zxsi0r 
     74      REAL(wp) ::   zchl 
    7475      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 
     76      REAL(wp), POINTER, DIMENSION(:,:  ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4, zqsr100 
     77      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpar, ze0, ze1, ze2, ze3 
    7778      !!--------------------------------------------------------------------- 
    7879      ! 
     
    8081      ! 
    8182      ! 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 ) 
     83      CALL wrk_alloc( jpi, jpj,      zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 ) 
     84      CALL wrk_alloc( jpi, jpj, jpk, zpar, ze0, ze1, ze2, ze3 ) 
     85 
     86      IF( knt == 1 .AND. ln_varpar ) CALL p4z_opt_sbc( kt ) 
    8687 
    8788      !     Initialisation of variables used to compute PAR 
    8889      !     ----------------------------------------------- 
    89       ze1(:,:,jpk) = 0._wp 
    90       ze2(:,:,jpk) = 0._wp 
    91       ze3(:,:,jpk) = 0._wp 
    92  
     90      ze1(:,:,:) = 0._wp 
     91      ze2(:,:,:) = 0._wp 
     92      ze3(:,:,:) = 0._wp 
    9393      !                                        !* attenuation coef. function of Chlorophyll and wavelength (Red-Green-Blue) 
    9494      DO jk = 1, jpkm1                         !  -------------------------------------------------------- 
    95 !CDIR NOVERRCHK 
    9695         DO jj = 1, jpj 
    97 !CDIR NOVERRCHK 
    9896            DO ji = 1, jpi 
    99                zchl = ( trn(ji,jj,jk,jpnch) + trn(ji,jj,jk,jpdch) + rtrn ) * 1.e6 
     97               zchl = ( trb(ji,jj,jk,jpnch) + trb(ji,jj,jk,jpdch) + rtrn ) * 1.e6 
    10098               zchl = MIN(  10. , MAX( 0.05, zchl )  ) 
    10199               irgb = NINT( 41 + 20.* LOG10( zchl ) + rtrn ) 
    102100               !                                                          
    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) 
     101               ekb(ji,jj,jk) = xkrgb(1,irgb) * e3t_n(ji,jj,jk) 
     102               ekg(ji,jj,jk) = xkrgb(2,irgb) * e3t_n(ji,jj,jk) 
     103               ekr(ji,jj,jk) = xkrgb(3,irgb) * e3t_n(ji,jj,jk) 
    106104            END DO 
    107105         END DO 
    108106      END DO 
    109  
    110  
    111107      !                                        !* Photosynthetically Available Radiation (PAR) 
    112108      !                                        !  -------------------------------------- 
    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) ) 
     109      IF( l_trcdm2dc ) THEN                     !  diurnal cycle 
     110         ! 1% of qsr to compute euphotic layer 
     111         zqsr100(:,:) = 0.01 * qsr_mean(:,:)     !  daily mean qsr 
     112         ! 
     113         CALL p4z_opt_par( kt, qsr_mean, ze1, ze2, ze3 )  
     114         ! 
     115         DO jk = 1, nksrp       
     116            etot_ndcy(:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk) 
     117            enano    (:,:,jk) =  2.1 * ze1(:,:,jk) + 0.42 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk) 
     118            ediat    (:,:,jk) =  1.6 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.7 * ze3(:,:,jk) 
     119         END DO 
     120         ! 
     121         CALL p4z_opt_par( kt, qsr, ze1, ze2, ze3 )  
     122         ! 
     123         DO jk = 1, nksrp       
     124            etot(:,:,jk) =  ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk) 
     125         END DO 
     126         ! 
    118127      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 
     128         ! 1% of qsr to compute euphotic layer 
     129         zqsr100(:,:) = 0.01 * qsr(:,:) 
     130         ! 
     131         CALL p4z_opt_par( kt, qsr, ze1, ze2, ze3 )  
     132         ! 
     133         DO jk = 1, nksrp       
     134            etot (:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk) 
     135            enano(:,:,jk) =  2.1 * ze1(:,:,jk) + 0.42 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk) 
     136            ediat(:,:,jk) =  1.6 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.7 * ze3(:,:,jk) 
     137         END DO 
     138         etot_ndcy(:,:,:) =  etot(:,:,:)  
     139      ENDIF 
     140 
    155141 
    156142      IF( ln_qsr_bio ) THEN                    !* heat flux accros w-level (used in the dynamics) 
    157143         !                                     !  ------------------------ 
    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 
     144         CALL p4z_opt_par( kt, qsr, ze1, ze2, ze3, pe0=ze0 ) 
     145         ! 
    173146         etot3(:,:,1) =  qsr(:,:) * tmask(:,:,1) 
    174          ! 
    175          ! 
    176147         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  
     148            etot3(:,:,jk) =  ( ze0(:,:,jk) + ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk) ) * tmask(:,:,jk) 
     149         END DO 
     150         !                                     !  ------------------------ 
     151      ENDIF 
    198152      !                                        !* Euphotic depth and level 
    199153      neln(:,:) = 1                            !  ------------------------ 
     
    203157         DO jj = 1, jpj 
    204158           DO ji = 1, jpi 
    205               IF( etot(ji,jj,jk) * tmask(ji,jj,jk) >= 0.0043 * qsr(ji,jj) )  THEN 
     159              IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.43 * zqsr100(ji,jj) )  THEN 
    206160                 neln(ji,jj) = jk+1                    ! Euphotic level : 1rst T-level strictly below Euphotic layer 
    207161                 !                                     ! 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 
     162                 heup(ji,jj) = gdepw_n(ji,jj,jk+1)     ! Euphotic layer depth 
    209163              ENDIF 
    210164           END DO 
    211165        END DO 
    212166      END DO 
    213   
     167      ! 
    214168      heup(:,:) = MIN( 300., heup(:,:) ) 
    215  
    216169      !                                        !* mean light over the mixed layer 
    217170      zdepmoy(:,:)   = 0.e0                    !  ------------------------------- 
    218       zetmp  (:,:)   = 0.e0 
    219171      zetmp1 (:,:)   = 0.e0 
    220172      zetmp2 (:,:)   = 0.e0 
     173      zetmp3 (:,:)   = 0.e0 
     174      zetmp4 (:,:)   = 0.e0 
    221175 
    222176      DO jk = 1, nksrp 
    223 !CDIR NOVERRCHK 
    224177         DO jj = 1, jpj 
    225 !CDIR NOVERRCHK 
    226178            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) 
     179               IF( gdepw_n(ji,jj,jk+1) <= hmld(ji,jj) ) THEN 
     180                  zetmp1 (ji,jj) = zetmp1 (ji,jj) + etot     (ji,jj,jk) * e3t_n(ji,jj,jk) ! remineralisation 
     181                  zetmp2 (ji,jj) = zetmp2 (ji,jj) + etot_ndcy(ji,jj,jk) * e3t_n(ji,jj,jk) ! production 
     182                  zetmp3 (ji,jj) = zetmp3 (ji,jj) + enano    (ji,jj,jk) * e3t_n(ji,jj,jk) ! production 
     183                  zetmp4 (ji,jj) = zetmp4 (ji,jj) + ediat    (ji,jj,jk) * e3t_n(ji,jj,jk) ! production 
     184                  zdepmoy(ji,jj) = zdepmoy(ji,jj) +                       e3t_n(ji,jj,jk) 
    232185               ENDIF 
    233186            END DO 
     
    235188      END DO 
    236189      ! 
    237       emoy(:,:,:) = etot(:,:,:) 
     190      emoy(:,:,:) = etot(:,:,:)       ! remineralisation 
     191      zpar(:,:,:) = etot_ndcy(:,:,:)  ! diagnostic : PAR with no diurnal cycle  
    238192      ! 
    239193      DO jk = 1, nksrp 
    240 !CDIR NOVERRCHK 
    241194         DO jj = 1, jpj 
    242 !CDIR NOVERRCHK 
    243195            DO ji = 1, jpi 
    244                IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN 
     196               IF( gdepw_n(ji,jj,jk+1) <= hmld(ji,jj) ) THEN 
    245197                  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 
     198                  emoy (ji,jj,jk) = zetmp1(ji,jj) * z1_dep 
     199                  zpar (ji,jj,jk) = zetmp2(ji,jj) * z1_dep 
     200                  enano(ji,jj,jk) = zetmp3(ji,jj) * z1_dep 
     201                  ediat(ji,jj,jk) = zetmp4(ji,jj) * z1_dep 
    249202               ENDIF 
    250203            END DO 
    251204         END DO 
    252205      END DO 
    253  
    254       IF( ln_diatrc ) THEN        ! save output diagnostics 
     206      ! 
     207      IF( lk_iomput ) THEN 
     208        IF( knt == nrdttrc ) THEN 
     209           IF( iom_use( "Heup"  ) ) CALL iom_put( "Heup" , heup(:,:  ) * tmask(:,:,1) )  ! euphotic layer deptht 
     210           IF( iom_use( "PARDM" ) ) CALL iom_put( "PARDM", zpar(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation 
     211           IF( iom_use( "PAR"   ) ) CALL iom_put( "PAR"  , emoy(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation 
     212        ENDIF 
     213      ELSE 
     214         IF( ln_diatrc ) THEN        ! save output diagnostics 
     215            trc2d(:,:,  jp_pcs0_2d + 10) = heup(:,:  ) * tmask(:,:,1) 
     216            trc3d(:,:,:,jp_pcs0_3d + 3)  = etot(:,:,:) * tmask(:,:,:) 
     217         ENDIF 
     218      ENDIF 
     219      ! 
     220      CALL wrk_dealloc( jpi, jpj,      zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 ) 
     221      CALL wrk_dealloc( jpi, jpj, jpk, zpar,  ze0, ze1, ze2, ze3 ) 
     222      ! 
     223      IF( nn_timing == 1 )  CALL timing_stop('p4z_opt') 
     224      ! 
     225   END SUBROUTINE p4z_opt 
     226 
     227   SUBROUTINE p4z_opt_par( kt, pqsr, pe1, pe2, pe3, pe0 )  
     228      !!---------------------------------------------------------------------- 
     229      !!                  ***  routine p4z_opt_par  *** 
     230      !! 
     231      !! ** purpose :   compute PAR of each wavelength (Red-Green-Blue) 
     232      !!                for a given shortwave radiation 
     233      !! 
     234      !!---------------------------------------------------------------------- 
     235      !! * arguments 
     236      INTEGER, INTENT(in)                                       ::  kt            !   ocean time-step 
     237      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in)              ::  pqsr          !   shortwave 
     238      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)           ::  pe1 , pe2 , pe3   !  PAR ( R-G-B) 
     239      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout), OPTIONAL ::  pe0   
     240      !! * local variables 
     241      INTEGER    ::   ji, jj, jk     ! dummy loop indices 
     242      REAL(wp), DIMENSION(jpi,jpj)     ::  zqsr          !   shortwave 
     243      !!---------------------------------------------------------------------- 
     244 
     245      !  Real shortwave 
     246      IF( ln_varpar ) THEN  ;  zqsr(:,:) = par_varsw(:,:) * pqsr(:,:) 
     247      ELSE                  ;  zqsr(:,:) = xparsw         * pqsr(:,:) 
     248      ENDIF 
     249      ! 
     250      IF( PRESENT( pe0 ) ) THEN     !  W-level 
     251         ! 
     252         pe0(:,:,1) = pqsr(:,:) - 3. * zqsr(:,:)    !   ( 1 - 3 * alpha ) * q 
     253         pe1(:,:,1) = zqsr(:,:)          
     254         pe2(:,:,1) = zqsr(:,:) 
     255         pe3(:,:,1) = zqsr(:,:) 
     256         ! 
     257         DO jk = 2, nksrp + 1 
     258            DO jj = 1, jpj 
     259               DO ji = 1, jpi 
     260                  pe0(ji,jj,jk) = pe0(ji,jj,jk-1) * EXP( -e3t_n(ji,jj,jk-1) * xsi0r ) 
     261                  pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -ekb(ji,jj,jk-1 ) ) 
     262                  pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -ekg(ji,jj,jk-1 ) ) 
     263                  pe3(ji,jj,jk) = pe3(ji,jj,jk-1) * EXP( -ekr(ji,jj,jk-1 ) ) 
     264               END DO 
     265              ! 
     266            END DO 
     267            ! 
     268         END DO 
    255269        ! 
    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 
     270      ELSE   ! T- level 
    265271        ! 
    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  *** 
     272        pe1(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekb(:,:,1) ) 
     273        pe2(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekg(:,:,1) ) 
     274        pe3(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekr(:,:,1) ) 
     275        ! 
     276        DO jk = 2, nksrp       
     277           DO jj = 1, jpj 
     278              DO ji = 1, jpi 
     279                 pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -0.5 * ( ekb(ji,jj,jk-1) + ekb(ji,jj,jk) ) ) 
     280                 pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -0.5 * ( ekg(ji,jj,jk-1) + ekg(ji,jj,jk) ) ) 
     281                 pe3(ji,jj,jk) = pe3(ji,jj,jk-1) * EXP( -0.5 * ( ekr(ji,jj,jk-1) + ekr(ji,jj,jk) ) ) 
     282              END DO 
     283           END DO 
     284        END DO     
     285        ! 
     286      ENDIF 
     287      !  
     288   END SUBROUTINE p4z_opt_par 
     289 
     290 
     291   SUBROUTINE p4z_opt_sbc( kt ) 
     292      !!---------------------------------------------------------------------- 
     293      !!                  ***  routine p4z_opt_sbc  *** 
    278294      !! 
    279295      !! ** purpose :   read and interpolate the variable PAR fraction 
     
    286302      !!---------------------------------------------------------------------- 
    287303      !! * arguments 
    288       INTEGER, INTENT( in  ) ::   kt   ! ocean time step 
     304      INTEGER ,                INTENT(in) ::   kt     ! ocean time step 
    289305 
    290306      !! * local declarations 
     
    299315         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_par > 1 ) ) THEN 
    300316            CALL fld_read( kt, 1, sf_par ) 
    301             par_varsw(:,:) = ( sf_par(1)%fnow(:,:,1) )/3.0 
     317            par_varsw(:,:) = ( sf_par(1)%fnow(:,:,1) ) / 3.0 
    302318         ENDIF 
    303319      ENDIF 
     
    305321      IF( nn_timing == 1 )  CALL timing_stop('p4z_optsbc') 
    306322      ! 
    307    END SUBROUTINE p4z_optsbc 
     323   END SUBROUTINE p4z_opt_sbc 
    308324 
    309325   SUBROUTINE p4z_opt_init 
     
    349365      ! 
    350366      xparsw = parlux / 3.0 
     367      xsi0r  = 1.e0 / rn_si0 
    351368      ! 
    352369      ! Variable PAR at the surface of the ocean 
     
    374391      IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksrp, ' ref depth = ', gdepw_1d(nksrp+1), ' m' 
    375392      ! 
    376                          etot (:,:,:) = 0._wp 
    377                          enano(:,:,:) = 0._wp 
    378                          ediat(:,:,:) = 0._wp 
    379       IF( ln_qsr_bio )   etot3(:,:,:) = 0._wp 
     393                         ekr      (:,:,:) = 0._wp 
     394                         ekb      (:,:,:) = 0._wp 
     395                         ekg      (:,:,:) = 0._wp 
     396                         etot     (:,:,:) = 0._wp 
     397                         etot_ndcy(:,:,:) = 0._wp 
     398                         enano    (:,:,:) = 0._wp 
     399                         ediat    (:,:,:) = 0._wp 
     400      IF( ln_qsr_bio )   etot3    (:,:,:) = 0._wp 
    380401      !  
    381402      IF( nn_timing == 1 )  CALL timing_stop('p4z_opt_init') 
     
    388409      !!                     ***  ROUTINE p4z_opt_alloc  *** 
    389410      !!---------------------------------------------------------------------- 
    390       ALLOCATE( enano(jpi,jpj,jpk), ediat(jpi,jpj,jpk), emoy (jpi,jpj,jpk), STAT=p4z_opt_alloc )  
     411      ALLOCATE( ekb(jpi,jpj,jpk)      , ekr(jpi,jpj,jpk), ekg(jpi,jpj,jpk),   & 
     412        &       enano(jpi,jpj,jpk)    , ediat(jpi,jpj,jpk), & 
     413        &       etot_ndcy(jpi,jpj,jpk), emoy (jpi,jpj,jpk), STAT=p4z_opt_alloc )  
    391414         ! 
    392415      IF( p4z_opt_alloc /= 0 ) CALL ctl_warn('p4z_opt_alloc : failed to allocate arrays.') 
     
    404427 
    405428   !!====================================================================== 
    406 END MODULE  p4zopt 
     429END MODULE p4zopt 
Note: See TracChangeset for help on using the changeset viewer.