Ignore:
Timestamp:
2018-11-21T11:25:53+01:00 (2 years ago)
Author:
smasson
Message:

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10344, see #2133

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/SED/sedini.F90

    r5215 r10345  
    44   !! Sediment : define sediment variables 
    55   !!===================================================================== 
    6 #if defined key_sed 
    76 
    87   !!---------------------------------------------------------------------- 
     
    1110   !! * Modules used 
    1211   USE sed     ! sediment global variable 
    13    USE seddta 
    14    USE sedrst 
    15    USE sedco3 
    16    USE sedchem 
     12   USE sed_oce 
    1713   USE sedarr 
     14   USE sedadv 
     15   USE trc_oce, ONLY : nn_dttrc 
     16   USE trcdmp_sed 
     17   USE trcdta 
    1818   USE iom 
     19   USE lib_mpp         ! distribued memory computing library 
    1920 
    2021 
     
    2425   !! Module variables 
    2526   REAL(wp)    ::  & 
    26       sisat   =  800.  ,  &  !: saturation for si    [ mol.l-1] 
    27       claysat =    0.  ,  &  !: saturation for clay  [ mol.l-1] 
     27      sedzmin = 0.3    ,  &  !: Minimum vertical spacing 
     28      sedhmax = 10.0   ,  &  !: Maximum depth of the sediment 
     29      sedkth  = 5.0    ,  &  !: Default parameters 
     30      sedacr  = 3.0          !: Default parameters 
     31       
     32   REAL(wp)    ::  & 
     33      porsurf =  0.95  ,  &  !: Porosity at the surface 
     34      porinf  =  0.75  ,  &  !: Porosity at infinite depth 
     35      rhox    =  2.0         !: Vertical length scale of porosity variation  
     36 
     37   REAL(wp)    ::  & 
    2838      rcopal  =   40.  ,  &  !: reactivity for si    [l.mol-1.an-1] 
    29       rcclay  =    0.  ,  &  !: reactivity for clay  [l.mol-1.an-1] 
    3039      dcoef   =  8.e-6       !: diffusion coefficient (*por)   [cm**2/s] 
    3140 
     41   REAL(wp), PUBLIC    ::  & 
     42      redO2    =  172.  ,  &  !: Redfield coef for Oxygen 
     43      redNo3   =   16.  ,  &  !: Redfield coef for Nitrate 
     44      redPo4   =    1.  ,  &  !: Redfield coef for Phosphate 
     45      redC     =  122.  ,  &  !: Redfield coef for Carbon 
     46      redfep   =  0.175 ,  &  !: Ratio for iron bound phosphorus 
     47      rcorgl   =   50.  ,  &  !: reactivity for POC/O2 [l.mol-1.an-1]     
     48      rcorgs   =   0.5  ,  &  !: reactivity of the semi-labile component 
     49      rcorgr   =   1E-4 ,  &  !: reactivity of the refractory component 
     50      rcnh4    =   10E6 ,  &  !: reactivity for O2/NH4 [l.mol-1.an-1]   
     51      rch2s    =   1.E5 ,  &  !: reactivity for O2/ODU [l.mol-1.an-1]  
     52      rcfe2    =   5.E8 ,  &  !: reactivity for O2/Fe2+ [l.mol-1.an-1] 
     53      rcfeh2s  =   1.E4 ,  &  !: Reactivity for FEOH/H2S [l.mol-1.an-1] 
     54      rcfes    =   1.E5 ,  &  !: Reactivity for FE2+/H2S [l.mol-1.an-1] 
     55      rcfeso   =   3.E5 ,  &  !: Reactivity for FES/O2 [l.mol-1.an-1] 
     56      xksedo2  =   5E-6 ,  &  !: half-sturation constant for oxic remin. 
     57      xksedno3 =   5E-6 ,  &  !: half-saturation constant for denitrification 
     58      xksedfeo =   0.6  ,  & !: half-saturation constant for iron remin 
     59      xksedso4 =   2E-3       !: half-saturation constant for SO4 remin 
     60 
    3261   REAL(wp)    ::  & 
    33       redO2   =  172.  ,  &  !: Redfield coef for Oxygene 
    34       redNo3  =   16.  ,  &  !: Redfield coef for Nitrates 
    35       redPo4  =    1.  ,  &  !: Redfield coef for Phosphates 
    36       redC    =  122.  ,  &  !: Redfield coef for Carbone 
    37       redDnit =  97.6  ,  &  !: Redfield coef for denitrification     
    38       rcorg   =   50.  ,  &  !: reactivity for POC/O2 [l.mol-1.an-1]     
    39       o2seuil =    1.  ,  &  !: threshold O2 concentration for [mol.l-1]      
    40       rcorgN  =   50.       !: reactivity for POC/No3 [l.mol-1.an-1] 
    41  
    42    REAL(wp)    ::  & 
    43       rccal   = 1000.        !: reactivity for calcite         [l.mol-1.an-1] 
    44  
    45    REAL(wp)    ::  & 
    46       dbiot   = 15.          !: coefficient for bioturbation    [cm**2.(1000an-1)] 
    47  
    48    LOGICAL     :: & 
    49       ln_rst_sed = .TRUE.    !: initialisation from a restart file or not 
    50  
     62      rccal   = 1000.,      & !: reactivity for calcite         [l.mol-1.an-1] 
     63      rcligc  = 1.E-4         !: L/C ratio in POC 
     64 
     65   REAL(wp), PUBLIC    ::  dbiot   = 15. , &  !: coefficient for bioturbation    [cm**2.(n-1)] 
     66      dbtbzsc =  10.0  ,    &  !: Vertical scale of variation. If no variation, mixed layer in the sed [cm] 
     67      xirrzsc = 2.0            !: Vertical scale of irrigation variation. 
    5168   REAL(wp)    ::  & 
    5269      ryear = 365. * 24. * 3600. !:  1 year converted in second 
     70 
     71   REAL(wp), DIMENSION(jpwat), PUBLIC  :: diff1 
     72   DATA diff1/4.59E-6, 1.104E-5, 4.81E-6 , 9.78E-6, 3.58E-6, 4.01E-6, 9.8E-6, 9.73E-6, 5.0E-6, 3.31E-6 / 
     73 
     74   REAL(wp), DIMENSION(jpwat), PUBLIC  :: diff2 
     75   DATA diff2/1.74E-7, 4.47E-7, 2.51E-7, 3.89E-7, 1.77E-7, 2.5E-7, 3.89E-7, 3.06E-7, 2.5E-7, 1.5E-7 / 
     76 
     77 
    5378 
    5479   !! *  Routine accessibility 
     
    76101      !!        !  06-07  (C. Ethe)  Re-organization 
    77102      !!---------------------------------------------------------------------- 
    78       INTEGER :: ji, jj, ikt 
    79 #if defined key_sed_off 
    80       INTEGER  :: numblt          
    81       INTEGER  :: nummsh    
    82       REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdta 
    83 #endif 
     103      INTEGER :: ji, jj, ikt, ierr 
    84104      !!---------------------------------------------------------------------- 
    85105 
     
    90110      CALL ctl_opn( numsed, 'sediment.output', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 
    91111 
    92       WRITE(numsed,*) 
    93       WRITE(numsed,*) '                 L S C E - C E A' 
    94       WRITE(numsed,*) '                 SEDIMENT model' 
    95       WRITE(numsed,*) '                version 2.0  (2007) ' 
    96       WRITE(numsed,*) 
    97       WRITE(numsed,*) 
    98  
    99      
    100       WRITE(numsed,*) ' sed_init : Initialization of sediment module  ' 
    101       WRITE(numsed,*) ' ' 
    102  
    103       !                                ! Allocate LOBSTER arrays 
    104       IF( sed_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sed_ini: unable to allocate sediment model arrays' ) 
    105  
     112      IF (lwp) THEN 
     113         WRITE(numsed,*) 
     114         WRITE(numsed,*) '                 PISCES framework' 
     115         WRITE(numsed,*) '                 SEDIMENT model' 
     116         WRITE(numsed,*) '                version 3.0  (2018) ' 
     117         WRITE(numsed,*) 
     118         WRITE(numsed,*) 
     119      ENDIF 
     120 
     121      IF(lwp) WRITE(numsed,*) ' sed_init : Initialization of sediment module  ' 
     122      IF(lwp) WRITE(numsed,*) ' ' 
     123 
     124      ! Read sediment Namelist 
     125      !------------------------- 
     126      CALL sed_init_nam 
     127 
     128      ! Allocate SEDIMENT arrays 
     129      ierr =        sed_alloc() 
     130      ierr = ierr + sed_oce_alloc() 
     131      ierr = ierr + sed_adv_alloc()  
     132      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'sed_ini: unable to allocate sediment model arrays' ) 
    106133 
    107134      ! Determination of sediments number of points and allocate global variables 
    108 #if defined key_sed_off 
    109       ! Reading Netcdf Pisces file for deepest water layer thickness [m] 
    110       ! This data will be used as mask to merdge space in 1D vector 
    111       !---------------------------------------------------------------- 
    112  
    113       CALL iom_open ( 'mesh_mask'  , nummsh )    
    114       CALL iom_open ( 'e3tbot'     , numblt ) 
    115  
    116       ! mask sediment points for outputs 
    117       CALL iom_get( nummsh, jpdom_data, 'tmask' , tmask ) 
    118       CALL iom_get( nummsh, jpdom_data, 'mbathy', sbathy ) 
    119        
    120       ! longitude/latitude values 
    121       CALL iom_get ( nummsh, jpdom_data,'nav_lon', glamt(:,:) ) 
    122       CALL iom_get ( nummsh, jpdom_data,'nav_lat', gphit(:,:) ) 
    123      
    124       ! bottom layer thickness 
    125       ALLOCATE( zdta(jpi,jpj) ) 
    126       CALL iom_get  ( numblt, jpdom_data, 'E3TBOT', zdta(:,:) ) 
    127       epkbot(:,:) = 0. 
    128       DO jj = 1, jpj 
    129          DO ji = 1, jpi 
    130             ikt = MAX( INT( sbathy(ji,jj) ), 1 ) 
    131             IF( tmask(ji,jj,ikt) == 1 ) epkbot(ji,jj) = zdta(ji,jj) 
    132          ENDDO 
    133       ENDDO 
    134  
    135       CALL iom_close( nummsh )   
    136       CALL iom_close( numblt )  
    137  
    138       DEALLOCATE( zdta ) 
    139 #else 
    140  
    141135      epkbot(:,:) = 0. 
    142136      DO jj = 1, jpj 
     
    144138            ikt = mbkt(ji,jj)  
    145139            IF( tmask(ji,jj,ikt) == 1 ) epkbot(ji,jj) = e3t_1d(ikt) 
     140            gdepbot(ji,jj) = gdepw_0(ji,jj,ikt) 
    146141         ENDDO 
    147142      ENDDO 
    148 #endif      
    149  
    150143 
    151144      ! computation of total number of ocean points 
    152145      !-------------------------------------------- 
    153       jpoce  = COUNT( epkbot(:,:) > 0. ) 
    154  
     146      sedmask = 0. 
     147      IF ( COUNT( epkbot(:,:) > 0. ) == 0 ) THEN  
     148          sedmask = 0. 
     149      ELSE 
     150          sedmask = 1. 
     151      ENDIF  
     152      jpoce  = MAX( COUNT( epkbot(:,:) > 0. ) , 1 ) 
    155153 
    156154      ! Allocate memory size of global variables 
     
    159157      ALLOCATE( rainrm(jpoce,jpsol) )        ;  ALLOCATE( rainrg(jpoce,jpsol) )        ;  ALLOCATE( raintg(jpoce) )  
    160158      ALLOCATE( dzdep(jpoce) )               ;  ALLOCATE( iarroce(jpoce) )             ;  ALLOCATE( dzkbot(jpoce) ) 
    161       ALLOCATE( temp(jpoce) )                ;  ALLOCATE( salt(jpoce) )                 
     159      ALLOCATE( zkbot(jpoce) )               ;  ALLOCATE( db(jpoce,jpksed) ) 
     160      ALLOCATE( temp(jpoce) )                ;  ALLOCATE( salt(jpoce) )   
     161      ALLOCATE( diff(jpoce,jpksed,jpwat ) )  ;  ALLOCATE( irrig(jpoce, jpksed) ) 
     162      ALLOCATE( wacc(jpoce) )                ;  ALLOCATE( fecratio(jpoce) ) 
    162163      ALLOCATE( press(jpoce) )               ;  ALLOCATE( densSW(jpoce) )  
    163       ALLOCATE( hipor(jpoce,jpksed) )        ;  ALLOCATE( co3por(jpoce,jpksed) )  
     164      ALLOCATE( hipor(jpoce,jpksed) )        ;  ALLOCATE( co3por(jpoce,jpksed) ) 
    164165      ALLOCATE( dz3d(jpoce,jpksed) )         ;  ALLOCATE( volw3d(jpoce,jpksed) )       ;  ALLOCATE( vols3d(jpoce,jpksed) ) 
     166      ALLOCATE( sedligand(jpoce, jpksed) ) 
    165167 
    166168      ! Initialization of global variables 
    167       pwcp  (:,:,:) = 0.  ;  pwcp0 (:,:,:) = 0.  ; pwcp_dta  (:,:) = 0.   
    168       solcp (:,:,:) = 0.  ;  solcp0(:,:,:) = 0.  ; rainrm_dta(:,:) = 0. 
    169       rainrm(:,:  ) = 0.  ;  rainrg(:,:  ) = 0.  ; raintg    (:  ) = 0.  
    170       dzdep (:    ) = 0.  ;  iarroce(:   ) = 0   ; dzkbot    (:  ) = 0. 
    171       temp  (:    ) = 0.  ;  salt   (:   ) = 0.   
    172       press (:    ) = 0.  ;  densSW (:   ) = 0.  
    173       hipor (:,:  ) = 0.  ;  co3por (:,: ) = 0.   
    174       dz3d  (:,:  ) = 0.  ;  volw3d (:,: ) = 0.  ;  vols3d   (:,:) = 0.  
     169      pwcp  (:,:,:) = 0.   ;  pwcp0 (:,:,:) = 0.  ; pwcp_dta  (:,:) = 0.   
     170      solcp (:,:,:) = 0.   ;  solcp0(:,:,:) = 0.  ; rainrm_dta(:,:) = 0. 
     171      rainrm(:,:  ) = 0.   ;  rainrg(:,:  ) = 0.  ; raintg    (:  ) = 0.  
     172      dzdep (:    ) = 0.   ;  iarroce(:   ) = 0   ; dzkbot    (:  ) = 0. 
     173      temp  (:    ) = 0.   ;  salt   (:   ) = 0.  ; zkbot     (:  ) = 0. 
     174      press (:    ) = 0.   ;  densSW (:   ) = 0.  ; db        (:,:) = 0.  
     175      hipor (:,:  ) = 0.   ;  co3por (:,: ) = 0.  ; irrig     (:,:) = 0.  
     176      dz3d  (:,:  ) = 0.   ;  volw3d (:,: ) = 0.  ; vols3d    (:,:) = 0.  
     177      fecratio(:)   = 1E-5  
     178      sedligand(:,:) = 0.6E-9 
    175179 
    176180      ! Chemical variables       
     
    178182      ALLOCATE( ak1ps (jpoce) )  ;  ALLOCATE( ak2ps  (jpoce) )  ;  ALLOCATE( ak3ps (jpoce) ) ;  ALLOCATE( aksis (jpoce) )     
    179183      ALLOCATE( aksps (jpoce) )  ;  ALLOCATE( ak12s  (jpoce) )  ;  ALLOCATE( ak12ps(jpoce) ) ;  ALLOCATE( ak123ps(jpoce) )     
    180       ALLOCATE( borats(jpoce) )  ;  ALLOCATE( calcon2(jpoce) )     
     184      ALLOCATE( borats(jpoce) )  ;  ALLOCATE( calcon2(jpoce) )  ;  ALLOCATE( sieqs (jpoce) )  
     185      ALLOCATE( aks3s(jpoce) )   ;  ALLOCATE( akf3s(jpoce) )    ;  ALLOCATE( sulfats(jpoce) ) 
     186      ALLOCATE( fluorids(jpoce) ) 
    181187 
    182188      akbs  (:) = 0. ;   ak1s   (:) = 0. ;  ak2s  (:) = 0. ;   akws   (:) = 0. 
    183189      ak1ps (:) = 0. ;   ak2ps  (:) = 0. ;  ak3ps (:) = 0. ;   aksis  (:) = 0. 
    184190      aksps (:) = 0. ;   ak12s  (:) = 0. ;  ak12ps(:) = 0. ;   ak123ps(:) = 0. 
    185       borats(:) = 0. ;   calcon2(:) = 0. 
     191      borats(:) = 0. ;   calcon2(:) = 0. ;  sieqs (:) = 0. 
     192      aks3s(:)  = 0. ;   akf3s(:)   = 0. ;  sulfats(:) = 0. ;  fluorids(:) = 0. 
    186193 
    187194      ! Mass balance calculation   
     
    191198      fromsed(:,:) = 0.    ;   tosed(:,:) = 0. ;  rloss(:,:) = 0.  ;   tokbot(:,:) = 0.  
    192199 
    193       ! Read sediment Namelist 
    194       !------------------------- 
    195       CALL sed_init_nam 
    196      
    197200      ! Initialization of sediment geometry 
    198201      !------------------------------------ 
    199202      CALL sed_init_geom 
    200203 
    201  
    202       ! sets initial sediment composition 
    203       ! ( only clay or reading restart file ) 
    204       !--------------------------------------- 
    205       CALL sed_init_data 
    206  
    207  
    208       CALL sed_init_wri 
    209  
     204      ! Offline specific mode 
     205      ! --------------------- 
     206      ln_sediment_offline = .FALSE. 
     207 
     208#if defined key_sed_off 
     209      ln_sediment_offline = .TRUE. 
     210      IF (lwp) write(numsed,*) 'Sediment module is run in offline mode' 
     211      IF (lwp) write(numsed,*) 'key_sed_off is activated at compilation time' 
     212      IF (lwp) write(numsed,*) 'ln_sed_2way is forced to false' 
     213      IF (lwp) write(numsed,*) '--------------------------------------------' 
     214      ln_sed_2way = .FALSE. 
     215#endif 
     216      ! Initialisation of tracer damping 
     217      ! -------------------------------- 
     218      IF (ln_sediment_offline) THEN 
     219         CALL trc_dta_ini(jptra) 
     220         CALL trc_dmp_sed_ini 
     221      ENDIF 
    210222 
    211223   END SUBROUTINE sed_init 
    212  
    213224 
    214225   SUBROUTINE sed_init_geom 
     
    227238      !! * Modules used 
    228239      !! * local declarations 
    229       INTEGER ::  & 
    230         ji, jj, jk 
    231   
    232 #if defined key_sed_off 
    233       REAL(wp) , DIMENSION (jpi,jpj) :: zdta 
    234       INTEGER  ::  numpres 
    235 #endif 
     240      INTEGER  :: ji, jj, jk, jn 
     241      REAL(wp) :: za0, za1, zt, zw, zsum, zsur, zprof, zprofw 
     242      REAL(wp) :: ztmp1, ztmp2 
    236243      !----------------------------------------------------------  
    237244 
    238       WRITE(numsed,*) ' sed_init_geom : Initialization of sediment geometry ' 
    239       WRITE(numsed,*) ' ' 
     245      IF(lwp) WRITE(numsed,*) ' sed_init_geom : Initialization of sediment geometry ' 
     246      IF(lwp) WRITE(numsed,*) ' ' 
    240247 
    241248      ! Computation of 1D array of sediments points 
     
    250257      END DO 
    251258 
     259      IF ( indoce .EQ. 0 ) THEN 
     260         indoce = 1 
     261         iarroce(indoce) = 1 
     262      ENDIF 
     263 
    252264      IF( indoce .NE. jpoce ) THEN 
    253          WRITE(numsed,*) ' ' 
    254          WRITE(numsed,*) 'Warning : number of ocean points indoce = ', indoce, & 
    255             &     ' doesn''t match  number of point where epkbot>0  jpoce = ', jpoce 
    256          WRITE(numsed,*) ' ' 
    257          WRITE(numsed,*) ' ' 
    258          STOP 
     265         CALL ctl_stop( 'STOP', 'sed_ini: number of ocean points indoce doesn''t match  number of point' ) 
    259266      ELSE 
    260          WRITE(numsed,*) ' ' 
    261          WRITE(numsed,*) ' total number of ocean points jpoce =  ',jpoce 
    262          WRITE(numsed,*) ' ' 
    263       ENDIF 
    264  
    265 #if defined key_sed_off 
    266  
    267       ! Reading Netcdf Pisces file for deepest water layer thickness [m] 
    268       ! This data will be used as mask to merdge space in 1D vector 
    269       !---------------------------------------------------------------- 
    270       CALL iom_open ( 'pressbot'   , numpres ) 
    271        
    272       ! pressure  in  bars       
    273       CALL iom_get ( numpres, jpdom_data,'BATH', zdta(:,:) ) 
    274       CALL pack_arr ( jpoce, press(1:jpoce), zdta(1:jpi,1:jpj), iarroce(1:jpoce) ) 
    275       press(1:jpoce) = press(1:jpoce) * 1.025e-1  
    276  
    277       CALL iom_close ( numpres ) 
    278 #endif 
    279  
    280  
    281       ! mask sediment points for outputs 
    282       DO jk = 1, jpksed 
    283          tmasksed(:,:,jk) = tmask(:,:,1) 
    284       ENDDO 
     267         IF (lwp) WRITE(numsed,*) ' ' 
     268         IF (lwp) WRITE(numsed,*) ' total number of ocean points jpoce =  ',jpoce 
     269         IF (lwp) WRITE(numsed,*) ' ' 
     270      ENDIF 
    285271 
    286272      ! initialization of dzkbot in [cm] 
     
    288274      CALL pack_arr ( jpoce, dzkbot(1:jpoce), epkbot(1:jpi,1:jpj), iarroce(1:jpoce) ) 
    289275      dzkbot(1:jpoce) = dzkbot(1:jpoce) * 1.e+2  
     276      CALL pack_arr ( jpoce, zkbot(1:jpoce), gdepbot(1:jpi,1:jpj), iarroce(1:jpoce) ) 
    290277 
    291278      ! Geometry and  constants  
     
    293280      ! (1st layer= diffusive layer = pur water)  
    294281      !------------------------------------------ 
    295       dz(1)  = 0.1 
    296       dz(2)  = 0.3 
    297       dz(3)  = 0.3 
    298       dz(4)  = 0.5 
    299       dz(5)  = 0.5 
    300       dz(6)  = 0.5 
    301       dz(7)  = 1. 
    302       dz(8)  = 1. 
    303       dz(9)  = 1. 
    304       dz(10) = 2.45 
    305       dz(11) = 2.45 
     282      za1  = (  sedzmin - sedhmax / FLOAT(jpksed-1)  )                                                      & 
     283         & / ( TANH((1-sedkth)/sedacr) - sedacr/FLOAT(jpksed-1) * (  LOG( COSH( (jpksed - sedkth) / sedacr) )      & 
     284         &                                                   - LOG( COSH( ( 1  - sedkth) / sedacr) )  )  ) 
     285      za0  = sedzmin - za1 * TANH( (1-sedkth) / sedacr ) 
     286      zsur = - za0 - za1 * sedacr * LOG( COSH( (1-sedkth) / sedacr )  ) 
     287 
     288      profsedw(1) = 0.0 
     289      profsed(1) = -dz(1) / 2. 
     290      DO jk = 2, jpksed 
     291         zw = REAL( jk , wp ) 
     292         zt = REAL( jk , wp ) - 0.5_wp 
     293         profsed(jk)  = ( zsur + za0 * zt + za1 * sedacr * LOG ( COSH( (zt-sedkth) / sedacr ) )  )  
     294         profsedw(jk) = ( zsur + za0 * zw + za1 * sedacr * LOG ( COSH( (zw-sedkth) / sedacr ) )  ) 
     295      END DO 
     296 
     297      dz(1) = 0.1 
     298      DO jk = 2, jpksed 
     299         dz(jk) = profsedw(jk) - profsedw(jk-1) 
     300      END DO 
    306301 
    307302      DO jk = 1, jpksed 
     
    311306      ENDDO 
    312307 
    313       ! Depth(jk)= depth of middle of each layer 
    314       !---------------------------------------- 
    315       profsed(1) = -dz(1)/ 2. 
    316       DO jk = 2, jpksed 
    317          profsed(jk) = profsed(jk-1) + dz(jk-1) / 2. + dz(jk) / 2. 
    318       ENDDO 
    319  
    320308      !  Porosity profile [0] 
    321309      !--------------------- 
    322       por(1)  = 1. 
    323       por(2)  = 0.95 
    324       por(3)  = 0.9 
    325       por(4)  = 0.85 
    326       por(5)  = 0.81 
    327       por(6)  = 0.75 
    328       por(7)  = 0.75       
    329       por(8)  = 0.75       
    330       por(9)  = 0.75       
    331       por(10) = 0.75 
    332       por(11) = 0.75 
     310      por(1) = 1.0 
     311      DO jk = 2, jpksed 
     312         por(jk) = porinf + ( porsurf-porinf) * exp(-rhox * (profsed(jk) ) ) 
     313      END DO 
    333314  
    334315      ! inverse of  Porosity profile 
     
    353334      dz3d(1:jpoce,1) = dz(1) 
    354335 
    355  
    356336      !--------------------------------------------- 
    357337      ! Molecular weight [g/mol] for solid species 
    358338      !--------------------------------------------- 
    359  
    360339 
    361340      ! opal=sio2*0.4(h20)=28+2*16+0.4*(2+16) 
     
    371350         &              0.59 * 27. + 10. * 16. 
    372351 
     352      mol_wgt(jsfeo)  = 55.0 + 3.0 * ( 16.0 + 1.0) 
     353 
     354      mol_wgt(jsfes)  = 55.0 + 32.0 
     355 
    373356      ! for chemistry Poc : C(122)H(244)O(86)N(16)P(1) 
    374357      ! But den sity of Poc is an Hydrated material (= POC + 30H2O) 
     
    377360      mol_wgt(jspoc) = ( 122. * 12. + 355. + 120. * 16.+ & 
    378361         &                16. * 14. + 31. ) / 122. 
     362      mol_wgt(jspos) = mol_wgt(jspoc) 
     363      mol_wgt(jspor) = mol_wgt(jspoc) 
    379364 
    380365      ! CaCO3 
     
    384369      ! Density of solid material in sediment [g/cm**3] 
    385370      !------------------------------------------------ 
    386       dens = 2.6 
     371      denssol = 2.6 
    387372 
    388373      ! Initialization of diffusion coefficient as function of porosity [cm**2/s] 
    389374      !-------------------------------------------------------------------- 
    390       diff(:) = dcoef * por(:) 
     375!      DO jn = 1, jpsol 
     376!         DO jk = 1, jpksed 
     377!            DO ji = 1, jpoce 
     378!               diff(ji,jk,jn) = dcoef / ( 1.0 - 2.0 * log(por(jk)) ) 
     379!            END DO 
     380!         END DO 
     381!      END DO 
     382 
     383      ! Accumulation rate from Burwicz et al. (2011). This is used to 
     384      ! compute the flux of clays and minerals 
     385      ! -------------------------------------------------------------- 
     386      DO ji = 1, jpoce 
     387          ztmp1 = 0.117 / ( 1.0 + ( zkbot(ji) / 200.)**3 ) 
     388          ztmp2 = 0.006 / ( 1.0 + ( zkbot(ji) / 4000.)**10 ) 
     389          wacc(ji) = ztmp1 + ztmp2 
     390      END DO 
    391391 
    392392 
    393393      ! Initialization of time step as function of porosity [cm**2/s] 
    394394      !------------------------------------------------------------------ 
    395       rdtsed(2:jpksed) = dtsed / ( dens * por1(2:jpksed) ) 
    396       
    397395   END SUBROUTINE sed_init_geom 
    398396 
     
    408406      !!---------------------------------------------------------------------- 
    409407 
    410       INTEGER :: & 
    411          numnamsed = 28 
     408      INTEGER ::   numnamsed_ref = -1           !! Logical units for namelist sediment 
     409      INTEGER ::   numnamsed_cfg = -1           !! Logical units for namelist sediment 
     410      INTEGER :: ios                 ! Local integer output status for namelist read 
     411      CHARACTER(LEN=20)   ::   clname 
    412412 
    413413      TYPE PSED 
     
    422422      TYPE(PSED) , DIMENSION(jpdia2dsed) :: seddiag2d 
    423423 
    424       NAMELIST/nam_time/nfreq 
     424      NAMELIST/nam_run/nrseddt,ln_sed_2way 
     425      NAMELIST/nam_geom/jpksed, sedzmin, sedhmax, sedkth, sedacr, porsurf, porinf, rhox 
    425426      NAMELIST/nam_trased/sedsol, sedwat 
    426427      NAMELIST/nam_diased/seddiag3d, seddiag2d 
    427       NAMELIST/nam_reac/sisat, claysat, rcopal, rcclay, dcoef  
    428       NAMELIST/nam_poc/redO2, redNo3, redPo4, redC, redDnit, & 
    429          &           rcorg, o2seuil, rcorgN 
    430       NAMELIST/nam_cal/rccal 
    431       NAMELIST/nam_dc13/pdb, rc13P, rc13Ca 
    432       NAMELIST/nam_btb/dbiot 
    433       NAMELIST/nam_rst/ln_rst_sed 
    434  
    435       INTEGER :: jn, jn1 
     428      NAMELIST/nam_inorg/rcopal, dcoef, rccal, ratligc, rcligc 
     429      NAMELIST/nam_poc/redO2, redNo3, redPo4, redC, redfep, rcorgl, rcorgs,  & 
     430         &             rcorgr, rcnh4, rch2s, rcfe2, rcfeh2s, rcfes, rcfeso,  & 
     431         &             xksedo2, xksedno3, xksedfeo, xksedso4 
     432      NAMELIST/nam_btb/dbiot, ln_btbz, dbtbzsc, adsnh4, ln_irrig, xirrzsc 
     433      NAMELIST/nam_rst/ln_rst_sed, cn_sedrst_indir, cn_sedrst_outdir, cn_sedrst_in, cn_sedrst_out 
     434 
     435      INTEGER :: ji, jn, jn1 
    436436      !------------------------------------------------------- 
    437437 
    438       WRITE(numsed,*) ' sed_init_nam : Read namelists ' 
    439       WRITE(numsed,*) ' ' 
     438      IF(lwp) WRITE(numsed,*) ' sed_init_nam : Read namelists ' 
     439      IF(lwp) WRITE(numsed,*) ' ' 
    440440 
    441441      ! ryear = 1 year converted in second 
    442442      !------------------------------------ 
    443       WRITE(numsed,*) ' ' 
    444       WRITE(numsed,*) 'number of seconds in one year : ryear = ', ryear 
    445       WRITE(numsed,*) ' '      
     443      IF (lwp) THEN 
     444         WRITE(numsed,*) ' ' 
     445         WRITE(numsed,*) 'number of seconds in one year : ryear = ', ryear 
     446         WRITE(numsed,*) ' '      
     447      ENDIF 
    446448 
    447449      ! Reading namelist.sed variables 
    448450      !--------------------------------- 
    449       CALL ctl_opn( numnamsed, 'namelist.sediment', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 
    450  
    451       dtsed = rdt 
     451      clname = 'namelist_sediment' 
     452      IF(lwp) WRITE(numsed,*) ' sed_init_nam : read SEDIMENT namelist' 
     453      IF(lwp) WRITE(numsed,*) ' ~~~~~~~~~~~~~~' 
     454      CALL ctl_opn( numnamsed_ref, TRIM( clname )//'_ref', 'OLD'    , 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 
     455      CALL ctl_opn( numnamsed_cfg, TRIM( clname )//'_cfg', 'OLD'    , 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 
     456 
    452457      nitsed000 = nittrc000 
    453458      nitsedend = nitend 
    454 #if ! defined key_sed_off 
    455       nwrised   = nwritetrc 
    456 #else 
    457       nwrised   = nwrite 
    458 #endif 
    459      ! Diffraction/reaction parameters 
    460       !---------------------------------- 
    461       READ( numnamsed, nam_time ) 
    462       WRITE(numsed,*) ' namelist nam_time' 
    463  
    464 #if ! defined key_sed_off 
    465       nfreq     = 1 
    466 #endif 
    467  
    468       WRITE(numsed,*) ' sedimentation time step              dtsed      = ', dtsed 
    469       WRITE(numsed,*) ' 1st time step for sediment.          nitsed000  = ', nitsed000 
    470       WRITE(numsed,*) ' last time step for sediment.         nitsedend  = ', nitsedend 
    471       WRITE(numsed,*) ' frequency of sediment outputs        nwrised    = ', nwrised 
    472       WRITE(numsed,*) ' frequency of restoring inputs data   nfreq      = ', nfreq 
    473       WRITE(numsed,*) ' ' 
    474  
    475       REWIND( numnamsed )               ! read nattrc 
    476       READ  ( numnamsed, nam_trased ) 
     459      ! Namelist nam_run 
     460      REWIND( numnamsed_ref )              ! Namelist nam_run in reference namelist : Pisces variables 
     461      READ  ( numnamsed_ref, nam_run, IOSTAT = ios, ERR = 901) 
     462901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_run in reference namelist', lwp ) 
     463 
     464      REWIND( numnamsed_cfg )              ! Namelist nam_run in reference namelist : Pisces variables 
     465      READ  ( numnamsed_cfg, nam_run, IOSTAT = ios, ERR = 902) 
     466902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_run in configuration namelist', lwp ) 
     467 
     468      IF (lwp) THEN 
     469         WRITE(numsed,*) ' namelist nam_run' 
     470         WRITE(numsed,*) ' Nb of iterations for fast species    nrseddt = ', nrseddt 
     471         WRITE(numsed,*) ' 2-way coupling between PISCES and Sed ln_sed_2way = ', ln_sed_2way 
     472      ENDIF 
     473 
     474      REWIND( numnamsed_ref )              ! Namelist nam_geom in reference namelist : Pisces variables 
     475      READ  ( numnamsed_ref, nam_geom, IOSTAT = ios, ERR = 903) 
     476903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_geom in reference namelist', lwp ) 
     477 
     478      REWIND( numnamsed_cfg )              ! Namelist nam_geom in reference namelist : Pisces variables 
     479      READ  ( numnamsed_cfg, nam_geom, IOSTAT = ios, ERR = 904) 
     480904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_geom in configuration namelist', lwp ) 
     481 
     482      IF (lwp) THEN  
     483         WRITE(numsed,*) ' namelist nam_geom' 
     484         WRITE(numsed,*) ' Number of vertical layers            jpksed  = ', jpksed 
     485         WRITE(numsed,*) ' Minimum vertical spacing             sedzmin = ', sedzmin 
     486         WRITE(numsed,*) ' Maximum depth of the sediment        sedhmax = ', sedhmax 
     487         WRITE(numsed,*) ' Default parameter                    sedkth  = ', sedkth 
     488         WRITE(numsed,*) ' Default parameter                    sedacr  = ', sedacr 
     489         WRITE(numsed,*) ' Sediment porosity at the surface     porsurf = ', porsurf 
     490         WRITE(numsed,*) ' Sediment porosity at infinite depth  porinf  = ', porinf 
     491         WRITE(numsed,*) ' Length scale of porosity variation   rhox    = ', rhox 
     492      ENDIF 
     493 
     494      jpksedm1  = jpksed - 1 
     495      dtsed = r2dttrc 
     496 
     497      REWIND( numnamsed_ref )              ! Namelist nam_trased in reference namelist : Pisces variables 
     498      READ  ( numnamsed_ref, nam_trased, IOSTAT = ios, ERR = 905) 
     499905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_trased in reference namelist', lwp ) 
     500 
     501      REWIND( numnamsed_cfg )              ! Namelist nam_trased in reference namelist : Pisces variables 
     502      READ  ( numnamsed_cfg, nam_trased, IOSTAT = ios, ERR = 906) 
     503906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_trased in configuration namelist', lwp ) 
    477504 
    478505      DO jn = 1, jpsol 
     
    489516      END DO 
    490517 
    491       WRITE(numsed,*) ' namelist nam_trased' 
    492       WRITE(numsed,*) ' ' 
    493       DO jn = 1, jptrased 
    494          WRITE(numsed,*) 'name of 3d output sediment field number :',jn,' : ',TRIM(sedtrcd(jn)) 
    495          WRITE(numsed,*) 'long name ', TRIM(sedtrcl(jn)) 
    496          WRITE(numsed,*) ' in unit = ', TRIM(sedtrcu(jn)) 
    497          WRITE(numsed,*) ' ' 
    498       END DO 
    499       WRITE(numsed,*) ' ' 
    500  
     518      IF (lwp) THEN 
     519         WRITE(numsed,*) ' namelist nam_trased' 
     520         WRITE(numsed,*) ' ' 
     521         DO jn = 1, jptrased 
     522            WRITE(numsed,*) 'name of 3d output sediment field number :',jn,' : ',TRIM(sedtrcd(jn)) 
     523            WRITE(numsed,*) 'long name ', TRIM(sedtrcl(jn)) 
     524            WRITE(numsed,*) ' in unit = ', TRIM(sedtrcu(jn)) 
     525            WRITE(numsed,*) ' ' 
     526         END DO 
     527         WRITE(numsed,*) ' ' 
     528      ENDIF 
     529 
     530      REWIND( numnamsed_ref )              ! Namelist nam_diased in reference namelist : Pisces variables 
     531      READ  ( numnamsed_ref, nam_diased, IOSTAT = ios, ERR = 907) 
     532907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diased in reference namelist', lwp ) 
     533 
     534      REWIND( numnamsed_cfg )              ! Namelist nam_diased in reference namelist : Pisces variables 
     535      READ  ( numnamsed_cfg, nam_diased, IOSTAT = ios, ERR = 908) 
     536908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diased in configuration namelist', lwp ) 
    501537       
    502       REWIND( numnamsed ) 
    503       READ( numnamsed, nam_diased ) 
    504  
    505538      DO jn = 1, jpdia3dsed 
    506539         seddia3d(jn) = seddiag3d(jn)%snamesed 
     
    515548      END DO 
    516549 
    517       WRITE(numsed,*) ' namelist nam_diased' 
    518       WRITE(numsed,*) ' ' 
    519       DO jn = 1, jpdia3dsed 
    520          WRITE(numsed,*) 'name of 3D output diag number :',jn, ' : ', TRIM(seddia3d(jn)) 
    521          WRITE(numsed,*) 'long name ', TRIM(seddia3l(jn)) 
    522          WRITE(numsed,*) ' in unit = ',TRIM(seddia3u(jn)) 
    523          WRITE(numsed,*) ' ' 
    524       END DO 
    525  
    526       DO jn = 1, jpdia2dsed 
    527          WRITE(numsed,*) 'name of 2D output diag number :',jn, ' : ', TRIM(seddia2d(jn)) 
    528          WRITE(numsed,*) 'long name ', TRIM(seddia2l(jn)) 
    529          WRITE(numsed,*) ' in unit = ',TRIM(seddia2u(jn)) 
    530          WRITE(numsed,*) ' ' 
    531       END DO 
    532  
    533       WRITE(numsed,*) ' ' 
    534  
    535  
    536       ! Diffraction/reaction parameters 
     550      IF (lwp) THEN 
     551         WRITE(numsed,*) ' namelist nam_diased' 
     552         WRITE(numsed,*) ' ' 
     553         DO jn = 1, jpdia3dsed 
     554            WRITE(numsed,*) 'name of 3D output diag number :',jn, ' : ', TRIM(seddia3d(jn)) 
     555            WRITE(numsed,*) 'long name ', TRIM(seddia3l(jn)) 
     556            WRITE(numsed,*) ' in unit = ',TRIM(seddia3u(jn)) 
     557            WRITE(numsed,*) ' ' 
     558         END DO 
     559 
     560         DO jn = 1, jpdia2dsed 
     561            WRITE(numsed,*) 'name of 2D output diag number :',jn, ' : ', TRIM(seddia2d(jn)) 
     562            WRITE(numsed,*) 'long name ', TRIM(seddia2l(jn)) 
     563            WRITE(numsed,*) ' in unit = ',TRIM(seddia2u(jn)) 
     564            WRITE(numsed,*) ' ' 
     565         END DO 
     566 
     567         WRITE(numsed,*) ' ' 
     568      ENDIF 
     569 
     570      ! Inorganic chemistry parameters 
    537571      !---------------------------------- 
    538       REWIND( numnamsed ) 
    539       READ( numnamsed, nam_reac ) 
    540       WRITE(numsed,*) ' namelist nam_reac' 
    541       WRITE(numsed,*) ' saturation for si    sisat   = ', sisat 
    542       WRITE(numsed,*) ' saturation for clay  claysat = ', claysat 
    543       WRITE(numsed,*) ' reactivity for Si    rcopal  = ', rcopal 
    544       WRITE(numsed,*) ' reactivity for clay  rcclay  = ', rcclay 
    545       WRITE(numsed,*) ' diff. coef for por.  dcoef   = ', dcoef 
    546       WRITE(numsed,*) ' ' 
    547  
     572      REWIND( numnamsed_ref )              ! Namelist nam_inorg in reference namelist : Pisces variables 
     573      READ  ( numnamsed_ref, nam_inorg, IOSTAT = ios, ERR = 909) 
     574909   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_inorg in reference namelist', lwp ) 
     575 
     576      REWIND( numnamsed_cfg )              ! Namelist nam_inorg in reference namelist : Pisces variables 
     577      READ  ( numnamsed_cfg, nam_inorg, IOSTAT = ios, ERR = 910) 
     578910   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_inorg in configuration namelist', lwp ) 
     579 
     580      IF (lwp) THEN 
     581         WRITE(numsed,*) ' namelist nam_inorg' 
     582         WRITE(numsed,*) ' reactivity for Si      rcopal  = ', rcopal 
     583         WRITE(numsed,*) ' diff. coef for por.    dcoef   = ', dcoef 
     584         WRITE(numsed,*) ' reactivity for calcite rccal   = ', rccal 
     585         WRITE(numsed,*) ' L/C ratio in POC       ratligc = ', ratligc 
     586         WRITE(numsed,*) ' reactivity for ligands rcligc  = ', rcligc 
     587         WRITE(numsed,*) ' ' 
     588      ENDIF 
    548589 
    549590      ! Unity conversion to get saturation conc. psat in [mol.l-1] 
    550591      ! and reactivity rc in  [l.mol-1.s-1] 
    551592      !---------------------------------------------------------- 
    552       sat_sil   = sisat   * 1.e-6    
    553       sat_clay  = claysat * 1.e-6 
    554  
    555       reac_sil  = rcopal / ryear      
    556       reac_clay  = rcclay / ryear 
    557  
     593      reac_sil   = rcopal / ryear      
     594      reac_ligc  = rcligc / ryear 
    558595 
    559596      ! Additional parameter linked to POC/O2/No3/Po4 
    560597      !---------------------------------------------- 
    561       REWIND( numnamsed ) 
    562       READ( numnamsed, nam_poc ) 
    563       WRITE(numsed,*) ' namelist nam_poc' 
    564       WRITE(numsed,*) ' Redfield coef for oxy     redO2   = ', redO2 
    565       WRITE(numsed,*) ' Redfield coef for no3     redNo3  = ', redNo3 
    566       WRITE(numsed,*) ' Redfield coef for po4     redPo4  = ', redPo4 
    567       WRITE(numsed,*) ' Redfield coef for carbone redC    = ', redC 
    568       WRITE(numsed,*) ' Redfield coef for denitri redDnit = ', redDnit 
    569       WRITE(numsed,*) ' reactivity for POC/O2     rcorg   = ', rcorg 
    570       WRITE(numsed,*) ' threshold O2 concen.      o2seuil = ', o2seuil 
    571       WRITE(numsed,*) ' reactivity for POC/NO3    rcorgN  = ', rcorgN 
    572       WRITE(numsed,*) ' ' 
    573  
    574  
    575       so2ut  = redO2   / redC 
    576       srno3  = redNo3  / redC 
    577       spo4r  = redPo4  / redC 
    578       srDnit = redDnit / redC 
    579       sthro2 = o2seuil * 1.e-6 ! threshold O2 concen. in [mol.l-1] 
     598      REWIND( numnamsed_ref )              ! Namelist nam_poc in reference namelist : Pisces variables 
     599      READ  ( numnamsed_ref, nam_poc, IOSTAT = ios, ERR = 911) 
     600911   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_poc in reference namelist', lwp ) 
     601 
     602      REWIND( numnamsed_cfg )              ! Namelist nam_poc in reference namelist : Pisces variables 
     603      READ  ( numnamsed_cfg, nam_poc, IOSTAT = ios, ERR = 912) 
     604912   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_poc in configuration namelist', lwp ) 
     605 
     606      IF (lwp) THEN 
     607         WRITE(numsed,*) ' namelist nam_poc' 
     608         WRITE(numsed,*) ' Redfield coef for oxy            redO2    = ', redO2 
     609         WRITE(numsed,*) ' Redfield coef for no3            redNo3   = ', redNo3 
     610         WRITE(numsed,*) ' Redfield coef for po4            redPo4   = ', redPo4 
     611         WRITE(numsed,*) ' Redfield coef for carbon         redC     = ', redC 
     612         WRITE(numsed,*) ' Ration for iron bound P          redfep   = ', redfep 
     613         WRITE(numsed,*) ' reactivity for labile POC        rcorgl   = ', rcorgl 
     614         WRITE(numsed,*) ' reactivity for semi-refract. POC rcorgs   = ', rcorgs 
     615         WRITE(numsed,*) ' reactivity for refractory POC    rcorgr   = ', rcorgr 
     616         WRITE(numsed,*) ' reactivity for NH4               rcnh4    = ', rcnh4 
     617         WRITE(numsed,*) ' reactivity for H2S               rch2s    = ', rch2s 
     618         WRITE(numsed,*) ' reactivity for Fe2+              rcfe2    = ', rcfe2 
     619         WRITE(numsed,*) ' reactivity for FeOH/H2S          rcfeh2s  = ', rcfeh2s 
     620         WRITE(numsed,*) ' reactivity for Fe2+/H2S          rcfes    = ', rcfes 
     621         WRITE(numsed,*) ' reactivity for FeS/O2            rcfeso   = ', rcfeso 
     622         WRITE(numsed,*) ' Half-sat. cste for oxic remin    xksedo2  = ', xksedo2 
     623         WRITE(numsed,*) ' Half-sat. cste for denit.        xksedno3 = ', xksedno3 
     624         WRITE(numsed,*) ' Half-sat. cste for iron remin    xksedfeo = ', xksedfeo 
     625         WRITE(numsed,*) ' Half-sat. cste for SO4 remin     xksedso4 = ', xksedso4 
     626         WRITE(numsed,*) ' ' 
     627      ENDIF 
     628 
     629 
     630      so2ut  = redO2    / redC 
     631      srno3  = redNo3   / redC 
     632      spo4r  = redPo4   / redC 
     633      srDnit = ( (redO2 + 32. ) * 0.8 - redNo3 - redNo3 * 0.6 ) / redC 
    580634      ! reactivity rc in  [l.mol-1.s-1] 
    581       reac_poc  = rcorg / ryear 
    582       reac_no3 = rcorgN  / ryear 
    583  
    584  
    585       ! Carbonate parameters 
    586       !--------------------- 
    587       READ( numnamsed, nam_cal ) 
    588       WRITE(numsed,*) ' namelist  nam_cal' 
    589       WRITE(numsed,*) ' reactivity for calcite rccal = ', rccal 
    590       WRITE(numsed,*) ' ' 
     635      reac_pocl  = rcorgl / ryear 
     636      reac_pocs  = rcorgs / ryear 
     637      reac_pocr  = rcorgr / ryear 
     638      reac_nh4   = rcnh4  / ryear 
     639      reac_h2s   = rch2s  / ryear 
     640      reac_fe2   = rcfe2  / ryear 
     641      reac_feh2s = rcfeh2s/ ryear 
     642      reac_fes   = rcfes  / ryear 
     643      reac_feso  = rcfeso / ryear 
    591644 
    592645      ! reactivity rc in  [l.mol-1.s-1]       
    593646      reac_cal = rccal / ryear 
    594647 
    595  
    596       ! C13  parameters 
    597       !---------------- 
    598       READ( numnamsed, nam_dc13 ) 
    599       WRITE(numsed,*) ' namelist nam_dc13 '  
    600       WRITE(numsed,*) ' 13C/12C in PD Belemnite        PDB    = ', pdb 
    601       WRITE(numsed,*) ' 13C/12C in POC = rc13P*PDB     rc13P  = ', rc13P 
    602       WRITE(numsed,*) ' 13C/12C in CaCO3 = rc13Ca*PDB  rc13Ca = ', rc13Ca 
    603       WRITE(numsed,*) ' ' 
    604  
    605        
    606648      ! Bioturbation parameter 
    607649      !------------------------ 
    608       READ( numnamsed, nam_btb ) 
    609       WRITE(numsed,*) ' namelist nam_btb '  
    610       WRITE(numsed,*) ' coefficient for bioturbation   dbiot    = ', dbiot 
    611       WRITE(numsed,*) ' ' 
    612  
    613       ! Unity convertion to get bioturb coefficient in [cm2.s-1] 
    614       db = dbiot / ( ryear * 1000. ) 
     650      REWIND( numnamsed_ref )              ! Namelist nam_btb in reference namelist : Pisces variables 
     651      READ  ( numnamsed_ref, nam_btb, IOSTAT = ios, ERR = 913) 
     652913   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_btb in reference namelist', lwp ) 
     653 
     654      REWIND( numnamsed_cfg )              ! Namelist nam_btb in reference namelist : Pisces variables 
     655      READ  ( numnamsed_cfg, nam_btb, IOSTAT = ios, ERR = 914) 
     656914   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_btb in configuration namelist', lwp ) 
     657 
     658      IF (lwp) THEN 
     659         WRITE(numsed,*) ' namelist nam_btb '  
     660         WRITE(numsed,*) ' coefficient for bioturbation      dbiot    = ', dbiot 
     661         WRITE(numsed,*) ' Depth varying bioturbation        ln_btbz  = ', ln_btbz 
     662         WRITE(numsed,*) ' coefficient for btb attenuation   dbtbzsc  = ', dbtbzsc 
     663         WRITE(numsed,*) ' Adsorption coefficient of NH4     adsnh4   = ', adsnh4 
     664         WRITE(numsed,*) ' Bioirrigation in sediment         ln_irrig = ', ln_irrig 
     665         WRITE(numsed,*) ' coefficient for irrig attenuation xirrzsc  = ', xirrzsc 
     666         WRITE(numsed,*) ' ' 
     667      ENDIF 
    615668 
    616669      ! Initial value (t=0) for sediment pore water and solid components 
    617670      !---------------------------------------------------------------- 
    618       READ( numnamsed, nam_rst ) 
    619       WRITE(numsed,*) ' namelist  nam_rst '  
    620       WRITE(numsed,*) '  boolean term for restart (T or F) ln_rst_sed = ', ln_rst_sed  
    621       WRITE(numsed,*) ' ' 
    622  
    623       CLOSE( numnamsed ) 
     671      REWIND( numnamsed_ref )              ! Namelist nam_rst in reference namelist : Pisces variables 
     672      READ  ( numnamsed_ref, nam_rst, IOSTAT = ios, ERR = 915) 
     673915   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_rst in reference namelist', lwp ) 
     674 
     675      REWIND( numnamsed_cfg )              ! Namelist nam_rst in reference namelist : Pisces variables 
     676      READ  ( numnamsed_cfg, nam_rst, IOSTAT = ios, ERR = 916) 
     677916   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_rst in configuration namelist', lwp ) 
     678 
     679      IF (lwp) THEN 
     680         WRITE(numsed,*) ' namelist  nam_rst '  
     681         WRITE(numsed,*) '  boolean term for restart (T or F) ln_rst_sed = ', ln_rst_sed  
     682         WRITE(numsed,*) ' ' 
     683      ENDIF 
     684      nn_dtsed = nn_dttrc 
     685 
     686      CLOSE( numnamsed_cfg ) 
     687      CLOSE( numnamsed_ref ) 
    624688 
    625689   END SUBROUTINE sed_init_nam 
    626690 
    627    SUBROUTINE sed_init_data 
    628       !!---------------------------------------------------------------------- 
    629       !!                   ***  ROUTINE sed_init_data  *** 
    630       !! 
    631       !! ** Purpose :  Initialization of sediment module 
    632       !!               - sets initial sediment composition 
    633       !!                 ( only clay or reading restart file ) 
    634       !! 
    635       !!   History : 
    636       !!        !  06-07  (C. Ethe)  original 
    637       !!---------------------------------------------------------------------- 
    638   
    639       ! local variables 
    640       INTEGER :: & 
    641          ji, jk, zhipor 
    642  
    643       !-------------------------------------------------------------------- 
    644   
    645  
    646       IF( .NOT. ln_rst_sed ) THEN 
    647  
    648           WRITE(numsed,*) ' Initilization of default values of sediment components' 
    649  
    650          ! default values for initial pore water concentrations [mol/l] 
    651          pwcp(:,:,:) = 0. 
    652          ! default value for initial solid component (fraction of dry weight dim=[0]) 
    653          ! clay 
    654          solcp(:,:,:) = 0. 
    655          solcp(:,2:jpksed,jsclay) = 1. 
    656  
    657          ! Initialization of [h+] and [co3--] 
    658  
    659          zhipor = 8. 
    660          ! Initialization of [h+] in mol/kg 
    661          DO jk = 1, jpksed 
    662             DO ji = 1, jpoce 
    663                hipor (ji,jk) = 10.**( -1. * zhipor ) 
    664             ENDDO 
    665          ENDDO 
    666  
    667          co3por(:,:) = 0. 
    668  
    669       ELSE    
    670    
    671          WRITE(numsed,*) ' Initilization of Sediment components from restart' 
    672  
    673          CALL sed_rst_read 
    674  
    675       ENDIF 
    676  
    677  
    678       ! Load initial Pisces Data for bot. wat. Chem and fluxes 
    679       CALL sed_dta ( nitsed000 )  
    680  
    681       ! Initialization of chemical constants 
    682       CALL sed_chem ( nitsed000 ) 
    683  
    684       ! Stores initial sediment data for mass balance calculation 
    685       pwcp0 (1:jpoce,1:jpksed,1:jpwat ) = pwcp (1:jpoce,1:jpksed,1:jpwat )  
    686       solcp0(1:jpoce,1:jpksed,1:jpsol ) = solcp(1:jpoce,1:jpksed,1:jpsol)  
    687  
    688       ! Conversion of [h+] in mol/Kg to get it in mol/l ( multiplication by density) 
    689       DO jk = 1, jpksed 
    690          hipor(1:jpoce,jk) = hipor(1:jpoce,jk) * densSW(1:jpoce) 
    691       ENDDO 
    692  
    693  
    694       ! In default case - no restart - sedco3 is run to initiate [h+] and [co32-] 
    695       ! Otherwise initiate values of pH and co3 read in restart 
    696       IF( .NOT. ln_rst_sed ) THEN 
    697          ! sedco3 is run to initiate[h+] [co32-] in mol/l of solution 
    698          CALL sed_co3 ( nitsed000 ) 
    699  
    700       ENDIF 
    701              
    702    END SUBROUTINE sed_init_data 
    703  
    704    SUBROUTINE sed_init_wri 
    705  
    706       INTEGER :: jk 
    707  
    708       WRITE(numsed,*)' ' 
    709       WRITE(numsed,*)'======== Write summary of initial state   ============' 
    710       WRITE(numsed,*)' ' 
    711       WRITE(numsed,*)' ' 
    712       WRITE(numsed,*)'-------------------------------------------------------------------' 
    713       WRITE(numsed,*)' Initial Conditions ' 
    714       WRITE(numsed,*)'-------------------------------------------------------------------' 
    715       WRITE(numsed,*)'dzm = dzkbot minimum to calculate ', 0. 
    716       WRITE(numsed,*)'Local zone : jpi, jpj : ',jpi, jpj 
    717       WRITE(numsed,*)'jpoce = ',jpoce,' nbtot pts = ',jpij,' nb earth pts = ',jpij - jpoce 
    718       WRITE(numsed,*)'sublayer thickness dz(1) [cm] : ', dz(1) 
    719       WRITE(numsed,*)'Coeff diff for k=1 (cm2/s) : ',diff(1) 
    720       WRITE(numsed,*)' nb solid comp : ',jpsol 
    721       WRITE(numsed,*)'(1=opal,2=clay,3=POC,4=CaCO3)' 
    722       WRITE(numsed,*)'weight mol 1,2,3,4' 
    723       WRITE(numsed,'(4(F0.2,3X))')mol_wgt(jsopal),mol_wgt(jsclay),mol_wgt(jspoc),mol_wgt(jscal) 
    724       WRITE(numsed,*)'nb dissolved comp',jpwat 
    725       WRITE(numsed,*)'(1=silicic acid,2="dissolved" clay,3=O2,4=DIC,5=Nitrate,& 
    726          &6=Phosphates,7=Alk))' 
    727       WRITE(numsed,*)'Psat (umol/l) for silicic Acid and "dissolved" clay' 
    728       WRITE(numsed,'(2(F0.2,3X))') sat_sil * 1e+6,  sat_clay * 1e+6 
    729       WRITE(numsed,*)'reaction rate rc for Op/si,Clay,POC/O2,caco3, POC/No3 (an-1)' 
    730       WRITE(numsed,'(5(F0.2,3X))') reac_sil * ryear, reac_clay * ryear, reac_poc * ryear, & 
    731                                    reac_cal * ryear, reac_no3 * ryear 
    732       WRITE(numsed,*)'redfield coef C,O,N P Dit ' 
    733       WRITE(numsed,'(5(F0.2,3X))')1./spo4r,so2ut/spo4r,srno3/spo4r,spo4r/spo4r,srDnit/spo4r 
    734       WRITE(numsed,*)'threshold for stating denitrification [mol/l]' 
    735       WRITE(numsed,'(1PE8.2)') sthrO2 
    736       WRITE(numsed,*)'-------------------------------------------------------------------' 
    737       WRITE(numsed,*)'Min-Max-Mean' 
    738       WRITE(numsed,*)'For each variable : min, max, moy value observed on selected local zone' 
    739       WRITE(numsed,*)'-------------------------------------------------------------------' 
    740       WRITE(numsed,*)'thickness of the last wet layer dzkbot [m]' 
    741       WRITE(numsed,'(3(F0.2,3X))') MINVAL(dzkbot(1:jpoce))/100.,MAXVAL(dzkbot(1:jpoce))/100.,& 
    742          &SUM(dzkbot(1:jpoce))/jpoce/100. 
    743       WRITE(numsed,*)'temp [°C]' 
    744       WRITE(numsed,'(3(F0.2,3X))') MINVAL(temp(1:jpoce)),MAXVAL(temp(1:jpoce)),& 
    745          &                         SUM(temp(1:jpoce))/jpoce 
    746       WRITE(numsed,*)'salt o/oo' 
    747       WRITE(numsed,'(3(F0.2,3X))')MINVAL(salt(1:jpoce)),MAXVAL(salt(1:jpoce)),& 
    748          &                        SUM(salt(1:jpoce))/jpoce 
    749 #if defined key_sed_off 
    750       WRITE(numsed,*)'pressure [bar] (depth in m is about 10*pressure)' 
    751       WRITE(numsed,'(3(F0.2,3X))') MINVAL(press(1:jpoce)),MAXVAL(press(1:jpoce)),& 
    752          &                         SUM(press(1:jpoce))/jpoce 
    753 #endif 
    754       WRITE(numsed,*)'density of Sea Water' 
    755       WRITE(numsed,'(3(F0.2,3X))') MINVAL(densSW(1:jpoce)), MAXVAL(densSW(1:jpoce)),& 
    756          &                         SUM(densSW(1:jpoce))/jpoce 
    757       WRITE(numsed,*)'' 
    758       WRITE(numsed,*)'     Dissolved Components ' 
    759       WRITE(numsed,*)'     =====================' 
    760       WRITE(numsed,*)'[Si(OH)4] dissolved (1)(k=1)(µmol/l)(and min value in mol/kg  of solution)' 
    761       WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwsil))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwsil))*1.e+6,& 
    762          &                         SUM(pwcp(1:jpoce,1,jwsil))*1.e+6/jpoce,& 
    763          &                         MINVAL(pwcp(1:jpoce,1,jwsil)*1.e+6/densSW(1:jpoce)) 
    764       WRITE(numsed,*)'[O2]    dissolved (3) (k=1)(µmol/l)(and min value in mol/kg  of solution)' 
    765       WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwoxy))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwoxy))*1.e+6,& 
    766          &SUM(pwcp(1:jpoce,1,jwoxy))*1.e+6/jpoce,& 
    767          &MINVAL(pwcp(1:jpoce,1,jwoxy)*1.e+6/densSW(1:jpoce)) 
    768       WRITE(numsed,*)'[DIC]    dissolved (4) (k=1)(µmol/l)(and min value in mol/kg  of solution)' 
    769       WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwdic))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwdic))*1.e+6,& 
    770          &SUM(pwcp(1:jpoce,1,jwdic))*1.e+6/jpoce,& 
    771          &MINVAL(pwcp(1:jpoce,1,jwdic)*1.e+6/densSW(1:jpoce)) 
    772       WRITE(numsed,*)'[NO3]    dissolved (5) (k=1)(µmol/l)(and min value in mol/kg  of solution)' 
    773       WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwno3))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwno3))*1.e+6,& 
    774          &SUM(pwcp(1:jpoce,1,jwno3))*1.e+6/jpoce,& 
    775          &MINVAL(pwcp(1:jpoce,1,jwno3)*1.e+6/densSW(1:jpoce)) 
    776       WRITE(numsed,*)'[PO4]    dissolved (6) (k=1)(µmol/l)(and min value in mol/kg  of solution)' 
    777       WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwpo4))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwpo4))*1.e+6,& 
    778          &SUM(pwcp(1:jpoce,1,jwpo4))*1.e+6/jpoce,& 
    779          &MINVAL(pwcp(1:jpoce,1,jwpo4)*1.e+6/densSW(1:jpoce)) 
    780       WRITE(numsed,*)'[Alk]    dissolved (7) (k=1)(µequi)(and min value in mol/kg  of solution)' 
    781       WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwalk))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwalk))*1.e+6,& 
    782          &SUM(pwcp(1:jpoce,1,jwalk))*1.e+6/jpoce,& 
    783          &MINVAL(pwcp(1:jpoce,1,jwalk)*1.e+6/densSW(1:jpoce)) 
    784       WRITE(numsed,*)'[DIC13]  dissolved (8) (k=1)(µmol/l)(and min value in mol/kg  of solution)' 
    785       WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwc13))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwc13))*1.e+6,& 
    786          &SUM(pwcp(1:jpoce,1,jwc13))*1.e+6/jpoce,& 
    787          &MINVAL(pwcp(1:jpoce,1,jwc13)*1.e+6/densSW(1:jpoce)) 
    788       WRITE(numsed,*)'' 
    789       WRITE(numsed,*)'     Solid Components ' 
    790       WRITE(numsed,*)'     =====================' 
    791       WRITE(numsed,*)'nmole of Opale rained per dt' 
    792       WRITE(numsed,'(3(1PE9.3,2X))') MINVAL(rainrm(1:jpoce,jsopal))*dtsed,MAXVAL(rainrm(1:jpoce,jsopal))*dtsed,& 
    793          &SUM(rainrm(1:jpoce,1))*dtsed/jpoce 
    794       WRITE(numsed,*)'nmole of Clay rained per dt' 
    795       WRITE(numsed,'(3(1PE9.3,2X))') MINVAL(rainrm(1:jpoce,jsclay))*dtsed,MAXVAL(rainrm(1:jpoce,jsclay))*dtsed,& 
    796          &SUM(rainrm(1:jpoce,jsclay))*dtsed/jpoce 
    797       WRITE(numsed,*)'nmole of POC rained per dt' 
    798       WRITE(numsed,'(3(1PE9.3,2X))') MINVAL(rainrm(1:jpoce,jspoc))*dtsed,MAXVAL(rainrm(1:jpoce,jspoc))*dtsed,& 
    799          &SUM(rainrm(1:jpoce,jspoc))*dtsed/jpoce 
    800       WRITE(numsed,*)'nmole of CaCO3 rained per dt' 
    801       WRITE(numsed,'(3(1PE9.3,2X))') MINVAL(rainrm(1:jpoce,jscal))*dtsed,MAXVAL(rainrm(1:jpoce,jscal))*dtsed,& 
    802          &SUM(rainrm(1:jpoce,jscal))*dtsed/jpoce 
    803       WRITE(numsed,*)' ' 
    804       WRITE(numsed,*)'Weight frac of opal rained (%) ' 
    805       WRITE(numsed,'(3(F0.3,7X))') MINVAL(rainrg(1:jpoce,jsopal)/raintg(1:jpoce))*100.,& 
    806          &MAXVAL(rainrg(1:jpoce,jsopal)/raintg(1:jpoce))*100.,& 
    807          & SUM(rainrg(1:jpoce,jsopal)/raintg(1:jpoce))*100./jpoce 
    808       WRITE(numsed,*)'Weight frac of clay  rained (%) ' 
    809       WRITE(numsed,'(3(F0.3,7X))') MINVAL(rainrg(1:jpoce,jsclay)/raintg(1:jpoce))*100.,& 
    810          &MAXVAL(rainrg(1:jpoce,jsclay)/raintg(1:jpoce))*100.,& 
    811          &SUM(rainrg(1:jpoce,jsclay)/raintg(1:jpoce))*100./jpoce 
    812       WRITE(numsed,*)'Weight frac of POC rained (%)' 
    813       WRITE(numsed,'(3(F0.3,7X))') MINVAL(rainrg(1:jpoce,jspoc)/raintg(1:jpoce))*100.,& 
    814          &MAXVAL(rainrg(1:jpoce,jspoc)/raintg(1:jpoce))*100.,& 
    815          &SUM(rainrg(1:jpoce,jspoc)/raintg(1:jpoce))*100./jpoce 
    816       WRITE(numsed,*)'Weight frac of CaCO3  rained (%)' 
    817       WRITE(numsed,'(3(F0.3,7X))') MINVAL(rainrg(1:jpoce,jscal)/raintg(1:jpoce))*100.,& 
    818          &MAXVAL(rainrg(1:jpoce,jscal)/raintg(1:jpoce))*100.,& 
    819          &SUM(rainrg(1:jpoce,jscal)/raintg(1:jpoce))*100./jpoce 
    820       WRITE(numsed,*)'' 
    821       WRITE(numsed,*)' Thickness of sediment layer added by particule rain, dzdep cm ' 
    822       WRITE(numsed,*)' ==============================================================' 
    823       WRITE(numsed,'(3(F0.5,2X))') MINVAL(dzdep(1:jpoce)),MAXVAL(dzdep(1:jpoce)),SUM(dzdep(1:jpoce))/jpoce 
    824       WRITE(numsed,*)'' 
    825       WRITE(numsed,*)' chemical constants K1,pK1,K2,pK2,Kw,pKw and Kb pKb (min max) [mol/kgsol] or [mol/kgsol]2 ' 
    826       WRITE(numsed,*)' =========================================================================================' 
    827       WRITE(numsed,'(4(1PE10.3,2X))')MINVAL(ak1s(1:jpoce)),MAXVAL(ak1s(1:jpoce)),-LOG10(MINVAL(ak1s(1:jpoce))),& 
    828          &-LOG10(MAXVAL(ak1s(1:jpoce))) 
    829       WRITE(numsed,'(4(1PE10.3,2X))')MINVAL(ak2s(1:jpoce)),MAXVAL(ak2s(1:jpoce)),-LOG10(MINVAL(ak2s(1:jpoce))),& 
    830          &-LOG10(MAXVAL(ak2s(1:jpoce))) 
    831       WRITE(numsed,'(4(1PE10.3,2X))')MINVAL(akws(1:jpoce)),MAXVAL(akws(1:jpoce)),-LOG10(MINVAL(akws(1:jpoce))),& 
    832          &-LOG10(MAXVAL(akws(1:jpoce))) 
    833       WRITE(numsed,'(4(1PE10.3,2X))')MINVAL(akbs(1:jpoce)),MAXVAL(akbs(1:jpoce)),-LOG10(MINVAL(akbs(1:jpoce))),& 
    834          &-LOG10(MAXVAL(akbs(1:jpoce))) 
    835       WRITE(numsed,*)'and boron' 
    836       WRITE(numsed,'(2(1PE10.3,2X))')MINVAL(borats(1:jpoce)),MAXVAL(borats(1:jpoce)) 
    837       WRITE(numsed,*)'' 
    838       WRITE(numsed,*)' Compo of initial sediment for point jpoce' 
    839       WRITE(numsed,*)' =========================================' 
    840       WRITE(numsed,*)'solcp(1),   solcp(2),   solcp(3),   solcp(4),   hipor,      pH,         co3por' 
    841       DO jk = 1,jpksed 
    842          WRITE(numsed,'(7(1PE10.3,2X))')solcp(jpoce,jk,jsopal),solcp(jpoce,jk,jsclay),solcp(jpoce,jk,jspoc),solcp(jpoce,jk,jscal),& 
    843             &       hipor(jpoce,jk),-LOG10(hipor(jpoce,jk)/densSW(jpoce)),co3por(jpoce,jk) 
    844       ENDDO 
    845       WRITE(numsed,'(A82)')'pwcp(1),    pwcp(2),    pwcp(3),    pwcp(4),    pwcp(5),    pwcp(6),    pwcp(7)' 
    846       DO jk = 1, jpksed 
    847          WRITE(numsed,'(7(1PE10.3,2X))')pwcp(jpoce,jk,jwsil),pwcp(jpoce,jk,jwoxy),pwcp(jpoce,jk,jwdic),& 
    848             & pwcp(jpoce,jk,jwno3),pwcp(jpoce,jk,jwpo4),pwcp(jpoce,jk,jwalk),pwcp(jpoce,jk,jwc13) 
    849       ENDDO 
    850       WRITE(numsed,*) ' ' 
    851       WRITE(numsed,*) ' End Of Initialization ' 
    852       WRITE(numsed,*) ' ' 
    853 ! 
    854    END SUBROUTINE sed_init_wri 
    855 #else 
    856    !!---------------------------------------------------------------------- 
    857    !!   Dummy module :                      NO Sediment model 
    858    !!---------------------------------------------------------------------- 
    859    !! $Id$ 
    860 CONTAINS 
    861    SUBROUTINE sed_ini              ! Empty routine 
    862    END SUBROUTINE sed_ini 
    863 #endif 
    864  
    865  
    866691END MODULE sedini 
Note: See TracChangeset for help on using the changeset viewer.