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

Ignore:
Timestamp:
2012-01-28T17:44:18+01:00 (12 years ago)
Author:
rblod
Message:

Merge of 3.4beta into the trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zsed.F90

    r2774 r3294  
    66   !! History :   1.0  !  2004-03 (O. Aumont) Original code 
    77   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90 
     8   !!             3.4  !  2011-06  (O. Aumont, C. Ethe) USE of fldread 
    89   !!---------------------------------------------------------------------- 
    910#if defined key_pisces 
     
    1516   !!   p4z_sed_init   :  Initialization of p4z_sed 
    1617   !!---------------------------------------------------------------------- 
    17    USE trc 
    18    USE oce_trc         ! 
    19    USE sms_pisces 
    20    USE prtctl_trc 
    21    USE p4zbio 
    22    USE p4zint 
    23    USE p4zopt 
    24    USE p4zsink 
    25    USE p4zrem 
    26    USE p4zlim 
    27    USE iom 
    28  
     18   USE oce_trc         !  shared variables between ocean and passive tracers 
     19   USE trc             !  passive tracers common variables  
     20   USE sms_pisces      !  PISCES Source Minus Sink variables 
     21   USE p4zsink         !  vertical flux of particulate matter due to sinking 
     22   USE p4zopt          !  optical model 
     23   USE p4zlim          !  Co-limitations of differents nutrients 
     24   USE p4zrem          !  Remineralisation of organic matter 
     25   USE p4zint          !  interpolation and computation of various fields 
     26   USE iom             !  I/O manager 
     27   USE fldread         !  time interpolation 
     28   USE prtctl_trc      !  print control for debugging 
    2929 
    3030   IMPLICIT NONE 
     
    3636 
    3737   !! * Shared module variables 
    38    LOGICAL, PUBLIC :: ln_dustfer  = .FALSE.    !: boolean for dust input from the atmosphere 
    39    LOGICAL, PUBLIC :: ln_river    = .FALSE.    !: boolean for river input of nutrients 
    40    LOGICAL, PUBLIC :: ln_ndepo    = .FALSE.    !: boolean for atmospheric deposition of N 
    41    LOGICAL, PUBLIC :: ln_sedinput = .FALSE.    !: boolean for Fe input from sediments 
    42  
    43    REAL(wp), PUBLIC :: sedfeinput = 1.E-9_wp   !: Coastal release of Iron 
    44    REAL(wp), PUBLIC :: dustsolub  = 0.014_wp   !: Solubility of the dust 
     38   LOGICAL  :: ln_dust     = .FALSE.    !: boolean for dust input from the atmosphere 
     39   LOGICAL  :: ln_river    = .FALSE.    !: boolean for river input of nutrients 
     40   LOGICAL  :: ln_ndepo    = .FALSE.    !: boolean for atmospheric deposition of N 
     41   LOGICAL  :: ln_ironsed  = .FALSE.    !: boolean for Fe input from sediments 
     42 
     43   REAL(wp) :: sedfeinput  = 1.E-9_wp   !: Coastal release of Iron 
     44   REAL(wp) :: dustsolub   = 0.014_wp   !: Solubility of the dust 
     45   REAL(wp) :: wdust       = 2.0_wp     !: Sinking speed of the dust  
     46   REAL(wp) :: nitrfix     = 1E-7_wp    !: Nitrogen fixation rate    
     47   REAL(wp) :: diazolight  = 50._wp     !: Nitrogen fixation sensitivty to light  
     48   REAL(wp) :: concfediaz  = 1.E-10_wp  !: Fe half-saturation Cste for diazotrophs  
     49 
    4550 
    4651   !! * Module variables 
    4752   REAL(wp) :: ryyss                  !: number of seconds per year  
    48    REAL(wp) :: ryyss1                 !: inverse of ryyss 
     53   REAL(wp) :: r1_ryyss                 !: inverse of ryyss 
    4954   REAL(wp) :: rmtss                  !: number of seconds per month 
    50    REAL(wp) :: rday1                  !: inverse of rday 
    51  
    52    INTEGER , PARAMETER :: jpmth = 12  !: number of months per year 
    53    INTEGER , PARAMETER :: jpyr  = 1   !: one year 
    54  
    55    INTEGER ::  numdust                !: logical unit for surface fluxes data 
    56    INTEGER ::  nflx1 , nflx2          !: first and second record used 
    57    INTEGER ::  nflx11, nflx12 
    58  
    59    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dustmo    !: set of dust fields 
     55   REAL(wp) :: r1_rday                  !: inverse of rday 
     56   LOGICAL  :: ll_sbc 
     57 
     58   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_dust      ! structure of input dust 
     59   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_riverdic  ! structure of input riverdic 
     60   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_riverdoc  ! structure of input riverdoc 
     61   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_ndepo     ! structure of input nitrogen deposition 
     62   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_ironsed   ! structure of input iron from sediment 
     63 
     64   INTEGER , PARAMETER :: nbtimes = 365  !: maximum number of times record in a file 
     65   INTEGER  :: ntimes_dust, ntimes_riv, ntimes_ndep       ! number of time steps in a file 
     66 
    6067   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:) :: dust      !: dust fields 
    6168   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:) :: rivinp, cotdep    !: river input fields 
     
    8693      !! ** Method  : - ??? 
    8794      !!--------------------------------------------------------------------- 
    88       USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 
    89       USE wrk_nemo, ONLY: zsidep => wrk_2d_1, zwork => wrk_2d_2, zwork1 => wrk_2d_3 
    90       USE wrk_nemo, ONLY: znitrpot => wrk_3d_2, zirondep => wrk_3d_3 
    9195      ! 
    9296      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step 
     
    96100      REAL(wp) ::   zrivalk, zrivsil, zrivpo4 
    97101#endif 
    98       REAL(wp) ::   zdenitot, znitrpottot, zlim, zfact 
    99       REAL(wp) ::   zwsbio3, zwsbio4, zwscal 
     102      REAL(wp) ::   zdenitot, znitrpottot, zlim, zfact, zfactcal 
     103      REAL(wp) ::   zsiloss, zcaloss, zwsbio3, zwsbio4, zwscal, zdep 
    100104      CHARACTER (len=25) :: charout 
     105      REAL(wp), POINTER, DIMENSION(:,:  ) :: zsidep, zwork1, zwork2, zwork3 
     106      REAL(wp), POINTER, DIMENSION(:,:,:) :: znitrpot, zirondep 
    101107      !!--------------------------------------------------------------------- 
    102  
    103       IF( ( wrk_in_use(2, 1,2,3) ) .OR. ( wrk_in_use(3, 2,3) ) ) THEN 
    104          CALL ctl_stop('p4z_sed: requested workspace arrays unavailable')  ;  RETURN 
    105       END IF 
    106  
    107       IF( jnt == 1  .AND.  ln_dustfer  )  CALL p4z_sbc( kt ) 
     108      ! 
     109      IF( nn_timing == 1 )  CALL timing_start('p4z_sed') 
     110      ! 
     111      ! Allocate temporary workspace 
     112      CALL wrk_alloc( jpi, jpj,      zsidep, zwork1, zwork2, zwork3 ) 
     113      CALL wrk_alloc( jpi, jpj, jpk, znitrpot, zirondep             ) 
     114 
     115      IF( jnt == 1 .AND. ll_sbc ) CALL p4z_sbc( kt ) 
     116 
     117      zirondep(:,:,:) = 0.e0          ! Initialisation of variables USEd to compute deposition 
     118      zsidep  (:,:)   = 0.e0 
    108119 
    109120      ! Iron and Si deposition at the surface 
    110121      ! ------------------------------------- 
    111  
    112122      DO jj = 1, jpj 
    113123         DO ji = 1, jpi 
    114             zirondep(ji,jj,1) = ( dustsolub * dust(ji,jj) / ( 55.85 * rmtss ) + 3.e-10 * ryyss1 )   & 
    115                &             * rfact2 / fse3t(ji,jj,1) 
    116             zsidep  (ji,jj)   = 8.8 * 0.075 * dust(ji,jj) * rfact2 / ( fse3t(ji,jj,1) * 28.1 * rmtss ) 
     124            zdep  = rfact2 / fse3t(ji,jj,1) 
     125            zirondep(ji,jj,1) = ( dustsolub * dust(ji,jj) / ( 55.85 * rmtss ) + 3.e-10 * r1_ryyss ) * zdep 
     126            zsidep  (ji,jj)   = 8.8 * 0.075 * dust(ji,jj) * zdep / ( 28.1 * rmtss ) 
    117127         END DO 
    118128      END DO 
     
    120130      ! Iron solubilization of particles in the water column 
    121131      ! ---------------------------------------------------- 
    122  
    123132      DO jk = 2, jpkm1 
    124          zirondep(:,:,jk) = dust(:,:) / ( 10. * 55.85 * rmtss ) * rfact2 * 1.e-4 
     133         zirondep(:,:,jk) = dust(:,:) / ( wdust * 55.85 * rmtss ) * rfact2 * 1.e-4 * EXP( -fsdept(:,:,jk) / 1000. ) 
    125134      END DO 
    126135 
    127136      ! Add the external input of nutrients, carbon and alkalinity 
    128137      ! ---------------------------------------------------------- 
    129  
    130138      trn(:,:,1,jppo4) = trn(:,:,1,jppo4) + rivinp(:,:) * rfact2  
    131139      trn(:,:,1,jpno3) = trn(:,:,1,jpno3) + (rivinp(:,:) + nitdep(:,:)) * rfact2 
     
    139147      ! (dust, river and sediment mobilization) 
    140148      ! ------------------------------------------------------ 
    141  
    142149      DO jk = 1, jpkm1 
    143150         trn(:,:,jk,jpfer) = trn(:,:,jk,jpfer) + zirondep(:,:,jk) + ironsed(:,:,jk) * rfact2 
    144151      END DO 
    145  
    146152 
    147153#if ! defined key_sed 
     
    154160            ikt = mbkt(ji,jj)  
    155161# if defined key_kriest 
    156             zwork (ji,jj) = trn(ji,jj,ikt,jpdsi) * wscal (ji,jj,ikt) 
    157             zwork1(ji,jj) = trn(ji,jj,ikt,jppoc) * wsbio3(ji,jj,ikt) 
     162            zwork1(ji,jj) = trn(ji,jj,ikt,jpdsi) * wscal (ji,jj,ikt) 
     163            zwork2(ji,jj) = trn(ji,jj,ikt,jppoc) * wsbio3(ji,jj,ikt) 
    158164# else 
    159             zwork (ji,jj) = trn(ji,jj,ikt,jpdsi) * wsbio4(ji,jj,ikt) 
    160             zwork1(ji,jj) = trn(ji,jj,ikt,jpgoc) * wsbio4(ji,jj,ikt) + trn(ji,jj,ikt,jppoc) * wsbio3(ji,jj,ikt)  
     165            zwork1(ji,jj) = trn(ji,jj,ikt,jpdsi) * wsbio4(ji,jj,ikt) 
     166            zwork2(ji,jj) = trn(ji,jj,ikt,jpgoc) * wsbio4(ji,jj,ikt) + trn(ji,jj,ikt,jppoc) * wsbio3(ji,jj,ikt)  
    161167# endif 
    162          END DO 
    163       END DO 
    164       zsumsedsi  = glob_sum( zwork (:,:) * e1e2t(:,:) ) * rday1 
    165       zsumsedpo4 = glob_sum( zwork1(:,:) * e1e2t(:,:) ) * rday1 
    166       DO jj = 1, jpj 
    167          DO ji = 1, jpi 
    168             ikt = mbkt(ji,jj)  
    169             zwork (ji,jj) = trn(ji,jj,ikt,jpcal) * wscal (ji,jj,ikt) 
    170          END DO 
    171       END DO 
    172       zsumsedcal = glob_sum( zwork (:,:) * e1e2t(:,:) ) * 2.0 * rday1 
     168            ! For calcite, burial efficiency is made a function of saturation 
     169            zfactcal      = MIN( excess(ji,jj,ikt), 0.2 ) 
     170            zfactcal      = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) ) 
     171            zwork3(ji,jj) = trn(ji,jj,ikt,jpcal) * wscal (ji,jj,ikt) * 2.e0 * zfactcal 
     172         END DO 
     173      END DO 
     174      zsumsedsi  = glob_sum( zwork1(:,:) * e1e2t(:,:) ) * r1_rday 
     175      zsumsedpo4 = glob_sum( zwork2(:,:) * e1e2t(:,:) ) * r1_rday 
     176      zsumsedcal = glob_sum( zwork3(:,:) * e1e2t(:,:) ) * r1_rday 
    173177#endif 
    174178 
    175       ! Then this loss is scaled at each bottom grid cell for 
     179      ! THEN this loss is scaled at each bottom grid cell for 
    176180      ! equilibrating the total budget of silica in the ocean. 
    177181      ! Thus, the amount of silica lost in the sediments equal 
    178182      ! the supply at the surface (dust+rivers) 
    179183      ! ------------------------------------------------------ 
     184#if ! defined key_sed 
     185      zrivsil =  1._wp - ( sumdepsi + rivalkinput * r1_ryyss / 6. ) / zsumsedsi  
     186      zrivpo4 =  1._wp - ( rivpo4input * r1_ryyss ) / zsumsedpo4  
     187#endif 
    180188 
    181189      DO jj = 1, jpj 
    182190         DO ji = 1, jpi 
    183             ikt = mbkt(ji,jj) 
    184             zfact = xstep / fse3t(ji,jj,ikt) 
    185             zwsbio3 = 1._wp - zfact * wsbio3(ji,jj,ikt) 
    186             zwsbio4 = 1._wp - zfact * wsbio4(ji,jj,ikt) 
    187             zwscal  = 1._wp - zfact * wscal (ji,jj,ikt) 
     191            ikt  = mbkt(ji,jj) 
     192            zdep = xstep / fse3t(ji,jj,ikt) 
     193            zwsbio4 = wsbio4(ji,jj,ikt) * zdep 
     194            zwscal  = wscal (ji,jj,ikt) * zdep 
     195# if defined key_kriest 
     196            zsiloss = trn(ji,jj,ikt,jpdsi) * zwsbio4 
     197# else 
     198            zsiloss = trn(ji,jj,ikt,jpdsi) * zwscal 
     199# endif 
     200            zcaloss = trn(ji,jj,ikt,jpcal) * zwscal 
    188201            ! 
    189 # if defined key_kriest 
    190             trn(ji,jj,ikt,jpdsi) = trn(ji,jj,ikt,jpdsi) * zwsbio4 
    191             trn(ji,jj,ikt,jpnum) = trn(ji,jj,ikt,jpnum) * zwsbio4 
    192             trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc) * zwsbio3 
    193             trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe) * zwsbio3 
    194 # else 
    195             trn(ji,jj,ikt,jpdsi) = trn(ji,jj,ikt,jpdsi) * zwscal  
    196             trn(ji,jj,ikt,jpgoc) = trn(ji,jj,ikt,jpgoc) * zwsbio4 
    197             trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc) * zwsbio3 
    198             trn(ji,jj,ikt,jpbfe) = trn(ji,jj,ikt,jpbfe) * zwsbio4 
    199             trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe) * zwsbio3 
    200 # endif 
    201             trn(ji,jj,ikt,jpcal) = trn(ji,jj,ikt,jpcal) * zwscal 
    202          END DO 
    203       END DO 
    204  
     202            trn(ji,jj,ikt,jpdsi) = trn(ji,jj,ikt,jpdsi) - zsiloss 
     203            trn(ji,jj,ikt,jpcal) = trn(ji,jj,ikt,jpcal) - zcaloss 
    205204#if ! defined key_sed 
    206       zrivsil =  1._wp - ( sumdepsi + rivalkinput * ryyss1 / 6. ) / zsumsedsi  
    207       zrivalk =  1._wp - ( rivalkinput * ryyss1 ) / zsumsedcal  
    208       zrivpo4 =  1._wp - ( rivpo4input * ryyss1 ) / zsumsedpo4  
     205            trn(ji,jj,ikt,jpsil) = trn(ji,jj,ikt,jpsil) + zsiloss * zrivsil  
     206            zfactcal = MIN( excess(ji,jj,ikt), 0.2 ) 
     207            zfactcal = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) ) 
     208            zrivalk  =  1._wp - ( rivalkinput * r1_ryyss ) * zfactcal / zsumsedcal  
     209            trn(ji,jj,ikt,jptal) =  trn(ji,jj,ikt,jptal) + zcaloss * zrivalk * 2.0 
     210            trn(ji,jj,ikt,jpdic) =  trn(ji,jj,ikt,jpdic) + zcaloss * zrivalk 
     211#endif 
     212         END DO 
     213      END DO 
     214 
    209215      DO jj = 1, jpj 
    210216         DO ji = 1, jpi 
    211             ikt = mbkt(ji,jj) 
    212             zfact = xstep / fse3t(ji,jj,ikt) 
    213             zwsbio3 = zfact * wsbio3(ji,jj,ikt) 
    214             zwsbio4 = zfact * wsbio4(ji,jj,ikt) 
    215             zwscal  = zfact * wscal (ji,jj,ikt) 
    216             trn(ji,jj,ikt,jptal) =  trn(ji,jj,ikt,jptal) + trn(ji,jj,ikt,jpcal) * zwscal  * zrivalk * 2.0 
    217             trn(ji,jj,ikt,jpdic) =  trn(ji,jj,ikt,jpdic) + trn(ji,jj,ikt,jpcal) * zwscal  * zrivalk 
    218 # if defined key_kriest 
    219             trn(ji,jj,ikt,jpsil) =  trn(ji,jj,ikt,jpsil) + trn(ji,jj,ikt,jpdsi) * zwsbio4 * zrivsil  
    220             trn(ji,jj,ikt,jpdoc) =  trn(ji,jj,ikt,jpdoc) + trn(ji,jj,ikt,jppoc) * zwsbio3 * zrivpo4  
     217            ikt  = mbkt(ji,jj) 
     218            zdep = xstep / fse3t(ji,jj,ikt) 
     219            zwsbio4 = wsbio4(ji,jj,ikt) * zdep 
     220            zwsbio3 = wsbio3(ji,jj,ikt) * zdep 
     221# if ! defined key_kriest 
     222            trn(ji,jj,ikt,jpgoc) = trn(ji,jj,ikt,jpgoc) - trn(ji,jj,ikt,jpgoc) * zwsbio4 
     223            trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc) - trn(ji,jj,ikt,jppoc) * zwsbio3 
     224            trn(ji,jj,ikt,jpbfe) = trn(ji,jj,ikt,jpbfe) - trn(ji,jj,ikt,jpbfe) * zwsbio4 
     225            trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe) - trn(ji,jj,ikt,jpsfe) * zwsbio3 
     226#if ! defined key_sed 
     227            trn(ji,jj,ikt,jpdoc) = trn(ji,jj,ikt,jpdoc) & 
     228               &               + ( trn(ji,jj,ikt,jpgoc) * zwsbio4 + trn(ji,jj,ikt,jppoc) * zwsbio3 ) * zrivpo4 
     229#endif 
     230 
    221231# else 
    222             trn(ji,jj,ikt,jpsil) =  trn(ji,jj,ikt,jpsil) + trn(ji,jj,ikt,jpdsi) * zwscal  * zrivsil  
    223             trn(ji,jj,ikt,jpdoc) =  trn(ji,jj,ikt,jpdoc)   & 
    224             &                     + ( trn(ji,jj,ikt,jppoc) * zwsbio3 + trn(ji,jj,ikt,jpgoc) * zwsbio4 ) * zrivpo4 
     232            trn(ji,jj,ikt,jpnum) = trn(ji,jj,ikt,jpnum) - trn(ji,jj,ikt,jpnum) * zwsbio4 
     233            trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc) - trn(ji,jj,ikt,jppoc) * zwsbio3 
     234            trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe) - trn(ji,jj,ikt,jpsfe) * zwsbio3 
     235#if ! defined key_sed 
     236            trn(ji,jj,ikt,jpdoc) = trn(ji,jj,ikt,jpdoc) & 
     237               &               + ( trn(ji,jj,ikt,jpnum) * zwsbio4 + trn(ji,jj,ikt,jppoc) * zwsbio3 ) * zrivpo4 
     238#endif 
     239 
    225240# endif 
    226241         END DO 
    227242      END DO 
    228 # endif 
     243 
    229244 
    230245      ! Nitrogen fixation (simple parameterization). The total gain 
     
    233248      ! ------------------------------------------------------------- 
    234249 
    235       zdenitot = glob_sum( denitr(:,:,:)  * cvol(:,:,:) * xnegtr(:,:,:) ) * rdenit 
     250      zdenitot = glob_sum(  ( denitr(:,:,:) * rdenit + denitnh4(:,:,:) * rdenita ) * cvol(:,:,:) )  
    236251 
    237252      ! Potential nitrogen fixation dependant on temperature and iron 
     
    246261               zlim = ( 1.- xnanono3(ji,jj,jk) - xnanonh4(ji,jj,jk) ) 
    247262               IF( zlim <= 0.2 )   zlim = 0.01 
    248                znitrpot(ji,jj,jk) = MAX( 0.e0, ( 0.6 * tgfunc(ji,jj,jk) - 2.15 ) * rday1 )   & 
    249 # if defined key_degrad 
    250                &                  * facvol(ji,jj,jk)   & 
    251 # endif 
    252                &                  * zlim * rfact2 * trn(ji,jj,jk,jpfer)   & 
    253                &                  / ( conc3 + trn(ji,jj,jk,jpfer) ) * ( 1.- EXP( -etot(ji,jj,jk) / 50.) ) 
     263#if defined key_degrad 
     264               zfact = zlim * rfact2 * facvol(ji,jj,jk) 
     265#else 
     266               zfact = zlim * rfact2  
     267#endif 
     268               znitrpot(ji,jj,jk) =  MAX( 0.e0, ( 0.6 * tgfunc(ji,jj,jk) - 2.15 ) * r1_rday )   & 
     269                 &                 *  zfact * trn(ji,jj,jk,jpfer) / ( concfediaz + trn(ji,jj,jk,jpfer) ) & 
     270                 &                 * ( 1.- EXP( -etot(ji,jj,jk) / diazolight ) ) 
    254271            END DO 
    255272         END DO  
     
    260277      ! Nitrogen change due to nitrogen fixation 
    261278      ! ---------------------------------------- 
    262  
    263279      DO jk = 1, jpk 
    264280         DO jj = 1, jpj 
    265281            DO ji = 1, jpi 
    266                zfact = znitrpot(ji,jj,jk) * 1.e-7 
     282               zfact = znitrpot(ji,jj,jk) * nitrfix 
    267283               trn(ji,jj,jk,jpnh4) = trn(ji,jj,jk,jpnh4) + zfact 
     284               trn(ji,jj,jk,jptal) = trn(ji,jj,jk,jptal) + rno3 * zfact 
    268285               trn(ji,jj,jk,jpoxy) = trn(ji,jj,jk,jpoxy) + zfact   * o2nit 
    269                trn(ji,jj,jk,jppo4) = trn(ji,jj,jk,jppo4) + 30./ 46.* zfact 
    270             END DO 
    271          END DO 
    272       END DO 
    273  
    274 #if defined key_diatrc 
    275       zfact = 1.e+3 * rfact2r 
    276 #  if  ! defined key_iomput 
    277       trc2d(:,:,jp_pcs0_2d + 11) = zirondep(:,:,1)         * zfact * fse3t(:,:,1) * tmask(:,:,1) 
    278       trc2d(:,:,jp_pcs0_2d + 12) = znitrpot(:,:,1) * 1.e-7 * zfact * fse3t(:,:,1) * tmask(:,:,1) 
    279 #  else 
    280       zwork (:,:)  =  ( zirondep(:,:,1) + ironsed(:,:,1) * rfact2 ) * zfact * fse3t(:,:,1) * tmask(:,:,1)  
    281       zwork1(:,:)  =  znitrpot(:,:,1) * 1.e-7                       * zfact * fse3t(:,:,1) * tmask(:,:,1) 
    282       IF( jnt == nrdttrc ) THEN 
    283          CALL iom_put( "Irondep", zwork  )  ! surface downward net flux of iron 
    284          CALL iom_put( "Nfix"   , zwork1 )  ! nitrogen fixation at surface 
    285       ENDIF 
    286 #  endif 
    287 #endif 
    288       ! 
    289        IF(ln_ctl)   THEN  ! print mean trends (used for debugging) 
    290          WRITE(charout, FMT="('sed ')") 
     286               trn(ji,jj,jk,jppo4) = trn(ji,jj,jk,jppo4) + 30. / 46. * zfact 
     287           !    trn(ji,jj,jk,jppo4) = trn(ji,jj,jk,jppo4) + zfact 
     288           END DO 
     289         END DO  
     290      END DO 
     291      ! 
     292      IF( ln_diatrc ) THEN 
     293         zfact = 1.e+3 * rfact2r 
     294         IF( lk_iomput ) THEN 
     295            zwork1(:,:)  =  ( zirondep(:,:,1) + ironsed(:,:,1) * rfact2 ) * zfact * fse3t(:,:,1) * tmask(:,:,1)  
     296            zwork2(:,:)  =    znitrpot(:,:,1) * nitrfix                   * zfact * fse3t(:,:,1) * tmask(:,:,1) 
     297            IF( jnt == nrdttrc ) THEN 
     298               CALL iom_put( "Irondep", zwork1  )  ! surface downward net flux of iron 
     299               CALL iom_put( "Nfix"   , zwork2 )  ! nitrogen fixation at surface 
     300            ENDIF 
     301         ELSE 
     302            trc2d(:,:,jp_pcs0_2d + 11) = zirondep(:,:,1)           * zfact * fse3t(:,:,1) * tmask(:,:,1) 
     303            trc2d(:,:,jp_pcs0_2d + 12) = znitrpot(:,:,1) * nitrfix * zfact * fse3t(:,:,1) * tmask(:,:,1) 
     304         ENDIF 
     305      ENDIF 
     306      ! 
     307      IF(ln_ctl) THEN  ! print mean trends (USEd for debugging) 
     308         WRITE(charout, fmt="('sed ')") 
    291309         CALL prt_ctl_trc_info(charout) 
    292310         CALL prt_ctl_trc(tab4d=trn, mask=tmask, clinfo=ctrcnm) 
    293        ENDIF 
    294  
    295       IF( ( wrk_not_released(2, 1,2,3) ) .OR. ( wrk_not_released(3, 2,3) ) )   & 
    296         &         CALL ctl_stop('p4z_sed: failed to release workspace arrays') 
    297  
     311      ENDIF 
     312      ! 
     313      CALL wrk_dealloc( jpi, jpj,      zsidep, zwork1, zwork2, zwork3 ) 
     314      CALL wrk_dealloc( jpi, jpj, jpk, znitrpot, zirondep             ) 
     315      ! 
     316      IF( nn_timing == 1 )  CALL timing_stop('p4z_sed') 
     317      ! 
    298318   END SUBROUTINE p4z_sed 
    299319 
    300320   SUBROUTINE p4z_sbc( kt ) 
    301  
    302321      !!---------------------------------------------------------------------- 
    303       !!                  ***  ROUTINE p4z_sbc  *** 
    304       !! 
    305       !! ** Purpose :   Read and interpolate the external sources of  
     322      !!                  ***  routine p4z_sbc  *** 
     323      !! 
     324      !! ** purpose :   read and interpolate the external sources of  
    306325      !!                nutrients 
    307326      !! 
    308       !! ** Method  :   Read the files and interpolate the appropriate variables 
     327      !! ** method  :   read the files and interpolate the appropriate variables 
    309328      !! 
    310329      !! ** input   :   external netcdf files 
     
    314333      INTEGER, INTENT( in  ) ::   kt   ! ocean time step 
    315334 
    316       !! * Local declarations 
    317       INTEGER :: imois, i15, iman  
    318       REAL(wp) :: zxy 
    319  
     335      !! * local declarations 
     336      INTEGER  :: ji,jj  
     337      REAL(wp) :: zcoef 
    320338      !!--------------------------------------------------------------------- 
    321  
    322       ! Initialization 
    323       ! -------------- 
    324  
    325       i15 = nday / 16 
    326       iman  = INT( raamo ) 
    327       imois = nmonth + i15 - 1 
    328       IF( imois == 0 ) imois = iman 
    329  
    330       ! Calendar computation 
    331       IF( kt == nit000 .OR. imois /= nflx1 ) THEN 
    332  
    333          IF( kt == nit000 )  nflx1  = 0 
    334  
    335          ! nflx1 number of the first file record used in the simulation 
    336          ! nflx2 number of the last  file record 
    337  
    338          nflx1 = imois 
    339          nflx2 = nflx1 + 1 
    340          nflx1 = MOD( nflx1, iman ) 
    341          nflx2 = MOD( nflx2, iman ) 
    342          IF( nflx1 == 0 )   nflx1 = iman 
    343          IF( nflx2 == 0 )   nflx2 = iman 
    344          IF(lwp) WRITE(numout,*)  
    345          IF(lwp) WRITE(numout,*) ' p4z_sbc : first record file used nflx1 ',nflx1 
    346          IF(lwp) WRITE(numout,*) ' p4z_sbc : last  record file used nflx2 ',nflx2 
    347  
    348       ENDIF 
    349  
    350       ! 3. at every time step interpolation of fluxes 
    351       ! --------------------------------------------- 
    352  
    353       zxy = FLOAT( nday + 15 - 30 * i15 ) / 30 
    354       dust(:,:) = ( (1.-zxy) * dustmo(:,:,nflx1) + zxy * dustmo(:,:,nflx2) ) 
    355  
     339      ! 
     340      IF( nn_timing == 1 )  CALL timing_start('p4z_sbc') 
     341      ! 
     342      ! Compute dust at nit000 or only if there is more than 1 time record in dust file 
     343      IF( ln_dust ) THEN 
     344         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_dust > 1 ) ) THEN 
     345            CALL fld_read( kt, 1, sf_dust ) 
     346            dust(:,:) = sf_dust(1)%fnow(:,:,1) 
     347         ENDIF 
     348      ENDIF 
     349 
     350      ! N/P and Si releases due to coastal rivers 
     351      ! Compute river at nit000 or only if there is more than 1 time record in river file 
     352      ! ----------------------------------------- 
     353      IF( ln_river ) THEN 
     354         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_riv > 1 ) ) THEN 
     355            CALL fld_read( kt, 1, sf_riverdic ) 
     356            CALL fld_read( kt, 1, sf_riverdoc ) 
     357            DO jj = 1, jpj 
     358               DO ji = 1, jpi 
     359                  zcoef = ryyss * cvol(ji,jj,1)  
     360                  cotdep(ji,jj) =   sf_riverdic(1)%fnow(ji,jj,1)                                  * 1E9 / ( 12. * zcoef + rtrn ) 
     361                  rivinp(ji,jj) = ( sf_riverdic(1)%fnow(ji,jj,1) + sf_riverdoc(1)%fnow(ji,jj,1) ) * 1E9 / ( 31.6* zcoef + rtrn ) 
     362               END DO 
     363            END DO 
     364         ENDIF 
     365      ENDIF 
     366 
     367      ! Compute N deposition at nit000 or only if there is more than 1 time record in N deposition file 
     368      IF( ln_ndepo ) THEN 
     369         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_ndep > 1 ) ) THEN 
     370            CALL fld_read( kt, 1, sf_ndepo ) 
     371            DO jj = 1, jpj 
     372               DO ji = 1, jpi 
     373                  nitdep(ji,jj) = 7.6 * sf_ndepo(1)%fnow(ji,jj,1) / ( 14E6 * ryyss * fse3t(ji,jj,1) + rtrn ) 
     374               END DO 
     375            END DO 
     376         ENDIF 
     377      ENDIF 
     378      ! 
     379      IF( nn_timing == 1 )  CALL timing_stop('p4z_sbc') 
     380      ! 
    356381   END SUBROUTINE p4z_sbc 
    357382 
    358  
    359383   SUBROUTINE p4z_sed_init 
    360384 
    361385      !!---------------------------------------------------------------------- 
    362       !!                  ***  ROUTINE p4z_sed_init  *** 
    363       !! 
    364       !! ** Purpose :   Initialization of the external sources of nutrients 
    365       !! 
    366       !! ** Method  :   Read the files and compute the budget 
    367       !!      called at the first timestep (nit000) 
     386      !!                  ***  routine p4z_sed_init  *** 
     387      !! 
     388      !! ** purpose :   initialization of the external sources of nutrients 
     389      !! 
     390      !! ** method  :   read the files and compute the budget 
     391      !!                called at the first timestep (nittrc000) 
    368392      !! 
    369393      !! ** input   :   external netcdf files 
    370394      !! 
    371395      !!---------------------------------------------------------------------- 
    372       USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 
    373       USE wrk_nemo, ONLY: zriverdoc => wrk_2d_1, zriver => wrk_2d_2, zndepo => wrk_2d_3 
    374       USE wrk_nemo, ONLY: zcmask => wrk_3d_2 
    375       ! 
    376       INTEGER :: ji, jj, jk, jm 
    377       INTEGER :: numriv, numbath, numdep 
    378       REAL(wp) ::   zcoef 
    379       REAL(wp) ::   expide, denitide,zmaskt 
    380       ! 
    381       NAMELIST/nampissed/ ln_dustfer, ln_river, ln_ndepo, ln_sedinput, sedfeinput, dustsolub 
     396      ! 
     397      INTEGER  :: ji, jj, jk, jm 
     398      INTEGER  :: numdust, numriv, numiron, numdepo 
     399      INTEGER  :: ierr, ierr1, ierr2, ierr3 
     400      REAL(wp) :: zexpide, zdenitide, zmaskt 
     401      REAL(wp), DIMENSION(nbtimes) :: zsteps                 ! times records 
     402      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zdust, zndepo, zriverdic, zriverdoc, zcmask 
     403      ! 
     404      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files 
     405      TYPE(FLD_N) ::   sn_dust, sn_riverdoc, sn_riverdic, sn_ndepo, sn_ironsed        ! informations about the fields to be read 
     406      NAMELIST/nampissed/cn_dir, sn_dust, sn_riverdic, sn_riverdoc, sn_ndepo, sn_ironsed, & 
     407        &                ln_dust, ln_river, ln_ndepo, ln_ironsed,         & 
     408        &                sedfeinput, dustsolub, wdust, nitrfix, diazolight, concfediaz  
    382409      !!---------------------------------------------------------------------- 
    383  
    384       IF( ( wrk_in_use(2, 1,2,3) ) .OR. ( wrk_in_use(3, 2) ) ) THEN 
    385          CALL ctl_stop('p4z_sed_init: requested workspace arrays unavailable')  ;  RETURN 
    386       END IF 
    387       ! 
    388       REWIND( numnat )                     ! read numnat 
    389       READ  ( numnat, nampissed ) 
     410      ! 
     411      IF( nn_timing == 1 )  CALL timing_start('p4z_sed_init') 
     412      ! 
     413      !                                    ! number of seconds per year and per month 
     414      ryyss    = nyear_len(1) * rday 
     415      rmtss    = ryyss / raamo 
     416      r1_rday  = 1. / rday 
     417      r1_ryyss = 1. / ryyss 
     418      !                            !* set file information 
     419      cn_dir  = './'            ! directory in which the model is executed 
     420      ! ... default values (NB: frequency positive => hours, negative => months) 
     421      !                  !   file       ! frequency !  variable   ! time intep !  clim   ! 'yearly' or ! weights  ! rotation   ! 
     422      !                  !   name       !  (hours)  !   name      !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs      ! 
     423      sn_dust     = FLD_N( 'dust'       ,    -1     ,  'dust'     ,  .true.    , .true.  ,   'yearly'  , ''       , ''         ) 
     424      sn_riverdic = FLD_N( 'river'      ,   -12     ,  'riverdic' ,  .false.   , .true.  ,   'yearly'  , ''       , ''         ) 
     425      sn_riverdoc = FLD_N( 'river'      ,   -12     ,  'riverdoc' ,  .false.   , .true.  ,   'yearly'  , ''       , ''         ) 
     426      sn_ndepo    = FLD_N( 'ndeposition',   -12     ,  'ndep'     ,  .false.   , .true.  ,   'yearly'  , ''       , ''         ) 
     427      sn_ironsed  = FLD_N( 'ironsed'    ,   -12     ,  'bathy'    ,  .false.   , .true.  ,   'yearly'  , ''       , ''         ) 
     428 
     429      REWIND( numnatp )                     ! read numnatp 
     430      READ  ( numnatp, nampissed ) 
    390431 
    391432      IF(lwp) THEN 
    392433         WRITE(numout,*) ' ' 
    393          WRITE(numout,*) ' Namelist : nampissed ' 
     434         WRITE(numout,*) ' namelist : nampissed ' 
    394435         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~ ' 
    395          WRITE(numout,*) '    Dust input from the atmosphere           ln_dustfer  = ', ln_dustfer 
    396          WRITE(numout,*) '    River input of nutrients                 ln_river    = ', ln_river 
    397          WRITE(numout,*) '    Atmospheric deposition of N              ln_ndepo    = ', ln_ndepo 
    398          WRITE(numout,*) '    Fe input from sediments                  ln_sedinput = ', ln_sedinput 
    399          WRITE(numout,*) '    Coastal release of Iron                  sedfeinput  =', sedfeinput 
    400          WRITE(numout,*) '    Solubility of the dust                   dustsolub   =', dustsolub 
    401       ENDIF 
    402  
    403       ! Dust input from the atmosphere 
     436         WRITE(numout,*) '    dust input from the atmosphere           ln_dust     = ', ln_dust 
     437         WRITE(numout,*) '    river input of nutrients                 ln_river    = ', ln_river 
     438         WRITE(numout,*) '    atmospheric deposition of n              ln_ndepo    = ', ln_ndepo 
     439         WRITE(numout,*) '    fe input from sediments                  ln_sedinput = ', ln_ironsed 
     440         WRITE(numout,*) '    coastal release of iron                  sedfeinput  = ', sedfeinput 
     441         WRITE(numout,*) '    solubility of the dust                   dustsolub   = ', dustsolub 
     442         WRITE(numout,*) '    sinking speed of the dust                wdust       = ', wdust 
     443         WRITE(numout,*) '    nitrogen fixation rate                   nitrfix     = ', nitrfix 
     444         WRITE(numout,*) '    nitrogen fixation sensitivty to light    diazolight  = ', diazolight 
     445         WRITE(numout,*) '    fe half-saturation cste for diazotrophs  concfediaz  = ', concfediaz 
     446       END IF 
     447 
     448      IF( ln_dust .OR. ln_river .OR. ln_ndepo ) THEN 
     449          ll_sbc = .TRUE. 
     450      ELSE 
     451          ll_sbc = .FALSE. 
     452      ENDIF 
     453 
     454      ! dust input from the atmosphere 
    404455      ! ------------------------------ 
    405       IF( ln_dustfer ) THEN  
    406          IF(lwp) WRITE(numout,*) '    Initialize dust input from atmosphere ' 
     456      IF( ln_dust ) THEN  
     457         IF(lwp) WRITE(numout,*) '    initialize dust input from atmosphere ' 
    407458         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ' 
    408          CALL iom_open ( 'dust.orca.nc', numdust ) 
    409          DO jm = 1, jpmth 
    410             CALL iom_get( numdust, jpdom_data, 'dust', dustmo(:,:,jm), jm ) 
     459         ! 
     460         ALLOCATE( sf_dust(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst 
     461         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_apr structure' ) 
     462         ! 
     463         CALL fld_fill( sf_dust, (/ sn_dust /), cn_dir, 'p4z_sed_init', 'Iron from sediment ', 'nampissed' ) 
     464                                   ALLOCATE( sf_dust(1)%fnow(jpi,jpj,1)   ) 
     465         IF( sn_dust%ln_tint )     ALLOCATE( sf_dust(1)%fdta(jpi,jpj,1,2) ) 
     466         ! 
     467         ! Get total input dust ; need to compute total atmospheric supply of Si in a year 
     468         CALL iom_open (  TRIM( sn_dust%clname ) , numdust ) 
     469         CALL iom_gettime( numdust, zsteps, kntime=ntimes_dust)  ! get number of record in file 
     470         ALLOCATE( zdust(jpi,jpj,ntimes_dust) ) 
     471         DO jm = 1, ntimes_dust 
     472            CALL iom_get( numdust, jpdom_data, TRIM( sn_dust%clvar ), zdust(:,:,jm), jm ) 
    411473         END DO 
    412474         CALL iom_close( numdust ) 
     475         sumdepsi = 0.e0 
     476         DO jm = 1, ntimes_dust 
     477            sumdepsi = sumdepsi + glob_sum( zdust(:,:,jm) * e1e2t(:,:) * tmask(:,:,1) )  
     478         ENDDO 
     479         sumdepsi = sumdepsi * r1_ryyss * 8.8 * 0.075 / 28.1  
     480         DEALLOCATE( zdust) 
    413481      ELSE 
    414          dustmo(:,:,:) = 0.e0 
    415          dust(:,:) = 0.0 
    416       ENDIF 
    417  
    418       ! Nutrient input from rivers 
     482         dust(:,:) = 0._wp 
     483         sumdepsi  = 0._wp 
     484      END IF 
     485 
     486      ! nutrient input from rivers 
    419487      ! -------------------------- 
    420488      IF( ln_river ) THEN 
    421          IF(lwp) WRITE(numout,*) '    Initialize the nutrient input by rivers from river.orca.nc file' 
    422          IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' 
    423          CALL iom_open ( 'river.orca.nc', numriv ) 
    424          CALL iom_get  ( numriv, jpdom_data, 'riverdic', zriver   (:,:), jpyr ) 
    425          CALL iom_get  ( numriv, jpdom_data, 'riverdoc', zriverdoc(:,:), jpyr ) 
     489         ALLOCATE( sf_riverdic(1), STAT=ierr1 )           !* allocate and fill sf_sst (forcing structure) with sn_sst 
     490         ALLOCATE( sf_riverdoc(1), STAT=ierr2 )           !* allocate and fill sf_sst (forcing structure) with sn_sst 
     491         IF( ierr1 + ierr2 > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_apr structure' ) 
     492         ! 
     493         CALL fld_fill( sf_riverdic, (/ sn_riverdic /), cn_dir, 'p4z_sed_init', 'Input DOC from river ', 'nampissed' ) 
     494         CALL fld_fill( sf_riverdoc, (/ sn_riverdoc /), cn_dir, 'p4z_sed_init', 'Input DOC from river ', 'nampissed' ) 
     495                                   ALLOCATE( sf_riverdic(1)%fnow(jpi,jpj,1)   ) 
     496                                   ALLOCATE( sf_riverdoc(1)%fnow(jpi,jpj,1)   ) 
     497         IF( sn_riverdic%ln_tint ) ALLOCATE( sf_riverdic(1)%fdta(jpi,jpj,1,2) ) 
     498         IF( sn_riverdoc%ln_tint ) ALLOCATE( sf_riverdoc(1)%fdta(jpi,jpj,1,2) ) 
     499         ! Get total input rivers ; need to compute total river supply in a year 
     500         CALL iom_open ( TRIM( sn_riverdic%clname ), numriv ) 
     501         CALL iom_gettime( numriv, zsteps, kntime=ntimes_riv) 
     502         ALLOCATE( zriverdic(jpi,jpj,ntimes_riv) )   ;     ALLOCATE( zriverdoc(jpi,jpj,ntimes_riv) ) 
     503         DO jm = 1, ntimes_riv 
     504            CALL iom_get( numriv, jpdom_data, TRIM( sn_riverdic%clvar ), zriverdic(:,:,jm), jm ) 
     505            CALL iom_get( numriv, jpdom_data, TRIM( sn_riverdoc%clvar ), zriverdoc(:,:,jm), jm ) 
     506         END DO 
    426507         CALL iom_close( numriv ) 
     508         ! N/P and Si releases due to coastal rivers 
     509         ! ----------------------------------------- 
     510         rivpo4input = 0._wp  
     511         rivalkinput = 0._wp  
     512         DO jm = 1, ntimes_riv 
     513            rivpo4input = rivpo4input + glob_sum( ( zriverdic(:,:,jm) + zriverdoc(:,:,jm) ) * tmask(:,:,1) )  
     514            rivalkinput = rivalkinput + glob_sum(   zriverdic(:,:,jm)                       * tmask(:,:,1) )  
     515         END DO 
     516         rivpo4input = rivpo4input * 1E9 / 31.6_wp 
     517         rivalkinput = rivalkinput * 1E9 / 12._wp  
     518         DEALLOCATE( zriverdic)   ;    DEALLOCATE( zriverdoc)  
    427519      ELSE 
    428          zriver   (:,:) = 0.e0 
    429          zriverdoc(:,:) = 0.e0 
    430       endif 
    431  
    432       ! Nutrient input from dust 
     520         rivinp(:,:) = 0._wp 
     521         cotdep(:,:) = 0._wp 
     522         rivpo4input = 0._wp 
     523         rivalkinput = 0._wp 
     524      END IF  
     525 
     526      ! nutrient input from dust 
    433527      ! ------------------------ 
    434528      IF( ln_ndepo ) THEN 
    435          IF(lwp) WRITE(numout,*) '    Initialize the nutrient input by dust from ndeposition.orca.nc' 
     529         IF(lwp) WRITE(numout,*) '    initialize the nutrient input by dust from ndeposition.orca.nc' 
    436530         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' 
    437          CALL iom_open ( 'ndeposition.orca.nc', numdep ) 
    438          CALL iom_get  ( numdep, jpdom_data, 'ndep', zndepo(:,:), jpyr ) 
    439          CALL iom_close( numdep ) 
     531         ALLOCATE( sf_ndepo(1), STAT=ierr3 )           !* allocate and fill sf_sst (forcing structure) with sn_sst 
     532         IF( ierr3 > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_apr structure' ) 
     533         ! 
     534         CALL fld_fill( sf_ndepo, (/ sn_ndepo /), cn_dir, 'p4z_sed_init', 'Iron from sediment ', 'nampissed' ) 
     535                                   ALLOCATE( sf_ndepo(1)%fnow(jpi,jpj,1)   ) 
     536         IF( sn_ndepo%ln_tint )    ALLOCATE( sf_ndepo(1)%fdta(jpi,jpj,1,2) ) 
     537         ! 
     538         ! Get total input dust ; need to compute total atmospheric supply of N in a year 
     539         CALL iom_open ( TRIM( sn_ndepo%clname ), numdepo ) 
     540         CALL iom_gettime( numdepo, zsteps, kntime=ntimes_ndep) 
     541         ALLOCATE( zndepo(jpi,jpj,ntimes_ndep) ) 
     542         DO jm = 1, ntimes_ndep 
     543            CALL iom_get( numdepo, jpdom_data, TRIM( sn_ndepo%clvar ), zndepo(:,:,jm), jm ) 
     544         END DO 
     545         CALL iom_close( numdepo ) 
     546         nitdepinput = 0._wp 
     547         DO jm = 1, ntimes_ndep 
     548           nitdepinput = nitdepinput + glob_sum( zndepo(:,:,jm) * e1e2t(:,:) * tmask(:,:,1) )  
     549         ENDDO 
     550         nitdepinput = nitdepinput * 7.6 / 14E6  
     551         DEALLOCATE( zndepo) 
    440552      ELSE 
    441          zndepo(:,:) = 0.e0 
    442       ENDIF 
    443  
    444       ! Coastal and island masks 
     553         nitdep(:,:) = 0._wp 
     554         nitdepinput = 0._wp 
     555      ENDIF 
     556 
     557      ! coastal and island masks 
    445558      ! ------------------------ 
    446       IF( ln_sedinput ) THEN      
    447          IF(lwp) WRITE(numout,*) '    Computation of an island mask to enhance coastal supply of iron' 
     559      IF( ln_ironsed ) THEN      
     560         IF(lwp) WRITE(numout,*) '    computation of an island mask to enhance coastal supply of iron' 
    448561         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' 
    449          IF(lwp) WRITE(numout,*) '       from bathy.orca.nc file ' 
    450          CALL iom_open ( 'bathy.orca.nc', numbath ) 
    451          CALL iom_get  ( numbath, jpdom_data, 'bathy', zcmask(:,:,:), jpyr ) 
    452          CALL iom_close( numbath ) 
     562         CALL iom_open ( TRIM( sn_ironsed%clname ), numiron ) 
     563         ALLOCATE( zcmask(jpi,jpj,jpk) ) 
     564         CALL iom_get  ( numiron, jpdom_data, TRIM( sn_ironsed%clvar ), zcmask(:,:,:), 1 ) 
     565         CALL iom_close( numiron ) 
    453566         ! 
    454567         DO jk = 1, 5 
     
    459572                        &                       * tmask(ji,jj-1,jk) * tmask(ji,jj,jk+1) 
    460573                     IF( zmaskt == 0. )   zcmask(ji,jj,jk ) = MAX( 0.1, zcmask(ji,jj,jk) )  
    461                   ENDIF 
     574                  END IF 
    462575               END DO 
    463576            END DO 
    464577         END DO 
     578         CALL lbc_lnk( zcmask , 'T', 1. )      ! lateral boundary conditions on cmask   (sign unchanged) 
    465579         DO jk = 1, jpk 
    466580            DO jj = 1, jpj 
    467581               DO ji = 1, jpi 
    468                   expide   = MIN( 8.,( fsdept(ji,jj,jk) / 500. )**(-1.5) ) 
    469                   denitide = -0.9543 + 0.7662 * LOG( expide ) - 0.235 * LOG( expide )**2 
    470                   zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( denitide ) / 0.5 ) 
     582                  zexpide   = MIN( 8.,( fsdept(ji,jj,jk) / 500. )**(-1.5) ) 
     583                  zdenitide = -0.9543 + 0.7662 * LOG( zexpide ) - 0.235 * LOG( zexpide )**2 
     584                  zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( zdenitide ) / 0.5 ) 
    471585               END DO 
    472586            END DO 
    473587         END DO 
     588         ! Coastal supply of iron 
     589         ! ------------------------- 
     590         ironsed(:,:,jpk) = 0._wp 
     591         DO jk = 1, jpkm1 
     592            ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( fse3t(:,:,jk) * rday ) 
     593         END DO 
     594         DEALLOCATE( zcmask) 
    474595      ELSE 
    475          zcmask(:,:,:) = 0.e0 
    476       ENDIF 
    477  
    478       CALL lbc_lnk( zcmask , 'T', 1. )      ! Lateral boundary conditions on zcmask   (sign unchanged) 
    479  
    480  
    481       !                                    ! Number of seconds per year and per month 
    482       ryyss  = nyear_len(1) * rday 
    483       rmtss  = ryyss / raamo 
    484       rday1  = 1. / rday 
    485       ryyss1 = 1. / ryyss 
    486       !                                    ! ocean surface cell 
    487  
    488       ! total atmospheric supply of Si 
    489       ! ------------------------------ 
    490       sumdepsi = 0.e0 
    491       DO jm = 1, jpmth 
    492          zcoef = 1. / ( 12. * rmtss ) * 8.8 * 0.075 / 28.1         
    493          sumdepsi = sumdepsi + glob_sum( dustmo(:,:,jm) * e1e2t(:,:) ) * zcoef 
    494       ENDDO 
    495  
    496       ! N/P and Si releases due to coastal rivers 
    497       ! ----------------------------------------- 
    498       DO jj = 1, jpj 
    499          DO ji = 1, jpi 
    500             zcoef = ryyss * e1e2t(ji,jj)  * fse3t(ji,jj,1) * tmask(ji,jj,1)  
    501             cotdep(ji,jj) =  zriver(ji,jj)                  *1E9 / ( 12. * zcoef + rtrn ) 
    502             rivinp(ji,jj) = (zriver(ji,jj)+zriverdoc(ji,jj)) *1E9 / ( 31.6* zcoef + rtrn ) 
    503             nitdep(ji,jj) = 7.6 * zndepo(ji,jj)                  / ( 14E6*ryyss*fse3t(ji,jj,1) + rtrn ) 
    504          END DO 
    505       END DO 
    506       ! Lateral boundary conditions on ( cotdep, rivinp, nitdep )   (sign unchanged) 
    507       CALL lbc_lnk( cotdep , 'T', 1. )  ;  CALL lbc_lnk( rivinp , 'T', 1. )  ;  CALL lbc_lnk( nitdep , 'T', 1. ) 
    508  
    509       rivpo4input = glob_sum( rivinp(:,:) * cvol(:,:,1) ) * ryyss 
    510       rivalkinput = glob_sum( cotdep(:,:) * cvol(:,:,1) ) * ryyss 
    511       nitdepinput = glob_sum( nitdep(:,:) * cvol(:,:,1) ) * ryyss 
    512  
    513  
    514       ! Coastal supply of iron 
    515       ! ------------------------- 
    516       DO jk = 1, jpkm1 
    517          ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( fse3t(:,:,jk) * rday ) 
    518       END DO 
    519       CALL lbc_lnk( ironsed , 'T', 1. )      ! Lateral boundary conditions on ( ironsed )   (sign unchanged) 
    520  
    521       IF( ( wrk_not_released(2, 1,2,3) ) .OR. ( wrk_not_released(3, 2) ) )   & 
    522         &         CALL ctl_stop('p4z_sed_init: failed to release workspace arrays') 
    523  
     596         ironsed(:,:,:) = 0._wp 
     597      ENDIF 
     598      ! 
     599      IF( ll_sbc ) CALL p4z_sbc( nit000 )  
     600      ! 
     601      IF(lwp) THEN  
     602         WRITE(numout,*) 
     603         WRITE(numout,*) '    Total input of elements from river supply' 
     604         WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' 
     605         WRITE(numout,*) '    N Supply   : ', rivpo4input/7.6*1E3/1E12*14.,' TgN/yr' 
     606         WRITE(numout,*) '    Si Supply  : ', rivalkinput/6.*1E3/1E12*32.,' TgSi/yr' 
     607         WRITE(numout,*) '    Alk Supply : ', rivalkinput*1E3/1E12,' Teq/yr' 
     608         WRITE(numout,*) '    DIC Supply : ', rivpo4input*2.631*1E3*12./1E12,'TgC/yr' 
     609         WRITE(numout,*)  
     610         WRITE(numout,*) '    Total input of elements from atmospheric supply' 
     611         WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' 
     612         WRITE(numout,*) '    N Supply   : ', nitdepinput/7.6*1E3/1E12*14.,' TgN/yr' 
     613         WRITE(numout,*)  
     614      ENDIF 
     615      ! 
     616      IF( nn_timing == 1 )  CALL timing_stop('p4z_sed_init') 
     617      ! 
    524618   END SUBROUTINE p4z_sed_init 
    525619 
     
    529623      !!---------------------------------------------------------------------- 
    530624 
    531       ALLOCATE( dustmo(jpi,jpj,jpmth), dust(jpi,jpj)       ,     & 
    532         &       rivinp(jpi,jpj)      , cotdep(jpi,jpj)     ,     & 
    533         &       nitdep(jpi,jpj)      , ironsed(jpi,jpj,jpk), STAT=p4z_sed_alloc )   
     625      ALLOCATE( dust  (jpi,jpj), rivinp(jpi,jpj)     , cotdep(jpi,jpj),      & 
     626        &       nitdep(jpi,jpj), ironsed(jpi,jpj,jpk), STAT=p4z_sed_alloc )   
    534627 
    535628      IF( p4z_sed_alloc /= 0 ) CALL ctl_warn('p4z_sed_alloc : failed to allocate arrays.') 
Note: See TracChangeset for help on using the changeset viewer.