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 2823 for branches/2011/dev_r2787_PISCES_improvment/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zsed.F90 – NEMO

Ignore:
Timestamp:
2011-08-09T13:11:24+02:00 (13 years ago)
Author:
cetlod
Message:

Add new parameterisation in PISCES, see ticket #854

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2011/dev_r2787_PISCES_improvment/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zsed.F90

    r2774 r2823  
    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 
     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 
     29 
    2830 
    2931 
     
    3638 
    3739   !! * 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 
     40   LOGICAL  :: ln_dust     = .FALSE.    !: boolean for dust input from the atmosphere 
     41   LOGICAL  :: ln_river    = .FALSE.    !: boolean for river input of nutrients 
     42   LOGICAL  :: ln_ndepo    = .FALSE.    !: boolean for atmospheric deposition of N 
     43   LOGICAL  :: ln_ironsed  = .FALSE.    !: boolean for Fe input from sediments 
     44 
     45   REAL(wp) :: sedfeinput  = 1.E-9_wp   !: Coastal release of Iron 
     46   REAL(wp) :: dustsolub   = 0.014_wp   !: Solubility of the dust 
     47   REAL(wp) :: wdust       = 2.0_wp     !: Sinking speed of the dust  
     48   REAL(wp) :: nitrfix     = 1E-7_wp    !: Nitrogen fixation rate    
     49   REAL(wp) :: diazolight  = 50._wp     !: Nitrogen fixation sensitivty to light  
     50   REAL(wp) :: concfediaz  = 1.E-10_wp  !: Fe half-saturation Cste for diazotrophs  
     51 
    4552 
    4653   !! * Module variables 
    4754   REAL(wp) :: ryyss                  !: number of seconds per year  
    48    REAL(wp) :: ryyss1                 !: inverse of ryyss 
     55   REAL(wp) :: r1_ryyss                 !: inverse of ryyss 
    4956   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 
     57   REAL(wp) :: r1_rday                  !: inverse of rday 
     58   LOGICAL  :: ll_sbc 
     59 
     60   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_dust      ! structure of input dust 
     61   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_riverdic  ! structure of input riverdic 
     62   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_riverdoc  ! structure of input riverdoc 
     63   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_ndepo     ! structure of input nitrogen deposition 
     64   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_ironsed   ! structure of input iron from sediment 
     65 
     66   INTEGER , PARAMETER :: nbtimes = 365  !: maximum number of times record in a file 
     67   INTEGER  :: ntimes_dust, ntimes_riv, ntimes_ndep       ! number of time steps in a file 
     68 
    6069   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:) :: dust      !: dust fields 
    6170   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:) :: rivinp, cotdep    !: river input fields 
     
    8695      !! ** Method  : - ??? 
    8796      !!--------------------------------------------------------------------- 
    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 
     97      USE wrk_nemo, ONLY: wrk_in_USE, wrk_not_released 
     98      USE wrk_nemo, ONLY: zsidep   => wrk_2d_11 
     99      USE wrk_nemo, ONLY: zwork1   => wrk_2d_12, zwork2 => wrk_2d_13, zwork3 => wrk_2d_14 
    90100      USE wrk_nemo, ONLY: znitrpot => wrk_3d_2, zirondep => wrk_3d_3 
    91101      ! 
     
    96106      REAL(wp) ::   zrivalk, zrivsil, zrivpo4 
    97107#endif 
    98       REAL(wp) ::   zdenitot, znitrpottot, zlim, zfact 
    99       REAL(wp) ::   zwsbio3, zwsbio4, zwscal 
     108      REAL(wp) ::   zdenitot, znitrpottot, zlim, zfact, zfactcal 
     109      REAL(wp) ::   zsiloss, zcaloss, zwsbio3, zwsbio4, zwscal, zdep 
    100110      CHARACTER (len=25) :: charout 
    101111      !!--------------------------------------------------------------------- 
    102112 
    103       IF( ( wrk_in_use(2, 1,2,3) ) .OR. ( wrk_in_use(3, 2,3) ) ) THEN 
     113      IF( ( wrk_in_USE(2, 11,12,13,14) ) .OR. ( wrk_in_USE(3, 2,3) ) ) THEN 
    104114         CALL ctl_stop('p4z_sed: requested workspace arrays unavailable')  ;  RETURN 
    105115      END IF 
    106116 
    107       IF( jnt == 1  .AND.  ln_dustfer  )  CALL p4z_sbc( kt ) 
     117      IF( jnt == 1 .AND. ll_sbc ) CALL p4z_sbc( kt ) 
     118 
     119      zirondep(:,:,:) = 0.e0          ! Initialisation of variables USEd to compute deposition 
     120      zsidep  (:,:)   = 0.e0 
    108121 
    109122      ! Iron and Si deposition at the surface 
    110123      ! ------------------------------------- 
    111  
    112124      DO jj = 1, jpj 
    113125         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 ) 
     126            zdep  = rfact2 / fse3t(ji,jj,1) 
     127            zirondep(ji,jj,1) = ( dustsolub * dust(ji,jj) / ( 55.85 * rmtss ) + 3.e-10 * r1_ryyss ) * zdep 
     128            zsidep  (ji,jj)   = 8.8 * 0.075 * dust(ji,jj) * zdep / ( 28.1 * rmtss ) 
    117129         END DO 
    118130      END DO 
     
    120132      ! Iron solubilization of particles in the water column 
    121133      ! ---------------------------------------------------- 
    122  
    123134      DO jk = 2, jpkm1 
    124          zirondep(:,:,jk) = dust(:,:) / ( 10. * 55.85 * rmtss ) * rfact2 * 1.e-4 
     135         zirondep(:,:,jk) = dust(:,:) / ( wdust * 55.85 * rmtss ) * rfact2 * 1.e-4 * EXP( -fsdept(:,:,jk) / 1000. ) 
    125136      END DO 
    126137 
    127138      ! Add the external input of nutrients, carbon and alkalinity 
    128139      ! ---------------------------------------------------------- 
    129  
    130140      trn(:,:,1,jppo4) = trn(:,:,1,jppo4) + rivinp(:,:) * rfact2  
    131141      trn(:,:,1,jpno3) = trn(:,:,1,jpno3) + (rivinp(:,:) + nitdep(:,:)) * rfact2 
     
    139149      ! (dust, river and sediment mobilization) 
    140150      ! ------------------------------------------------------ 
    141  
    142151      DO jk = 1, jpkm1 
    143152         trn(:,:,jk,jpfer) = trn(:,:,jk,jpfer) + zirondep(:,:,jk) + ironsed(:,:,jk) * rfact2 
    144153      END DO 
    145  
    146154 
    147155#if ! defined key_sed 
     
    154162            ikt = mbkt(ji,jj)  
    155163# 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) 
     164            zwork1(ji,jj) = trn(ji,jj,ikt,jpdsi) * wscal (ji,jj,ikt) 
     165            zwork2(ji,jj) = trn(ji,jj,ikt,jppoc) * wsbio3(ji,jj,ikt) 
    158166# 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)  
     167            zwork1(ji,jj) = trn(ji,jj,ikt,jpdsi) * wsbio4(ji,jj,ikt) 
     168            zwork2(ji,jj) = trn(ji,jj,ikt,jpgoc) * wsbio4(ji,jj,ikt) + trn(ji,jj,ikt,jppoc) * wsbio3(ji,jj,ikt)  
    161169# 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 
    173 #endif 
    174  
    175       ! Then this loss is scaled at each bottom grid cell for 
     170            ! For calcite, burial efficiency is made a function of saturation 
     171            zfactcal      = MIN( excess(ji,jj,ikt), 0.2 ) 
     172            zfactcal      = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) ) 
     173            zwork3(ji,jj) = trn(ji,jj,ikt,jpcal) * wscal (ji,jj,ikt) * 2.e0 * zfactcal 
     174         END DO 
     175      END DO 
     176      zsumsedsi  = glob_sum( zwork1(:,:) * e1e2t(:,:) ) * r1_rday 
     177      zsumsedpo4 = glob_sum( zwork2(:,:) * e1e2t(:,:) ) * r1_rday 
     178      zsumsedcal = glob_sum( zwork3(:,:) * e1e2t(:,:) ) * r1_rday 
     179#endif 
     180 
     181      ! THEN this loss is scaled at each bottom grid cell for 
    176182      ! equilibrating the total budget of silica in the ocean. 
    177183      ! Thus, the amount of silica lost in the sediments equal 
    178184      ! the supply at the surface (dust+rivers) 
    179185      ! ------------------------------------------------------ 
     186#if ! defined key_sed 
     187      zrivsil =  1._wp - ( sumdepsi + rivalkinput * r1_ryyss / 6. ) / zsumsedsi  
     188      zrivpo4 =  1._wp - ( rivpo4input * r1_ryyss ) / zsumsedpo4  
     189#endif 
    180190 
    181191      DO jj = 1, jpj 
    182192         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) 
     193            ikt  = mbkt(ji,jj) 
     194            zdep = xstep / fse3t(ji,jj,ikt) 
     195            zwsbio4 = wsbio4(ji,jj,ikt) * zdep 
     196            zwscal  = wscal (ji,jj,ikt) * zdep 
     197# if defined key_kriest 
     198            zsiloss = trn(ji,jj,ikt,jpdsi) * zwsbio4 
     199# else 
     200            zsiloss = trn(ji,jj,ikt,jpdsi) * zwscal 
     201# endif 
     202            zcaloss = trn(ji,jj,ikt,jpcal) * zwscal 
    188203            ! 
    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  
     204            trn(ji,jj,ikt,jpdsi) = trn(ji,jj,ikt,jpdsi) - zsiloss 
     205            trn(ji,jj,ikt,jpcal) = trn(ji,jj,ikt,jpcal) - zcaloss 
    205206#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  
     207            trn(ji,jj,ikt,jpsil) = trn(ji,jj,ikt,jpsil) + zsiloss * zrivsil  
     208            zfactcal = MIN( excess(ji,jj,ikt), 0.2 ) 
     209            zfactcal = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) ) 
     210            zrivalk  =  1._wp - ( rivalkinput * r1_ryyss ) * zfactcal / zsumsedcal  
     211            trn(ji,jj,ikt,jptal) =  trn(ji,jj,ikt,jptal) + zcaloss * zrivalk * 2.0 
     212            trn(ji,jj,ikt,jpdic) =  trn(ji,jj,ikt,jpdic) + zcaloss * zrivalk 
     213#endif 
     214         END DO 
     215      END DO 
     216 
    209217      DO jj = 1, jpj 
    210218         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  
     219            ikt  = mbkt(ji,jj) 
     220            zdep = xstep / fse3t(ji,jj,ikt) 
     221            zwsbio4 = wsbio4(ji,jj,ikt) * zdep 
     222            zwsbio3 = wsbio3(ji,jj,ikt) * zdep 
     223# if ! defined key_kriest 
     224            trn(ji,jj,ikt,jpgoc) = trn(ji,jj,ikt,jpgoc) - trn(ji,jj,ikt,jpgoc) * zwsbio4 
     225            trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc) - trn(ji,jj,ikt,jppoc) * zwsbio3 
     226            trn(ji,jj,ikt,jpbfe) = trn(ji,jj,ikt,jpbfe) - trn(ji,jj,ikt,jpbfe) * zwsbio4 
     227            trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe) - trn(ji,jj,ikt,jpsfe) * zwsbio3 
     228#if ! defined key_sed 
     229            trn(ji,jj,ikt,jpdoc) = trn(ji,jj,ikt,jpdoc) & 
     230               &               + ( trn(ji,jj,ikt,jpgoc) * zwsbio4 + trn(ji,jj,ikt,jppoc) * zwsbio3 ) * zrivpo4 
     231#endif 
     232 
    221233# 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 
     234            trn(ji,jj,ikt,jpnum) = trn(ji,jj,ikt,jpnum) - trn(ji,jj,ikt,jpnum) * zwsbio4 
     235            trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc) - trn(ji,jj,ikt,jppoc) * zwsbio3 
     236            trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe) - trn(ji,jj,ikt,jpsfe) * zwsbio3 
     237#if ! defined key_sed 
     238            trn(ji,jj,ikt,jpdoc) = trn(ji,jj,ikt,jpdoc) & 
     239               &               + ( trn(ji,jj,ikt,jpnum) * zwsbio4 + trn(ji,jj,ikt,jppoc) * zwsbio3 ) * zrivpo4 
     240#endif 
     241 
    225242# endif 
    226243         END DO 
    227244      END DO 
    228 # endif 
     245 
    229246 
    230247      ! Nitrogen fixation (simple parameterization). The total gain 
     
    233250      ! ------------------------------------------------------------- 
    234251 
    235       zdenitot = glob_sum( denitr(:,:,:)  * cvol(:,:,:) * xnegtr(:,:,:) ) * rdenit 
     252      zdenitot = glob_sum(  ( denitr(:,:,:) * rdenit + denitnh4(:,:,:) * rdenita ) * cvol(:,:,:) * xnegtr(:,:,:) ) 
    236253 
    237254      ! Potential nitrogen fixation dependant on temperature and iron 
     
    246263               zlim = ( 1.- xnanono3(ji,jj,jk) - xnanonh4(ji,jj,jk) ) 
    247264               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.) ) 
     265#if defined key_degrad 
     266               zfact = zlim * rfact2 * facvol(ji,jj,jk) 
     267#else 
     268               zfact = zlim * rfact2  
     269#endif 
     270               znitrpot(ji,jj,jk) =  MAX( 0.e0, ( 0.6 * tgfunc(ji,jj,jk) - 2.15 ) * r1_rday )   & 
     271                 &                 *  zfact * trn(ji,jj,jk,jpfer) / ( concfediaz + trn(ji,jj,jk,jpfer) ) & 
     272                 &                 * ( 1.- EXP( -etot(ji,jj,jk) / diazolight ) ) 
    254273            END DO 
    255274         END DO  
     
    260279      ! Nitrogen change due to nitrogen fixation 
    261280      ! ---------------------------------------- 
    262  
    263281      DO jk = 1, jpk 
    264282         DO jj = 1, jpj 
    265283            DO ji = 1, jpi 
    266                zfact = znitrpot(ji,jj,jk) * 1.e-7 
     284               zfact = znitrpot(ji,jj,jk) * nitrfix 
    267285               trn(ji,jj,jk,jpnh4) = trn(ji,jj,jk,jpnh4) + zfact 
     286               trn(ji,jj,jk,jptal) = trn(ji,jj,jk,jptal) + rno3 * zfact 
    268287               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 
     288               trn(ji,jj,jk,jppo4) = trn(ji,jj,jk,jppo4) + 30. / 46. * zfact 
     289           !    trn(ji,jj,jk,jppo4) = trn(ji,jj,jk,jppo4) + zfact 
     290           END DO 
     291         END DO  
    272292      END DO 
    273293 
    274294#if defined key_diatrc 
    275295      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) 
     296#if defined key_iomput 
     297      zwork1(:,:)  =  ( zirondep(:,:,1) + ironsed(:,:,1) * rfact2 ) * zfact * fse3t(:,:,1) * tmask(:,:,1)  
     298      zwork2(:,:)  =    znitrpot(:,:,1) * nitrif                    * zfact * fse3t(:,:,1) * tmask(:,:,1) 
    282299      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 
     300          CALL iom_put( "Irondep", zwork1  )  ! surface downward net flux of iron 
     301          CALL iom_put( "Nfix"   , zwork2 )  ! nitrogen fixation at surface 
     302      ENDIF 
     303#else 
     304      trc2d(:,:,jp_pcs0_2d + 11) = zirondep(:,:,1)          * zfact * fse3t(:,:,1) * tmask(:,:,1) 
     305      trc2d(:,:,jp_pcs0_2d + 12) = znitrpot(:,:,1) * nitrif * zfact * fse3t(:,:,1) * tmask(:,:,1) 
    287306#endif 
    288307      ! 
    289        IF(ln_ctl)   THEN  ! print mean trends (used for debugging) 
    290          WRITE(charout, FMT="('sed ')") 
     308      IF(ln_ctl) THEN  ! print mean trends (USEd for debugging) 
     309         WRITE(charout, fmt="('sed ')") 
    291310         CALL prt_ctl_trc_info(charout) 
    292311         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) ) )   & 
     312      ENDIF 
     313 
     314      IF( ( wrk_not_released(2, 11,12,13,14) ) .OR. ( wrk_not_released(3, 2,3) ) )   & 
    296315        &         CALL ctl_stop('p4z_sed: failed to release workspace arrays') 
    297316 
     
    299318 
    300319   SUBROUTINE p4z_sbc( kt ) 
    301  
    302320      !!---------------------------------------------------------------------- 
    303       !!                  ***  ROUTINE p4z_sbc  *** 
    304       !! 
    305       !! ** Purpose :   Read and interpolate the external sources of  
     321      !!                  ***  routine p4z_sbc  *** 
     322      !! 
     323      !! ** purpose :   read and interpolate the external sources of  
    306324      !!                nutrients 
    307325      !! 
    308       !! ** Method  :   Read the files and interpolate the appropriate variables 
     326      !! ** method  :   read the files and interpolate the appropriate variables 
    309327      !! 
    310328      !! ** input   :   external netcdf files 
     
    314332      INTEGER, INTENT( in  ) ::   kt   ! ocean time step 
    315333 
    316       !! * Local declarations 
    317       INTEGER :: imois, i15, iman  
    318       REAL(wp) :: zxy 
     334      !! * local declarations 
     335      INTEGER  :: ji,jj  
     336      REAL(wp) :: zcoef 
    319337 
    320338      !!--------------------------------------------------------------------- 
    321339 
    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  
     340      ! Compute dust at nit000 or only if there is more than 1 time record in dust file 
     341      IF( ln_dust ) THEN 
     342         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_dust > 1 ) ) THEN 
     343            CALL fld_read( kt, 1, sf_dust ) 
     344            dust(:,:) = sf_dust(1)%fnow(:,:,1) 
     345         ENDIF 
     346      ENDIF 
     347 
     348      ! N/P and Si releases due to coastal rivers 
     349      ! Compute river at nit000 or only if there is more than 1 time record in river file 
     350      ! ----------------------------------------- 
     351      IF( ln_river ) THEN 
     352         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_riv > 1 ) ) THEN 
     353            CALL fld_read( kt, 1, sf_riverdic ) 
     354            CALL fld_read( kt, 1, sf_riverdoc ) 
     355            DO jj = 1, jpj 
     356               DO ji = 1, jpi 
     357                  zcoef = ryyss * cvol(ji,jj,1)  
     358                  cotdep(ji,jj) =   sf_riverdic(1)%fnow(ji,jj,1)                                  * 1E9 / ( 12. * zcoef + rtrn ) 
     359                  rivinp(ji,jj) = ( sf_riverdic(1)%fnow(ji,jj,1) + sf_riverdoc(1)%fnow(ji,jj,1) ) * 1E9 / ( 31.6* zcoef + rtrn ) 
     360               END DO 
     361            END DO 
     362         ENDIF 
     363      ENDIF 
     364 
     365      ! Compute N deposition at nit000 or only if there is more than 1 time record in N deposition file 
     366      IF( ln_ndepo ) THEN 
     367         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_ndep > 1 ) ) THEN 
     368            CALL fld_read( kt, 1, sf_ndepo ) 
     369            DO jj = 1, jpj 
     370               DO ji = 1, jpi 
     371                  nitdep(ji,jj) = 7.6 * sf_ndepo(1)%fnow(ji,jj,1) / ( 14E6 * ryyss * fse3t(ji,jj,1) + rtrn ) 
     372               END DO 
     373            END DO 
     374         ENDIF 
     375      ENDIF 
     376      ! 
    356377   END SUBROUTINE p4z_sbc 
    357378 
     
    360381 
    361382      !!---------------------------------------------------------------------- 
    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) 
     383      !!                  ***  routine p4z_sed_init  *** 
     384      !! 
     385      !! ** purpose :   initialization of the external sources of nutrients 
     386      !! 
     387      !! ** method  :   read the files and compute the budget 
     388      !!                called at the first timestep (nit000) 
    368389      !! 
    369390      !! ** input   :   external netcdf files 
    370391      !! 
    371392      !!---------------------------------------------------------------------- 
    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 
    375393      ! 
    376       INTEGER :: ji, jj, jk, jm 
    377       INTEGER :: numriv, numbath, numdep 
    378       REAL(wp) ::   zcoef 
    379       REAL(wp) ::   expide, denitide,zmaskt 
     394      INTEGER  :: ji, jj, jk, jm 
     395      INTEGER  :: numdust, numriv, numiron, numdepo 
     396      INTEGER  :: ierr, ierr1, ierr2, ierr3 
     397      REAL(wp) :: zexpide, zdenitide, zmaskt 
     398      REAL(wp), DIMENSION(nbtimes) :: zsteps                 ! times records 
     399      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zdust, zndepo, zriverdic, zriverdoc, zcmask 
    380400      ! 
    381       NAMELIST/nampissed/ ln_dustfer, ln_river, ln_ndepo, ln_sedinput, sedfeinput, dustsolub 
     401      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files 
     402      TYPE(FLD_N) ::   sn_dust, sn_riverdoc, sn_riverdic, sn_ndepo, sn_ironsed        ! informations about the fields to be read 
     403      NAMELIST/nampissed/cn_dir, sn_dust, sn_riverdic, sn_riverdoc, sn_ndepo, sn_ironsed, & 
     404        &                ln_dust, ln_river, ln_ndepo, ln_ironsed,         & 
     405        &                sedfeinput, dustsolub, wdust, nitrfix, diazolight, concfediaz  
    382406      !!---------------------------------------------------------------------- 
    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 ) 
     407      !                                    ! number of seconds per year and per month 
     408      ryyss    = nyear_len(1) * rday 
     409      rmtss    = ryyss / raamo 
     410      r1_rday  = 1. / rday 
     411      r1_ryyss = 1. / ryyss 
     412      !                            !* set file information 
     413      cn_dir  = './'            ! directory in which the model is executed 
     414      ! ... default values (NB: frequency positive => hours, negative => months) 
     415      !                  !   file       ! frequency !  variable   ! time intep !  clim   ! 'yearly' or ! weights  ! rotation   ! 
     416      !                  !   name       !  (hours)  !   name      !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs      ! 
     417      sn_dust     = FLD_N( 'dust'       ,    -1     ,  'dust'     ,  .true.    , .true.  ,   'yearly'  , ''       , ''         ) 
     418      sn_riverdic = FLD_N( 'river'      ,   -12     ,  'riverdic' ,  .false.   , .true.  ,   'yearly'  , ''       , ''         ) 
     419      sn_riverdoc = FLD_N( 'river'      ,   -12     ,  'riverdoc' ,  .false.   , .true.  ,   'yearly'  , ''       , ''         ) 
     420      sn_ndepo    = FLD_N( 'ndeposition',   -12     ,  'ndep'     ,  .false.   , .true.  ,   'yearly'  , ''       , ''         ) 
     421      sn_ironsed  = FLD_N( 'ironsed'    ,   -12     ,  'bathy'    ,  .false.   , .true.  ,   'yearly'  , ''       , ''         ) 
     422 
     423 
     424      REWIND( numnatp )                     ! read numnat 
     425      READ  ( numnatp, nampissed ) 
    390426 
    391427      IF(lwp) THEN 
    392428         WRITE(numout,*) ' ' 
    393          WRITE(numout,*) ' Namelist : nampissed ' 
     429         WRITE(numout,*) ' namelist : nampissed ' 
    394430         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 
     431         WRITE(numout,*) '    dust input from the atmosphere           ln_dust     = ', ln_dust 
     432         WRITE(numout,*) '    river input of nutrients                 ln_river    = ', ln_river 
     433         WRITE(numout,*) '    atmospheric deposition of n              ln_ndepo    = ', ln_ndepo 
     434         WRITE(numout,*) '    fe input from sediments                  ln_sedinput = ', ln_ironsed 
     435         WRITE(numout,*) '    coastal release of iron                  sedfeinput  = ', sedfeinput 
     436         WRITE(numout,*) '    solubility of the dust                   dustsolub   = ', dustsolub 
     437         WRITE(numout,*) '    sinking speed of the dust                wdust       = ', wdust 
     438         WRITE(numout,*) '    nitrogen fixation rate                   nitrfix     = ', nitrfix 
     439         WRITE(numout,*) '    nitrogen fixation sensitivty to light    diazolight  = ', diazolight 
     440         WRITE(numout,*) '    fe half-saturation cste for diazotrophs  concfediaz  = ', concfediaz 
     441       END IF 
     442 
     443      IF( ln_dust .OR. ln_river .OR. ln_ndepo ) THEN 
     444          ll_sbc = .TRUE. 
     445      ELSE 
     446          ll_sbc = .FALSE. 
     447      ENDIF 
     448 
     449      ! dust input from the atmosphere 
    404450      ! ------------------------------ 
    405       IF( ln_dustfer ) THEN  
    406          IF(lwp) WRITE(numout,*) '    Initialize dust input from atmosphere ' 
     451      IF( ln_dust ) THEN  
     452         IF(lwp) WRITE(numout,*) '    initialize dust input from atmosphere ' 
    407453         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 ) 
     454         ! 
     455         ALLOCATE( sf_dust(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst 
     456         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_apr structure' ) 
     457         ! 
     458         CALL fld_fill( sf_dust, (/ sn_dust /), cn_dir, 'p4z_sed_init', 'Iron from sediment ', 'nampissed' ) 
     459                                   ALLOCATE( sf_dust(1)%fnow(jpi,jpj,1)   ) 
     460         IF( sn_dust%ln_tint )     ALLOCATE( sf_dust(1)%fdta(jpi,jpj,1,2) ) 
     461         ! 
     462         ! Get total input dust ; need to compute total atmospheric supply of Si in a year 
     463         CALL iom_open (  TRIM( sn_dust%clname ) , numdust ) 
     464         CALL iom_gettime( numdust, zsteps, kntime=ntimes_dust)  ! get number of record in file 
     465         ALLOCATE( zdust(jpi,jpj,ntimes_dust) ) 
     466         DO jm = 1, ntimes_dust 
     467            CALL iom_get( numdust, jpdom_data, TRIM( sn_dust%clvar ), zdust(:,:,jm), jm ) 
    411468         END DO 
    412469         CALL iom_close( numdust ) 
     470         sumdepsi = 0.e0 
     471         DO jm = 1, ntimes_dust 
     472            sumdepsi = sumdepsi + glob_sum( zdust(:,:,jm) * e1e2t(:,:) * tmask(:,:,1) )  
     473         ENDDO 
     474         sumdepsi = sumdepsi * r1_ryyss * 8.8 * 0.075 / 28.1  
     475         DEALLOCATE( zdust) 
    413476      ELSE 
    414          dustmo(:,:,:) = 0.e0 
    415          dust(:,:) = 0.0 
    416       ENDIF 
    417  
    418       ! Nutrient input from rivers 
     477         dust(:,:) = 0._wp 
     478         sumdepsi  = 0._wp 
     479      END IF 
     480 
     481      ! nutrient input from rivers 
    419482      ! -------------------------- 
    420483      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 ) 
     484         ALLOCATE( sf_riverdic(1), STAT=ierr1 )           !* allocate and fill sf_sst (forcing structure) with sn_sst 
     485         ALLOCATE( sf_riverdoc(1), STAT=ierr2 )           !* allocate and fill sf_sst (forcing structure) with sn_sst 
     486         IF( ierr1 + ierr2 > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_apr structure' ) 
     487         ! 
     488         CALL fld_fill( sf_riverdic, (/ sn_riverdic /), cn_dir, 'p4z_sed_init', 'Input DOC from river ', 'nampissed' ) 
     489         CALL fld_fill( sf_riverdoc, (/ sn_riverdoc /), cn_dir, 'p4z_sed_init', 'Input DOC from river ', 'nampissed' ) 
     490                                   ALLOCATE( sf_riverdic(1)%fnow(jpi,jpj,1)   ) 
     491                                   ALLOCATE( sf_riverdoc(1)%fnow(jpi,jpj,1)   ) 
     492         IF( sn_riverdic%ln_tint ) ALLOCATE( sf_riverdic(1)%fdta(jpi,jpj,1,2) ) 
     493         IF( sn_riverdoc%ln_tint ) ALLOCATE( sf_riverdoc(1)%fdta(jpi,jpj,1,2) ) 
     494         ! Get total input rivers ; need to compute total river supply in a year 
     495         CALL iom_open ( TRIM( sn_riverdic%clname ), numriv ) 
     496         CALL iom_gettime( numriv, zsteps, kntime=ntimes_riv) 
     497         ALLOCATE( zriverdic(jpi,jpj,ntimes_riv) )   ;     ALLOCATE( zriverdoc(jpi,jpj,ntimes_riv) ) 
     498         DO jm = 1, ntimes_riv 
     499            CALL iom_get( numriv, jpdom_data, TRIM( sn_riverdic%clvar ), zriverdic(:,:,jm), jm ) 
     500            CALL iom_get( numriv, jpdom_data, TRIM( sn_riverdoc%clvar ), zriverdoc(:,:,jm), jm ) 
     501         END DO 
    426502         CALL iom_close( numriv ) 
     503         ! N/P and Si releases due to coastal rivers 
     504         ! ----------------------------------------- 
     505         rivpo4input = 0._wp  
     506         rivalkinput = 0._wp  
     507         DO jm = 1, ntimes_riv 
     508            rivpo4input = rivpo4input + glob_sum( ( zriverdic(:,:,jm) + zriverdoc(:,:,jm) ) * tmask(:,:,1) )  
     509            rivalkinput = rivalkinput + glob_sum(   zriverdic(:,:,jm)                       * tmask(:,:,1) )  
     510         END DO 
     511         rivpo4input = rivpo4input * 1E9 / 31.6_wp 
     512         rivalkinput = rivalkinput * 1E9 / 12._wp  
     513         DEALLOCATE( zriverdic)   ;    DEALLOCATE( zriverdoc)  
    427514      ELSE 
    428          zriver   (:,:) = 0.e0 
    429          zriverdoc(:,:) = 0.e0 
    430       endif 
    431  
    432       ! Nutrient input from dust 
     515         rivinp(:,:) = 0._wp 
     516         cotdep(:,:) = 0._wp 
     517         rivpo4input = 0._wp 
     518         rivalkinput = 0._wp 
     519      END IF  
     520 
     521      ! nutrient input from dust 
    433522      ! ------------------------ 
    434523      IF( ln_ndepo ) THEN 
    435          IF(lwp) WRITE(numout,*) '    Initialize the nutrient input by dust from ndeposition.orca.nc' 
     524         IF(lwp) WRITE(numout,*) '    initialize the nutrient input by dust from ndeposition.orca.nc' 
    436525         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 ) 
     526         ALLOCATE( sf_ndepo(1), STAT=ierr3 )           !* allocate and fill sf_sst (forcing structure) with sn_sst 
     527         IF( ierr3 > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_apr structure' ) 
     528         ! 
     529         CALL fld_fill( sf_ndepo, (/ sn_ndepo /), cn_dir, 'p4z_sed_init', 'Iron from sediment ', 'nampissed' ) 
     530                                   ALLOCATE( sf_ndepo(1)%fnow(jpi,jpj,1)   ) 
     531         IF( sn_ndepo%ln_tint )    ALLOCATE( sf_ndepo(1)%fdta(jpi,jpj,1,2) ) 
     532         ! 
     533         ! Get total input dust ; need to compute total atmospheric supply of N in a year 
     534         CALL iom_open ( TRIM( sn_ndepo%clname ), numdepo ) 
     535         CALL iom_gettime( numdepo, zsteps, kntime=ntimes_ndep) 
     536         ALLOCATE( zndepo(jpi,jpj,ntimes_ndep) ) 
     537         DO jm = 1, ntimes_ndep 
     538            CALL iom_get( numdepo, jpdom_data, TRIM( sn_ndepo%clvar ), zndepo(:,:,jm), jm ) 
     539         END DO 
     540         CALL iom_close( numdepo ) 
     541         nitdepinput = 0._wp 
     542         DO jm = 1, ntimes_ndep 
     543           nitdepinput = nitdepinput + glob_sum( zndepo(:,:,jm) * e1e2t(:,:) * tmask(:,:,1) )  
     544         ENDDO 
     545         nitdepinput = nitdepinput * 7.6 / 14E6  
     546         DEALLOCATE( zndepo) 
    440547      ELSE 
    441          zndepo(:,:) = 0.e0 
    442       ENDIF 
    443  
    444       ! Coastal and island masks 
     548         nitdep(:,:) = 0._wp 
     549         nitdepinput = 0._wp 
     550      ENDIF 
     551 
     552      ! coastal and island masks 
    445553      ! ------------------------ 
    446       IF( ln_sedinput ) THEN      
    447          IF(lwp) WRITE(numout,*) '    Computation of an island mask to enhance coastal supply of iron' 
     554      IF( ln_ironsed ) THEN      
     555         IF(lwp) WRITE(numout,*) '    computation of an island mask to enhance coastal supply of iron' 
    448556         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 ) 
     557         CALL iom_open ( TRIM( sn_ironsed%clname ), numiron ) 
     558         ALLOCATE( zcmask(jpi,jpj,jpk) ) 
     559         CALL iom_get  ( numiron, jpdom_data, TRIM( sn_ironsed%clvar ), zcmask(:,:,:), 1 ) 
     560         CALL iom_close( numiron ) 
    453561         ! 
    454562         DO jk = 1, 5 
     
    459567                        &                       * tmask(ji,jj-1,jk) * tmask(ji,jj,jk+1) 
    460568                     IF( zmaskt == 0. )   zcmask(ji,jj,jk ) = MAX( 0.1, zcmask(ji,jj,jk) )  
    461                   ENDIF 
     569                  END IF 
    462570               END DO 
    463571            END DO 
    464572         END DO 
     573         CALL lbc_lnk( zcmask , 'T', 1. )      ! lateral boundary conditions on cmask   (sign unchanged) 
    465574         DO jk = 1, jpk 
    466575            DO jj = 1, jpj 
    467576               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 ) 
     577                  zexpide   = MIN( 8.,( fsdept(ji,jj,jk) / 500. )**(-1.5) ) 
     578                  zdenitide = -0.9543 + 0.7662 * LOG( zexpide ) - 0.235 * LOG( zexpide )**2 
     579                  zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( zdenitide ) / 0.5 ) 
    471580               END DO 
    472581            END DO 
    473582         END DO 
     583         ! Coastal supply of iron 
     584         ! ------------------------- 
     585         ironsed(:,:,jpk) = 0._wp 
     586         DO jk = 1, jpkm1 
     587            ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( fse3t(:,:,jk) * rday ) 
     588         END DO 
     589         DEALLOCATE( zcmask) 
    474590      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  
     591         ironsed(:,:,:) = 0._wp 
     592      ENDIF 
     593      ! 
     594      IF(lwp) THEN  
     595         WRITE(numout,*) 
     596         WRITE(numout,*) '    Total input of elements from river supply' 
     597         WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' 
     598         WRITE(numout,*) '    N Supply   : ', rivpo4input/7.6*1E3/1E12*14.,' TgN/yr' 
     599         WRITE(numout,*) '    Si Supply  : ', rivalkinput/6.*1E3/1E12*32.,' TgSi/yr' 
     600         WRITE(numout,*) '    Alk Supply : ', rivalkinput*1E3/1E12,' Teq/yr' 
     601         WRITE(numout,*) '    DIC Supply : ', rivpo4input*2.631*1E3*12./1E12,'TgC/yr' 
     602         WRITE(numout,*)  
     603         WRITE(numout,*) '    Total input of elements from atmospheric supply' 
     604         WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' 
     605         WRITE(numout,*) '    N Supply   : ', nitdepinput/7.6*1E3/1E12*14.,' TgN/yr' 
     606         WRITE(numout,*)  
     607      ENDIF 
     608       ! 
    524609   END SUBROUTINE p4z_sed_init 
    525610 
     
    529614      !!---------------------------------------------------------------------- 
    530615 
    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 )   
     616      ALLOCATE( dust  (jpi,jpj), rivinp(jpi,jpj)     , cotdep(jpi,jpj),      & 
     617        &       nitdep(jpi,jpj), ironsed(jpi,jpj,jpk), STAT=p4z_sed_alloc )   
    534618 
    535619      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.