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.
sedini.F90 in NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/TOP/PISCES/SED – NEMO

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/TOP/PISCES/SED/sedini.F90 @ 11624

Last change on this file since 11624 was 11624, checked in by acc, 5 years ago

Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Substantive changes required to replace all namelists with internal files. These are the key changes only; to compile and run tests all REWIND and CLOSE operations on the (no longer) units have to be removed. These changes affect many more files but can be scripted so are not included here in order to make a later merge easier. The scripts used to prepare code for testing are included on: wiki:2019WP/ENHANCE-04_AndrewC-reporting/Internal_Namelists. With these additional changes this code passes most SETTE tests but the AGRIF preprocessor does not currently accept the new allocatable character strings. To be investigated.
  • Property svn:keywords set to Id
File size: 28.3 KB
RevLine 
[3443]1MODULE sedini
2   !!======================================================================
3   !!              ***  MODULE  sedini  ***
4   !! Sediment : define sediment variables
5   !!=====================================================================
6
7   !!----------------------------------------------------------------------
8   !!   sed_init    : initialization, namelist read, and parameters control
9   !!----------------------------------------------------------------------
10   !! * Modules used
11   USE sed     ! sediment global variable
[10222]12   USE sed_oce
[3443]13   USE sedarr
[10222]14   USE sedadv
[10333]15   USE trc_oce, ONLY : nn_dttrc
[10222]16   USE trcdmp_sed
17   USE trcdta
[3443]18   USE iom
[10222]19   USE lib_mpp         ! distribued memory computing library
[3443]20
21
22   IMPLICIT NONE
23   PRIVATE
24
25   !! Module variables
26   REAL(wp)    ::  &
[10222]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)    ::  &
[3443]38      rcopal  =   40.  ,  &  !: reactivity for si    [l.mol-1.an-1]
39      dcoef   =  8.e-6       !: diffusion coefficient (*por)   [cm**2/s]
40
[10222]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
[3443]60
61   REAL(wp)    ::  &
[10222]62      rccal   = 1000.,      & !: reactivity for calcite         [l.mol-1.an-1]
63      rcligc  = 1.E-4         !: L/C ratio in POC
[3443]64
[10222]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.
[3443]68   REAL(wp)    ::  &
[10222]69      ryear = 365. * 24. * 3600. !:  1 year converted in second
[3443]70
[10222]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 /
[3443]73
[10222]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 /
[3443]76
[10222]77
78
[3443]79   !! *  Routine accessibility
80   PUBLIC sed_init          ! routine called by opa.F90
81
[5215]82   !! $Id$
[3443]83CONTAINS
84
85
86   SUBROUTINE sed_init
87      !!----------------------------------------------------------------------
88      !!                   ***  ROUTINE sed_init  ***
89      !!
90      !! ** Purpose :  Initialization of sediment module
91      !!               - Reading namelist
92      !!               - Read the deepest water layer thickness
93      !!                 ( using as mask ) in Netcdf file
94      !!               - Convert unity if necessary
95      !!               - sets initial sediment composition
96      !!                 ( only clay or reading restart file )
97      !!               - sets sediment grid, porosity and others constants
98      !!
99      !!   History :
100      !!        !  04-10  (N. Emprin, M. Gehlen )  Original code
101      !!        !  06-07  (C. Ethe)  Re-organization
102      !!----------------------------------------------------------------------
[10222]103      INTEGER :: ji, jj, ikt, ierr
[3443]104      !!----------------------------------------------------------------------
105
106
107      ! Reading namelist.sed variables
108      !---------------------------------------
109
110      CALL ctl_opn( numsed, 'sediment.output', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
111
[10222]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
[3443]120
[10222]121      IF(lwp) WRITE(numsed,*) ' sed_init : Initialization of sediment module  '
122      IF(lwp) WRITE(numsed,*) ' '
[3443]123
[10222]124      ! Read sediment Namelist
125      !-------------------------
126      CALL sed_init_nam
[3443]127
[10222]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' )
[3443]133
134      ! Determination of sediments number of points and allocate global variables
135      epkbot(:,:) = 0.
136      DO jj = 1, jpj
137         DO ji = 1, jpi
138            ikt = mbkt(ji,jj) 
[4292]139            IF( tmask(ji,jj,ikt) == 1 ) epkbot(ji,jj) = e3t_1d(ikt)
[10222]140            gdepbot(ji,jj) = gdepw_0(ji,jj,ikt)
[3443]141         ENDDO
142      ENDDO
143
144      ! computation of total number of ocean points
145      !--------------------------------------------
[10222]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 )
[3443]153
154      ! Allocate memory size of global variables
155      ALLOCATE( pwcp (jpoce,jpksed,jpwat) )  ;  ALLOCATE( pwcp0 (jpoce,jpksed,jpwat) ) ;  ALLOCATE( pwcp_dta  (jpoce,jpwat) )
156      ALLOCATE( solcp(jpoce,jpksed,jpsol) )  ;  ALLOCATE( solcp0(jpoce,jpksed,jpsol) ) ;  ALLOCATE( rainrm_dta(jpoce,jpsol) )
157      ALLOCATE( rainrm(jpoce,jpsol) )        ;  ALLOCATE( rainrg(jpoce,jpsol) )        ;  ALLOCATE( raintg(jpoce) ) 
158      ALLOCATE( dzdep(jpoce) )               ;  ALLOCATE( iarroce(jpoce) )             ;  ALLOCATE( dzkbot(jpoce) )
[10222]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) )
[3443]163      ALLOCATE( press(jpoce) )               ;  ALLOCATE( densSW(jpoce) ) 
[10222]164      ALLOCATE( hipor(jpoce,jpksed) )        ;  ALLOCATE( co3por(jpoce,jpksed) )
[3443]165      ALLOCATE( dz3d(jpoce,jpksed) )         ;  ALLOCATE( volw3d(jpoce,jpksed) )       ;  ALLOCATE( vols3d(jpoce,jpksed) )
[10222]166      ALLOCATE( sedligand(jpoce, jpksed) )
[3443]167
168      ! Initialization of global variables
[10222]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
[3443]179
180      ! Chemical variables     
181      ALLOCATE( akbs  (jpoce) )  ;  ALLOCATE( ak1s   (jpoce) )  ;  ALLOCATE( ak2s  (jpoce) ) ;  ALLOCATE( akws  (jpoce) )     
182      ALLOCATE( ak1ps (jpoce) )  ;  ALLOCATE( ak2ps  (jpoce) )  ;  ALLOCATE( ak3ps (jpoce) ) ;  ALLOCATE( aksis (jpoce) )   
183      ALLOCATE( aksps (jpoce) )  ;  ALLOCATE( ak12s  (jpoce) )  ;  ALLOCATE( ak12ps(jpoce) ) ;  ALLOCATE( ak123ps(jpoce) )   
[10222]184      ALLOCATE( borats(jpoce) )  ;  ALLOCATE( calcon2(jpoce) )  ;  ALLOCATE( sieqs (jpoce) ) 
185      ALLOCATE( aks3s(jpoce) )   ;  ALLOCATE( akf3s(jpoce) )    ;  ALLOCATE( sulfats(jpoce) )
186      ALLOCATE( fluorids(jpoce) )
[3443]187
188      akbs  (:) = 0. ;   ak1s   (:) = 0. ;  ak2s  (:) = 0. ;   akws   (:) = 0.
189      ak1ps (:) = 0. ;   ak2ps  (:) = 0. ;  ak3ps (:) = 0. ;   aksis  (:) = 0.
190      aksps (:) = 0. ;   ak12s  (:) = 0. ;  ak12ps(:) = 0. ;   ak123ps(:) = 0.
[10222]191      borats(:) = 0. ;   calcon2(:) = 0. ;  sieqs (:) = 0.
192      aks3s(:)  = 0. ;   akf3s(:)   = 0. ;  sulfats(:) = 0. ;  fluorids(:) = 0.
[3443]193
194      ! Mass balance calculation 
195      ALLOCATE( fromsed(jpoce, jpsol) ) ; ALLOCATE( tosed(jpoce, jpsol) ) ;  ALLOCATE( rloss(jpoce, jpsol) )
196      ALLOCATE( tokbot (jpoce, jpwat) ) 
197
198      fromsed(:,:) = 0.    ;   tosed(:,:) = 0. ;  rloss(:,:) = 0.  ;   tokbot(:,:) = 0. 
199
200      ! Initialization of sediment geometry
201      !------------------------------------
202      CALL sed_init_geom
203
[10222]204      ! Offline specific mode
205      ! ---------------------
206      ln_sediment_offline = .FALSE.
[3443]207
[10222]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
[3443]222
223   END SUBROUTINE sed_init
224
225   SUBROUTINE sed_init_geom
226      !!----------------------------------------------------------------------
227      !!                   ***  ROUTINE sed_init_geom  ***
228      !!
229      !! ** Purpose :  Initialization of sediment geometry
230      !!               - Read the deepest water layer thickness
231      !!                 ( using as mask ) in Netcdf file
232      !!               - sets sediment grid, porosity and molecular weight
233      !!                 and others constants
234      !!
235      !!   History :
236      !!        !  06-07  (C. Ethe)  Original
237      !!----------------------------------------------------------------------
238      !! * Modules used
239      !! * local declarations
[10222]240      INTEGER  :: ji, jj, jk, jn
241      REAL(wp) :: za0, za1, zt, zw, zsum, zsur, zprof, zprofw
242      REAL(wp) :: ztmp1, ztmp2
[3443]243      !----------------------------------------------------------
244
[10222]245      IF(lwp) WRITE(numsed,*) ' sed_init_geom : Initialization of sediment geometry '
246      IF(lwp) WRITE(numsed,*) ' '
[3443]247
248      ! Computation of 1D array of sediments points
249      indoce = 0
250      DO jj = 1, jpj
251         DO ji = 1, jpi
252            IF (  epkbot(ji,jj) > 0. ) THEN
253               indoce          = indoce + 1
254               iarroce(indoce) = (jj - 1) * jpi + ji
255            ENDIF
256         END DO
257      END DO
258
[10222]259      IF ( indoce .EQ. 0 ) THEN
260         indoce = 1
261         iarroce(indoce) = 1
262      ENDIF
263
[3443]264      IF( indoce .NE. jpoce ) THEN
[10222]265         CALL ctl_stop( 'STOP', 'sed_ini: number of ocean points indoce doesn''t match  number of point' )
[3443]266      ELSE
[10222]267         IF (lwp) WRITE(numsed,*) ' '
268         IF (lwp) WRITE(numsed,*) ' total number of ocean points jpoce =  ',jpoce
269         IF (lwp) WRITE(numsed,*) ' '
[3443]270      ENDIF
271
272      ! initialization of dzkbot in [cm]
273      !------------------------------------------------   
274      CALL pack_arr ( jpoce, dzkbot(1:jpoce), epkbot(1:jpi,1:jpj), iarroce(1:jpoce) )
275      dzkbot(1:jpoce) = dzkbot(1:jpoce) * 1.e+2 
[10222]276      CALL pack_arr ( jpoce, zkbot(1:jpoce), gdepbot(1:jpi,1:jpj), iarroce(1:jpoce) )
[3443]277
278      ! Geometry and  constants
279      ! sediment layer thickness [cm]
280      ! (1st layer= diffusive layer = pur water)
281      !------------------------------------------
[10222]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 )  )
[3443]287
[10222]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
301
[3443]302      DO jk = 1, jpksed
303         DO ji = 1, jpoce
304            dz3d(ji,jk) = dz(jk)
305         END DO
306      ENDDO
307
308      !  Porosity profile [0]
309      !---------------------
[10222]310      por(1) = 1.0
311      DO jk = 2, jpksed
312         por(jk) = porinf + ( porsurf-porinf) * exp(-rhox * (profsed(jk) ) )
313      END DO
[3443]314 
315      ! inverse of  Porosity profile
316      !-----------------------------
317      por1(:) = 1. - por(:)
318
319      ! Volumes of pore water and solid fractions (vector and array)
320      !     WARNING : volw(1) and vols(1) are sublayer volums
321      volw(:) = dz(:) * por(:)
322      vols(:) = dz(:) * por1(:)
323
324      ! temporary new value for dz3d(:,1)
325      dz3d(1:jpoce,1) = dzkbot(1:jpoce)
326
327      ! WARNING : volw3d(:,1) and vols3d(:,1) are deepest water column volums
328      DO jk = 1, jpksed
329         volw3d(1:jpoce,jk) = dz3d(1:jpoce,jk) * por (jk)
330         vols3d(1:jpoce,jk) = dz3d(1:jpoce,jk) * por1(jk)
331      ENDDO
332
333      ! Back to the old sublayer vlaue for dz3d(:,1)
334      dz3d(1:jpoce,1) = dz(1)
335
336      !---------------------------------------------
337      ! Molecular weight [g/mol] for solid species
338      !---------------------------------------------
339
340      ! opal=sio2*0.4(h20)=28+2*16+0.4*(2+16)
341      !---------------------------------------
342      mol_wgt(jsopal) = 28. + 2. * 16. + 0.4 * ( 2. + 16. ) 
343
344      !  clay
345      !  some kind of Illit (according to Pape)
346      !  K0.58(Al 1.38 Fe(III)0.37Fe(II)0.04Mg0.34)[(OH)2|(Si3.41Al0.59)O10]
347      !--------------------------------------------------------------------
348      mol_wgt(jsclay) = 0.58 * 39. + 1.38 * 27. + ( 0.37 + 0.04 ) * 56.+ &
349         &              0.34 * 24. + 2. * ( 16. + 1. ) + 3.41 * 38. +    &
350         &              0.59 * 27. + 10. * 16.
351
[10222]352      mol_wgt(jsfeo)  = 55.0 + 3.0 * ( 16.0 + 1.0)
353
354      mol_wgt(jsfes)  = 55.0 + 32.0
355
[3443]356      ! for chemistry Poc : C(122)H(244)O(86)N(16)P(1)
357      ! But den sity of Poc is an Hydrated material (= POC + 30H2O)
358      ! So C(122)H(355)O(120)N(16)P(1)
359      !------------------------------------------------------------
360      mol_wgt(jspoc) = ( 122. * 12. + 355. + 120. * 16.+ &
361         &                16. * 14. + 31. ) / 122.
[10222]362      mol_wgt(jspos) = mol_wgt(jspoc)
363      mol_wgt(jspor) = mol_wgt(jspoc)
[3443]364
365      ! CaCO3
366      !---------
367      mol_wgt(jscal) = 40. + 12. + 3. * 16.
368
369      ! Density of solid material in sediment [g/cm**3]
370      !------------------------------------------------
[10222]371      denssol = 2.6
[3443]372
373      ! Initialization of diffusion coefficient as function of porosity [cm**2/s]
374      !--------------------------------------------------------------------
[10222]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
[3443]382
[10222]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
[3443]391
[10222]392
[3443]393      ! Initialization of time step as function of porosity [cm**2/s]
394      !------------------------------------------------------------------
395   END SUBROUTINE sed_init_geom
396
397   SUBROUTINE sed_init_nam
398      !!----------------------------------------------------------------------
399      !!                   ***  ROUTINE sed_init_nam  ***
400      !!
401      !! ** Purpose :  Initialization of sediment geometry
402      !!               - Reading namelist and defines constants variables
403      !!
404      !!   History :
405      !!        !  06-07  (C. Ethe)  Original
406      !!----------------------------------------------------------------------
407
[11624]408      CHARACTER(:), ALLOCATABLE ::   numnamsed_ref           !! Character buffer for reference namelist sediment
409      CHARACTER(:), ALLOCATABLE ::   numnamsed_cfg           !! Character buffer for configuration namelist sediment
[10222]410      INTEGER :: ios                 ! Local integer output status for namelist read
411      CHARACTER(LEN=20)   ::   clname
[3443]412
413      TYPE PSED
414         CHARACTER(len = 20)  :: snamesed   !: short name
415         CHARACTER(len = 80 ) :: lnamesed   !: long name
416         CHARACTER(len = 20 ) :: unitsed    !: unit
417      END TYPE PSED
418
419      TYPE(PSED) , DIMENSION(jpsol     ) :: sedsol
420      TYPE(PSED) , DIMENSION(jpwat     ) :: sedwat
421      TYPE(PSED) , DIMENSION(jpdia3dsed) :: seddiag3d
422      TYPE(PSED) , DIMENSION(jpdia2dsed) :: seddiag2d
423
[10222]424      NAMELIST/nam_run/nrseddt,ln_sed_2way
425      NAMELIST/nam_geom/jpksed, sedzmin, sedhmax, sedkth, sedacr, porsurf, porinf, rhox
[3443]426      NAMELIST/nam_trased/sedsol, sedwat
427      NAMELIST/nam_diased/seddiag3d, seddiag2d
[10222]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
[3443]434
[10222]435      INTEGER :: ji, jn, jn1
[3443]436      !-------------------------------------------------------
437
[10222]438      IF(lwp) WRITE(numsed,*) ' sed_init_nam : Read namelists '
439      IF(lwp) WRITE(numsed,*) ' '
[3443]440
441      ! ryear = 1 year converted in second
442      !------------------------------------
[10222]443      IF (lwp) THEN
444         WRITE(numsed,*) ' '
445         WRITE(numsed,*) 'number of seconds in one year : ryear = ', ryear
446         WRITE(numsed,*) ' '     
447      ENDIF
[3443]448
449      ! Reading namelist.sed variables
450      !---------------------------------
[10222]451      clname = 'namelist_sediment'
452      IF(lwp) WRITE(numsed,*) ' sed_init_nam : read SEDIMENT namelist'
453      IF(lwp) WRITE(numsed,*) ' ~~~~~~~~~~~~~~'
[11624]454      CALL load_nml( numnamsed_ref, TRIM( clname )//'_ref', numout, .FALSE. )
455      CALL load_nml( numnamsed_cfg, TRIM( clname )//'_cfg', numout, .FALSE. )
[3443]456
457      nitsed000 = nittrc000
458      nitsedend = nitend
[10222]459      ! Namelist nam_run
460      READ  ( numnamsed_ref, nam_run, IOSTAT = ios, ERR = 901)
[11536]461901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_run in reference namelist' )
[3443]462
[10222]463      READ  ( numnamsed_cfg, nam_run, IOSTAT = ios, ERR = 902)
[11536]464902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_run in configuration namelist' )
[3443]465
[10222]466      IF (lwp) THEN
467         WRITE(numsed,*) ' namelist nam_run'
468         WRITE(numsed,*) ' Nb of iterations for fast species    nrseddt = ', nrseddt
469         WRITE(numsed,*) ' 2-way coupling between PISCES and Sed ln_sed_2way = ', ln_sed_2way
470      ENDIF
[3443]471
[10362]472      IF ( ln_p5z .AND. ln_sed_2way ) CALL ctl_stop( '2 ways coupling with sediment cannot be activated with PISCES-QUOTA' )
473
[10222]474      READ  ( numnamsed_ref, nam_geom, IOSTAT = ios, ERR = 903)
[11536]475903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_geom in reference namelist' )
[3443]476
[10222]477      READ  ( numnamsed_cfg, nam_geom, IOSTAT = ios, ERR = 904)
[11536]478904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_geom in configuration namelist' )
[10222]479
480      IF (lwp) THEN
481         WRITE(numsed,*) ' namelist nam_geom'
482         WRITE(numsed,*) ' Number of vertical layers            jpksed  = ', jpksed
483         WRITE(numsed,*) ' Minimum vertical spacing             sedzmin = ', sedzmin
484         WRITE(numsed,*) ' Maximum depth of the sediment        sedhmax = ', sedhmax
485         WRITE(numsed,*) ' Default parameter                    sedkth  = ', sedkth
486         WRITE(numsed,*) ' Default parameter                    sedacr  = ', sedacr
487         WRITE(numsed,*) ' Sediment porosity at the surface     porsurf = ', porsurf
488         WRITE(numsed,*) ' Sediment porosity at infinite depth  porinf  = ', porinf
489         WRITE(numsed,*) ' Length scale of porosity variation   rhox    = ', rhox
490      ENDIF
491
492      jpksedm1  = jpksed - 1
493      dtsed = r2dttrc
494
495      READ  ( numnamsed_ref, nam_trased, IOSTAT = ios, ERR = 905)
[11536]496905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_trased in reference namelist' )
[10222]497
498      READ  ( numnamsed_cfg, nam_trased, IOSTAT = ios, ERR = 906)
[11536]499906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_trased in configuration namelist' )
[10222]500
[3443]501      DO jn = 1, jpsol
502         sedtrcd(jn) = sedsol(jn)%snamesed
503         sedtrcl(jn) = sedsol(jn)%lnamesed
504         sedtrcu(jn) = sedsol(jn)%unitsed
505      END DO
506
507      DO jn = 1, jpwat
508         jn1 = jn + jpsol
509         sedtrcd(jn1) = sedwat(jn)%snamesed
510         sedtrcl(jn1) = sedwat(jn)%lnamesed
511         sedtrcu(jn1) = sedwat(jn)%unitsed
512      END DO
513
[10222]514      IF (lwp) THEN
515         WRITE(numsed,*) ' namelist nam_trased'
[3443]516         WRITE(numsed,*) ' '
[10222]517         DO jn = 1, jptrased
518            WRITE(numsed,*) 'name of 3d output sediment field number :',jn,' : ',TRIM(sedtrcd(jn))
519            WRITE(numsed,*) 'long name ', TRIM(sedtrcl(jn))
520            WRITE(numsed,*) ' in unit = ', TRIM(sedtrcu(jn))
521            WRITE(numsed,*) ' '
522         END DO
523         WRITE(numsed,*) ' '
524      ENDIF
[3443]525
[10222]526      READ  ( numnamsed_ref, nam_diased, IOSTAT = ios, ERR = 907)
[11536]527907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diased in reference namelist' )
[10222]528
529      READ  ( numnamsed_cfg, nam_diased, IOSTAT = ios, ERR = 908)
[11536]530908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diased in configuration namelist' )
[3443]531     
532      DO jn = 1, jpdia3dsed
533         seddia3d(jn) = seddiag3d(jn)%snamesed
534         seddia3l(jn) = seddiag3d(jn)%lnamesed
535         seddia3u(jn) = seddiag3d(jn)%unitsed
536      END DO
537
538      DO jn = 1, jpdia2dsed
539         seddia2d(jn) = seddiag2d(jn)%snamesed
540         seddia2l(jn) = seddiag2d(jn)%lnamesed
541         seddia2u(jn) = seddiag2d(jn)%unitsed
542      END DO
543
[10222]544      IF (lwp) THEN
545         WRITE(numsed,*) ' namelist nam_diased'
[3443]546         WRITE(numsed,*) ' '
[10222]547         DO jn = 1, jpdia3dsed
548            WRITE(numsed,*) 'name of 3D output diag number :',jn, ' : ', TRIM(seddia3d(jn))
549            WRITE(numsed,*) 'long name ', TRIM(seddia3l(jn))
550            WRITE(numsed,*) ' in unit = ',TRIM(seddia3u(jn))
551            WRITE(numsed,*) ' '
552         END DO
[3443]553
[10222]554         DO jn = 1, jpdia2dsed
555            WRITE(numsed,*) 'name of 2D output diag number :',jn, ' : ', TRIM(seddia2d(jn))
556            WRITE(numsed,*) 'long name ', TRIM(seddia2l(jn))
557            WRITE(numsed,*) ' in unit = ',TRIM(seddia2u(jn))
558            WRITE(numsed,*) ' '
559         END DO
560
[3443]561         WRITE(numsed,*) ' '
[10222]562      ENDIF
[3443]563
[10222]564      ! Inorganic chemistry parameters
565      !----------------------------------
566      READ  ( numnamsed_ref, nam_inorg, IOSTAT = ios, ERR = 909)
[11536]567909   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_inorg in reference namelist' )
[3443]568
[10222]569      READ  ( numnamsed_cfg, nam_inorg, IOSTAT = ios, ERR = 910)
[11536]570910   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_inorg in configuration namelist' )
[3443]571
[10222]572      IF (lwp) THEN
573         WRITE(numsed,*) ' namelist nam_inorg'
574         WRITE(numsed,*) ' reactivity for Si      rcopal  = ', rcopal
575         WRITE(numsed,*) ' diff. coef for por.    dcoef   = ', dcoef
576         WRITE(numsed,*) ' reactivity for calcite rccal   = ', rccal
577         WRITE(numsed,*) ' L/C ratio in POC       ratligc = ', ratligc
578         WRITE(numsed,*) ' reactivity for ligands rcligc  = ', rcligc
579         WRITE(numsed,*) ' '
580      ENDIF
[3443]581
582      ! Unity conversion to get saturation conc. psat in [mol.l-1]
583      ! and reactivity rc in  [l.mol-1.s-1]
584      !----------------------------------------------------------
[10222]585      reac_sil   = rcopal / ryear     
586      reac_ligc  = rcligc / ryear
[3443]587
588      ! Additional parameter linked to POC/O2/No3/Po4
589      !----------------------------------------------
[10222]590      READ  ( numnamsed_ref, nam_poc, IOSTAT = ios, ERR = 911)
[11536]591911   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_poc in reference namelist' )
[3443]592
[10222]593      READ  ( numnamsed_cfg, nam_poc, IOSTAT = ios, ERR = 912)
[11536]594912   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_poc in configuration namelist' )
[3443]595
[10222]596      IF (lwp) THEN
597         WRITE(numsed,*) ' namelist nam_poc'
598         WRITE(numsed,*) ' Redfield coef for oxy            redO2    = ', redO2
599         WRITE(numsed,*) ' Redfield coef for no3            redNo3   = ', redNo3
600         WRITE(numsed,*) ' Redfield coef for po4            redPo4   = ', redPo4
601         WRITE(numsed,*) ' Redfield coef for carbon         redC     = ', redC
602         WRITE(numsed,*) ' Ration for iron bound P          redfep   = ', redfep
603         WRITE(numsed,*) ' reactivity for labile POC        rcorgl   = ', rcorgl
604         WRITE(numsed,*) ' reactivity for semi-refract. POC rcorgs   = ', rcorgs
605         WRITE(numsed,*) ' reactivity for refractory POC    rcorgr   = ', rcorgr
606         WRITE(numsed,*) ' reactivity for NH4               rcnh4    = ', rcnh4
607         WRITE(numsed,*) ' reactivity for H2S               rch2s    = ', rch2s
608         WRITE(numsed,*) ' reactivity for Fe2+              rcfe2    = ', rcfe2
609         WRITE(numsed,*) ' reactivity for FeOH/H2S          rcfeh2s  = ', rcfeh2s
610         WRITE(numsed,*) ' reactivity for Fe2+/H2S          rcfes    = ', rcfes
611         WRITE(numsed,*) ' reactivity for FeS/O2            rcfeso   = ', rcfeso
612         WRITE(numsed,*) ' Half-sat. cste for oxic remin    xksedo2  = ', xksedo2
613         WRITE(numsed,*) ' Half-sat. cste for denit.        xksedno3 = ', xksedno3
614         WRITE(numsed,*) ' Half-sat. cste for iron remin    xksedfeo = ', xksedfeo
615         WRITE(numsed,*) ' Half-sat. cste for SO4 remin     xksedso4 = ', xksedso4
616         WRITE(numsed,*) ' '
617      ENDIF
[3443]618
619
[10222]620      so2ut  = redO2    / redC
621      srno3  = redNo3   / redC
622      spo4r  = redPo4   / redC
623      srDnit = ( (redO2 + 32. ) * 0.8 - redNo3 - redNo3 * 0.6 ) / redC
624      ! reactivity rc in  [l.mol-1.s-1]
625      reac_pocl  = rcorgl / ryear
626      reac_pocs  = rcorgs / ryear
627      reac_pocr  = rcorgr / ryear
628      reac_nh4   = rcnh4  / ryear
629      reac_h2s   = rch2s  / ryear
630      reac_fe2   = rcfe2  / ryear
631      reac_feh2s = rcfeh2s/ ryear
632      reac_fes   = rcfes  / ryear
633      reac_feso  = rcfeso / ryear
[3443]634
635      ! reactivity rc in  [l.mol-1.s-1]     
636      reac_cal = rccal / ryear
637
638      ! Bioturbation parameter
639      !------------------------
[10222]640      READ  ( numnamsed_ref, nam_btb, IOSTAT = ios, ERR = 913)
[11536]641913   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_btb in reference namelist' )
[3443]642
[10222]643      READ  ( numnamsed_cfg, nam_btb, IOSTAT = ios, ERR = 914)
[11536]644914   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_btb in configuration namelist' )
[3443]645
[10222]646      IF (lwp) THEN
647         WRITE(numsed,*) ' namelist nam_btb ' 
648         WRITE(numsed,*) ' coefficient for bioturbation      dbiot    = ', dbiot
649         WRITE(numsed,*) ' Depth varying bioturbation        ln_btbz  = ', ln_btbz
650         WRITE(numsed,*) ' coefficient for btb attenuation   dbtbzsc  = ', dbtbzsc
651         WRITE(numsed,*) ' Adsorption coefficient of NH4     adsnh4   = ', adsnh4
652         WRITE(numsed,*) ' Bioirrigation in sediment         ln_irrig = ', ln_irrig
653         WRITE(numsed,*) ' coefficient for irrig attenuation xirrzsc  = ', xirrzsc
654         WRITE(numsed,*) ' '
655      ENDIF
656
[3443]657      ! Initial value (t=0) for sediment pore water and solid components
658      !----------------------------------------------------------------
[10222]659      READ  ( numnamsed_ref, nam_rst, IOSTAT = ios, ERR = 915)
[11536]660915   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_rst in reference namelist' )
[3443]661
[10222]662      READ  ( numnamsed_cfg, nam_rst, IOSTAT = ios, ERR = 916)
[11536]663916   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_rst in configuration namelist' )
[3443]664
[10222]665      IF (lwp) THEN
666         WRITE(numsed,*) ' namelist  nam_rst ' 
667         WRITE(numsed,*) '  boolean term for restart (T or F) ln_rst_sed = ', ln_rst_sed 
668         WRITE(numsed,*) ' '
[3443]669      ENDIF
[10222]670      nn_dtsed = nn_dttrc
[3443]671
672
[10222]673   END SUBROUTINE sed_init_nam
[3443]674
675END MODULE sedini
Note: See TracBrowser for help on using the repository browser.