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 15459 for NEMO/trunk/src/TOP/PISCES/P4Z/p4zprod.F90 – NEMO

Ignore:
Timestamp:
2021-10-29T10:19:18+02:00 (2 years ago)
Author:
cetlod
Message:

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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/TOP/PISCES/P4Z/p4zprod.F90

    r15090 r15459  
    22   !!====================================================================== 
    33   !!                         ***  MODULE p4zprod  *** 
    4    !! TOP :  Growth Rate of the two phytoplanktons groups  
     4   !! TOP :  Growth Rate of the two phytoplankton groups of PISCES  
    55   !!====================================================================== 
    66   !! History :   1.0  !  2004     (O. Aumont) Original code 
     
    2424   PUBLIC   p4z_prod         ! called in p4zbio.F90 
    2525   PUBLIC   p4z_prod_init    ! called in trcsms_pisces.F90 
    26    PUBLIC   p4z_prod_alloc 
    27  
    28    REAL(wp), PUBLIC ::   pislopen     !: 
    29    REAL(wp), PUBLIC ::   pisloped     !: 
    30    REAL(wp), PUBLIC ::   xadap        !: 
    31    REAL(wp), PUBLIC ::   excretn      !: 
    32    REAL(wp), PUBLIC ::   excretd      !: 
    33    REAL(wp), PUBLIC ::   bresp        !: 
    34    REAL(wp), PUBLIC ::   chlcnm       !: 
    35    REAL(wp), PUBLIC ::   chlcdm       !: 
    36    REAL(wp), PUBLIC ::   chlcmin      !: 
    37    REAL(wp), PUBLIC ::   fecnm        !: 
    38    REAL(wp), PUBLIC ::   fecdm        !: 
    39    REAL(wp), PUBLIC ::   grosip       !: 
     26   PUBLIC   p4z_prod_alloc   ! called in trcini_pisces.F90 
     27 
     28   REAL(wp), PUBLIC ::   pislopen     !:  P-I slope of nanophytoplankton 
     29   REAL(wp), PUBLIC ::   pisloped     !:  P-I slope of diatoms 
     30   REAL(wp), PUBLIC ::   xadap        !:  Adaptation factor to low light  
     31   REAL(wp), PUBLIC ::   excretn      !:  Excretion ratio of nanophyto 
     32   REAL(wp), PUBLIC ::   excretd      !:  Excretion ratio of diatoms 
     33   REAL(wp), PUBLIC ::   bresp        !:  Basal respiration rate 
     34   REAL(wp), PUBLIC ::   chlcnm       !:  Maximum Chl/C ratio of nano 
     35   REAL(wp), PUBLIC ::   chlcdm       !:  Maximum Chl/C ratio of diatoms 
     36   REAL(wp), PUBLIC ::   chlcmin      !:  Minimum Chl/C ratio of phytoplankton 
     37   REAL(wp), PUBLIC ::   fecnm        !:  Maximum Fe/C ratio of nano 
     38   REAL(wp), PUBLIC ::   fecdm        !:  Maximum Fe/C ratio of diatoms 
     39   REAL(wp), PUBLIC ::   grosip       !:  Mean Si/C ratio of diatoms 
    4040 
    4141   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   quotan   !: proxy of N quota in Nanophyto 
    42    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   quotad   !: proxy of N quota in diatomee 
     42   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   quotad   !: proxy of N quota in diatoms 
    4343    
    4444   REAL(wp) ::   r1_rday    ! 1 / rday 
     
    6060      !!                     ***  ROUTINE p4z_prod  *** 
    6161      !! 
    62       !! ** Purpose :   Compute the phytoplankton production depending on 
    63       !!              light, temperature and nutrient availability 
    64       !! 
    65       !! ** Method  : - ??? 
     62      !! ** Purpose :   Computes phytoplankton production depending on 
     63      !!                light, temperature and nutrient availability 
     64      !!                Computes also the uptake of Iron and Si as well  
     65      !!                as the chlorophyll content of the cells 
     66      !!                PISCES relies on a mixed Monod-Quota formalism  
    6667      !!--------------------------------------------------------------------- 
    6768      INTEGER, INTENT(in) ::   kt, knt   ! 
     
    7071      INTEGER  ::   ji, jj, jk 
    7172      REAL(wp) ::   zsilfac, znanotot, zdiattot, zconctemp, zconctemp2 
    72       REAL(wp) ::   zratio, zmax, zsilim, ztn, zadap, zlim, zsilfac2, zsiborn 
    73       REAL(wp) ::   zprod, zproreg, zproreg2, zprochln, zprochld 
    74       REAL(wp) ::   zmaxday, zdocprod, zpislopen, zpisloped 
    75       REAL(wp) ::   zmxltst, zmxlday 
    76       REAL(wp) ::   zrum, zcodel, zargu, zval, zfeup, chlcnm_n, chlcdm_n 
    77       REAL(wp) ::   zfact 
     73      REAL(wp) ::   zratio, zmax, zsilim, ztn, zadap, zlim, zsiborn 
     74      REAL(wp) ::   zpptot, zpnewtot, zpregtot, zprochln, zprochld 
     75      REAL(wp) ::   zproddoc, zprodsil, zprodfer, zprodlig 
     76      REAL(wp) ::   zpislopen, zpisloped, zfact 
     77      REAL(wp) ::   zratiosi, zmaxsi, zlimfac, zsizetmp, zfecnm, zfecdm 
     78      REAL(wp) ::   zprod, zval 
    7879      CHARACTER (len=25) :: charout 
    79       REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zw2d 
    80       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d 
    81       REAL(wp), DIMENSION(jpi,jpj    ) :: zstrn, zmixnano, zmixdiat 
    8280      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprmaxn,zprmaxd 
    8381      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpislopeadn, zpislopeadd, zysopt   
    84       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprdia, zprbio, zprdch, zprnch    
     82      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprdia, zprbio, zprchld, zprchln    
    8583      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprorcan, zprorcad, zprofed, zprofen 
    8684      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpronewn, zpronewd 
    8785      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmxl_fac, zmxl_chl 
    88       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpligprod1, zpligprod2 
     86      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zpligprod 
    8987      !!--------------------------------------------------------------------- 
    9088      ! 
     
    9391      !  Allocate temporary workspace 
    9492      ! 
    95       zprorcan  (:,:,:) = 0._wp ; zprorcad  (:,:,:) = 0._wp ; zprofed (:,:,:) = 0._wp 
    96       zprofen   (:,:,:) = 0._wp ; zysopt    (:,:,:) = 0._wp 
    97       zpronewn  (:,:,:) = 0._wp ; zpronewd  (:,:,:) = 0._wp ; zprdia  (:,:,:) = 0._wp 
    98       zprbio    (:,:,:) = 0._wp ; zprdch    (:,:,:) = 0._wp ; zprnch  (:,:,:) = 0._wp  
    99       zmxl_fac  (:,:,:) = 0._wp ; zmxl_chl  (:,:,:) = 0._wp  
    100       zpligprod1(:,:,:) = 0._wp ; zpligprod2(:,:,:) = 0._wp  
    101  
    102       ! Computation of the optimal production 
    103       zprmaxn(:,:,:) = 0.8_wp * r1_rday * tgfunc(:,:,:) 
     93      zprorcan(:,:,:) = 0._wp ; zprorcad(:,:,:) = 0._wp ; zprofed(:,:,:) = 0._wp 
     94      zprofen (:,:,:) = 0._wp ; zysopt  (:,:,:) = 0._wp 
     95      zpronewn(:,:,:) = 0._wp ; zpronewd(:,:,:) = 0._wp ; zprdia(:,:,:) = 0._wp 
     96      zprbio  (:,:,:) = 0._wp ; zprchld (:,:,:) = 0._wp ; zprchln(:,:,:) = 0._wp  
     97      zmxl_fac(:,:,:) = 0._wp ; zmxl_chl(:,:,:) = 0._wp  
     98 
     99      ! Computation of the maximimum production. Based on a Q10 description 
     100      ! of the thermal dependency 
     101      ! Parameters are taken from Bissinger et al. (2008) 
     102      zprmaxn(:,:,:) = 0.65_wp * r1_rday * tgfunc(:,:,:) 
    104103      zprmaxd(:,:,:) = zprmaxn(:,:,:) 
    105104 
    106       ! compute the day length depending on latitude and the day 
    107       zrum = REAL( nday_year - 80, wp ) / REAL( nyear_len(1), wp ) 
    108       zcodel = ASIN(  SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp )  ) 
    109  
    110       ! day length in hours 
    111       zstrn(:,:) = 0. 
    112       DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    113          zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad ) 
    114          zargu = MAX( -1., MIN(  1., zargu ) ) 
    115          zstrn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. ) 
    116       END_2D 
    117  
    118       ! Impact of the day duration and light intermittency on phytoplankton growth 
    119       DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1) 
    120          IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 
    121             zval = MAX( 1., zstrn(ji,jj) ) 
    122             IF( gdept(ji,jj,jk,Kmm) <= hmld(ji,jj) ) THEN 
    123                zval = zval * MIN(1., heup_01(ji,jj) / ( hmld(ji,jj) + rtrn )) 
     105      ! Intermittency is supposed to have a similar effect on production as  
     106      ! day length (Shatwell et al., 2012). The correcting factor is zmxl_fac.  
     107      ! zmxl_chl is the fractional day length and is used to compute the mean 
     108      ! PAR during daytime. The effect of mixing is computed using the  
     109      ! absolute light level definition of the euphotic zone 
     110      ! -------------------------------------------------------------------------  
     111      IF ( ln_p4z_dcyc ) THEN    ! Diurnal cycle in PISCES 
     112 
     113         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1) 
     114            IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 
     115               IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN 
     116                  zval = MIN(1., heup_01(ji,jj) / ( hmld(ji,jj) + rtrn )) 
     117               ENDIF 
     118               zmxl_chl(ji,jj,jk) = zval / 24. 
     119               zmxl_fac(ji,jj,jk) = 1.0 - exp( -0.26 * zval ) 
    124120            ENDIF 
    125             zmxl_chl(ji,jj,jk) = zval / 24. 
    126             zmxl_fac(ji,jj,jk) = 1.5 * zval / ( 12. + zval ) 
    127          ENDIF 
    128       END_3D 
     121         END_3D 
     122  
     123      ELSE ! No diurnal cycle in PISCES 
     124         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1) 
     125            IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 
     126               zval = MAX( 1., strn(ji,jj) ) 
     127               IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN 
     128                  zval = zval * MIN(1., heup_01(ji,jj) / ( hmld(ji,jj) + rtrn )) 
     129               ENDIF 
     130               zmxl_chl(ji,jj,jk) = zval / 24. 
     131               zmxl_fac(ji,jj,jk) = 1.0 - exp( -0.26 * zval ) 
     132            ENDIF 
     133         END_3D 
     134 
     135      ENDIF 
    129136 
    130137      zprbio(:,:,:) = zprmaxn(:,:,:) * zmxl_fac(:,:,:) 
    131138      zprdia(:,:,:) = zprmaxd(:,:,:) * zmxl_fac(:,:,:) 
    132139 
    133       ! Maximum light intensity 
    134       WHERE( zstrn(:,:) < 1.e0 ) zstrn(:,:) = 24. 
    135  
    136       ! Computation of the P-I slope for nanos and diatoms 
     140      ! The formulation proposed by Geider et al. (1997) has been modified  
     141      ! to exclude the effect of nutrient limitation and temperature in the PI 
     142      ! curve following Vichi et al. (2007) 
     143      ! ----------------------------------------------------------------------- 
    137144      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1) 
    138145         IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 
     
    142149            zconctemp2  = tr(ji,jj,jk,jpdia,Kbb) - zconctemp 
    143150            ! 
     151            ! The initial slope of the PI curve can be increased for nano 
     152            ! to account for photadaptation, for instance in the DCM 
     153            ! This parameterization is adhoc and should be either  
     154            ! improved or removed in future versions of the model 
     155 
     156            ! Nanophytoplankton 
    144157            zpislopeadn(ji,jj,jk) = pislopen * ( 1.+ zadap  * EXP( -0.25 * enano(ji,jj,jk) ) )  & 
    145158            &                   * tr(ji,jj,jk,jpnch,Kbb) /( tr(ji,jj,jk,jpphy,Kbb) * 12. + rtrn) 
    146             ! 
     159 
     160            ! Diatoms 
    147161            zpislopeadd(ji,jj,jk) = (pislopen * zconctemp2 + pisloped * zconctemp) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn )   & 
    148162            &                   * tr(ji,jj,jk,jpdch,Kbb) /( tr(ji,jj,jk,jpdia,Kbb) * 12. + rtrn) 
     
    153167         IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 
    154168             ! Computation of production function for Carbon 
    155              !  --------------------------------------------- 
     169             ! Actual light levels are used here  
     170             ! ---------------------------------------------- 
    156171             zpislopen = zpislopeadn(ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) & 
    157172             &            * zmxl_fac(ji,jj,jk) * rday + rtrn) 
     
    160175             zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) )  ) 
    161176             zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) )  ) 
     177 
    162178             !  Computation of production function for Chlorophyll 
    163              !-------------------------------------------------- 
     179             !  Mean light level in the mixed layer (when appropriate) 
     180             !  is used here (acclimation is in general slower than  
     181             !  the characteristic time scales of vertical mixing) 
     182             !  ------------------------------------------------------ 
    164183             zpislopen = zpislopeadn(ji,jj,jk) / ( zprmaxn(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn ) 
    165184             zpisloped = zpislopeadd(ji,jj,jk) / ( zprmaxd(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn ) 
    166              zprnch(ji,jj,jk) = zprmaxn(ji,jj,jk) * ( 1.- EXP( -zpislopen * enanom(ji,jj,jk) ) ) 
    167              zprdch(ji,jj,jk) = zprmaxd(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediatm(ji,jj,jk) ) ) 
     185             zprchln(ji,jj,jk) = zprmaxn(ji,jj,jk) * ( 1.- EXP( -zpislopen * enanom(ji,jj,jk) ) ) 
     186             zprchld(ji,jj,jk) = zprmaxd(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediatm(ji,jj,jk) ) ) 
    168187         ENDIF 
    169188      END_3D 
    170189 
    171       !  Computation of a proxy of the N/C ratio 
    172       !  --------------------------------------- 
     190      !  Computation of a proxy of the N/C quota from nutrient limitation  
     191      !  and light limitation. Steady state is assumed to allow the computation 
     192      !  ---------------------------------------------------------------------- 
    173193      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1) 
    174194          zval = MIN( xnanopo4(ji,jj,jk), ( xnanonh4(ji,jj,jk) + xnanono3(ji,jj,jk) ) )   & 
    175195          &      * zprmaxn(ji,jj,jk) / ( zprbio(ji,jj,jk) + rtrn ) 
    176           quotan(ji,jj,jk) = MIN( 1., 0.2 + 0.8 * zval ) 
     196          quotan(ji,jj,jk) = MIN( 1., 0.3 + 0.7 * zval ) 
    177197          zval = MIN( xdiatpo4(ji,jj,jk), ( xdiatnh4(ji,jj,jk) + xdiatno3(ji,jj,jk) ) )   & 
    178198          &      * zprmaxd(ji,jj,jk) / ( zprdia(ji,jj,jk) + rtrn ) 
    179           quotad(ji,jj,jk) = MIN( 1., 0.2 + 0.8 * zval ) 
     199          quotad(ji,jj,jk) = MIN( 1., 0.3 + 0.7 * zval ) 
    180200      END_3D 
    181201 
     
    184204 
    185205          IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 
    186              !    Si/C of diatoms 
    187              !    ------------------------ 
    188              !    Si/C increases with iron stress and silicate availability 
    189              !    Si/C is arbitrariliy increased for very high Si concentrations 
    190              !    to mimic the very high ratios observed in the Southern Ocean (silpot2) 
     206             ! Si/C of diatoms 
     207             ! ------------------------ 
     208             ! Si/C increases with iron stress and silicate availability 
     209             ! Si/C is arbitrariliy increased for very high Si concentrations 
     210             ! to mimic the very high ratios observed in the Southern Ocean (zsilfac) 
     211             ! A parameterization derived from Flynn (2003) is used for the control 
     212             ! when Si is not limiting which is similar to the parameterisation 
     213             ! proposed by Gurney and Davidson (1999). 
     214             ! ----------------------------------------------------------------------- 
    191215            zlim  = tr(ji,jj,jk,jpsil,Kbb) / ( tr(ji,jj,jk,jpsil,Kbb) + xksi1 ) 
    192             zsilim = MIN( zprdia(ji,jj,jk) / ( zprmaxd(ji,jj,jk) + rtrn ), xlimsi(ji,jj,jk) ) 
    193             zsilfac = 4.4 * EXP( -4.23 * zsilim ) * MAX( 0.e0, MIN( 1., 2.2 * ( zlim - 0.5 ) )  ) + 1.e0 
     216            zsilim = xlimdia(ji,jj,jk) * zprdia(ji,jj,jk) / ( zprmaxd(ji,jj,jk) + rtrn ) 
    194217            zsiborn = tr(ji,jj,jk,jpsil,Kbb) * tr(ji,jj,jk,jpsil,Kbb) * tr(ji,jj,jk,jpsil,Kbb) 
    195218            IF (gphit(ji,jj) < -30 ) THEN 
    196               zsilfac2 = 1. + 2. * zsiborn / ( zsiborn + xksi2**3 ) 
     219              zsilfac = 1. + 2. * zsiborn / ( zsiborn + xksi2**3 ) 
    197220            ELSE 
    198               zsilfac2 = 1. +      zsiborn / ( zsiborn + xksi2**3 ) 
     221              zsilfac = 1. +      zsiborn / ( zsiborn + xksi2**3 ) 
    199222            ENDIF 
    200             zysopt(ji,jj,jk) = grosip * zlim * zsilfac * zsilfac2 
     223            zratiosi = 1.0 - tr(ji,jj,jk,jpdsi,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn ) / ( zsilfac * grosip * 3.0 + rtrn ) 
     224            zratiosi = MAX(0., MIN(1.0, zratiosi) ) 
     225            zmaxsi  = (1.0 + 0.1**4) * zratiosi**4 / ( zratiosi**4 + 0.1**4 ) 
     226            IF( xlimsi(ji,jj,jk) /= xlimdia(ji,jj,jk) ) THEN 
     227               zysopt(ji,jj,jk) = zlim * zsilfac * grosip * 1.0 * zmaxsi 
     228            ELSE 
     229               zysopt(ji,jj,jk) = zlim * zsilfac * grosip * 1.0 * zsilim**0.7 * zmaxsi 
     230            ENDIF 
    201231        ENDIF 
    202232      END_3D 
    203233 
    204       !  Mixed-layer effect on production  
    205       !  Sea-ice effect on production 
    206  
     234      ! Sea-ice effect on production 
     235      ! No production is assumed below sea ice 
     236      ! --------------------------------------  
    207237      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1) 
    208238         zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1. - fr_i(ji,jj) ) 
     
    210240      END_3D 
    211241 
    212       ! Computation of the various production terms  
     242      ! Computation of the various production  and nutrient uptake terms 
     243      ! --------------------------------------------------------------- 
    213244      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1) 
    214245         IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 
    215246            !  production terms for nanophyto. (C) 
    216247            zprorcan(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * tr(ji,jj,jk,jpphy,Kbb) * rfact2 
     248 
     249            !  New production (uptake of NO3) 
    217250            zpronewn(ji,jj,jk)  = zprorcan(ji,jj,jk)* xnanono3(ji,jj,jk) / ( xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) + rtrn ) 
    218251            ! 
    219             zratio = tr(ji,jj,jk,jpnfe,Kbb) / ( tr(ji,jj,jk,jpphy,Kbb) * fecnm + rtrn ) 
    220             zmax   = MAX( 0., ( 1. - zratio ) / ABS( 1.05 - zratio ) )  
    221             zprofen(ji,jj,jk) = fecnm * zprmaxn(ji,jj,jk) * ( 1.0 - fr_i(ji,jj) )  & 
    222             &             * ( 4. - 4.5 * xlimnfe(ji,jj,jk) / ( xlimnfe(ji,jj,jk) + 0.5 ) )    & 
    223             &             * biron(ji,jj,jk) / ( biron(ji,jj,jk) + concnfe(ji,jj,jk) )  & 
    224             &             * zmax * tr(ji,jj,jk,jpphy,Kbb) * rfact2 
    225             !  production terms for diatoms (C) 
     252            ! Size computation 
     253            ! Size is made a function of the limitation of of phytoplankton growth 
     254            ! Strongly limited cells are supposed to be smaller. sizena is the  
     255            ! size at time step t+1 and is thus updated at the end of the  
     256            ! current time step 
     257            ! -------------------------------------------------------------------- 
     258            zlimfac = xlimphy(ji,jj,jk) * zprchln(ji,jj,jk) / ( zprmaxn(ji,jj,jk) + rtrn ) 
     259            zsizetmp = 1.0 + 1.3 * ( xsizern - 1.0 ) * zlimfac**3/(0.3 + zlimfac**3) 
     260            sizena(ji,jj,jk) = min(xsizern, max( sizena(ji,jj,jk), zsizetmp ) ) 
     261 
     262            ! Iron uptake rates of nanophytoplankton. Upregulation is   
     263            ! not parameterized at low iron concentrations as observations 
     264            ! do not suggest it for accimated cells. Uptake is 
     265            ! downregulated when the quota is close to the maximum quota 
     266            zfecnm = xqfuncfecn(ji,jj,jk) + ( fecnm - xqfuncfecn(ji,jj,jk) ) * ( xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) ) 
     267            zratio = 1.0 - MIN(1.0,tr(ji,jj,jk,jpnfe,Kbb) / ( tr(ji,jj,jk,jpphy,Kbb) * zfecnm + rtrn ) ) 
     268            zmax   = MAX( 0., MIN( 1.0, zratio**2/ (0.05**2+zratio**2) ) )  
     269            zprofen(ji,jj,jk) = zfecnm * zprmaxn(ji,jj,jk) * ( 1.0 - fr_i(ji,jj) )  & 
     270            &          * (1. + 0.8 * xnanono3(ji,jj,jk) / ( rtrn + xnanono3(ji,jj,jk)  & 
     271            &          + xnanonh4(ji,jj,jk) ) * (1. - xnanofer(ji,jj,jk) ) )   & 
     272            &          * xnanofer(ji,jj,jk) * zmax * tr(ji,jj,jk,jpphy,Kbb) * rfact2 
     273            ! production terms of diatoms (C) 
    226274            zprorcad(ji,jj,jk) = zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * tr(ji,jj,jk,jpdia,Kbb) * rfact2 
     275 
     276            ! New production (uptake of NO3) 
    227277            zpronewd(ji,jj,jk) = zprorcad(ji,jj,jk) * xdiatno3(ji,jj,jk) / ( xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) + rtrn ) 
    228             ! 
    229             zratio = tr(ji,jj,jk,jpdfe,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) * fecdm + rtrn ) 
    230             zmax   = MAX( 0., ( 1. - zratio ) / ABS( 1.05 - zratio ) )  
    231             zprofed(ji,jj,jk) = fecdm * zprmaxd(ji,jj,jk) * ( 1.0 - fr_i(ji,jj) )  & 
    232             &             * ( 4. - 4.5 * xlimdfe(ji,jj,jk) / ( xlimdfe(ji,jj,jk) + 0.5 ) )    & 
    233             &             * biron(ji,jj,jk) / ( biron(ji,jj,jk) + concdfe(ji,jj,jk) )  & 
    234             &             * zmax * tr(ji,jj,jk,jpdia,Kbb) * rfact2 
     278 
     279            ! Size computation 
     280            ! Size is made a function of the limitation of of phytoplankton growth 
     281            ! Strongly limited cells are supposed to be smaller. sizeda is 
     282            ! size at time step t+1 and is thus updated at the end of the  
     283            ! current time step.  
     284            ! -------------------------------------------------------------------- 
     285            zlimfac = zprchld(ji,jj,jk) * xlimdia(ji,jj,jk) / ( zprmaxd(ji,jj,jk) + rtrn ) 
     286            zsizetmp = 1.0 + 1.3 * ( xsizerd - 1.0 ) * zlimfac**3/(0.3 + zlimfac**3) 
     287            sizeda(ji,jj,jk) = min(xsizerd, max( sizeda(ji,jj,jk), zsizetmp ) ) 
     288 
     289            ! Iron uptake rates of diatoms. Upregulation is   
     290            ! not parameterized at low iron concentrations as observations 
     291            ! do not suggest it for accimated cells. Uptake is 
     292            ! downregulated when the quota is close to the maximum quota 
     293            zfecdm = xqfuncfecd(ji,jj,jk) + ( fecdm - xqfuncfecd(ji,jj,jk) ) * ( xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) ) 
     294            zratio = 1.0 - MIN(1.0, tr(ji,jj,jk,jpdfe,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) * zfecdm + rtrn ) ) 
     295            zmax   = MAX( 0., MIN( 1.0, zratio**2/ (0.05**2+zratio**2) ) )  
     296            zprofed(ji,jj,jk) = zfecdm * zprmaxd(ji,jj,jk) * (1.0 - fr_i(ji,jj) )  & 
     297            &          * (1. + 0.8 * xdiatno3(ji,jj,jk) / ( rtrn + xdiatno3(ji,jj,jk)  & 
     298            &          + xdiatnh4(ji,jj,jk) ) * (1. - xdiatfer(ji,jj,jk) ) )   & 
     299            &          * xdiatfer(ji,jj,jk) * zmax * tr(ji,jj,jk,jpdia,Kbb) * rfact2 
    235300         ENDIF 
    236301      END_3D 
    237302 
    238303      ! Computation of the chlorophyll production terms 
     304      ! The parameterization is taken from Geider et al. (1997) 
     305      ! ------------------------------------------------------- 
    239306      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1) 
    240307         IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 
    241308            !  production terms for nanophyto. ( chlorophyll ) 
    242309            znanotot = enanom(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn ) 
    243             zprod    = rday * zprorcan(ji,jj,jk) * zprnch(ji,jj,jk) * xlimphy(ji,jj,jk) 
     310            zprod    = rday * zprorcan(ji,jj,jk) * zprchln(ji,jj,jk) * xlimphy(ji,jj,jk) 
    244311            zprochln = chlcmin * 12. * zprorcan (ji,jj,jk) 
    245             chlcnm_n   = MIN ( chlcnm, ( chlcnm / (1. - 1.14 / 43.4 *ts(ji,jj,jk,jp_tem,Kmm))) * (1. - 1.14 / 43.4 * 20.)) 
    246             zprochln = zprochln + (chlcnm_n-chlcmin) * 12. * zprod / & 
     312            zprochln = zprochln + (chlcnm - chlcmin) * 12. * zprod / & 
    247313                                  & (  zpislopeadn(ji,jj,jk) * znanotot +rtrn) 
     314 
    248315            !  production terms for diatoms ( chlorophyll ) 
    249316            zdiattot = ediatm(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn ) 
    250             zprod    = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * xlimdia(ji,jj,jk) 
     317            zprod    = rday * zprorcad(ji,jj,jk) * zprchld(ji,jj,jk) * xlimdia(ji,jj,jk) 
    251318            zprochld = chlcmin * 12. * zprorcad(ji,jj,jk) 
    252             chlcdm_n   = MIN ( chlcdm, ( chlcdm / (1. - 1.14 / 43.4 * ts(ji,jj,jk,jp_tem,Kmm))) * (1. - 1.14 / 43.4 * 20.)) 
    253             zprochld = zprochld + (chlcdm_n-chlcmin) * 12. * zprod / & 
     319            zprochld = zprochld + (chlcdm - chlcmin) * 12. * zprod / & 
    254320                                  & ( zpislopeadd(ji,jj,jk) * zdiattot +rtrn ) 
     321 
    255322            !   Update the arrays TRA which contain the Chla sources and sinks 
    256323            tr(ji,jj,jk,jpnch,Krhs) = tr(ji,jj,jk,jpnch,Krhs) + zprochln * texcretn 
     
    262329      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1) 
    263330        IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 
    264            zproreg  = zprorcan(ji,jj,jk) - zpronewn(ji,jj,jk) 
    265            zproreg2 = zprorcad(ji,jj,jk) - zpronewd(ji,jj,jk) 
    266            zdocprod = excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk) 
    267            tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) - zprorcan(ji,jj,jk) - zprorcad(ji,jj,jk) 
    268            tr(ji,jj,jk,jpno3,Krhs) = tr(ji,jj,jk,jpno3,Krhs) - zpronewn(ji,jj,jk) - zpronewd(ji,jj,jk) 
    269            tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) - zproreg - zproreg2 
     331           zpptot   = zprorcan(ji,jj,jk) + zprorcad(ji,jj,jk) 
     332           zpnewtot = zpronewn(ji,jj,jk) + zpronewd(ji,jj,jk) 
     333           zpregtot = zpptot - zpnewtot 
     334           zprodsil  = zprmaxd(ji,jj,jk) * zysopt(ji,jj,jk) * rfact2 * tr(ji,jj,jk,jpdia,Kbb) 
     335           zproddoc  = excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk) 
     336           zprodfer  = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk) 
     337           ! 
     338           tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) - zpptot 
     339           tr(ji,jj,jk,jpno3,Krhs) = tr(ji,jj,jk,jpno3,Krhs) - zpnewtot 
     340           tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) - zpregtot 
    270341           tr(ji,jj,jk,jpphy,Krhs) = tr(ji,jj,jk,jpphy,Krhs) + zprorcan(ji,jj,jk) * texcretn 
    271            tr(ji,jj,jk,jpnfe,Krhs) = tr(ji,jj,jk,jpnfe,Krhs) + zprofen(ji,jj,jk) * texcretn 
     342           tr(ji,jj,jk,jpnfe,Krhs) = tr(ji,jj,jk,jpnfe,Krhs) + zprofen(ji,jj,jk)  * texcretn 
    272343           tr(ji,jj,jk,jpdia,Krhs) = tr(ji,jj,jk,jpdia,Krhs) + zprorcad(ji,jj,jk) * texcretd 
    273            tr(ji,jj,jk,jpdfe,Krhs) = tr(ji,jj,jk,jpdfe,Krhs) + zprofed(ji,jj,jk) * texcretd 
    274            tr(ji,jj,jk,jpdsi,Krhs) = tr(ji,jj,jk,jpdsi,Krhs) + zprorcad(ji,jj,jk) * zysopt(ji,jj,jk) * texcretd 
    275            tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zdocprod 
    276            tr(ji,jj,jk,jpoxy,Krhs) = tr(ji,jj,jk,jpoxy,Krhs) + o2ut * ( zproreg + zproreg2) & 
    277            &                   + ( o2ut + o2nit ) * ( zpronewn(ji,jj,jk) + zpronewd(ji,jj,jk) ) 
     344           tr(ji,jj,jk,jpdfe,Krhs) = tr(ji,jj,jk,jpdfe,Krhs) + zprofed(ji,jj,jk)  * texcretd 
     345           tr(ji,jj,jk,jpdsi,Krhs) = tr(ji,jj,jk,jpdsi,Krhs) + zprodsil 
     346           tr(ji,jj,jk,jpsil,Krhs) = tr(ji,jj,jk,jpsil,Krhs) - zprodsil 
     347           tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zproddoc 
     348           tr(ji,jj,jk,jpoxy,Krhs) = tr(ji,jj,jk,jpoxy,Krhs) + o2ut * zpregtot + ( o2ut + o2nit ) * zpnewtot 
    278349           ! 
    279            zfeup = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk) 
    280            tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) - zfeup 
    281            tr(ji,jj,jk,jpsil,Krhs) = tr(ji,jj,jk,jpsil,Krhs) - texcretd * zprorcad(ji,jj,jk) * zysopt(ji,jj,jk) 
    282            tr(ji,jj,jk,jpdic,Krhs) = tr(ji,jj,jk,jpdic,Krhs) - zprorcan(ji,jj,jk) - zprorcad(ji,jj,jk) 
    283            tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + rno3 * ( zpronewn(ji,jj,jk) + zpronewd(ji,jj,jk) ) & 
    284            &                                         - rno3 * ( zproreg + zproreg2 ) 
     350           tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) - zprodfer 
     351           consfe3(ji,jj,jk)   = zprodfer * 75.0 / ( rtrn + ( plig(ji,jj,jk) + 75.0 * (1.0 - plig(ji,jj,jk) ) )   & 
     352           &                   * tr(ji,jj,jk,jpfer,Kbb) ) / rfact2 
     353           ! 
     354           tr(ji,jj,jk,jpdic,Krhs) = tr(ji,jj,jk,jpdic,Krhs) - zpptot 
     355           tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + rno3 * ( zpnewtot - zpregtot ) 
    285356        ENDIF 
    286357      END_3D 
    287      ! 
     358 
     359     ! Production and uptake of ligands by phytoplankton. This part is activated  
     360     ! when ln_ligand is set to .true. in the namelist. Ligand uptake is small  
     361     ! and based on the FeL model by Morel et al. (2008) and on the study of 
     362     ! Shaked et al. (2020) 
     363     ! ------------------------------------------------------------------------- 
    288364     IF( ln_ligand ) THEN 
    289          zpligprod1(:,:,:) = 0._wp    ;    zpligprod2(:,:,:) = 0._wp 
    290365         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1) 
    291366           IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 
    292               zdocprod = excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk) 
    293               zfeup    = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk) 
    294               tr(ji,jj,jk,jplgw,Krhs) = tr(ji,jj,jk,jplgw,Krhs) + zdocprod * ldocp - zfeup * plig(ji,jj,jk) * lthet 
    295               zpligprod1(ji,jj,jk) = zdocprod * ldocp 
    296               zpligprod2(ji,jj,jk) = zfeup * plig(ji,jj,jk) * lthet 
     367              zproddoc = excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk) 
     368              zprodfer = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk) 
     369              zprodlig = plig(ji,jj,jk) / ( rtrn + plig(ji,jj,jk) + 75.0 * (1.0 - plig(ji,jj,jk) ) ) * lthet 
     370              ! 
     371              tr(ji,jj,jk,jplgw,Krhs) = tr(ji,jj,jk,jplgw,Krhs) + zproddoc * ldocp - zprodfer * zprodlig 
    297372           ENDIF 
    298373         END_3D 
     
    300375 
    301376 
     377    ! Output of the diagnostics 
    302378    ! Total primary production per year 
    303379    IF( iom_use( "tintpp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  )  & 
     
    312388       CALL iom_put( "PPNEWD"  , zpronewd(:,:,:) * zfact * tmask(:,:,:)   ) ! new primary production by diatomes 
    313389       CALL iom_put( "PBSi"    , zprorcad(:,:,:) * zfact * tmask(:,:,:) * zysopt(:,:,:)  ) ! biogenic silica production 
    314        CALL iom_put( "PFeN"    , zprofen(:,:,:) * zfact * tmask(:,:,:)  ) ! biogenic iron production by nanophyto 
    315        CALL iom_put( "PFeD"    , zprofed(:,:,:) * zfact * tmask(:,:,:)  ) ! biogenic iron production by  diatomes 
    316        IF( ln_ligand ) THEN 
    317          CALL iom_put( "LPRODP"  , zpligprod1(:,:,:) * 1e9 * zfact * tmask(:,:,:) ) 
    318          CALL iom_put( "LDETP"   , zpligprod2(:,:,:) * 1e9 * zfact * tmask(:,:,:) ) 
     390       CALL iom_put( "PFeN"    , zprofen(:,:,:)  * zfact * tmask(:,:,:)  ) ! biogenic iron production by nanophyto 
     391       CALL iom_put( "PFeD"    , zprofed(:,:,:)  * zfact * tmask(:,:,:)  ) ! biogenic iron production by  diatomes 
     392       IF( ln_ligand .AND. ( iom_use( "LPRODP" ) .OR. iom_use( "LDETP" ) ) ) THEN 
     393           ALLOCATE(  zpligprod(jpi,jpj,jpk) ) 
     394           zpligprod(:,:,:) = excretd * zprorcad(:,:,:) + excretn * zprorcan(:,:,:) 
     395           CALL iom_put( "LPRODP"  , zpligprod(:,:,:) * ldocp * 1e9 * zfact * tmask(:,:,:) ) 
     396           ! 
     397           zpligprod(:,:,:) = ( texcretn * zprofen(:,:,:) + texcretd * zprofed(:,:,:) ) &  
     398             &                  * plig(:,:,:) / ( rtrn + plig(:,:,:) + 75.0 * (1.0 - plig(:,:,:) ) ) 
     399           CALL iom_put( "LDETP"   , zpligprod(:,:,:)  * lthet * 1e9 * zfact * tmask(:,:,:) ) 
     400           DEALLOCATE(  zpligprod ) 
    319401       ENDIF 
    320402       CALL iom_put( "Mumax"   , zprmaxn(:,:,:) * tmask(:,:,:)  ) ! Maximum growth rate 
     
    346428      !! ** Purpose :   Initialization of phytoplankton production parameters 
    347429      !! 
    348       !! ** Method  :   Read the nampisprod namelist and check the parameters 
     430      !! ** Method  :   Read the namp4zprod namelist and check the parameters 
    349431      !!      called at the first timestep (nittrc000) 
    350432      !! 
    351       !! ** input   :   Namelist nampisprod 
     433      !! ** input   :   Namelist namp4zprod 
    352434      !!---------------------------------------------------------------------- 
    353435      INTEGER ::   ios   ! Local integer 
    354436      ! 
     437      ! Namelist block 
    355438      NAMELIST/namp4zprod/ pislopen, pisloped, xadap, bresp, excretn, excretd,  & 
    356439         &                 chlcnm, chlcdm, chlcmin, fecnm, fecdm, grosip 
     
    365448      READ  ( numnatp_ref, namp4zprod, IOSTAT = ios, ERR = 901) 
    366449901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namp4zprod in reference namelist' ) 
     450 
    367451      READ  ( numnatp_cfg, namp4zprod, IOSTAT = ios, ERR = 902 ) 
    368452902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namp4zprod in configuration namelist' ) 
Note: See TracChangeset for help on using the changeset viewer.