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 5965 for branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90 – NEMO

Ignore:
Timestamp:
2015-12-01T16:35:30+01:00 (8 years ago)
Author:
timgraham
Message:

Upgraded branch to r5518 of trunk (v3.6 stable revision)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90

    • Property svn:keywords set to Id
    r4641 r5965  
    2121   USE p4zopt          !  optical model 
    2222   USE p4zlim          !  Co-limitations of differents nutrients 
    23    USE p4zrem          !  Remineralisation of organic matter 
    2423   USE p4zsbc          !  External source of nutrients  
    2524   USE p4zint          !  interpolation and computation of various fields 
     
    3029   PRIVATE 
    3130 
    32    PUBLIC   p4z_sed    
     31   PUBLIC   p4z_sed   
     32   PUBLIC   p4z_sed_alloc 
     33  
    3334 
    3435   !! * Module variables 
    35    REAL(wp) :: ryyss                    !: number of seconds per year  
    36    REAL(wp) :: r1_ryyss                 !: inverse of ryyss 
    37    REAL(wp) :: rmtss                    !: number of seconds per month 
     36   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: nitrpot    !: Nitrogen fixation  
     37   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:  ) :: sdenit     !: Nitrate reduction in the sediments 
    3838   REAL(wp) :: r1_rday                  !: inverse of rday 
    39  
    40    INTEGER ::  numnit   
    41  
    4239 
    4340   !!* Substitution 
     
    4542   !!---------------------------------------------------------------------- 
    4643   !! NEMO/TOP 3.3 , NEMO Consortium (2010) 
    47    !! $Header:$  
     44   !! $Id$  
    4845   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    4946   !!---------------------------------------------------------------------- 
    5047CONTAINS 
    5148 
    52    SUBROUTINE p4z_sed( kt, jnt ) 
     49   SUBROUTINE p4z_sed( kt, knt ) 
    5350      !!--------------------------------------------------------------------- 
    5451      !!                     ***  ROUTINE p4z_sed  *** 
     
    6158      !!--------------------------------------------------------------------- 
    6259      ! 
    63       INTEGER, INTENT(in) ::   kt, jnt ! ocean time step 
     60      INTEGER, INTENT(in) ::   kt, knt ! ocean time step 
    6461      INTEGER  ::   ji, jj, jk, ikt 
    6562#if ! defined key_sed 
     
    7269      REAL(wp) ::  zsiloss, zcaloss, zws3, zws4, zwsc, zdep, zwstpoc 
    7370      REAL(wp) ::  ztrfer, ztrpo4, zwdust, zlight 
    74       REAL(wp) ::  zrdenittot, zsdenittot, znitrpottot 
    7571      ! 
    7672      CHARACTER (len=25) :: charout 
    77       REAL(wp), POINTER, DIMENSION(:,:  ) :: zpdep, zsidep, zwork1, zwork2, zwork3, zwork4 
     73      REAL(wp), POINTER, DIMENSION(:,:  ) :: zpdep, zsidep, zwork1, zwork2, zwork3 
    7874      REAL(wp), POINTER, DIMENSION(:,:  ) :: zdenit2d, zironice, zbureff 
    7975      REAL(wp), POINTER, DIMENSION(:,:  ) :: zwsbio3, zwsbio4, zwscal 
    80       REAL(wp), POINTER, DIMENSION(:,:,:) :: znitrpot, zirondep, zsoufer 
     76      REAL(wp), POINTER, DIMENSION(:,:,:) :: zirondep, zsoufer 
    8177      !!--------------------------------------------------------------------- 
    8278      ! 
    8379      IF( nn_timing == 1 )  CALL timing_start('p4z_sed') 
    8480      ! 
    85       IF( kt == nittrc000 .AND. jnt == 1 )  THEN 
    86          ryyss    = nyear_len(1) * rday    ! number of seconds per year and per month 
    87          rmtss    = ryyss / raamo 
    88          r1_rday  = 1. / rday 
    89          r1_ryyss = 1. / ryyss 
    90          IF( ln_check_mass .AND. lwp)  & 
    91            &  CALL ctl_opn( numnit, 'nitrogen.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea ) 
    92       ENDIF 
     81      IF( kt == nittrc000 .AND. knt == 1 )   r1_rday  = 1. / rday 
    9382      ! 
    9483      ! Allocate temporary workspace 
    95       CALL wrk_alloc( jpi, jpj, zdenit2d, zwork1, zwork2, zwork3, zwork4, zbureff ) 
     84      CALL wrk_alloc( jpi, jpj, zdenit2d, zwork1, zwork2, zwork3, zbureff ) 
    9685      CALL wrk_alloc( jpi, jpj, zwsbio3, zwsbio4, zwscal ) 
    97       CALL wrk_alloc( jpi, jpj, jpk, znitrpot, zsoufer ) 
     86      CALL wrk_alloc( jpi, jpj, jpk, zsoufer ) 
    9887 
    9988      zdenit2d(:,:) = 0.e0 
    10089      zbureff (:,:) = 0.e0 
     90      zwork1  (:,:) = 0.e0 
     91      zwork2  (:,:) = 0.e0 
     92      zwork3  (:,:) = 0.e0 
    10193 
    10294      ! Iron input/uptake due to sea ice : Crude parameterization based on Lancelot et al. 
     
    110102               zdep    = rfact2 / fse3t(ji,jj,1) 
    111103               zwflux  = fmmflx(ji,jj) / 1000._wp 
    112                zfminus = MIN( 0._wp, -zwflux ) * trn(ji,jj,1,jpfer) * zdep 
     104               zfminus = MIN( 0._wp, -zwflux ) * trb(ji,jj,1,jpfer) * zdep 
    113105               zfplus  = MAX( 0._wp, -zwflux ) * icefeinput * zdep 
    114106               zironice(ji,jj) =  zfplus + zfminus 
     
    116108         END DO 
    117109         ! 
    118          trn(:,:,1,jpfer) = trn(:,:,1,jpfer) + zironice(:,:)  
    119          !                                               
    120          IF( ln_diatrc .AND. lk_iomput .AND. jnt == nrdttrc )   & 
     110         tra(:,:,1,jpfer) = tra(:,:,1,jpfer) + zironice(:,:)  
     111         !  
     112         IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "Ironice" ) )   & 
    121113            &   CALL iom_put( "Ironice", zironice(:,:) * 1.e+3 * rfact2r * fse3t(:,:,1) * tmask(:,:,1) ) ! iron flux from ice 
     114         ! 
    122115         CALL wrk_dealloc( jpi, jpj, zironice ) 
    123116         !                                               
     
    132125         !                                              ! Iron and Si deposition at the surface 
    133126         IF( ln_solub ) THEN 
    134             zirondep(:,:,1) = solub(:,:) * dust(:,:) * mfrac * rfact2 / fse3t(:,:,1) / ( 55.85 * rmtss ) + 3.e-10 * r1_ryyss  
     127            zirondep(:,:,1) = solub(:,:) * dust(:,:) * mfrac * rfact2 / fse3t(:,:,1) / 55.85 + 3.e-10 * r1_ryyss  
    135128         ELSE 
    136             zirondep(:,:,1) = dustsolub  * dust(:,:) * mfrac * rfact2 / fse3t(:,:,1) / ( 55.85 * rmtss ) + 3.e-10 * r1_ryyss  
     129            zirondep(:,:,1) = dustsolub  * dust(:,:) * mfrac * rfact2 / fse3t(:,:,1) / 55.85 + 3.e-10 * r1_ryyss  
    137130         ENDIF 
    138          zsidep(:,:) = 8.8 * 0.075 * dust(:,:) * mfrac * rfact2 / fse3t(:,:,1) / ( 28.1  * rmtss ) 
    139          zpdep (:,:) = 0.1 * 0.021 * dust(:,:) * mfrac * rfact2 / fse3t(:,:,1) / ( 31.   * rmtss ) / po4r  
     131         zsidep(:,:) = 8.8 * 0.075 * dust(:,:) * mfrac * rfact2 / fse3t(:,:,1) / 28.1  
     132         zpdep (:,:) = 0.1 * 0.021 * dust(:,:) * mfrac * rfact2 / fse3t(:,:,1) / 31. / po4r  
    140133         !                                              ! Iron solubilization of particles in the water column 
    141134         !                                              ! dust in kg/m2/s ---> 1/55.85 to put in mol/Fe ;  wdust in m/j 
     
    145138         END DO 
    146139         !                                              ! Iron solubilization of particles in the water column 
    147          trn(:,:,1,jppo4) = trn(:,:,1,jppo4) + zpdep   (:,:) 
    148          trn(:,:,1,jpsil) = trn(:,:,1,jpsil) + zsidep  (:,:) 
    149          trn(:,:,:,jpfer) = trn(:,:,:,jpfer) + zirondep(:,:,:)  
    150          !                                               
    151          IF( ln_diatrc ) THEN 
    152             zfact = 1.e+3 * rfact2r 
    153             IF( lk_iomput ) THEN 
    154                IF( jnt == nrdttrc ) THEN 
    155                   CALL iom_put( "Irondep", zirondep(:,:,1) * zfact * fse3t(:,:,1) * tmask(:,:,1) ) ! surface downward dust depo of iron 
    156                   CALL iom_put( "pdust"  , dust(:,:) / ( wdust * rday )  * tmask(:,:,1) ) ! dust concentration at surface 
    157                ENDIF 
    158             ELSE 
    159                trc2d(:,:,jp_pcs0_2d + 11) = zirondep(:,:,1) * zfact * fse3t(:,:,1) * tmask(:,:,1) 
     140         tra(:,:,1,jppo4) = tra(:,:,1,jppo4) + zpdep   (:,:) 
     141         tra(:,:,1,jpsil) = tra(:,:,1,jpsil) + zsidep  (:,:) 
     142         tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + zirondep(:,:,:)  
     143         !  
     144         IF( lk_iomput ) THEN 
     145            IF( knt == nrdttrc ) THEN 
     146                IF( iom_use( "Irondep" ) )   & 
     147                &  CALL iom_put( "Irondep", zirondep(:,:,1) * 1.e+3 * rfact2r * fse3t(:,:,1) * tmask(:,:,1) ) ! surface downward dust depo of iron 
     148                IF( iom_use( "pdust" ) )   & 
     149                &  CALL iom_put( "pdust"  , dust(:,:) / ( wdust * rday )  * tmask(:,:,1) ) ! dust concentration at surface 
    160150            ENDIF 
     151         ELSE                                     
     152            IF( ln_diatrc )  & 
     153              &  trc2d(:,:,jp_pcs0_2d + 11) = zirondep(:,:,1) * 1.e+3 * rfact2r * fse3t(:,:,1) * tmask(:,:,1) 
    161154         ENDIF 
    162155         CALL wrk_dealloc( jpi, jpj,      zpdep, zsidep ) 
     
    168161      ! ---------------------------------------------------------- 
    169162      IF( ln_river ) THEN 
    170          trn(:,:,1,jppo4) = trn(:,:,1,jppo4) + rivdip(:,:) * rfact2 
    171          trn(:,:,1,jpno3) = trn(:,:,1,jpno3) + rivdin(:,:) * rfact2 
    172          trn(:,:,1,jpfer) = trn(:,:,1,jpfer) + rivdic(:,:) * 5.e-5 * rfact2 
    173          trn(:,:,1,jpsil) = trn(:,:,1,jpsil) + rivdsi(:,:) * rfact2 
    174          trn(:,:,1,jpdic) = trn(:,:,1,jpdic) + rivdic(:,:) * rfact2 
    175          trn(:,:,1,jptal) = trn(:,:,1,jptal) + ( rivalk(:,:) - rno3 * rivdin(:,:) ) * rfact2 
     163         DO jj = 1, jpj 
     164            DO ji = 1, jpi 
     165               DO jk = 1, nk_rnf(ji,jj) 
     166                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) +  rivdip(ji,jj) * rfact2 
     167                  tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) +  rivdin(ji,jj) * rfact2 
     168                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) +  rivdic(ji,jj) * 5.e-5 * rfact2 
     169                  tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) +  rivdsi(ji,jj) * rfact2 
     170                  tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) +  rivdic(ji,jj) * rfact2 
     171                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) +  ( rivalk(ji,jj) - rno3 * rivdin(ji,jj) ) * rfact2 
     172               ENDDO 
     173            ENDDO 
     174         ENDDO 
    176175      ENDIF 
    177176       
     
    179178      ! ---------------------------------------------------------- 
    180179      IF( ln_ndepo ) THEN 
    181          trn(:,:,1,jpno3) = trn(:,:,1,jpno3) + nitdep(:,:) * rfact2 
    182          trn(:,:,1,jptal) = trn(:,:,1,jptal) - rno3 * nitdep(:,:) * rfact2 
     180         tra(:,:,1,jpno3) = tra(:,:,1,jpno3) + nitdep(:,:) * rfact2 
     181         tra(:,:,1,jptal) = tra(:,:,1,jptal) - rno3 * nitdep(:,:) * rfact2 
    183182      ENDIF 
    184183 
     
    186185      ! ------------------------------------------------------ 
    187186      IF( ln_ironsed ) THEN 
    188          trn(:,:,:,jpfer) = trn(:,:,:,jpfer) + ironsed(:,:,:) * rfact2 
     187         tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + ironsed(:,:,:) * rfact2 
    189188         ! 
    190          IF( ln_diatrc .AND. lk_iomput .AND. jnt == nrdttrc )   & 
     189         IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "Ironsed" ) )   & 
    191190            &   CALL iom_put( "Ironsed", ironsed(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! iron inputs from sediments 
    192191      ENDIF 
     
    195194      ! ------------------------------------------------------ 
    196195      IF( ln_hydrofe ) THEN 
    197          trn(:,:,:,jpfer) = trn(:,:,:,jpfer) + hydrofe(:,:,:) * rfact2 
     196         tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + hydrofe(:,:,:) * rfact2 
    198197         ! 
    199          IF( ln_diatrc .AND. lk_iomput .AND. jnt == nrdttrc )   & 
     198         IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "HYDR" ) )   & 
    200199            &   CALL iom_put( "HYDR", hydrofe(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! hydrothermal iron input 
    201200      ENDIF 
    202  
    203201 
    204202      ! OA: Warning, the following part is necessary, especially with Kriest 
     
    224222              ikt = mbkt(ji,jj) 
    225223# if defined key_kriest 
    226               zflx =    trn(ji,jj,ikt,jppoc) * zwsbio3(ji,jj)    * 1E3 * 1E6 / 1E4 
     224              zflx =    trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj)    * 1E3 * 1E6 / 1E4 
    227225# else 
    228               zflx = (  trn(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   & 
    229                 &     + trn(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) )  * 1E3 * 1E6 / 1E4 
     226              zflx = (  trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   & 
     227                &     + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) )  * 1E3 * 1E6 / 1E4 
    230228#endif 
    231229              zflx  = LOG10( MAX( 1E-3, zflx ) ) 
    232               zo2   = LOG10( MAX( 10. , trn(ji,jj,ikt,jpoxy) * 1E6 ) ) 
    233               zno3  = LOG10( MAX( 1.  , trn(ji,jj,ikt,jpno3) * 1E6 * rno3 ) ) 
     230              zo2   = LOG10( MAX( 10. , trb(ji,jj,ikt,jpoxy) * 1E6 ) ) 
     231              zno3  = LOG10( MAX( 1.  , trb(ji,jj,ikt,jpno3) * 1E6 * rno3 ) ) 
    234232              zdep  = LOG10( fsdepw(ji,jj,ikt+1) ) 
    235233              zdenit2d(ji,jj) = -2.2567 - 1.185 * zflx - 0.221 * zflx**2 - 0.3995 * zno3 * zo2 + 1.25 * zno3    & 
     
    237235              zdenit2d(ji,jj) = 10.0**( zdenit2d(ji,jj) ) 
    238236              ! 
    239               zflx = (  trn(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   & 
    240                 &     + trn(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) ) * 1E6 
     237              zflx = (  trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   & 
     238                &     + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) ) * 1E6 
    241239              zbureff(ji,jj) = 0.013 + 0.53 * zflx**2 / ( 7.0 + zflx )**2 
    242240           ENDIF 
     
    250248      DO jj = 1, jpj 
    251249         DO ji = 1, jpi 
    252             ikt = mbkt(ji,jj)  
     250            IF( tmask(ji,jj,1) == 1 ) THEN 
     251               ikt = mbkt(ji,jj)  
    253252# if defined key_kriest 
    254             zwork1(ji,jj) = trn(ji,jj,ikt,jpgsi) * zwscal (ji,jj) 
    255             zwork2(ji,jj) = trn(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) 
     253               zwork1(ji,jj) = trb(ji,jj,ikt,jpgsi) * zwscal (ji,jj) 
     254               zwork2(ji,jj) = trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) 
    256255# else 
    257             zwork1(ji,jj) = trn(ji,jj,ikt,jpgsi) * zwsbio4(ji,jj) 
    258             zwork2(ji,jj) = trn(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj) + trn(ji,jj,ikt,jppoc) * zwsbio3(ji,jj)  
     256               zwork1(ji,jj) = trb(ji,jj,ikt,jpgsi) * zwsbio4(ji,jj) 
     257               zwork2(ji,jj) = trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj) + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj)  
    259258# endif 
    260             ! For calcite, burial efficiency is made a function of saturation 
    261             zfactcal      = MIN( excess(ji,jj,ikt), 0.2 ) 
    262             zfactcal      = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) ) 
    263             zwork3(ji,jj) = trn(ji,jj,ikt,jpcal) * zwscal(ji,jj) * 2.e0 * zfactcal 
     259               ! For calcite, burial efficiency is made a function of saturation 
     260               zfactcal      = MIN( excess(ji,jj,ikt), 0.2 ) 
     261               zfactcal      = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) ) 
     262               zwork3(ji,jj) = trb(ji,jj,ikt,jpcal) * zwscal(ji,jj) * 2.e0 * zfactcal 
     263            ENDIF 
    264264         END DO 
    265265      END DO 
     
    279279         DO ji = 1, jpi 
    280280            ikt  = mbkt(ji,jj) 
    281             zdep = xstep / fse3t(ji,jj,ikt) 
     281            zdep = xstep / fse3t(ji,jj,ikt)  
    282282            zws4 = zwsbio4(ji,jj) * zdep 
    283283            zwsc = zwscal (ji,jj) * zdep 
    284284# if defined key_kriest 
    285             zsiloss = trn(ji,jj,ikt,jpgsi) * zws4 
     285            zsiloss = trb(ji,jj,ikt,jpgsi) * zws4 
    286286# else 
    287             zsiloss = trn(ji,jj,ikt,jpgsi) * zwsc 
     287            zsiloss = trb(ji,jj,ikt,jpgsi) * zwsc 
    288288# endif 
    289             zcaloss = trn(ji,jj,ikt,jpcal) * zwsc 
     289            zcaloss = trb(ji,jj,ikt,jpcal) * zwsc 
    290290            ! 
    291             trn(ji,jj,ikt,jpgsi) = trn(ji,jj,ikt,jpgsi) - zsiloss 
    292             trn(ji,jj,ikt,jpcal) = trn(ji,jj,ikt,jpcal) - zcaloss 
     291            tra(ji,jj,ikt,jpgsi) = tra(ji,jj,ikt,jpgsi) - zsiloss 
     292            tra(ji,jj,ikt,jpcal) = tra(ji,jj,ikt,jpcal) - zcaloss 
    293293#if ! defined key_sed 
    294             trn(ji,jj,ikt,jpsil) = trn(ji,jj,ikt,jpsil) + zsiloss * zrivsil  
     294            tra(ji,jj,ikt,jpsil) = tra(ji,jj,ikt,jpsil) + zsiloss * zrivsil  
    295295            zfactcal = MIN( excess(ji,jj,ikt), 0.2 ) 
    296296            zfactcal = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) ) 
    297297            zrivalk  =  1._wp - ( rivalkinput * r1_ryyss ) * zfactcal / ( zsumsedcal + rtrn ) 
    298             trn(ji,jj,ikt,jptal) =  trn(ji,jj,ikt,jptal) + zcaloss * zrivalk * 2.0 
    299             trn(ji,jj,ikt,jpdic) =  trn(ji,jj,ikt,jpdic) + zcaloss * zrivalk 
     298            tra(ji,jj,ikt,jptal) =  tra(ji,jj,ikt,jptal) + zcaloss * zrivalk * 2.0 
     299            tra(ji,jj,ikt,jpdic) =  tra(ji,jj,ikt,jpdic) + zcaloss * zrivalk 
    300300#endif 
    301301         END DO 
     
    304304      DO jj = 1, jpj 
    305305         DO ji = 1, jpi 
    306             ikt     = mbkt(ji,jj) 
    307             zdep    = xstep / fse3t(ji,jj,ikt) 
     306            ikt  = mbkt(ji,jj) 
     307            zdep = xstep / fse3t(ji,jj,ikt)  
    308308            zws4 = zwsbio4(ji,jj) * zdep 
    309309            zws3 = zwsbio3(ji,jj) * zdep 
    310310            zrivno3 = 1. - zbureff(ji,jj) 
    311311# if ! defined key_kriest 
    312             trn(ji,jj,ikt,jpgoc) = trn(ji,jj,ikt,jpgoc) - trn(ji,jj,ikt,jpgoc) * zws4 
    313             trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc) - trn(ji,jj,ikt,jppoc) * zws3 
    314             trn(ji,jj,ikt,jpbfe) = trn(ji,jj,ikt,jpbfe) - trn(ji,jj,ikt,jpbfe) * zws4 
    315             trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe) - trn(ji,jj,ikt,jpsfe) * zws3 
    316             zwstpoc              =  trn(ji,jj,ikt,jpgoc) * zws4 + trn(ji,jj,ikt,jppoc) * zws3  
     312            tra(ji,jj,ikt,jpgoc) = tra(ji,jj,ikt,jpgoc) - trb(ji,jj,ikt,jpgoc) * zws4  
     313            tra(ji,jj,ikt,jppoc) = tra(ji,jj,ikt,jppoc) - trb(ji,jj,ikt,jppoc) * zws3 
     314            tra(ji,jj,ikt,jpbfe) = tra(ji,jj,ikt,jpbfe) - trb(ji,jj,ikt,jpbfe) * zws4 
     315            tra(ji,jj,ikt,jpsfe) = tra(ji,jj,ikt,jpsfe) - trb(ji,jj,ikt,jpsfe) * zws3 
     316            zwstpoc              = trb(ji,jj,ikt,jpgoc) * zws4 + trb(ji,jj,ikt,jppoc) * zws3 
    317317# else 
    318             trn(ji,jj,ikt,jpnum) = trn(ji,jj,ikt,jpnum) - trn(ji,jj,ikt,jpnum) * zws4 
    319             trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc) - trn(ji,jj,ikt,jppoc) * zws3 
    320             trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe) - trn(ji,jj,ikt,jpsfe) * zws3 
    321             zwstpoc = trn(ji,jj,ikt,jppoc) * zws3  
     318            tra(ji,jj,ikt,jpnum) = tra(ji,jj,ikt,jpnum) - trb(ji,jj,ikt,jpnum) * zws4  
     319            tra(ji,jj,ikt,jppoc) = tra(ji,jj,ikt,jppoc) - trb(ji,jj,ikt,jppoc) * zws3 
     320            tra(ji,jj,ikt,jpsfe) = tra(ji,jj,ikt,jpsfe) - trb(ji,jj,ikt,jpsfe) * zws3 
     321            zwstpoc = trb(ji,jj,ikt,jppoc) * zws3  
    322322# endif 
    323323 
     
    325325            ! The 0.5 factor in zpdenit and zdenitt is to avoid negative NO3 concentration after both denitrification 
    326326            ! in the sediments and just above the sediments. Not very clever, but simpliest option. 
    327             zpdenit  = MIN( 0.5 * ( trn(ji,jj,ikt,jpno3) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3 ) 
     327            zpdenit  = MIN( 0.5 * ( trb(ji,jj,ikt,jpno3) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3 ) 
    328328            z1pdenit = zwstpoc * zrivno3 - zpdenit 
    329             zolimit = MIN( ( trn(ji,jj,ikt,jpoxy) - rtrn ) / o2ut, z1pdenit * ( 1.- nitrfac(ji,jj,ikt) ) ) 
    330             zdenitt = MIN(  0.5 * ( trn(ji,jj,ikt,jpno3) - rtrn ) / rdenit, z1pdenit * nitrfac(ji,jj,ikt) ) 
    331             trn(ji,jj,ikt,jpdoc) = trn(ji,jj,ikt,jpdoc) + z1pdenit - zolimit - zdenitt 
    332             trn(ji,jj,ikt,jppo4) = trn(ji,jj,ikt,jppo4) + zpdenit + zolimit + zdenitt 
    333             trn(ji,jj,ikt,jpnh4) = trn(ji,jj,ikt,jpnh4) + zpdenit + zolimit + zdenitt 
    334             trn(ji,jj,ikt,jpno3) = trn(ji,jj,ikt,jpno3) - rdenit * (zpdenit + zdenitt) 
    335             trn(ji,jj,ikt,jpoxy) = trn(ji,jj,ikt,jpoxy) - zolimit * o2ut 
    336             trn(ji,jj,ikt,jptal) = trn(ji,jj,ikt,jptal) + rno3 * (zolimit + (1.+rdenit) * (zpdenit + zdenitt) ) 
    337             trn(ji,jj,ikt,jpdic) = trn(ji,jj,ikt,jpdic) + zpdenit + zolimit + zdenitt 
    338             zwork4(ji,jj) = rdenit * zpdenit * fse3t(ji,jj,ikt) 
     329            zolimit = MIN( ( trb(ji,jj,ikt,jpoxy) - rtrn ) / o2ut, z1pdenit * ( 1.- nitrfac(ji,jj,ikt) ) ) 
     330            zdenitt = MIN(  0.5 * ( trb(ji,jj,ikt,jpno3) - rtrn ) / rdenit, z1pdenit * nitrfac(ji,jj,ikt) ) 
     331            tra(ji,jj,ikt,jpdoc) = tra(ji,jj,ikt,jpdoc) + z1pdenit - zolimit - zdenitt 
     332            tra(ji,jj,ikt,jppo4) = tra(ji,jj,ikt,jppo4) + zpdenit + zolimit + zdenitt 
     333            tra(ji,jj,ikt,jpnh4) = tra(ji,jj,ikt,jpnh4) + zpdenit + zolimit + zdenitt 
     334            tra(ji,jj,ikt,jpno3) = tra(ji,jj,ikt,jpno3) - rdenit * (zpdenit + zdenitt) 
     335            tra(ji,jj,ikt,jpoxy) = tra(ji,jj,ikt,jpoxy) - zolimit * o2ut 
     336            tra(ji,jj,ikt,jptal) = tra(ji,jj,ikt,jptal) + rno3 * (zolimit + (1.+rdenit) * (zpdenit + zdenitt) ) 
     337            tra(ji,jj,ikt,jpdic) = tra(ji,jj,ikt,jpdic) + zpdenit + zolimit + zdenitt 
     338            sdenit(ji,jj) = rdenit * zpdenit * fse3t(ji,jj,ikt) 
    339339#endif 
    340340         END DO 
     
    356356#endif 
    357357               ztrfer = biron(ji,jj,jk)       / ( concfediaz + biron(ji,jj,jk)       ) 
    358                ztrpo4 = trn  (ji,jj,jk,jppo4) / ( concnnh4   + trn  (ji,jj,jk,jppo4) )  
    359                zlight =  ( 1.- EXP( -etot(ji,jj,jk) / diazolight ) )  
    360                znitrpot(ji,jj,jk) =  MAX( 0.e0, ( 0.6 * tgfunc(ji,jj,jk) - 2.15 ) * r1_rday )   & 
     358               ztrpo4 = trb  (ji,jj,jk,jppo4) / ( concnnh4   + trb  (ji,jj,jk,jppo4) )  
     359               zlight =  ( 1.- EXP( -etot_ndcy(ji,jj,jk) / diazolight ) )  
     360               nitrpot(ji,jj,jk) =  MAX( 0.e0, ( 0.6 * tgfunc(ji,jj,jk) - 2.15 ) * r1_rday )   & 
    361361                 &         *  zfact * MIN( ztrfer, ztrpo4 ) * zlight 
    362362               zsoufer(ji,jj,jk) = zlight * 2E-11 / (2E-11 + biron(ji,jj,jk)) 
     
    370370         DO jj = 1, jpj 
    371371            DO ji = 1, jpi 
    372                zfact = znitrpot(ji,jj,jk) * nitrfix 
    373                trn(ji,jj,jk,jpnh4) = trn(ji,jj,jk,jpnh4) +             zfact 
    374                trn(ji,jj,jk,jptal) = trn(ji,jj,jk,jptal) + rno3      * zfact 
    375                trn(ji,jj,jk,jpoxy) = trn(ji,jj,jk,jpoxy) + o2nit     * zfact  
    376                trn(ji,jj,jk,jppo4) = trn(ji,jj,jk,jppo4) + concdnh4 / ( concdnh4 + trn(ji,jj,jk,jppo4) ) & 
    377                &                     * 0.002 * trn(ji,jj,jk,jpdoc) * rfact2 / rday 
    378                trn(ji,jj,jk,jpfer) = trn(ji,jj,jk,jpfer) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday 
     372               zfact = nitrpot(ji,jj,jk) * nitrfix 
     373               tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) +             zfact 
     374               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3      * zfact 
     375               tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + o2nit     * zfact  
     376               tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + concdnh4 / ( concdnh4 + trb(ji,jj,jk,jppo4) ) & 
     377               &                     * 0.002 * trb(ji,jj,jk,jpdoc) * xstep 
     378               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * xstep 
    379379           END DO 
    380380         END DO  
    381381      END DO 
    382382 
    383   
    384       IF( ln_check_mass ) THEN 
    385         ! Global budget of N SMS : denitrification in the water column and in the sediment 
    386          !                          nitrogen fixation by the diazotrophs 
    387          ! -------------------------------------------------------------------------------- 
    388          zrdenittot   = glob_sum ( denitr(:,:,:) * rdenit * xnegtr(:,:,:) * cvol(:,:,:) ) 
    389          zsdenittot   = glob_sum ( zwork4(:,:)   * e1e2t(:,:) ) 
    390          znitrpottot  = glob_sum ( znitrpot(:,:,:) * nitrfix              * cvol(:,:,:) ) 
    391          IF( kt == nitend .AND. jnt == nrdttrc ) THEN 
    392             zfact = 1.e+3 * rfact2r * rno3 * ryyss * 14. / 1e12 
    393             IF(lwp) WRITE(numnit,9100) ndastp, znitrpottot * nitrfix * zfact, zrdenittot * zfact , zsdenittot * zfact 
     383      IF( lk_iomput ) THEN 
     384         IF( knt == nrdttrc ) THEN 
     385            zfact = 1.e+3 * rfact2r * rno3  !  conversion from molC/l/kt  to molN/m3/s 
     386            IF( iom_use("Nfix"   ) ) CALL iom_put( "Nfix", nitrpot(:,:,:) * nitrfix * zfact * tmask(:,:,:) )  ! nitrogen fixation  
     387            IF( iom_use("INTNFIX") ) THEN   ! nitrogen fixation rate in ocean ( vertically integrated ) 
     388               zwork1(:,:) = 0. 
     389               DO jk = 1, jpkm1 
     390                 zwork1(:,:) = zwork1(:,:) + nitrpot(:,:,jk) * nitrfix * zfact * fse3t(:,:,jk) * tmask(:,:,jk) 
     391               ENDDO 
     392               CALL iom_put( "INTNFIX" , zwork1 )  
     393            ENDIF 
    394394         ENDIF 
    395        ENDIF 
    396       ! 
    397       IF( ln_diatrc ) THEN 
    398          zfact = 1.e+3 * rfact2r 
    399          IF( lk_iomput ) THEN 
    400             IF( jnt == nrdttrc ) THEN 
    401                CALL iom_put( "Nfix"  , znitrpot(:,:,:) * nitrfix * rno3 * zfact * tmask(:,:,:) )  ! nitrogen fixation  
    402                CALL iom_put( "Sdenit", zwork4(:,:)               * rno3 * zfact * tmask(:,:,1) )  ! Nitrate reduction in the sediments 
    403             ENDIF 
    404          ELSE 
    405             trc2d(:,:,jp_pcs0_2d + 12) = znitrpot(:,:,1) * nitrfix * zfact * fse3t(:,:,1) * tmask(:,:,1) 
    406          ENDIF 
     395      ELSE 
     396         IF( ln_diatrc )  & 
     397            &  trc2d(:,:,jp_pcs0_2d + 12) = nitrpot(:,:,1) * nitrfix * rno3 * 1.e+3 * rfact2r * fse3t(:,:,1) * tmask(:,:,1) 
    407398      ENDIF 
    408399      ! 
     
    410401         WRITE(charout, fmt="('sed ')") 
    411402         CALL prt_ctl_trc_info(charout) 
    412          CALL prt_ctl_trc(tab4d=trn, mask=tmask, clinfo=ctrcnm) 
    413       ENDIF 
    414       ! 
    415       CALL wrk_dealloc( jpi, jpj, zdenit2d, zwork1, zwork2, zwork3, zwork4, zbureff ) 
     403         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm) 
     404      ENDIF 
     405      ! 
     406      CALL wrk_dealloc( jpi, jpj, zdenit2d, zwork1, zwork2, zwork3, zbureff ) 
    416407      CALL wrk_dealloc( jpi, jpj, zwsbio3, zwsbio4, zwscal ) 
    417       CALL wrk_dealloc( jpi, jpj, jpk, znitrpot, zsoufer ) 
     408      CALL wrk_dealloc( jpi, jpj, jpk, zsoufer ) 
    418409      ! 
    419410      IF( nn_timing == 1 )  CALL timing_stop('p4z_sed') 
     
    422413      ! 
    423414   END SUBROUTINE p4z_sed 
     415 
     416 
     417   INTEGER FUNCTION p4z_sed_alloc() 
     418      !!---------------------------------------------------------------------- 
     419      !!                     ***  ROUTINE p4z_sed_alloc  *** 
     420      !!---------------------------------------------------------------------- 
     421      ALLOCATE( nitrpot(jpi,jpj,jpk), sdenit(jpi,jpj), STAT=p4z_sed_alloc ) 
     422      ! 
     423      IF( p4z_sed_alloc /= 0 )   CALL ctl_warn('p4z_sed_alloc: failed to allocate arrays') 
     424      ! 
     425   END FUNCTION p4z_sed_alloc 
     426 
    424427 
    425428#else 
Note: See TracChangeset for help on using the changeset viewer.