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 12759 for NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zprod.F90 – NEMO

Ignore:
Timestamp:
2020-04-17T00:17:29+02:00 (4 years ago)
Author:
aumont
Message:

make parameterizations in PISCES-operationnal more similar to thos of PISCES-QUOTA (prey switching, optimal allocation, size, ...)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zprod.F90

    r12537 r12759  
    6868      REAL(wp) ::   zratio, zmax, zsilim, ztn, zadap, zlim, zsilfac2, zsiborn 
    6969      REAL(wp) ::   zprod, zproreg, zproreg2, zprochln, zprochld 
    70       REAL(wp) ::   zmaxday, zdocprod, zpislopen, zpisloped 
    71       REAL(wp) ::   zmxltst, zmxlday 
     70      REAL(wp) ::   zdocprod, zpislopen, zpisloped, zfact 
     71      REAL(wp) ::   zratiosi, zmaxsi, zlimfac, zsizetmp 
    7272      REAL(wp) ::   zrum, zcodel, zargu, zval, zfeup, chlcnm_n, chlcdm_n 
    73       REAL(wp) ::   zfact 
    7473      CHARACTER (len=25) :: charout 
    7574      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zw2d 
    7675      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d 
    77       REAL(wp), DIMENSION(jpi,jpj    ) :: zstrn, zmixnano, zmixdiat 
     76      REAL(wp), DIMENSION(jpi,jpj    ) :: zstrn 
    7877      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprmaxn,zprmaxd 
    7978      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpislopeadn, zpislopeadd, zysopt   
     
    129128                  ENDIF 
    130129                  zmxl_chl(ji,jj,jk) = zval / 24. 
    131                   zmxl_fac(ji,jj,jk) = 1.5 * zval / ( 12. + zval ) 
     130                  zmxl_fac(ji,jj,jk) = 1.0 - exp( -0.26 * zval ) 
    132131               ENDIF 
    133132            END DO 
     
    153152                  zconctemp   = MAX( 0.e0 , trb(ji,jj,jk,jpdia) - xsizedia ) 
    154153                  zconctemp2  = trb(ji,jj,jk,jpdia) - zconctemp 
    155                   ! 
     154 
    156155                  ! The initial slope of the PI curve can be increased for nano 
    157156                  ! to account for photadaptation, for instance in the DCM 
     
    197196                zval = MIN( xnanopo4(ji,jj,jk), ( xnanonh4(ji,jj,jk) + xnanono3(ji,jj,jk) ) )   & 
    198197                &      * zprmaxn(ji,jj,jk) / ( zprbio(ji,jj,jk) + rtrn ) 
    199                 quotan(ji,jj,jk) = MIN( 1., 0.2 + 0.8 * zval ) 
     198                quotan(ji,jj,jk) = MIN( 1., 0.3 + 0.7 * zval ) 
    200199                zval = MIN( xdiatpo4(ji,jj,jk), ( xdiatnh4(ji,jj,jk) + xdiatno3(ji,jj,jk) ) )   & 
    201200                &      * zprmaxd(ji,jj,jk) / ( zprdia(ji,jj,jk) + rtrn ) 
    202                 quotad(ji,jj,jk) = MIN( 1., 0.2 + 0.8 * zval ) 
    203             END DO 
    204          END DO 
    205       END DO 
    206  
    207  
    208       DO jk = 1, jpkm1 
    209          DO jj = 1, jpj 
    210             DO ji = 1, jpi 
    211  
    212                 IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 
     201                quotad(ji,jj,jk) = MIN( 1., 0.3 + 0.7 * zval ) 
     202            END DO 
     203         END DO 
     204      END DO 
     205 
     206 
     207      DO jk = 1, jpkm1 
     208         DO jj = 1, jpj 
     209            DO ji = 1, jpi 
     210 
     211               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 
    213212                   ! Si/C of diatoms 
    214213                   ! ------------------------ 
     
    217216                   ! to mimic the very high ratios observed in the Southern Ocean (zsilfac2) 
    218217                   ! ----------------------------------------------------------------------- 
    219                   zlim  = trb(ji,jj,jk,jpsil) / ( trb(ji,jj,jk,jpsil) + xksi1 ) 
    220                   zsilim = MIN( zprdia(ji,jj,jk) / ( zprmaxd(ji,jj,jk) + rtrn ), xlimsi(ji,jj,jk) ) 
    221                   zsilfac = 4.4 * EXP( -4.23 * zsilim ) * MAX( 0.e0, MIN( 1., 2.2 * ( zlim - 0.5 ) )  ) + 1.e0 
    222                   zsiborn = trb(ji,jj,jk,jpsil) * trb(ji,jj,jk,jpsil) * trb(ji,jj,jk,jpsil) 
    223                   IF (gphit(ji,jj) < -30 ) THEN 
    224                     zsilfac2 = 1. + 2. * zsiborn / ( zsiborn + xksi2**3 ) 
    225                   ELSE 
    226                     zsilfac2 = 1. +      zsiborn / ( zsiborn + xksi2**3 ) 
    227                   ENDIF 
    228                   zysopt(ji,jj,jk) = grosip * zlim * zsilfac * zsilfac2 
     218                 zlim  = trb(ji,jj,jk,jpsil) / ( trb(ji,jj,jk,jpsil) + xksi1 ) 
     219                 zsilim = xlimdia(ji,jj,jk) * zprdia(ji,jj,jk) / ( zprmaxd(ji,jj,jk) + rtrn ) 
     220                 zsiborn = trb(ji,jj,1,jpsil) * trb(ji,jj,1,jpsil) * trb(ji,jj,1,jpsil) 
     221                 IF (gphit(ji,jj) < -30 ) THEN 
     222                   zsilfac2 = 1. + 2. * zsiborn / ( zsiborn + xksi2**3 ) 
     223                 ELSE 
     224                   zsilfac2 = 1. +      zsiborn / ( zsiborn + xksi2**3 ) 
     225                 ENDIF 
     226                 zratiosi = 1.0 - trb(ji,jj,jk,jpdsi) / ( trb(ji,jj,jk,jpdia) + rtrn ) / ( zsilfac2 * grosip * 3.0 + rtrn ) 
     227                 zratiosi = MAX(0., MIN(1.0, zratiosi) ) 
     228                 zmaxsi  = (1.0 + 0.1**4) * zratiosi**4 / ( zratiosi**4 + 0.1**4 ) 
     229                 IF ( xlimsi(ji,jj,jk) /= xlimdia(ji,jj,jk) ) THEN 
     230                    zysopt(ji,jj,jk) = zlim * zsilfac2 * grosip * 1.0 * zmaxsi 
     231                 ELSE 
     232                    zysopt(ji,jj,jk) = zlim * zsilfac2 * grosip * 1.0 * zsilim**0.7 * zmaxsi 
     233                 ENDIF 
    229234              ENDIF 
    230235            END DO 
     
    255260                  !  New production (uptake of NO3) 
    256261                  zpronewn(ji,jj,jk)  = zprorcan(ji,jj,jk)* xnanono3(ji,jj,jk) / ( xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) + rtrn ) 
    257                   ! 
    258                   zratio = trb(ji,jj,jk,jpnfe) / ( trb(ji,jj,jk,jpphy) * fecnm + rtrn ) 
     262 
     263                  ! Size computation 
     264                  ! Size is made a function of the limitation of of phytoplankton growth 
     265                  ! Strongly limited cells are supposed to be smaller. sizena is the  
     266                  ! size at time step t+1 and is thus updated at the end of the  
     267                  ! current time step 
     268                  ! -------------------------------------------------------------------- 
     269                  zlimfac = xlimphy(ji,jj,jk) * zprchln(ji,jj,jk) / ( zprmaxn(ji,jj,jk) + rtrn ) 
     270                  zsizetmp = 1.0 + 1.3 * ( xsizern - 1.0 ) * zlimfac**3/(0.3 + zlimfac**3) 
     271                  sizena(ji,jj,jk) = min(xsizern, max( sizena(ji,jj,jk), zsizetmp ) ) 
    259272 
    260273                  !  Iron uptake rates of nanophytoplankton. Upregulation   
     
    262275                  !  formulation used in quota formulations. Uptake is downregulated 
    263276                  !  when the quota is close to the maximum quota  
    264                   zmax   = MAX( 0., ( 1. - zratio ) / ABS( 1.05 - zratio ) )  
     277                  zratio = 1.0 - MIN(1.0,trb(ji,jj,jk,jpnfe) / ( trb(ji,jj,jk,jpphy) * fecnm + rtrn ) ) 
     278                  zmax   = MAX( 0., MIN( 1.0, zratio**2/ (0.05**2+zratio**2) ) )  
    265279                  zprofen(ji,jj,jk) = fecnm * zprmaxn(ji,jj,jk) * ( 1.0 - fr_i(ji,jj) )  & 
    266                   &             * ( 4. - 4.5 * xlimnfe(ji,jj,jk) / ( xlimnfe(ji,jj,jk) + 0.5 ) )    & 
    267                   &             * biron(ji,jj,jk) / ( biron(ji,jj,jk) + concnfe(ji,jj,jk) )  & 
    268                   &            * zmax * trb(ji,jj,jk,jpphy) * rfact2 
     280                  &          * (1. + 0.8 * xnanono3(ji,jj,jk) / ( rtrn + xnanono3(ji,jj,jk)  & 
     281                  &          + xnanonh4(ji,jj,jk) ) * (1. - xnanofer(ji,jj,jk) ) )   & 
     282                  &          * xnanofer(ji,jj,jk) * zmax * trb(ji,jj,jk,jpphy) * rfact2 
    269283                  !  production terms of diatoms (C) 
    270284                  zprorcad(ji,jj,jk) = zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * trb(ji,jj,jk,jpdia) * rfact2 
     
    272286                  ! New production (uptake of NO3) 
    273287                  zpronewd(ji,jj,jk) = zprorcad(ji,jj,jk) * xdiatno3(ji,jj,jk) / ( xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) + rtrn ) 
     288 
     289                  ! Size computation 
     290                  ! Size is made a function of the limitation of of phytoplankton growth 
     291                  ! Strongly limited cells are supposed to be smaller. sizeda is 
     292                  ! size at time step t+1 and is thus updated at the end of the  
     293                  ! current time step.  
     294                  ! -------------------------------------------------------------------- 
     295                  zlimfac = zprchld(ji,jj,jk) * xlimdia(ji,jj,jk) / ( zprmaxd(ji,jj,jk) + rtrn ) 
     296                  zsizetmp = 1.0 + 1.3 * ( xsizerd - 1.0 ) * zlimfac**3/(0.3 + zlimfac**3) 
     297                  sizeda(ji,jj,jk) = min(xsizerd, max( sizeda(ji,jj,jk), zsizetmp ) ) 
    274298 
    275299                  !  Iron uptake rates of nanophytoplankton. Upregulation   
     
    277301                  !  formulation used in quota formulations. Uptake is downregulated 
    278302                  !  when the quota is close to the maximum quota  
    279                   zratio = trb(ji,jj,jk,jpdfe) / ( trb(ji,jj,jk,jpdia) * fecdm + rtrn ) 
    280                   zmax   = MAX( 0., ( 1. - zratio ) / ABS( 1.05 - zratio ) )  
    281                   zprofed(ji,jj,jk) = fecdm * zprmaxd(ji,jj,jk) * ( 1.0 - fr_i(ji,jj) )  & 
    282                   &             * ( 4. - 4.5 * xlimdfe(ji,jj,jk) / ( xlimdfe(ji,jj,jk) + 0.5 ) )    & 
    283                   &             * biron(ji,jj,jk) / ( biron(ji,jj,jk) + concdfe(ji,jj,jk) )  & 
    284                   &            * zmax * trb(ji,jj,jk,jpdia) * rfact2 
     303                  zratio = 1.0 - MIN(1.0, trb(ji,jj,jk,jpdfe) / ( trb(ji,jj,jk,jpdia) * fecdm + rtrn ) ) 
     304                  zmax   = MAX( 0., MIN( 1.0, zratio**2/ (0.05**2+zratio**2) ) )  
     305                  zprofed(ji,jj,jk) = fecdm * zprmaxd(ji,jj,jk) * (1.0 - fr_i(ji,jj) )  & 
     306                  &          * (1. + 0.8 * xdiatno3(ji,jj,jk) / ( rtrn + xdiatno3(ji,jj,jk)  & 
     307                  &          + xdiatnh4(ji,jj,jk) ) * (1. - xdiatfer(ji,jj,jk) ) )   & 
     308                  &          * xdiatfer(ji,jj,jk) * zmax * trb(ji,jj,jk,jpdia) * rfact2 
    285309               ENDIF 
    286310            END DO 
     
    336360                 tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) + zprorcad(ji,jj,jk) * texcretd 
    337361                 tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) + zprofed(ji,jj,jk) * texcretd 
    338                  tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) + zprorcad(ji,jj,jk) * zysopt(ji,jj,jk) * texcretd 
     362                 tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) + zprmaxd(ji,jj,jk) * zysopt(ji,jj,jk)   & 
     363                 &                   * rfact2 * trb(ji,jj,jk,jpdia) 
    339364                 tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zdocprod 
    340365                 tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + o2ut * ( zproreg + zproreg2) & 
     
    343368                 zfeup = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk) 
    344369                 tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zfeup 
    345                  tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) - texcretd * zprorcad(ji,jj,jk) * zysopt(ji,jj,jk) 
     370                 tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) - zprmaxd(ji,jj,jk) * zysopt(ji,jj,jk)   & 
     371                 &                   * rfact2 * trb(ji,jj,jk,jpdia) 
    346372                 tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprorcan(ji,jj,jk) - zprorcad(ji,jj,jk) 
    347373                 tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * ( zpronewn(ji,jj,jk) + zpronewd(ji,jj,jk) ) & 
     
    401427          ENDIF 
    402428          IF( iom_use( "PBSi" ) )  THEN 
    403               zw3d(:,:,:) = zprorcad(:,:,:) * zfact * tmask(:,:,:) * zysopt(:,:,:) ! biogenic silica production 
     429              zw3d(:,:,:) = zprmaxd(:,:,:) * 1.E3 * tmask(:,:,:) * zysopt(:,:,:) * trb(:,:,:,jpdia) ! biogenic silica production 
    404430              CALL iom_put( "PBSi"  , zw3d ) 
    405431          ENDIF 
Note: See TracChangeset for help on using the changeset viewer.