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 2823 for branches/2011/dev_r2787_PISCES_improvment/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zprod.F90 – NEMO

Ignore:
Timestamp:
2011-08-09T13:11:24+02:00 (13 years ago)
Author:
cetlod
Message:

Add new parameterisation in PISCES, see ticket #854

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2011/dev_r2787_PISCES_improvment/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zprod.F90

    r2730 r2823  
    66   !! History :   1.0  !  2004     (O. Aumont) Original code 
    77   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90 
     8   !!             3.4  !  2011-05  (O. Aumont, C. Ethe) New parameterization of light limitation 
    89   !!---------------------------------------------------------------------- 
    910#if defined key_pisces 
     
    1112   !!   'key_pisces'                                       PISCES bio-model 
    1213   !!---------------------------------------------------------------------- 
    13    !!   p4z_prod       :   
    14    !!---------------------------------------------------------------------- 
    15    USE trc 
    16    USE oce_trc         ! 
    17    USE sms_pisces      !  
    18    USE prtctl_trc 
    19    USE p4zopt 
    20    USE p4zint 
    21    USE p4zlim 
    22    USE iom 
     14   !!---------------------------------------------------------------------- 
     15   !!   p4z_prod       :   Compute the growth Rate of the two phytoplanktons groups 
     16   !!   p4z_prod_init  :   Initialization of the parameters for growth 
     17   !!   p4z_prod_alloc :   Allocate variables for growth 
     18   !!---------------------------------------------------------------------- 
     19   USE oce_trc         !  shared variables between ocean and passive tracers 
     20   USE trc             !  passive tracers common variables  
     21   USE sms_pisces      !  PISCES Source Minus Sink variables 
     22   USE p4zopt          !  optical model 
     23   USE p4zlim          !  Co-limitations of differents nutrients 
     24   USE prtctl_trc      !  print control for debugging 
     25   USE iom             !  I/O manager 
    2326 
    2427   IMPLICIT NONE 
     
    2932   PUBLIC   p4z_prod_alloc 
    3033 
    31    REAL(wp), PUBLIC ::   & 
    32      pislope   = 3.0_wp          ,  &  !: 
    33      pislope2  = 3.0_wp          ,  &  !: 
    34      excret    = 10.e-5_wp       , &   !: 
    35      excret2   = 0.05_wp         , &   !: 
    36      chlcnm    = 0.033_wp        , &   !: 
    37      chlcdm    = 0.05_wp         , &   !: 
    38      fecnm     = 10.E-6_wp       , &   !: 
    39      fecdm     = 15.E-6_wp       , &   !: 
    40      grosip    = 0.151_wp 
    41  
    42    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prmax   !: 
     34   !! * Shared module variables 
     35   REAL(wp), PUBLIC ::  pislope   = 3.0_wp            !: 
     36   REAL(wp), PUBLIC ::  pislope2  = 3.0_wp            !: 
     37   REAL(wp), PUBLIC ::  excret    = 10.e-5_wp         !: 
     38   REAL(wp), PUBLIC ::  excret2   = 0.05_wp           !: 
     39   REAL(wp), PUBLIC ::  bresp     = 0.00333_wp        !: 
     40   REAL(wp), PUBLIC ::  chlcnm    = 0.033_wp          !: 
     41   REAL(wp), PUBLIC ::  chlcdm    = 0.05_wp           !: 
     42   REAL(wp), PUBLIC ::  chlcmin   = 0.00333_wp        !: 
     43   REAL(wp), PUBLIC ::  fecnm     = 10.E-6_wp         !: 
     44   REAL(wp), PUBLIC ::  fecdm     = 15.E-6_wp         !: 
     45   REAL(wp), PUBLIC ::  grosip    = 0.151_wp          !: 
     46 
     47   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prmax     !: optimal prduction = f(temperature) 
     48   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   quotan   !: proxy of N quota in Nanophyto 
     49   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   quotad   !: proxy of N quota in diatomee 
    4350    
    44    REAL(wp) ::   & 
    45       rday1                      ,  &  !: 0.6 / rday 
    46       texcret                    ,  &  !: 1 - excret  
    47       texcret2                   ,  &  !: 1 - excret2         
    48       tpp                              !: Total primary production 
     51   REAL(wp) :: r1_rday                !: 1 / rday 
     52   REAL(wp) :: r1_bresp               !: 1 / bresp  
     53   REAL(wp) :: texcret                !: 1 - excret  
     54   REAL(wp) :: texcret2               !: 1 - excret2         
     55   REAL(wp) :: tpp                    !: Total primary production 
     56 
    4957 
    5058   !!* Substitution 
     
    6775      !!--------------------------------------------------------------------- 
    6876      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    69       USE wrk_nemo, ONLY:   zmixnano   => wrk_2d_1  , zmixdiat    => wrk_2d_2  , zstrn  => wrk_2d_3 
    70       USE wrk_nemo, ONLY:   zpislopead => wrk_3d_2  , zpislopead2 => wrk_3d_3 
    71       USE wrk_nemo, ONLY:   zprdia     => wrk_3d_4  , zprbio      => wrk_3d_5  , zysopt => wrk_3d_6 
    72       USE wrk_nemo, ONLY:   zprorca    => wrk_3d_7  , zprorcad    => wrk_3d_8 
    73       USE wrk_nemo, ONLY:   zprofed    => wrk_3d_9  , zprofen     => wrk_3d_10 
    74       USE wrk_nemo, ONLY:   zprochln   => wrk_3d_11 , zprochld    => wrk_3d_12 
    75       USE wrk_nemo, ONLY:   zpronew    => wrk_3d_13 , zpronewd    => wrk_3d_14 
     77      USE wrk_nemo, ONLY:   zmixnano   => wrk_2d_1 , zmixdiat    => wrk_2d_2, zstrn => wrk_2d_3 
     78      USE wrk_nemo, ONLY:   zpislopead => wrk_3d_2 , zpislopead2 => wrk_3d_3 
     79      USE wrk_nemo, ONLY:   zprdia     => wrk_3d_4 , zprbio      => wrk_3d_5  
     80      USE wrk_nemo, ONLY:   zprdch     => wrk_3d_6 , zprnch      => wrk_3d_7 
     81      USE wrk_nemo, ONLY:   zprorca    => wrk_3d_8 , zprorcad    => wrk_3d_9 
     82      USE wrk_nemo, ONLY:   zprofed    => wrk_3d_10, zprofen     => wrk_3d_11 
     83      USE wrk_nemo, ONLY:   zprochln   => wrk_3d_12, zprochld    => wrk_3d_13 
     84      USE wrk_nemo, ONLY:   zpronew    => wrk_3d_14, zpronewd    => wrk_3d_15 
    7685      ! 
    7786      INTEGER, INTENT(in) :: kt, jnt 
    7887      ! 
    7988      INTEGER  ::   ji, jj, jk 
    80       REAL(wp) ::   zsilfac, zfact 
    81       REAL(wp) ::   zprdiachl, zprbiochl, zsilim, ztn, zadap, zadap2 
    82       REAL(wp) ::   zlim, zsilfac2, zsiborn, zprod, zetot2, zmax, zproreg, zproreg2 
    83       REAL(wp) ::   zmxltst, zmxlday, zlim1 
     89      REAL(wp) ::   zsilfac, zfact, znanotot, zdiattot, zconctemp, zconctemp2 
     90      REAL(wp) ::   zratio, zmax, zsilim, ztn, zadap 
     91      REAL(wp) ::   zlim, zsilfac2, zsiborn, zprod, zetot2, zproreg, zproreg2 
     92      REAL(wp) ::   zmxltst, zmxlday 
    8493      REAL(wp) ::   zpislopen  , zpislope2n 
    85       REAL(wp) ::   zrum, zcodel, zargu, zval, zvol 
    86 #if defined key_diatrc 
     94      REAL(wp) ::   zrum, zcodel, zargu, zval 
     95      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zysopt  
    8796      REAL(wp) ::   zrfact2 
    88 #endif 
    8997      CHARACTER (len=25) :: charout 
    9098      !!--------------------------------------------------------------------- 
    9199 
    92100      IF( wrk_in_use(2, 1,2,3)                             .OR.  & 
    93           wrk_in_use(3, 2,3,4,5,6,7,8,9,10,11,12,13,14)  ) THEN 
     101          wrk_in_use(3, 2,3,4,5,6,7,8,9,10,11,12,13,14,15)  ) THEN 
    94102          CALL ctl_stop('p4z_prod: requested workspace arrays unavailable')   ;   RETURN 
    95103      ENDIF 
     104    
     105      ALLOCATE( zysopt(jpi,jpj,jpk) ) 
    96106 
    97107      zprorca (:,:,:) = 0._wp 
     
    105115      zprdia  (:,:,:) = 0._wp 
    106116      zprbio  (:,:,:) = 0._wp 
     117      zprdch  (:,:,:) = 0._wp 
     118      zprnch  (:,:,:) = 0._wp 
    107119      zysopt  (:,:,:) = 0._wp 
    108120 
    109121      ! Computation of the optimal production 
    110 # if defined key_degrad 
    111       prmax(:,:,:) = rday1 * tgfunc(:,:,:) * facvol(:,:,:) 
    112 # else 
    113       prmax(:,:,:) = rday1 * tgfunc(:,:,:) 
    114 # endif 
     122      prmax(:,:,:) = 0.6_wp * r1_rday * tgfunc(:,:,:)  
     123      IF( lk_degrad ) THEN 
     124        prmax(:,:,:) = prmax(:,:,:) * facvol(:,:,:) 
     125      ENDIF 
    115126 
    116127      ! compute the day length depending on latitude and the day 
     
    119130 
    120131      ! day length in hours 
    121       zstrn(:,:) = 0._wp 
     132      zstrn(:,:) = 0. 
    122133      DO jj = 1, jpj 
    123134         DO ji = 1, jpi 
    124135            zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad ) 
    125136            zargu = MAX( -1., MIN(  1., zargu ) ) 
    126             zval  = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. ) 
    127             IF( zval < 1.e0 )   zval = 24. 
    128             zstrn(ji,jj) = 24. / zval 
    129          END DO 
    130       END DO 
    131  
     137            zstrn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. ) 
     138         END DO 
     139      END DO 
     140 
     141     ! Impact of the day duration on phytoplankton growth 
     142      DO jk = 1, jpkm1 
     143         DO jj = 1 ,jpj 
     144            DO ji = 1, jpi 
     145               zval = MAX( 1., zstrn(ji,jj) ) 
     146               zval = 1.5 * zval / ( 12. + zval ) 
     147               zprbio(ji,jj,jk) = prmax(ji,jj,jk) * zval 
     148               zprdia(ji,jj,jk) = zprbio(ji,jj,jk) 
     149            END DO 
     150         END DO 
     151      END DO 
     152 
     153      WHERE( zstrn(:,:) < 1.e0 ) zstrn(:,:) = 24. 
     154      zstrn(:,:) = 24. / zstrn(:,:) 
    132155 
    133156!CDIR NOVERRCHK 
     
    141164               IF( etot(ji,jj,jk) > 1.E-3 ) THEN 
    142165                   ztn    = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. ) 
    143                    zadap  = 0.+ 1.* ztn / ( 2.+ ztn ) 
    144                    zadap2 = 0.e0 
    145  
    146                    zfact  = EXP( -0.21 * emoy(ji,jj,jk) ) 
    147  
    148                    zpislopead (ji,jj,jk) = pislope  * ( 1.+ zadap  * zfact ) 
    149                    zpislopead2(ji,jj,jk) = pislope2 * ( 1.+ zadap2 * zfact ) 
    150  
    151                    zpislopen = zpislopead(ji,jj,jk) * trn(ji,jj,jk,jpnch)                 & 
    152                      &         / ( trn(ji,jj,jk,jpphy) * 12.                   + rtrn )   & 
    153                      &         / ( prmax(ji,jj,jk) * rday * xlimphy(ji,jj,jk) + rtrn ) 
    154  
    155                    zpislope2n = zpislopead2(ji,jj,jk) * trn(ji,jj,jk,jpdch)                & 
    156                      &          / ( trn(ji,jj,jk,jpdia) * 12.                   + rtrn )   & 
    157                      &          / ( prmax(ji,jj,jk) * rday * xlimdia(ji,jj,jk) + rtrn ) 
    158  
    159                    ! Computation of production function 
    160                    zprbio(ji,jj,jk) = prmax(ji,jj,jk) * & 
    161                      &                (  1.- EXP( -zpislopen * enano(ji,jj,jk) )  ) 
    162                    zprdia(ji,jj,jk) = prmax(ji,jj,jk) * & 
    163                      &                (  1.- EXP( -zpislope2n * ediat(ji,jj,jk) )  ) 
     166                   zadap  = ztn / ( 2.+ ztn ) 
     167 
     168                   zconctemp   = MAX( 0.e0 , trn(ji,jj,jk,jpdia) - 5e-7 ) 
     169                   zconctemp2  = trn(ji,jj,jk,jpdia) - zconctemp 
     170 
     171                   znanotot = enano(ji,jj,jk) * zstrn(ji,jj) 
     172                   zdiattot = ediat(ji,jj,jk) * zstrn(ji,jj) 
     173                   zfact  = EXP( -0.21 * znanotot ) 
     174 
     175                   zpislopead (ji,jj,jk) = pislope  * ( 1.+ zadap  * zfact )  & 
     176                   &                       * trn(ji,jj,jk,jpnch) /( trn(ji,jj,jk,jpphy) * 12. + rtrn) 
     177 
     178                   zpislopead2(ji,jj,jk) = (pislope * zconctemp2 + pislope2 * zconctemp)   & 
     179                     &                     / ( trn(ji,jj,jk,jpdia) + rtrn )   & 
     180                     &                     * trn(ji,jj,jk,jpdch) /( trn(ji,jj,jk,jpdia) * 12. + rtrn) 
     181 
     182                  !  Computation of production function for Chlorophyll 
     183                  !-------------------------------------------------- 
     184                   zpislopen  = zpislopead (ji,jj,jk) / ( prmax(ji,jj,jk) * rday + rtrn ) 
     185                   zpislope2n = zpislopead2(ji,jj,jk) / ( prmax(ji,jj,jk) * rday + rtrn ) 
     186                   zprnch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopen * znanotot ) ) 
     187                   zprdch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislope2n * zdiattot ) ) 
     188 
     189                   ! Computation of production function for Carbon 
     190                   !  --------------------------------------------- 
     191                   zpislopen  = zpislopead (ji,jj,jk) / ( ( r1_rday + r1_bresp / chlcnm ) * rday + rtrn) 
     192                   zpislope2n = zpislopead2(ji,jj,jk) / ( ( r1_rday + r1_bresp / chlcdm ) * rday + rtrn) 
     193                   zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen  * znanotot )  ) 
     194                   zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpislope2n * zdiattot )  ) 
    164195               ENDIF 
     196            END DO 
     197         END DO 
     198      END DO 
     199 
     200      !  Computation of a proxy of the N/C ratio 
     201      !  --------------------------------------- 
     202!CDIR NOVERRCHK 
     203      DO jk = 1, jpkm1 
     204!CDIR NOVERRCHK 
     205         DO jj = 1, jpj 
     206!CDIR NOVERRCHK 
     207            DO ji = 1, jpi 
     208                zval = ( xnanonh4(ji,jj,jk) + xnanono3(ji,jj,jk) ) * prmax(ji,jj,jk) / ( zprbio(ji,jj,jk) + rtrn ) 
     209                quotan(ji,jj,jk) = MIN( 1., 0.5 + 0.5 * zval ) 
     210                zval = ( xdiatnh4(ji,jj,jk) + xdiatno3(ji,jj,jk) ) * prmax(ji,jj,jk) / ( zprdia(ji,jj,jk) + rtrn ) 
     211                quotad(ji,jj,jk) = MIN( 1., 0.5 + 0.5 * zval ) 
    165212            END DO 
    166213         END DO 
     
    178225                   !    Si/C is arbitrariliy increased for very high Si concentrations 
    179226                   !    to mimic the very high ratios observed in the Southern Ocean (silpot2) 
    180  
    181                   zlim1  = trn(ji,jj,jk,jpsil) / ( trn(ji,jj,jk,jpsil) + xksi1 ) 
    182                   zlim   = xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) 
    183  
    184                   zsilim = MIN( zprdia(ji,jj,jk)    / ( rtrn + prmax(ji,jj,jk) ),                 & 
    185                   &          trn(ji,jj,jk,jpfer) / ( concdfe(ji,jj,jk) + trn(ji,jj,jk,jpfer) ),   & 
    186                   &          trn(ji,jj,jk,jppo4) / ( concdnh4 + trn(ji,jj,jk,jppo4) ),            & 
    187                   &          zlim ) 
    188                   zsilfac = 5.4 * EXP( -4.23 * zsilim ) * MAX( 0.e0, MIN( 1., 2.2 * ( zlim1 - 0.5 ) )  ) + 1.e0 
     227                  zlim  = trn(ji,jj,jk,jpsil) / ( trn(ji,jj,jk,jpsil) + xksi1 ) 
     228                  zsilim = MIN( zprdia(ji,jj,jk) / ( prmax(ji,jj,jk) + rtrn ), xlimsi(ji,jj,jk) ) 
     229                  zsilfac = 4.4 * EXP( -4.23 * zsilim ) * MAX( 0.e0, MIN( 1., 2.2 * ( zlim - 0.5 ) )  ) + 1.e0 
    189230                  zsiborn = MAX( 0.e0, ( trn(ji,jj,jk,jpsil) - 15.e-6 ) ) 
    190                   zsilfac2 = 1.+ 3.* zsiborn / ( zsiborn + xksi2 ) 
    191                   zsilfac = MIN( 6.4,zsilfac * zsilfac2) 
    192                   zysopt(ji,jj,jk) = grosip * zlim1 * zsilfac 
     231                  zsilfac2 = 1.+ 2.* zsiborn / ( zsiborn + xksi2 ) 
     232                  zsilfac = MIN( 5.4, zsilfac * zsilfac2) 
     233                  zysopt(ji,jj,jk) = grosip * zlim * zsilfac 
    193234              ENDIF 
    194235            END DO 
     
    201242         DO ji = 1, jpi 
    202243            zmxltst = MAX( 0.e0, hmld(ji,jj) - heup(ji,jj) ) 
    203             zmxlday = zmxltst**2 / rday 
    204             zmixnano(ji,jj) = 1.- zmxlday / ( 1.+ zmxlday ) 
    205             zmixdiat(ji,jj) = 1.- zmxlday / ( 3.+ zmxlday ) 
     244            zmxlday = zmxltst * zmxltst * r1_rday 
     245            zmixnano(ji,jj) = 1. - zmxlday / ( 3. + zmxlday ) 
     246            zmixdiat(ji,jj) = 1. - zmxlday / ( 4. + zmxlday ) 
    206247         END DO 
    207248      END DO 
     
    214255                  zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * zmixnano(ji,jj) 
    215256                  zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * zmixdiat(ji,jj) 
     257                  zprnch(ji,jj,jk) = zprnch(ji,jj,jk) * zmixnano(ji,jj) 
     258                  zprdch(ji,jj,jk) = zprdch(ji,jj,jk) * zmixdiat(ji,jj) 
    216259               ENDIF 
    217260            END DO 
    218261         END DO 
    219262      END DO 
    220  
    221263 
    222264!CDIR NOVERRCHK 
     
    230272                  !     Computation of the various production terms for nanophyto. 
    231273                  zetot2 = enano(ji,jj,jk) * zstrn(ji,jj) 
    232                   zmax = MAX( 0.1, xlimphy(ji,jj,jk) ) 
    233                   zpislopen = zpislopead(ji,jj,jk)          & 
    234                   &         * trn(ji,jj,jk,jpnch) / ( rtrn + trn(ji,jj,jk,jpphy) * 12.)         & 
    235                   &         / ( prmax(ji,jj,jk) * rday * zmax + rtrn ) 
    236  
    237                   zprbiochl = prmax(ji,jj,jk) * (  1.- EXP( -zpislopen * zetot2 )  ) 
    238  
    239274                  zprorca(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * trn(ji,jj,jk,jpphy) * rfact2 
    240  
    241                   zpronew(ji,jj,jk) = zprorca(ji,jj,jk) * xnanono3(ji,jj,jk)    & 
    242                   &             / ( xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) + rtrn ) 
    243                   zprod = rday * zprorca(ji,jj,jk) * zprbiochl * trn(ji,jj,jk,jpphy) *zmax 
    244  
    245                   zprofen(ji,jj,jk) = (fecnm)**2 * zprod / chlcnm            & 
    246                   &              / (  zpislopead(ji,jj,jk) * zetot2 * trn(ji,jj,jk,jpnfe) + rtrn  ) 
    247  
    248                   zprochln(ji,jj,jk) = chlcnm * 144. * zprod                  & 
    249                   &              / (  zpislopead(ji,jj,jk) * zetot2 * trn(ji,jj,jk,jpnch) + rtrn  ) 
     275                  zpronew(ji,jj,jk) = zprorca(ji,jj,jk) * xnanono3(ji,jj,jk) / ( xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) + rtrn ) 
     276                  ! 
     277                  zprod = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * xlimphy(ji,jj,jk) 
     278                  zprochln(ji,jj,jk) = chlcmin * 12. * zprorca(ji,jj,jk) 
     279                  zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + chlcnm * 12. * zprod / (  zpislopead(ji,jj,jk) * zetot2 +rtrn) 
     280                  ! 
     281                  zratio = trn(ji,jj,jk,jpnfe) / ( trn(ji,jj,jk,jpphy) + rtrn ) 
     282                  zratio = zratio / fecnm  
     283                  zmax   = MAX( 0., ( 1. - zratio ) / ABS( 1.05 - zratio ) )  
     284                  zprofen(ji,jj,jk) = fecnm * prmax(ji,jj,jk)  & 
     285                  &             * ( 4. - 4.5 * xlimnfe(ji,jj,jk) / ( xlimnfe(ji,jj,jk) + 0.5 ) )    & 
     286                  &             * trn(ji,jj,jk,jpfer) / ( trn(ji,jj,jk,jpfer) + concnfe(ji,jj,jk) )  & 
     287                  &             * zmax * trn(ji,jj,jk,jpphy) * rfact2 
    250288               ENDIF 
    251289            END DO 
     
    262300                  !  Computation of the various production terms for diatoms 
    263301                  zetot2 = ediat(ji,jj,jk) * zstrn(ji,jj) 
    264                   zmax = MAX( 0.1, xlimdia(ji,jj,jk) ) 
    265                   zpislope2n = zpislopead2(ji,jj,jk) * trn(ji,jj,jk,jpdch)        & 
    266                   &           / ( rtrn + trn(ji,jj,jk,jpdia) * 12.)        & 
    267                   &           / ( prmax(ji,jj,jk) * rday * zmax + rtrn ) 
    268  
    269                   zprdiachl = prmax(ji,jj,jk) * (  1.- EXP( -zetot2 * zpislope2n )  ) 
    270  
    271302                  zprorcad(ji,jj,jk) = zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * trn(ji,jj,jk,jpdia) * rfact2 
    272  
    273                   zpronewd(ji,jj,jk) = zprorcad(ji,jj,jk) * xdiatno3(ji,jj,jk)     & 
    274                   &              / ( xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) + rtrn ) 
    275  
    276                   zprod = rday * zprorcad(ji,jj,jk) * zprdiachl * trn(ji,jj,jk,jpdia) * zmax 
    277  
    278                   zprofed(ji,jj,jk) = (fecdm)**2 * zprod / chlcdm                   & 
    279                   &              / ( zpislopead2(ji,jj,jk) * zetot2 * trn(ji,jj,jk,jpdfe) + rtrn ) 
    280  
    281                   zprochld(ji,jj,jk) = chlcdm * 144. * zprod       & 
    282                   &              / ( zpislopead2(ji,jj,jk) * zetot2 * trn(ji,jj,jk,jpdch) + rtrn ) 
    283  
     303                  zpronewd(ji,jj,jk) = zprorcad(ji,jj,jk) * xdiatno3(ji,jj,jk) / ( xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) + rtrn ) 
     304                  ! 
     305                  zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * xlimdia(ji,jj,jk) 
     306                  zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk) 
     307                  zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + chlcdm * 12. * zprod / ( zpislopead2(ji,jj,jk) * zetot2 +rtrn ) 
     308                  ! 
     309                  zratio = trn(ji,jj,jk,jpdfe) / ( trn(ji,jj,jk,jpdia) + rtrn ) 
     310                  zratio = zratio / fecdm  
     311                  zmax   = MAX( 0., ( 1. - zratio ) / ABS( 1.05 - zratio ) )  
     312                  zprofed(ji,jj,jk) = fecdm * prmax(ji,jj,jk)  & 
     313                  &             * ( 4. - 4.5 * xlimdfe(ji,jj,jk) / ( xlimdfe(ji,jj,jk) + 0.5 ) )    & 
     314                  &             * trn(ji,jj,jk,jpfer) / ( trn(ji,jj,jk,jpfer) + concdfe(ji,jj,jk) )  & 
     315                  &             * zmax * trn(ji,jj,jk,jpdia) * rfact2 
    284316               ENDIF 
    285317            END DO 
    286318         END DO 
    287319      END DO 
    288       ! 
    289320 
    290321      !   Update the arrays TRA which contain the biological sources and sinks 
     
    304335              tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) + zprofed(ji,jj,jk) * texcret2 
    305336              tra(ji,jj,jk,jpbsi) = tra(ji,jj,jk,jpbsi) + zprorcad(ji,jj,jk) * zysopt(ji,jj,jk) * texcret2 
    306               tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + & 
    307               &                     excret2 * zprorcad(ji,jj,jk) + excret * zprorca(ji,jj,jk) 
     337              tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + excret2 * zprorcad(ji,jj,jk) + excret * zprorca(ji,jj,jk) 
    308338              tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + o2ut * ( zproreg + zproreg2) & 
    309               &                    + ( o2ut + o2nit ) * ( zpronew(ji,jj,jk) + zpronewd(ji,jj,jk) ) 
    310               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) & 
    311               &                     - texcret * zprofen(ji,jj,jk) - texcret2 * zprofed(ji,jj,jk) 
    312               tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) & 
    313               &                     - texcret2 * zprorcad(ji,jj,jk) * zysopt(ji,jj,jk) 
     339                 &                + ( o2ut + o2nit ) * ( zpronew(ji,jj,jk) + zpronewd(ji,jj,jk) ) 
     340              tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - texcret * zprofen(ji,jj,jk) - texcret2 * zprofed(ji,jj,jk) 
     341              tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) - texcret2 * zprorcad(ji,jj,jk) * zysopt(ji,jj,jk) 
    314342              tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprorca(ji,jj,jk) - zprorcad(ji,jj,jk) 
    315               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) & 
    316               &                    + rno3 * ( zpronew(ji,jj,jk) + zpronewd(ji,jj,jk) ) 
     343              tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * ( zpronew(ji,jj,jk) + zpronewd(ji,jj,jk) ) & 
     344                 &                                      - rno3 * ( zproreg + zproreg2 ) 
    317345          END DO 
    318346        END DO 
     
    320348 
    321349     ! Total primary production per year 
    322  
    323 #if defined key_degrad 
    324      tpp = tpp + glob_sum( ( zprorca(:,:,:) + zprorcad(:,:,:) ) * cvol(:,:,:) * facvol(:,:,:) ) 
    325 #else 
    326      tpp = tpp + glob_sum( ( zprorca(:,:,:) + zprorcad(:,:,:) ) * cvol(:,:,:) ) 
    327 #endif 
     350     IF( lk_degrad )  THEN 
     351        tpp = tpp + glob_sum( ( zprorca(:,:,:) + zprorcad(:,:,:) ) * cvol(:,:,:) * facvol(:,:,:) ) 
     352     ELSE 
     353        tpp = tpp + glob_sum( ( zprorca(:,:,:) + zprorcad(:,:,:) ) * cvol(:,:,:) ) 
     354     ENDIF 
    328355 
    329356     IF( kt == nitend .AND. jnt == nrdttrc .AND. lwp ) THEN 
     
    333360      ENDIF 
    334361 
    335 #if defined key_diatrc && ! defined key_iomput 
    336       !   Supplementary diagnostics 
     362#if defined key_diatrc  
    337363      zrfact2 = 1.e3 * rfact2r 
    338       trc3d(:,:,:,jp_pcs0_3d + 4)  = zprorca (:,:,:) * zrfact2 * tmask(:,:,:) 
    339       trc3d(:,:,:,jp_pcs0_3d + 5)  = zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) 
    340       trc3d(:,:,:,jp_pcs0_3d + 6)  = zpronew (:,:,:) * zrfact2 * tmask(:,:,:) 
    341       trc3d(:,:,:,jp_pcs0_3d + 7)  = zpronewd(:,:,:) * zrfact2 * tmask(:,:,:) 
    342       trc3d(:,:,:,jp_pcs0_3d + 8)  = zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) * zysopt(:,:,:) 
    343       trc3d(:,:,:,jp_pcs0_3d + 9)  = zprofed (:,:,:) * zrfact2 * tmask(:,:,:) 
    344 #  if ! defined key_kriest 
    345       trc3d(:,:,:,jp_pcs0_3d + 10) = zprofen (:,:,:) * zrfact2 * tmask(:,:,:) 
    346 #  endif 
    347 #endif 
    348  
    349 #if defined key_diatrc && defined key_iomput 
    350       zrfact2 = 1.e3 * rfact2r 
    351       IF ( jnt == nrdttrc ) then 
     364#if defined key_iomput 
     365      IF( jnt == nrdttrc ) THEN 
    352366         CALL iom_put( "PPPHY" , zprorca (:,:,:) * zrfact2 * tmask(:,:,:) )  ! primary production by nanophyto 
    353367         CALL iom_put( "PPPHY2", zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) )  ! primary production by diatom 
     
    357371         CALL iom_put( "PFeD"  , zprofed (:,:,:) * zrfact2 * tmask(:,:,:) )  ! biogenic iron production by diatom 
    358372         CALL iom_put( "PFeN"  , zprofen (:,:,:) * zrfact2 * tmask(:,:,:) )  ! biogenic iron production by nanophyto 
    359       ENDIF 
     373#else 
     374         trc3d(:,:,:,jp_pcs0_3d + 4)  = zprorca (:,:,:) * zrfact2 * tmask(:,:,:) 
     375         trc3d(:,:,:,jp_pcs0_3d + 5)  = zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) 
     376         trc3d(:,:,:,jp_pcs0_3d + 6)  = zpronew (:,:,:) * zrfact2 * tmask(:,:,:) 
     377         trc3d(:,:,:,jp_pcs0_3d + 7)  = zpronewd(:,:,:) * zrfact2 * tmask(:,:,:) 
     378         trc3d(:,:,:,jp_pcs0_3d + 8)  = zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) * zysopt(:,:,:) 
     379         trc3d(:,:,:,jp_pcs0_3d + 9)  = zprofed (:,:,:) * zrfact2 * tmask(:,:,:) 
     380#  if ! defined key_kriest 
     381         trc3d(:,:,:,jp_pcs0_3d + 10) = zprofen (:,:,:) * zrfact2 * tmask(:,:,:) 
     382#  endif 
     383#endif 
    360384#endif 
    361385 
     
    367391 
    368392      IF(  wrk_not_released(2, 1,2,3)                          .OR.  & 
    369            wrk_not_released(3, 2,3,4,5,6,7,8,9,10,11,12,13,14)   )   & 
     393           wrk_not_released(3, 2,3,4,5,6,7,8,9,10,11,12,13,14,15)   )   & 
    370394           CALL ctl_stop('p4z_prod: failed to release workspace arrays') 
     395      ! 
     396      DEALLOCATE( zysopt ) 
    371397      ! 
    372398   END SUBROUTINE p4z_prod 
     
    384410      !! ** input   :   Namelist nampisprod 
    385411      !!---------------------------------------------------------------------- 
    386       NAMELIST/nampisprod/ pislope, pislope2, excret, excret2, chlcnm, chlcdm,   & 
    387          &              fecnm, fecdm, grosip 
     412      ! 
     413      NAMELIST/nampisprod/ pislope, pislope2, bresp, excret, excret2,  & 
     414         &                 chlcnm, chlcdm, chlcmin, fecnm, fecdm, grosip 
    388415      !!---------------------------------------------------------------------- 
    389416 
    390       REWIND( numnat )                     ! read numnat 
    391       READ  ( numnat, nampisprod ) 
     417      REWIND( numnatp )                     ! read numnat 
     418      READ  ( numnatp, nampisprod ) 
    392419 
    393420      IF(lwp) THEN                         ! control print 
     
    399426         WRITE(numout,*) '    excretion ratio of nanophytoplankton      excret    =', excret 
    400427         WRITE(numout,*) '    excretion ratio of diatoms                excret2   =', excret2 
     428         WRITE(numout,*) '    basal respiration in phytoplankton        bresp     =', bresp 
    401429         WRITE(numout,*) '    P-I slope  for diatoms                    pislope2  =', pislope2 
    402430         WRITE(numout,*) '    Minimum Chl/C in nanophytoplankton        chlcnm    =', chlcnm 
    403431         WRITE(numout,*) '    Minimum Chl/C in diatoms                  chlcdm    =', chlcdm 
     432         WRITE(numout,*) '    Maximum Chl/C in phytoplankton            chlcmin   =', chlcmin 
    404433         WRITE(numout,*) '    Maximum Fe/C in nanophytoplankton         fecnm     =', fecnm 
    405434         WRITE(numout,*) '    Minimum Fe/C in diatoms                   fecdm     =', fecdm 
    406435      ENDIF 
    407436      ! 
    408       rday1     = 0.6 / rday  
    409       texcret   = 1.0 - excret 
    410       texcret2  = 1.0 - excret2 
    411       tpp       = 0. 
     437      r1_rday   = 1._wp / rday  
     438      r1_bresp  = bresp * r1_rday 
     439      texcret   = 1._wp - excret 
     440      texcret2  = 1._wp - excret2 
     441      tpp       = 0._wp 
    412442      ! 
    413443   END SUBROUTINE p4z_prod_init 
     
    418448      !!                     ***  ROUTINE p4z_prod_alloc  *** 
    419449      !!---------------------------------------------------------------------- 
    420       ALLOCATE( prmax(jpi,jpj,jpk), STAT=p4z_prod_alloc ) 
     450      ALLOCATE( prmax(jpi,jpj,jpk), quotan(jpi,jpj,jpk), quotad(jpi,jpj,jpk), STAT = p4z_prod_alloc ) 
    421451      ! 
    422452      IF( p4z_prod_alloc /= 0 ) CALL ctl_warn('p4z_prod_alloc : failed to allocate arrays.') 
Note: See TracChangeset for help on using the changeset viewer.