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 7646 for trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90 – NEMO

Ignore:
Timestamp:
2017-02-06T10:25:03+01:00 (7 years ago)
Author:
timgraham
Message:

Merge of dev_merge_2016 into trunk. UPDATE TO ARCHFILES NEEDED for XIOS2.
LIM_SRC_s/limrhg.F90 to follow in next commit due to change of kind (I'm unable to do it in this commit).
Merged using the following steps:

1) svn merge --reintegrate ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk .
2) Resolve minor conflicts in sette.sh and namelist_cfg for ORCA2LIM3 (due to a change in trunk after branch was created)
3) svn commit
4) svn switch ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
5) svn merge ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2016/dev_merge_2016 .
6) At this stage I checked out a clean copy of the branch to compare against what is about to be committed to the trunk.
6) svn commit #Commit code to the trunk

In this commit I have also reverted a change to Fcheck_archfile.sh which was causing problems on the Paris machine.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90

    r6140 r7646  
    99   !!             3.5  !  2012-07 (O. Aumont) improvment of river input of nutrients  
    1010   !!---------------------------------------------------------------------- 
    11 #if defined key_pisces 
    12    !!---------------------------------------------------------------------- 
    13    !!   'key_pisces'                                       PISCES bio-model 
    14    !!---------------------------------------------------------------------- 
    1511   !!   p4z_sed        :  Compute loss of organic matter in the sediments 
    1612   !!---------------------------------------------------------------------- 
    … …  
    1814   USE trc             !  passive tracers common variables  
    1915   USE sms_pisces      !  PISCES Source Minus Sink variables 
    20    USE p4zsink         !  vertical flux of particulate matter due to sinking 
    21    USE p4zopt          !  optical model 
    2216   USE p4zlim          !  Co-limitations of differents nutrients 
    2317   USE p4zsbc          !  External source of nutrients  
    … …  
    5650      INTEGER, INTENT(in) ::   kt, knt ! ocean time step 
    5751      INTEGER  ::   ji, jj, jk, ikt 
    58 #if ! defined key_sed 
    5952      REAL(wp) ::   zsumsedsi, zsumsedpo4, zsumsedcal 
    6053      REAL(wp) ::   zrivalk, zrivsil, zrivno3 
    61 #endif 
    6254      REAL(wp) ::  zwflux, zfminus, zfplus 
    6355      REAL(wp) ::  zlim, zfact, zfactcal 
    6456      REAL(wp) ::  zo2, zno3, zflx, zpdenit, z1pdenit, zdenitt, zolimit 
    65       REAL(wp) ::  zsiloss, zcaloss, zws3, zws4, zwsc, zdep, zwstpoc 
    66       REAL(wp) ::  ztrfer, ztrpo4, zwdust, zlight 
     57      REAL(wp) ::  zsiloss, zcaloss, zws3, zws4, zwsc, zdep 
     58      REAL(wp) ::  zwstpoc, zwstpon, zwstpop 
     59      REAL(wp) ::  ztrfer, ztrpo4s, ztrdp, zwdust, zmudia, ztemp 
     60      REAL(wp) ::  xdiano3, xdianh4 
     61      REAL(wp) ::  zwssfep 
    6762      ! 
    6863      CHARACTER (len=25) :: charout 
    69       REAL(wp), POINTER, DIMENSION(:,:  ) :: zpdep, zsidep, zwork1, zwork2, zwork3 
     64      REAL(wp), POINTER, DIMENSION(:,:  ) :: zsidep, zwork1, zwork2, zwork3 
    7065      REAL(wp), POINTER, DIMENSION(:,:  ) :: zdenit2d, zironice, zbureff 
    7166      REAL(wp), POINTER, DIMENSION(:,:  ) :: zwsbio3, zwsbio4, zwscal 
    72       REAL(wp), POINTER, DIMENSION(:,:,:) :: zirondep, zsoufer 
     67      REAL(wp), POINTER, DIMENSION(:,:  ) :: zsedcal, zsedsi, zsedc 
     68      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrpo4, ztrdop, zirondep, zsoufer, zpdep, zlight 
     69      REAL(wp), POINTER, DIMENSION(:,:  ) :: zwsfep 
     70 
    7371      !!--------------------------------------------------------------------- 
    7472      ! 
    … …  
    7876      ! 
    7977      ! Allocate temporary workspace 
    80       CALL wrk_alloc( jpi, jpj, zdenit2d, zwork1, zwork2, zwork3, zbureff ) 
    81       CALL wrk_alloc( jpi, jpj, zwsbio3, zwsbio4, zwscal ) 
    82       CALL wrk_alloc( jpi, jpj, jpk, zsoufer ) 
     78                      CALL wrk_alloc( jpi, jpj, zdenit2d, zwork1, zwork2, zwork3, zbureff ) 
     79                      CALL wrk_alloc( jpi, jpj, zwsbio3, zwsbio4, zwscal ) 
     80                      CALL wrk_alloc( jpi, jpj, zsedcal,  zsedsi, zsedc ) 
     81                      CALL wrk_alloc( jpi, jpj, jpk, zlight, zsoufer ) 
     82      IF( ln_p5z )    CALL wrk_alloc( jpi, jpj, jpk, ztrpo4, ztrdop ) 
     83      IF( ln_ligand ) CALL wrk_alloc( jpi, jpj, zwsfep ) 
     84 
    8385 
    8486      zdenit2d(:,:) = 0.e0 
    … …  
    8789      zwork2  (:,:) = 0.e0 
    8890      zwork3  (:,:) = 0.e0 
     91      zsedsi  (:,:) = 0.e0 
     92      zsedcal (:,:) = 0.e0 
     93      zsedc   (:,:) = 0.e0 
     94 
    8995 
    9096      ! Iron input/uptake due to sea ice : Crude parameterization based on Lancelot et al. 
    … …  
    117123      IF( ln_dust ) THEN 
    118124         !                                               
    119          CALL wrk_alloc( jpi, jpj,      zpdep, zsidep ) 
    120          CALL wrk_alloc( jpi, jpj, jpk, zirondep      ) 
     125         CALL wrk_alloc( jpi, jpj,      zsidep ) 
     126         CALL wrk_alloc( jpi, jpj, jpk, zpdep, zirondep      ) 
    121127         !                                              ! Iron and Si deposition at the surface 
    122128         IF( ln_solub ) THEN 
    … …  
    125131            zirondep(:,:,1) = dustsolub  * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 55.85 + 3.e-10 * r1_ryyss  
    126132         ENDIF 
    127          zsidep(:,:) = 8.8 * 0.075 * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 28.1  
    128          zpdep (:,:) = 0.1 * 0.021 * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 31. / po4r  
     133         zsidep(:,:)   = 8.8 * 0.075 * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 28.1  
     134         zpdep (:,:,1) = 0.1 * 0.021 * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 31. / po4r  
    129135         !                                              ! Iron solubilization of particles in the water column 
    130136         !                                              ! dust in kg/m2/s ---> 1/55.85 to put in mol/Fe ;  wdust in m/j 
    … …  
    132138         DO jk = 2, jpkm1 
    133139            zirondep(:,:,jk) = dust(:,:) * mfrac * zwdust * rfact2 * EXP( -gdept_n(:,:,jk) / 540. ) 
     140            zpdep   (:,:,jk) = zirondep(:,:,jk) * 0.023 
    134141         END DO 
    135142         !                                              ! Iron solubilization of particles in the water column 
    136          tra(:,:,1,jppo4) = tra(:,:,1,jppo4) + zpdep   (:,:) 
    137143         tra(:,:,1,jpsil) = tra(:,:,1,jpsil) + zsidep  (:,:) 
     144         tra(:,:,:,jppo4) = tra(:,:,:,jppo4) + zpdep   (:,:,:) 
    138145         tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + zirondep(:,:,:)  
    139146         !  
    … …  
    145152                &  CALL iom_put( "pdust"  , dust(:,:) / ( wdust * rday )  * tmask(:,:,1) ) ! dust concentration at surface 
    146153            ENDIF 
    147          ELSE                                     
    148             IF( ln_diatrc )  & 
    149               &  trc2d(:,:,jp_pcs0_2d + 11) = zirondep(:,:,1) * 1.e+3 * rfact2r * e3t_n(:,:,1) * tmask(:,:,1) 
    150154         ENDIF 
    151          CALL wrk_dealloc( jpi, jpj,      zpdep, zsidep ) 
    152          CALL wrk_dealloc( jpi, jpj, jpk, zirondep      ) 
     155         CALL wrk_dealloc( jpi, jpj,      zsidep ) 
     156         CALL wrk_dealloc( jpi, jpj, jpk, zpdep, zirondep      ) 
    153157         !                                               
    154158      ENDIF 
    … …  
    169173            ENDDO 
    170174         ENDDO 
     175         IF( ln_p5z ) THEN 
     176            DO jj = 1, jpj 
     177               DO ji = 1, jpi 
     178                  DO jk = 1, nk_rnf(ji,jj) 
     179                     tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + rivdop(ji,jj) * rfact2 
     180                     tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + rivdon(ji,jj) * rfact2 
     181                     tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + rivdoc(ji,jj) * rfact2 
     182                  ENDDO 
     183               ENDDO 
     184            ENDDO 
     185         ENDIF 
    171186      ENDIF 
    172187       
    … …  
    181196      ! ------------------------------------------------------ 
    182197      IF( ln_ironsed ) THEN 
    183          tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + ironsed(:,:,:) * rfact2 
     198                         tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + ironsed(:,:,:) * rfact2 
     199         IF( ln_ligand ) tra(:,:,:,jpfep) = tra(:,:,:,jpfep) + ( ironsed(:,:,:) * fep_rats ) * rfact2 
    184200         ! 
    185201         IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "Ironsed" ) )   & 
    … …  
    190206      ! ------------------------------------------------------ 
    191207      IF( ln_hydrofe ) THEN 
    192          tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + hydrofe(:,:,:) * rfact2 
     208            tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + hydrofe(:,:,:) * rfact2 
     209         IF( ln_ligand ) THEN 
     210            tra(:,:,:,jpfep) = tra(:,:,:,jpfep) + ( hydrofe(:,:,:) * fep_rath ) * rfact2 
     211            tra(:,:,:,jplgw) = tra(:,:,:,jplgw) + ( hydrofe(:,:,:) * lgw_rath ) * rfact2 
     212         ENDIF 
    193213         ! 
    194214         IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "HYDR" ) )   & 
    … …  
    196216      ENDIF 
    197217 
    198       ! OA: Warning, the following part is necessary, especially with Kriest 
    199       ! to avoid CFL problems above the sediments 
     218      ! OA: Warning, the following part is necessary to avoid CFL problems above the sediments 
    200219      ! -------------------------------------------------------------------- 
    201220      DO jj = 1, jpj 
    … …  
    208227         END DO 
    209228      END DO 
    210  
    211 #if ! defined key_sed 
    212       ! Computation of the sediment denitrification proportion: The metamodel from midlleburg (2006) is being used 
    213       ! Computation of the fraction of organic matter that is permanently buried from Dunne's model 
    214       ! ------------------------------------------------------- 
    215       DO jj = 1, jpj 
    216          DO ji = 1, jpi 
    217            IF( tmask(ji,jj,1) == 1 ) THEN 
    218               ikt = mbkt(ji,jj) 
    219 # if defined key_kriest 
    220               zflx =    trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj)    * 1E3 * 1E6 / 1E4 
    221 # else 
    222               zflx = (  trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   & 
    223                 &     + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) )  * 1E3 * 1E6 / 1E4 
    224 #endif 
    225               zflx  = LOG10( MAX( 1E-3, zflx ) ) 
    226               zo2   = LOG10( MAX( 10. , trb(ji,jj,ikt,jpoxy) * 1E6 ) ) 
    227               zno3  = LOG10( MAX( 1.  , trb(ji,jj,ikt,jpno3) * 1E6 * rno3 ) ) 
    228               zdep  = LOG10( gdepw_n(ji,jj,ikt+1) ) 
    229               zdenit2d(ji,jj) = -2.2567 - 1.185 * zflx - 0.221 * zflx**2 - 0.3995 * zno3 * zo2 + 1.25 * zno3    & 
    230               &                + 0.4721 * zo2 - 0.0996 * zdep + 0.4256 * zflx * zo2 
    231               zdenit2d(ji,jj) = 10.0**( zdenit2d(ji,jj) ) 
    232               ! 
    233               zflx = (  trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   & 
    234                 &     + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) ) * 1E6 
    235               zbureff(ji,jj) = 0.013 + 0.53 * zflx**2 / ( 7.0 + zflx )**2 
    236            ENDIF 
    237          END DO 
    238       END DO  
    239  
    240       ! Loss of biogenic silicon, Caco3 organic carbon in the sediments.  
    241       ! First, the total loss is computed. 
    242       ! The factor for calcite comes from the alkalinity effect 
    243       ! ------------------------------------------------------------- 
    244       DO jj = 1, jpj 
    245          DO ji = 1, jpi 
    246             IF( tmask(ji,jj,1) == 1 ) THEN 
    247                ikt = mbkt(ji,jj)  
    248 # if defined key_kriest 
    249                zwork1(ji,jj) = trb(ji,jj,ikt,jpgsi) * zwscal (ji,jj) 
    250                zwork2(ji,jj) = trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) 
    251 # else 
    252                zwork1(ji,jj) = trb(ji,jj,ikt,jpgsi) * zwsbio4(ji,jj) 
    253                zwork2(ji,jj) = trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj) + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj)  
    254 # endif 
    255                ! For calcite, burial efficiency is made a function of saturation 
    256                zfactcal      = MIN( excess(ji,jj,ikt), 0.2 ) 
    257                zfactcal      = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) ) 
    258                zwork3(ji,jj) = trb(ji,jj,ikt,jpcal) * zwscal(ji,jj) * 2.e0 * zfactcal 
    259             ENDIF 
    260          END DO 
    261       END DO 
    262       zsumsedsi  = glob_sum( zwork1(:,:) * e1e2t(:,:) ) * r1_rday 
    263       zsumsedpo4 = glob_sum( zwork2(:,:) * e1e2t(:,:) ) * r1_rday 
    264       zsumsedcal = glob_sum( zwork3(:,:) * e1e2t(:,:) ) * r1_rday 
    265 #endif 
     229      ! 
     230      IF( ln_ligand ) THEN 
     231         DO jj = 1, jpj 
     232            DO ji = 1, jpi 
     233               ikt  = mbkt(ji,jj) 
     234               zdep = e3t_n(ji,jj,ikt) / xstep 
     235               zwsfep(ji,jj)  = MIN( 0.99 * zdep, wsfep(ji,jj,ikt)  ) 
     236            END DO 
     237         ENDDO 
     238      ENDIF 
     239 
     240      IF( .NOT.lk_sed ) THEN 
     241         ! Computation of the sediment denitrification proportion: The metamodel from midlleburg (2006) is being used 
     242         ! Computation of the fraction of organic matter that is permanently buried from Dunne's model 
     243         ! ------------------------------------------------------- 
     244         DO jj = 1, jpj 
     245            DO ji = 1, jpi 
     246              IF( tmask(ji,jj,1) == 1 ) THEN 
     247                 ikt = mbkt(ji,jj) 
     248                 zflx = (  trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   & 
     249                   &     + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) )  * 1E3 * 1E6 / 1E4 
     250                 zflx  = LOG10( MAX( 1E-3, zflx ) ) 
     251                 zo2   = LOG10( MAX( 10. , trb(ji,jj,ikt,jpoxy) * 1E6 ) ) 
     252                 zno3  = LOG10( MAX( 1.  , trb(ji,jj,ikt,jpno3) * 1E6 * rno3 ) ) 
     253                 zdep  = LOG10( gdepw_n(ji,jj,ikt+1) ) 
     254                 zdenit2d(ji,jj) = -2.2567 - 1.185 * zflx - 0.221 * zflx**2 - 0.3995 * zno3 * zo2 + 1.25 * zno3    & 
     255                   &                + 0.4721 * zo2 - 0.0996 * zdep + 0.4256 * zflx * zo2 
     256                 zdenit2d(ji,jj) = 10.0**( zdenit2d(ji,jj) ) 
     257                   ! 
     258                 zflx = (  trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   & 
     259                   &     + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) ) * 1E6 
     260                 zbureff(ji,jj) = 0.013 + 0.53 * zflx**2 / ( 7.0 + zflx )**2 
     261                ENDIF 
     262              END DO 
     263           END DO  
     264 
     265           ! Loss of biogenic silicon, Caco3 organic carbon in the sediments.  
     266           ! First, the total loss is computed. 
     267           ! The factor for calcite comes from the alkalinity effect 
     268           ! ------------------------------------------------------------- 
     269           DO jj = 1, jpj 
     270              DO ji = 1, jpi 
     271                 IF( tmask(ji,jj,1) == 1 ) THEN 
     272                    ikt = mbkt(ji,jj)  
     273                    zwork1(ji,jj) = trb(ji,jj,ikt,jpgsi) * zwsbio4(ji,jj) 
     274                    zwork2(ji,jj) = trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj) + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj)  
     275                    ! For calcite, burial efficiency is made a function of saturation 
     276                    zfactcal      = MIN( excess(ji,jj,ikt), 0.2 ) 
     277                    zfactcal      = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) ) 
     278                    zwork3(ji,jj) = trb(ji,jj,ikt,jpcal) * zwscal(ji,jj) * 2.e0 * zfactcal 
     279                ENDIF 
     280            END DO 
     281         END DO 
     282         zsumsedsi  = glob_sum( zwork1(:,:) * e1e2t(:,:) ) * r1_rday 
     283         zsumsedpo4 = glob_sum( zwork2(:,:) * e1e2t(:,:) ) * r1_rday 
     284         zsumsedcal = glob_sum( zwork3(:,:) * e1e2t(:,:) ) * r1_rday 
     285         ! 
     286      ENDIF 
    266287 
    267288      ! This loss is scaled at each bottom grid cell for equilibrating the total budget of silica in the ocean. 
    268289      ! Thus, the amount of silica lost in the sediments equal the supply at the surface (dust+rivers) 
    269290      ! ------------------------------------------------------ 
    270 #if ! defined key_sed 
    271       zrivsil =  1._wp - ( sumdepsi + rivdsiinput * r1_ryyss ) / ( zsumsedsi + rtrn ) 
    272 #endif 
     291      IF( .NOT.lk_sed )  zrivsil =  1._wp - ( sumdepsi + rivdsiinput * r1_ryyss ) / ( zsumsedsi + rtrn ) 
    273292 
    274293      DO jj = 1, jpj 
    … …  
    276295            ikt  = mbkt(ji,jj) 
    277296            zdep = xstep / e3t_n(ji,jj,ikt)  
    278             zws4 = zwsbio4(ji,jj) * zdep 
    279297            zwsc = zwscal (ji,jj) * zdep 
    280 # if defined key_kriest 
    281             zsiloss = trb(ji,jj,ikt,jpgsi) * zws4 
    282 # else 
    283298            zsiloss = trb(ji,jj,ikt,jpgsi) * zwsc 
    284 # endif 
    285299            zcaloss = trb(ji,jj,ikt,jpcal) * zwsc 
    286300            ! 
    287301            tra(ji,jj,ikt,jpgsi) = tra(ji,jj,ikt,jpgsi) - zsiloss 
    288302            tra(ji,jj,ikt,jpcal) = tra(ji,jj,ikt,jpcal) - zcaloss 
    289 #if ! defined key_sed 
    290             tra(ji,jj,ikt,jpsil) = tra(ji,jj,ikt,jpsil) + zsiloss * zrivsil  
    291             zfactcal = MIN( excess(ji,jj,ikt), 0.2 ) 
    292             zfactcal = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) ) 
    293             zrivalk  =  1._wp - ( rivalkinput * r1_ryyss ) * zfactcal / ( zsumsedcal + rtrn ) 
    294             tra(ji,jj,ikt,jptal) =  tra(ji,jj,ikt,jptal) + zcaloss * zrivalk * 2.0 
    295             tra(ji,jj,ikt,jpdic) =  tra(ji,jj,ikt,jpdic) + zcaloss * zrivalk 
    296 #endif 
    297303         END DO 
    298304      END DO 
    299  
     305      ! 
     306      IF( .NOT.lk_sed ) THEN 
     307         DO jj = 1, jpj 
     308            DO ji = 1, jpi 
     309               ikt  = mbkt(ji,jj) 
     310               zdep = xstep / e3t_n(ji,jj,ikt)  
     311               zwsc = zwscal (ji,jj) * zdep 
     312               zsiloss = trb(ji,jj,ikt,jpgsi) * zwsc 
     313               zcaloss = trb(ji,jj,ikt,jpcal) * zwsc 
     314               tra(ji,jj,ikt,jpsil) = tra(ji,jj,ikt,jpsil) + zsiloss * zrivsil  
     315               ! 
     316               zfactcal = MIN( excess(ji,jj,ikt), 0.2 ) 
     317               zfactcal = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) ) 
     318               zrivalk  =  1._wp - ( rivalkinput * r1_ryyss ) * zfactcal / ( zsumsedcal + rtrn ) 
     319               tra(ji,jj,ikt,jptal) =  tra(ji,jj,ikt,jptal) + zcaloss * zrivalk * 2.0 
     320               tra(ji,jj,ikt,jpdic) =  tra(ji,jj,ikt,jpdic) + zcaloss * zrivalk 
     321               zsedcal(ji,jj) = (1.0 - zrivalk) * zcaloss / zdep 
     322               zsedsi (ji,jj) = (1.0 - zrivsil) * zsiloss / zdep 
     323            END DO 
     324         END DO 
     325      ENDIF 
     326      ! 
    300327      DO jj = 1, jpj 
    301328         DO ji = 1, jpi 
    … …  
    304331            zws4 = zwsbio4(ji,jj) * zdep 
    305332            zws3 = zwsbio3(ji,jj) * zdep 
    306             zrivno3 = 1. - zbureff(ji,jj) 
    307 # if ! defined key_kriest 
    308333            tra(ji,jj,ikt,jpgoc) = tra(ji,jj,ikt,jpgoc) - trb(ji,jj,ikt,jpgoc) * zws4  
    309334            tra(ji,jj,ikt,jppoc) = tra(ji,jj,ikt,jppoc) - trb(ji,jj,ikt,jppoc) * zws3 
    310335            tra(ji,jj,ikt,jpbfe) = tra(ji,jj,ikt,jpbfe) - trb(ji,jj,ikt,jpbfe) * zws4 
    311336            tra(ji,jj,ikt,jpsfe) = tra(ji,jj,ikt,jpsfe) - trb(ji,jj,ikt,jpsfe) * zws3 
    312             zwstpoc              = trb(ji,jj,ikt,jpgoc) * zws4 + trb(ji,jj,ikt,jppoc) * zws3 
    313 # else 
    314             tra(ji,jj,ikt,jpnum) = tra(ji,jj,ikt,jpnum) - trb(ji,jj,ikt,jpnum) * zws4  
    315             tra(ji,jj,ikt,jppoc) = tra(ji,jj,ikt,jppoc) - trb(ji,jj,ikt,jppoc) * zws3 
    316             tra(ji,jj,ikt,jpsfe) = tra(ji,jj,ikt,jpsfe) - trb(ji,jj,ikt,jpsfe) * zws3 
    317             zwstpoc = trb(ji,jj,ikt,jppoc) * zws3  
    318 # endif 
    319  
    320 #if ! defined key_sed 
    321             ! The 0.5 factor in zpdenit and zdenitt is to avoid negative NO3 concentration after both denitrification 
    322             ! in the sediments and just above the sediments. Not very clever, but simpliest option. 
    323             zpdenit  = MIN( 0.5 * ( trb(ji,jj,ikt,jpno3) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3 ) 
    324             z1pdenit = zwstpoc * zrivno3 - zpdenit 
    325             zolimit = MIN( ( trb(ji,jj,ikt,jpoxy) - rtrn ) / o2ut, z1pdenit * ( 1.- nitrfac(ji,jj,ikt) ) ) 
    326             zdenitt = MIN(  0.5 * ( trb(ji,jj,ikt,jpno3) - rtrn ) / rdenit, z1pdenit * nitrfac(ji,jj,ikt) ) 
    327             tra(ji,jj,ikt,jpdoc) = tra(ji,jj,ikt,jpdoc) + z1pdenit - zolimit - zdenitt 
    328             tra(ji,jj,ikt,jppo4) = tra(ji,jj,ikt,jppo4) + zpdenit + zolimit + zdenitt 
    329             tra(ji,jj,ikt,jpnh4) = tra(ji,jj,ikt,jpnh4) + zpdenit + zolimit + zdenitt 
    330             tra(ji,jj,ikt,jpno3) = tra(ji,jj,ikt,jpno3) - rdenit * (zpdenit + zdenitt) 
    331             tra(ji,jj,ikt,jpoxy) = tra(ji,jj,ikt,jpoxy) - zolimit * o2ut 
    332             tra(ji,jj,ikt,jptal) = tra(ji,jj,ikt,jptal) + rno3 * (zolimit + (1.+rdenit) * (zpdenit + zdenitt) ) 
    333             tra(ji,jj,ikt,jpdic) = tra(ji,jj,ikt,jpdic) + zpdenit + zolimit + zdenitt 
    334             sdenit(ji,jj) = rdenit * zpdenit * e3t_n(ji,jj,ikt) 
    335 #endif 
    336337         END DO 
    337338      END DO 
     339      ! 
     340      IF( ln_ligand ) THEN 
     341         DO jj = 1, jpj 
     342            DO ji = 1, jpi 
     343               ikt     = mbkt(ji,jj) 
     344               zdep    = xstep / e3t_n(ji,jj,ikt)  
     345               zwssfep = zwsfep(ji,jj) * zdep 
     346               tra(ji,jj,ikt,jpfep) = tra(ji,jj,ikt,jpfep) - trb(ji,jj,ikt,jpfep) * zwssfep 
     347            END DO 
     348         END DO 
     349      ENDIF 
     350      ! 
     351      IF( ln_p5z ) THEN 
     352         DO jj = 1, jpj 
     353            DO ji = 1, jpi 
     354               ikt  = mbkt(ji,jj) 
     355               zdep = xstep / e3t_n(ji,jj,ikt)  
     356               zws4 = zwsbio4(ji,jj) * zdep 
     357               zws3 = zwsbio3(ji,jj) * zdep 
     358               tra(ji,jj,ikt,jpgon) = tra(ji,jj,ikt,jpgon) - trb(ji,jj,ikt,jpgon) * zws4 
     359               tra(ji,jj,ikt,jppon) = tra(ji,jj,ikt,jppon) - trb(ji,jj,ikt,jppon) * zws3 
     360               tra(ji,jj,ikt,jpgop) = tra(ji,jj,ikt,jpgop) - trb(ji,jj,ikt,jpgop) * zws4 
     361               tra(ji,jj,ikt,jppop) = tra(ji,jj,ikt,jppop) - trb(ji,jj,ikt,jppop) * zws3 
     362            END DO 
     363         END DO 
     364      ENDIF 
     365 
     366      IF( .NOT.lk_sed ) THEN 
     367         ! The 0.5 factor in zpdenit and zdenitt is to avoid negative NO3 concentration after both denitrification 
     368         ! in the sediments and just above the sediments. Not very clever, but simpliest option. 
     369         DO jj = 1, jpj 
     370            DO ji = 1, jpi 
     371               ikt  = mbkt(ji,jj) 
     372               zdep = xstep / e3t_n(ji,jj,ikt)  
     373               zws4 = zwsbio4(ji,jj) * zdep 
     374               zws3 = zwsbio3(ji,jj) * zdep 
     375               zrivno3 = 1. - zbureff(ji,jj) 
     376               zwstpoc = trb(ji,jj,ikt,jpgoc) * zws4 + trb(ji,jj,ikt,jppoc) * zws3 
     377               zpdenit  = MIN( 0.5 * ( trb(ji,jj,ikt,jpno3) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3 ) 
     378               z1pdenit = zwstpoc * zrivno3 - zpdenit 
     379               zolimit = MIN( ( trb(ji,jj,ikt,jpoxy) - rtrn ) / o2ut, z1pdenit * ( 1.- nitrfac(ji,jj,ikt) ) ) 
     380               zdenitt = MIN(  0.5 * ( trb(ji,jj,ikt,jpno3) - rtrn ) / rdenit, z1pdenit * nitrfac(ji,jj,ikt) ) 
     381               tra(ji,jj,ikt,jpdoc) = tra(ji,jj,ikt,jpdoc) + z1pdenit - zolimit - zdenitt 
     382               tra(ji,jj,ikt,jppo4) = tra(ji,jj,ikt,jppo4) + zpdenit + zolimit + zdenitt 
     383               tra(ji,jj,ikt,jpnh4) = tra(ji,jj,ikt,jpnh4) + zpdenit + zolimit + zdenitt 
     384               tra(ji,jj,ikt,jpno3) = tra(ji,jj,ikt,jpno3) - rdenit * (zpdenit + zdenitt) 
     385               tra(ji,jj,ikt,jpoxy) = tra(ji,jj,ikt,jpoxy) - zolimit * o2ut 
     386               tra(ji,jj,ikt,jptal) = tra(ji,jj,ikt,jptal) + rno3 * (zolimit + (1.+rdenit) * (zpdenit + zdenitt) ) 
     387               tra(ji,jj,ikt,jpdic) = tra(ji,jj,ikt,jpdic) + zpdenit + zolimit + zdenitt 
     388               sdenit(ji,jj) = rdenit * zpdenit * e3t_n(ji,jj,ikt) 
     389               zsedc(ji,jj)   = (1. - zrivno3) * zwstpoc / zdep 
     390               IF( ln_p5z ) THEN 
     391                  zwstpop              = trb(ji,jj,ikt,jpgop) * zws4 + trb(ji,jj,ikt,jppop) * zws3 
     392                  zwstpon              = trb(ji,jj,ikt,jpgon) * zws4 + trb(ji,jj,ikt,jppon) * zws3 
     393                  tra(ji,jj,ikt,jpdon) = tra(ji,jj,ikt,jpdon) + (z1pdenit - zolimit - zdenitt) * zwstpon / (zwstpoc + rtrn) 
     394                  tra(ji,jj,ikt,jpdop) = tra(ji,jj,ikt,jpdop) + (z1pdenit - zolimit - zdenitt) * zwstpop / (zwstpoc + rtrn) 
     395               ENDIF 
     396            END DO 
     397         END DO 
     398       ENDIF 
     399 
    338400 
    339401      ! Nitrogen fixation process 
    … …  
    341403      !----------------------------------- 
    342404      DO jk = 1, jpkm1 
    343          DO jj = 1, jpj 
    344             DO ji = 1, jpi 
    345                !                      ! Potential nitrogen fixation dependant on temperature and iron 
    346                zlim = ( 1.- xnanono3(ji,jj,jk) - xnanonh4(ji,jj,jk) ) 
    347                IF( zlim <= 0.2 )   zlim = 0.01 
    348 #if defined key_degrad 
    349                zfact = zlim * rfact2 * facvol(ji,jj,jk) 
    350 #else 
    351                zfact = zlim * rfact2 
    352 #endif 
    353                ztrfer = biron(ji,jj,jk)       / ( concfediaz + biron(ji,jj,jk)       ) 
    354                ztrpo4 = trb  (ji,jj,jk,jppo4) / ( concnnh4   + trb  (ji,jj,jk,jppo4) )  
    355                zlight =  ( 1.- EXP( -etot_ndcy(ji,jj,jk) / diazolight ) )  
    356                nitrpot(ji,jj,jk) =  MAX( 0.e0, ( 0.6 * tgfunc(ji,jj,jk) - 2.15 ) * r1_rday )   & 
    357                  &         *  zfact * MIN( ztrfer, ztrpo4 ) * zlight 
    358                zsoufer(ji,jj,jk) = zlight * 2E-11 / (2E-11 + biron(ji,jj,jk)) 
    359             END DO 
    360          END DO 
    361       END DO 
     405         zlight (:,:,jk) =  ( 1.- EXP( -etot_ndcy(:,:,jk) / diazolight ) ) * ( 1. - fr_i(:,:) )  
     406         zsoufer(:,:,jk) = zlight(:,:,jk) * 2E-11 / ( 2E-11 + biron(:,:,jk) ) 
     407      ENDDO 
     408      IF( ln_p4z ) THEN 
     409         DO jk = 1, jpkm1 
     410            DO jj = 1, jpj 
     411               DO ji = 1, jpi 
     412                  !                      ! Potential nitrogen fixation dependant on temperature and iron 
     413                  zlim = ( 1.- xnanono3(ji,jj,jk) - xnanonh4(ji,jj,jk) ) 
     414                  IF( zlim <= 0.2 )   zlim = 0.01 
     415                  zfact = zlim * rfact2 
     416 
     417                  ztrfer  = biron(ji,jj,jk)       / ( concfediaz + biron(ji,jj,jk)       ) 
     418                  ztrpo4s = trb  (ji,jj,jk,jppo4) / ( concnnh4   + trb  (ji,jj,jk,jppo4) )  
     419                  nitrpot(ji,jj,jk) =  MAX( 0.e0, ( 0.6 * tgfunc(ji,jj,jk) - 2.15 ) * r1_rday ) & 
     420                    &                *  zfact * MIN( ztrfer, ztrpo4s ) * zlight(ji,jj,jk) 
     421               END DO 
     422            END DO 
     423         END DO 
     424      ELSE       ! p5z 
     425         DO jk = 1, jpkm1 
     426            DO jj = 1, jpj 
     427               DO ji = 1, jpi 
     428                  !                      ! Potential nitrogen fixation dependant on temperature and iron 
     429                  ztemp = tsn(ji,jj,jk,jp_tem) 
     430                  zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625 
     431                  !       Potential nitrogen fixation dependant on temperature and iron 
     432                  xdianh4 = trb(ji,jj,jk,jpnh4) / ( concnnh4 + trb(ji,jj,jk,jpnh4) ) 
     433                  xdiano3 = trb(ji,jj,jk,jpno3) / ( concnno3 + trb(ji,jj,jk,jpno3) ) * (1. - xdianh4) 
     434                  zlim = ( 1.- xdiano3 - xdianh4 ) 
     435                  IF( zlim <= 0.1 )   zlim = 0.01 
     436                  zfact = zlim * rfact2 
     437                  ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) ) 
     438                  ztrpo4(ji,jj,jk) = trb(ji,jj,jk,jppo4) / ( 1E-6 + trb(ji,jj,jk,jppo4) ) 
     439                  ztrdop(ji,jj,jk) = trb(ji,jj,jk,jpdop) / ( 1E-6 + trb(ji,jj,jk,jpdop) ) * (1. - ztrpo4(ji,jj,jk)) 
     440                  ztrdp = ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) 
     441                  nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk) 
     442               END DO 
     443            END DO 
     444         END DO 
     445      ENDIF 
    362446 
    363447      ! Nitrogen change due to nitrogen fixation 
    364448      ! ---------------------------------------- 
    365       DO jk = 1, jpkm1 
    366          DO jj = 1, jpj 
    367             DO ji = 1, jpi 
    368                zfact = nitrpot(ji,jj,jk) * nitrfix 
    369                tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) +             zfact 
    370                tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3      * zfact 
    371                tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + o2nit     * zfact  
    372                tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + concdnh4 / ( concdnh4 + trb(ji,jj,jk,jppo4) ) & 
    373                &                     * 0.002 * trb(ji,jj,jk,jpdoc) * xstep 
    374                tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * xstep 
    375            END DO 
    376          END DO  
    377       END DO 
     449      IF( ln_p4z ) THEN 
     450         DO jk = 1, jpkm1 
     451            DO jj = 1, jpj 
     452               DO ji = 1, jpi 
     453                  zfact = nitrpot(ji,jj,jk) * nitrfix 
     454                  tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) +             zfact 
     455                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3      * zfact 
     456                  tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + o2nit     * zfact  
     457                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + concdnh4 / ( concdnh4 + trb(ji,jj,jk,jppo4) ) & 
     458                  &                     * 0.002 * trb(ji,jj,jk,jpdoc) * xstep 
     459                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * xstep 
     460              END DO 
     461            END DO  
     462         END DO 
     463      ELSE    ! p5z 
     464         DO jk = 1, jpkm1 
     465            DO jj = 1, jpj 
     466               DO ji = 1, jpi 
     467                  zfact = nitrpot(ji,jj,jk) * nitrfix 
     468                  tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zfact / 3.0 
     469                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * zfact / 3.0 
     470                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) - 16.0 / 46.0 * zfact * ( 1.0 - 1.0 / 3.0 ) & 
     471                  &                     * ztrpo4(ji,jj,jk) / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn) 
     472                  tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + zfact * 1.0 / 3.0 
     473                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zfact * 1.0 / 3.0 
     474                  tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + 16.0 / 46.0 * zfact / 3.0  & 
     475                  &                     - 16.0 / 46.0 * zfact * ztrdop(ji,jj,jk)   & 
     476                  &                     / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn) 
     477                  tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zfact * 1.0 / 3.0 * 2.0 / 3.0 
     478                  tra(ji,jj,jk,jppon) = tra(ji,jj,jk,jppon) + zfact * 1.0 / 3.0 * 2.0 /3.0 
     479                  tra(ji,jj,jk,jppop) = tra(ji,jj,jk,jppop) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 2.0 /3.0 
     480                  tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zfact * 1.0 / 3.0 * 1.0 / 3.0 
     481                  tra(ji,jj,jk,jpgon) = tra(ji,jj,jk,jpgon) + zfact * 1.0 / 3.0 * 1.0 /3.0 
     482                  tra(ji,jj,jk,jpgop) = tra(ji,jj,jk,jpgop) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 1.0 /3.0 
     483                  tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0 
     484                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - 30E-6 * zfact * 1.0 / 3.0  
     485                  tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0 
     486                  tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0 
     487                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday 
     488              END DO 
     489            END DO  
     490         END DO 
     491         ! 
     492      ENDIF 
    378493 
    379494      IF( lk_iomput ) THEN 
    … …  
    388503               CALL iom_put( "INTNFIX" , zwork1 )  
    389504            ENDIF 
     505            IF( iom_use("SedCal" ) ) CALL iom_put( "SedCal", zsedcal(:,:) * 1.e+3 ) 
     506            IF( iom_use("SedSi" ) )  CALL iom_put( "SedSi",  zsedsi (:,:) * 1.e+3 ) 
     507            IF( iom_use("SedC" ) )   CALL iom_put( "SedC",   zsedc  (:,:) * 1.e+3 ) 
     508            IF( iom_use("Sdenit" ) ) CALL iom_put( "Sdenit", sdenit (:,:) * 1.e+3 * rno3 ) 
    390509         ENDIF 
    391       ELSE 
    392          IF( ln_diatrc )  & 
    393             &  trc2d(:,:,jp_pcs0_2d + 12) = nitrpot(:,:,1) * nitrfix * rno3 * 1.e+3 * rfact2r * e3t_n(:,:,1) * tmask(:,:,1) 
    394510      ENDIF 
    395511      ! 
    … …  
    400516      ENDIF 
    401517      ! 
    402       CALL wrk_dealloc( jpi, jpj, zdenit2d, zwork1, zwork2, zwork3, zbureff ) 
    403       CALL wrk_dealloc( jpi, jpj, zwsbio3, zwsbio4, zwscal ) 
    404       CALL wrk_dealloc( jpi, jpj, jpk, zsoufer ) 
     518                      CALL wrk_dealloc( jpi, jpj, zdenit2d, zwork1, zwork2, zwork3, zbureff ) 
     519                      CALL wrk_dealloc( jpi, jpj, zwsbio3, zwsbio4, zwscal ) 
     520                      CALL wrk_dealloc( jpi, jpj, zsedcal,  zsedsi, zsedc ) 
     521                      CALL wrk_dealloc( jpi, jpj, jpk, zlight, zsoufer ) 
     522      IF( ln_p5z )    CALL wrk_dealloc( jpi, jpj, jpk, ztrpo4, ztrdop ) 
     523      IF( ln_ligand ) CALL wrk_dealloc( jpi, jpj, zwsfep ) 
    405524      ! 
    406525      IF( nn_timing == 1 )  CALL timing_stop('p4z_sed') 
    407       ! 
    408  9100  FORMAT(i8,3f10.5) 
    409526      ! 
    410527   END SUBROUTINE p4z_sed 
    … …  
    422539 
    423540 
    424 #else 
    425    !!====================================================================== 
    426    !!  Dummy module :                                   No PISCES bio-model 
    427    !!====================================================================== 
    428 CONTAINS 
    429    SUBROUTINE p4z_sed                         ! Empty routine 
    430    END SUBROUTINE p4z_sed 
    431 #endif  
    432  
    433541   !!====================================================================== 
    434542END MODULE p4zsed 
Note: See TracChangeset for help on using the changeset viewer.